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

Olivier Thomas, Jean-François Deü, Julien Ducarne

► To cite this version:

Olivier Thomas, Jean-François Deü, Julien Ducarne. Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients. International Journal for Numerical Methods in Engineering, 2009, 80 (2), pp.235-268. 10.1002/nme.2632 . hal-01572491

HAL Id: hal-01572491

https://hal.science/hal-01572491

Submitted on 7 Aug 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients

O. Thomas, J. F. Deü and J. Ducarne

*Structural Mechanics and Coupled Systems Laboratory, Conservatoire National des Arts et Métiers, 2 rue Conté, 75003 Paris, France*

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.

KEY WORDS: piezoelectric materials; vibration reduction; electric boundary conditions; shunt; electromechanical coupling
# 1. 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–10]. For more references, one can refer to the following review articles [11–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–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
operation requires a maximization of the effective electromechanical coupling factors (EEMCF). In the case of passive resistive or resonant shunt techniques, whose purpose is the structural vibration reduction, the EEMCFs of the targeted vibration modes have to be maximized, as shown in [14, 22, 23]. Moreover, it is proved in [24, 25] that the EEMCF is the only free parameter of the optimization, the optimal value of the electric circuit parameters being known as functions of the EEMCF and the system structural characteristics, and evaluated in a second step. The same results are also obtained in the case of switch techniques [26].

The first goal of this article is to write a discretized formulation of the dynamics of an elastic structure of arbitrary geometry with piezoelectric patches, the latter having a thin shape and being polarized in their transverse direction. The main advantage of our formulation is that provided a few non-restrictive electrical hypotheses, the system electrical state is defined by very few discrete global unknowns: on the one hand, the free electric charges contained in the electrodes of the piezoelectric patches and on the other hand the associated terminal potential differences. This formulation enables to take into account both the direct and indirect piezoelectric effects, separately or at the same time, while there is no need to compute the whole electric potential field. Concerning the system mechanical state, no restrictions are formulated so that it appears in a standard finite element form in the formulation, with the displacement d.o.fs at the mesh nodes and the associated generalized forces.

The present problem formulation can be used to simulate the dynamics of the system in any standard situation of a smart structure connected to an electrical network: a part of the piezoelectric patches can be sensors and the remaining part actuators. However, this study especially addresses the simulation of the dynamics of the system coupled to shunted electric circuits, since each piezoelectric patch appears in the formulation with a pair of global intuitive unknowns, its potential difference and its electric charge, that can be introduced without restrictions in the evolution equation of the electrical circuit to which it is plugged. This is an originality of this work, since most other studies are restricted to piezoelectric elements used either as sensor or as actuators (see above cited references). Another advantage of our formulation is that realistic electric boundary conditions on the piezoelectric patches—i.e. on the one hand equipotentiality in the electrodes and on the other hand a global free electric charge in the electrodes—are naturally imposed and included in our formulation. This is another originality since classical finite element formulations of the literature often use local electric variables—the electric potential and displacement—that are not naturally adapted to the above cited global boundary conditions. The equipotentiality in the electrodes is often enforced by assigning a single d.o.f. for voltage on all their nodes [27–29].

The second goal of this work is to introduce an efficient manner of evaluating the effective electromechanical coupling factors (EEMCF) of the system. As stated above, most of the literature works on the subject show that the EEMCFs are the major parameters that govern, on the one hand, all the performances of the vibration dampers made with shunted piezoelectric elements, and on the other hand, all the optimizations procedures of those devices [24–26]. Thus, evaluating the EEMCF for a given system is crucial when designing its coupling with an electrical circuit. To do so, a given EEMCF is associated in our study with one system vibration eigenmode and one piezoelectric patch. The general discretized finite element formulation introduced above is expanded onto the vibration eigenmodes of the system with all its piezoelectric elements short-circuited. A modal discretized formulation is obtained, where all the EEMCFs appear naturally. The latter are known with analytical expressions as functions of the matrices of the finite element discretization and the short-circuit eigenmodes. This result can be useful in an engineering context, since the exhibited EEMCFs numerical computation involves only the assembling of the finite
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.

## 2. 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.

### 2.1. Electro-mechanical local equations

We consider a piezoelectric structure occupying a domain $\Omega_p$ at the equilibrium. The structure is subjected to a prescribed displacement $u_i^d$ on a part $\Gamma_u$ and to a prescribed surface force density $t_i^d$ on the complementary part $\Gamma_t$ of its external boundary. The electric boundary conditions are defined by a prescribed electric potential $\psi^d$ on $\Gamma_\psi$ and a surface density of free electric charges $q^d$ on the remaining part $\Gamma_q$. Thus, the total structure boundary, denoted $\partial\Omega_p$, is such that $\partial\Omega_p = \Gamma_u \cup \Gamma_t = \Gamma_\psi \cup \Gamma_q$ with $\Gamma_u \cap \Gamma_t = \Gamma_\psi \cap \Gamma_q = \emptyset$. In addition, $\Omega_p$ is subjected to prescribed body forces $f_i^d$.

The linearized deformation tensor and the stress tensor are denoted, respectively, by $\varepsilon_{ij}$ and $\sigma_{ij}$. Moreover, $D_i$ denotes the electric displacement vector components and $E_i$ the electric field vector components. $\rho$ is the mass density, $n_i$ is the unit normal external to $\Omega_p$ and $t$ is the time.

The local equations of the electro-mechanical coupled problem are:

$$ \sigma_{ij,j} + f_i^d = \rho \frac{\partial^2 u_i}{\partial t^2} \quad \text{in } \Omega_p \quad (1a) $$

$$ \sigma_{ij} n_j = t_i^d \quad \text{on } \Gamma_t \quad (1b) $$

$$ u_i = u_i^d \quad \text{on } \Gamma_u \quad (1c) $$
and

$$D_{i,i} = 0 \quad \text{in } \Omega_p \qquad (2a)$$

$$D_i n_i = -q^d \quad \text{on } \Gamma_q \qquad (2b)$$

$$\psi = \psi^d \quad \text{on } \Gamma_{\psi} \qquad (2c)$$

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 $\Omega_p$ can be neglected with respect to its value inside $\Omega_p$; otherwise, the jump of $D_i n_i$ through $\partial\Omega_p$ would have appeared in (2b). Moreover, no volume free electric charges are prescribed in $\Omega_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 $\sigma_{ij}$ and electric displacement $D_i$ are related to the linear strain tensor $\epsilon_{kl}$ and electric field $E_k$ through the converse and direct linear piezoelectric constitutive equations [31, 32]:

$$\sigma_{ij} = c_{ijkl}\epsilon_{kl} - e_{kij}E_k \qquad (3)$$

$$D_i = e_{ikl}\epsilon_{kl} + \epsilon_{ik}E_k \qquad (4)$$

