Abstract

A three-dimensional spectral element method (SEM) was developed for analysis of Lamb wave propagation in composite laminates containing a delamination. SEM is more efficient in simulating wave propagation in structures than conventional finite element method (FEM) because of its unique diagonal form of the mass matrix. Three types of composite laminates, namely, unidirectional-ply laminates, cross-ply laminates, and angle-ply laminates are modeled using three-dimensional spectral finite elements. Wave propagation characteristics in intact composite laminates are investigated, and the effectiveness of the method is validated by comparison of the simulation results with analytical solutions based on transfer matrix method. Different Lamb wave mode interactions with delamination are evaluated, and it is demonstrated that symmetric Lamb wave mode may be insensitive to delamination at certain interfaces of laminates while the antisymmetric mode is more suited for identification of delamination in composite structures.

1. Introduction

Owing to its superior mechanical properties and light weight, composite materials are finding more and more applications especially in aerospace industries [1]. However, composite structures still run a high risk of suffering from damage due to abrupt impact or growth of fatigue defects, which can result in catastrophic failure during their service life. It is, therefore, essential to develop techniques to inspect integrity and improve safety, reliability, and operational life of structures [27].

Traditionally, nondestructive evaluation (NDE) techniques, such as C-scan and radiographic inspection, are used to evaluate the integrity and degradation of structures on a periodic basis. Now online structural health monitoring (SHM) techniques, for example, vibration-based and Lamb-wave-based techniques, are being developed to provide early warning and assure the performance of structures. Among them, Lamb-wave-based techniques for damage detection have received a considerable attention in the past decades due to its ability of long-distance propagation and sensitivity to a variety of defects.

For Lamb waves-based damage detection techniques, understanding wave propagation characteristics in structures is essential for their successful application. Therefore, a number of numerical methods have been applied to analyze propagation of elastic waves, such as finite difference method (FDM) [8, 9], finite element method (FEM) [1012], boundary element method (BEM) [13, 14], finite strip elements (FSE) [15, 16], mass-spring lattice model (MSLM) [1721], and local interaction simulation approach (LISA) [22, 23].

In the literature, two kinds of spectral element method (SEM) were applied to wave propagation modeling, namely, fast-Fourier-transform- (FFT-) based SEM and the orthogonal polynomials-based SEM [2426]. In the FFT-based SEM [24], a single element is sufficient to model wave propagation in large uniform structures, which is suited for simple 1D and 2D problems [27]. On the other hand, the orthogonal-polynomials- (e.g., Legendre and Cheybysev polynomials) based SEM [25] is much more suitable for analyzing wave propagation in structures with complex geometry. The formulation of the spectral finite element (SFE) is similar to FE when assembly of element matrices and solution of equations are considered. Hence, the SEM can be used to handle the wave propagation in structures with complex geometry, and various types of defects can be modeled. The difference between the SEM and FEM comes in that orthogonal polynomials are used as approximation functions in SEM and, therefore, calculation can be more efficient because of the diagonal mass matrix. This method has been successfully applied to many fields, such as fluids, seismology and acoustics [28, 29]. More recently, the SEM was used to simulate wave propagation in structures for damage detection, for example, wave propagation in 1D structures, such as rod and beam [30, 31]. Numerical simulation results of the elastic wave propagation in a composite plate were presented by Kudela et al. [32]. A 2D spectral membrane finite element-based model, developed by Zak et al. [33], was used to analyze wave propagation in a cracked isotropic panel. Wave propagation in 2D plate structures using a 3D SEM for damage detection was also discussed by Peng et al. [34].

A wave propagation analysis in composite structures using 3D SEM for the purpose of damage detection has not been widely reported in the literature so far. Although characteristics of Lamb wave propagation in composite laminates and the damage evaluation using numerical simulation have been rigorously explored for a couple of years, in most of the related work, structures were simplified by either one- or two-dimensional models, resulting in approximate results especially for laminates of complicated layup or with damage. The SEM combines accuracy with flexibility in describing problems with complex geometries, which is highly desirable for modeling of elastic wave propagation. In this paper, multilayered composite laminates are modeled using the Legendre polynomials-based spectral finite element, elastic wave propagation characteristics are analyzed, and wave interaction with delamination is discussed.

