Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients

This article is devoted to the numerical simulation of the vibrations of an elastic mechanical structure equipped with several piezoelectric patches, with applications for the control, sensing and reduction of vibrations. At first, a finite element formulation of the coupled electromechanical problem is introduced, whose originality is that provided a set of non-restrictive assumptions, the system’s electrical state is fully described by very few global discrete unknowns: only a couple of variables per piezoelectric patches, namely (1) the electric charge contained in the electrodes and (2) the voltage between the electrodes. The main advantages are (1) since the electrical state is fully discretized at the weak formulation step, any standard (elastic only) finite element formulation can be easily modified to include the piezoelectric patches (2) realistic electrical boundary conditions such that equipotentiality on the electrodes and prescribed global charges naturally appear (3) the global charge/voltage variables are intrinsically adapted to include any external electrical circuit into the electromechanical problem and to simulate shunted piezoelectric patches. The second part of the article is devoted to the introduction of a reduced-order model (ROM) of the problem, by means of a modal expansion. This leads to show that the classical efficient electromechanical coupling factors (EEMCF) naturally appear as the main parameters that master the electromechanical coupling in the ROM. Finally, all the above results are applied in the case of a cantilever beam whose vibrations are reduced by means of a resonant shunt. A finite element formulation of this problem is described. It enables to compute the system EEMCF as well as its frequency response, which are compared with experimental results, showing an excellent agreement. Copyright q 2009 John Wiley & Sons, Ltd.


INTRODUCTION
Piezoelectric materials are proposed for many applications, most of the time to couple the vibrations of a structure including the piezoelectric material with an electric circuit. This can lead to several applications like oscillators for electronic circuits, acceleration sensors and gyroscopes, sound transducers, vibration monitoring, control or energy harvesting and reduction of structural vibrations. Both direct and indirect piezoelectric effects can be used, with the piezoelectric element being used as sensors or actuators, or even both simultaneously. Except for very few applications where the structure's geometry-including the piezoelectric elements-is simple and allows an analytical solution, numerical methods, such as the finite element method (FEM), have to be used to compute the vibration response and simulate the system electromechanical behavior.
The present study proposes a finite element formulation adapted to elastic structures to which are attached piezoelectric plate-like patches with electrodes. Its originality lies in the use of only one pair of electrical variables per piezoelectric patch: the terminal voltage and the charge contained in one of the electrodes. By essence, this choice of variable is naturally adapted to simulating the behavior of the structure coupled to any electrical circuit and thus to practical engineering applications.
A pioneering work on the subject of piezoelectric finite element [1] proposed a formulation of the coupled problem between mechanical displacement and the electric potential field, taking into account both direct and inverse effects. If the piezoelectric elements are used either as sensors or actuators for active vibration control, it is worth remarking that the electrical unknowns can be condensed so that the problem to solve has the form of a standard elastic vibration problem. In particular, the actuators action on the system appears as an external forcing proportional to the applied voltage, whereas the system's stiffness is purely elastic, without any modifications due to the electric state of the patches. In the case of sensors, the problem appears with an added stiffness term due to the open-circuit condition of the piezoelectric patches, whose terminal voltage (the sensor output) is proportional to the mechanical displacement unknowns. This formulation is the basis of most numerical studies involving elastic structure with piezoelectric elements [2][3][4][5][6][7][8][9][10]. For more references, one can refer to the following review articles [11][12][13].
However, the need for modeling both sensing and actuation at the same time is required by applications where a passive (or semi-passive) electric circuit is connected to the piezoelectric elements. These so-called shunt techniques, proposed in [14], are reviewed in [15] in the case of vibration reduction. The initial models used for simulating these shunted systems were lumped, with only two degrees of freedom (d.o.f.) (one mechanical and one electrical, in most cases one vibration mode and the electric-free charge). On the contrary, richer models can be obtained via piezoelectric FEM commercial softwares, by simulating the shunt as an active controller with simultaneous sensing and actuation [16]. This latter technique is often numerically intensive. A new category of shunts has recently become of interest, often known as switch techniques, where the electrical circuit impedance is switched periodically between two values [17][18][19]. In this case, the same distinction between the lumped model, characterized by fast computation, and the full finite element simulations, more accurate but slower, are found [20]. A good compromise is to use the modal basis to project the problem and reduce the number of variables without losing the accuracy in the frequency domain [21].
In addition to simulating the vibratory behavior of the mechanical structure coupled to an electric circuit, a key issue is the optimization of the whole system, in terms of size/shape/location of the piezoelectric patches as well as the choice of the electric circuit components. In most cases, this 238 O. THOMAS, J.-F. DEÜ AND J. DUCARNE element matrices and the computation of the system short-circuit eigenmodes, which take the form of a standard elastic eigenvalue problem. To conclude, the EEMCFs for a given elastic structure with piezoelectric patches can potentially be computed in post-processing, after a modal analysis done with a standard finite element elastic code.
The present article is divided into six main parts. After the present introduction, Section 2 recalls the main equation of a standard piezoelectric vibration problem and the associated variational formulation. Section 3 gives the main electrical hypotheses that lead to introduce the global electrical variables-the voltage/charge pair associated with each piezoelectric patch-in the variational formulation and derives the associated general discretized finite element formulation. The modal expansion of the general finite element formulation is proposed in Section 4, thus, enabling to introduce and compute the problem EEMCFs. Section 5 exposes the details of a finite element discretization of a laminated beam with piezoelectric patches. To conclude, Section 6 applies all the above formulations to a beam with two collocated piezoelectric patches connected to a resonant shunt circuit. The EEMCFs are numerically evaluated as well as the vibratory response of the beam with and without the shunt. Experiments are described and an excellent agreement is obtained, thus, validating the introduced formulations.

