Simulation of wave propagation in plate structures by using new spectral element with piezoelectric coupling

1. Introduction

It is well known that detecting small damage in structures is an important and challenging task [1-2]. Recently, the use of Lamb wave based technology in structural health monitoring has made significant progress [3-5]. It is observed that the behavior of Lamb wave is very complicated. For example, the wave occurs in multiple modes which can convert into each other under special conditions during propagation [6]. These phenomena make the Lamb wave based monitoring technology more complicated in order to detect and locate the failures. Therefore, it is important to understand thoroughly the complex behavior of the Lamb wave interacting with different kinds of damage in structures to reduce the cost and time in the design process of such structural health monitoring systems. Besides the experimental investigations, numerical simulations will definitely play an important role in design of efficient Lamb wave based monitoring systems.

Although different higher order finite element (FE) schemes are available for simulations of Lamb waves [7], however, the time domain spectral element method is a more promising method for simulating wave propagation in complex structures [8-18]. The time domain spectral element method combines the accuracy of the global pseudospectral method with the flexibility of the local finite element method thus can accurately simulate the wave propagation in structures in terms of both phase and amplitude.

For real Lamb wave based structural health monitoring systems, the wave is usually generated by piezoelectric (PZT) actuators and picked up by PZT sensors. Therefore, the piezoelectric coupling effect should be considered in the numerical simulations. Currently most numerical simulations of Lamb wave propagations [8-18] have not considered the coupling effect to reduce the modeling effort. Only fewer papers [19, 20] have considered the piezoelectric coupling effect. In doing so there-dimensional (3-D) spectral elements [19] have been used to take into account the piezoelectric coupling effect, however, the computing time and memory may become impractical for real complex structures. To reduce the computational resource and simplify the problem, a new approach is proposed to incorporate the coupling effect within the spectral element framework [20]. Although the method reported in [20] seems feasible, it may still suffer some shortages. For example, the piezoelectric effect in the thickness direction is not considered, a non-diagonal mass matrix may be resulted since PZT elements are usually mounted only on one face of the host plate in practice. Therefore, it is still necessary to investigate other efficient ways.

For plate structures, it is desirable to use a two dimensional (2D) plate element to reduce the modeling effort as well as the computational resources. In this paper, a new time domain 2D spectral element is proposed to incorporate the coupling of piezoelectric transducers. Detailed formulations are presented. To validate the formulations, a benchmark problem presented in [19] is analyzed. Numerical results are compared with data obtained by ABAQUS with fine meshes. Conclusions are drawn based on the results reported herein.

2. Formulations of the new spectral plate element

2.1. Displacement and strain fields

To model the 3-D behavior of the Lamb wave propagation in plates by using a 2-D spectral plate element, seven variables, $u^T$, $u^B$, $v^T$, $v^B$, $w^T$, $w^M$, $w^B$, are introduced. $u$, $v$, $w$ are displacement components in the x, y, and z direction and superscripts $B$, $M$, $T$ denote the bottom, middle and top surfaces, respectively. The reason to choose these variables will be demonstrated in Section 3. The proposed method is also applicable if more variables are required. Fig. 1 shows a typical 6×6-node line (or simply node) 2-D spectral plate element with (Fig. 1a) and without (Fig. 1b) a PZT layer. It is seen that the modeling effort will be reduced since only 2-D elements are involved.

In the global $x$-$y$-$z$ coordinate system set in the middle plane of the host plate shown in Fig. 1, the 3-D displacement field of the plate can be written as:

1

$u\left(x,y,z\right)=\left(0.5-\frac{z}{h}\right)u^B+\left(0.5+\frac{z}{h}\right)u^T,$
$v\left(x,y,z\right)=\left(0.5-\frac{z}{h}\right)v^B+\left(0.5+\frac{z}{h}\right)v^T,$
$\left(x,y,z\right)=\left(\frac{2z^2}{h^2}-\frac{z}{h}\right)w^B+\left(1-\frac{4z^2}{h^2}\right)w^M+\left(\frac{2z^2}{h^2}+\frac{z}{h}\right)w^T,$