2. Wave Propagation in Composite Plate

Composite laminates are commonly fabricated by stacking unidirectional lamina with a certain layup configuration. After a composite is properly cured, a multilayered structure is formed with all the layers bonded together. For analysis of wave propagation, each lamina can be regarded as a quasihomogeneous orthotropic or transversely isotropic material with the main principal axis parallel to the fibres. Therefore, in numerical calculations, a composite laminate is modeled as a multilayered medium with different elastic and anisotropic material properties [35].

In order to form a numerical model of composite structure, a Cartesian coordinate system is firstly defined by the 𝑧-axis normal to the central plane of a composite laminate spanned by the 𝑥-axis and the 𝑦-axis for modeling of the wave propagation, as shown in Figure 1.

For each layer of the composite laminate, the stress-strain relations for arbitrary direction have the following form [36]:𝝈=𝐷𝜺,(2.1) where 𝝈 is stress vector, 𝜺 is strain vector, and𝐷 is the flexibility matrix. In order to study wave propagation, the elastic constants of all the layers must be expressed in the global coordinate system. For those layers, the principal material coordinate system does not coincide with the global coordinate system; this can be achieved by using a transformation matrix method.

In a homogeneous media, the elastic wave propagation is described by the governing equation [37]: 𝜌𝜕2𝑡𝐮=𝜎+𝐟,(2.2) where 𝐮 is the displacement vector. 𝜌 is mass density and 𝐟is external force vector.

3. Formulation of 3D Spectral Element Method

For the Legendre polynomials-based 3D spectral finite element, it requires that the domain Ω in three dimensions is decomposed into a number of nonoverlapping hexahedrons,Ω𝑒, as in the conventional FE method. In SEM, the equations of wave propagation is [37] Ω𝜌𝐰𝜕2𝑡𝐮𝑑Ω+Ω𝐰𝐂𝐮𝑑Ω=Ω𝐰𝐟𝑑Ω,(3.1) where the Ω denotes the physical region of interest and 𝐰is an arbitrary test vector.

In SEM, the nodes are defined into two steps: (1) each element in its physical domain is mapped to a reference domain Λ=[1,1]3 using an invertible local mapping 𝑓; a set of basis functions consisting of Legendre polynomials of degree 𝑁 are introduced; (2) a set of nodes are defined as 𝜉𝑖[1,1], 𝑖1,,𝑁+1. These Gauss-Lobatto-Legendre (GLL) points are the (𝑁+1) roots of [29]1𝜉2𝑃𝑁(𝜉)=0,(3.2) where 𝑃𝑁(𝜉) is the derivative of the Legendre polynomial of degree 𝑁. The 𝜉𝑖 are different from the classical FE method in which the nodes are uniformly spaced. As an example, a 108-node spectral element in the local coordinate is shown in Figure 2.

The Lagrange interpolation function,𝑢𝑒𝑁, supported by the GLL points can be expressed as𝑢𝑒𝑁(𝜉,𝜂,𝛾)=𝑛1𝑛𝑚=12𝑛𝑛=13𝑟=1𝑢𝑒𝑁𝜉𝑚,𝜂𝑛,𝛾𝑟𝑚(𝜉)𝑛(𝜂)𝑟=(𝛾)𝑛1𝑛𝑚=12𝑛𝑛=13𝑟=1𝑢𝑒𝑁𝜉𝑚,𝜂𝑛,𝛾𝑟Ψ𝑚𝑛𝑟,(3.3) where Ψ𝑚𝑛𝑟 is defined as the orthogonal shape functions in 3-D 𝑚(𝜉) denotes the 𝑚th 1D Lagrange interpolation at the (𝑁+1) GLL points defined above. The property of 𝑚(𝜉) is𝑚𝜉𝑛=𝛿𝑚𝑛,(3.4) where 𝛿𝑚𝑛 denotes the Kronecker delta and 𝑛𝑖, 𝑖=1,2,3, are the numbers of GLL points in each direction in the local coordinate.