AN ARBITRARY PIEZOELECTRIC MEDIUM
This section is devoted to the general formulation of the equations that govern the mechanical and electric state of an arbitrary piezoelectric medium, on their strong form as well as on a weak form, suitable for a finite element formulation. Standard indicial notations are adopted throughout the article: subscripts i, j, k,l denote the three-dimensional vectors and tensor components and repeated subscripts imply summation. In addition, a comma indicates a partial derivative.

Electro-mechanical local equations
We consider a piezoelectric structure occupying a domain p at the equilibrium. The structure is subjected to a prescribed displacement u d i on a part u and to a prescribed surface force density t d i on the complementary part t of its external boundary. The electric boundary conditions are defined by a prescribed electric potential d on and a surface density of free electric charges q d on the remaining part q . Thus, the total structure boundary, denoted * p , is such that * p = u ∪ t = ∪ q with u ∩ t = ∩ q = ∅. In addition, p is subjected to prescribed body forces f d i . The linearized deformation tensor and the stress tensor are denoted, respectively, by ε i j and i j . Moreover, D i denotes the electric displacement vector components and E i the electric field vector components. is the mass density, n i is the unit normal external to p and t is the time.
The local equations of the electro-mechanical coupled problem are: with appropriate initial conditions. Equation (1a) corresponds to the elastodynamic equation, Equations (1b) and (1c) are the prescribed mechanical boundary conditions, Equation (2a) corresponds to Gauss's law, Equations (2b) and (2c) are the prescribed electric conditions [30]. One can remark that Equation (2b) is valid as long as the electric displacement in the medium that surrounds p can be neglected with respect to its value inside p ; otherwise, the jump of D i n i through * p would have appeared in (2b). Moreover, no volume free electric charges are prescribed in p , since at the electric equilibrium, for classical homogeneous piezoelectric, dielectric or conducting media, the free charges are always concentrated on the boundaries, in the form of surface densities [30].
The stress tensor i j and electric displacement D i are related to the linear strain tensor ε kl and electric field E k through the converse and direct linear piezoelectric constitutive equations [31,32]: where c i jkl denote the elastic moduli at constant electric field, e ki j the piezoelectric constants and ik the dielectric permittivities at constant strain. This expression of the constitutive equation is well suited to finite element approaches. Moreover, we have the following gradient relations between the linearized strain tensor ε kl and the displacement u k , and between the electric field E k and the electric potential : Equation (6) is justified under the electrostatic approximation (the curl of E k is zero), that is valid because the characteristic time of the electric phenomena is much smaller than the one of the mechanical phenomena, which means that the eventual electromagnetic waves are neglected [32].
In this case, the curl of E k vanishes, which justifies the existence of in Equation (6) [30, 32].

Variational formulation in terms of (u i , )
The local equations in the previous section are expressed in terms of the chosen unknown fields of the boundary value problem, i.e. the mechanical displacement u i and the electric potential . In order to obtain the variational formulation associated with the local equations of the coupled system, the test-function method is applied. We introduce the spaces C u of sufficiently regular functions u i defined in p , C * u = {u i ∈ C u | u i = 0 on u } and C of sufficiently regular functions in p . ∈ C , integrating over p , applying Green's formula, taking into account Equations (1b) and (2b) as well as the constitutive laws (3) and (4) leads to: and where q r = −D i n i is the reaction free charge density that appears on . The variational formulation of the coupled electro-mechanical problem, equivalent to the strong form (Equations (1a)-(6)), consists in finding u i ∈ C u such that u i = u d i on u , ∈ C such that = d on and the reaction electric free charge density q r on , satisfying Equations (7) and (8) for all u i ∈ C * u and for all ∈ C , with appropriate initial conditions.

AN ELASTIC STRUCTURE WITH PIEZOELECTRIC PATCHES
In this section, the general formulation introduced in the previous section is applied to the special case of an elastic structure equipped with thin piezoelectric patches, either glued on its outer surface or embedded inside. It allows, with some electrical hypotheses, to replace the local electric variables and D i , defined on the whole system domain, by discrete variables, defined on each piezoelectric patch: on the one hand the electric potential difference between its electrodes and on the other hand the global charge contained in one of its electrodes.

