Three-dimensional BEM for transient dynamic analysis of piezoelectric and anisotropic elastic solids

A boundary element approach based on the Green’s functions in integral representation and the convolution quadrature method is presented. Proposed approach is designed for analyzing 3D initial boundary-value problems of the dynamics of general anisotropic elastic and piezoelectric linear homogeneous solids with mixed boundary conditions. Numerical modelling of transient dynamics of elastic anisotropic and piezoelectric three-dimensional solids is carried out to demonstrate the potential of the developed boundary element software. Obtained solutions are compared with the corresponding FEM results and results of the dynamic experiment. A numerical technique based on the exact Laplace-domain boundary integral equations for the direct approach of 3D linear theories of anisotropic elasticity and piezoelectricity is employed. The BEM scheme is constructed using the collocation method and the convolution quadrature method in the form of a stepping method for numerical inversion of integral Laplace transform. Results of the stepped BE-modelling of the problems when a transient force is acting on 3D piezoelectric and anisotropic elastic homogeneous solids are presented.


Introduction
Three-dimensional transient problems of linear anisotropic elasticity and piezoelectricity can be successfully treated by the direct Boundary Element Method (BEM).Yet such BEM formulations hardly reported in the literature due to the lack of the exact closed-form expressions of dynamic fundamental solutions (often referred to as Green's functions) and high complexity of their integral representations [1,2].The Dual Reciprocity BEM (DRBEM) formulations which require only static fundamental solutions were proposed [3][4][5] to avoid this complication.In this study we present a Laplace domain direct BEM formulation based on exact singular boundary integral equation (BIE) and on integral expressions of corresponding dynamic fundamental solutions.

Basic equations 2.1. Linear anisotropic elasticity
First we consider a three-dimensional linear and homogeneous anisotropic elastic solid of volume and boundary .In the absence of body forces, equations of motion are given by where C E i jkl , u i and ρ represent the elasticity tensor, displacement vector and mass density, respectively.
The initial and boundary condition are given as follows a Corresponding author: igumnov@mech.unn.ru where = u ∪ t and t is the traction vector, defined as where n is the unit outward normal vector to the boundary and σ i j is the Cauchy stress tensor.

Linear piezoelectricity
Under quasi-electrostatic assumption and in the absence of body forces, the coupled equations of motion for a homogeneous linear piezoelectric solid take the following form [5,6]: where e i jk and κ ik are piezoelectric and dielectric permittivity tensors, respectively.In addition, φ represents the electric potential and it has been taken into account that piezoelectric materials possess no free electric charges.
In order to combine elastic and electric variables into corresponding generalized vectors and tensors we employ the following conventional contracted notations [5] EPJ Web of Conferences where q denotes the surface charge density, defined as q = D j n j (11) and D j is the electric displacements.
We now can rewrite the equations of motion ( 6)-( 7) in compacted form as follows Initial conditions for the elastic quantities have to be prescribed The following boundary conditions are assumed for the elastic field for the electrical field where φ and q are the parts of the boundary with the prescribed electric potential and surface charge density, respectively.

Laplace transform formulation
As is can be observed, Eq. ( 12) has essentially the same form as Eq. ( 1) and for the sake of simplicity we now formulate the generalized initial-boundary value problem for three-dimensional linear and homogeneous anisotropic elastic or piezoelectric solid of volume and boundary , assuming zero body forces and zero initial conditions: where N denotes the dimension of the problem; N = 3 for elasticity, N = 4 for piezoelectricity and we assume that the proper boundary conditions are prescribed.
Applying the Laplace transform in time variable to Eq. ( 19) we get where a bar over a variable hereinafter represents its Laplace transform and s is the transform variable.Equation (20) along with corresponding boundary conditions represents the generalized boundary value problem in Laplace domain.

Boundary integral equation and BEM formulation
Besides the case of elasticity, the generalized Betti-Rayleigh reciprocal theorem also holds for piezoelectricity.Thus, assuming zero body forces and zero initial conditions a boundary integral representation of the solution of Eq. ( 20) can be obtained in the form where c jk are the free term coefficients, p.v. indicates that the integral is treated in the sense of its Cauchy principal value, ḡij (x, y, s) and hij (x, y, s) represent the corresponding (elastic or piezoelectric) displacement and traction fundamental solutions in Laplace domain.
In this study we employ matched boundary element approach [7]: the boundary is approximated by eight-node quadratic quadrilateral elements, the boundary variables Ūk and Tk are approximated by linear and constant elements, respectively.
Once Eq. ( 21) is discretized and boundary conditions are taken into account, the resulting system of linear equations is solved for unknown boundary generalized displacements and tractions for the given value of Laplace transform variable s.In order to obtain solution in timedomain we use convolution quadrature method in the form of stepping method for numerical inversion of Laplace transforms [8][9][10].

Dynamic anisotropic elastic and piezoelectric fundamental solutions
Laplace domain displacement fundamental solutions for anisotropic elastic or piezoelectric solids are obtained by applying the Laplace transform to the corresponding timedomain expressions [1,2]: where ḡS jk and ḡR jk denote the static and dynamic parts of the fundamental solution, respectively.
Decomposition of dynamic fundamental solution into two parts Eq. ( 22) allows using the interpolation scheme for calculation of static part first proposed by Wilson and Cruse [11].In this study we employ linear Lagrange interpolation on the 800 × 800 grid.
For 3D anisotropic elastic problems, ḡS jk and ḡR jk are expressed as follows [12]: where λ m , E jm are eigenvalues and corresponding eigenvectors of Christoffel tensor jk (d); c m and k m are the phase velocities and wave numbers, which are functions of the propagation vector n.
And for 3D piezoelectric problems ḡS jk and ḡR j p are expressed as follows [13]: where A m j p is the adjugate of matrix L j p − ρc 2 m δ j p and it was taken into consideration that the eigenvectors need not to be evaluated explicitly, hence λ m are the eigenvalues of matrix L j p and Q is number of distinct eigenvalues.
In addition, the traction fundamental solution h jk is obtained using the following relation h jm (x, y, s) = C i jkl ḡkm,l n i . (36)

Example A
First we consider 3D anisotropic elastic rod which is clamped at one end, and subjected to Heavisidetype traction t 2 = t * 2 H (t), t * 2 = −1N/m 2 at the opposite free end, as illustrated in Fig. 1.The remaining surfaces are traction free.The material considered is an orthotropic 65% graphite-35% epoxy [14] with density ρ = 1600 kg/m 3 and whose material principal axes are rotated sequentially by 30 • around x 2 and then by the same angle around x 3 .This yields the following fully populated elasticity tensor, which is given in the Voigt notation as: GPa.
(37) The employed BEM mesh consists of 504 boundary elements and 506 nodes.The transient response of displacements u 2 (t) at boundary point (0, 3, 0) m is shown in Fig. 2 where BEM solution are given in comparison with the corresponding FEM results.It can be observed from   this figure that the agreement between the BEM and FEM results are very good.

Example B
As another example, consider the problem depicted in Fig. 3: an anisotropic elastic prismatic bar is simply supported at one side on two areas with width b = 1.5 mm each, that are located symmetric with respect to the x 3 axis at a distance c = 45 mm from each other.Bar is subjected to a uniformly distributed normal load t 3 = t * 3 f (t), t * 3 = −9.724• 10 7 /m 2 acting on the area with width a = 1.4 mm located exactly in the middle on the opposite side.Function f (t) is given in Fig. 4.
In the first case we consider an orthotropic material for which we also conducted dynamic experiment.The employed BEM mesh consists of 540 boundary elements and 542 nodes.The results from the BEM analysis and dynamic experiment for the vertical displacements −u 3 (t) at boundary point O (0, 0, 0) m are in the good agreement as it shown in Fig. 5.
In the second case we consider the same triclinic material as in Example A. Obtained BEM results for the vertical displacements −u 3 (t) at boundary point O (0, 0, 0) m are given in Fig. 6.

Example C
For the final example we consider a -shaped piezoelectric solid [5] under Heaviside-type mechanical load t 1 (t) = t * 1 H (t), t * 1 = −10 8 N m 2 at one face and applied voltage φ(t) = φ * = 0 V at another face, which is clamped    Mesh used in computations is regular and consists of 350 boundary elements and 352 nodes.Obtained BEM solutions for u 1 (t) and φ(t) at boundary point (0, 0, 0.1) m are given in Fig. 8 and Fig. 9, respectively.

Concluding remarks
The proposed Laplace domain direct BEM formulation for transient piezoelectric and anisotropic elastic analyses produce highly accurate results, which was demonstrated comparing the BEM solutions with the corresponding FEM results and with the results of dynamic experiment.

Figure 3 .
Figure 3. Example B: geometry and boundary conditions.