Therefore, the element matrices,𝐌𝑒 (mass matrix), 𝐊𝑒 (stiffness matrix), and 𝐅𝑒 (force vectors), are calculated numerically on the reference coordinate:𝐌𝑒=𝜌Ω𝑒𝚿𝑒(𝑥,𝑦,𝑧)𝑇𝚿𝑒(𝑥,𝑦,𝑧)𝑑Ω𝑒=𝜌+11+11+11𝚿𝑒(𝜉,𝜂,𝛾)𝑇𝚿𝑒𝐉(𝜉,𝜂,𝛾)det𝑒𝑑𝜉𝑑𝜂𝑑𝛾=𝜌𝑛1𝑖=1𝜔𝑖𝑛2𝑗=1𝜔𝑗𝑛3𝑘=1𝜔𝑘𝚿𝑒𝜉𝑖,𝜂𝑗,𝛾𝑘𝑇𝚿𝑒𝜉𝑖,𝜂𝑗,𝛾𝑘𝐉det𝑒,𝐊𝑒=Ω𝑒(𝐁𝑒)(𝑥,𝑦,𝑧)𝑇𝐃𝑒(𝐁𝑒(𝑥,𝑦,𝑧))𝑑Ω𝑒=+11+11+11(𝐁𝑒(𝜉,𝜂,𝛾))𝑇𝐃𝑒(𝐁𝑒𝐉(𝜉,𝜂,𝛾))det𝑒=𝑑𝜉𝑑𝜂𝑑𝛾𝑛1𝑖=1𝜔𝑖𝑛2𝑗=1𝜔𝑗𝑛3𝑘=1𝜔𝑘𝐁𝑒𝜉𝑖,𝜂𝑗,𝛾𝑘𝑇𝐃𝑒𝐁𝑒𝜉𝑖,𝜂𝑗,𝛾𝑘𝐉det𝑒,𝐅𝑒=Ω𝑒𝚿𝑒(𝑥,𝑦,𝑧)𝑇𝐏𝑑Ω𝑒=+11+11+11𝚿𝑒𝐉(𝜉,𝜂,𝛾)𝐏(𝜉,𝜂,𝛾)det𝑒=𝑑𝜉𝑑𝜂𝑑𝛾𝑛1𝑖=1𝜔𝑖𝑛2𝑗=1𝜔𝑗𝑛3𝑘=1𝜔𝑘𝚿𝑒𝜉𝑖,𝜂𝑗,𝛾𝑘𝑇𝐏𝜉𝑖,𝜂𝑗,𝛾𝑘𝐉det𝑒,(3.5) where 𝜌is the mass density, 𝐃𝑒is termed material stiffness matrix, and 𝐏 is a distributed load. Ψe are the shape functions based on the Legendre polynomials. The matrix 𝐁𝑒 is the strain-displacement matrix calculated by 𝐁𝑒=𝐋𝚿𝑒(𝑥,𝑦,𝑧),(3.6) where 𝐋 denotes a differential operator matrix:𝜕𝐋=𝑥00𝜕𝑦0𝜕𝑧0𝜕𝑦0𝜕𝑥𝜕𝑧000𝜕𝑧0𝜕𝑦𝜕𝑥𝑇.(3.7)𝐉𝑒is the Jacobian matrix associated with the mapping 𝑓 from the physical domain Ωe to the reference domainΛ:𝐉𝑒=𝜕𝜉𝑥𝜕𝜉𝑦𝜕𝜉𝑧𝜕𝜂𝑥𝜕𝜂𝑦𝜕𝜂𝑧𝜕𝛾𝑥𝜕𝛾𝑦𝜕𝛾𝑧.(3.8) The weights𝜔𝑖 are defined by𝜔𝑖=2𝑛𝑃(𝑛1)𝑛1𝜉𝑖2,𝑖1,,𝑛,𝑛=𝑁+1.(3.9) Therefore, wave equation (3.1) can be rewritten into matrix form and wave propagation modeling is transformed to an ordinary differential equation in time. Let 𝐔 denotes the global vector of unknown displacement in the medium. Then the ordinary differential equation can be written as𝐌̈𝐔+𝐊𝐔=𝐅,(3.10) where 𝐌 denotes the global mass matrix, 𝐊 is the global stiffness matrix, and 𝐅 is the vector of time-dependent excitation force.