General problem and hypotheses
We consider an elastic structure equipped with P ∈ N piezoelectric patches ( Figure 1). The domain occupied by the elastic part of the system is denoted by s . Each piezoelectric patch has the shape of a plate or a slightly curved shell, with its upper and lower surfaces covered with a very thin layer of conducting material to obtain electrodes. The pth patch, p ∈{1, . . . , P}, occupies a domain ( p) such that ( s , (1) , . . . , (P) ) is a partition of the whole domain . The pth patch upper and lower electrode surfaces are denoted, respectively, by − and its lateral surface is In the remaining part of the article, superscript ( p) will denote a quantity related to the pth piezoelectric patch.
A set of hypotheses, that apply to a wide spectrum of practical applications, as it will be seen in Section 3.4, is now formulated.
1. Only the piezoelectric patches are made of piezoelectric material. Consequently, the piezoelectric material constants e i jk vanish in s . However, we can mention that the dielectric constants ik do not vanish in s , in general [30]. 2. The piezoelectric patches are thin, with a constant thickness, denoted h ( p) for the pth patch, smaller than its characteristic longitudinal length (see Figure 1). 3. The thickness of the electrodes is much smaller than h ( p) and is thus neglected. 4. The piezoelectric patches are polarized in their transverse direction (i.e. the direction normal to the electrodes). 5. The piezoelectric patch materials are transverse isotropic, hence, they have the same mechanical and electrical properties in any longitudinal direction. 6. No free charges are localized on the lateral boundary ( p) 0 of the piezoelectric patches. 7. The electric end effects for the piezoelectric patches are neglected. 8. The electric field vector, of components E k , is normal to the electrodes and its intensity is uniform in the piezoelectric patch, so that for all p ∈{1, . . . , P}: where − is the potential difference between − and n k is the kth component of the normal unit vector to the surface of the electrodes (see Figure 1). In this case, two main couplings can be obtained: a '33' coupling (between the transverse electric field and the stress/strain components in the same direction) and a '31' coupling (between the transverse electric field and the membrane stresses/strains). The present hypothesis is justified by hypotheses 2, 5 and 7. 9. The electric displacement vector of component D i is neglected in the elastic domain s , as compared with its value in any of the piezoelectric patches ( p) . The justification of this hypothesis will be discussed in Section 3.4. However, the electric field in s does not vanish, in general.

Variational formulation in terms of
A variational formulation of the electromechanical problem associated with the system can be obtained by considering the above hypotheses 1-9 and applying Equations (7) and (8) successively to each of the P +1 subdomains of , namely the elastic part s and the piezoelectric patches ( p) , for p ∈{1, . . . , P}.

Mechanical equation.
We first apply Equation (7) successively to each subdomain of . Then, the P +1 obtained equations are added term by term considering that the integrals over each subdomain combine themselves to give a single integral on the whole domain: Consequently to the equilibrium of forces and the equality of displacement through each interface between two subdomains of , all the surface integrals that are related to the subdomain boundaries inside cancel by pair. Then, we redefine u and t to be a partition of the boundary * of where the displacements and the forces are prescribed, respectively. We also redefine C u as the space of sufficiently regular functions u i defined in the whole and C * u = {u i ∈ C u | u i = 0 on u }. Then, the piezoelectric integral term in Equation (7), with p replaced by , can be rewritten as where V ( p) is the pth piezoelectric patch potential difference, introduced in Equation (9). To prove the above identity, one has first to consider Equation (10) together with assumption 1, that implies that e i jk = 0 in s so that the corresponding integral vanishes. Finally, assumption 8 (together with Equation (9)) leads to the result.

Electrical equation.
In a second step, we apply Equation (8) to each piezoelectric patch. We chose here to prescribe the electrical potential on all the piezoelectric patches electrodes  0 . The piezoelectric term in Equation (8) is processed in the same manner as in Section 3.2.1, so that one obtains, with p replaced by ( p) : where V ( p) =  The dielectric term in Equation (8) can be rewritten as, with again p replaced by ( p) : The above identity has been obtained in the same manner as for Equations (11) and (12), by considering hypothesis 8.
With assumption 2, the volume of the pth piezoelectric patch is very close to h ( p) S ( p) , with S ( p) the average area of the patch, so that Equation (13) can be further simplified: the capacitance of the pth piezoelectric patch and 33 = ik n i n k the piezoelectric material permittivity at constant strain, in the direction normal to the electrodes. The surface integral terms in the right-hand side of Equation (8)  − . Moreover, we apply the integral form of Gauss's law to a closed surface that surrounds the pth patch. It stands that, at the electric equilibrium, where Q is the free electric charge contained inside . In our particular case, lies in s and a part of it can be inside the system surrounding medium. However, in any case, thanks to hypothesis 9, D i is neglected on , which means that the left-hand side of Equation (15) vanishes (the electric displacement in the electric cables that eventually cross is zero because they are conducting media). Consequently, the whole charge Q contained inside is zero, which means that the charge contained in the piezoelectric patch and its electrodes vanishes. As the free charges are concentrated in the electrodes, Gauss's law finally shows that the charge in the upper electrode is the opposite of the one in the lower electrode (see Figure 1). If Q ( p) is the electric charge contained in the upper electrode ( p) + of the pth piezoelectric patch, one obtains: Finally, the remaining term of the right-hand side of Equation (8) writes, by substituting for

Variational formulation.
The variational formulation of the present problem can now be written, by considering Equations (7) and (8) for each of the P patches, together with Equations (11), (12), (14) and (17).
The electromechanical problem now consists in finding u i ∈C u such that u i =u d i on u and the P electric charges Q ( p) ( p∈{1, . . . , P}), such that the associated potential differences V ( p) are and with appropriate initial conditions.

