Abstract

Be-Fe Coupled Method for the Analysis of 3D Elastodynamic Problems in the Time Domain

Francisco Célio de Araújo

Department of Civil Engineering, Escola de Minas/UFOP

CEP: 35400-000, Ouro Preto-MG, Brazil

fcelio@em.ufop.br

Webe João Mansur

COPPE/UFRJ, CXP 68506

CEP: 21945-970, Rio de Janeiro-RJ, Brazil

webe@coc.ufrj.br

Introduction

Nowadays the well-established Finite and Boundary Element Methods are among the most commonly used numerical techniques to treat Engineering problems. These techniques are applicable to solve general real physical problems that occur in Engineering, but is not possible to affirm, in general, which one performs better. According to the problem, each one of the methods can offer advantages or disadvantages.

The Finite Element Method (FEM) is especially attractive because of the good properties of the resulting matrices and of the facility to consider in its formulation non-linear and nonhomogeneous features of the problem definition domain. Analyses with FE models distort however the response of problems defined in infinite or semi-infinite regions, as e.g. in dynamic soil-structure interaction problems. In this case it is not possible to reproduce correctly the radiation of the elastic waves in the soil, which are reflected at the boundaries of the finite elements model. Such problems are more conveniently analysed with the Boundary Element Method (BEM), in which this phenomenon is perfectly accounted for. As a consequence of using boundary integral equations as a starting point for BE formulations one concludes also that the data generation modules in the related computer codes can be implemented in a much easier way than for FE codes. Some disadvantages of this numerical procedure arise however if non-linear problems or problems defined in nonhomogeneous regions are to be solved. Also the determination of the state variables at internal points of domains discretized with boundary elements is rather cumbersome on the numerical view point.

Thus, facing problems as the above mentioned, some engineers and applied mathematicians were led to develop algorithms based on the coupled process. In elastodynamics these works were developed mainly for two- and three-dimensional analyses of frequency-dependent problems (Goto, Matsumoto and Urayama,1985), (Kobayashi and Kawakami, 1985), (Gaitanaros and Karabalis, 1986), (Kobayashi and Mori, 1986). Among the first papers in which coupling formulations for direct time domain BEM analyses are presented are the contributions of Karabalis and Beskos (1985) and Spyrakos and Beskos (1986), who investigated the dynamic response of flexible foundations. Based on the works of these authors von Estorff and Prabucki (1990) developed a general algorithm for the solution of two-dimensional problems in the time domain. Araújo (1994) presented general ideas for developing BE ^_ FE coupled processes for solving 3D wave propagation problems direct in the time domain. Coda and Venturini (1995) presented a BE ^_ FE coupling formulation in which special fundamental solutions are introduced while Araújo, Nishikava and Mansur (1997) considered the so-called enclosing elements in a BE ^_ FE algorithm, which enables then the use of the rigid body displacement criterion for determination of the leading diagonal block matrices (LDBM) of the H-matrices in problems involving semi-infinite domains; later, Araújo, Mansur Nishikava (1998) applied this coupled algorithm for solving problems defined in layered soils. The present work discusses details of a BE-FE coupling algorithm for treatment of general three-dimensional transient problems.

The Boundary Element Formulation

Wave propagation problems in a linear elastic medium W are governed by the Navier equations, which when only inertial forces are considered, can be written as

(1)

where and are respectively the irrotational and equivoluminal wave propagation velocities, r is the material density and b(x,t) represents the body forces. In addition to the partial differential Eqs. (1), which must be satisfied in the interior of W, one must also consider the boundary conditions

(2)

(3)

and the initial conditions

(4)

The elastodynamic problem is completely defined by Eqs. (1) to (4). Alternatively it can be described by the boundary integral equation (Eringen and Suhubi, 1975 and Wheeler and Sternberg, 1968)

(5)

Nomenclature

B = body force vector

c₁ = irrotational wave propagation velocity

c₂ = equivoluminal wave propagation velocity

c_ik = integral free term of the boundary integral

H, G = BE coefficient matrices originated of the traction and displacement fundamental kernels respectively

H = Heaviside function

= effective stiffness matrix

M_a^(m) = time interpolation functions

= elastodynamic traction kernel

= elastostatic traction kernel

P = boundary force field

Q = inverse of A

= effective equivalent nodal load vector

R = radius module

T = time point

= elastodynamic displacement kernel

U = displacement field

x = field point

d_ik = Kronecker's symbol

= stress tensor

= time integration variable

r = material density

W = problem domain

x = source point

B = subscript indicating the BE domain

F = subscript indicating the FE domain

i = subscript indicating the interface

i, k = indexes assuming values 1, 2 and 3

N = Nth time step

n = indicates derivatives in relation to the normal direction

o = subscript indicating the outer boundary

where for 3D problems

(6)

(7)

with and

(8)

The integral free term in (5) is given by

(9)

where are the elastostatic fundamental tractions in G and is a spherical surface centred at x introduced in order to exclude the singular point of the integration surface G (see Fig. 1). With this procedure it becomes clear that the improper boundary integrals indicated in Eq. (5) are actually to be computed in the Cauchy principal value sense. Eq. (9) is correctly demonstrated by Araújo (1994).

