# A boundary element formulation in time domain for viscoelastic solids

Martin Schanz

<figure><img src="image_1.png" /></figure>

## To cite this version:

Martin Schanz. A boundary element formulation in time domain for viscoelastic solids. International Journal for Numerical Methods in Biomedical Engineering, 1999, 15 (11), pp.799-809. 10.1002/(SICI)1099-0887(199911)15:113.0.CO;2-F. hal-01401766

HAL Id: hal-01401766

https://hal.science/hal-01401766

Submitted on 25 Nov 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A BOUNDARY ELEMENT FORMULATION IN TIME DOMAIN FOR VISCOELASTIC SOLIDS

MARTIN SCHANZ

Technical University Braunschweig, Institute of Applied Mechanics, P.O. Box 3329, D-38023 Braunschweig, Germany

## SUMMARY

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.

**KEY WORDS** viscoelasticity; boundary integral equations; time domain; transform methods

## 1. INTRODUCTION

The boundary element method (BEM) has become a widely used numerical tool in statics and dynamics. A review about the efforts in dynamics is published by Beskos.¹,² The main advantages of the method are reduction of the problem dimension by one and the implicit fulfilment of the radiation condition for unbounded domains. Due to the last advantage the BEM seems to be preferable for the calculation of infinite or semi-infinite 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 differential or integral equation of hereditary type.³ Improved curve fitting of measured material properties by constitutive equations with fewer parameters is achieved with the concept of fractional differ-integration.⁴,⁵ A review about fractional calculus applied to dynamic problems has recently been given by Rossikhin and Shitikova.⁶
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,¹³ 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¹⁴ and has been implemented in a boundary element formulation.¹⁵ Based on the frequency domain fundamental solution with subsequent inverse transformation, a 1D solution has been proposed by Wolf and Dabre.¹⁶ In the 3D case, Gaul and Schanz¹⁷ 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¹⁸ published a viscoelastic formulation based on the so-called ‘Convolution Quadrature Method’ proposed by Lubich.¹⁹,²⁰ 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.²¹ 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, first, 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 efficiency of the proposed method. Finally, as a more realistic example, an elastic foundation slab resting on a viscoelastic half space is studied.

## 2. CONSTITUTIVE EQUATION

Decomposing the stress and strain tensors $\sigma_{ij}$, $\varepsilon_{ij}$ into the hydrostatic parts $\delta_{ij}\sigma_{kk}/3$, $\delta_{ij}\varepsilon_{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:

$$ (1 + p^D s^{\alpha^D}) \hat{s}_{ij} = \frac{E}{1 + \nu} \hat{e}_{ij} (1 + q^D s^{\alpha^D}) (1 + p^H s^{\alpha^H}) \hat{\sigma}_{kk} = \frac{E}{1 - 2\nu} \hat{e}_{kk} (1 + q^H s^{\alpha^H}) \quad (1) $$

when $(^D) = \mathcal{L}\{\}$ 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 $\nu$ 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 *α* can be different 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⁴ and finite wave velocities. Therefore, the parameters are constrained by²²

$$0 \le p^D < q^D \quad 0 < \alpha^D < 2 \quad 0 \le p^H < q^H \quad 0 < \alpha^H < 2 \qquad (2)$$

Comparing the viscoelastic constitutive equations (1) with Hook's law, the elastic viscoelastic correspondence principle is obtained.

$$\begin{align}
E &\Leftrightarrow \frac{3E(1 + q^H s^{\alpha^H})(1 + q^D s^{\alpha^D})}{2(1 + \nu)(1 + q^H s^{\alpha^H})(1 + p^D s^{\alpha^D}) + (1 - 2\nu)(1 + q^D s^{\alpha^D})(1 + p^H s^{\alpha^H})} \nonumber \\
\nu &\Leftrightarrow \frac{(1 + \nu)(1 + q^H s^{\alpha^H})(1 + p^D s^{\alpha^D}) - (1 - 2\nu)(1 + q^D s^{\alpha^D})(1 + p^H s^{\alpha^H})}{2(1 + \nu)(1 + q^H s^{\alpha^H})(1 + p^D s^{\alpha^D}) + (1 - 2\nu)(1 + q^D s^{\alpha^D})(1 + p^H s^{\alpha^H})} \tag{3}
\end{align}$$

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).³

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