Finite element formulation
Using any finite element procedure to discretize the mechanical part of Equations (18) and (19) leads to introduce U, the vector of nodal values of u i . By introducing Q = (Q (1) Q (2) . . . Q (P) ) T and V = (V (1) V (2) . . . V (P) ) T , the column vectors of electric charges and potential differences, one finally obtains the general finite element formulation of the electromechanical problem: In the above equation, M m and K m are the mechanical mass and stiffness matrices, of size N × N , if N is the number of mechanical unknowns. K e = diag(C (1) C (2) . . .C (P) ) is a diagonal matrix filled with the P capacitances of the piezoelectric patches. K c is the electromechanical coupling matrix, of size N × P. F is the column vector of mechanical forcing, of length N . The above discretized formulation (20) can be used for a wide range of applications of a mechanical structure associated with piezoelectric patches. It is particularly adapted to the case where the piezoelectric patches are 'shunted', that is to stay connected to a passive electrical network. In this case, neither V nor Q are prescribed by the electrical network but the latter imposes only a relation between them, that can be written as f(V,V, V,Q,Q, Q) = 0. As an example, this additional relation writes, in the case of a resonant shunt connected to the pth patch and composed of a resistance R and an inductance L [24,25]: Back to the general case, the problem consists in finding U, V and Q so that: It is sometimes convenient to rewrite Equation (22a) in a (U, Q) form, using Equation (22b). One thus obtains:

245
In this latter formulation, the electromechanical coupling effect on the mechanical structure appears with an added-stiffness term, that depends on the coupling matrix K c and the patches capacitances gathered in K e [33]. This (U, Q) formulation is well suited for switch shunting applications, where the electric circuit is left open most of the time [26]. It will be used in Section 6.

Some practical system configurations
The practical system configurations that verify the hypotheses of Section 3.1 must present the following two characteristics.
• The piezoelectric patches have to be thin. It leads to hypotheses 2 to 8.
• The system has to be such that the dielectric effects of the piezoelectric patches are equivalent to an electrically independent set of capacitors. This leads to hypotheses 6 and 9.
In the case of the above first item, hypothesis 8 neglects the linear through-the-thickness dependence of the electric field, associated with flexural deformations of the patch (the so-called 'induced electric field' [34]). This hypothesis is justified if the piezoelectric patches are mainly subjected to membrane deformations, which is the case if they are thin and associated with a main elastic structure. This hypothesis has been verified for some simple geometries, by comparison with three-dimensional solutions [35,36]. For the above second item, hypothesis 9 stands that the value of the electric displacement D i can be neglected in s as compared with its value in the piezoelectric patches. We first recall that thanks to hypothesis 1, D i = ik E k in s so that D i depends only on the medium permittivity ik and the electric field E k . Most of the engineering materials, conductors as well as dielectrics, have their permittivity constants ik of the order of magnitude of 3 0 , a value much smaller than those of common piezoelectric materials, of the order of 1000 0 [30] ( 0 is the permittivity of free space). Consequently, hypothesis 9 is verified if the value of the electric field E k in s is of the same order of magnitude than in the piezoelectric patches. Since E k is the gradient of the electric potential in , its value is of the order of magnitude of the potential difference between two conductors, divided by the distance that separates them. Consequently, large electric fields are obtained if two conductors with different electric potentials are situated very close to each other. Thus, in order to obtain an electrical field in s of the same order of magnitude than in the piezoelectric patches, the electrodes of a given piezoelectric patch have to be situated as far as possible from the other patches electrodes and from the conducting parts of s of different potentials. Otherwise, one has to electrically connect all inside conductors, so that they share the same potential.
A common system configuration that leads to large values of E k in s is obtained if one wants to glue a given piezoelectric patch onto a conducting part of s without imposing the same potential in s and in the glued electrode. In this case, a very thin film of dielectric glue is used and constitutes a very thin condenser, whose capacitance cannot be neglected as compared with the capacitance of the piezoelectric patches. In this situation, the formulation of Section 3.3 is not correct, as the thin dielectric capacitance would have to be included in matrix K e .
However, this situation has to be avoided as much as possible in a practical context, because a part of the electrical energy associated with the dynamics of the system is stored in s , since neither D i nor E k vanish. This in general reduces the amount of energy that can be exchanged between the mechanical deformations and the electric circuits connected to the piezoelectric patches and generally reduces the efficiency of the whole device (in particular, the values of the electromechanical coupling coefficients, introduced in Section 4, can be significantly reduced). Finally, hypothesis 6 is verified as long as the lateral boundary ( p) 0 of the piezoelectric patches coincides with a dielectric medium that contains no free charges.
Some practical system configurations, which verify the hypotheses of Section 3.1 and that lead to efficient systems, are introduced below and sketched in Figure 2.
• If s is made only of a dielectric material, one has just to verify that its permittivity ik is small and that the piezoelectric patches are not situated too close to one another. Some patches can be embedded in s (Figure 2(a)). This would be the case of a non-conductive composite structure. • If s is made only of a conducting material (this is the case for most of the metallic materials), one has to connect all the electrodes glued on s to s , to impose the same electric potential, and avoid unwanted capacitive effects of the glue film (Figure 2(b)). • If some parts of s are dielectric and other parts are conductors, if possible, one has to impose the same potential in all the conducting parts of s (Figure 2(c)).

MODAL EXPANSION AND COUPLING COEFFICIENTS
In this section, a reduced-order formulation of the discretized problem obtained in Section 3.3 is introduced, by expanding the mechanical displacement unknown vector onto the short-circuit eigenmode basis. The main motivation of choosing this particular basis is that it can be computed with a classical elastic mechanical problem, as it will be seen, whereas open-circuit modes depend also on the piezoelectric system properties. This basis could be enriched by other functions (like static modes, for example), which is out of the scope of the present article, that focuses on the introduction of the coupling coefficients.