Under consideration of homogeneous initial conditions, no volume forces, and by adopting a linear time variation for the field variables, one obtains then from (5) a time discretized integral equation of the form

(10)

where the time interpolation functions are given by

(11)

It is admitted in the development of the time domain BE formulation that the time is subdivided in intervals of equal length Dt and in order to better handle the time variables in Eq. (10), we rewrite this equation in the fashion

(12)

with

(13)

(14)

It should be noticed that N in the expressions (13) and (14) refers to the current time step, m to the time step in which the time integrations are carried out, and a equal 1 and 2 are associated respectively with the beginning and the end of the the time steps.

By using the time translation property of the fundamental displacements and tractions it can be seen that

(15)

(16)

and thus and can be substituted in the expressions above simply by and , where

(17)

(18)

and

(19)

Clearly the time translation property of the fundamental solution together with the choice of uniform time step lengths is of great importance for the efficiency of the time domain BE algorithm. In this way it is possible to express the time discretized boundary integral equation with time integrations restricted to the first time step, and consequently only two new matrices need to be computed for each new time step. Furthermore only one inversion of a coefficient matrix is required in the whole analysis. In the case of using nonuniform time increments the number of matrices to be calculated per time step grows linearly with the number m of time steps and in each of them a new coefficient matrix must be inverted.

For the evaluation of the problem response at the time t_N = ND t only the integrals

(20)

(21)

with a = 1, 2 are therefore necessary, as the other needed time integrals at this time point have already been determined along the previous time steps.

Because of the aspect of the elastodynamic kernels the following three different types of integrals must be considered in Eqs. (20) and (21):

(i) 1^st integration type

(22)

With the aid of the graphic of the Heaviside function in fig. 2, it can be seen that ^(a)I_l exists only within the interval, ( N - 1)c₂Dt < r < Nc₁Dt, which is the region between two spherical surfaces centred in x and with radii ( N ^_ 1)c₂Dt (the internal) and Nc₁Dt (the external).

It follows:

(23)

with

(24)

By making a = 1 and a = 2 in (23) one obtains

(25)

and

(26)

By plotting t_I and t_F versus the radius r (see fig. 3), one can in a simple way determine the integration limits in Eq. (23) for the various possible cases. We conclude from Fig. 3 that the sought limits depend not only on the distance between source and field point but also on the time step number, time step length and wave propagation velocities. The two graphics in Fig. 3 are necessary because not always the inequality ( N ^-1)c₁ > Nc₂ is valid. For normal values of the Poisson's ratio( v = 0.25 ^_0.33}the integration limits for the first two time iterations are obtained from the graphic in fig. 3(a).

In these graphics r_aI and r_aF are give respectively by r_aI = ( N ^-1)c_aD t and r_aF = Nc_aD t , with a = 1, 2.

(ii) 2^nd integration type

(27)

By using the properties of the Dirac's d function one has

(28)

where

(29)

(iii) 3^rd integration type

(30)

where

(31)

The complete expressions for and can be found in Araújo (1994).

By using the 8-nodes isoparametric elements for discretization of the boundary G of the problem domain and after consideration of the prescribed boundary conditions in G the following algebraic equations system is obtained:

(32)

where N_I = max(1, N – N_max + 1). The matrix Eq. (32) furnishes the problem response at the time ^t N.

The Couling Process

Below it is described a time dependent coupling process that consists in inserting the influence of the BE subdomain into the system of FE equations. For development of this process, the continuous medium is subdivided in two subdomains W_BE and W_FE as shown in Fig. 4, which are respectively discretized with boundary and finite elements.

By considering the BE Eq. (32) and the effective FE equilibrium equation obtained after the Newmark integration process (see Bathe, 1982) for the correspondent subdomains in Fig. 4 and under suitable organisation of the state variables it follows:

(i) for the BE subdomain

(33)

where is obtained by reorganisation of the vector that, according to the global nodes numbering from the BE model, is calculated by

(34)

(34)

with N_{I = max (1,N - Nmax + 1).}

(ii) for the FE subdomain

(35)

where

(36)

(37)

with

The coupling conditions on the interface given by

(38)

(39)

can then be introduced into the system by means of the expression for evaluating the equivalent nodal load vector of the FE model on the interface. One has:

(40)

where nfei is the number of finite elements on the interface and it is assumed in (33) that . It results following equations for treatment of the coupled BE-FE system:

(41)

where

(42)

and the submatrices and are determined by inverting the matrix in (33), i.e.,

Displacements in the FE subdomain can be obtained by solving Eq. (41). Other quantities such as velocities, accelerations, support reactions, stresses, etc. in that subdomain can be determined following standard procedures for the Finite Element Method. By observing that u_Bi = u_Fi we can also calculate by means of (33) the boundary unknowns in the region discretized with boundary elements, and then, the internal state variables.

Applications

The developed algorithm is used to analyse the transient response of the problems described next.

Vibration analysis of a cantilever beam