In this study, the differential equation (3.10) is solved using a central difference time integration scheme. Initial conditions of displacement and velocity are 𝐔=0 anḋ𝐔=0 at the time 𝑡=0, and the central difference time integration scheme is implemented as1Δ𝑡2𝐌𝐔𝑡+Δ𝑡=𝐅𝑡2𝐊Δ𝑡2𝐌𝐔𝑡1Δ𝑡2𝐌𝐔𝑡Δ𝑡,(3.11) where the symbol 𝑡 denotes time and Δ𝑡 denotes the time step of integration. In the central difference time integration scheme, the stable condition is Δ𝑡Δ𝑡𝑐𝑟=𝐿/𝑐, where 𝐿 is the minimum distance between two adjacent nodes and 𝑐 is the wave speed in elastic medium.

In comparison with FE, less computation effort is required for SFE because of the choice of Lagrange interpolation supported on the GLL points in conjunction with the GLL integration rule. The efficiency of SEM was demonstrated using two 3D models based on SFEs and FEs with the same degrees of freedom, and a reduction of about 65% in CPU time for calculation is used in comparison with the FEM.

4. Numerical Calculation

Lamb wave propagation in 8-ply T300/F593 composite laminates is analyzed in this study. The material properties of unidirectional lamina are listed in Table 1. Three types of composite laminates, namely, unidirectional [08], symmetric cross-ply [02/902]𝑠, and quasi-isotropic [+45/45/0/90]s laminates are investigated. All the laminates have the same geometric configurations of 500 mm × 500 mm × 1.72 mm, as shown in Figure 3. The composite plate was meshed using the three-dimensional spectral finite elements with 108-node (shown in Figure 2).

In the literature, piezoelectric (PZT) wafer is one of the frequently applied transducers to excite and capture Lamb waves propagating in structures for damage detection. The PZT transducers can be modeled by adding force to the incident points or using electromechanical coupling for elastic wave modeling in structures [38, 39]. In the present study, in order to obtain relatively simple Lamb wave mode, two forces perpendicular to the plane of the plate are applied at point A in Figure 3 on the upper and the lower surfaces of the plate, respectively. When the two forces are in phase, antisymmetric modes are activated, and when the two forces are out of phase, symmetric modes are excited. The excitation force is a 5-cycle sinusoidal signal modulated by Hanning window with a center frequency of 100 kHz and absolute maximum magnitude is 1𝑁, as shown in Figure 4. The displacement responses at point B and C on the upper surface of the composite laminate are used to investigate wave propagation characteristics.

4.1. Wave Propagation in a Multilayered Composite Plate