$$(1 + q^H s^{\alpha^H})(1 + p^D s^{\alpha^D}) = (1 + q^D s^{\alpha^D})(1 + p^H s^{\alpha^H}) \qquad (4)$$

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 $\Omega$ with boundary $\Gamma$ is governed by the boundary integral equation²³ (* denotes a convolution in time)

$$c_{ij}(\mathbf{y}, t) * u_j(\mathbf{y}, t) = \int_{\Gamma} U_{ij}(\mathbf{x}, \mathbf{y}, t) * p_j(\mathbf{x}, t) \, d\Gamma - \oint_{\Gamma} P_{ij}(\mathbf{x}, \mathbf{y}, t) * u_j(\mathbf{x}, t) \, d\Gamma \quad \text{for } \mathbf{y} \in \Omega \cup \Gamma \quad (5)$$

where $u_j = u_j(\mathbf{x},t)$ and $p_j = p_j(\mathbf{x},t)$ is the displacement and the traction vector at point $\mathbf{x}$ and time $t$, respectively. The fundamental solutions of the displacements and tractions are denoted by $U_{ij}(\mathbf{x},\mathbf{y},t)$ and $P_{ij}(\mathbf{x},\mathbf{y},t)$, respectively. They are available in the Laplace domain from the elastodynamic fundamental solutions (see Reference 24) via the elastic viscoelastic correspondence principle.¹⁸