As first example consider the cantilever beam shown in Fig. 5 that is submitted to the given rectangular impulsive shearing load of duration time t₀ = 0.52 [s].

The load acting on the free end cross section is given by:

where p₀ = 10³ [N / m² ]. The used material constants are E = 0.52 x 10⁸ [N / m² ], v = 0.30 and r = 2.00 ´ 10³ [Kg / m³ ]. Time step Dt = 0.72 ´ 10^-2 [s] and the mesh shown in Fig. 6, which is composed of 16 boundary elements and 12 finite elements, are adopted for the coupled time domain FE-BE analysis.

As it can be seen, the numerical solutions agree quite well, although the BE/FE solution after the second vibration period suffers a little damping. Damping can be reduced if a shorter time step is used, but simply by handling the problem simultaneously with two different numerical methods little differences between the BE/FE solution and the FE or BE solution are expected. In fact the coupled algorithm causes the body to behave as non-homogeneous. Thus, surface waves with propagation velocities other than C₁ and C₂ (Timoshenko, 1979) which appear in the neighbourhood of the interface influences the response time history.

Dynamic soil - stiff foundation interaction

The dynamic interaction between the soil - considered as a homogeneous half-space - and a stiff foundation is studied now. The load consists of a concentrated trapezoidal impulse as shown Fig. 8.

The system is characterised by the following physical properties:

In order to simulate the stiff foundation a relatively high elasticity module value is used. The subdomain corresponding to the soil is discretized with 32 boundary elements while in the foundation mesh 32 finite elements are used (see Fig. 9 below).

Two time domain analyses were carried out: considering the foundation having mass and being massless.

case a: foundation with mass

In this case the values characterising the trapezoidal impulse (see Fig. 8) are t₀ = 2.42 ´ 10^-4 [s] and p₀ = 444 [N]. Time step lengths Dt = 0.4033 ´ 10^-4 [s] and Dt = 0.6875 ´ 10^-4 [ s ] were considered in this transient analysis. It is made a comparison of the solution obtained with the BE/FE coupling algorithm developed in this paper and a BEM solution obtained with the BEMFIS algorithm (Latz, 1993). The mesh employed for determination of the BEMFIS solution is more refined in the neighbourhood of the edges of the contact surface than that depicted in Fig. 9. As a consequence, it is more suitable for reproducing the singular behaviour of the tractions on such contact surface than the mesh employed in the BE/FE analysis. Even though, results for vertical displacements for both analyses were quite close (see Fig. 10a). Fig. 10(a) may give the wrong impression that the time step Dt = 0.4033 ´ 10^-4 [ s ] resulted in an amplification of the vertical displacement of the foundation, which is not true. In fact the amplitude is bigger in this case because the impulse on the foundation (see Fig. 8) is I_I = 0.1075[ Ns ], whereas for the other two analyses the impulse is I₂ = 0.0915 [Ns].

The dynamic load is characterised by t₀ = 2.75 ´ 10^-4 [ s ] and p₀ = 444 [ N ]. BE/FE solutions determined with time steps Dt = 0.4583 ´ 10^-4 [ s ] and Dt = 0.6875 ´ 10^-4 [ s ] are again compared against a BEMFIS solution (Latz, 1993). As depicted in Fig. 10(b), time history of the vertical displacement of the stiff foundation for the three analyses agree quite well. It is worth mentioning that the time variation of the impulsive load considered by the BEMFIS computer code is square whereas for the BE/FE analysis it is trapezoidal. If the same time variation had been considered, results of BEMFIS and BE/FE would have been even closer.

Conclusions

The examples analysed in this paper lead to the conclusion that the BE/FE algorithm presented has a good performance. Semi-infinite regions can be modelled in a very simple way and radiation conditions are correctly satisfied, as waves are not reflected at artificial boundaries as happens for finite element models. An accurate response for a soil-foundation interaction problem was obtained for the half-space discretized as a square whose side was 8 times the side of the foundation block. Numerical instability which often occurs in time domain analyses did not appear for the time interval of interest.

Full advantage of the coupled BE/FE procedure presented here can also be taken for problems involving stress concentration, non-linear physical behaviour of materials, non-linear structural behaviour in soil-structure interaction, etc.

Lastly, it should be mentioned that the evaluation of field variables at internal points of BE domains in 3D time domain BEM formulations can be expensive and thus it must be restricted to the minimum. Integral equations to compute internal stresses and velocities have not been presented here, however they must be included in BE codes, as at least stresses are required in most designs. In order to determine internal stresses in a consistent way, by means of a BE time domain algorithm, time interpolation functions for the boundary nodal displacements of at least 2^nd order must be used (Araújo, 1994). But this causes the rank of matrices generated at each time step to be twice the rank of the matrices obtained when linear interpolation functions are considered (Brebbia, Telles and Wrobel, 1984). With the BE/FE coupled algorithm presented here stresses and other field variables can be computed in the FE discretized region in a more suitable way.

Manuscript received: June 1998. Technical Editor: Angela Ourívio Nieckele.

Publication Dates

History