where $c_{ijkl}$ denote the elastic moduli at constant electric field, $e_{kij}$ the piezoelectric constants and $\epsilon_{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 $\epsilon_{kl}$ and the displacement $u_k$, and between the electric field $E_k$ and the electric potential $\psi$:

$$\epsilon_{kl}(u) = \frac{1}{2}(u_{k,l} + u_{l,k}) \qquad (5)$$

$$E_k(\psi) = -\psi_{,k} \qquad (6)$$

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 $\psi$ in Equation (6) [30, 32].

## 2.2. Variational formulation in terms of $(u_i, \psi)$

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 $\psi$. 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 $\Omega_p$, $C_u^* = \{u_i \in C_u | u_i = 0 \text{ on } \Gamma_u\}$ and $C_\psi$ of sufficiently regular functions $\psi$ in $\Omega_p$.
Multiplying Equations (1a) and (2a), respectively, by any time-independent test functions $\delta u_i \in \mathcal{C}_u^*$ and $\delta \psi \in \mathcal{C}_\psi$, integrating over $\Omega_p$, applying Green's formula, taking into account Equations (1b) and (2b) as well as the constitutive laws (3) and (4) leads to:

$$ \int_{\Omega_p} c_{ijkl} \varepsilon_{kl}(u) \varepsilon_{ij}(\delta u) \, d\Omega - \int_{\Omega_p} e_{kij} E_k(\psi) \varepsilon_{ij}(\delta u) \, d\Omega + \int_{\Omega_p} \rho \frac{\partial^2 u_i}{\partial t^2} \delta u_i \, d\Omega = \int_{\Omega_p} f_i^d \delta u_i \, d\Omega + \int_{\Gamma_t} t_i^d \delta u_i \, dS \quad (7) $$

and

$$ \int_{\Omega_p} e_{ikl} \varepsilon_{kl}(u) E_i(\delta \psi) \, d\Omega + \int_{\Omega_p} \varepsilon_{ik} E_k(\psi) E_i(\delta \psi) \, d\Omega = \int_{\Gamma_q} q^d \delta \psi \, dS + \int_{\Gamma_\psi} q^r \delta \psi \, dS \quad (8) $$

where $q^r = -D_i n_i$ is the reaction free charge density that appears on $\Gamma_\psi$.

The variational formulation of the coupled electro-mechanical problem, equivalent to the strong form (Equations (1a)–(6)), consists in finding $u_i \in \mathcal{C}_u$ such that $u_i = u_i^d$ on $\Gamma_u$, $\psi \in \mathcal{C}_\psi$ such that $\psi = \psi^d$ on $\Gamma_\psi$ and the reaction electric free charge density $q^r$ on $\Gamma_\psi$, satisfying Equations (7) and (8) for all $\delta u_i \in \mathcal{C}_u^*$ and for all $\delta \psi \in \mathcal{C}_\psi$, with appropriate initial conditions.

### 3. 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 $\psi$ 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.

#### 3.1. General problem and hypotheses

We consider an elastic structure equipped with $P \in \mathbb{N}$ piezoelectric patches (Figure 1). The domain occupied by the elastic part of the system is denoted by $\Omega_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 $p$th patch, $p \in \{1, \dots, P\}$, occupies a domain $\Omega^{(p)}$ such that $(\Omega_s, \Omega^{(1)}, \dots, \Omega^{(P)})$ is a partition of the whole domain $\Omega$. The $p$th patch upper and lower electrode surfaces are denoted, respectively, by $\Gamma_+^{(p)}$ and $\Gamma_-^{(p)}$ and its lateral surface is $\Gamma_0^{(p)}$. In the remaining part of the article, superscript $(p)$ will denote a quantity related to the $p$th 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_{ijk}$ vanish in $\Omega_s$. However, we can mention that the dielectric constants $\epsilon_{ik}$ do not vanish in $\Omega_s$, in general [30].
<figure><img src="image_2.png" /><figcaption>Figure 1. Left: an elastic structure with two piezoelectric patches. Right: the pth piezoelectric patch submitted to the potential difference $V^{(p)}$, with its electrodes containing the free electric charges $Q^{(p)}$ and $-Q^{(p)}$ on their surface in contact with the patch.</figcaption></figure>

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 $\Gamma_0^{(p)}$ 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 \in \{1, \dots, P\}$:

$$E_k = -\frac{\psi_+^{(p)} - \psi_-^{(p)}}{h^{(p)}} n_k = -\frac{V^{(p)}}{h^{(p)}} n_k \quad \text{in } \Omega^{(p)} \qquad (9)$$

where $V^{(p)} = \psi_+^{(p)} - \psi_-^{(p)}$ is the potential difference between $\Gamma_+^{(p)}$ and $\Gamma_-^{(p)}$ 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 $\Omega_s$, as compared with its value in any of the piezoelectric patches $\Omega^{(p)}$. The justification of this hypothesis will be discussed in Section 3.4. However, the electric field in $\Omega_s$ does not vanish, in general.
### 3.2. Variational formulation in terms of $(u_i, V^{(p)})$

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 $\Omega$, namely the elastic part $\Omega_s$ and the piezoelectric patches $\Omega^{(p)}$, for $p \in \{1, \dots, P\}$.

#### 3.2.1. Mechanical equation.
We first apply Equation (7) successively to each subdomain of $\Omega$. 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:

$$ \int_{\Omega} = \int_{\Omega_s} + \sum_{p=1}^{P} \int_{\Omega^{(p)}} \qquad (10) $$

Consequently to the equilibrium of forces and the equality of displacement through each interface between two subdomains of $\Omega$, all the surface integrals that are related to the subdomain boundaries inside $\Omega$ cancel by pair. Then, we redefine $\Gamma_u$ and $\Gamma_t$ to be a partition of the boundary $\partial\Omega$ of $\Omega$ where the displacements and the forces are prescribed, respectively. We also redefine $\mathcal{C}_u$ as the space of sufficiently regular functions $u_i$ defined in the whole $\Omega$ and $\mathcal{C}_u^* = \{u_i \in \mathcal{C}_u | u_i = 0 \text{ on } \Gamma_u\}$. Then, the piezoelectric integral term in Equation (7), with $\Omega_p$ replaced by $\Omega$, can be rewritten as

$$ \int_{\Omega} e_{kij} E_k(\psi) \varepsilon_{ij}(\delta u) \, d\Omega = \sum_{p=1}^{P} -\frac{V^{(p)}}{h^{(p)}} \int_{\Omega^{(p)}} e_{kij} n_k \varepsilon_{ij}(\delta u) \, d\Omega \qquad (11) $$

where $V^{(p)}$ is the $p$th 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_{ijk} = 0$ in $\Omega_s$ so that the corresponding integral vanishes. Finally, assumption 8 (together with Equation (9)) leads to the result.

#### 3.2.2. Electrical equation.
In a second step, we apply Equation (8) to each piezoelectric patch. We chose here to prescribe the electrical potential $\psi$ on all the piezoelectric patches electrodes $\Gamma_+^{(p)}$ and $\Gamma_-^{(p)}$. Moreover, following hypothesis 6, a zero surface free charge density $q^d = 0$ is prescribed on the lateral boundaries $\Gamma_0^{(p)}$ of the piezoelectric patches. Consequently, the electrical boundary conditions on the $p$th piezoelectric patch are defined by $\Gamma_\psi$ that coincides with $\Gamma_+^{(p)} \cup \Gamma_-^{(p)}$ and $\Gamma_q$ that coincides with $\Gamma_0^{(p)}$.

The piezoelectric term in Equation (8) is processed in the same manner as in Section 3.2.1, so that one obtains, with $\Omega_p$ replaced by $\Omega^{(p)}$:

$$ \int_{\Omega^{(p)}} e_{ikl} \varepsilon_{kl}(u) E_i(\delta\psi) \, d\Omega = -\frac{\delta V^{(p)}}{h^{(p)}} \int_{\Omega^{(p)}} e_{ikl} \varepsilon_{kl}(u) n_i \, d\Omega \qquad (12) $$

where $\delta V^{(p)} = \delta\psi_+^{(p)} - \delta\psi_-^{(p)}$ with $\delta\psi_+^{(p)}$ and $\delta\psi_-^{(p)}$ defined as the uniform values of $\delta\psi$ on the $p$th patch upper and lower electrodes.
The dielectric term in Equation (8) can be rewritten as, with again $\Omega_p$ replaced by $\Omega^{(p)}$:

$$ \int_{\Omega^{(p)}} \epsilon_{ik} E_k(\psi) E_i(\delta\psi) d\Omega = \frac{V^{(p)} \delta V^{(p)}}{(h^{(p)})^2} \int_{\Omega^{(p)}} \epsilon_{ik} n_i n_k d\Omega \quad (13) $$

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:

$$ \int_{\Omega^{(p)}} \epsilon_{ik} E_k(\psi) E_i(\delta\psi) d\Omega = C^{(p)} V^{(p)} \delta V^{(p)} \quad \text{with} \quad C^{(p)} = \frac{\epsilon_{33} S^{(p)}}{h^{(p)}} \quad (14) $$

the capacitance of the pth piezoelectric patch and $\epsilon_{33} = \epsilon_{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) are now addressed. First, as explained above, $q^d = 0$. Then, considering that $\delta\psi$ is uniform on the upper and lower electrode with values $\delta\psi_+^{(p)}$ and $\delta\psi_-^{(p)}$, the remaining integral corresponds to the sums of all the free charges contained in the electrodes, named $Q_+^{(p)}$ and $Q_-^{(p)}$. Moreover, we apply the integral form of Gauss's law to a closed surface $\Sigma$ that surrounds the pth patch. It stands that, at the electric equilibrium,

$$ \oint_{\Sigma} D_i n_i dS = Q \quad (15) $$

where $Q$ is the free electric charge contained inside $\Sigma$. In our particular case, $\Sigma$ lies in $\Omega_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 $\Sigma$, which means that the left-hand side of Equation (15) vanishes (the electric displacement in the electric cables that eventually cross $\Sigma$ is zero because they are conducting media). Consequently, the whole charge $Q$ contained inside $\Sigma$ 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 $\Gamma_+^{(p)}$ of the pth piezoelectric patch, one obtains:

$$ Q_+^{(p)} = -Q_-^{(p)} = Q^{(p)} \quad (16) $$

Finally, the remaining term of the right-hand side of Equation (8) writes, by substituting $\Gamma_\psi$ for $\Gamma_+^{(p)} \cup \Gamma_-^{(p)}$:

$$ \int_{\Gamma_{+}^{(p)} \cup \Gamma_{-}^{(p)}} q^r \delta\psi \, dS = Q_{+}^{(p)} \delta\psi_{+}^{(p)} + Q_{-}^{(p)} \delta\psi_{-}^{(p)} = Q^{(p)} (\delta\psi_{+}^{(p)} - \delta\psi_{-}^{(p)}) = Q^{(p)} \delta V^{(p)} \quad (17) $$

### 3.2.3. 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 \in \mathcal{C}_u$ such that $u_i = u_i^d$ on $\Gamma_u$ and the P electric charges $Q^{(p)}$ ($p \in \{1, \dots, P\}$), such that the associated potential differences $V^{(p)}$ are
prescribed, satisfying:

$$ \int_{\Omega} c_{ijkl} \varepsilon_{kl}(u) \varepsilon_{ij}(\delta u) \, d\Omega + \sum_{p=1}^{P} \frac{V^{(p)}}{h^{(p)}} \int_{\Omega^{(p)}} e_{kij} n_k \varepsilon_{ij}(\delta u) \, d\Omega \\ + \int_{\Omega} \rho \frac{\partial^2 u_i}{\partial t^2} \, \delta u_i \, d\Omega = \int_{\Omega} f_i^d \, \delta u_i \, d\Omega + \int_{\Gamma_i} t_i^d \, \delta u_i \, dS \quad \forall \delta u_i \in \mathcal{C}_u^* \quad (18) $$

and

$$ -\sum_{p=1}^{P} \frac{\delta V^{(p)}}{h^{(p)}} \int_{\Omega^{(p)}} e_{ikl} \varepsilon_{kl}(u) n_i \, d\Omega + \sum_{p=1}^{P} \delta V^{(p)} C^{(p)} V^{(p)} = \sum_{p=1}^{P} \delta V^{(p)} Q^{(p)} \quad \forall \delta V^{(p)} \in \mathbb{R} \quad (19) $$

with appropriate initial conditions.

### 3.3. 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:

$$ \begin{pmatrix} \mathbf{M}_m & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \ddot{\mathbf{U}} \\ \ddot{\mathbf{V}} \end{pmatrix} + \begin{pmatrix} \mathbf{K}_m & \mathbf{K}_c \\ -\mathbf{K}_c^\mathrm{T} & \mathbf{K}_e \end{pmatrix} \begin{pmatrix} \mathbf{U} \\ \mathbf{V} \end{pmatrix} = \begin{pmatrix} \mathbf{F} \\ \mathbf{Q} \end{pmatrix} \quad (20) $$

In the above equation, $\mathbf{M}_m$ and $\mathbf{K}_m$ are the mechanical mass and stiffness matrices, of size $N \times N$, if $N$ is the number of mechanical unknowns. $\mathbf{K}_e = \mathrm{diag}(C^{(1)} \ C^{(2)} \ ... \ C^{(P)})$ is a diagonal matrix filled with the $P$ capacitances of the piezoelectric patches. $\mathbf{K}_c$ is the electromechanical coupling matrix, of size $N \times P$. $\mathbf{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 $\mathbf{V}$ nor $\mathbf{Q}$ are prescribed by the electrical network but the latter imposes only a relation between them, that can be written as $\mathbf{f}(\ddot{\mathbf{V}}, \dot{\mathbf{V}}, \mathbf{V}, \ddot{\mathbf{Q}}, \dot{\mathbf{Q}}, \mathbf{Q}) = \mathbf{0}$. As an example, this additional relation writes, in the case of a resonant shunt connected to the $p$th patch and composed of a resistance $R$ and an inductance $L$ [24, 25]:

$$ V^{(p)} - L \ddot{\mathbf{Q}}^{(p)} - R \dot{\mathbf{Q}}^{(p)} = 0 \quad (21) $$

Back to the general case, the problem consists in finding $U$, $V$ and $Q$ so that:

$$ \mathbf{M}_m \ddot{\mathbf{U}} + \mathbf{K}_m \mathbf{U} + \mathbf{K}_c \mathbf{V} = \mathbf{F} \quad (22a) $$

$$ \mathbf{K}_e \mathbf{V} - \mathbf{Q} - \mathbf{K}_c^\mathrm{T} \mathbf{U} = \mathbf{0} \quad (22b) $$

$$ \mathbf{f}(\ddot{\mathbf{V}}, \dot{\mathbf{V}}, \mathbf{V}, \ddot{\mathbf{Q}}, \dot{\mathbf{Q}}, \mathbf{Q}) = \mathbf{0} \quad (22c) $$

It is sometimes convenient to rewrite Equation (22a) in a $(\mathbf{U}, \mathbf{Q})$ form, using Equation (22b). One thus obtains:

$$ \mathbf{M}_m \ddot{\mathbf{U}} + (\mathbf{K}_m + \mathbf{K}_c \mathbf{K}_e^{-1} \mathbf{K}_c^\mathrm{T}) \mathbf{U} + \mathbf{K}_c \mathbf{K}_e^{-1} \mathbf{Q} = \mathbf{F} \quad (23) $$
In this latter formulation, the electromechanical coupling effect on the mechanical structure appears with an added-stiffness term, that depends on the coupling matrix $\mathbf{K}_e$ and the patches capacitances gathered in $\mathbf{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.

### 3.4. 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 $\Omega_s$ as compared with its value in the piezoelectric patches. We first recall that thanks to hypothesis 1, $D_i = \epsilon_{ik} E_k$ in $\Omega_s$ so that $D_i$ depends only on the medium permittivity $\epsilon_{ik}$ and the electric field $E_k$. Most of the engineering materials, conductors as well as dielectrics, have their permittivity constants $\epsilon_{ik}$ of the order of magnitude of $3\epsilon_0$, a value much smaller than those of common piezoelectric materials, of the order of $1000\epsilon_0$ [30] ($\epsilon_0$ is the permittivity of free space). Consequently, hypothesis 9 is verified if the value of the electric field $E_k$ in $\Omega_s$ is of the same order of magnitude than in the piezoelectric patches. Since $E_k$ is the gradient of the electric potential in $\Omega$, 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 $\Omega_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 $\Omega_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 $\Omega_s$ is obtained if one wants to glue a given piezoelectric patch onto a conducting part of $\Omega_s$ without imposing the same potential in $\Omega_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 $\mathbf{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 $\Omega_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).
<figure><img src="image_3.png" /><figcaption>Figure 2. Examples of system configurations that verify the hypotheses of Section 3.1. (a) The elastic structure $\Omega_s$ is made of a dielectric material; (b) the elastic structure $\Omega_s$ is made of a conducting material; and (c) parts of the elastic structure $\Omega_s$ are dielectric, others are conductors.</figcaption></figure>

Finally, hypothesis 6 is verified as long as the lateral boundary $\Gamma_0^{(p)}$ 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 $\Omega_s$ is made only of a dielectric material, one has just to verify that its permittivity $\epsilon_{ik}$ is small and that the piezoelectric patches are not situated too close to one another. Some patches can be embedded in $\Omega_s$ (Figure 2(a)). This would be the case of a non-conductive composite structure.

* If $\Omega_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 $\Omega_s$ to $\Omega_s$, to impose the same electric potential, and avoid unwanted capacitive effects of the glue film (Figure 2(b)).

* If some parts of $\Omega_s$ are dielectric and other parts are conductors, if possible, one has to impose the same potential in all the conducting parts of $\Omega_s$ (Figure 2(c)).

## 4. 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.

### 4.1. Short-circuit normal modes

The system short-circuit normal modes are solutions of Equation (22a) with $V=0$ and $F=0$. The natural frequencies $\omega_i$ and mode shapes $\Phi_i$ (each $\Phi_i$ is a column vectors of length $N$) are the eigensolutions of the following problem:

$$ K_m \Phi - \omega^2 M_m \Phi = 0 \qquad (24) $$

which depends only on the mechanical properties of the system.
These modes verify the following orthogonality properties:

$$ \forall (i, j), \quad \mathbf{\Phi}_i^\mathrm{T} \mathbf{M}_m \mathbf{\Phi}_j = \delta_{ij}, \quad \mathbf{\Phi}_i^\mathrm{T} \mathbf{K}_m \mathbf{\Phi}_j = \omega_i^2 \delta_{ij} \qquad (25) $$

where $\delta_{ij}$ is the Kronecker delta and $\Phi_j$ have been normalized with respect to the mass matrix:

$$ \forall i, \quad \mathbf{\Phi}_i^\mathrm{T} \mathbf{M}_m \mathbf{\Phi}_i = 1 \qquad (26) $$

One can note that the above equation imposes that the mode shapes $\mathbf{\Phi}_i$ units are $kg^{-1/2}$ for the elements linked to d.o.f. $u$ and $w$ and $kg^{-1/2}m^{-1}$ for the elements linked to d.o.f. $\theta$ (see Equation (A3) for the units of $\mathbf{M}_m$).

## 4.2. Modal expansion of the general problem

The displacement vector is sought as:

$$ \mathbf{U}(t) = \sum_{i=1}^{N} \mathbf{\Phi}_i q_i(t) \qquad (27) $$

By inserting the above equation in Equations (22a), (22b), multiplying the first obtained equation by $\mathbf{\Phi}_i^\mathrm{T}$ and using the orthogonality properties of Equations (25), the problem writes, for all $i \in \{1, \dots, N\}$:

$$ \ddot{q}_i + \omega_i^2 q_i + \sum_{p=1}^{P} \chi_i^{(p)} V^{(p)} = F_i, \quad \forall i \in \{1 \dots N\} \qquad (28a) $$

$$ C^{(p)} V^{(p)} - Q^{(p)} - \sum_{i=1}^{N} \chi_i^{(p)} q_i = 0, \quad \forall p \in \{1 \dots P\} \qquad (28b) $$

where $\chi_i^{(p)}$ is the coupling coefficient of the $i$th mode and the $p$th piezoelectric patch, that is defined by:

$$ (\chi_i^{(1)} \chi_i^{(2)} \dots \chi_i^{(P)}) = \mathbf{\Phi}_i^\mathrm{T} \mathbf{K}_c, \quad \forall i \in \{1 \dots N\} \qquad (29) $$

and $F_i(t) = \mathbf{\Phi}_i^\mathrm{T} \mathbf{F}(t)$ is the external forcing of the $i$th 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:

$$ \ddot{q}_i + \omega_i^2 q_i + \sum_{p=1}^{P} \sum_{k=1}^{N} \frac{\chi_i^{(p)} \chi_k^{(p)}}{C^{(p)}} q_k + \sum_{p=1}^{P} \frac{\chi_i^{(p)}}{C^{(p)}} Q^{(p)} = F_i, \quad \forall i \in \{1 \dots N\} \qquad (30) $$

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}sA$) for $\chi_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.
To recall, one has to:

* Calculate all matrices of the initial problem (20), that is to say $M_m, K_m, K_c$.

* Calculate the $\tilde{N}$ natural frequencies $\omega_i$ under interest in short circuit as well as the associated mode shapes $\Phi_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 $\chi_i^{(p)}$.

### 4.3. Effective coupling factors

This section aims at showing that the above introduced coupling coefficients $\chi_i^{(p)}$ 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 $i$th mode, by [31]:

$$k_{\text{eff},i}^2 = \frac{\hat{\omega}_i^2 - \omega_i^2}{\omega_i^2} \qquad (31)$$

where $\omega_i$ and $\hat{\omega}_i$ are, respectively, the short-circuit and open-circuit $i$th system natural frequencies (with all piezoelectric patches, respectively, short circuited or in open circuit).

The set (30) truncated to the $i$th equation only—namely $q_j(t) \equiv 0$ for all $j \neq i$—writes:

$$\ddot{q}_i + \left[ \omega_i^2 + \sum_{p=1}^{P} \frac{(\chi_i^{(p)})^2}{C^{(p)}} \right] q_i + \sum_{p=1}^{P} \frac{\chi_i^{(p)}}{C^{(p)}} Q^{(p)} = F_i \qquad (32)$$

In open circuit ($Q^{(p)}=0$ for all $p$), the problem reads

$$\ddot{q}_i + \hat{\omega}_i^2 q_i = F_i \qquad (33)$$

so that, by comparing Equations (32) and (33),

$$\hat{\omega}_i \approx \omega_i^2 + \sum_{p=1}^{P} \frac{(\chi_i^{(p)})^2}{C^{(p)}} \qquad (34)$$

The above expression gives an approximation of the system $i$th mode natural frequency in open circuit, denoted by $\hat{\omega}_i$, provided that the above introduced one-mode truncation is valid. As a consequence, by mixing Equations (31) and (34), the systems $i$th mode EEMCF is

$$k_{\text{eff},i}^2 \approx \sum_{p=1}^{P} (k_i^{(p)})^2 \qquad (35)$$

where $k_i^{(p)}$, defined by

$$k_i^{(p)} = \frac{\chi_i^{(p)}}{\sqrt{C^{(p)}} \omega_i} \qquad (36)$$

is the $i$th mode and the $p$th piezoelectric patch dimensionless coupling factor. By considering the definitions of $\chi_i^{(p)}, \omega_i$ and $C^{(p)}$ as a function of the matrices of the initial problem (20), one
obtains another expression for $k_{\text{eff},i}$ stemming from Equation (35):

$$k_{\text{eff},i}^2 \approx \frac{\Phi_i^T \mathbf{K}_c \mathbf{K}_c^{-1} \mathbf{K}_c^T \Phi_i}{\Phi_i^T \mathbf{K}_m \Phi_i} \quad (37)$$

If only one piezoelectric patch is connected to the structure, one can observe that $k_{\text{eff},i}$ is directly linked to the only coupling coefficient $\chi_i$:

$$k_{\text{eff},i}^2 \approx k_i^2 = \frac{\chi_i^2}{C \omega_i^2} \quad (38)$$

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_{\text{eff},i}^{(p)}$, which is associated with the $i$th mode and the $p$th piezoelectric patch, by:

$$(k_{\text{eff},i}^{(p)})^2 = \frac{(\hat{\omega}_i^{(p)})^2 - \omega_i^2}{\omega_i^2} \quad (39)$$

In the above equation, $\hat{\omega}_i^{(p)}$ denotes the $i$th system natural frequency with the $p$th piezoelectric patch in open-circuit condition and all the others short circuited. Under these conditions, all the voltages but the $p$th are zero and $Q^{(p)}=0$. By truncating Equations (28a), (28b) to the $i$th mode, one gets the following approximation of $\hat{\omega}_i^{(p)}$:

$$(\hat{\omega}_i^{(p)})^2 \approx \omega_i^2 + \frac{(\chi_i^{(p)})^2}{C^{(p)}} \quad (40)$$

so that $k_{\text{eff},i}^{(p)}$ is found close to the coupling factor defined by Equation (36):

$$k_{\text{eff},i}^{(p)} \approx k_i^{(p)} = \frac{\chi_i^{(p)}}{\sqrt{C^{(p)}} \omega_i} \quad (41)$$

The above equation shows that the coupling coefficient $\chi_i^{(p)}$ is directly linked to the EEMCF $k_{\text{eff},i}^{(p)}$. Coupling factor $k_i^{(p)}$ has been previously introduced in [37].

The above calculations have introduced several coupling coefficients and factors, whose main features are summarized now. The modal coupling factors $k_i^{(p)}$ (defined by Equation (36)) play the same role as $\chi_i^{(p)}$ (defined by Equation (29)), but have the advantage of being dimensionless, which thus gives them a larger relevance. Moreover, they are directly linked to the EEMCF $k_{\text{eff},i}^{(p)}$ (defined by Equation (39)), and thus to $\hat{\omega}_i^{(p)}$ and $\omega_i$. Finally, observing Equation (35) leads to remark that the global EEMCF $k_{\text{eff},i}$ groups all the contributions of the individual piezoelectric patches, since $k_i^{(p)}$ is linked to the $p$th 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–26]. In this context, the approximate relation (41), in addition to giving a physical
interpretation of $k_i^{(p)}$ and $\chi_i^{(p)}$, gives a useful practical way, alternative to Equation (29), for evaluating the coupling factors $k_i^{(p)}$ (and $\chi_i^{(p)}$), provided that the natural frequencies $\omega_i$ and $\hat{\omega}_i^{(p)}$ are known. It can be used in two practical cases:

* In an experimental situation, each $k_i^{(p)}$ can be obtained by simply measuring $\omega_i$ and the $\hat{\omega}_i^{(p)}$, $p=1, \dots, P$.

* If a given finite element code is able to compute each $\hat{\omega}_i^{(p)}$, $\omega_i$ being given by a standard elastic modal analysis of the whole structure.

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_i^{(p)}$.

## 5. 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.

### 5.1. 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

$$u_x(x, z, t) = u(x, t) + z\theta(x, t) \qquad (42a)$$

$$u_z(x, z, t) = w(x, t) \qquad (42b)$$

<figure><img src="image_4.png" /><figcaption>Figure 3. Piezoelectric laminated beam.</figcaption></figure>
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 $\theta$ the fiber rotation defined by

$$ \theta = -\frac{\partial w}{\partial x} \qquad (43) $$

From Equations (5) and the above-described kinematics, axial strain can be written as follows:

$$ \varepsilon_1 = e + z\kappa \qquad (44) $$

where membrane strain $e$ and curvature $\kappa$ are defined by

$$ e = \frac{\partial u}{\partial x} \quad \text{and} \quad \kappa = -\frac{\partial^2 w}{\partial x^2} \qquad (45) $$

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 $p$th patch. The only non-vanishing component of the electric field is:

$$ E_3^{(p)} = -\frac{V^{(p)}}{h^{(p)}} \qquad (46) $$

## 5.2. 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 $\varepsilon_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 ($\sigma_2 = \sigma_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:

$$ \sigma_1^k = \bar{c}_{11}^k \varepsilon_1 - \bar{e}_{31}^k E_3 \qquad (47) $$

$$ D_3^k = \bar{e}_{31}^k \varepsilon_1 + \bar{e}_{33}^k E_3 \qquad (47) $$

where $\sigma_1^k$ and $D_3^k$ are the axial stress and the transverse electrical displacement within the $k$th layer. $\bar{c}_{11}^k$ is the modified elastic stiffness at constant electric field, $\bar{e}_{31}^k$ is the modified piezoelectric constant and $\bar{e}_{33}^k$ is the modified dielectric permittivity at constant strain. They are defined by:

$$ \bar{c}_{11}^k = Y_1^k = c_{11} + \frac{c_{12}^2 c_{33} + c_{13}^2 c_{11} - 2 c_{12} c_{13} c_{23}}{\Delta} \qquad (48a) $$

$$ \bar{e}_{31}^k = e_{31} - \frac{c_{12}c_{33} - c_{13}c_{23}}{\Delta} e_{31} - \frac{c_{13}c_{11} - c_{12}c_{23}}{\Delta} e_{33} \qquad (48b) $$

$$ \bar{e}_{33}^k = \epsilon_{33}^k + \frac{c_{33}e_{31}^2 + c_{11}e_{33}^2 - 2c_{23}e_{31}e_{33}}{\Delta} \qquad (48c) $$
with

$$ \Delta = c_{11}c_{33} - c_{23}^2 \quad (49) $$

and $Y_1^k$ 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.

### 5.3. 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:

$$ \sum_{k=1}^{K} \int_{\Omega^k} \rho^k (\ddot{u}_x \delta u_x + \ddot{u}_z \delta u_z) d\Omega + \sum_{k=1}^{K} \int_{\Omega^k} \bar{c}_{11}^k \varepsilon_1 \delta \varepsilon_1 d\Omega + \sum_{p=1}^{P} \frac{V^{(p)}}{h^{(p)}} \int_{\Omega^{(p)}} \bar{e}_{31}^{(p)} \delta \varepsilon_1 d\Omega \\ = \sum_{k=1}^{K} \int_{\Gamma_t^k} (t_x^k \delta u_x + t_z^k \delta u_z) dS + \sum_{k=1}^{K} \int_{\Omega^k} (f_x^k \delta u_x + f_z^k \delta u_z) d\Omega \quad (50) $$

$$ -\sum_{p=1}^{P} \frac{\delta V^{(p)}}{h^{(p)}} \int_{\Omega^{(p)}} \bar{e}_{31}^{(p)} \varepsilon_1 d\Omega + \sum_{p=1}^{P} V^{(p)} C^{(p)} \delta V^{(p)} = \sum_{p=1}^{P} Q^{(p)} \delta V^{(p)} \quad (51) $$

where $\rho^k$ and $\Omega^k$ are the mass density and the domain occupied by the kth layer and the pth patch capacitance is (Equation (14)):

$$ C^{(p)} = \frac{bl^{(p)}}{h^{(p)}} \bar{e}_{33}^{(p)} \quad (52) $$

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_x^k, t_z^k$) and to body ones ($f_x^k, f_z^k$). 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 \times 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:

$$ K_c = (K_c^{(1)} K_c^{(2)} \dots K_c^{(P)}) \quad (53) $$