As aforementioned, three types of the composite laminates are analyzed in this study. Firstly, Lamb wave propagation in the unidirectional composite laminate [0]8 is investigated. The plate is meshed to 50 × 50 × 1 spectral elements. Under the two excitation of in-phase forces, the fundamental antisymmetric mode A0 is excited. According to the simulation results, the displacement component in the 𝑧-direction is dominant. Therefore, only the displacement component in the 𝑧-direction is analyzed here. Responses of the laminate at 0.097 ms are plotted in Figure 5(a). It can be observed that the wave front of the A0 mode has ellipse-like shape, and the group velocity in the fibre direction is greater than that in the direction perpendicular to the fibre, indicating that group velocity of this type of wave mode in the composite depends on the orientation of wave propagation. Generally, group velocity in composite laminate is a function of the direction of wave propagation and direction of the fibres. The group velocity in composite laminates can be calculated analytically using the transfer matrix method (In the appendix). For such an 8-ply unidirectional laminate, the analytical group velocity is plotted in a polar coordinate, as shown in Figure 5(b). It presents similar feature as Figure 5(a). Group velocities in the direction along the fibers (𝑥-axis) and perpendicular to the fibres (𝑦-axis) are further calculated from displacement responses of SEM simulations at the point B and C, as shown in Figure 6. Based on the peak of the received responses, time of flight (ToF) from point A to B can be defined. The calculated group velocity 𝑐𝑔 of the A0 mode of Lamb waves propagating in the 𝑥 direction is 1794 m/s, which is 2.6% smaller than the analytical value. In a similar way, the group velocity of the A0 mode in the 𝑦-direction is 1319 m/s and the one calculated analytically is 1245 m/s, which gives a relative error of 5.9%. It can be concluded that the simulation results of SEM model, the proposed model, are in good agreement with the analytical results thus validating the effectiveness of the model.

Under the two excitations of out-of-phase forces, the fundamental symmetric wave modes, the S0 mode and the SH0 mode, are excited. According to the simulation results, the displacement components in the 𝑥-direction and in the 𝑦-direction are dominated. Hence, only those two displacement components are plotted, as shown in Figures 7(a) and 7(b). It can be seen that S0 mode and the SH0 mode are excited simultaneously, which is because the S0 mode and the SH0 mode are coupled in multilayered composite laminate. Analytical group velocities of the S0 mode and the SH0 mode are also calculated, as shown in Figure 7(c), and the simulation results of SEM modeling are in good agreement of the analytical results.

In case of cross-ply laminates [02/902]𝑠, the laminate was meshed to 50 × 50 × 4 spectral finite elements. The symmetric mode and antisymmetric modes are excited using the above-mentioned method. In case of antisymmetricmode, the displacement component in the 𝑧-direction is shown in Figure 8. On the other hand, in case of symmetric mode the displacement in the 𝑥-direction and the 𝑦-direction are plotted in Figure 9. The complexity of Lamb wave propagation in composite is also demonstrated, and the effectiveness of the proposed model is validated by comparison of simulation results and the analytical results.

In case of angle-ply laminates [45/45/0/90]𝑠, the composite laminate is meshed by using 50 × 50 × 8 spectral finite elements. The displacement components in the 𝑧-direction at 0.086 ms are plotted in Figure 10(a), while displacement components in the 𝑥-direction and in the 𝑦-direction at 0.069 ms are plotted in Figures 11(a) and 11(b). The group velocities of analytical solutions for those Lamb modes are also plotted in a polar coordinate for comparison, as shown in Figures 10(b) and 11(c). A good agreement found in both symmetric and antisymmetricmodes. It can be observed that the angular dependence of Lamb nodes in the laminates [45/45/0/90]𝑠 becomes weaker because of its quasi-isotropic layup. Under the S0 and the SH0 modes, the group velocities are approximately independent of direction of wave propagation, but it still can be discerned that the A0 mode has the maximum in the 45° (or 225°) directions because outer lamina are orientated in these directions, which dominates the bending properties related to A0 mode.

In case of angle-ply laminate, the wavelength of the A0 mode is shorter than that of the S0 and the SH0 modes in composite laminate from the simulation results and analytical results. It can be expected that the A0 mode is possibly more sensitive to small damage than the S0 mode and the SH0 mode because of its short wavelength.

It is demonstrated that characteristics of the wave propagation in multilayered composite is complex due to the nature of anisotropic of the material constants and the multilayered configurations, which leads to the group velocity of Lamb waves depending on the laminate layup and the direction of wave propagation. Good agreement between the simulation results based on SEM and the analytical results demonstrates that the proposed model provides an effective tool to investigate the wave propagation in composite structures.

4.2. Wave Propagation in a Composite Plate Containing Delamination

