A boundary element formulation in time domain for viscoelastic solids

Viscoelastic solids may be effectively treated by the boundary element method (BEM) in the Laplace domain. However, calculation of transient response via the Laplace domain requires the inverse transform. Since all numerical inversion formulas depend heavily on a proper choice of their parameters, a direct evaluation in the time domain seems to be preferable. On the other hand, direct calculation of viscoelastic solids in the time domain requires the knowledge of viscoelastic fundamental solutions. 
 
 
 
Such solutions are simply obtained in the Laplace domain with the elastic–viscoelastic correspondence principle, but not in the time domain. Due to this, a quadrature rule for convolution integrals, the ‘convolution quadrature method' proposed by Lubich, is applied. This numerical quadrature formula determines their integration weights from the Laplace transformed fundamental solution and a linear multistep method. Finally, a boundary element formulation in the time domain using all the advantages of the Laplace domain formulation is obtained. Even materials with complex Poisson ratio, leading to time-dependent integral free terms in the boundary integral equation, can be treated by this formulation. 
 
 
 
Two numerical examples, a 3D rod and an elastic concrete slab resting on a viscoelastic half-space, are presented in order to assess the accuracy and the parameter choice of the proposed method. Copyright © 1999 John Wiley & Sons, Ltd.


INTRODUCTION
The boundary element method (BEM) has become a widely used numerical tool in statics and dynamics. A review about the eorts in dynamics is published by Beskos. 1,2 The main advantages of the method are reduction of the problem dimension by one and the implicit ful®lment of the radiation condition for unbounded domains. Due to the last advantage the BEM seems to be preferable for the calculation of in®nite or semi-in®nite domains, e.g. soil.
Modelling soil as an elastic solid can only be a crude approximation of the real material. A better description of the material properties can be achieved using viscoelastic constitutive equations. In these equations, normally, the stress strain relation is a dierential or integral equation of hereditary type. 3 Improved curve ®tting of measured material properties by constitutive equations with fewer parameters is achieved with the concept of fractional dierintegration. 4,5 A review about fractional calculus applied to dynamic problems has recently been given by Rossikhin and Shitikova. 6 Viscoelastic boundary element formulations are mostly published for the quasi-static case (see, for example, References 7 9), or in dynamics using a frequency or Laplace domain representation of the governing integral equation. These formulations are developed by applying the elastic viscoelastic correspondence principle to the elastodynamic boundary element formulation, e.g. Reference 10 for a frequency domain or Reference 11 for a Laplace domain formulation. In these formulations, the complex moduli concept is used which can be extended to fractional operator viscoelasticity allowing not only integer powers of the frequency, e.g. Reference 12. Calculation of transient response, however, requires the inverse transformation. Since all numerical inversion formulas depend on a proper choice of their parameters, 13 a direct evaluation in the time domain seems to be preferable. However, formulations directly in the time domain require the knowledge of viscoelastic fundamental solutions which are not yet known for the general visco elastodynamic case. Only for a simple Maxwell model, a solution has been obtained analytically by Gaul and Schanz 14 and has been implemented in a boundary element formulation. 15 Based on the frequency domain fundamental solution with subsequent inverse transformation, a 1D solution has been proposed by Wolf and Dabre. 16 In the 3D case, Gaul and Schanz 17 developed a formulation for a generalized (with fractional derivatives) 3-parameter model using the Laplace transformed fundamental solution which is inverted within each time step. Recently, Schanz and Antes 18 published a viscoelastic formulation based on the so-called Convolution Quadrature Method' proposed by Lubich. 19,20 In this formulation, the convolution integral is numerically approximated by a quadrature formula whose weights are determined by the Laplace transform of the fundamental solution and a linear multistep method. A comparison of the above-mentioned two time domain formulations and the Laplace transformed formulation with a subsequent inverse transformation was presented by Gaul and Schanz. 21 However, all published formulations assumed the same viscoelastic behaviour of the deviatoric and the hydrostatic part of the stress strain relations, leading to a constant, real numbered Poisson ratio which does not represent all viscoelastic materials. Now, here, the formulation based on the convolution quadrature will be extended to the case of a complex Poisson ratio. For this, ®rst, the constitutive equations are recalled and also, in the second step, the governing integral equation for linear visco elastodynamics. With a spatial discretization and applying the`Convolution Quadrature Method' for temporal discretization a time-stepping procedure is achieved. Some tests of the method with respect to the spatial and temporal discretization will show the eciency of the proposed method. Finally, as a more realistic example, an elastic foundation slab resting on a viscoelastic half space is studied.