where $h$ is the thickness of host plate. Each node, denoted by the filled points in Fig. 1, has seven degrees of freedom $q_k(k=\mathrm{1,2},...,7)$, i.e., $u^T({\xi }_i,{\eta }_j)$, $u^B({\xi }_i,{\eta }_j)$, $v^T({\xi }_i,{\eta }_j)$, $v^B({\xi }_i,{\eta }_j)$, $w^T({\xi }_i,{\eta }_j)$, $w^M({\xi }_i,{\eta }_j)$, $w^B({\xi }_i,{\eta }_j)$. Each node in the PZT layer, denoted by the hollow circles shown in Fig. 1a, has only one degree of freedom $\varphi$, the electrical potential. Note that the electrical potential at the interface is assumed to be zero identically.

Fig. 1Sketches of spectral plate element: a) with a PZT layer, b) without a PZT layer

a)

b)

For small strains, the 3-D strain-displacement relations can be written by:

2

${\epsilon }_x=\left(\frac{1}{2}-\frac{z}{h}\right)\frac{\partial u^B}{\partial x}+\left(\frac{1}{2}+\frac{z}{h}\right)\frac{\partial u^T}{\partial x},$
${\epsilon }_y=\left(\frac{1}{2}-\frac{z}{h}\right)\frac{\partial v^B}{\partial y}+\left(\frac{1}{2}+\frac{z}{h}\right)\frac{\partial v^T}{\partial y},$
${\epsilon }_z=\left(\frac{4z}{h^2}-\frac{1}{h}\right)w^B-\frac{8z}{h^2}{w}^M+\left(\frac{4z}{h^2}+\frac{1}{h}\right)w^T,$
${\gamma }_{yz}=\left(\frac{2z^2}{h^2}-\frac{z}{h}\right)\frac{\partial w^B}{\partial y}+\left(1-\frac{4z^2}{h^2}\right)\frac{\partial w^M}{\partial y}+\left(\frac{2z^2}{h^2}+\frac{z}{h}\right)\frac{\partial w^T}{\partial y}-\frac{v^B}{h}+\frac{v^T}{h},$
${\gamma }_{zx}=\left(\frac{2z^2}{h^2}-\frac{z}{h}\right)\frac{\partial w^B}{\partial x}+\left(1-\frac{4z^2}{h^2}\right)\frac{\partial w^M}{\partial x}+\left(\frac{2z^2}{h^2}+\frac{z}{h}\right)\frac{\partial w^T}{\partial x}-\frac{u^B}{h}+\frac{u^T}{h},$
${\gamma }_{xy}=\left(\frac{1}{2}-\frac{z}{h}\right)\frac{\partial u^B}{\partial y}+\left(\frac{1}{2}+\frac{z}{h}\right)\frac{\partial u^T}{\partial y}+\left(\frac{1}{2}-\frac{z}{h}\right)\frac{\partial v^B}{\partial x}+\left(\frac{1}{2}+\frac{z}{h}\right)\frac{\partial v^T}{\partial x},$

or simply:

3

$\left\{\gamma \right\}=\left[\frac{\partial }{\partial x}\right]\left\{q\right\},$

where $\left\{q\right\}={\left[u^T,u^B,v^T,v^B,w^T,w^M,w^B\right]}^T$.

For an $N\times N$-node 2-D spectral element, the displacement field in the element can be assumed as:

4

$q_k\left(\xi ,\eta \right)=\sum _{i=1}^N\sum _{j=1}^N{N}_i\left(\xi \right)N_j\left(\eta \right)q_k\left({\xi }_i,{\eta }_j\right)=\sum _{i=1}^N\sum _{j=1}^N{N}_{ij}\left(\xi ,\eta \right)q_k\left({\xi }_i,{\eta }_j\right),(k=\mathrm{1,2},...,7),$