Short-circuit normal modes
The system short-circuit normal modes are solutions of Equation (22a) with V = 0 and F = 0. The natural frequencies i and mode shapes U i (each U i is a column vectors of length N ) are the N eigensolutions of the following problem: which depends only on the mechanical properties of the system.

247
These modes verify the following orthogonality properties: where i j is the Kronecker delta and j have been normalized with respect to the mass matrix: One can note that the above equation imposes that the mode shapes U i units are kg

Modal expansion of the general problem
The displacement vector is sought as: By inserting the above equation in Equations (22a), (22b), multiplying the first obtained equation by U T i and using the orthogonality properties of Equations (25), the problem writes, for all i ∈{1, . . . , N }:q is the coupling coefficient of the ith mode and the pth piezoelectric patch, that is defined by: and F i (t) = U T i F(t) is the external forcing of the ith mode. One can get another writing of the above problem (28a), (28b) by substituting V ( p) as a function of Q ( p) (with Equation (28b)) in Equation (28a), so that the latter is replaced by: The units of the above introduced parameters/variables are: kg 1/2 m for q i , Nkg −1/2 for F i and kg 1/2 m s −2 V −1 (or in S.I. units kg −1/2 m −1 s A) for i . The initial finite element formulation introduced in Section 3.3 has been replaced by the modal formulation of Equations (28a), (28b) or Equations (30), (28b), whose unknowns are the N modal coordinates q i and the P voltage/charge pairs (V ( p) , Q ( p) ) associated with the P piezoelectric patches. Its major interest, and especially the choice of the short-circuit eigenmodes as the expansion basis, is that the above computations of the parameters necessitate only a single modal analysis of the elastic problem. This operation can thus be done by any standard finite elements code. • Calculate all matrices of the initial problem (20), that is to say M m , K m , K c . • Calculate theÑ natural frequencies i under interest in short circuit as well as the associated mode shapes U i , by solving problem (24). This is a classical eigenvalues/vectors computation, that involves only the elastic mass and stiffness matrices M m and K m . • Apply the matrix product of Equation (29) to obtain ( p) i .

Effective coupling factors
This section aims at showing that the above introduced coupling coefficients ( p) i are directly linked to the well-known effective electromechanical modal coupling factor (EEMCF). This parameter characterizes the energy exchanges between the mechanical structure and the piezoelectric patches. It is usually defined, for the system ith mode, by [31]: where i andˆ i are, respectively, the short-circuit and open-circuit ith system natural frequencies (with all piezoelectric patches, respectively, short circuited or in open circuit).
The set (30) truncated to the ith equation only-namely q j (t) ≡ 0 for all j = i-writes: In open circuit (Q ( p) = 0 for all p), the problem reads so that, by comparing Equations (32) and (33), The above expression gives an approximation of the system ith mode natural frequency in open circuit, denoted byˆ i , provided that the above introduced one-mode truncation is valid. As a consequence, by mixing Equations (31) and (34), the systems ith mode EEMCF is where k is the ith mode and the pth piezoelectric patch dimensionless coupling factor. By considering the definitions of obtains another expression for k eff,i stemming from Equation (35): If only one piezoelectric patch is connected to the structure, one can observe that k eff,i is directly linked to the only coupling coefficient i : where C is the patch capacitance. This relation can be generalized to the case of a system with several piezoelectric patches, by defining a new EEMCF, denoted by k ( p) eff,i , which is associated with the ith mode and the pth piezoelectric patch, by: In the above equation,ˆ so that k eff,i is found close to the coupling factor defined by Equation (36): The above equation shows that the coupling coefficient   (29)), but have the advantage of being dimensionless, which thus gives them a larger relevance. Moreover, they are directly linked to the EEMCF k  is linked to the pth patch only. As stated in Section 1, the coupling factors have a major importance in practice since they often constitute the main parameter of the system optimization, especially when considering shunts [24][25][26]. In this context, the approximate relation (41), in addition to giving a physical However, one must bear in mind that the above formula (41) can be used only if the one-mode truncation related to Equation (40) is valid, that is to say only if the natural frequencies are widely spaced. Otherwise, one can define and evaluate the problem EEMCFs with Equations (31) and (39), but they would have different values than the corresponding k ( p) i .

APPLICATION TO A PIEZOELECTRIC LAMINATED BEAM
This section aims at applying the general formulation introduced in Section 3 to a laminated beam composed of elastic and piezoelectric layers ( Figure 3). The main purpose of this section is to validate the general formulation of Section 3. The beam system has been chosen because it offers possible comparisons with experiments (see Section 6), even if its basic geometry can be thought as restrictive as compared with one of the systems considered in the general formulation of Section 3. Further validations on more complex three-dimensional systems are under study and will be published in the future.

Mechanical and electrical field assumptions
In the following, the explicit Voigt notations will be used rather than the tensor (indicial) notations. The beam is modeled using the classical laminate theory based on Euler-Bernoulli assumptions [38]. The mechanical displacement field can then be written as where u x and u z are the axial and transverse displacements; u is the axial displacement of the center line of the beam, w the transverse displacement, and the fiber rotation defined by From Equations (5) and the above-described kinematics, axial strain can be written as follows: where membrane strain e and curvature are defined by Concerning the electrostatic aspects, all hypotheses of Section 3.1 are kept so that the electric field is normal to the electrodes and uniform in the pth patch. The only non-vanishing component of the electric field is:

Piezoelectric constitutive equations
Following the hypotheses of Section 3.1, the piezoelectric layers of the laminated beam are poled in the thickness (z, 3)-direction with an electrical field applied parallel to this polarization. Such a configuration is characterized in particular by the electromechanical coupling between the axial strain ε 1 and the transverse electrical field E 3 . Even if neither plane-strain nor plane-stress assumption in the width direction are physically justified [27], we consider in this work, for the sake of simplicity, the classical plane-stress assumptions of beam theories ( 2 = 3 = 0). In this case, the 3D linear piezoelectric constitutive equations (3), (4) (using here a standard Voigt notation) are restricted to the mechanical axial components only: where k 1 and D k 3 are the axial stress and the transverse electrical displacement within the kth layer.c k 11 is the modified elastic stiffness at constant electric field,ē k 31 is the modified piezoelectric constant and¯ k 33 is the modified dielectric permittivity at constant strain. They are defined by: and Y k 1 the kth layer Young's modulus in the (1, 2) plane. The above equations are obtained by introducing the plane stress hypothesis in Equations (3), (4) together with the material transverse isotropy in the (1, 2) plane.

Finite elements discretization
Considering a particular axial part of the laminated beam (of length x + − x − and with a total of K layers including P piezoelectric patches), Equations (18) and (19) of the coupled electro-mechanical variational formulation can be rewritten in the following form: where k and k are the mass density and the domain occupied by the kth layer and the pth patch capacitance is (Equation (14)): with b and l ( p) are the width of the beam and the pth patch length. Notice that the beam is subjected to surface axial and transversal forces at the boundaries of each face sublayer (t k x , t k z ) and to body ones ( f k x , f k z ). This variational formulation is discretized using a two node finite element and the matrix equation (20) is obtained. All details that lead to matrices M m , K m , K e and F are written in the Appendix.
The case of coupling matrix K c is discussed here. We recall that this matrix, of size N × P, couples the N mechanical d.o.f. contained in U to the P electrical potential differences contained in V. Thus, K c can be written in the following form: c is a column vector that couples the mechanical d.o.f. to the potential difference of the pth patch only.
One single patch is generally discretized by more than one finite element such that K ( p) c is simply obtained by assembling the elementary vectors K e( p) c of size (6×1), whose explicit expression is given in the Appendix (Equation (A7)). Owing to their particular form, during the assembling operation, all components of K between the potential difference V ( p) and the mechanical dofs (axial translation u i and rotation i ) of the first and the last patch's nodes only. As expected, the piezoelectric effect only produces localized moments and forces at the patch edges [39,40].

NUMERICAL AND EXPERIMENTAL RESULTS
The purpose of this section is to apply the theoretical developments introduced before to the vibration reduction of a cantilevered beam, by using two piezoelectric patches and a resonant shunt.

System under study
A cantilever beam, already used in [24,26], is partially covered with two collocated piezoelectric elements, polarized in opposite directions ( Figure 4). The electrodes are connected in series to a passive electrical circuit composed of a resistor R and an inductor L, thus constituting a resonant shunt [14]. We use the '31' coupling of the piezoelectric elements that couple longitudinal (x-direction) deformations with the electric field in the transverse (z-direction), normal to the electrodes. Since the polarization directions are opposite in the elements, only the flexural deformations of the beam are coupled to the electric circuit. The beam material is elastic, homogeneous and isotropic, of density b and Young's modulus Y b . The piezoelectric material constants are its density p , Young's modulus Y p , coupling piezoelectric constantē 31 and dielectric constant¯ 33 . The system geometry is defined by the length/thickness of the beam and the piezoelectric patches, respectively, denoted l b / h b and l p / h p ; b denotes their common width (in the y-direction). Finally, the piezoelectric patches limits are defined by x − and x + , with x + − x − =l p . All numerical values of the studied system are gathered in Table I  to a constitutive law ( , D) as a function of (ε, E). The following relations have been used to obtain the parameters of In our case, k 31 = 0.372 andk 31 = 0.401.

Problem formulation
The finite element formulation introduced in Section 5 is here applied to the system sketched on Figure 4. The beam geometry is divided into three regions (Figure 4). Region I (for x ∈[0, x − ]) and region III (x ∈[x + ,l b ]) are composed of one layer of beam material of thickness h b . Region II corresponds to the piezoelectric patches (for x ∈[x − , x + ]); it is symmetric with respect to the transverse direction and is composed of three layers (two piezoelectric layers of thickness h p and one beam layer of thickness h b , in the middle). Regions I, II and III are discretized, respectively, with 1, 5 and 35 beam finite element, so that the total number of elements is N e = 41 and the total number of node is N n = 42. Because the finite element formulation uses three d.o.f. per node ((u, w, ), see Appendix), the number of mechanical dofs is N = 3N n = 126.
The above finite element discretization enables to evaluate the matrices of formulation (22a), (22b), which is reproduced here: In where C is the patches capacitance, identical for both patches, defined either by Equation (14) or by Equation (52). K c is a N ×2 matrix composed of two column vectors K 1 = K (1) c and K 2 = K Connecting the piezoelectric patches in series (see Figure 4) leads to impose equal charges in the piezoelectric patches: Q 1 = Q 2 = Q. By defining V = V 1 + V 2 the shunt terminal voltage (Figure 4), considering Equation (57) and summing up the two scalar equations of (56b), one obtains: By considering the RL electrical circuit of Figure 4, one can write: Finally, by gathering Equations (56a), (57) and (58) and eliminating V , V 1 and V 2 , one obtains the following matrix formulation of the electromechanical problem: with the following notations: where the stiffness matrix in open circuit-i.e. with Q = 0-is: (62) Figure 5 shows the location of non-zero elements of matrices K m , K c and K. An interesting feature is that the coupling matrix K c (and consequently the electromechanical coupling part of matrix K ) has non-zero elements only in the rows corresponding to the d.o.f. related to the boundary nodes of the piezoelectric patches. This stems from the electromechanical action of the piezoelectric patches on the mechanical structure, which is equivalent to concentrated moments localized on the boundaries of the patches, when the piezoelectric patches are perfectly bonded to the structure (see [39,40]).
Another practical way of connecting the shunt to the structure is to plug the piezoelectric patches in parallel and not in series. In this case, the systems behavior is similar to the series plugging: imposing V 1 = V 2 instead of Q 1 = Q 2 leads to the same general formulation of (60). However, the obtained expressions of matrices K andK m as functions of C, K 1 and K 2 are slightly different.

