Simulation of wave propagation in plate structures by using new spectral element with piezoelectric coupling
Feng Wang^{1} , Xinwei Wang^{2} , Zhaohui Feng^{3}
^{1, 2, 3}State Key Laboratory of Mechanics and Control of Mechanical Structures Nanjing University of Aeronautics and Astronautics, No. 29 Yudao Street, Nanjing 210016, China
Journal of Vibroengineering, Vol. 15, Issue 1, 2013, p. 214222.
Received 23 December 2012; accepted 28 February 2013; published 31 March 2013
This paper presents an efficient approach to simulate Lamb wave propagations in thin plate structures by using new time domain spectral plate elements. A novel approach is proposed to incorporate the coupling of piezoelectric transducers within the twodimensional plate element. The diagonal mass matrix is obtained by using a simple method with less computational effort. Detailed formulations are given. The benchmark problem of a thin aluminum plate with two surfacemounted piezoelectric transducers was investigated in detail. Comparisons are made with results obtained by using ABAQUS to verify the developed spectral plate element. It is shown that the proposed element can efficiently simulate the propagation of Lamb wave generated by a piezoelectric actuator and picked up by a piezoelectric sensor.
Keywords: piezoelectric transducers, spectral plate element, wave propagation, Lamb wave.
1. Introduction
It is well known that detecting small damage in structures is an important and challenging task [12]. Recently, the use of Lamb wave based technology in structural health monitoring has made significant progress [35]. 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 [818]. 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 [818] 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 theredimensional (3D) 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 nondiagonal 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 3D behavior of the Lamb wave propagation in plates by using a 2D 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×6node line (or simply node) 2D 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 2D 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 3D displacement field of the plate can be written as:
$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{2{z}^{2}}{{h}^{2}}\frac{z}{h}\right){w}^{B}+\left(1\frac{4{z}^{2}}{{h}^{2}}\right){w}^{M}+\left(\frac{2{z}^{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. 1. Sketches of spectral plate element: a) with a PZT layer, b) without a PZT layer
a)
b)
For small strains, the 3D straindisplacement relations can be written by:
${\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},\text{\hspace{0.17em}\hspace{0.17em}}$
${\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{2{z}^{2}}{{h}^{2}}\frac{z}{h}\right)\frac{\partial {w}^{B}}{\partial y}+\left(1\frac{4{z}^{2}}{{h}^{2}}\right)\frac{\partial {w}^{M}}{\partial y}+\left(\frac{2{z}^{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{2{z}^{2}}{{h}^{2}}\frac{z}{h}\right)\frac{\partial {w}^{B}}{\partial x}+\left(1\frac{4{z}^{2}}{{h}^{2}}\right)\frac{\partial {w}^{M}}{\partial x}+\left(\frac{2{z}^{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}\text{\hspace{0.17em}}+\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:
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 2D spectral element, the displacement field in the element can be assumed as:
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 onedimensional Lagrange interpolation functions defined by:
where the distribution of nodes, ${\xi}_{i}$, ${\eta}_{i}$, can be the Gauss–Lobatto–Legendre (GLL) points, or the Chebyshev points, or the approximate Lebesgueoptimal grid points, i.e., the expandedChebyshev grid points [2122]. 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 expandedChebyshev 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:
where ${H}_{Vi}$ is the shape functions and defined by:
2.2. Constitutive relations
For the host plate, the 3D stressstrain relation is given by:
where $\left[C\right]$ is the elasticity matrix of the material.
For the isotropic piezoelectric element, the adopted constitutive relations are given by:
where $\left[e\right]$ and $\left[\epsilon \right]$ are matrices of the piezoelectric constant and the dielectric constant measured at zero strain, $\left\{D\right\}$ and $\left\{E\right\}$ are vectors of the electric displacement and electric field, respectively.
2.3. Mass and stiffness matrices
Since the general equations can be found in many reference papers [1520] and textbooks [23], therefore, the detailed formulations of the mass matrix and stiffness matrix for element with and without a PZT layer are omitted.
Due to $\sum _{i=1}^{N}\sum _{j=1}^{N}{N}_{ij}\left(\xi \text{,}\eta \right)=1$, the diagonal mass matrix of the element can be obtained by row summation and integrated by using Gaussian quadrature as:
$k,p=\mathrm{1,2},...,N;L=1,...,7;I=(k1)N+7(p1)+L,$
where $NN=(N+1)/2$, ${H}_{i}$, ${H}_{j}$ and ${\xi}_{i}$, ${\eta}_{j}$ are weights and abscissas of Gaussian quadrature, $\leftJ\right(\xi ,\eta \left)\right$ is the determinant of the Jacobian matrix, and ${\mu}_{L}$ for plates with and without a PZT layer are respectively given by:
and:
where $\rho $ and ${\rho}_{1}$ are the mass density of the host plate and PZT layer. It was found that the computational effort in obtaining the mass matrix is only about 25 % of that by the existing GLL quadrature rule [17].
The stiffness matrix of the element with a PZT layer can be partitioned by:
$\left[{K}_{UV}\right]={\left[{K}_{VU}\right]}^{T}={\int}_{{V}_{P}}{\left[B\right]}^{T}{\left[e\right]}^{T}\left[Bv\right]dV,$
$\left[{K}_{VV}\right]={\int}_{{V}_{P}}{\left[Bv\right]}^{T}\left[{\epsilon}^{s}\right]\left[Bv\right]dV,$
in which $V$ and ${V}_{p}$ are the volume of the host plate and PZT layer, $\left[B\right]$ and $\left[Bv\right]$ are strain matrix and electric field matrix, $\left[{K}_{uu}\right]$ and $\left[{K}_{vv}\right]$ are the stiffness matrix related to displacement and electric potential, and $\left[{K}_{uv}\right]$ is the coupled stiffness matrix, respectively. The stiffness matrix can be obtained numerically by using Gaussian quadrature or GLL rule if GLL points are used as the nodes. Detailed derivations are ommitted and may be found in reference [24].
2.4. The time integration scheme
In terms of the spectral plate element with piezoelectric coupling, the wellknown governing differential equations for the wave propagation in plate structures can be written by:
$\sum {K}_{VU}U\left(t\right)+\sum {K}_{VV}\Phi \left(t\right)=Q\left(t\right),$
where $M$ is the structural mass matrix, $U\left(t\right)$ is the generalized displacement vector and the over double dots denote the secondorder derivative with respect to time $t$, $\Phi \left(t\right)$, $F\left(t\right)$ and $Q\left(t\right)$ are vectors of electric potential, external applied force, and electric charge, respectively.
To simplify the solution procedures, the dynamic interface interaction between the PZT layer and the host structure is neglected. Thus, Eq. (14) is uncoupled. Due to the fact of that $M$ is a diagonal matrix, a crucial reduction of the complexity and the cost of the numerical time integration can be achieved by using the central finite difference method [25] to solve the first part of Eq. (14). The second part of Eq. (14) is merely used to compute the output voltage.
3. Results and discussion
To verify the proposed formulations, a benchmark problem presented in [19] is analyzed. A thin aluminum plate with two surfacemounted 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 fivepeak toneburst 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^{3}, 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. 2. Geometry of an aluminum plate
Fig. 3. Meshes by (a) proposed spectral plate element and (b) finite plain strain element
a)
b)
Table 1. Material properties of the PZT element
${C}_{11}^{E}$

${C}_{12}^{E}$

${C}_{13}^{E}$

${C}_{33}^{E}$

${C}_{44}^{E}={C}_{55}^{E}$

${C}_{66}^{E}$

147 GPa

105 GPa

93.7 GPa

113 GPa

23 GPa

21.2 GPa

${e}_{31}$

${e}_{33}$

${e}_{15}$

${\epsilon}_{11,r}^{}$

${\epsilon}_{33,r}^{}$

$\rho $

–3.09 C/m^{2}

16.0 C/m^{2}

11.64 C/m^{2}

1.0^{8}^{}F/m

8.04^{9}^{}F/m

7700 kg/m^{3}

Fig. 4. Excitation 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 3D 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. 5. Comparisons of the output voltage with data by ABAQUS
Fig. 6. Comparisons of the output voltage with data by ABAQUS
Fig. 7. Comparisons 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 sixvariable or sixdegree 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×6node with six DOFs per node spectral plate element considering the piezoelectric effect can be derived.
The same problem is then reanalyzed by using the 6×6node 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×6node 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 threedimensinal 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 2D 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.
Acknowledgements
The work is partially supported by the National Natural Science Foundation of China (50830201) and by the Priority Academic Program Development of Jiangsu Higher Education Institutions.
References
 Boller C., Staszewski W., Tomlinson G. Health Monitoring of Aerospace Structures. John Wiley & Sons, 2004. [Search CrossRef]
 Giurgiutiu V. Structural Health Monitoring with Piezoelectric Wafer Active Sensors. Academic Press, Elsevier, 2008. [Search CrossRef]
 Kessler S. S., Spearing S. M., Soutis C. Damage detection in composite materials using Lamb wave methods. Smart Materials and Structures, Vol. 11, Issue 2, 2002, p. 269278. [Search CrossRef]
 Ihn J. B., Chang F. K. Detection and monitoring of hidden fatigue crack growth using a builtin piezoelectric sensor/actuator network, 1: Diagnostics. Smart Materials and Structures, Vol. 13, Issue 3, 2004, p. 609620. [Search CrossRef]
 Ihn J. B., Chang F. K. Detection and monitoring of hidden fatigue crack growth using a builtin piezoelectric sensor/actuator network, 2: Validation using riveted joints and repair patches. Smart Materials and Structures, Vol. 13, Issue 3, 2004, p. 621630. [Search CrossRef]
 Willberg C., Koch S., Mook G., et al. Continuous mode conversion of Lamb waves in CFRP plates. Smart Materials and Structures, Vol. 21, Issue 7, 2012, p. 19. [Search CrossRef]
 Willberg C., Duczek S., Vivar Perez J. M., et al. Comparison of different higher order finite element schemes for the simulation of Lamb waves. Computer Methods in Applied Mechanics and Engineering, Vol. 241244, 2012, p. 246261. [Search CrossRef]
 Komatitsch D., Barnes C. H., Tromp J. Simulation of anisotropic wave propagation based upon a spectral element method. Geophysics, Vol. 4, 2000, p. 12511260. [Search CrossRef]
 Igawa H., Komatsu K., Yamaguchi I., et al. Wave propagation analysis of frame structures using the spectral element method. Journal of Sound and Vibration, Vol. 277, 2004, p. 10711081. [Search CrossRef]
 Palacz M., Krawczuk M., Ostachowicz W. The spectral finite element model for analysis of flexuralshear coupled wave propagation, Part 1: Laminated multilayer composite. Composite Structures, Vol. 68, 2005, p. 3744. [Search CrossRef]
 Sun H., Zhou L. Analysis of damage characteristics for cracked composite structures using spectral element method. Journal of Vibroengineering, Vol. 14, Issue 1, 2012, p. 430439. [Search CrossRef]
 Seriani G. 3D large scale wave propagation modeling by spectral element method on Cray T3E multiprocessor. Computer Methods in Applied Mechanics and Engineering, Vol. 164, Issue 12, 1998, p. 235247. [Search CrossRef]
 Kudela P., Krawczuk M., Ostachowicz W. Wave propagation modeling in 1D structures using spectral finite elements. Journal of Sound and Vibration, Vol. 300, 2007, p. 88100. [Search CrossRef]
 Zak A., Krawczuk M., Ostachowicz W. Propagation of inplane waves in an isotropic panel with a crack. Finite Elements in Analysis and Design, Vol. 42, 2006, p. 929941. [Search CrossRef]
 Peng H. K., Meng G., Li F. C. Modeling of wave propagation in plate structures using threedimensional spectral element method for damage detection. Journal of Sound and Vibration, Vol. 320, 2009, p. 942954. [Search CrossRef]
 Kudela P., Zak A., Krawczuk M., et al. Modeling of wave propagation in composite plates using the time domain spectral element method. Journal of Sound and Vibration, Vol. 302, 2007, p. 728745. [Search CrossRef]
 Xu C., Wang X. Efficient modeling and simulations of Lamb wave propagation in thin plates by using a new spectral plate element. Journal of Vibroengineering, Vol. 14, Issue 3, 2012, p. 11871199. [Search CrossRef]
 Wang X., Wang F., Xu C., et al. New spectral plate element for simulating Lamb wave propagations in plate structures. Journal of Nanjing University of Aeronautics & Astronautics, Vol. 44, Issue 5, 2012, p. 645651. [Search CrossRef]
 Kim Y., Ha S., Chang F. K. Timedomain spectral element method for builtin piezoelectricactuatorinduced Lamb wave propagation analysis. AIAA Journal, Vol. 46, Issue 3, 2008, p. 591600. [Search CrossRef]
 Schulte R. T., Fritzen C. P. Simulation of wave propagation in damped composite structures with piezoelectric coupling. Journal of Theoretical and Applied Mechanics, Vol. 49, Issue 3, 2011, p. 879903. [Search CrossRef]
 Boyd J. P. A numerical comparison of seven grids for polynomial interpolation on the interval. Computers and Mathematics with Applications, Vol. 38, 1999, p. 3550. [Search CrossRef]
 Brutman L. On the Lebesgue function for polynomial interpolation. SIAM Journal of Numerical Analysis, Vol. 15, 1978, p. 694704. [Search CrossRef]
 Ikeda T. Fundamentals of Piezoelectricity, Oxford University Press, New York, 1990. [Search CrossRef]
 Feng Z. Numerical Simulation of Piezoelectric Actuator Generated Lamb Wave Propagation in Plates. Master Thesis, Nanjing University of Aeronautics & Astronautics, Nanjing, China, March, 2013. [Search CrossRef]
 Subbaraj K., Dokainish M. A. A survey of direct timeintegration methods in computational structural dynamics – I. Explicit methods. Computers & Structures, Vol. 32, Issue 6, 1989, p. 13711386. [Search CrossRef]