where $N_{ij}(\xi ,\eta )=N_i\left(\xi \right)N_j\left(\eta \right)$ are shape functions, $q_k\left({\xi }_i,{\eta }_j\right),(i,j=\mathrm{1,2},...,N,k=\mathrm{1,2},...,7)$ denote the nodal degrees of freedom, and $N_i\left(\xi \right)$ and $N_j\left(\eta \right)$ are the one-dimensional Lagrange interpolation functions defined by:

5

$N_i\left(\xi \right)=\prod _{k=1,k\ne i}^N\frac{\xi -{\xi }_k}{{\xi }_i-{\xi }_k};N_j\left(\eta \right)=\prod _{k=1,k\ne j}^N\frac{\eta -{\eta }_k}{{\eta }_j-{\eta }_k},(i,j=\mathrm{1,2},...,N),$

where the distribution of nodes, ${\xi }_i$, ${\eta }_i$, can be the Gauss–Lobatto–Legendre (GLL) points, or the Chebyshev points, or the approximate Lebesgue-optimal grid points, i.e., the expanded-Chebyshev grid points [21-22]. It is shown [18] that the critical time step to ensure stable time integration by using the central finite difference method is the largest if the expanded-Chebyshev grid points are used as the nodes of the spectral element.

Let $z^{*}\in [-h_1/2,h_1/2]$, where $h_1$ is the thickness of the PZT layer. Since the electric potential on the lower surface of the PZT layer, connected with the host plate, is zero identically, thus the electric voltage $V$ in the PZT layer is given by:

6

$V=\sum _{i=1}^4{H}_{Vi}{\varphi }_i,$

where $H_{Vi}$ is the shape functions and defined by:

7

$H_{Vi}=\frac{1}{8}\left(1+{\xi }_i\xi \right)\left(1+{\eta }_i\eta \right)\left(1+\frac{2z^{*}}{h_1}\right),\left(i=1...4\right).$

3. Results and discussion

To verify the proposed formulations, a benchmark problem presented in [19] is analyzed. A thin aluminum plate with two surface-mounted piezoelectric transducers shown in Fig. 2 is considered. The dimension of the plate and the two piezoelectric transducers are 508 mm×∞ mm×1.02 mm and 6.35 mm×∞ mm 0.25 mm. The poling direction is parallel to the z direction. One piezoelectric transducer on the left side is used as the actuator and the other as the sensor. A uniform voltage is input on the piezoelectric element’s upper surface. Since this problem can be regarded as a plane strain problem, thus the displacement in the y direction is set to zero identically. This simplifies the problem for investigation of the effect of spatial resolution in numerical methods without loss of generality [19].

Fig. 3 shows the mesh used in the analysis. Fig. 3(a) is the model by using proposed 6×6 spectral plate element with seven degrees of freedom (DOFs) per node. GLL nodes are used. Fig. 3(b) is the mesh used in the analysis by ABAQUS. In the numerical simulations, the excitation signal is a five-peak tone-burst wave with a center frequency of 100 kHz and kept a peak voltage at 50 V, shown in Fig. 4. The material properties of the host plate are: elasticity modulus $E=$ 70GPa, mass density $\rho =$ 2700 kg/m³, and Poisson’s ratio $v=$ 0.3. The material properties of the PZT element (Navy Type II), provided by APC company, are listed in Table 1.

Fig. 2Geometry of an aluminum plate

Fig. 3Meshes by (a) proposed spectral plate element and (b) finite plain strain element

a)

b)

Fig. 4Excitation signal