Modal formulation and coupling coefficients
The alternative modal reduced-order formulation of Section 4 can also be used. Since we are interested in the flexural motion of the beam, only the short-circuit flexural beam normal modes are considered here. They are denoted by ( i , U i ) and can be calculated with Equation (24). By considering Equations (29) and (57), the coupling coefficients Then, since U i are flexural modes, they have zero elements on the lines corresponding to d.o.f. u, which proves that i = i . Then, by eliminating V , V 1 and V 2 in Equations (28b) and (59), in the same way than in the previous section, the problem modal formulation truncated toÑ flexural modes becomes in this particular case: The present reduced-order formulation depends only on theÑ retained short-circuit modes ( i , U i ), theÑ associated coupling coefficients i and the patches capacitance C. To apply Section 4.3 results, we can define the following coupling factors: whereˆ i is the ith system open-circuit natural frequency, i.e with an unplugged shunt that imposes Q = 0 and thus zero electric charges on both piezoelectric patches. Then, Equation (35) shows that: which gives another way to evaluate i , provided that the one-mode approximation of Equation (32) is valid. Table II gives the frequencies and coupling coefficients numerical values obtained from the finite element model, for the first three flexural modes of the beam. The natural frequencies (in Hz) are defined as functions of the angular frequencies (in rad/s) by: The short-circuit natural frequencies ( f i , second column) and mode shapes (U i , Figure 6) have been calculated by numerically solving the eigenvalue problem of Equation (24). The open-circuit natural frequencies (f i , fourth column) and mode shapes (Û i , Figure 6) have been obtained by solving the same eigenvalue problem, but withK m (Equation (62)) as the stiffness matrix. The sixth column gives the coupling factors k i , evaluated by Equation (65), with C = 18.30 nF (Equation (57)) and i obtained with Equation (63). Finally, the seventh column gives the EEMCF k eff,i , obtained with Equation (65). The values obtained for k i and k eff,i are very close to each other, proving that the one-mode approximation that justifies Equation (66) is valid for the present system. In order to validate the present finite element formulation, experiments have been performed on a beam similar to the one sketched in Figure 4. All details are written in [24,25] and only the main features are recalled here. The beam is clamped with a vice. An non-contact electromagnetic driving system is used, composed of a small magnet glued on the structure with bees-wax, subjected to the magnetic field created by a coil, fed by a sine electrical signal. This exciter is fully described in [43]. The force acting on the beam is estimated by measuring the current intensity in the coil, proportional to the force. The beam motion is obtained with a laser Doppler vibrometer, that measures the velocity of the beam in one point. Two modal analyses are performed, the first one with the piezoelectric patches short circuited with a wire and the second one with unplugged patches. The first four natural frequencies (corresponding to three flexural modes and one torsional mode) are measured, in the above-described short-circuit ( f i ) and open-circuit conditions (f i ). The corresponding values are gathered in Table II, together with an estimation of k eff,i with Equation (65). The good agreement between the FE model and the experiments has been obtained by adding a concentrated mass of 4.2 g at the closest node to the magnet location (at the tip), that slightly lowers the natural frequency as compared with a naked beam and that equally modifies the mode shapes ( Figure 6), thus, validating the present formulation.
One can remark that the torsional mode has a low coupling coefficient as compared with the ones of the flexural modes. It would be zero in theory, since the extensional motion of the piezoelectric patches is not activated when the beam undergoes torsion. Moreover, since the FE model does not include torsion d.o.f., no FE values are available in Table II. Another remark is that for the first mode, the measured value of the effective coupling coefficient k eff,i is higher that the finite-elements one. One would expect the opposite, for example, because of non-perfect bonding, as found for the second and third flexural modes. A difference between model and experiments is the wrapped electrode, used to easily solder the electric wire to the lower electrode on the piezoelectric patch upper face. It locally modifies the electric field intensity and direction near the edge of the piezoelectric patch. This may explain the observed increase of coupling coefficient in the experiments, even if the exact effect of the wrapped electrode is difficult to characterize without a three-dimensional electromechanical simulation. This is beyond the scope of the present article.