The typical damage forms in composite laminate are transverse microcracking, fiber-breakage, and delamination. Typically, the transverse microcracking through the thickness of the ply occurs as the first-ply failure, and then delamination damage follows. The fiber breakage usually happens at the last stage of the failure. However, a catastrophic failure can occur only with the microcracking and delamination damage without the fiber breakage. Delamination is known to happen because of excessive interlamina normal and shear stress at the ply boundaries, which not only causes reduction in stiffness, but also affects the strength and integrity of the structure, leading to failure.

The wave propagation in composite laminates containing a delamination is investigated in this study. The delamination in the laminates is modeled using nodes separation method. Laminate without delamination is initially meshed. At the interface between the adjacent elements where the delamination occurs, nodes that are affected by the delamination are separated, as shown in Figure 12.

Lamb wave modes interaction with delamination is analyzed in quasi-isotropic laminates [45/−45/0/90]s using the proposed 3D SEM model. The size of the delamination is 30 mm × 30 mm in a square shape, as shown in Figure 3. The wave interaction of the symmetric mode and antisymmetricmode with delamination is investigated.

The effects of the delamination at different interfaces in the composite laminate are addressed. Under the symmetric mode, the displacement responses at the point B in the 𝑥- and the 𝑦-directions are plotted in Figure 13. Responses of the intact composite laminate are also provided for comparison. It is evident that the scattered waves from the delamination are not so obvious in comparison with the intact laminate, on the captured responses at the point B, when the delamination is located in different interfaces. The displacement responses at 0.069 ms in the 𝑥- and the 𝑦-directions of the composite laminates are plotted in Figure 14, when there is a delamination at the interfaces between 3-4 layers. However, when the delamination is at the interface between 4-5 layers, responses of the point B are same as those from intact composite laminate, indicating that the symmetric mode is insensitive to the delamination in the symmetric plane of the plate, attributed to the reason that the layup of the composite laminate is symmetric about the central plane between 4-5 layers. Therefore, the symmetric wave mode related to extensional wave mode can travel through the delamination without any scattering from it.

Under the antisymmetricmode, the displacements at the point B in the 𝑧-direction are plotted in Figure 15. It can be seen that the antisymmetricmode is more sensitive to delamination located at all interfaces of the composite laminate. The displacement responses of the laminate in the 𝑧-direction at 0.173 ms are plotted in Figure 16.

Under the same excitation frequency, the scattered waves from delamination become clearer in the captured response at the point B under the antisymmetricmode than that under symmetric mode. Hence, the antisymmetricmode can be used to detect smaller delamination than the symmetric mode, since the wavelength of the antisymmetricmode is less than that of the symmetric mode. In addition, according to the simulation results, the symmetric mode may be less sensitive to delamination located at certain interfaces. It can be expected that the antisymmetricmode is more suitable for identification of the delamination.

The scattered waves are very weak compared with the incident wave, resulting in that the reflected wave packet is difficult to identify delamination in such composite materials with a high attenuation ratio. Fortunately, the transmitted wave has been affected a lot, as shown in Figure 16.

5. Conclusions

A three-dimensional spectral element method is developed to investigate the wave propagation characteristics in composite laminates in the present study. Three types of laminates, namely, unidirectional [0]8, cross-ply [02/902]𝑠, and quasi-isotropic [45/45/0/90]𝑠 laminates are modeled using 3D spectral elements, and the laminates are excited using two forces in-phase and out-of-phase to generate the symmetric mode and antisymmetricmode respectively. It is demonstrated that the proposed 3D spectral element method can be efficiently and effectively used to simulate wave propagation in composite laminates. Complexity of wave propagation characteristic is also demonstrated even for a single Lamb wave mode in composite laminates. Finally, interactions between the Lamb wave mode and a delamination are analyzed. It is concluded that symmetric mode of Lamb wave may be insensitive to delamination in certain interfaces in the laminate. And, therefore, it is essential to understand wave propagation characteristics in composite laminate when Lamb-wave-based structural health monitoring strategy is carried out.

