BEM/FEM non-linear model applied to transient analysis with viscous damping

Abstract

In this article a two-dimensional transient boundary element formulation based on the mass matrix approach is discussed. The implicit formulation of the method to deal with elastoplastic analysis is considered, as well as the way to deal with viscous damping effects. The time integration processes are based on the Newmark rhoand Houbolt methods, while the domain integrals for mass, elastoplastic and damping effects are carried out by the well known cell approximation technique. The boundary element algebraic relations are also coupled with finite element frame relations to solve stiffened domains. Some examples to illustrate the accuracy and efficiency of the proposed formulation are also presented.

BEM/FEM Non-Linear Model Applied to Transient Analysis with Viscous Damping

Humberto Breves Coda

Wilson Sérgio Venturini

Departamento de Engenharia de Estruturas

Escola de Engenharia de São Carlos

Universidade de São Paulo

13560-250 São Carlos, SP, Brazil,

[e-mail]

[e-mail]

Abstract

In this article a two-dimensional transient boundary element formulation based on the mass matrix approach is discussed. The implicit formulation of the method to deal with elastoplastic analysis is considered, as well as the way to deal with viscous damping effects. The time integration processes are based on the Newmark r and Houbolt methods, while the domain integrals for mass, elastoplastic and damping effects are carried out by the well known cell approximation technique. The boundary element algebraic relations are also coupled with finite element frame relations to solve stiffened domains. Some examples to illustrate the accuracy and efficiency of the proposed formulation are also presented.

Introduction

In recent years much progress has been made regarding BEM formulations applied to dynamic problems. Many works have been published on the developments of time and frequency domain elastodynamic formulations. As the amount of works already carried out on this subject is rather large, it is convenient to point out some specific reviews, in which the main improvements on that technique are reported. In 1988, Beskos (1988) published a very important review reporting the main developments known so far. Recently the same author published another review(1997) to complete the missing works and to provide the community with the references issued between 1986 and 1997. In the context of the present work the authors wish to point out some relevant works. Coupling techniques, in which BEM equations are combined with algebraic relations coming from other numerical methods, are the natural choice that engineers can follow to solve practical problems. Several works have already shown the applicability and efficiency of this kind of model, e.g. Beskos and Stamos(1995), Beer (1992) and Coda and Venturini (1995). It is also interesting to mention several works that describe mass matrix approach, in which a static fundamental solution is adopted to derive the required integral representations, leading to domain integral terms to take into account the body density, (Wrong, and Hutchinson, 1981; Niwa et al., 1982 and Brebbia, 1985).

Analysis of solids, that exhibit non-linear behaviour subjected to dynamic loads, represent another promising class of problems for which the boundary element methods may be applied. Only recently a few works based on BEM formulations have been published in this context. For most of those works the BEM mass matrix approach has been the scheme adopted, (Telles et al. 1994; Kontoni et al. 1993; Carrer et al. 1991 and Kontoni et al. 1992).

In this paper, the non-linear approach with transient effects is studied again, now including other terms to analyse more complex problems. The static fundamental solution is adopted to derive the integral representations, leading to the standard mass matrix approach for dynamics. The Newmark r and Houbolt schemes are employed to perform the time integrals. Initial strain and stress fields are considered to model non-linear behaviour, viscous damping effects, as well as other temperature like domain loads. Quadrilateral isoparametric cells are adopted to approach all the domain integrals appearing in the integral representations. It is also shown a simple technique to connect finite element bars or a frame structure with the two-dimensional medium modelled by the boundary element. Several examples were chosen to emphasise the performance of the BEM formulation, combined or not with finite elements, to analyse two-dimensional elasto-plastic bodies. The elasto-plastic examples, presented in this paper, were solved by using appropriate implicit and explicit techniques.

Dynamic BEM Mass Approach Formulation

The starting point of the BEM dynamic mass approach for elastic and homogeneous bodies is the generalised Maxwell-Betti reciprocal theorem applied to two independent displacement states {u, s }and {u’, s ’} written as, (Nardini and Brebbia, 1985):

(1)