The sensing signal obtained by the present element is shown in Fig. 5. The data obtained by ABAQUS with plane strain elements are also included in the figure for comparisons. It is seen that the solution obtained by using the proposed spectral element is quantitatively in accord with data obtained by the FEM (ABAQUS). The comparison verifies the formulation of the proposed element and the developed program. Although the solutions match quite well for the $S_0$ mode, however, a small discrepancy is seen in the slower $A_0$ mode. It is pointed out [19] that due to the numerical dispersion error, errors in phase/group velocities exist for both FEM and spectral element method (SPE). For FEM, the solution accuracy would be further improved if smaller size of element than the one shown in Fig. 3(b) is used. For the SEM, an exponential convergence rate is observed and results are more accurate [19]. Due to restriction of the displacement field adopted in the proposed element, however, numerical dispersion error may be larger than the existing 3-D spectral element [19] in the numerical simulation of Lamb wave propagations. Thus, the accuracy of the solution obtained by the proposed element is similar to the result obtained by the FEM. Further increase in the order of the displacement field along $z$ direction or refinement of the mesh will definitely reduce the numerical dispersion error, but the computational efficiency will also be decreased.

Fig. 5Comparisons of the output voltage with data by ABAQUS

Fig. 6Comparisons of the output voltage with data by ABAQUS

Fig. 7Comparisons of the output voltages obtained by different elements

To investigate the effect of different variation of displacement w along the thickness direction, $w$ is assumed a linear function of $z$, the same as the other two displacement components $u$ and $v$ shown in Eq. (1). Therefore, a six-variable or six-degree of freedom per node spectral plate element can be established. The element is similar to the spectral plate element presented in [18]. Following a similar way, a 6×6-node with six DOFs per node spectral plate element considering the piezoelectric effect can be derived.

The same problem is then re-analyzed by using the 6×6-node with six DOFs per node spectral plate element. Again GLL nodes are used and the same mesh shown in Fig. 3a is adopted in the analysis. The sensing signal obtained by the spectral plate element is shown in Fig. 6. The data obtained by ABAQUS is also included in the figure for comparisons. It is seen that although the agreement is quite well for $S_0$ mode, the magnitude of $A_0$ mode is much smaller than the data obtained by ABAQUS. It suggests that higher order variation of displacement in the thickness direction is necessary. This is the main reason why Eq. (1) is used as the displacement field for the proposed spectral plate element.

Fig. 7 shows the comparisons of output signals obtained by the 6×6-node spectral plate elements with six DOFs or seven DOFs per node. It is seen that if the piezoelectric effect in the poling ($z$) direction is not considered, namely, $E_{33}=$ 0, the magnitudes of both $A_0$ mode and $S_0$ mode will be reduced greatly. In other words, keeping the piezoelectric effect in the poling ($z$) direction is important to obtain correct amplitude of the sensing signals by PZT sensors.

4. Conclusions

In this paper, a new spectral plate element is proposed to take into account the coupling effect of piezoelectric actuators and sensors. A simple way to formulate the diagonal mass matrix is used to reduce the computational effort. Detailed formulations are given. A benchmark problem is studied. It is shown that without considering the coupling effect of piezoelectric actuators and sensors in the thickness direction may result much smaller amplitude of signals picked up by PZT sensors. The simulated results by the proposed element are well compared to data obtained by ABAQUS with fine meshes. In modeling Lamb wave propagation in plate structures, the proposed element has advantages over the three-dimensinal spectral element in terms of modeling effort, computational time and memory space.

It should be pointed out that the proposed method can also be used for variables higher than seven. In other words, the displacement fields $u$, $v$, $w$ along the thickness direction can be assumed even more complicated if additional variables are introduced to increase the accuracy. Besides, variables can be located at any position along the thickness direction, since nodes for a 2-D spectral plate element are actually nodal lines. In this way, the numbering effort will be greatly reduced in establishing a spectral element model for plate structures.

In the present paper, the dynamic interface interaction between the PZT actuator/sensor and the host structure as well as the dynamic behavior of the adhesive interface has not been considered. Further researches are necessary and underway to study these effects.