CONSTITUTIVE EQUATION
Decomposing the stress and strain tensors s ij , e ij into the hydrostatic parts d ij s kk /3, d ij e kk /3 and the deviatoric parts s ij , e ij and assuming for both parts a linear isothermal viscoelastic material, e.g. a generalized three parameter model, two independent sets of constitutive equations exist in the Laplace domain: when ^ v{} denotes the Laplace transformed functions and s the Laplace variable. In equation (1) the elastic material constants Young's modulus E and Poisson's ratio are used. In this equation, and in the following, Latin indices receive the values 1, 2 or 1, 2, 3 in 2D or in 3D, respectively, where summation convention is implied over repeated indices. The index H or D shows that the constitutive parameters p, q and the fractional power a can be dierent for the hydrostatic part H and the deviatoric part D of the stress strain relation. A plausible model of the viscoelastic phenomenon should predict non-negative internal work, a non-negative rate of energy dissipation 4 and ®nite wave velocities. Therefore, the parameters are constrained by 22 Comparing the viscoelastic constitutive equations (1) with Hook's law, the elastic viscoelastic correspondence principle is obtained.
This means that every elastodynamic solution of a distinct problem can be converted to the solution of the related viscoelastic problem by inserting the correspondence (3). 3 In equation (3), the dependence of Poisson's ratio from the Laplace parameter s is clearly observed, which corresponds to a time dependent Poisson ratio. If, however, the parameters satisfy the relation which means that the deviatoric and the hydrostatic part of the stress strain relation behave similarly, then Poisson's ratio is also for viscoelasticity no function of s, i.e. time invariant, and equal to the elastic case. 3