As in the elastostatic case, the first integral in equation (5) is weakly singular, and the second integral has to be defined 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 $\nu(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 field variables in both space and time. Therefore, the boundary surface $\Gamma$ is discretized by $E$ isoparametric elements $\Gamma_e$, where $F$ polynomial shape functions $N_e^f(\mathbf{x})$ are defined. Hence, with the time dependent nodal values $u_j^e(t)$ and $p_j^e(t)$ the following representation is adapted:

$$u_j(\mathbf{x}, t) = \sum_{e=1}^{E} \sum_{f=1}^{F} N_e^f(\mathbf{x}) u_j^e(t), \quad p_j(\mathbf{x}, t) = \sum_{e=1}^{E} \sum_{f=1}^{F} N_e^f(\mathbf{x}) p_j^e(t) \quad (6)$$

Inserting these 'ansatz' functions in equation (5) results in

$$c_{ij}(\mathbf{y}, t) * u_j(\mathbf{y}, t) = \sum_{e=1}^{E} \sum_{f=1}^{F} \left\{ \int_{\Gamma} U_{ij}(\mathbf{x}, \mathbf{y}, t) N_e^f(\mathbf{x}) \, d\Gamma_e * p_j^e(t) - \oint_{\Gamma} P_{ij}(\mathbf{x}, \mathbf{y}, t) N_e^f(\mathbf{x}) \, d\Gamma_e * u_j^e(t) \right\} \quad (7)$$

When the time $t$ is divided into $N$ intervals of equal duration $\Delta t$ so that $t = N\Delta t$, the convolution integral between the fundamental solutions $U_{ij}(\mathbf{x},\mathbf{y},t)$ or $P_{ij}(\mathbf{x},\mathbf{y},t)$ and the nodal values $p_i^e(t)$ or $u_i^e(t)$, respectively, may be performed by the convolution quadrature method proposed by Lubich.¹⁰ This quadrature formula numerically approximates a convolution integral $y(t) = f(t)*g(t)$ by the finite sum

$$y(n\Delta t) = \sum_{k=0}^{n} \omega_{n-k}(\Delta t)g(k\Delta t), \quad n = 0, 1, \dots, N \quad (8)$$

The integration weights $\omega_n(\Delta t)$ are the coefficients of the power series for the function $\hat{f}(\gamma(z)/\Delta t)$ at the point $\gamma(z)/\Delta t$. Herein, $\gamma(z)$ is the quotient of the characteristic polynomials of a linear $A(\alpha)$-stable multistep method, e.g. $\gamma(z) = \frac{3}{2} - 2z + \frac{1}{2}z^2$ for the backward differential formula of second order (BDF 2). The coefficients of this power series are calculated by the integral

$$\omega_n(\Delta t) = \frac{1}{2\pi i} \int_{|z|=R} \hat{f}\left(\frac{\gamma(z)}{\Delta t}\right) z^{-n-1} dz \approx \frac{1}{L} \sum_{l=0}^{L-1} \hat{f}\left(\frac{\gamma(\xi_l)}{\Delta t}\right) \xi_l^{-n} \quad (9)$$

with $\xi_l = e^{il2\pi/L}$, $R$ being the radius of a circle in the domain of analyticity of $\hat{f}(z)$. The integral in equation (9) is approximated by a trapezoidal rule with $L$ equal steps $2\pi/L$, after a transformation to polar co-ordinates. Further details are described in Reference 18 as well as in the original work of Lubich.¹⁹ The characteristic advantage of the quadrature rule (8) is that only the Laplace transformed function $\hat{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, \dots, N$.

$$\sum_{k=0}^{n} \omega_{n-k}(\hat{c}_{ij}, \mathbf{y}, \Delta t) u_j(\mathbf{y}, k\Delta t) = \sum_{e=1}^{E} \sum_{f=1}^{F} \sum_{k=0}^{n} \left\{ \omega_{n-k}(\hat{U}_{ij}, \mathbf{y}, \Delta t) p_j^e(k\Delta t) - \omega_{n-k}(\hat{P}_{ij}, \mathbf{y}, \Delta t) u_j^e(k\Delta t) \right\} \quad (10)$$
with the weights corresponding to equation (9),

$$ \omega_{n-k}(\hat{c}_{ij}, \mathbf{y}, \Delta t) = \frac{1}{L} \sum_{l=0}^{L-1} \hat{c}_{ij} \left( \mathbf{y}, \frac{\gamma(\xi_l)}{\Delta t} \right) \xi_l^{-n} \quad (11) $$

$$ \omega_{n-k}(\hat{U}_{ij}, \mathbf{y}, \Delta t) = \frac{1}{L} \sum_{l=0}^{L-1} \int_{\Gamma} \hat{U}_{ij} \left( \mathbf{x}, \mathbf{y}, \frac{\gamma(\xi_l)}{\Delta t} \right) N_e^f(\mathbf{x}) \, d\Gamma_e \xi_l^{-n} \quad (12) $$

$$ \omega_{n-k}(\hat{P}_{ij}, \mathbf{y}, \Delta t) = \frac{1}{L} \sum_{l=0}^{L-1} \oint_{\Gamma} \hat{P}_{ij} \left( \mathbf{x}, \mathbf{y}, \frac{\gamma(\xi_l)}{\Delta t} \right) N_e^f(\mathbf{x}) \, d\Gamma_e \xi_l^{-n} \quad (13) $$

Obviously, although applying only the Laplace transformed fundamental solutions $\hat{U}_{ij}$ and $\hat{P}_{ij}$, a time stepping procedure directly in the time domain is achieved, i.e. a viscoelastic boundary element formulation in the time domain without the knowledge of the time dependent fundamental solutions is obtained. Further on, because the integral over each boundary element involves only the Laplace transformed fundamental solution (see equations (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.²⁵ In order to arrive in equation (10) at systems of algebraic equations, collocation is used at every node of the shape functions $N_e^f(\mathbf{x})$, and, finally, 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 coefficient $c_{13}$ of a rectangular corner

$$ \hat{c}_{13}(s) = \frac{1}{8\pi(1 - \nu(s))} $$

($\nu(s)$ denotes the lengthy expression in equation (3)) is convoluted with a unit step function $u_3 = H(t) \cdot 1m$. The related elastic Poisson ratio in this test $\nu = 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 different. 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.

## 4. 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 influence of the time step size, the underlying multistep method and the complex
<figure><img src="image_3.png" /><figcaption>Figure 1. Convolution $c_{13} * u_3$ vs. time for a viscoelastic material</figcaption></figure>

test data:

$$v = 0.3$$

hydrostatic damping:

$$p^D = 0 \text{ s}^{-1} \qquad p^H = 0.8 \text{ s}^{-1}$$

$$q^D = 0 \text{ s}^{-1} \qquad q^H = 1 \text{ s}^{-1}$$

equal damping:

$$p^D = 0.8 \text{ s}^{-1} \qquad p^H = 0.8 \text{ s}^{-1}$$

$$q^D = 1 \text{ s}^{-1} \qquad q^H = 1 \text{ s}^{-1}$$

deviatoric damping:

$$p^D = 0.8 \text{ s}^{-1} \qquad p^H = 0 \text{ s}^{-1}$$

$$q^D = 1 \text{ s}^{-1} \qquad q^H = 0 \text{ s}^{-1}$$

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:²⁰ $L = N$ and $\mathcal{R}^N = \sqrt{\epsilon}$ with $\epsilon = 10^{-10}$. For a better comparison with results of other example, the dimensionless parameter $\beta = c_1 \Delta t / r_e$ is used instead of the time step size. This parameter is calculated using the viscoelastic compression wave velocity²¹ $c_1$ and the mean value of a characteristic element length $r_e$.

### 4.1. Three-dimensional rod

The problem geometry and material data of the 3D rod are shown in Figure 2. The rod is taken to be fixed on one end, and is excited by a pressure jump according to a unit step function

$$p_y(\mathbf{x}, t) = 1 \frac{N}{m^2 H(t)}$$

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 $\beta$ is restricted to a very small range where stable and satisfactory results are achieved. In Figure 3 the displacements for different $\beta$ 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.²¹ Mesh 1 is the mesh shown in Figure 2 and mesh 2 denotes a finer mesh (every triangle in Figure 2 is bisected). Both meshes show a critical time step size, corresponding to $\beta \le 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 $\beta$, i.e. the solution deviates more from the 1D solution if $\beta$ increases.
<figure><img src="image_4.png" /><figcaption>Figure 2. Step function excitation of a free-fixed rod</figcaption></figure>

<table><thead><tr><th>material data:</th><th>geometry data:</th></tr></thead><tbody><tr><td><i>E</i> = 1 <sup><i>N</i></sup>/<sub><i>m</i><sup>2</sup></sub></td><td>length = 3 m</td></tr><tr><td><i>v</i> = 0</td><td>height = 1 m</td></tr><tr><td><i>ρ</i> = 1 <sup><i>kg</i></sup>/<sub><i>m</i><sup>3</sup></sub></td><td>width = 1 m</td></tr><tr><td><i>α</i><sup>D</sup> = <i>α</i><sup>H</sup> = 1</td><td></td></tr><tr><td><i>q</i><sup>D</sup> = <i>q</i><sup>H</sup> = 1 <i>s</i><sup>-1</sup></td><td></td></tr><tr><td><i>p</i><sup>D</sup> = <i>p</i><sup>H</sup> = 0.8 <i>s</i><sup>-1</sup></td><td></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>(a) mesh 1</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) mesh 2</figcaption></figure>

<caption>Figure 3. Longitudinal displacement of the free end vs. time: different time step sizes, i.e. different $\beta$</caption>

Compared to mesh 2, there the results for all $\beta$, except $\beta = 0.1$, are close to the 1D solution. The stronger dependency of the results from $\beta$, respectively from the time step size $\Delta t$, 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 $\beta$ for mesh 1 and mesh 2 is achieved by half the time step size for mesh 2 compared to mesh 1. Because a finer 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 fine enough mesh the results are not dependent on the time step size if $0.1 < \beta < 1$ is regarded.

Not only the time step size but also the underlying multistep method influences the results. Therefore, in Figure 4 the longitudinal displacements and tractions vs. time for different multistep methods using mesh 1 are presented. The linear multistep methods backward differential
<figure><img src="image_7.png" /><figcaption>Figure 4. Longitudinal displacements and tractions of the free and fixed end of the rod, respectively, vs. time: different multistep methods</figcaption></figure>

formulation of first, 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 different 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* (α)-stable method. The trapezoidal rule is certainly *A*-stable but not stable at infinity. These results might indicate that the usable multistep methods for visco elastodynamics should be restricted to *A*-stable methods which are also stable at infinity. 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 final experiment, different 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 influence 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 influence is observed.

## 4.2. 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 effects are neglected.
<figure><img src="image_8.png" /><figcaption>Figure 5. Longitudinal displacement of the free end vs. time: different behaviour of the hydrostatic and the deviatoric part of the stress-strain relation</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6. Elastic concrete slab on viscoelastic clay half-space: boundary element discretization, material data and loading function</figcaption></figure>

test data:

$$E = 1 \frac{N}{m^2} \qquad \rho = 1 \frac{kg}{m^3}$$

$$v = 0.3 \qquad \beta = 0.28$$

hydrostatic damping:

$$p^D = 0 s^{-1} \qquad p^H = 0.8 s^{-1}$$

$$q^D = 0 s^{-1} \qquad q^H = 1 s^{-1}$$

equal damping:

$$p^D = 0.8 s^{-1} \qquad p^H = 0.8 s^{-1}$$

$$q^D = 1 s^{-1} \qquad q^H = 1 s^{-1}$$

deviatoric damping:

$$p^D = 0.8 s^{-1} \qquad p^H = 0 s^{-1}$$

$$q^D = 1 s^{-1} \qquad q^H = 0 s^{-1}$$

clay:

$$E = 2 \cdot 10^6 \frac{N}{m^2}$$

$$v = 0.35$$

$$\rho = 1800 \frac{kg}{m^3}$$

$$\alpha^D = \alpha^H = 1$$

$$p^D = p s^{-1} \qquad p^H = 1 s^{-1}$$

$$q^D = 1 s^{-1} \qquad q^H = 1 s^{-1}$$

concrete:

$$E = 3 \cdot 10^8 \frac{N}{m^2}$$

$$v = 0.20$$

$$\rho = 2000 \frac{kg}{m^3}$$

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 different wave velocities of the two domains and varying element sizes, the value $\beta$
<figure><img src="image_10.png" /><figcaption>Figure 7. Displacement $u_z$ at point A vs. time for different values of $p$</figcaption></figure>

cannot be chosen optimal for both domains. The experience shows that the unbounded half-space is not so sensitive to an optimal choice of $\beta$ 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 $\Delta t = 0.0025$, and kept constant during all tests.

In Figure 7, the vertical displacements $u_z$ at point A vs. time are depicted for different viscous damping values $p$. The great influence of the viscosity is observed. Lower values of $p$, i.e. increasing viscosity, leads to a stiffening 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.

## 5. CONCLUSIONS

In the paper, a novel visco elastodynamic boundary element formulation in time domain is presented. The constitutive equations are formulated in the differential 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’,¹⁹ 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 finally 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.
REFERENCES

1. D. E. Beskos, 'Boundary element methods in dynamic analysis', *Appl. Mech. Rev.*, **40**(1), 1 23 (1987).

2. D. E. Beskos, 'Boundary element methods in dynamic analysis: Part II (1986–1996)', *Appl. Mech. Rev.*, **50**(3), 149 197 (1997).

3. R. M. Christensen, *Theory of Viscoelasticity*, Academic Press, New York, 1971.

4. R. L. Bagley and P. J. Torvik, 'On the fractional calculus model of viscoelastic behaviour', *J. Rheol.*, **30**(1), 133 155 (1986).

5. A. D. Drozdov, 'Fractional differential models in finite viscoelasticity', *Acta Mech.*, **124**, 155 180 (1997).

6. Y. A. Rossikhin and M. V. Shitikova, 'Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of solids', *Appl. Mech. Rev.*, **50**(1), 15 67 (1997).

7. M. Tanaka, *New Integral Equation Approach to Viscoelastic Problems*, Vol. 3: Computational Aspects, Topics in Boundary Element Research, Chapter 2, 25 35, Springer, Berlin, Heidelberg, 1987.

8. W. J. Sim and B. M. Kwak, 'Linear viscoelastic analysis in time domain by boundary element method', *Comput. Struct.*, **29**(4), 531 539 (1988).

9. A. Carini, M. Diligenti and G. Maier, 'Boundary integral equation analysis in linear viscoelasticity: Variational and saddle point formulations', *Comput. Mech.*, **8**, 87 98 (1991).

10. S. Kobayashi and T. Kawakami, 'Application of BE-FE combined method to analysis of dynamic interactions between structures and viscoelastic soil', in C. A. Brebbia and G. Maier (eds), *Boundary Elements VII*, Vol. I, Springer, Berlin, Heidelberg, 6-3 6-12, 1985.

11. G. D. Manolis and D. E. Beskos, 'Dynamic stress concentration studies by boundary integrals and Laplace transform', *Int. J. Numer. Methods Eng.*, **17**, 573 599 (1981).

12. L. Gaul, P. Klein and M. Plenge, 'Dynamic boundary element analysis of foundation slabs on layered soil', in C. A. Brebbia (ed.), *Boundary Elements X*, Vol. 4, Computational Mechanics Publications, Southampton, Boston, 29 44, 1988.

13. G. V. Narayanan and D. E. Beskos, 'Numerical operational methods for time-dependent linear problems', *Int. J. Num. Methods Eng.*, **18**, 1829 1854 (1982).

14. L. Gaul and M. Schanz, 'BEM formulation in time domain for viscoelastic media based on analytical time integration', in C. A. Brebbia, J. Dominguez and F. Paris (eds), *Boundary Elements XIV*, Vol. II, Computational Mechanics Publications, Southampton, 223 234, 1992.

15. M. Schanz and L. Gaul, 'Implementation of viscoelastic behaviour in a time domain boundary element formulation', *Appl. Mech. Rev.*, **46**(11, Part 2), S41 S46 (1993).

16. J. P. Wolf and G. R. Dabre, 'Time-domain boundary element method in viscoelasticity with application to a spherical cavity', *Soil Dyn. Earthq. Eng.*, **5**, 138 148 (1986).

17. L. Gaul and M. Schanz, 'A viscoelastic boundary element formulation in time domain', *Arch. Mech.*, **46**(4), 583 594 (1994).

18. M. Schanz and H. Antes, 'A new visco- and elastodynamic time domain boundary element formulation', *Comput. Mech.*, **20**(5), 452 459 (1997).

19. C. Lubich, 'Convolution quadrature and discretized operational calculus. I', *Numer. Math.*, **52**, 129 145 (1988).

20. C. Lubich, 'Convolution quadrature and discretized operational calculus. II', *Numer. Math.*, **52**, 413 425 (1988).

21. L. Gaul and M. Schanz, 'A comparative study of three boundary element approaches to calculate the transient response of viscoelastic solids with unbounded domains', *Comput. Methods Appl. Mech. Eng.*, (accepted for publication).

22. L. Gaul and M. Schanz, 'Dynamics of viscoelastic solids treated by boundary element approaches in time domain', *Eur. J. Mech. A/Solids*, **13**(4-suppl.), 43 59 (1994).

23. A. T. de Hoop, 'An elastodynamic reciprocity theorem for linear, viscoelastic media', *Appl. Sci. Res.*, **16**, 39 45 (1966).

24. T. A. Cruse and F. J. Rizzo, 'A direct formulation and numerical solution of the general transient elastodynamic problem', *J. Math. Anal. Appl.*, **22**, 244 259 (1968).

25. M. Guiggiani and A. Gigante, 'A general algorithm for multidimensional cauchy principal value integrals in the boundary element method', *J. Appl. Mech.*, **ASME**, **57**, 906 915 (1990).