Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 49, 3, pp. 879-903, Warsaw 2011 SIMULATION OF WAVE PROPAGATION IN DAMPED COMPOSITE STRUCTURES WITH PIEZOELECTRIC COUPLING Rolf T. Schulte Claus-Peter Fritzen University of Siegen, Institute of Mechanics and Automatic Control – Mechatronics, Siegen, Germany; e-mail: fritzen@imr.mb.uni-siegen.de This paper presents an efficient and accurate simulation approach to shorten time and cost of the necessary pre-tests in the design process of structural health monitoring (SHM) systems. The simulation is perfor- med using the time domain spectral elementmethod, which leads to an optimally concentratedmassmatrix and results in a crucial reduction of complexity of the time integration algorithm. The theoretical backgro- und of the method and a spectral element for flat shells are presented. Newapproaches to incorporate the anisotropicmaterial damping and an efficient coupling of piezoelectric elements within the spectral element framework are developed. Some numerical calculations are performed showing both the accura- cy of this methodology, by comparing to experimental values, and the applicability to more complex structures like stiffened curved panels. Key words: spectral element method, composites, anisotropic attenu- ation 1. Introduction During the last decades,many structuralhealthmonitoring (SHM)approaches based on activating and sensing elastic waves using piezoelectric patches have been developed (Balageas et al., 2006; Giurgiutiu, 2007). Themajority of the- se systems have shown their capabilities in the laboratory environment, and their technology readiness level (TRL) is in the range of two to five, compare Mankis (1995). So far, there are few examples of the transfer to more com- plex real world structures. The reason for this can be found inmany obstacles 880 R.T. Schulte, C.-P. Fritzen arising from different sectors: complexity of the technology, manufacturing, certification, regulation authorities and somemore. This paper concentrates on the first of these points, the complexity of the technology itself. Many wave propagation phenomena and especially the in- teractions with different kinds of damage are difficult to understand. This is particularly true with regard to received sensor data because of the presence of different wave modes that often cannot be distinguished. However, only a profound comprehension of wave propagation leads to the ability to configure optimized SHM systems for a desired application. The system performance strongly depends on an adequate choice of parameters like excitation signals, damage evaluation algorithms, actuator/sensor types and actuator/sensor di- stribution. For successful monitoring of complex real world structures, the influence of all these parameters has to be studied in detail. For this reason there is a strong interest in accurate simulation and visualisation tools to give insight into the details of thewavemotion. Precise and efficient numericalmo- delling techniques can help to cut time and cost of the parameter optimization process, and thereby can assist the challenging step to real world structures and higher TRLs. Generally, for the modelling of wave propagation phenomena a variety of methods can be utilised. Well-established are the finite difference method (FDM) (Graves, 1996), the pseudospectral method (PSM) (Fornberg, 1987; Ku et al., 1987), the finite element method (FEM) (Moser et al., 1999; Zien- kiewicz et al., 2005) and the boundary elementmethod (BEM) (Beskos, 1997; Cho and Rose, 1996). Because of the high frequency signals, which are requ- ired to detect small damage, a very dense finite element mesh is necessary for accurate simulation of wave propagation and possible wave scattering at defects. For that reason, the conventional finite element simulation becomes computationally very inefficient. Finite difference methods suffer from nume- rical dispersion and difficulties arise when implementing boundary conditions (Graves, 1996). Other approaches use mass-spring-lattice models (MSLM), see e.g. Yim and Sohn (2000). Developments in this field include the so called local interaction simulation approach (LISA), in which wave propagation is simulated without using finite difference equations, but directly from physical phenomena and properties (Delsanto et al., 1997). A quite different methodology is the utilisation of so called spectral ele- ments. Using this terminology, it is important to clearly distinguish between two different techniques knownunder the same term.TheFourier transforma- tionbased spectral elements, seeDoyle (1997), start fromexact solutions of the partial differential equations in the frequency domain. The excitation signals Simulation of wave propagation... 881 are transferred to the frequency domain using Fast Fourier Transformation (FFT) and after solving the inverse FFT is used to calculate the response in the time domain. This technique shows very good results for one-dimensional structures (beams, rods) but difficulties arise in handling more complex geo- metries, because exact solutions of the problem are often not available. Amore promisingmethod for complex structures is the time domain spec- tral elementmethod (SEM) that was first proposed by Patera (1984). It com- bines the accuracy of the global pseudospectral method with the flexibility of the FEM.However, the spectral elementmethod is notwidespread in the con- text of wave propagation in thin-walled structures. InKudela et al.(2007) and Kudela and Ostachowicz (2009) spectral elements are used to simulate trans- verse wave propagation in undamped composite plates. In Ostachowicz(2008) different elements (rods, beams, plates) based on the spectral elementmethod including damage models are derived. In Peng et al. (2009) waves in planar shells are calculated bymeans of spectral volume elements. The SEM is used in recent publications to investigate the wave propagation and interaction with damages in complex structures like a welded L-joint (Rucka, 2010) or a half-pipe structure (Ostachowicz andKudela, 2010). In this paper, the formulation of a flat shell spectral element is presen- ted, including consideration of material damping. The construction of mass-, stiffness- anddampingmatrices, andthe time integration algorithmarepresen- tedbriefly.Amodelof the electromechanical couplingofpiezoelectric actuators and sensors is adapted and integrated into the spectral element framework. By construction, this leads to an efficient tool to simulate wave propagation phenomena in thin-walled structures. An example of propagating waves in a unidirectional composite plate is presented including the comparison of simu- lation results with measured data. Moreover, a model of a curved panel with stiffeners is analysed to demonstrate the applicability tomore complex-shaped geometries. The simulation results of the presented model can be used to de- velop advanced signal processing strategies and to performdetailed parameter studies. 2. Nodal base and shape functions in the spectral element method For the numerical treatment of high-frequency wave propagation one impor- tant aspect to be kept inmind is the necessity of a densemesh-grid to account for the shortwavelength. This requires to consider the interpolation properties 882 R.T. Schulte, C.-P. Fritzen of the used element type to be able to capture the wave motion with as few degrees of freedom as possible. Theoretical analysis shows that the highest interpolation accuracy is achieved by distributing the interior nodes of an ele- ment corresponding to the roots of certain orthogonal polynomials, namely Chebyshev or Lobatto polynomials (Pozrikidis, 2005). In contrast to the Che- byshev polynomials, the N +1 nodes of the Lobatto nodal base, the roots of equation (2.1) (1− ξ2)LoN−1(ξ) (2.1) in conjunction with Lagrange interpolation polynomials leads by construc- tion to an optimally concentrated mass matrix as will be seen subsequently. LoN−1(ξ) denotes the Lobatto polynomial of (N −1)-th order. The Lagrange interpolation polynomials fulfil the discrete orthogonality ψi(ξj) = δij where δij denotes the Kronecker delta. In Fig.1 these interpolation polynomials – also called shape functions – are depicted within the one-dimensional referen- ce domain for a nodal base of 7 Gauss-Lobatto-Legendre (GLL) nodes. The discrete orthogonality can be clearly seen as each polynomial is exactly one at its correspondingGLL node and zero at all other nodes. Fig. 1. Lagrange interpolation polynomials on the reference interval [−1,1] through 7 Gauss-Lobatto-Legendre (GLL) points Moreover, the excellent interpolation properties of this kind of elements can be seen in Fig.1, because the shape functions vary smoothly between approximately −0.2 and 1 and any kind of Runge oscillations, that occur at the ends of a higher order element with evenly spaced nodes when the polynomial order is raised, are totally suppressed. For the extension to two-dimensional elements, as the flat shell elements presented within this paper, quadrilateral elements based on a tensorial pro- duct of the one-dimensional constituents should be favoured over triangular elements. The tensorial product guarantees that the discrete orthogonality Simulation of wave propagation... 883 also holds for the two-dimensional shape functions. If triangular elements ha- ve to be used, Fekete points should be used as the nodal base (Komatitsch et al., 2001; Pozrikidis, 2005), but with respect to accuracy and numerical sta- bility quadrangles are preferable, as it is stated in Komatitsch et al. (2001). According to Seriani, only about 5 to 6 nodes (exact value depends on the de- gree of the interpolation polynomial) per shortest wavelength are necessary to correctly simulate the wave propagation behaviour in contrast to 15-30 which are needed using lower order finite elements (Seriani and Priolo, 1994). 3. Flat shell spectral element for composite material The flat shell spectral element presented here is based on the first order shear deformation theory (FSDT) for plates that was developed by Mindlin. The three independent displacements ũ, ṽ and w̃ are expressed using a displace- ment field of the form ũ(x,y,z,t) = u(x,y,t)+ ∂z0 ∂x w(x,y,z,t)+zθy(x,y,t) ṽ(x,y,z,t) = v(x,y,t)+ ∂z0 ∂y w(x,y,z,t)−zθx(x,y,t) (3.1) w̃(x,y,z,t) =w(x,y,z,t) where (u,v,w,θx,θy) are unknown functions to be determined. u, v and w are the displacement of the plane z = 0, θx and θy denote right-hand-rule rotations, see Fig.2. Fig. 2. Kinematic parameters used for the description of a flat shell element, the dashed line is the midface of the shell Arotational dof θz about the local z-axis is not usedwithin the kinematic description. This inhibits arbitrary transformations of the elements in spa- ce between local and global coordinates, because general transformations in 884 R.T. Schulte, C.-P. Fritzen space account for 6 dofs. In principle, it is possible to use an in-plane formula- tion with explicit θz-dof but this would complicate the problemunnecessarily. Instead, an artificial θz-dof that is not used within the kinematic formula- tion, can be introduced. To regularise the resulting system of equations, small artificial values for the stiffness andmass of this dof are defined following Ba- the (1986). The terms (∂z0/∂x)w(x,y,z,t) and (∂z0/∂y)w(x,y,z,t) take into account the curvature of the shell in form of a pre-deformation, see Wagner (1985). For planar shells, these terms are dropped. The displacement field on the element level can be approximated in the following form    w(ξ,η) θx(ξ,η) θy(ξ,η) u(ξ,η) v(ξ,η)    = N+1∑ i=1 N+1∑ j=1 Ψij(ξ,η)q (e) ij = N+1∑ i=1 N+1∑ j=1 ψi(ξ)ψj(η)    ŵ(ξi,ηj) θ̂x(ξi,ηj) θ̂y(ξi,ηj) û(ξi,ηj) v̂(ξi,ηj)    (3.2) where ξ and η denote the local element coordinates, Ψ is a two-dimensional shape function and ψi is the i-th one-dimensional shape function defined above. q(e) represents the vector of nodal variables and the hat indica- tes nodal degrees of freedom that are the unknowns of the resulting sys- tem of equations. Examples of the GLL nodal distribution within the two- dimensional reference element are depicted in Fig.3 for a 25-node and a 81-node element. Raising the degree of the interpolation polynomial – that is increasing the number of nodes per element – leads to clustering of nodes at the element corners. The implication of this effect is discussed subsequently. Fig. 3. Nodal distribution within the reference element for 5 GLL nodes (left) and 9 GLL nodes (right) per element edge Simulation of wave propagation... 885 Some examples of the two-dimensional shape functions Ψ on a grid using 25 element nodes are given in Fig.4. Fig. 4. Examples of shape functions Ψij of a spectral-shellelement with 5GLL-nodes per element edge The derivation of the variational equation and the weak form follows the standard FE procedures, see e.g. Hughes (1987). The element matrices are calculated using the Gauss-Lobatto integration rule. For the stiffness matrix this leads to K (e) = ∫∫ Ωe [B(x,y)]⊤DB(x,y)det(J) dΩ (3.3) ≈ N+1∑ i=1 N+1∑ j=1 wiwj[B(xij,yij)] ⊤ DB(xij,yij)det(J) 886 R.T. Schulte, C.-P. Fritzen where w are quadrature weights and D is the material stiffness matrix that is defined as D=   D11 D12 D16 0 0 B11 B12 B16 D12 D22 D26 0 0 B12 B22 B26 D16 D26 D66 0 0 B16 B26 B66 0 0 0 κA55 κA45 0 0 0 0 0 0 κA45 κA44 0 0 0 B11 B12 B16 0 0 A11 A12 A16 B12 B22 B26 0 0 A12 A22 A26 B16 B26 B66 0 0 A16 A26 A66   (3.4) where κ is a shear correction factor. Themass matrix can be calculated as M (e) = ∫∫ Ωe [Ψ(x,y)]⊤HΨ(x,y)det(J) dΩ (3.5) ≈ N+1∑ i=1 N+1∑ j=1 wiwj[Ψ(xij,yij)] ⊤ HΨ(xij,yij)det(J) where H is given in the following equation H=   I0 0 0 0 0 0 I2 0 0 −I1 0 0 I2 I1 0 0 0 I1 I0 0 0 −I1 0 0 I0   (3.6) Themass moments of inertia Ii are given according to   I0 I1 I2  = ∫ he   1 z z2  ρ(z) dz (3.7) where ρ is the material density. Themassmatrix is optimally concentrated by construction because of the aforementioned discrete orthogonality of the shape functions. For a density distribution that is symmetric towards themidplane of the shell, the terms I1 in thematrix H vanish. Therefore, a completely diagonal massmatrix results in that case. Both Ke and Me can be assembled in the samemanner as it is done with the conventional finite elements to form global stiffness andmass matrices K and M. The spectral elements naturally lead to free-free boundaries, other boundary conditions can be implemented as for low-order finite elements. Simulation of wave propagation... 887 4. Incorporation of material damping For a realistic simulation model of a composite material, the influence of at- tenuation has to be considered. While in other papers this point is neglected for simplicity (e.g. Kudela and Ostachowicz, 2009) the incorporation of dam- ping into the spectral element framework for composite shells is an important aspect of this contribution. The viscous dampingmodel that is often applied in finite element analysis of dynamics of structures, assumes the equation of motion in the following form Mq̈+Cq̇+Kq=F(t) (4.1) The dampingmatrix C is usually assumed proportional to themass and stiff- ness matrix C= αM+βK (4.2) where the scalar coefficients α and β are determined from experimental me- asurements. Besides thesemodels, some other approaches are used, especially the application of fractional derivatives seems promising (Bagley and Torvik, 1985). Within this contribution, viscous damping is assumed leading to an equ- ation of motion in form of Eq. (4.1). Indeed, in contrast to Eq. (4.2), the damping matrix is assembled in a different manner. As can be noticed from experiments on composite plates, thedampingbehaviour often shows an inver- se proportionality to the in-plane stiffness. For that reason, the second term in Eq. (4.2), the linear couplingwith the stiffnessmatrix, is completely dropped. However, the proportionality to themassmatrix alone is not able to take into account the direction-depending attenuation behaviour. Instead, the followingmaterial dampingmatrix is introduced for each layer C (k) mat =   1 2 (C (k) θf +C (k) θm ) 0 0 0 0 0 C (k) θm 0 0 0 0 0 C (k) θf 0 0 0 0 0 C (k) f 0 0 0 0 0 C (k) m   (4.3) providing different coefficients for the different degrees of freedom. Index f denotes the fibre direction, index m the transverse direction. This matrix can be formulated for each of the material layers separately in the material coordinate system, where the first axis is parallel to the fibre direction, see Fig.5. 888 R.T. Schulte, C.-P. Fritzen Fig. 5. Laminate coordinate system (x,y,z) and k-th layer with its material coordinate system (x(k),y(k),z(k)) Thematerial dampingmatrices fromEq. (4.3) can be transformed into the laminate coordinate system by using C (k) mat =Tc(−φ (k))C (k) matT ⊤ c (−φ (k)) (4.4) where the matrix Tc is an appropriate rotation matrix. Finally, the material dampingmatrix for the laminate Cmat = ∫ he Cmat dz = L∑ k=1 zk+1∫ zk C (k) mat dz (4.5) and the element dampingmatrix can be constructed C (e) = ∫∫ Ωe [Ψ(xij,yij)] ⊤ CmatΨ(xij,yij)det(J) dΩ (4.6) ≈ N+1∑ i=1 N+1∑ j=1 wiwj[Ψ(xij,yij)] ⊤ CmatΨ(xij,yij)det(J) If this approach is used on the element level, different attenuation coefficients can be used for the fibre direction and the transverse direction. With this methodology, two important properties are achieved: on the one hand, the direction-dependent attenuation behaviour can be realised and the damping coefficients of in-planeandout-of-planewaves canbechosen separately.Onthe other hand, the resulting dampingmatrix is diagonal. This property preserves the opportunity tomake use of a rapid time integration scheme, whichwill be explained during the next section. Simulation of wave propagation... 889 5. Time integration of the global system The equation of motion of the global system (4.1) can be discretised in time using either implicit or explicit time integration schemes. To fully utilise the advantage of the diagonal mass and damping matrices, the explicit central difference scheme is used here. Themain assumption of the central difference scheme is an acceleration of the form q̈ t = 1 ∆t2 (qt−∆t−2qt+qt+∆t) (5.1) The velocity can be assumed as q̇ t = 1 2∆t (−qt−∆t+qt+∆t) (5.2) Inserting (5.1) and (5.2) into the systemof ordinarydifferential equations (4.1) at the time instant t results in ( 1 ∆t2 M+ 1 2∆t C ) q t+∆t =Ft− ( K− 2 ∆t2 M ) q t − ( 1 ∆t2 M− 1 2∆t C ) q t−∆t (5.3) that can be solved for qt+∆t. The solution of Eq. (5.3) can be performed very efficiently if matrix inversions can be avoided, which is the case if M and C are diagonal matrices. This methodology allows rapid calculation of the time dependent solution of the system. If thesematrices are not diagonal, the term in brackets on the left hand side of Eq. (5.3) has to be decomposed, which in- creases the computational cost drastically. If instead of theGLL interpolation, Chebyshevpolynomials are used, themassmatrix is no longer diagonal. In this case, either an iterative time-stepping algorithm could be used (Seriani, 2004) or the mass- and damping matrices can be diagonalised by row-summing the contributions on the diagonal (Dauksher and Emert, 2000). The drawback of this scheme is the difficulty to set the timestep ∆t that has to be appropriately small to reach stability of the solution. It is very difficult to derive a generally applicable formula to compute this timestep for wave propagation problems. It depends onmany parameters like themaximal velocity of the travelling wave vmax, the excitation frequency, the distortion of elements andminimumdistance ofmesh-nodes ∆xmin. The optimal length of the time increment ∆t has to be adjusted to the particular problem of interest. It can be estimated by using the Courant-Friedrichs-Levy condition ∆tmax =CFL ∆xmin vmax (5.4) 890 R.T. Schulte, C.-P. Fritzen CFL is the Courant-number that is close to 0.85 for the one-dimensional case, but can only be estimated for more complex situations. The positive aspect of this particular problem is the fact that a numerical error resulting from a too large timestep ∆t makes the solution grow very fast to infinity. This behaviour enables the user to distinguish easily between the correct results and numerically distorted results by the too large timestep. AsEq. (5.4) indicates a direct relationship of the timestep to theminimum grid-spacing, this equation dictates theuseful range ofGLLnodesper element- edge.Asdiscussedabove anddepicted inFig.2, increasing theorder of spectral elements leads to clustering of nodes at the element edges and reduces ∆xmin. While on the one hand higher order elements lead to a better approximation, on theother hand the length of the time incrementhas tobe reduced (denoting more timesteps for the same period T , resulting in a higher computational cost). From a practical point of view, the elements with an approximation order of 5 to 9 seem an acceptable compromise. 6. Piezoelectric coupling Asmentioned before, an efficient coupling of piezoelectric patches is included into the spectral element framework.On the one hand these patches have local contributions to themass, stiffness anddampingproperties of the structure, on the other hand the coupling of electrical andmechanical quantities due to the piezoelectric effect has to be considered. Thewell-known piezoelectric effect is not to be explained in detail here because the general equations can be found in many textbooks, see e.g. Ikeda (1990). Assuming isotropic piezoelectric actuators, themagnitude of the induced line forces andmoments at the edges of a rectangular piezo actuator patch that is bondedonto the surface of a plate structure (actuator equation) can be expressed, following Banks et al. (1996), as Npztx = N pzt y = −Epzthpzt 1−νpzt d31 hpzt Vpzt (6.1) Mpztx =M pzt y =− 1 8 Epzt 1−νpzt [ 4 (h 2 +hpzt )2 −h2 ] d31 hpzt Vpzt Here, h denotes the plate thickness, hpzt is the thickness of the piezo element, d31 is a piezo-ceramic strain constant, Epzt and νpzt areYoung’smodulus and the Poisson’s ratio of the piezo patch, and Vpzt denotes the applied voltage. If Simulation of wave propagation... 891 instead of rectangular patches, disc-shaped pzt-elements are used, the forces and moments act at those nodes nearest to the physical edges of the pzt- element, as it is depicted in Fig.6. The line forces of the quadratic-shaped element (indicated in grey) are replaced by forces at those GLL nodes nearest to the radius of the round element (indicated in black). Mass, stiffness and damping contributions to the global system are also only considered for those areas physically covered by the pzt-element. Fig. 6. Alignment of equivalent piezoelectric nodal forces (black arrows) for a pzt-disc The sensor equation is developed for orthotropic situations in Yang and Ngoi (1999), and is used to calculate the output charge Q(t) Q(t)=− ∫∫ S (F +εX33E3) dxdy +zk ∫∫ S G dxdy (6.2) with F(x,y)= e31 ∂u0 ∂x +e32 ∂v0 ∂y G(x,y)= e31 ∂ϕy ∂x +e32 ∂ϕx ∂y (6.3) The charge depends on the strain in the piezo patch, as can be seen in equation (6.3), where e31 and e32 are the coupling coefficients. For circular patches, only the proportion of the element that is physically covered by the patch is considered for the strain calculation. Because of the high density of the interior mesh, the resulting output charge is accurate even for a small wavelength corresponding to high frequency waves. Using these equations is an efficient possibility to incorporate the electro- mechanical coupling because it avoids adding further degrees of freedom to the system, as has to be done using an approach including electrical degrees of 892 R.T. Schulte, C.-P. Fritzen freedom, see e.g. Lammering andMesecke-Rischmann (2003). It is straight for- ward to compute the generalised actuator forces and add themas external for- ces Fext to the global system. The strains for the sensing equation (Eq. (6.2)) can be easily calculated either analytically using the strain-displacement ma- trix B or numerically from the resulting displacements. Effects of shear-lag due to the glue-layer between piezo patch and host structure can easily be im- plemented by loss factors, as it is done in the application shown subsequently. It is also possible to consider more refinedmodels of the piezoelectric patches including eigenresonances and so on, but this is not in the scope of this paper and not necessary for the considered frequency range used here. 7. Application I, wave propagation in a GFRP plate The presentedmethodology is used to simulate the propagation of waves in a unidirectional GFRP plate of dimensions 800mm×800mm that is depicted in Fig.7. The thickness of the plate is approximately 1.5mm, the fibre direction is in the y-direction. As can be noticed in Fig.7, nine piezoelectric patches are attached to the structure. The piezo-discs have a diameter of 10mm and a thickness of 0.25mm. For the simulation, the plate is meshed with 68×68 elements using 36 nodes per element. Fig. 7. GFRP plate with piezo patches To demonstrate the capabilities of this simulation approach, the patch located in the centre (P5) is used as an actuator, and the other patches are Simulation of wave propagation... 893 used as sensors. A voltage pulse signal with a carrier frequency of 100kHz lasting for 5 cycles (windowed by an Hann-window) is used as the actuation signal, and the received sensor signals at the other patches are bothmeasured and calculated. In Fig.8, themeasured and simulated sensor voltage signals are compared for patchesP2,P3 andP6 corresponding to 90◦, 45◦ and 0◦. For this compari- son, a loss factor of 0.78 is introduced to compensate for the shear lag effect of the actuator and sensors, andYoung’smoduli, shearmoduli anddensity of the plate are optimized. The signals that can be seen at themeasurement voltage of P2 andP3 at the time before 0.05ms are an electrical coupling phenomenon with the actuator channel, because the contacts are not perfectly shielded. In general, Fig.8 shows a good agreement of the simulated and measured data including the time of flight of different wave-modes, amplitude and phase. Fig. 8. Simulated andmeasured sensor signals for P2, P3 and P6; excitation at 100kHz 894 R.T. Schulte, C.-P. Fritzen Comparison of the signals of P2 and P6 (Fig.8) shows the anisotropy of this plate, because both sensors have the same distance to the actuator, but the time of flight is significantly shorter to P2 that is placed in the direction of the fibres. The anisotropic characteristic becomes even more obvious by considering snapshots of thepropagatingwaves. InFigs.9 and10 someof these snapshots are depicted. As the excitation signal for these plots, a windowed pulse signal with a carrier frequency of 75kHz lasting for four cycles was used. In Fig.9 the resulting in-plane displacement, calculated from the x- and y- displacement is depicted (in-plane waves), while in Fig.10 the z-displacement is shown(out-of-planewaves). InFig.9, the samescale is used for all snapshots, because the damping of the in-plane waves is not as strong as for the out-of- planewaves, where different scales had to be used (see Fig.10). In both figures small disturbances of the wave field can be noticed, when the travelling waves reach the piezo sensors (Fig.9c,d and Fig.10b,c). In contrast to isotropic plates, where the SH0 mode is usually not exci- ted by bonded piezo actuators, in Fig.9 two different wave modes, namely S0 and SH0 can be clearly identified, because besides the S0 wave, there is another group of waves with a lower wave speed and wavelength. The SH0 wave propagates predominantly in a direction close to ±45◦ (see Fig.9d,e). In Figs.10c and 10d, the anisotropic damping characteristics can be noti- ced, because travelling thewave packet in the x-direction looses its amplitude much faster than in the y-direction. Moreover, the dispersive character of the out-of plane waves is obvious, because the pulse widens with increasing runtime, in contrast to the in-plane waves that are barely dispersive. In comparison to pure sensor data, the presented simulation methodology provides a deeper insight tomanydetails of thewave propagation phenomena. While on a single plate the behaviour of the travelling waves is comparatively simple to predict, matters get more complicated, if more complex structu- res are under investigation. This makes it particularly attractive to use this simulation approach for the analysis of more complex structures. 8. Application II, wave propagation in a stringer stiffened curved panel In the aerospace industry, typical sub-structures are curved panels with stif- feners. For that reason the presented simulation approach is demonstrated subsequently on suchapanelwith a radius of 2m,which is typical formedium- sized aircraft. The geometry of the panel with its element-mesh is shown in Simulation of wave propagation... 895 Fig. 9. Snapshots of propagating in-plane-waves at 75kHz in the GFPR plate; S0 and SH0 wave can be identified 896 R.T. Schulte, C.-P. Fritzen Fig. 10. Snapshots of propagating out-of-plane-waves at 75kHz in the GFPR plate Fig. 11. Curved panel with stiffeners and piezo patches Fig.11, where piezo elements are indicated in black. The excitation signal is a windowed pulse signal with a carrier frequency of 50kHz lasting for three cycles. Simulation of wave propagation... 897 The T-shape stiffeners are modelled close to reality as it can bee seen in Fig.12. For this model, an isotropic aluminium material is assumed. The connection of the base-shell to the stiffener is modelled using a five times higher attenuation factor to account for higher damping because of gluing. Fig. 12. Detailed view of the stiffener geometry The subsequent figures show snapshots of different waves travelling within the panel at different time instances. The out-of-plane waves excited by the piezopatch mainly remain in the area between the two neighbouring stiffeners (see Figs. 13 to 15). Fig. 13. y-displacement at t =0.10ms 898 R.T. Schulte, C.-P. Fritzen Fig. 14. y-displacement at t =0.20ms Fig. 15. y-displacement at t =0.30ms The in-planewaves can leave this areamuch easier, as can be seen in Figs. 16 and 17, where the resulting in plane displacement is depicted. The in-plane waves within the base-shell excites out-of plane waves in the stiffeners, as can be seen in Fig.17 from the shorter wavelength in the stiffeners. As in Fig.17 some contributions of waveswith a shorterwavelength also in the base-shell are noticeable, in Fig.18 the local in-plane displacements at the same time instance are plotted for comparison. This makes clear that those contributions arise just fromthe curvature of thebase-shell, because the global xz-displacement contains parts from the out-of-plane wave as well. Simulation of wave propagation... 899 Fig. 16. xz-displacement at t =0.10ms Fig. 17. xz-displacement at t =0.30ms Fig. 18. Local in-plane displacement at t =0.30ms 900 R.T. Schulte, C.-P. Fritzen 9. Conclusions The formulation of flat spectral shell elements has been presented in detail. The novelty of this contribution is, in particular, the capability of the presen- ted formulation to include direction dependent attenuationwithin the spectral element framework. The construction of the dampingmatrix and the efficient implementation within the central difference scheme for the time integration were discussed. Moreover, piezoelectric coupling is implemented within the spectral element framework. This is done using analytical equations and avo- idingadditional degrees of freedom.Theability to handledampingphenomena broadens the scope of application of this methodology towards realistic simu- lation models especially regarding composite materials. The proposed approach is demonstrated by two examples including a com- parison to experimental data for a unidirectional GFRP plate. Simulated and measured sensor data show a good correlation and affirm the necessity to con- sider the attenuation. The second, pure numerical example demonstrates the application to a stiffened curved panel, amore complex structure that is often used as sub-structure in the aircraft industry. The simulation can support a deeper understanding of the propagation behaviour of different wavemodes including direction-depending attenuation, reflections from various obstacles as the stiffeners, and somemore. In the next step, the presented approach can be used to optimize actu- ator/sensornetworks on structures includingextensiveparameter studies.This process can help to shorten time, cost and material consumption of pretests and may help to bring SHM a small step further to real life application and higher TRL levels. References 1. BagleyR.L., TorvikP.J., 1985,Fractional calculus in the transient analysis of viscoelastically damped structures,AIAA Journal, 23, 6, 918-925 2. Balageas D., Fritzen C.-P., Güemes A. (eds.), 2006, Structural Health Monitoring, Hermes Science Publishing, London UK 3. Banks H.T., Smith R.C., Wang Y., 1996, Smart Material Structures, Mo- delling, Estimation and Control, Wiley, Masson, Paris Simulation of wave propagation... 901 4. Bathe K.-J., 1986, Finite-Element-Methoden, Springer, Berlin, Heidelberg, NewYork 5. Beskos D.E., 1997, Boundary element methods in dynamic analysis: Part II (1986-1996),App. Mech. Rev., 50, 3, 149-198 6. Cho Y., Rose J.L., 1996, A boundary element solution for mode conversion study of the edge reflections of Lambwaves,J.Acoust. Soc. Am.,99, 2079-2109 7. Dauksher W., Emery A.F., 2000, The solution of elastostatic and elasto- dynamic problems with Chebyshev spectral finite elements, Comput. Methods Appl. Mech. Engrg., 188, 217-233 8. Delsanto P.P., Schechter R.S., Mignogna R.B., 1997, Connectionma- chine simulation of ultrasonic wave propagation in materials III: The three- dimensional case,Wave Motion, 26, 329-339 9. Doyle J.F., 1997,Wave Propagation in Structures, Springer, Berlin 10. Fornberg B., 1987, The pseudospectralmethod; Comparisonswith finite dif- ferences for the elastic wave equations,Geophys., 52, 483-501 11. Giurgiutiu V., 2007, Structural Health Monitoring with Piezoelectric Wafer Active Sensors, Academic Press, Elsevier, San Diego 12. Graves R.W., 1996, Simulating seismicwave propagation in 3D elasticmedia using staggered-gridfinite differences,Bull. Seismol. Soc.Am.,86, 4, 1091-1106 13. HughesT.J.R., 1987,TheFinite ElementMethod: Linear Static andDynamic Finite Element Analysis, Prentice-Hall EnglewoodCliffs, NJ 14. Ikeda T., 1990, Fundamentals of Piezoelectricity, Oxford University Press, NewYork 15. Komatitsch D., Martin R., Tromp J., Taylor M.A., Wingate B.A., 2001,Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles, J. Comp. Acoust., 9, 2, 703-718 16. Ku H., Hirsh R., Taylor T., 1987, A pseudospectral method for solution of the three-dimensional incompressible Navier-Stokes equations, J. Comput. Phys. 17. Kudela P., Ostachowicz W., 2009, A multilayer delaminated composite beamandplate elements: reflections ofLambwavesat delamination,Mechanics of Advanced Materials and Structures, 16, 3, 174-187 18. Kudela P., Zak A., Krawczuk M., Ostachowicz W., 2007,Modelling of wave propagation in composite plates using the time domain spectral element method, J. Sound and Vib., 302, 728-745 902 R.T. Schulte, C.-P. Fritzen 19. Lammering R., Mesecke-Rischmann S., 2003, Multi-field variational for- mulations and related finite elements for piezoelectric shells, Smart Mater. Struct., 12, 904-913 20. Mankins J.C., 1995, Technology readiness levels,White Paper, 6 21. Moser F., Jacobs L.J., Qu J., 1999, Modeling elastic wave propagation in waveguideswith the finite elementmethod,NDT&E International,32, 225-234 22. Ostachowicz W., 2008, Damage detection of structures using spectral finite element method,Computers and Structures, 86, 3/5, 454-462 23. Ostachowicz W., Kudela P., 2010, Wave propagation numerical models in damage detection based on the time domain spectral element method, IOP Conf. Series: Materials Science and Engineering, 10, 012068 24. Patera A.T., 1984, A spectral element method for fluid dynamics: Laminar flow in a channel expansion, J. Comput. Phys., 54, 468-488 25. PengH.,MengG.,Li F., 2009,Modelling ofwavepropagation inplate struc- tures using three-dimensional spectral element method for damage detection, Journal of Sound and Vibration, 320, 942-954 26. Pozrikidis C., 2005, Introduction to Finite and Spectral Element Methods using Matlab, Chapman&Hall/CRC, Boca Raton 27. Rucka M., 2010, Experimental and numerical study on damage detection in an L-joint using guided wave propagation, Journal of Sound and Vibration, 329, 1760-1779 28. Seriani G., 2004, Double-grid Chebyshev spectral elements for acoustic wave modelling,Wave Motion, 39, 351-360 29. Seriani G., Priolo E., 1994, Spectral element method for acoustic wave simulation in heterogeneousmedia,Finite Elem. Anal. Des., 16, 337-348 30. Wagner W., 1985, Eine geometrisch nichtlineare Theorie schubelastischer Schalen mit Anwendung auf Finite-Elemente-Berechnungen von Durchschlag- und Kontakt problemen, Dissertation, Universität Hannover, Hannover 31. Yang S., Ngoi B., 1999, General sensor equation and actuator equation for the theory of laminated piezoelectric plates, Smart Mater. Struct., 8, 411-415 32. YimH., SohnY., 2000,Numerical simulationandvisualizationof elasticwaves using mass-spring lattice model, IEEE Transactions on Ultrasonics, Ferroelec- trics and Frequency Control, 47, 549-558 33. Zienkiewicz O.W., Taylor R.L., Zhu J.Z., 2005,The Finite Element Me- thod, 6th ed., Elsevier Butterworth-Heinemann, Oxford Simulation of wave propagation... 903 Symulacja propagacji fali w kompozytach tłumionych efektem sprzężenia piezoelektrycznego Streszczenie Wpracy przedstawiono efektywną i dokładną procedurę symulacji układów z sys- temem monitorowania stanu (SHM) stosowaną na wstępnym etapie projektowania takich układów i pozwalającą na redukcję czasu i kosztów niezbędnych analiz wstęp- nych. Procedurę oparto na metodzie elementów spektralnych w dziedzinie czasu, co umożliwiło uzyskanie optymalnie określonychmacierzybezwładności i przez to zdecy- dowane skrócenie czasu całkowania. Przybliżono podstawy teoretycznemetody i opi- sano elementy spektralne dla płaskich powłok. Przedyskutowano nową metodologię badań symulacyjnychpozwalającychnawprowadzenie tłumienia anizotropowegooraz efektu sprzężenia piezoelektrycznegowujęciu elementów spektralnych.Na kilku przy- kładach zaprezentowano dokładność obliczeń numerycznych poprzez porównanie ich z wynikami eksperymentu oraz wykazano ich przydatność do analizy bardziej złożo- nych układów, jak np. zakrzywionych paneli z lokalnymi usztywnieniami. Manuscript received March 23, 2011; accepted for print April 19, 2011