VISCOELASTIC BOUNDARY ELEMENT FORMULATION
Assuming homogeneity and a linear strain displacement relation the dynamic behaviour of a viscoelastic structure as well as the wave propagation in a domain O with boundary G is governed by the boundary integral equation 23 (* denotes a convolution in time) where u j u j (x,t) and p j p j (x,t) is the displacement and the traction vector at point x and time t, respectively. The fundamental solutions of the displacements and tractions are denoted by U ij (x,y,t) and P ij (x,y,t), respectively. They are available in the Laplace domain from the elastodynamic fundamental solutions (see Reference 24) via the elastic viscoelastic correspondence principle. 18 As in the elastostatic case, the ®rst integral in equation (5) is weakly singular, and the second integral has to be de®ned in the sense of a Cauchy principal value, i.e. the integral is strongly singular. Further on, due to the time dependence of Poisson's ratio (t), now, in general, c ij (t) is also time dependent. However, the integral free terms c ij (t) are independent of time and equal to the elastostatic case, if Poisson's ratio is independent of time. This is true if equation (4) is valid, i.e. the hydrostatic and the deviatoric part of the stress strain relation behave similarly.
The numerical implementation of equation (5) requires approximations of boundary ®eld variables in both space and time. Therefore, the boundary surface G is discretized by E isoparametric elements G e , where F polynomial shape functions N f e (x) are de®ned. Hence, with the time dependent nodal values u ef j (t) and p ef j (t) the following representation is adapted: Inserting these`ansatz' functions in equation (5) results in When the time t is divided into N intervals of equal duration Dt so that t NDt, the convolution integral between the fundamental solutions U ij (x,y,t) or P ij (x,y,t) and the nodal values p ef i (t) or u ef i (t), respectively, may be performed by the convolution quadrature method proposed by Lubich. 10 This quadrature formula numerically approximates a convolution integral y(t) f(t)*g(t) by the ®nite sum ynDt X n k0 o nÀk DtgkDt; n 0; 1; . . . ; N 8 The integration weights o n (Dt) are the coecients of the power series for the function f Ã (g(z)/Dt) at the point g(z)/Dt. Herein, g(z) is the quotient of the characteristic polynomials of a linear A (a)stable multistep method, e.g. gz 3 2 À 2z 1 2 z 2 for the backward dierential formula of second order (BDF 2). The coecients of this power series are calculated by the integral with x l e il 2p L , being the radius of a circle in the domain of analyticity of f Ã (z). The integral in equation (9) is approximated by a trapezoidal rule with L equal steps 2p/L, after a transformation to polar co-ordinates. Further details are described in Reference 18 as well as in the original work of Lubich. 19 The characteristic advantage of the quadrature rule (8) is that only the Laplace transformed function f Ã is used.
Applying the quadrature formula (8) to the integral equation (7) results in the following boundary element time-stepping formulation for n 0, 1, . . . , N.  (12) and (13)), for the regular integrals a standard Gaussian quadrature rule and for the singular integrals the well known techniques from elastodynamics can be used. This means that the weakly singular integrals in (12) are regularized by co-ordinate transformation and the strongly singular integrals in (13) by the method suggested by Guiggiani and Gigante. 25 In order to arrive in equation (10) at systems of algebraic equations, collocation is used at every node of the shape functions N f e (x), and, ®nally, a direct equation solver is applied.
Before the method is applied to a numerical example, the time dependence of the integral free terms c ij (t) is shown in Figure 1. For example, the coecient c 13 of a rectangular corner c 13 s 1 8p1 À s ( (s) denotes the lengthy expression in equation (3)) is convoluted with a unit step function u 3 H(t) . 1m. The related elastic Poisson ratio in this test 0 . 3, leading to c 13 0 . 0568, and a BDF 2 as the underlying multistep method is used. Three cases are considered in Figure 1: (i) hydrostatic damping the deviatoric part of the stress strain relation has an elastic behaviour; (ii) equal damping mechanisms for the deviatoric and hydrostatic part of the stress strain relation; and (iii) deviatoric damping the hydrostatic part of the stress strain relation has elastic behaviour. The initial values of the three cases are quite dierent. With increasing time a relaxation or a creep behaviour of the results is observed. The case with the equal damping values gives, as expected, a constant value which is equal to the elastic case.

NUMERICAL EXAMPLES
The propagation of waves in a 3D viscoelastic continuum has been calculated by the presented viscoelastic boundary element formulation. The numerical tests are divided into two main parts. First, the displacements and tractions of a 3D rod are compared to a 1D solution with respect to the in¯uence of the time step size, the underlying multistep method and the complex Poisson ratio. In the second part, an elastic concrete foundation slab on a viscoelastic half-space is considered.
For all the tests, the parameters of the convolution quadrature are chosen as suggested by Lubich: 20 L N and N p e with e 10 À10 . For a better comparison with results of other example, the dimensionless parameter b c 1 Dt/r e is used instead of the time step size. This parameter is calculated using the viscoelastic compression wave velocity 21 c 1 and the mean value of a characteristic element length r e .

Three-dimensional rod
The problem geometry and material data of the 3D rod are shown in Figure 2. The rod is taken to be ®xed on one end, and is excited by a pressure jump according to a unit step function p y x; t 1 N m 2 Ht on the other free end. The remaining surfaces are traction free. Linear spatial shape functions on 56 triangles are used.
In most boundary element formulations in the time domain, the value b is restricted to a very small range where stable and satisfactory results are achieved. In Figure 3 the displacements for dierent b and BDF 2 as the underlying multistep method are depicted and compared to the 1D solution. The 1D solution is obtained with the elastic viscoelastic correspondence principle applied to the known solution of an elastic rod. 21 Mesh 1 is the mesh shown in Figure 2 and mesh 2 denotes a ®ner mesh (every triangle in Figure 2 is bisected). Both meshes show a critical time step size, corresponding to b 4 0 . 1, below which the results are unstable. For mesh 1 only the tendency, whereas for mesh 2 the instabilities, are clearly observed. Furthermore, the results of mesh 1 strongly depend on b, i.e. the solution deviates more from the 1D solution if b increases. Compared to mesh 2, there the results for all b, except b 0 . 1, are close to the 1D solution. The stronger dependency of the results from b, respectively from the time step size Dt, of mesh 1 has two reasons. First, the spatial integration on mesh 2 has twice the quality as on mesh 1 due to the smaller elements. Second, the same b for mesh 1 and mesh 2 is achieved by half the time step size for mesh 2 compared to mesh 1. Because a ®ner time discretization will always lead to a better approximation of the time history, the time dependency of the displacements must be better approximated by mesh 2. This gives reason to the conclusion, that for a ®ne enough mesh the results are not dependent on the time step size if 0 . 1 5 b 5 1 is regarded.
Not only the time step size but also the underlying multistep method in¯uences the results. Therefore, in Figure 4 the longitudinal displacements and tractions vs. time for dierent multistep methods using mesh 1 are presented. The linear multistep methods backward dierential Figure 2. Step function excitation of a free-®xed rod formulation of ®rst, second and third order (BDF 1, BDF 2 and BDF 3) are compared with the trapezoidal rule. In the test, the relevant optimal time step size for each of the dierent methods was chosen. The traction solution clearly indicates that the BDF 3 and the trapezoidal rule lead to not acceptable results. However, with BDF 1 and BDF 2 good solutions are achieved. The BDF 1 and 2 are A-stable methods, whilst the BDF 3 is an A (a)-stable method. The trapezoidal rule is certainly A-stable but not stable at in®nity. These results might indicate that the usable multistep methods for visco elastodynamics should be restricted to A-stable methods which are also stable at in®nity. However, note that for a proper solution with the BDF 1, a very small time step size has to be taken to avoid numerical damping. Finally, it should be mentioned that the BDF 2 seems to be the best choice. Therefore, in the following, all tests are performed with the BDF 2 as the underlying multistep method.
In the ®nal experiment, dierent viscoelastic behaviour of the stress strain relation is studied. First, the extreme, but not realistic, case of pure hydrostatic damping i.e. the deviatoric part of the stress strain relation is assumed to be elastic is compared to the more realistic cases of equal damping mechanisms and the pure deviatoric case, i.e. the hydrostatic part of the stress strain relation is assumed to be elastic. In Figure 5, the displacement at point P is plotted vs. time for these three cases. This experiment shows that the deviatoric part of the stress strain relation has more in¯uence on the damping behaviour than the hydrostatic part. This is clearly indicated by the nearly identical results for the deviatoric case as for the case where both parts are viscous. On the contrary, when only the hydrostatic part is damped, only a very small in¯uence is observed.

Elastic concrete slab foundation on a viscoelastic clay half-space
The propagation of waves in an elastic concrete foundation slab (1 m Â 1 m Â 0 . 5 m) bonded on a viscoelastic clay half-space will be analysed. Both domains are coupled by a substructure technique based on displacement-and traction-continuity at the interface. With this assumption, any uplifting or other non-linear contact eects are neglected. The problem geometry and the associated boundary discretization are shown in Figure 6. The viscoelastic half-space is modelled with the material data of clay, where only for the deviatoric part of the stress strain relation damping is assumed. The spatial discretization is done with linear shape functions on triangles, and a BDF 2 is used as the underlying multistep method. Due to the very dierent wave velocities of the two domains and varying element sizes, the value b  cannot be chosen optimal for both domains. The experience shows that the unbounded halfspace is not so sensitive to an optimal choice of b as the bounded domain of the slab. Therefore, the time step is oriented at the material and discretization data of the slab. It is chosen with Dt 0 . 0025, and kept constant during all tests. In Figure 7, the vertical displacements u z at point A vs. time are depicted for dierent viscous damping values p. The great in¯uence of the viscosity is observed. Lower values of p, i.e. increasing viscosity, leads to a stiening of the material, indicated by the reduction of the displacement amplitudes and higher wave velocities. This example shows the proposed method is not only applicable to academic problems as to the rod; it can also be applied to more realistic examples. The insensitivity to the optimal choice of the time step size makes this method feasible for realistic applications.

CONCLUSIONS
In the paper, a novel visco elastodynamic boundary element formulation in time domain is presented. The constitutive equations are formulated in the dierential operator form taking also fractional derivatives into account, and the viscoelastic boundary integral equation is recalled. Next, the spatial discretization with standard techniques is introduced. Due to the lack of time dependent fundamental solutions, the`convolution quadrature method', 19 a quadrature rule for convolution integrals, is used. This numerical quadrature formula determines their integration weights from the Laplace transformed fundamental solution and a linear multistep method. These fundamental solutions can easily be obtained from the elastodynamic ones by applying the elastic viscoelastic correspondence principle. This leads ®nally to a boundary element formulation in the time domain using all the advantages of the Laplace domain formulation. Even materials with complex Poisson ratio, leading to time-dependent integral free terms in the boundary integral equation, can be treated by this formulation.
In two examples, a 3D rod and an elastic concrete slab bonded to a viscoelastic half-space, numerical studies show the accuracy of the proposed formulation.