where u_i represents the displacement component in the i direction, s_ij is the stress tensor and p_i is the surface traction component given by p_i = s_ijn_j, being n the outward normal vector related to the surface G of the body W; Einstein notation is implied.

The dynamic equilibrium equation representing this problem is given by:

(2)

where b_i,and are body force, acceleration and velocity vectors at a field point in W, while mass density and viscous damping are given by r and ¶, respectively.

In order to derive the direct mass BEM formulation a particular known elastostatic state must be adopted as one of those specified states in Eq. (1). The Kelvin’s fundamental solution represents the most common state to be used. This fundamental solution can be derived from Eq. (2), considering an infinite domain W* subjected to a concentrated unit force at a source point s.

Taking into account that usually non-linear BEM formulations are derived by assuming an initial stress state, s_ik^p, applied over the domain W, similar to the ones used to model temperature effects, the body force term b_i has to be replaced by b_i -s_ik^P,k in Eq. (2).

The displacement integral representation, containing mass, viscous damping and initial stress terms can now be derived using appropriately Eqs. (1) and (2), as follows, (Nardini and Brebbia, 1985; Telles and Carrer, 1994; Kontoni and Breskos, 1993; Carrer and Telles, 1991; Kontoni and Breskos, 1992 and Brebbia et al., 1984.

(3)

where superscript (*) represents the adopted Kelvin fundamental values; C_ik is the well known free term, which is equal to d _ik (Kronecker delta) when the source point s is inside the domain, zero for outside source points and depends on the boundary geometry when s is defined along the boundary.

By differentiating Eq. (3) one can write the strain integral representation for internal points, given by,

(4)

where e_ij(s) represents a strain component at the source point s.

In order to obtain the system of time differential equations, appropriate to perform numerical analysis, the boundary of the body is discretized into elements, over which displacements and tractions are approximated. In this work, quadratic isoparametric boundary elements were adopted to transform the boundary integral terms into algebraic ones. The algebraic terms corresponding to the last three domain integrals are computed by adopting isoparametric quadratic cells with eight nodes. The body force integral, the first domain integral, exhibits a known density, therefore can be transformed to the boundary and then evaluated numerically or analytically.

Displacement and strain representations, Eqs. (3) and (4), can be conveniently modified to include internal line loads. Those line loads can represent actual loads distributed along internal lines or the internal forces resulting from the coupling between the two-dimensional body and the embedded reinforcement bars or frame structures. The interactive forces are assumed to act along either one single internal line or along a couple of lines, symmetrically located in relation to the bar middle axis, as shown in Fig. 1. The extra integrals due to the internal line loads are added to Eqs. (3) and (4); they are identical to the surface traction ones, i.e. they are performed similarly to boundary elements.

No displacement approximation is imposed along those internal load lines. As the distance between the two load lines (bar height) is too small in comparison with the internal line element size, an accurate algorithm must be adopted to compute the resulting quasi-singular integrals. In order to carry out accurately those quasi-singular integrals, a progressive element subdivision technique was adopted, (Mon-Ma, 1996).

After performing all integrals along the boundary elements, internal lines and internal cells, the following sets of algebraic Eqs. are achieved from Eqs. (3) and (4) respectively, (Coda et al. 1996)

(5)

(6)

where f represents the body force vector; M and C are the mass and the viscous damping matrices, respectively; while the other vectors and matrices are the standard ones usually employed in BEM formulations.

In this work, two well-known time integration procedures are applied to Eqs. (5) and (6): Newmark b and Houbolt schemes. The main algebraic relations of both procedures are summarised as follows:

Newmark b Non-linear Scheme

Writing the Newmark b scheme for Eq. (5) one can find:

(7)

where

(8)

and

(9)

After imposing boundary conditions and solving the system (7), results:

(10)

where Y_s+1^e is the elastic part of the solution written as:

(11)

while the non-linear part of the solution for displacements is given by:

(12)

In the above algebraic relations, matrix A is the usual BEM matrix when the columns of and bDt²A^-1Q are conveniently interchanged and the superscript ‘BC’ represents prescribed boundary values.

For the Newmark b method, acceleration and velocity representations are expressed in terms of displacements, as follows:

(13)

(14)

(15)

(16)

The values v_s and a_s , known from the previous time step, are linear uncoupled parts of velocity and acceleration, respectively. The remaining parts depend upon the actual displacement values, therefore they cannot be assumed linear as made in the previous procedures found in the literature. The proposed procedure is called here implicit acceleration and velocity evaluation.

Before finding the total solution of Eqs. (5) and (6) let us modify Eq. (6) right hand side by imposing boundary conditions, using the following relation:

(17)

where f'_s+1 and Y_s+1 are the independent vector and the computed displacements, respectively.

Introducing Eq. (10) into Eq. (6) and using expressions (13-17) results:

(18)

where

(19)

(20)

(21)

(22)

in which the superscript ‘d’ indicates that only displacement values should be taken into account, while the superscript ‘u’ indicates the unknown displacement columns to be considered.

It is important noting that the matrices M, C and B are different from each other only by a scalar product when constant density and viscous damping are assumed.

Equation (18) is solved by both explicit and implicit standard procedures, (Telles and Carrer, 1994 and Bonnet, 1995). If matrix W' is zero the process is called explicit in acceleration and velocity, as assumed in references^, (Telles and Carrer, 1995; Kontoni and Beskos, 1993; Carre and Telles, 1991 and Kontoni e Beskos, 1992). The accepted error, less than a prescribed tolerance, is given by:

(23)

where the superscript k represents an iteration inside the time step ‘s+1’.

Houbolt non-linear scheme:

Similar expressions can be written for the Houbolt scheme, as follows:

(24)

(25)

(26)

(27)

(28)

where

(29)

and

(30)

After taking into account the boundary conditions and solving Eq. (28), one finds:

(31)

where

(32)

(33)

The usual BEM matrix A was achieved by conveniently interchanging columns of and Dt²G.

Equation (31) can be introduced into Eq. (6) and them using expressions (17) and (24-27) leads to the same expressions obtained for the Newmark b method, i.e. Eqs. (18-21) with:

(34)

It is worth noting that to compute strains, for boundary nodes, we applied the BEM simplified expressions, (Mon-Ma et al. 1996), therefore the related columns of matrix W' are neglected.

Coupling

Let us consider now a body subjected to dynamic loads. For this problem using the D’Alembert principle, the Virtual Work Principle, one can find, (Clough and Penzien, 1975):

(35)

where e is the strain tensor, upper bars represent virtual displacements, subscripts give directions and superscripts the defined positions.

As usual, from Eq. (35) one can obtain the FEM algebraic system of equations by discretizing the body into elements and approaching appropriately the displacement field. After performing all integrations over the elements, the usual system of time differential equations for dynamic problems is obtained, as follows:

(36)

where K is the stiffness matrix, U and F represent the displacement and nodal force vectors, B gives the equivalent body force nodal values and M and C are the mass and damping matrices, respectively; and are the velocity and acceleration vectors.

The product GP, in Eq. (36), usually called equivalent nodal distributed force values, is here written in its extended form to make possible considering surface forces as unknowns. In this sense, G is called the consistent lumping matrix.

The time step integration schemes used to solve Eq. (36) are the same ones used for the BEM analysis, replacing H by K. In this sense, BEM and FEM algebraic systems of equations can be represented by the same reduced form, as follows:

(37)

The finite element adopted here is the Bernouilli frame element. As usual for that structural element, cubic and linear shape functions are adopted to approach transversal and longitudinal displacements, respectively. Similarly to the BEM approach the interactive forces are assumed linear along the longitudinal and transversal directions. The nodal values are two translations and one rotation at each node. As for the frame element case, stiffness, mass and damping matrices are easily found in the literature, they are omitted here. According to Eq. (36), one has also the lumping matrix, G, which is given by:

(38)

where L is the bar length.

In order to couple two-dimensional BEM equations with FEM frame structure relations, the sub-region technique is followed. Considering two sub-regions defined by W_i and < W_j , coupled together along an interface G_ij, one can apply Eq. (37) for both bodies resulting:

(39a)

(39b)

The equilibrium and geometrical compatibility conditions referred to interface nodes can be written as:

(40a)

(40b)

where the superscripts denote the first and the second sub-regions that define the contact.

The values U^ij and P^ij are interface displacement and contact forces, respectively. The exterior values for both sub-regions, i.e., values that do not belong to the interface, are denoted by U^ie and P^ie. When these values are introduced into Eqs. (39), we have:

(41a)

(41b)

Taking into account Eqs. (40), expressions (41) can be coupled together to give:

(42)

where P^ij denotes prescribed surface forces at the contact surface.

The above expression can be easily generalised for more than two sub-regions. As it can be seen, H and G matrices appearing in expression (42) contain full and null blocks. Several numerical procedures available in the specialised literature can be taken to solve the resulting system of equations.

Examples

In the first example, a finite element beam is coupled with a BEM rectangular domain as shown in Fig. 2. The material properties are the same for both sub-regions and given by the young modulus: E=210 10'kg/(dms²) and the Poisson's ratio:v=0.2. The moment of inertia I=2.25dm⁴ is adopted for the horizontal FEM element, for the vertical ones (or the contact elements) the moment of inertia varies from I=0.25dm⁴ to I=900dm⁴. The adopted discretization is shown in Fig. 2. The contact tractions, shear and normal components, are shown respectively in Figs. 3a -^b.

As it was expected, the contact tractions tend to the simplified Bernoulli hypothesis, as the contact becomes more rigid. For moment of inertia larger than 900 dm⁴, the normal stress profile is almost linear, therefore the results approach the rigid contact case.

The second example shows analysis of a clamped beam subjected to a suddenly applied transversal load, Fig. 4. The material properties are: E=21GPa, v=0.2, r=2500kg/m³ and ¶=100kg/(m³s). The adopted time step is Dt=0.0005s. This analysis is carried out to demonstrate the accuracy of the 2D dynamic BEM mass approach when compared with a simple bar finite element procedure. Moreover, the computed results demonstrated that the BEM solutions are not too sensitive to the mass discretization (see Fig. 5). Coarse domain discretizations give acceptable numerical solutions. The results for vertical displacements at the bar free end are shown in Fig. 6.

In the third example, the same beam is analysed, but inserting a finite element bar inside the domain in order to reinforce the continuous medium. This reinforcement is accomplished by removing the displacement restriction along the tensioned part of the supported bar end, see Fig. 7. The material properties assumed for the embedded bar are: E=210GPa and r=7000kg/m³. The bar moment of inertia for this analysis is neglected. Figure 8 shows the results computed for that reinforced beam together with the numerical solution obtained for the second example, indicated by the caption 'fixed continuum'. In Figure 9, the maximum computed values along the beam length for both, the reinforced and non-reinforced domains, are compared.

Figure 8 shows that as the reinforcement cross section becomes larger the level of displacements is reduced to the fixed continuum case. On the contrary, the frequencies observed are not affected, as expected, when the reinforcement increases. The displacement curves computed along the beam, exhibited in Fig. 9, confirm that the frequencies should not be the same. A larger displacement slope at the supported beam end is verified for the reinforced case. The reinforced beam exhibits more rigid shape, leading to show more mass influence than for the continuum case, resulting into a lower vibration frequency.

In the next example, the simply supported deep beam, solved by Carrer and Telles (1991) and Kontoni and Beskos (1994) using explicit acceleration and velocity evaluations, is studied again adopting the implicit scheme. That simply supported deep beam is subjected to a suddenly applied uniformly distributed load. The material is considered to be elastic perfectly plastic governed by the Von Mises yield criterion. The following physical parameters are adopted: E=100; s_y=0.16;n=0.333;r=1.5

The geometry and the adopted discretization are shown in Fig. 10. The uniformly distributed load is given by the value p=0.01. First, the problem is solved using the Houbolt scheme, for which we have adopted the time step Dt=0.56 and the tolerance tol=1%.

The values computed for this load level are practically the same ones found in references Kontoni and Beskos (1993) and Carrer and Telles (1991), therefore validating our code. The numerical solution obtained by using the implicit acceleration scheme does not show significant differences at this load level. At that level only few points of the solid have yielded, reducing therefore possible differences between both procedures.

Applying a rather large load (p=0.012) to increase the amount of yielding, a considerable difference between using the explicit and implicit acceleration evaluation schemes is observed, as shown in Fig. 12.

As it can be seen, Fig. 12, the implicit evaluation of acceleration reduces the amount of yielding. This behaviour is in agreement with expression (19). Higher load levels can be considered if the discretization is refined. Other complementary results obtained by using the Newmark b scheme can be taken from Figs. 13 and 14.

From Figure 13, one can see that, for the same level of loads, the Houbolt scheme leads to larger displacements in comparison with the Newmark b method. As a consequence, the amount of developed yielding is higher for the Houbolt scheme, by comparing Figs. 11 and 14, what seems to be more reliable for structural analysis. One can observe, as well, that the Newmark b method together with BEM is unable to solve high frequency problems. Thus, due to these two reasons, we can say that the Houbolt integration scheme is the appropriate method to be used together with BEM mass matrix for dynamic analysis.

In Figure 15, elastoplastic results computed using the Houbolt scheme with implicit acceleration forp=0.012 are shown, assuming the viscous damping ¶=0.5; 0.2; 0.1; 0.05; 0.01 and 0..

The same beam is now analysed assuming a vertical base excitation defined by making u_v=0.04sin(0.18t). The vertical displacements, computed at point A, assuming elastic and elastoplastic behaviours, are given in Fig. 16.

Figure 16

The last example was chosen particularly to illustrate the elastoplastic algorithm and the capability of the formulation to deal with higher frequency modes using the Houbolt scheme. A thick cylinder geometrically defined by internal and external radii, a=100mm and b=200mm, is subjected to a suddenly applied and maintained internal pressure, (Kontoni and Beskos, 1993). Perfect plasticity and plane stress conditions are assumed, together with the Von Mises yield criterion. In order to preserve the proper time dependence of the problem, the material characteristics and internal pressure p should be conveniently given as follows:

The tolerance adopted to solve this problem is tol=10^-5. Fig. 17 shows the discretization using quadratic isoparametric cells and boundary elements; double symmetry properties have been adopted as well, therefore no boundary element is required along radial directions. The radial displacements at point A(r, q)=A(100,45°) versus time, shown in Fig. 18, were obtained by following the implicit acceleration scheme, named DI-BEM for simplicity, and the non accelerated formulation, named NA-BEM (the acceleration term is neglecting for strain evaluation). The time step dependence of the first scheme is illustrated in Fig. 19, where several values have been tested showing the convergence behaviour.

The numerical solution obtained here, for this simple example, is in good agreement with the ones presented by Kontoni and Beskos (1993)using the dual reciprocity boundary elements (DR-BEM) and finite elements. Unfortunately, Kontoni and Beskos (1993) has not used smaller time step to capture higher frequencies. This has been provided recently by Frangi (1998), who has also solved the same example using the symmetric Galerkin BEM, for which the following parameters were adopted:

In order to compare Frangi's results with the solutions obtained by using the proposed scheme, a new discretization (Fig. 20) is adopted together with the time step Dt=3x10^-7s. For this comparison we took only the radial displacements computed at points A(r, q)=A(100,45°) and B(r, q)=B(200,45°).

Frangi has adopted a rather smaller time step to analyse this problem in comparison with the one taken by Kontoni and Beskos (1993). Due to this reason, Fig. 21 shows much more wave reflection details than one can see in Figs. 18 and 19. This example confirms once more the good agreement between the presented formulation and previous published SG-BEM (using 32 cells) and FEM (ABAQUS with 300 quadrilateral bilinear elements) solutions.

Conclusions

A general 2D BEM mass matrix technique applied to elastoplastic dynamic analysis of bodies has been successfully implemented. An approach to evaluate acceleration and velocity in an implicit way along the time increment scheme was developed. The coupling of a finite element bar with the 2D body modelled by the boundary element method was included in the formulation, as well. In addition, the BEM formulation was successfully extended to treat damping effects. The examples presented in this article demonstrate that the proposed formulations are general and useful for practical purposes. The numerical results obtained for elastoplastic and damping analyses show clearly the accuracy of the developed algorithms when compared with other BEM mass matrix formulations and finite element approaches.

Manuscript received: September 1998. Technical Editor: Hans Ingo Weber.

• Beer, G. and Watson, J. O., 1992 "Introduction to Finite and Boundary Element Methods for Engineers". John Wiley and Sons, New York.
• Beskos, D. E., 1988 "Boundary Element Methods in Dynamic Analysis". App. Mech. Rev. Vol. 40, p 1-23.
• Beskos, D. E., 1997 "Boundary Element Methods in Dynamic Analysis: Part II (1986-1996)". App. Mech. Rev., 50, p.149-197,
• Bonnet, M., 1995, "Équations intégrales et éléments de frontière. Sciences et Techniques de l'ingénieur", CNRS Éditions/Eyrolles.
• Brebbia, C.A.; Telles, J.C.F. and Wrobel, L.C., 1984, "Boundary Element Techniques. Springer Verlag, Berlin".
• Carrer J.A.M and Telles J.C.F., 1991, "Transient Dynamic Elastoplastic Analysis by the Boundary Element Method". In: Boundary Element Technology VI, Proceedings, ed. C.A. Brebbia. CMP, UK, and Elsevier, p.265-277.
• Clough, R.W. and Penzien, J., 1975, "Dynamics of structures". McGraw-Hill.
• Coda,H.B.; Venturini W.S and Aliabadi, M.H., 1996, "A Simple Coupling of 2D BEM and FEM bar Model Applied to Mass Matrix Elastodynamic Analysis". In: Boundary Elements XVIII, Proceedings, Brebbia et. al eds.,CMP. UK, p.363-372.
• Coda. H. B.and Venturini, W. S., 1995, "Three Dimensional Transient BEM Analysis". Computer of Structures, 56, n. 5, p.751-768, Pergamon.
• Frangi A., 1998, "Some Developments in the Symmetric Galerkin Boundary Element Method". PhD Theses, Politecnico di Milano, Milan.
• Kontoni, D.P.N. and Beskos D.E., 1992, "The Dual Reciprocity Boundary Element Method for the Transient Dynamic Analysis of Elastoplastic Problems". In: Boundary Element Technology VII, Proceedings, ed. C.A. Brebbia. CMP, UK, and Elsevier, p.259-272.
• Kontoni, D.P.N. and Beskos, D.E, 1993, "Transient Dynamic Elastoplastic Analysis by the Dual Reciprocity BEM". Engineering Analysis with Boundary Elements, 12, p.1-16, 1993.
• Mon-Ma, M.L.; Venturini, W.S. and Coda, H.B., 1996, "A Simple Technique to Evaluate Quasi-Singular Integrals". In: First Brazilian Seminar on the Boundary Element Method in Engineering, p. 44-53, Federal University of Rio de Janeiro.
• Nardini, D. and Brebbia, C.A., 1985, "Boundary Integral Formulation of Mass Matrices for Dynamic Analysis". In: Topics in Boundary Element Research, Ed. Brebbia, C.A., V2.
• Niwa, Y, Kobayashi and Kitarrara, M., 1982, "Applications of the Bem to Eigenvalue Problems in Eladtodynamics". In Boundary Element Methods in Engineering. Ed. C.A. Brebbia, Springer Verlag, Berlin.
• Stamos, A.A. and Beskos, D.E., 1995, "Dynamic Analysis of Large 3D Underground Structures by the BEM". Earthquake engineering and structural dynamics, 24, p.917-934.
• Telles J.C.F.and Carrer J.A.M, 1994, "Static and Dynamic Non-linear Stress Analysis by the Boundary Element Method With Implicit Techniques". Engineering Analysis with Boundary Elements, 14, p.65-74.
• Wong, G.I.K and Hutchinson, J.R., 1981, "An Improved Boundary Element Method for Plate Vibrations". In Boundary Element Methods. Ed. C.A., Brebbia, Springer Verlag, Berlin.

Publication Dates

• Publication in this collection
  20 Nov 2002
• Date of issue
  Sept 1999

History

• Received
  Sept 1998