Appendix

The analytical wave front can be calculated as follows.

Consider wave propagation solutions in the following form:{𝑢,𝑣,𝑤}=(1,𝑉,𝑊)𝑈𝑒[𝑖𝜉(𝑥1+𝛼𝑥3)𝜔𝑡],(A.1) where 𝜉is the wave number, 𝜔is the circular frequency, 𝛼is still an unknown parameter, and 𝑉and 𝑊 are ratios of the displacement amplitudes of 𝑣 and 𝑤, respectively. The choice of the solution leads to the three coupled equations that can be written as𝐾11𝐾12𝐾13𝐾21𝐾22𝐾23𝐾31𝐾32𝐾33𝑈1𝑈2𝑈3=0,(A.2) where the elements of the matrix 𝐾 are𝐾11=𝐶11+𝐶55𝜌𝑐2,𝐾12=𝐶16+𝐶45𝛼2,𝐾13=𝐶13+𝐶55𝐾𝛼,22=𝐶66𝜌𝑐2+𝐶44𝛼2,𝐾23=𝐶36+𝐶45𝐾𝛼,33=𝐶55𝜌𝑐2+𝐶33𝛼2.(A.3) The existence of nontrivial solutions for 𝑈1, 𝑈2, and 𝑈3 demands the vanishing of the determinant of the matrix 𝐾 and yields the sixth-degree polynomial equation:𝛼6+𝐴1𝛼4+𝐴2𝛼2+𝐴3=0.(A.4) There are six roots of this equation, which correspond to the three sets of mode pairs. For each 𝛼𝑞, 𝑞=1,2,,6, the displacement ratios𝑉𝑞=𝑈2𝑞/𝑈1𝑞 and 𝑊𝑞=𝑈3𝑞/𝑈1𝑞 can be expressed as𝑉𝑞=𝐾11𝛼𝑞𝐾23𝛼𝑞𝐾13𝛼𝑞𝐾12𝛼𝑞𝐾13𝛼𝑞𝐾22𝛼𝑞𝐾12𝛼𝑞𝐾23𝛼𝑞,𝑊𝑞=𝐾11𝛼𝑞𝐾23𝛼𝑞𝐾12𝛼𝑞𝐾13𝛼𝑞𝐾12𝛼𝑞𝐾33𝛼𝑞𝐾23𝛼𝑞𝐾13𝛼𝑞.(A.5) The formal solutions for the displacements and stresses in the expanded matrix form𝑢1𝑢2𝑢3𝜎33𝜎13𝜎23=𝑉1111111𝑉1𝑉3𝑉3𝑉5𝑉5𝑊1𝑊1𝑊3𝑊3𝑊5𝑊5𝐷11𝐷11𝐷13𝐷13𝐷15𝐷15𝐷21𝐷21𝐷23𝐷23𝐷25𝐷25𝐷31𝐷31𝐷33𝐷33𝐷35𝐷35𝑈11𝐸1𝑈12𝐸2𝑈13𝐸3𝑈14𝐸4𝑈15𝐸5𝑈16𝐸6,(A.6) where𝐸𝑞=𝑒𝑖𝜉𝛼𝑞𝑥3,𝐷1𝑞𝐶=𝑖𝜉13+𝐶36𝑉𝑞+𝐶33𝛼𝑞𝑊𝑞,𝐷2𝑞𝐶=𝑖𝜉55𝛼𝑞+𝑊𝑞+𝐶45𝛼𝑞𝑉𝑞,𝐷3𝑞𝐶=𝑖𝜉45𝛼𝑞+𝑊𝑞+𝐶44𝛼𝑞𝑉𝑞,𝑞=1,2,,6.(A.7)

Acknowledgments

The authors are grateful for the support received from the National Natural Science Foundation of China (NSFC nos. 11072148 and 11061160491), Research Project of State Key Laboratory of Mechanical System and Vibration (MSV201110), and the National High Technology Research and Development Program of China (no. 2009AA044800).