where $K_c^{(p)}$ 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_c^{(p)}$ is simply obtained by assembling the elementary vectors $K_c^{e(p)}$ of size $(6 \times 1)$, whose explicit expression is given in the Appendix (Equation (A7)). Owing to their particular form, during the assembling operation, all components of $K_c^{(p)}$ vanish except those associated with the end nodes of the pth patch. An example is given on Figure 5. It shows that the electromechanical coupling occurs
between the potential difference $V^{(p)}$ and the mechanical dofs (axial translation $u_i$ and rotation $\theta_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].

# 6. 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.

## 6.1. 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 $\rho_b$ and Young's modulus $Y_b$. The piezoelectric material constants are its density $\rho_p$, Young's modulus $Y_p$, coupling piezoelectric constant $\bar{e}_{31}$ and dielectric constant $\bar{\epsilon}_{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.

The piezoelectric material parameters are those of a PIC151 [41]. This reference gives the following constants: $s_{11}^E=15 \times 10^{-12} \text{ Pa}^{-1}$, $d_{31}=-210 \times 10^{-12} \text{ C/N}$ and $\epsilon_{33}^\sigma=2400\epsilon_0$ corresponding

<figure><img src="image_5.png" /><figcaption>Figure 4. Beam with two piezoelectric patches.</figcaption></figure>
<table><caption>Table I. Numerical values of system parameters.</caption><thead><tr><th>Parameters</th><th>Beam</th><th>Values</th><th>Piezo.</th><th>Values</th></tr></thead><tbody><tr><td>Lengths (mm)</td><td>l<sub>b</sub></td><td>170</td><td>l<sub>p</sub></td><td>25</td></tr><tr><td>Thicknesses (mm)</td><td>h<sub>b</sub></td><td>2</td><td>h<sub>p</sub></td><td>0.5</td></tr><tr><td>Width (mm)</td><td>b</td><td>20</td><td>b</td><td>20</td></tr><tr><td>Patches position (mm)</td><td></td><td></td><td>x<sub>-</sub></td><td>0.5</td></tr><tr><td>Densities (kg/m<sup>3</sup>)</td><td>ρ<sub>b</sub></td><td>2800</td><td>ρ<sub>p</sub></td><td>8500</td></tr><tr><td>Young's moduli (GPa)</td><td>Y<sub>b</sub></td><td>72</td><td>Y<sub>p</sub></td><td>66.7</td></tr><tr><td>Piezoelectric constant (C/m<sup>2</sup>)</td><td></td><td></td><td>ε&#773;<sub>31</sub></td><td>-14</td></tr><tr><td>Dielectric constant (F/m)</td><td></td><td></td><td>ε&#773;<sub>33</sub></td><td>2068ε<sub>0</sub></td></tr></tbody></table>

$ε_0 = 8.854 \times 10^{-12} \text{F/m}$ is the free space permittivity [31].

to a constitutive law ($\sigma$, $D$) as a function of ($\epsilon$, $E$). The following relations have been used to obtain the parameters of Table I, linked to Equations (47):

$$Y_p = \bar{c}_{11} = \frac{1}{s_{11}^E}, \quad \bar{e}_{31} = \frac{d_{31}}{s_{11}^E}, \quad \bar{\epsilon}_{33} = \epsilon_{33}^\sigma - \frac{d_{31}^2}{s_{11}^E} \qquad (54)$$

For the PIC151, the material coupling factors are [42]:

$$k_{31}^2 = \frac{d_{31}^2}{s_{11}^E \epsilon_{33}^\sigma}, \quad \tilde{k}_{31}^2 = \frac{\bar{e}_{31}^2}{\bar{c}_{11} \bar{\epsilon}_{33}}, \quad 1 + \tilde{k}_{31}^2 = \frac{1}{1 - k_{31}^2} \qquad (55)$$

In our case, $k_{31} = 0.372$ and $\tilde{k}_{31} = 0.401$.

## 6.2. 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 \in [0, x_-]$) and region III ($x \in [x_+, l_b]$) are composed of one layer of beam material of thickness $h_b$. Region II corresponds to the piezoelectric patches (for $x \in [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, \theta$), 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:

$$M_m \ddot{U} + K_m U + K_c V = F \qquad (56a)$$

$$K_e V - Q - K_c^T U = 0 \qquad (56b)$$

In the above equations, $U$ is of size $N \times 1$, $M_m$ and $K_m$ are of size $N \times N$. The electric d.o.f., associated with the two piezoelectric patches, lead to define the voltage vector $V = (V_1 V_2)^T$ and the
charge vector $Q=(Q_1 Q_2)^T$, both of size two, where $(V_1, Q_1)$ and $(V_2, Q_2)$ are the voltage/charge pairs associated, respectively, with the upper and the lower patch. Matrices $K_e$ and $K_c$ write:

$$ K_e = \begin{pmatrix} C & 0 \\ 0 & C \end{pmatrix}, \quad K_c = (K_1 \ K_2), \quad C = \frac{\bar{\epsilon}_{33}bl_p}{h_p} \qquad (57) $$

where $C$ is the patches capacitance, identical for both patches, defined either by Equation (14) or by Equation (52). $K_c$ is a $N \times 2$ matrix composed of two column vectors $K_1=K_c^{(1)}$ and $K_2=K_c^{(2)}$ (see Equation (53)) of length $N$. $K_1$ (resp. $K_2$) couples the mechanical d.o.f. to the upper (resp. lower) piezoelectric patches voltage. By considering the geometry of region II of the beam as well as the opposite polarization direction of the piezoelectric patches, one obtains $G_0^{(1)}=-G_0^{(2)}=b\bar{\epsilon}_{31}$ and $G_1^{(1)}=G_1^{(2)}=b\bar{\epsilon}_{31}(h_p+h_b)/2$ in Equation (A7). Consequently, it shows that (i) $K_1$ and $K_2$ are equal in absolute value, (ii) $K_1$ and $K_2$ elements corresponding to d.o.f. $u$ are opposite in sign and (iii) $K_1$ and $K_2$ elements corresponding to d.o.f. $\theta$ are equal. It means that if the piezoelectric patches are subjected to two identical voltages ($V_1=V_2$), only the flexural motion is activated; conversely, two opposite voltages ($V_1=-V_2$) activate only the axial motion.

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:

$$ CV - 2Q - (K_1^T + K_2^T)U = 0 \qquad (58) $$

By considering the RL electrical circuit of Figure 4, one can write:

$$ V = -R\dot{Q} - L\ddot{Q} \qquad (59) $$

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:

$$ M\ddot{X} + D\dot{X} + KX = \tilde{F} \qquad (60) $$

with the following notations:

$$ X = \begin{pmatrix} U \\ Q \end{pmatrix}, \quad M = \begin{pmatrix} M_m & 0 \\ 0^T & L \end{pmatrix}, \quad D = \begin{pmatrix} 0 & 0 \\ 0^T & R \end{pmatrix} \qquad (61a) $$

$$ K = \begin{pmatrix} \hat{K}_m & (K_1 + K_2)/C \\ (K_1^T + K_2^T)/C & 2/C \end{pmatrix}, \quad \tilde{F} = \begin{pmatrix} F \\ 0 \end{pmatrix} \qquad (61b) $$

where the stiffness matrix in open circuit—i.e. with $Q=0$—is:

$$ \hat{K}_m = K_m + \frac{K_1 K_1^T + K_2 K_2^T}{C} \qquad (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
<figure><img src="image_6.png" /><figcaption>Figure 5. Sparsity pattern of matrices $\mathbf{K}_m$, $\mathbf{K}_c$ and $\mathbf{K}$, from left to right.</figcaption></figure>

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 $\mathbf{K}$ and $\hat{\mathbf{K}}_m$ as functions of $C$, $\mathbf{K}_1$ and $\mathbf{K}_2$ are slightly different.

### 6.3. 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 $(\omega_i, \mathbf{\Phi}_i)$ and can be calculated with Equation (24). By considering Equations (29) and (57), the coupling coefficients $\chi_i^{(p)}$, $p=1, 2$, in Equations (28a), (28b) are

$$ \chi_i^{(p)} = \mathbf{\Phi}_i^T \mathbf{K}_p, \quad p=1,2 \qquad (63) $$

Then, since $\mathbf{\Phi}_i$ are flexural modes, they have zero elements on the lines corresponding to d.o.f. $u$, which proves that $\chi_i^{(1)} = \chi_i^{(2)} = \chi_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 $\tilde{N}$ flexural modes becomes in this particular case:

$$ \mathbf{U}(t) = \sum_{i=1}^{\tilde{N}} \mathbf{\Phi}_i q_i(t) \qquad (64a) $$

$$ \ddot{q}_i + \omega_i^2 q_i + \sum_{j=1}^{\tilde{N}} \frac{2\chi_i \chi_j}{C} q_j + \frac{2\chi_i}{C} Q = F_i, \quad \forall i \in \{1 \dots \tilde{N}\} \qquad (64b) $$

$$ L\ddot{Q} + R\dot{Q} + \frac{2}{C}Q + \sum_{j=1}^{\tilde{N}} \frac{2\chi_j}{C}q_j = 0 \qquad (64c) $$
<table><caption>Table II. Natural frequencies and coupling coefficients, from finite element computation (FE) and experiments (Exp.).</caption><thead><tr><th rowspan="2">Modes</th><th colspan="2">Short-circuit freq. (Hz)</th><th colspan="2">Open-circuit freq. (Hz)</th><th colspan="3">Coupling factors</th></tr><tr><th>f<sub>i</sub> FE</th><th>f<sub>i</sub> Exp.</th><th>f&#x0302;<sub>i</sub> FE</th><th>f&#x0302;<sub>i</sub> Exp.</th><th>k<sub>i</sub> FE</th><th>k<sub>eff,i</sub> FE</th><th>k<sub>eff,i</sub> Exp.</th></tr></thead><tbody><tr><td>1 flex. (F<sub>1</sub>)</td><td>48.96</td><td>51.64</td><td>49.42</td><td>52.17</td><td>0.143</td><td>0.137</td><td>0.144</td></tr><tr><td>2 flex. (F<sub>2</sub>)</td><td>337.1</td><td>337.0</td><td>340.7</td><td>340.2</td><td>0.150</td><td>0.145</td><td>0.138</td></tr><tr><td>1 tors. (T<sub>1</sub>)</td><td>—</td><td>853.0</td><td>—</td><td>854.0</td><td>—</td><td>—</td><td>0.050</td></tr><tr><td>3 flex. (F<sub>3</sub>)</td><td>951.8</td><td>936.3</td><td>960.6</td><td>940.0</td><td>0.140</td><td>0.137</td><td>0.089</td></tr></tbody></table>

The present reduced-order formulation depends only on the $\tilde{N}$ retained short-circuit modes ($\omega_i$, $\Phi_i$), the $\tilde{N}$ associated coupling coefficients $\chi_i$ and the patches capacitance $C$. To apply Section 4.3 results, we can define the following coupling factors:

$$k_{\text{eff},i}^2 = \frac{\hat{\omega}_i^2 - \omega_i^2}{\omega_i^2}, \quad k_i = \sqrt{\frac{2}{C}} \frac{\chi_i}{\omega_i} \qquad (65)$$

where $\hat{\omega}_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:

$$k_{\text{eff},i} \approx k_i \qquad (66)$$

which gives another way to evaluate $\chi_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:

$$f_i = \frac{\omega_i}{2\pi}, \quad \hat{f}_i = \frac{\hat{\omega}_i}{2\pi} \qquad (67)$$

The short-circuit natural frequencies ($f_i$, second column) and mode shapes ($\Phi_i$, Figure 6) have been calculated by numerically solving the eigenvalue problem of Equation (24). The open-circuit natural frequencies ($\hat{f}_i$, fourth column) and mode shapes ($\hat{\Phi}_i$, Figure 6) have been obtained by solving the same eigenvalue problem, but with $\hat{\mathbf{K}}_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 $\chi_i$ obtained with Equation (63). Finally, the seventh column gives the EEMCF $k_{\text{eff},i}$, obtained with Equation (65). The values obtained for $k_i$ and $k_{\text{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.
<figure><img src="image_7.png" /><figcaption>Figure 6. First three short-circuit (‘—’) and open-circuit (‘--’) mode shapes, from the finite element model, normalized with Equation (26). The piezoelectric patches are displayed by thick lines.</figcaption></figure>

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 ($\hat{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 than 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.

### 6.4. 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 \Omega t$. Defining the complex amplitude $\mathbf{X}$ by $\mathbf{X}(t)=\Re(\mathbf{X}e^{j\Omega t})$ leads to the frequency formulation of Equation (60):

$$Z(\Omega)\mathbf{X}=F_0 \qquad (68)$$
<figure><img src="image_8.png" /><figcaption>Figure 7. Beam tip frequency response functions. ‘-’: experimental; ‘--’: finite element. (a) Short-circuit FRFs and (b) short circuit and damped FRFs, around mode 2 ($F_2$) resonance.</figcaption></figure>

with the impedance matrix defined by:

$$Z(\Omega) = -\Omega^2 M + j\Omega D + K \quad (69)$$

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 $\Omega$:

$$H_{ij}(\Omega) = \frac{X_i}{F_{0j}} = Z_{ij}^{-1}(\Omega) \quad (70)$$

where $Z_{ij}^{-1}(\Omega)$ is the ith line and the jth column term of $Z(\Omega)$ inverse matrix.

Figure 7 shows the driving point FRF modulus ($|H_{ii}(\Omega)|$) 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 853Hz, 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\Omega$ and $L=21.8H$ tuned to lower the second flexural resonance, around 336Hz. In the experiments, an attenuation of 25dB 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.

### 6.5. Reduced-order model

The system vibratory response can also be obtained with the modal ROM defined by Equations (64), truncated to $\tilde{N}$ modal coordinates. The problem to solve takes the form of Equation (60), where
<figure><img src="image_9.png" /><figcaption>Figure 8. Beam tip frequency response functions with RL shunt tuned on mode $F_2$. ‘-’: finite element curve of Figure 7(b); ‘--’: modal model truncated to $\tilde{N}=1, 2, 5, 12$ modes.</figcaption></figure>

the matrices of Equations (61) are replaced by

$$ \mathbf{X} = (q_1 \dots q_{\tilde{N}} Q)^{\mathrm{T}}, \quad \mathbf{M} = \begin{pmatrix} \mathbf{1} & \mathbf{0} \\ \mathbf{0}^{\mathrm{T}} & L \end{pmatrix}, \quad \mathbf{D} = \begin{pmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0}^{\mathrm{T}} & R \end{pmatrix} \qquad (71a) $$

$$ \mathbf{K} = \begin{pmatrix} \hat{\mathbf{K}}_{\mathrm{m}} & 2\chi/C \\ 2\chi^{\mathrm{T}}/C & 2/C \end{pmatrix}, \quad \mathbf{F} = (F_1 \dots F_{\tilde{N}} Q)^{\mathrm{T}} \qquad (71b) $$

where $\chi=(\chi_1 \dots \chi_{\tilde{N}})^{\mathrm{T}}$ and the ith row and jth column component of $\hat{\mathbf{K}}_{\mathrm{m}}$ is $[\hat{\mathbf{K}}_{\mathrm{m}}]_{ij} = \omega_i^2 \delta_{ij} + 2\chi_i \chi_j / C$.

Figure 8 shows the driving point FRF modulus ($|H_{ii}(\Omega)|$), calculated with the full FE model (as explained in Section 6.4) as well as with the present modal model, truncated to $\tilde{N}=1, 2, 5$ and 12 modes. A good convergence is obtained with at least 5 modes at all points (even at the anti-resonances), except around the attenuated $F_2$ resonance, where convergence is attained with $\tilde{N}=12$ modes. The full model is composed of $N=126$ d.o.f. plus one electrical d.o.f. (charge $Q$). As compared with the 13 d.o.f. modal model, a 30 ratio in computation time is gained.

## 7. 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
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–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. Full finite
element simulations are possible on systems with a small number of d.o.f. [20]. However, for
systems with a large number of d.o.f., ROMs are essential to obtain the vibratory response within
realistic computation time and to perform optimizations.

# APPENDIX A

Details of the stratified beam finite element discretization of Section 5 are specified here, by starting
from the variational formulation of Equations (50) and (51). In the following, a superscript or a
subscript *e* refers to elementary quantities. The generalized displacements are discretized using linear (axial displacement) and cubic (transverse displacement) shape functions. Thus, axial and transverse displacements are related to the elementary d.o.f. vector **U**<sup>*e*</sup> by

$$ u^e = N_u U^e \quad \text{and} \quad w^e = N_w U^e \qquad (A1) $$

where

$$ \mathbf{U}^e = (u_1 \ w_1 \ \theta_1 \ | \ u_2 \ w_2 \ \theta_2)^T $$

$u_1$, $w_1$ and $\theta_1$ (respectively, $u_2$, $w_2$ and $\theta_2$) corresponding to the first (respectively, second) node of the beam element.

In equation (A1), the interpolation vectors are defined as follows:

$$ \mathbf{N}_u = (N_1 \ 0 \ 0 \ | \ N_2 \ 0 \ 0) $$

and

$$ \mathbf{N}_w = (0 \ N_3 \ N_4 \ | \ 0 \ N_5 \ N_6) $$

with, for all $x \in [0, L_e]$

$$ N_1 = 1 - \frac{x}{L_e}, \quad N_2 = \frac{x}{L_e}, \quad N_3 = \left(1 - \frac{x}{L_e}\right)^2 \left(1 + \frac{2x}{L_e}\right) $$

$$ N_4 = x \left(1 - \frac{x}{L_e}\right)^2, \quad N_5 = \frac{x^2}{L_e^2} \left(3 - \frac{2x}{L_e}\right), \quad N_6 = \frac{x^2}{L_e} \left(\frac{x}{L_e} - 1\right) $$

where $L_e$ is the length of the considered finite element. Note that Equations (43), (45) and (A1) lead to the following expression of elementary rotation, membrane strain and curvature:

$$ \theta^e = -\mathbf{N}'_w \mathbf{U}^e \quad \text{and} \quad e^e = \mathbf{N}'_u \mathbf{U}^e \quad \text{and} \quad \kappa^e = -\mathbf{N}''_w \mathbf{U}^e \qquad (A2) $$

where $(\cdot)' = \partial(\cdot)/\partial x$.

The various terms appearing in the variational formulation of Equations (50) and (51) are now successively discussed.

* The kinetic energy variation is:

$$ \begin{aligned} \sum_{k=1}^{K} \int_{\Omega^k} \rho^k (\ddot{u}_x \delta u_x + \ddot{u}_z \delta u_z) \, d\Omega &= \int_{x^{-}}^{x^{+}} [I_0(\ddot{u}\delta u + \ddot{w}\delta w) + I_1(\ddot{u}\delta\theta + \ddot{w}\delta u) + I_2\ddot{\theta}\delta\theta] \, dx \\ &= \int_{x^{-}}^{x^{+}} (\delta u \ \delta w \ \delta \theta) \begin{pmatrix} I_0 & 0 & I_1 \\ 0 & I_0 & 0 \\ I_1 & 0 & I_2 \end{pmatrix} \begin{pmatrix} \ddot{u} \\ \ddot{w} \\ \ddot{\theta} \end{pmatrix} \, dx \end{aligned} $$

where the zero-, first-, and second-order inertia moments are defined by

$$ [I_0, 1, I_2] = \sum_{k=1}^{K} \rho^k [J_0^k, J_1^k, J_2^k] $$
with the zero, first and second moment of area of each sublayer, respectively, defined by

$$[J_0^k, J_1^k, J_2^k] = b \int_{z_{k-1}}^{z_k} [1, z, z^2] dz = b \left[ h^k, \frac{z_k^2 - z_{k-1}^2}{2}, \frac{z_k^3 - z_{k-1}^3}{3} \right]$$

where b and $h^k$ are the width of the beam and the thickness of the kth layer. The elementary mass matrix is then:

$$\mathbf{M}_{\mathrm{m}}^{\mathrm{e}}=\int_{0}^{L_{\mathrm{e}}}\left(\mathbf{N}_{u}^{\mathrm{T}} \mathbf{N}_{w}^{\mathrm{T}}-\mathbf{N}_{w}^{\prime \mathrm{T}}\right)\left(\begin{array}{lll}
I_{0} & 0 & I_{1} \\
0 & I_{0} & 0 \\
I_{1} & 0 & I_{2}
\end{array}\right)\left(\begin{array}{l}
\mathbf{N}_{u} \\
\mathbf{N}_{w} \\
-\mathbf{N}_{w}^{\prime}
\end{array}\right) \mathrm{d} x$$

whose explicit expression is:

$$\mathbf{M}_{\mathrm{m}}^{\mathrm{e}}=\frac{I_{0} L_{\mathrm{e}}}{420}\left[
\begin{array}{cccccc}
140 & 0 & 0 & 70 & 0 & 0 \\
0 & 156 & 22 L_{\mathrm{e}} & 0 & 54 & -13 L_{\mathrm{e}} \\
0 & 22 L_{\mathrm{e}} & 4 L_{\mathrm{e}}^{2} & 0 & 13 L_{\mathrm{e}} & -3 L_{\mathrm{e}}^{2} \\
70 & 0 & 0 & 140 & 0 & 0 \\
0 & 54 & 13 L_{\mathrm{e}} & 0 & 156 & -22 L_{\mathrm{e}} \\
0 & -13 L_{\mathrm{e}} & -3 L_{\mathrm{e}}^{2} & 0 & -22 L_{\mathrm{e}} & 4 L_{\mathrm{e}}^{2}
\end{array}
\right] \\
+\frac{I_{1}}{12}\left[
\begin{array}{cccccc}
0 & 6 & -L_{\mathrm{e}} & 0 & -6 & L_{\mathrm{e}} \\
6 & 0 & 0 & 6 & 0 & 0 \\
-L_{\mathrm{e}} & 0 & 0 & L_{\mathrm{e}} & 0 & 0 \\
0 & 6 & L_{\mathrm{e}} & 0 & -6 & -L_{\mathrm{e}} \\
-6 & 0 & 0 & -6 & 0 & 0 \\
L_{\mathrm{e}} & 0 & 0 & -L_{\mathrm{e}} & 0 & 0
\end{array}
\right] \\
+\frac{I_{2}}{30 L_{\mathrm{e}}}\left[
\begin{array}{cccccc}
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 36 & 3 L_{\mathrm{e}} & 0 & -36 & 3 L_{\mathrm{e}} \\
0 & 3 L_{\mathrm{e}} & 4 L_{\mathrm{e}}^{2} & 0 & -3 L_{\mathrm{e}} & -L_{\mathrm{e}}^{2} \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & -36 & -3 L_{\mathrm{e}} & 0 & 36 & -3 L_{\mathrm{e}} \\
0 & 3 L_{\mathrm{e}} & -L_{\mathrm{e}}^{2} & 0 & -3 L_{\mathrm{e}} & 4 L_{\mathrm{e}}^{2}
\end{array}
\right] \quad (\mathrm{A}3)$$
• The mechanical contribution to the internal energy variation is:

$$
\begin{align*}
\sum_{k=1}^{K} \int_{\Omega^k} \bar{c}_{11}^k \varepsilon_1 \delta \varepsilon_1 \mathrm{d}\Omega &= \int_{x^{-}}^{x^{+}} [A \varepsilon \delta e + B(\kappa \delta e + e \delta \kappa) + D \kappa \delta \kappa] \mathrm{d}x \\
&= \int_{x^{-}}^{x^{+}} (\delta e \ \delta \kappa) \begin{pmatrix} A & B \\ B & D \end{pmatrix} \begin{pmatrix} e \\ \kappa \end{pmatrix} \mathrm{d}x
\end{align*}
$$

where A is the extensional stiffness, D is the bending stiffness, and B is the bending-
extensional coupling stiffness defined as

$$
[A, B, D] = \sum_{k=1}^{K} \bar{c}_{11}^{k} [J_{0}^{k}, J_{1}^{k}, J_{2}^{k}]
$$

The elementary stiffness matrix is then:

$$
\mathbf{K}_{\text{m}}^{\text{e}} = \int_{0}^{L_{e}} (\mathbf{N}_{u}^{\prime T} - \mathbf{N}_{w}^{\prime\prime T}) \begin{pmatrix} A & B \\ B & D \end{pmatrix} \begin{pmatrix} \mathbf{N}_{u}^{\prime} \\ -\mathbf{N}_{w}^{\prime\prime} \end{pmatrix} \mathrm{d}x
$$

whose explicit expression is:

$$
\mathbf{K}_m^e = \frac{A}{L_e} \begin{bmatrix} 1 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} + \frac{B}{L_e} \begin{bmatrix} 0 & 0 & -1 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ -1 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & -1 & 0 & 0 \end{bmatrix} \\
+ \frac{D}{L_e^3} \begin{bmatrix} 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 12 & 6L_e & 0 & -12 & 6L_e \\ 0 & 6L_e & 4L_e^2 & 0 & -6L_e & 2L_e^2 \\ 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -12 & -6L_e & 0 & 12 & -6L_e \\ 0 & 6L_e & 2L_e^2 & 0 & -6L_e & 4L_e^2 \end{bmatrix} \tag{A4}
$$

• Work done by the external mechanical forces

$$
\begin{align*}
& \sum_{k=1}^{K} \int_{\Gamma_t^k} (t_x^k \delta u_x + t_z^k \delta u_z) \, dS + \sum_{k=1}^{K} \int_{\Omega^k} (f_x^k \delta u_x + f_z^k \delta u_z) \, d\Omega \\
&= \sum_{k=1}^{K} \left\{ \left[ \int_{\Sigma^k} (t_x^k \delta u + t_z^k \delta w - z t_x^k \delta \theta) \, dS \right]_{x^{-}}^{x^{+}} + \int_{\Omega^k} (f_x^k \delta u + f_z^k \delta w - z f_x^k \delta \theta) \, d\Omega \right\} \\
&= [N \delta u + T \delta w - M \delta \theta]_{x^{-}}^{x^{+}} + \int_{x^{-}}^{x^{+}} (n \delta u + t \delta w - m \delta \theta) \, dx \tag{A5}
\end{align*}
$$
where $\Sigma^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:

$$[N, T, M] = \sum_{k=1}^{K} \left[ \int_{\Sigma^k} t_x^k dS, \int_{\Sigma^k} t_z^k dS, \int_{\Sigma^k} zt_z^k dS \right]$$

and

$$[n, t, m] = \sum_{k=1}^{K} \left[ \int_{\Sigma^k} f_x^k dS, \int_{\Sigma^k} f_z^k dS, \int_{\Sigma^k} zf_z^k dS \right]$$

The elementary vector of generalized mechanical forces is then:

$$\mathbf{F}^e = [NN_u^T + TN_w^T + MN_w^T]_0^{L_e} + \int_0^{L_e} (nN_u^T + tN_w^T + mN_w^T) dx$$

whose explicit expression is:

$$\mathbf{F}^e = \begin{bmatrix} N \\ T \\ M \\ 0 \\ 0 \\ 0 \end{bmatrix}_{x=0} + \begin{bmatrix} 0 \\ 0 \\ 0 \\ N \\ T \\ M \end{bmatrix}_{x=L_e} + \frac{nL_e}{2} \begin{bmatrix} 1 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{bmatrix} + \frac{tL_e}{12} \begin{bmatrix} 0 \\ 6 \\ L_e \\ 0 \\ 6 \\ -L_e \end{bmatrix} + m \begin{bmatrix} 0 \\ -1 \\ 0 \\ 0 \\ 1 \\ 0 \end{bmatrix} \quad (A6)$$

* The piezoelectric contributions of the internal energy variation, related to the inverse effect, are:

$$\sum_{p=1}^{P} \frac{V^{(p)}}{h^{(p)}} \int_{\Omega^{(p)}} \bar{e}_{31}^{(p)} \delta\epsilon_1 d\Omega = \sum_{p=1}^{P} \left[ \int_{x^{-}}^{x^{+}} G_{0}^{(p)} V^{(p)} \delta e \, dx + \int_{x^{-}}^{x^{+}} G_{1}^{(p)} V^{(p)} \delta \kappa \, dx \right] \\ = \sum_{p=1}^{P} \int_{x^{-}}^{x^{+}} (\delta e \ \delta \kappa) \begin{pmatrix} G_{0}^{(p)} \\ G_{1}^{(p)} \end{pmatrix} V^{(p)} \, dx$$

where $G_0^{(p)}$ and $G_1^{(p)}$ are defined by

$$[G_0^{(p)}, G_1^{(p)}] = \frac{\bar{e}_{31}^{(p)}}{h^{(p)}} [J_0^{(p)}, J_1^{(p)}]$$
Then, following Equation (53), the elementary coupling vector is:

$$ \mathbf{K}_e^{(p)} = \int_0^{L_e} (\mathbf{N}_u^T - \mathbf{N}_w^{\prime\prime T}) \begin{pmatrix} G_0^{(p)} \\ G_1^{(p)} \end{pmatrix} dx = G_0^{(p)} \begin{bmatrix} -1 \\ 0 \\ 0 \\ 1 \\ 0 \\ 0 \end{bmatrix} + G_1^{(p)} \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \\ 0 \\ -1 \end{bmatrix} \quad (A7) $$

The same matrix is obtained for the direct effect.

ACKNOWLEDGEMENTS

The authors are grateful to the reviewers for their careful reading of the manuscript and for their various relevant comments. Aurélien Sénéchal is also warmly thanked, also for his careful reading, and for the forgotten plus sign discovery. This research is partially funded under DGA contract number 05.43.063.004707565 (INSA Lyon/LGEF, CNAM Paris/LMSSc, UCBL Lyon/LENAC), for which the authors are grateful.

REFERENCES

1. Allik H, Hughes TJR. Finite element method for piezoelectric vibration. *International Journal for Numerical Methods in Engineering* 1970; **2**(2):151–157.

2. Tzou HS, Tseng CI. Distributed piezoelectric sensor/actuator design for dynamic measurement/control of distributed parameter systems: A piezoelectric finite element approach. *Journal of Sound and Vibration* 1990; **138**(1):17–34.

3. Hwang WS, Park HC. Finite element modeling of piezoelectric sensors and actuators. *AIAA Journal* 1993; **31**(5):930–937.

4. Vasques CMA, Rodrigues JD. Coupled three-layered analysis of smart piezoelectric beams with different electric boundary conditions. *International Journal for Numerical Methods in Engineering* 2005; **62**(11):1488–1518.

5. Kim J, Varadan VV, Varadan VK. Finite element modelling of structures including piezoelectric active devices. *International Journal for Numerical Methods in Engineering* 1997; **40**(5):817–832.

6. Chen S, Wang Z, Liu XH. Active vibration control and suppression for intelligent structures. *Journal of Sound and Vibration* 1997; **200**(2):167–177.

7. Saravanos DA, Heyliger PR, Hopkins DA. Layerwise mechanics and finite element for the dynamic analysis of piezoelectric composite plates. *International Journal of Solids and Structures* 1997; **34**(3):359–378.

8. Trindade MA, Benjeddou A, Ohayon R. Finite element modelling of hybrid active-passive vibration damping of multilayer piezoelectric sandwich beams—part I: formulation. *International Journal for Numerical Methods in Engineering* 2001; **51**(7):835–854.

9. Lammering R, Mesecke-Rischmann S. Multi-field variational formulations and related finite elements for piezoelectric shells. *Smart Materials and Structures* 2003; **12**(6):904–913.

10. Carrera E, Boscolo M. Classical and mixed finite elements for static and dynamic analysis of piezoelectric plates. *International Journal for Numerical Methods in Engineering* 2007; **70**(10):1135–1181.

11. Saravanos DA, Heyliger PR. Mechanics and computational models for laminated piezoelectric beams, plates, and shell. *Applied Mechanics Review* 1999; **52**(10):305–320.

12. Benjeddou A. Advances in piezoelectric finite element modeling of adaptive structural elements: a survey. *Computers and Structures* 2000; **76**(1–3):347–363.

13. Mackerle J. Smart materials and structures—a finite element approach—an addendum: a bibliography (1997–2002). *Modelling and Simulation in Materials Science and Engineering* 2003; **11**(5):707–744.
14. Hagood NW, Flotow AV. Damping of structural vibrations with piezoelectric materials and passive electrical networks. *Journal of Sound and Vibration* 1991; **146**(2):243–268.

15. Lesieutre GA. Vibration damping and control using shunted piezoelectric materials. *The Shock and Vibration Digest* 1998; **30**(3):187–195.

16. Nguyen CH, Pietrzko SJ. FE analysis of a PZT-actuated adaptive beam with vibration damping using a parallel R-L shunt circuit. *Finite Elements in Analysis and Design* 2006; **42**(14):1231–1239.

17. Richard C, Guyomar D, Audigier D, Bassaler H. Enhanced semi-passive damping using continuous switching of a piezoelectric device on an inductor. *Smart Structures and Materials: Passive Damping and Isolation, SPIE* 2000; **3989**:288–299.

18. Corr LR, Clark WW. A novel semi-active multi-modal vibration control law for a piezoceramic actuator. *Journal of Vibration and Acoustics* 2003; **125**(2):214–222.

19. Makihara K, Onoda J, Minesugi K. Low-energy-consumption hybrid vibration suppression based on an energy-recycling approach. *AIAA Journal* 2005; **43**(8):1706–1715.

20. Badel A, Lagache M, Guyomar D, Lefeuvre E, Richard C. Finite element and simple lumped modeling for flexural nonlinear semi-passive damping. *Journal of Intelligent Material Systems and Structures* 2007; **18**(7):727–742.

21. Collinger JC, Wickert JA. Adaptive piezoelectric vibration control with synchronized switching. *Proceedings of ASME International Mechanical Engineering Congress and Exposition, IMECE 2007-41427*, Seattle, Washington, U.S.A., 2007.

22. Davis CL, Lesieutre GA. A modal strain energy approach to the prediction of resistively shunted piezoceramic damping. *Journal of Sound and Vibration* 1995; **184**(1):129–139.

23. Becker J, Fein O, Maess M, Gaul L. Finite element-based analysis of shunted piezoelectric structures for vibration damping. *Computers and Structures* 2006; **84**(31–32):2340–2350.

24. Ducarne J, Thomas O, Deü JF. Optimisation de dispositif passif d'atténuation de vibration par shunt piézoélectrique. *Actes du 8ème Colloque National en Calcul de Structures*, vol. 2. Hermes: Giens, France, 2007; 519–524 (in French).

25. Ducarne J, Thomas O, Deü JF. Performance and optimization of piezoelectric shunts for vibration reduction. *Journal of Sound and Vibration* 2009; 2009; in preparation.

26. Ducarne J, Thomas O, Deü JF. Structural vibration reduction optimization by switch shunting of piezoelectric elements. *Proceedings of ASME International Mechanical Engineering Congress and Exposition, IMECE, 2007-43396*, Seattle, Washington, U.S.A., 2007.

27. Maurini C, Pouget J, dell'Isola F. Extension of the Euler-Bernoulli model of piezoelectric laminates to include 3D effects via a mixed approach. *Computers and Structures* 2006; **4**(22–23):1438–1458.

28. ANSYS, Inc. Coupling and constraint equations. *Basic Analysis Procedures Guide* (5.5 edn), Chapter 12. Notes on CP command, 1998. Available from: http://uic.edu/depts/accc/software/ansys/html/guide_55/g-mod/GMOD12.htm.

29. Piefort V, Preumont A. Finite element modeling of piezoelectric structures. *Samtech User's Conference*, Paris, France, 2001. Available from: http://www.ulb.ac.be/scmero/piezo.html.

30. Johnk CTA. *Engineering Electromagnetic Fields and Waves* (2nd edn). CRC Press: Boca Raton, 1988.

31. ANSI/IEEE Standard 176-1987. *IEEE Standard on Piezoelectricity*, 1988.

32. Tiersten H. *Linear Piezoelectric Plate Vibration*. Plenum Press: New York, 1969.

33. Deü JF, Larbi W, Ohayon R. Piezoelectric structural acoustic problems: symmetric variational formulations and finite element results. *Computer Methods in Applied Mechanics and Engineering* 2008; **197**(19–20):1715–1724.

34. Rahmoune M, Benjeddou A, Ohayon R. New thin piezoelectric plate models. *Journal of Intelligent Material Systems and Structures* 1998; **9**(12):1017–1029.

35. Benjeddou A, Deü JF. A two-dimensional closed-form solution for the free-vibrations analysis of piezoelectric sandwich plates. *International Journal of Solids and Structures* 2002; **39**(6):1463–1486.

36. Galucio AC, Deü JF, Ohayon R. A fractional derivative viscoelastic model for hybrid active–passive damping treatments in time domain—application to sandwich beams. *Journal of Intelligent Material Systems and Structures* 2005; **16**(1):33–45.

37. Porfiri M, Maurini C, Pouget J. Identification of electromechanical modal parameters of linear piezoelectric structures. *Smart Materials and Structures* 2007; **16**(2):323–331.

38. Reddy JN. *Mechanics of Laminated Composite Plates and Shells: Theory and Analysis*. Wiley: New York, 1997.

39. Crawley EF, de Luis J. Use of piezoelectric actuators as elements of intelligent structures. *AIAA Journal* 1987; **25**(10):1373–1385.

40. Sonti VR, Kim SJ, Jones JD. Equivalent forces and wavenumber spectra of shaped piezoelectric actuators. *Journal of Sound and Vibration* 1995; **187**(1):111–131.
41. Polytec PI. *Piezoelectric Ceramic Products*, 2008. Available from: http://www.piceramic.de/pdf/KATALOG_english.pdf.

42. Ikeda T. *Fundamentals of Piezoelectricity*. Oxford University Press: Oxford, 1990.

43. Thomas O, Touzé C, Chaigne A. Asymmetric non-linear forced vibrations of free-edge circular plates, part 2: experiments. *Journal of Sound and Vibration* 2003; **265**(5):1075–1101.

44. Sénéchal A, Thomas O, Deü JF, Jean P. Atténuation de vibrations de structures complexes par élément piézoélectrique. application une aube de turbomachine. *Actes du 9ème Colloque National en Calcul de Structures*, Giens, France, 2009 (in French).