Vibratory response
The vibratory response of the system subjected to a harmonic forcing can be evaluated with the formulation of Equation (60), when F = F 0 cos t. Defining the complex amplitude X by X(t) = (Xe j t ) leads to the frequency formulation of Equation (60):  with the impedance matrix defined by: The system frequency response function (FRF) at the ith d.o.f. to a forcing at the jth d.o.f. is then, for each forcing frequency : where Z −1 i j ( ) is the ith line and the jth column term of Z( ) inverse matrix. Figure 7 shows the driving point FRF modulus (|H ii ( )|) at the tip of the beam (i.e. the response at the tip for an excitation at the same point), from the present finite element formulation as well as from the experiments, with the setup described in Section 6.3. Again, an excellent agreement is obtained between simulations and experiments. The main difference lies in the structural damping that is not taken into account in the FE model, which leads to infinite responses, in short circuit, at the beam resonances (Figure 7(a)). Moreover, as announced in the previous section, the torsion mode, of frequency around 853 Hz, is not predicted by the finite element model.
Finally, the effect of the shunt is shown in Figure 7(b), with the values of R = 7900 and L = 21.8 H tuned to lower the second flexural resonance, around 336 Hz. In the experiments, an attenuation of 25 dB is obtained as compared with the short-circuit resonance. Again an excellent agreement is obtained between the simulated and experimental damped FRFs. As already said, no agreement is possible around the resonances for the short-circuit FRFs because of the absence of structural damping in the FE model.

Reduced-order model
The system vibratory response can also be obtained with the modal ROM defined by Equations (64), truncated toÑ modal coordinates. The problem to solve takes the form of Equation (60)

CONCLUSION
In this article, a finite element formulation of the coupled electromechanical problem of an elastic structure equipped with piezoelectric patches has been introduced. Its originality is that provided a set of non restrictive-assumptions, the system electrical state is fully described by very few global discrete unknowns: only a couple of variables per piezoelectric patches, namely (1) the electric charge contained in the electrodes and (2) the voltage between the electrodes.
It has been shown that this formulation has several important features. First of all, the electrical part of the problem is fully discretized at the weak formulation step, by introducing the above cited voltage/charge variables, without any restriction on the mechanical part of the problem. As a consequence, any standard (elastic only) finite element formulation can be easily modified to 261 include the piezoelectric patches and thus the effect of an external electrical action. A second advantage of this formulation is that since global electrical variables are used, realistic electrical boundary conditions are naturally introduced. First, the equipotentiality in any of the patches electrodes is exactly satisfied when introducing the potential difference variable. Second, using the global charge contained in the electrodes as the second electrical variable is realistic since plugging an external electrical circuit to the electrodes of the patches imposes only the global charge contained in the electrodes and not a local charge surface density.
The second advantage of using the global charge/voltage variables is that they are intrinsically adapted to include any external electrical circuit into the electromechanical problem. Thus, the present formulation can be used to efficiently simulate the system vibratory behavior in any classical situation and in particular in the field of vibration control, where some of the patches are used as sensors and the other as actuators. It is also especially adapted to simulate shunt systems, where a passive (or semi-passive) electrical network is connected to the patches and where both piezoelectric effects are used at the same time. In this case, one has only to add to the finite element formulation the equation(s) of the electric network, that link the patches charge and voltage variables. This last application has been thoroughly illustrated in the article in the case of a cantilever beam whose vibrations are reduced by means of a resonant shunt. A finite element formulation of the problem has been described as well as experiment and an excellent agreement was obtained.
Finally, a reduced-order model (ROM) of the problem has been introduced by expanding the solution onto the system's normal mode basis with all patches short circuited. It shows that the classical efficient electromechanical coupling factors (EEMCF) naturally appear in the ROM as the main parameters that master the electromechanical coupling. A simple analytical formula is obtained, so that the EEMCF depend only on the matrices (mass, stiffness and coupling) of the coupled finite element problem as well as on the short-circuit modes. Since those modes are solution of a modal analysis of the elastic problem, any standard elastic finite element code can be used to obtain them. This is an advantage, since the EEMCF is often the main parameter that influences the efficiency of vibration control systems, especially in the case of shunts. To illustrate this result, the EEMCF for the first three modes of the cantilever beam, obtained from the finite element study as well as from the experiments are compared, showing an excellent agreement.
The natural extension of this work is to apply the finite-elements formulation introduced in the present article to more complex three-dimensional structures. In particular, the simulation and the optimization of the reduction of turbomachinery blades vibration is under study [44]. Moreover, the ROM introduced in the present work will also be tested to simulate the behavior of a new category of shunts, often known as 'switch', where the electrical circuit impedance is switched periodically between two values [17][18][19]. Since the process is intrinsically non-linear, it specifically requires the use of direct time integration, over long time periods, to attain the steady state. where U e = (u 1 w 1 1 u 2 w 2 2 ) T u 1 , w 1 and 1 (respectively, u 2 , w 2 and 2 ) corresponding to the first (respectively, second) node of the beam element.
In equation (A1), the interpolation vectors are defined as follows: with, for all x ∈[0, L e ] where L e is the length of the considered finite element. Note that Equations (43) where (·) = *(·)/*x. The various terms appearing in the variational formulation of Equations (50) and (51) are now successively discussed.
• The kinetic energy variation is: where the zero-, first-, and second-order inertia moments are defined by with the zero, first and second moment of area of each sublayer, respectively, defined by • Work done by the external mechanical forces where k is the cross-section area of the layer k and where the boundary and distributed (i) normal resultant, (ii) shear resultant and (iii) bending moment are, respectively, defined by: The elementary vector of generalized mechanical forces is then: whose explicit expression is: • The piezoelectric contributions of the internal energy variation, related to the inverse effect, are: where G The same matrix is obtained for the direct effect.