Inelastic 2D analysis by adaptive iterative BEM–FEM coupling procedures

Computers and Structures 156 (2015) 134–148

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Inelastic 2D analysis by adaptive iterative BEM–FEM coupling procedures D. Soares a,⇑, L. Godinho b a b

Structural Engineering Department, Federal University of Juiz de Fora, CEP 36036-330 Juiz de Fora, MG, Brazil ISISE, Department of Civil Engineering, University of Coimbra, 3030-788 Coimbra, Portugal

a r t i c l e

i n f o

Article history: Received 20 January 2015 Accepted 5 May 2015 Available online 25 May 2015 Keywords: Iterative coupling Boundary elements Adaptive finite elements Elastoplasticity Non-matching discretizations Optimal relaxation parameters

a b s t r a c t In this work, two iterative BEM–FEM coupling procedures are discussed, taking into account nonlinear models and adaptive spatial discretizations. In both procedures, optimal relaxation parameters are computed (speeding up the convergence of the iterative coupling) and non-matching discretizations at common interfaces are allowed (enabling more versatile adaptive approaches). A single unified iterative loop is considered in order to deal with all the focused iterative solutions simultaneously (i.e., nonlinear analysis, adaptive analysis and coupling analysis), rendering a very efficient methodology. At the end of the paper, numerical results are presented, illustrating the good performance of the proposed technique. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The numerical simulation of arbitrarily shaped/loaded continuous linear/nonlinear mechanical models remains, despite much effort and progress over the last decades, a challenging area of research. In most cases, discrete techniques, such as the Boundary Element Method (BEM) and the Finite Element Method (FEM), have been employed and continuously further developed with respect to accuracy and efficiency. Both the BEM and the FEM have relative benefits and limitations. The Finite Element Method, for instance, is well suited for inhomogeneous and anisotropic materials as well as for dealing with nonlinear behavior. For systems with infinite extension, discontinuities and regions of high stress concentration, however, the use of the Boundary Element Method is by far more advantageous (for a deeper discussion about the advantages and disadvantages of the BEM, reference [1] is indicated). In fact, it did not take long until some researchers started to combine the BEM and the FEM in order to profit from their respective advantages, trying to evade their disadvantages, and nowadays several works dealing with BEM–FEM coupling are available [2–20]. However, standard coupling procedures of BEM/FEM can lead to several problems with respect to efficiency, ⇑ Corresponding author. E-mail addresses: delfi[email protected] (D. Soares), [email protected] (L. Godinho). http://dx.doi.org/10.1016/j.compstruc.2015.05.007 0045-7949/Ó 2015 Elsevier Ltd. All rights reserved.

accuracy and flexibility. First, the coupled system of equations has a banded symmetric structure only in the FEM part, while in the BEM part it is non-symmetric and fully populated; consequently, for its solution, the optimized solvers usually employed with the FEM can no longer be used, leading to rather expensive calculations with respect to computer time. Second, quite different physical properties may be involved in the coupled model, resulting in ill-conditioned matrices when standard coupling procedures are considered; this may affect the accuracy of the methodology, providing misleading results. Third, the standard coupling methodology does not allow independent discretization for each subdomain of the model, requiring matching nodes at common interfaces, which drastically affects the flexibility and versatility of the technique. Moreover, the efficiency of the technique may be further reduced if nonlinear models are considered. In this case, the larger and overly complex coupled system of equations must be dealt with several times within a single analysis, leading to excessive computational costs. In order to evade these drawbacks, iterative coupling procedures have been developed. Initially, static problems were studied considering iterative coupling approaches, and linear and nonlinear behavior have been simulated [5–11]. Later on, dynamic problems were focused, and time [12–15] and frequency domain [16–19] iterative analyses have been implemented (in this case, an overview is presented by Soares and Godinho [20]). Iterative coupling approaches allow BEM and FEM subdomains to be analyzed separately, leading to smaller and

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

better-conditioned systems of equations (different solvers, suitable for each subdomain, may be employed). Additionally, a small number of iterations is required for the algorithm to converge and the matrices related to the smaller governing systems of equations do not need to be treated (inverted, triangularized etc.) at each iterative step (unless it is necessary, as for instance for some nonlinear analysis etc.), providing an efficient methodology. This coupling technique also allows independent discretizations to be easily and efficiently employed for the boundary and finite element subdomains, without any requirement of matching nodes at common interfaces. This is particularly important when some sort of adaptive discretization technique is considered within one of the subdomains, allowing non matching nodes to occur at the common interfaces. In the present work, adaptive iterative BEM–FEM coupling procedures are proposed, in order to analyze nonlinear static models. More precisely, elastoplastic problems are focused here. In this context, an adaptive FEM approach is applied to the regions where the inelastic behavior is expected to occur, whereas the BEM is applied to the regions where the elastic behavior remains. Iterative BEM–FEM coupling techniques applied to elastoplastic problems have already been reported in the literature (see, for instance, the work of Elleithy et al. [7], Boumaiza and Aour [11], etc.); however, in most of these references, two iterative loops are considered, being the complete iterative loop of the nonlinear FEM analysis carried out within an iterative step of the iterative loop of the BEM–FEM coupling. This multiple iterative procedure is computationally very demanding. In the present work, a single iterative loop is considered, in which three iterative analyses are carried out together, at a common iterative step, namely: (i) the FEM nonlinear analysis; (ii) the FEM adaptive discretization; (iii) the BEM–FEM coupling. As it is reported in Section 5, this unified iterative approach does not significantly increase the number of iterations that would be required by the dominant iterative analysis of the model when acting isolated, rendering a very effective technique. Adaptive BEM–FEM coupled analyses, taking into account elastoplastic models, have been reported by Elleithy et al. [8,9]. In these works, an interesting adaptive procedure is presented, in which the FEM subdomains are expanded (and BEM subdomains are shrunk) as the plastic zones evolve. Here, a very different approach is considered. In the present work, the regions modeled by the BEM do not change along the analysis, and the matrices of the BEM subdomains are computed (and treated) just once, enhancing efficiency. In addition, a coarse discretization is initially adopted for the FEM subdomains, and this discretization is adaptively enriched as the solution evolves. Thus, at the end of the process, a refined FEM discretization is expected to be available at the regions where the plastic zones occur, providing an optimal FEM simulation linked to a very efficient BEM analysis. As one can notice, following this adaptive approach, it is essential to allow considering non-matching nodes at BEM–FEM common interfaces, otherwise the BEM discretizations would also have to adapt, increasing the computational cost of the analysis. Two iterative coupling approaches are discussed here. In these approaches, either displacements or tractions may be considered prescribed to the BEM common interface (and, as a consequence, either tractions or displacements may be considered prescribed to the FEM common interface, respectively). Once these two approaches are available, their use can be decided in accordance to the characteristics of the model (thus, avoiding singular systems of equations to be obtained due to the lack of essential boundary conditions on a subdomain, for instance). For both coupling approaches, optimal relaxation parameters are employed to speed up and/or to ensure the convergence of the iterative coupling

135

analysis. As it has been previously remarked [5–20], this procedure is highly important to guarantee the effectiveness of the technique. The paper is organized as follows: first, the governing equations of the elastoplastic model are briefly presented, and basic aspects of the BEM and FEM are described; in the sequence, the iterative BEM–FEM coupling algorithm is discussed, and numerical applications are considered, illustrating the versatility and good performance of the proposed techniques. 2. Governing equations The basic governing equations related to elastoplastic materials are given by:

rij;j ¼ ci

ð1Þ

drij ¼ Dep ijkl dekl

ð2Þ

deij ¼

1 ðdui;j þ duj;i Þ 2

ð3Þ

where Eq. (1) is the equilibrium equation and Eqs. (2) and (3) stand for incremental relations. The Cauchy stress, using the usual indicial notation for Cartesian axes, is represented by rij , and ui and ci stand for displacement and body force distribution components, respectively (inferior commas indicate partial space derivatives and the summation convention is considered, with the indices varying from 1 to 2 for 2D analyses, and from 1 to 3 for 3D analyses). Eq. (2) is the constitutive law, written incrementally. The incremental strain components deij are defined in the usual way from the displacements, as described by Eq. (3). In Eq. (2), Dep ijkl is a tangential tensor defined by suitable state variables and the direction of the increment. Within the context of associated isotropic work hardening theory [21,22], the tangent constitutive tensor is defined as:

Dep ijkl ¼ Dijkl  ð1=wÞDijmn amn aop Dopkl

ð4Þ

where

Dijkl ¼ 2lm=ð1  2mÞdij dkl þ lðdik djl þ dil djk Þ

ð5aÞ

 =@ rkl akl ¼ @ r

ð5bÞ

w ¼ aij Dijkl akl þ H

ð5cÞ

H ¼ @ r0 =@ ep

ð5dÞ

 and e are the equivalent (or effective) stress and plastic In Eq. (5), r strain, respectively; r0 is the uniaxial yield stress; H is the plastic-hardening modulus (for the case of a perfectly plastic material H ¼ 0); l and m stand for the shear modulus and the Poisson ratio, respectively; and dik is the Kronecker delta. In case of elastic analyses, the Cauchy stresses can be defined by rij ¼ Dijkl ekl , where Dijkl (see Eq. (5a)) is the elastic constitutive tensor (this linear relation is a particular case of Eq. (2)). In addition to Eqs. (1)–(5), boundary conditions are prescribed as follows, in order to completely define the problem: p

 i at Cu ui ¼ u

ð6aÞ

si ¼ rij nj ¼ si at Cs

ð6bÞ

where the prescribed values are indicated by overbars and si stands for traction components along the boundary whose unit outward normal vector is represented by ni . Following Eq. (6), the boundary of the model (C) is divided into an essential (Cu ) and a natural (Cs ) boundary, where Cu [ Cs ¼ C and Cu \ Cs ¼ 0.

136

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

In the next section, boundary and finite element techniques are briefly described. Here, elastoplastic regions are treated by the FEM, whereas elastic subdomains may be analyzed by the BEM or by the FEM.

For elastoplastic models, the vectors and matrices represented in Eq. (10a) may be described as follows, considering each finite element e: nþ1

fk

¼

[Z

In the present work, an iterative algorithm is considered to analyze coupled BEM–FEM elastoplastic models. In this iterative algorithm, nonlinear adaptive finite element analysis and BEM–FEM coupling are carried out together within the same iterative loop. Thus, although in this section BEM and FEM formulations are described independently (i.e., not considering a coupled BEM–FEM analysis), the variables described here follow the generic notation unk . In this notation, the superscript n stands for the current incremental step of the analysis and the subscript k stands for its iterative step. Considering this notation and assuming a boundary element discretization, elastostatic models may be analyzed by the following system of equations [23,24]: nþ1

nþ1 ðC þ HÞunþ1 kþ1 ¼ Gskþ1 þ dkþ1

ð7Þ

where u and s stand for displacement and traction vectors, respectively; C stands for the location matrix and G and H are influence matrices. Vector d accounts for domain terms, such as body force contributions. Matrices G and H can be computed taking into account each boundary element e of the model, and vector d can be evaluated considering each domain integration cell c:

H¼

[Z [Z e

nþ1 dkþ1

¼

ð8aÞ



U NdC

ð8bÞ

Ce

[Z c



nþ1 kþ1 d

Uc Xc



X

ð8cÞ



In Eq. (8), T and U stand for traction and displacement fundamental matrices, respectively, and N represents the adopted BEM interpolation matrix. By re-arranging the system of Eq. (7), taking into account the boundary conditions of the problem, just known and unknown variables can be disposed at the right and at the left hand side of the system of equations, respectively, allowing its solution. Eqs. (7–8) only intend to summarily describe the boundary element formulation considered here; for further details on the topic, the following references are suggested [23,24]. Considering a finite element discretization, the basic equation describing a nonlinear model is given by: nþ1

f ðunþ1 kþ1 Þ ¼ f kþ1

ð9Þ

where f ðuÞ and f stand for internal and external forces, respectively. Taking into account a linearized incremental approach, Eq. (9) can be rewritten as [25,26]: nþ1

Kk dukþ1 ¼ f k

 f ðunþ1 k Þ

nþ1 unþ1 þ dukþ1 kþ1 ¼ uk

ð10aÞ ð10bÞ

where Kk stands for the tangent nonlinear stiffness matrix, nþ1

[Z

f ðunþ1 k Þ ¼

Xe

e

Kk ¼

[Z e

Xe

Z Cse

nþ1 NT s k dC



BT rnþ1 k dX

BT Dep k BdX

ð11aÞ

ð11bÞ

ð11cÞ

where B stands for the strain matrix and Dep k is the nonlinear constitutive matrix (as described by Eq. (4)). N stands for the adopted FEM interpolation matrix. The stress state in Eq. (11b) can be evaluated ¼ rn þ drnþ1 ¼ following Eqs. (2) and (3); i.e., by considering rnþ1 k k nþ1 nþ1 rn þ Dep ¼ rn þ Dep , where dunþ1 ¼ unþ1  un . k k k dek k Bduk

If adaptive discretizations are considered, at each stage that the FEM mesh is actualized, the fields of the model (i.e., displacements, stresses etc.) have to be updated, taking into account the new element and node distributions. This updating can be carried out taking into account the adopted interpolation functions of the model. Eqs. (9–11) briefly describe the finite element formulation considered here; for further details on the topic, the following references are suggested [25–27]. Once Eqs. (7–11) are considered, subdomains modeled by the BEM or by the FEM can be analyzed. In the next section, the iterative coupling of these subdomains is discussed. 4. Coupled BEM–FEM analysis

Ce

e

G¼

T NdC

Xe

e

3. Boundary and finite element methods

NT cnþ1 k dX þ

Þ represents the nonlinear residual vector and dukþ1 is f k  f ðunþ1 k the variation of the incremental displacements, calculated at each iterative step.

For the coupled analysis in focus, the following continuity and equilibrium equations must hold at the interfaces between the BEM and FEM subdomains: F

ui ¼ B ui

ð12aÞ

F

si þ B si ¼ 0

ð12bÞ

where the superscripts B and F indicate if a variable is related to the BEM or to the FEM subdomain, respectively. In the present work, two different configurations can be considered, regarding the BEM–FEM interface: (i) in the first configuration (configuration 1), the common interface is considered as an essential boundary for the BEM, with prescribed displacements, and as a natural boundary for the FEM, with prescribed tractions; (ii) in the second configuration (configuration 2), the common interface is considered as a natural boundary for the BEM, with prescribed tractions, and as an essential boundary for the FEM, with prescribed displacements. For both cases, an iterative renewal of the variables at the common interface is carried out, taking into account Eq. (12). In the BEM–FEM iterative coupling algorithm considered here, the FEM subdomain is analyzed first, and the FEM displacements (F unþ1 ) or tractions (F snþ1 ) at the common interface are evaluated kþ1 kþ1 (in this case, tractions are computed as projections of the stress state, as indicated by Eq. (6b)). In the sequence, these values are applied to the BEM as prescribed boundary conditions, following Eq. (12) and the configuration in focus (i.e., configuration 1 or 2) for the common interface. This is carried out taking into account a relaxation parameter a, which is introduced here in order to ensure and/or to speed up convergence. Thus, the following variables may be computed: B  nþ1 ukþ1

F nþ1 B  nþ1 ¼ ðaÞIBF u ð ukþ1 Þ þ ð1  aÞ uk

ð13aÞ

137

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148 B nþ1 kþ1

s

F nþ1 B nþ1 ¼ ðaÞIBF s ð skþ1 Þ þ ð1  aÞ s k

ð13bÞ

where I stands for spatial interpolation functions; Eqs. (13a) and (13b) are related to configuration 1 and 2, respectively. In the present work, non-matching BEM and FEM nodes may be considered (this is particularly important when adaptive discretizations are employed). Thus, a routine to spatially link the FEM and BEM computed values must be taken into account. In IBF u

the FEM as prescribed boundary conditions, following Eq. (12). Thus, in order to obtain the FEM prescribed values, the following variables may be computed: F nþ1 kþ1

B nþ1 ¼ IFB s ð skþ1 Þ

ð14aÞ

F  nþ1 ukþ1

B nþ1 ¼ IFB u ð ukþ1 Þ

ð14bÞ

s

BF

Eq. (13), and Is link the displacements and tractions computed at the finite element nodes and faces, respectively, to their respective values at the boundary element functional nodes. These interpolation routines can be carried out based on the interpolation functions of the elements, or considering enriched approaches. In Fig. 1, for instance, a sketch illustrating linear interpolation procedures is depicted, indicating how a variable ‘‘v’’ at a common interface can be computed from the correspondent adjacent variables. Once the prescribed variables at the BEM common interface are computed, the BEM subdomains can be analyzed, allowing to comB nþ1 pute displacements (B unþ1 kþ1 , configuration 2) or tractions ( skþ1 , configuration 1) at the common interface. These values are applied to

vi di

v = I(vk)

dj vj

Fig. 1. Spatial interpolation procedures: example of an interpolation scheme of values in order to obtain v (linear interpolation: v ¼ ðv i dj þ v j di Þ=ðdj þ di Þ).

vk

where once again, I stands for spatial interpolation functions and Eqs. (14a) and (14b) are related to configuration 1 and 2, respecFB tively. In Eq. (14), IFB u and I s link the displacements and tractions computed at the boundary element functional nodes, respectively, to their respective values at the finite element nodes. Once the FEM prescribed boundary conditions are computed, the iterative cycle is reinitiated, being all the above described procedures repeated until convergence is achieved. It is important to highlight that just one iterative loop is considered here to deal with several types of iterative approaches, namely: (i) FEM nonlinear analysis; (ii) FEM adaptive discretization; (iii) BEM–FEM coupling. As it is illustrated in the next section, this seems to be an appropriate approach, being the number of iterative steps required by the isolated dominant iterative procedure not significantly increased by considering all the iterating procedures together, in the same iterative loop. The usual alternative to the single iterative approach is a multiple iterative algorithm. In this case, an entire iterative loop is carried out within each iterative step of the ‘‘host’’ iterative loop. Thus, the multiple iterative procedure may considerably increase the computational effort of the analysis, once too many iterative steps occur and/or iterative processes are considered. In Fig. 2, sketches for the multiple and the single iterative algorithms are depicted, taking into account the application in focus. As one can observe, not only the single iterative algorithm is a more efficient approach

Multiple iterative algorithm Single iterative algorithm Coupling iterative loop

Unified iterative loop Adaptive iterative loop FEM analysis

Nonlinear iterative loop FEM analysis

Check for convergence Interface treatment

Interface treatment

Check for convergence

BEM analysis

BEM analysis

Interface treatment

Interface treatment Check for convergence

Check for convergence

Fig. 2. Sketches for the multiple and the single iterative algorithms.

138

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

(as it is further described in the next section), but also it is more simple to implement. It is important to observe that the effectiveness of the present iterative coupling algorithm is intimately related to the relaxation parameter selection. An inappropriate selection for a can drastically increase the number of iterations in the analysis or, even worse, make convergence unfeasible. In the next sub-section, an expression for an optimal relaxation parameter is provided. 4.1. Optimal relaxation parameter In order to evaluate an optimal relaxation parameter, the following square error functional is here minimized: 2 B  nþ1  nþ1 f ðaÞ ¼ jjB v kþ1 ðaÞ  v k ðaÞjj

ð15Þ

where vector v stands for the traction or displacement vector, according to the configuration in focus. Taking into account the relaxation of vector v for the k þ 1 and k iterations, Eq. (16) may be written, regarding Eq. (13):

v

F nþ1 B  nþ1 ¼ ðaÞIBF v ð v kþ1 Þ þ ð1  aÞ v k

ð16aÞ

v

F nþ1 B  nþ1 ¼ ðaÞIBF v ð v k Þ þ ð1  aÞ v k1

ð16bÞ

B  nþ1 kþ1 B  nþ1 k

Substituting Eq. (16) into Eq. (15) yields:

f ðaÞ ¼ jja F w þ ð1  aÞB wjj2 ¼ ða2 jjF wjj2 þ 2að1  aÞðF w; B wÞ þ ð1  aÞ2 jjB wjj2 Þ where

the

inner

product

definition

is

employed

ð17Þ (e.g.,

ðw; wÞ ¼ jjwjj2 ) and new variables, as defined in Eq. (18), are considered. F

BF F nþ1 F nþ1 w ¼ IBF v ð v kþ1 Þ  Iv ð v k Þ

ð18aÞ

B

 nþ1  nþ1 w ¼ Bv  Bv k k1

ð18bÞ

To find the optimal a that minimizes the functional f ðaÞ, Eq. (17) is differentiated with respect to a and the result is set to zero, as described below:

ajjF wjj2 þ ð1  2aÞðF w; B wÞ þ ða  1ÞjjB wjj2 ¼ 0

ð19Þ

Re-arranging the terms in Eq. (19), yields:

a ¼ ðB w; B w  F wÞ=jjB w  F wjj2

ð20Þ

which is an easy to implement expression that provides an optimal value for the relaxation parameter a, at each iterative step. It is important to note that the relation 0 < a 6 1 should hold. In the present work, the optimal relaxation parameter is evaluated according to Eq. (20) and if a R ½0:01; 1:00 the previous iterative-step relaxation parameter is adopted. For the first iterative step, a ¼ 0:5 is selected.

5. Numerical applications In this section, three numerical examples are presented, illustrating the performance and potentialities of the discussed techniques. In the first example, a circular cavity is considered, and BEM–FEM coupled results are presented, taking into account a Mohr–Coulomb criterion and adaptive and fixed FEM discretizations. Analogously, in the second example, a half-space model is analyzed. In this case, the Drucker–Prager yield criterion is considered and adaptive and fixed BEM–FEM discretizations are carried out. Finally, in the third example, a cantilever beam is studied, following the von Mises yield criterion. Here, adaptive BEM–FEM coupled results are compared to those provided by an adaptive and fixed FEM approach. In the first two examples, infinite domain models are focused. In these cases, the so-called configuration 2 is considered, and prescribed tractions are applied at the BEM common interface, whereas prescribed displacements are applied at the FEM common interface. In the third application, a finite domain model is analyzed and configuration 1 is followed, considering prescribed displacement and tractions applied to the BEM and FEM common interfaces, respectively. For this third application, linear boundary elements are considered, whereas, for the other applications, constant boundary elements are employed. For all the analyses that follow, linear triangular finite elements are adopted, since discretizations considering this type of finite element are easier to adaptively refine. Adaptive finite element methods are nowadays widely used and a typical loop for an adaptive FEM through local refinement basically involves 4 steps, namely: (i) solve; (ii) estimate; (iii) mark; (iv) refine/coarsen. In this sense, FEM subdomains are here analyzed and the solution on the current triangulation is obtained (first step). The error is then estimated using the computed solution (second step) and it is used to mark (third step) a set of triangles that are to be refined, which is carried out keeping the triangulation shape regularity and conformity (fourth step). The adaptive procedure implemented here is based on the package provided by Chen and Zhang [27]. In this work, the FEM discretization is adapted just at the first iterative step of each incremental step. Moreover, a Newton–Raphson initial stress configuration is considered here for the nonlinear analysis, and a tight relative error tolerance (displacement and force residual) of 105 is adopted for the convergence of the iterative process. SI units are considered in the subsections that follow. 5.1. Circular cavity This plane strain problem consists of a circular cavity under a uniform internal pressure distribution. A sketch of the model is depicted in Fig. 3. The geometry of the problem is defined by

L

BEM

H

FEM

Fig. 3. Sketch of the cavity model and initial FEM discretization for the adaptive analyses.

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

L ¼ 2 m and H ¼ 4 m. The physical properties of the model are E ¼ 108 N=m2 (Young modulus) and m ¼ 0:3 (Poisson coefficient). A perfectly plastic material obeying the Mohr–Coulomb yield criterion is assumed, where c ¼ 7  102 N=m2 (cohesion) and / ¼ 20 (internal friction angle). Two different BEM–FEM discretizations are considered here. In the first discretization, a fixed FEM mesh is employed, which is composed of 2880 elements. In the second discretization, an adaptive FEM mesh is considered, being its initial

139

configuration depicted in Fig. 3 (320 elements). For both cases, an invariable BEM discretization with 64 constant elements is adopted. For the analyses in focus, 6 incremental steps are considered. In Fig. 4, displacement results are depicted considering both discretization procedures and linear and nonlinear analyses. The number of the incremental step, the number of iterative steps necessary for that incremental step, and the number of nodes and

Fig. 4. Computed displacements (modulus) along the FEM deformed meshes (scale factor of 104) for the cavity model: (a and b) elastic analysis; (c and d) elastoplastic analysis; (a) and (c) fixed FEM discretization; (b) and (d) adaptive FEM discretization.

140

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

elements in the current FEM discretization are also depicted in the figure. As one can observe, good agreement is observed among the results. In Fig. 5, the curves of the relaxation parameters computed along each incremental step of the analyses are depicted. As one can observe, the relaxation parameter is characterized by an intricate evolution. In fact, since expression (20) is based on residuals (see Eq. (18)), a very dynamic behavior is expected for this parameter. The evolution of the BEM–FEM adaptive elastoplastic analysis is depicted in Fig. 6. In this figure, the evolution of the FEM discretizations and the computed equivalent plastic strains, at the end of each incremental step, are illustrated. In Table 1, the number of iterations necessary for convergence, at each incremental step, is presented, considering elastic/elatoplastic and fixed/adaptive analyses. As one can observe, the number of iterations does not drastically increase by considering all iterative approaches in a unique iterative loop. In fact, here, this number is controlled by the most demanding iterative process being carried out, not significantly mattering if other iterative approaches are carried out together. Thus, the single iterative loop considered here provides a quite efficient methodology. In the present example, for instance, the iterative demand of an adaptive nonlinear BEM–FEM coupled analysis would be basically the same of a simple nonlinear FEM

analysis, since the nonlinear analysis is the most demanding iterative process here (one should keep in mind that a Newton– Raphson initial stress approach is being considered). For the present nonlinear model, the CPU time obtained for the single iterative adaptive BEM–FEM coupled analysis is approximately 58% of the CPU time obtained for the single iterative fixed BEM–FEM coupled analysis. In addition, the CPU time obtained for the single iterative fixed BEM–FEM coupled analysis is approximately 34% of the CPU time obtained for the multiple iterative fixed BEM–FEM coupled analysis. As one can observe, the proposed technique considerably enhances the performance of the BEM– FEM coupled analysis.

5.2. Half-space In this application, a half-space is analyzed, in which the region closer to the loaded area is discretized by the FEM, whereas the remaining domain is discretized by the BEM. A sketch of the model is depicted in Fig. 7. The geometry of the problem is defined by L ¼ 6:4 m and H ¼ 4:0 m. The physical properties of the model are E ¼ 109 N=m2 and m ¼ 0:3. A perfectly plastic material obeying the Drucker–Prager yield criterion is assumed, where

Fig. 5. Computed relaxation parameters along each incremental step for the cavity model (results related to higher incremental steps are depicted as lighter gray curves): (a and b) elastic analysis; (c and d) elastoplastic analysis; (a) and (c) fixed FEM discretization; (b) and (d) adaptive FEM discretization.

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

141

Fig. 6. Evolution of the adaptive elastoplastic analysis for the cavity model: FEM discretization and equivalent plastic strains at the end of each incremental step.

Table 1 Number of iterations per incremental step for the cavity model. Incremental step 1

2

3

4

5

Elastic

Fixed Adaptive

8 11

8 11

8 12

8 12

8 12

6 8 13

Elastoplastic

Fixed Adaptive

8 11

8 11

8 12

8 12

48 50

86 115

c ¼ 1:7  102 N=m2 and / ¼ 10 . Two different BEM–FEM discretizations are considered here. In the first discretization, a fixed FEM mesh is employed, which is composed of 5120 elements. In the second discretization, an adaptive FEM mesh is considered, being its initial configuration depicted in Fig. 7 (64 elements). For the first BEM–FEM discretization, 20 boundary elements of length 0.20 m are employed to discretize each vertical common interface, whereas 16 boundary elements of length 0.25 m are considered for the second discretization. For both discretizations, boundary elements of length 0.25 m and 0.20 m are applied to discretize

H FEM BEM

L Fig. 7. Sketch of the half-space model and initial FEM discretization for the adaptive analysis.

142

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

Fig. 8. Computed displacements (modulus) along the FEM deformed meshes (scale factor of 105) for the half-space model: (a and b) elastic analysis; (c and d) elastoplastic analysis; (a) and (c) fixed FEM discretization; (b) and (d) adaptive FEM discretization.

-1.75x10-6

Vertical displacement

-2.00x10-6 -2.25x10-6 -2.50x10-6 -2.75x10-6 -3.00x10-6 -3.25x10-6 Adaptive (single iterative) Fixed (single iterative) Fixed (multiple iterative)

-3.50x10-6 -3.75x10-6 -4.00x10-6 0.0

0.5

1.0

1.5

2.0

2.5

3.0

3.5

4.0

Vertical coordinate Fig. 9. Vertical displacements along the symmetrical vertical axis of the half-space model considering BEM–FEM elastoplastic analyses.

the half-space (which is sufficiently extended up to 25 m at each side of the model) and the horizontal common interface, respectively. For the analyses in focus, 12 incremental steps are considered. In Fig. 8, displacement results are depicted throughout the FEM meshes, considering both discretization procedures and linear and nonlinear analyses. In Fig. 9, the vertical displacements computed along the symmetrical vertical axis of the half-space model are depicted, considering the elastoplastic analyses. In this case, results obtained from a multiple iterative analysis are also depicted in the figure, for comparison. As one can observe, good agreement is observed among the results. In Fig. 10, the curves of the relaxation parameters computed along each incremental step of the analyses are depicted, once again illustrating their complex patterns. In Fig. 11, the evolution of the equivalent plastic strains, at the last three incremental steps, is described, considering the fixed and the adaptive BEM–FEM discretizations. Once again, good agreement is observed among the results. The evolution of the FEM discretization, taking into account the adaptive analysis, is

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

143

Fig. 10. Computed relaxation parameters along each incremental step for the half-space model (results related to higher incremental steps are depicted as lighter gray curves): (a and b) elastic analysis; (c and d) elastoplastic analysis; (a) and (c) fixed FEM discretization; (b) and (d) adaptive FEM discretization.

also illustrated in Fig. 11, considering these last three incremental steps. As one can observe, refinement is properly introduced into the analysis, being the region under the applied load progressively refined. In Fig. 12, computed displacements and tractions along the BEM discretization, at the last incremental step of the elastoplastic analysis, are illustrated, considering both discretization procedures. Since the so-called configuration 2 is considered here, the tractions depicted in Fig. 12 are obtained from the projections of the FEM computed stresses, whereas the displacements depicted in the figure are computed by the solution of the BEM system of equations. As one can observe, although quite a different number of finite elements is considered along the coupling interface for each discretization, providing different traction distributions as depicted in Fig. 12, the patterns of the displacements computed by the BEM are still very similar, illustrating the robustness of the technique. In Table 2, the number of iterations per incremental step for the half-space model is presented. As one can observe, once again the number of iterations does not significantly increase by considering several iterative procedures at a unique iterative loop, illustrating the good performance of the adopted technique. For the present nonlinear model, the CPU time obtained for the single iterative adaptive BEM–FEM coupled analysis is

approximately 19% of the CPU time obtained for the single iterative fixed BEM–FEM coupled analysis. In addition, the CPU time obtained for the single iterative fixed BEM–FEM coupled analysis is approximately 14% of the CPU time obtained for the multiple iterative fixed BEM–FEM coupled analysis. Thus, once again, the considerably superior performance of the proposed technique is observed.

5.3. Cantilever beam In this last application, a cantilever beam is analyzed, in which the first half of the model is discretized by the FEM, whereas its second half is discretized by the BEM. A sketch of the model is depicted in Fig. 13. The geometry of the problem is defined by L ¼ 1:0 m and H ¼ 0:5 m. The physical properties of the model are E ¼ 2  1011 N=m2 and m ¼ 0:3. A perfectly plastic material obeying the von Mises yield criterion is assumed, where the uniaxial yield stress is 1:5  104 N=m2 . For the adaptive BEM–FEM analysis, an initial FEM mesh with 400 elements is considered (as depicted in Fig. 13), as well as 60 boundary elements of equal length (i.e., 0.05 m) are employed (double nodes are also

144

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

Fig. 11. Evolution of the equivalent plastic strains at the last three incremental steps: (a) fixed FEM discretization; and (b) adaptive FEM discretization.

145

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

Fig. 12. Displacements (left) and tractions (right) along the BEM meshes at the last incremental step of the elastoplastic analysis: (a) fixed FEM discretization; and (b) adaptive FEM discretization.

Table 2 Number of iterations per incremental step for the half-space model. Incremental step 1

2

3

4

5

6

7

8

9

10

11

12

Elastic

Fixed Adaptive

18 23

17 27

20 26

18 20

18 21

18 20

18 21

18 20

18 25

21 16

18 21

18 21

Elastoplastic

Fixed Adaptive

18 23

17 27

20 26

18 20

18 21

18 20

18 21

18 20

18 25

23 23

40 37

65 65

H

FEM

BEM

L

L

Fig. 13. Sketch of the cantilever beam model and initial FEM discretization for the BEM–FEM adaptive analysis.

considered at the corners of the BEM subdomain). Results provided by this BEM–FEM coupled configuration are compared to those provided by an adaptive FEM solution, in which an initial mesh of 800 elements is considered to discretize the entire model, and by a fixed FEM solution, in which a mesh of 5200 elements is

employed. For all configurations, 6 incremental steps are considered. In Fig. 14, the vertical displacements computed along the anti-symmetrical horizontal axis of the beam are depicted, considering elastoplastic analyses. As one can observe, good agreement is

146

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

3.0x10 -8

Vertical displacement

0.0 -3.0x10-8 -6.0x10-8 -9.0x10-8 -1.2x10-7 -1.5x10-7

Adaptive BEM-FEM Adaptive FEM Fixed FEM

-1.8x10-7 -2.1x10-7 -2.4x10-7 0.0

0.2

0.4

0.6

0.8

1.0

Horizontal coordinate Fig. 14. Vertical displacements along the anti-symmetrical horizontal axis of the cantilever beam model considering elastoplastic analyses.

observed among the results (it is important to highlight that the considered fixed FEM mesh provides a poorer discretization in the regions were the nonlinear behavior occurs; thus, results provided by this approach are expected to be less accurate). In Fig. 15,

the equivalent plastic strains at the last incremental step of the analyses are plotted, considering the adaptive FEM and BEM– FEM discretizations. Again, good agreement among the results is observed. In Fig. 16, computed displacements and tractions along the BEM discretization, at the last incremental step of the elastoplastic analysis, are illustrated. Since the so-called configuration 1 is considered here, at the common interface, the displacements depicted in Fig. 16 are obtained from the FEM computed displacements, whereas the tractions depicted in the figure are computed by the solution of the BEM system of equations. In Fig. 17, the curves of the relaxation parameters computed along each incremental step of the analyses are depicted and, in Table 3, the number of iterations per incremental step for the beam model is presented. As one can observe, in the present application, convergence is very quickly obtained for the elastic model. Moreover, as it is described in the table, there is basically no additional iterative cost introduced into the analysis by the adopted BEM–FEM coupling approach (i.e., basically the same number of iterations are required by the FEM adaptive analysis and by the BEM–FEM adaptive analysis), once again highlighting the good performance of the adopted single iterative loop. For the present nonlinear model, the CPU time obtained for the adaptive BEM–FEM coupled analysis is approximately 55% of the CPU time obtained for the adaptive FEM analysis.

Fig. 15. Equivalent plastic strains along the FEM deformed meshes (scale factor of 105) for the cantilever beam model: (a) adaptive FEM; and (b) adaptive BEM–FEM.

Fig. 16. Displacements (left) and tractions (right) along the BEM mesh at the last incremental step of the adaptive elastoplastic analysis.

147

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

zoom

(a)

(b)

Fig. 17. Computed relaxation parameters along each incremental step for the cantilever beam model (results related to higher incremental steps are depicted as lighter gray curves): (a) elastic analysis; and (b) elastoplastic analysis.

Table 3 Number of iterations per incremental step for the cantilever beam model. Incremental step 1

2

3

Elastic

FEM BEM–FEM

4 4

4 4

4 5

4 4 4

5 4 5

6 4 5

Elastoplastic

FEM BEM–FEM

4 4

4 4

5 5

22 4

86 24

271 247

6. Conclusions In this work, two iterative BEM–FEM coupling approaches are discussed, according to the prescribed boundary conditions at the BEM–FEM common interface. In both coupling procedures, adaptive nonlinear analyses are enabled within the FEM subdomains, and independent discretizations at common interfaces are allowed. Thus, a quite powerful and versatile numerical methodology is provided, enabling the proper solution of a wide range of complex applications. Moreover, the proposed techniques are quite efficient. Here, optimal relaxation parameters are employed, speeding up the convergence of the iterative coupling, and a single iterative loop is considered, enabling all iterative solutions to be carried out at once. Thus, the computational cost of the analysis is greatly reduced, and the proposed iterative approach becomes highly competitive. As a matter of fact, as it is illustrated in the previous section, the number of iterative steps related to the adopted unified iterative solution is not significantly higher than that of the dominant iterative approach acting isolated. Thus, standard multiple iterative solutions can be avoided, eliminating overly excessive computational costs. Several numerical examples are discussed here, illustrating the versatility and effectiveness of the proposed techniques. In fact, the present work enables a wide range of analyses, providing effective numerical tools to properly deal with several complex problems. The main advantages of the proposed formulations can be summarized as follows: (i) different subdomains can be analyzed separately, leading to smaller and better-conditioned systems of equations (different solvers, suitable for each subdomain, may be employed);

(ii) coupled systems may be easily treated by separate program modules, taking full advantage of specialized features and disciplinary expertise; (iii) matching nodes at common interfaces are not required, greatly improving the flexibility and versatility of the coupled analyses, especially when adaptive discretizations are considered; (iv) nonlinear and/or adaptive analyses (as well as other iterative-based analyses) may be carried out in the same iterative loop of the iterative coupling, not introducing a relevant extra computational effort for the model; (v) more efficient analyses can be obtained, once the global model can be reduced to several subdomains with reduced size matrices, and adaptive techniques can be applied, enabling a further reduction of the model by optimal discretization; (vi) the advantages of the BEM and of the FEM can be fully explored, evading their drawbacks and characterizing an efficient union of different ‘‘worlds’’.

Acknowledgements The financial support by CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) and FAPEMIG (Fundação de Amparo à Pesquisa do Estado de Minas Gerais) is greatly acknowledged. References [1] Katsikadelis JT. Boundary elements: theory and applications. Elsevier; 2002. [2] Zienkiewicz OC, Kelly DM, Bettes P. The coupling of the finite element method and boundary solution procedures. Int J Numer Meth Eng 1977;11:355–76. [3] Pavlatos GD, Beskos DE. Dynamic elastoplastic analysis by BEM/FEM. Eng Anal Boundary Elem 1994;14:51–63. [4] von Estorff O, Firuziaan M. Coupled BEM/FEM approach for nonlinear soil/ structure interaction. Eng Anal Boundary Elem 2000;24:715–25. [5] Lin CC, Lawton EC, Caliendo JA, Anderson LR. An iterative finite element – boundary element algorithm. Comput Struct 1996;39:899–909. [6] Elleithy WM, Al-Gahtani HJ, El-Gebeily M. Iterative coupling of BE and FE methods in elastostatics. Eng Anal Boundary Elem 2001;25:685–95. [7] Elleithy WM, Tanaka M, Guzik A. Interface relaxation FEM–BEM coupling method for elasto-plastic analysis. Eng Anal Boundary Elem 2004;28:849–57. [8] Elleithy WM, Grzhibovskis R. An adaptive domain decomposition coupled finite element–boundary element method for solving problems in elastoplasticity. Int J Numer Meth Eng 2009;79:1019–104.

148

D. Soares, L. Godinho / Computers and Structures 156 (2015) 134–148

[9] Elleithy WM. Multi-region adaptive finite element-boundary element method for elasto-plastic analysis. Int J Comput Math 2012;89:1525–39. [10] Jahromi HZ, Izzuddin BA, Zdravkovic L. A domain decomposition approach for coupled modelling of nonlinear soil-structure interaction. Comput Methods Appl Mech Eng 2009;198:2738–49. [11] Boumaiza D, Aour B. On the efficiency of the iterative coupling FEM–BEM for solving the elasto-plastic problems. Eng Struct 2014;72:12–25. [12] Soares D, von Estorff O, Mansur WJ. Iterative coupling of BEM and FEM for nonlinear dynamic analyses. Comput Mech 2004;34:67–73. [13] Soares D. An optimised FEM–BEM time-domain iterative coupling algorithm for dynamic analyses. Comput Struct 2008;86:1839–44. [14] Soares D. Fluid-structure interaction analysis by optimised boundary element – finite element coupling procedures. J Sound Vib 2009;322:184–95. [15] Soares D. FEM–BEM iterative coupling procedures to analyze interacting wave propagation models: fluid–fluid, solid–solid and fluid–solid analyses. Coupled Syst Mech 2012;1:19–37. [16] Bendali A, Boubendir Y, Fares M. A FETI-like domain decomposition method for coupling finite elements and boundary elements in large-size problems of acoustic scattering. Comput Struct 2007;85:526–35. [17] Soares D, Godinho L. An optimized BEM–FEM iterative coupling algorithm for acoustic–elastodynamic interaction analyses in the frequency domain. Comput Struct 2012;106–107:68–80.

[18] Coulier P, François S, Lombaert G, Degrande G. Coupled finite element – hierarchical boundary element methods for dynamic soil–structure interaction in the frequency domain. Int J Numer Meth Eng 2014;97:505–30. [19] Godinho L, Soares D. Frequency domain analysis of interacting acoustic– elastodynamic models taking into account optimized iterative coupling of different numerical methods. Eng Anal Boundary Elem 2013;37:1074–88. [20] Soares D, Godinho L. An overview of recent advances in the iterative analysis of coupled models for wave propagation. J Appl Math 2014; Article ID 426283. [21] Chen WF, Han DJ. Plasticity for structural engineers. New York: Spring-Verlag; 1988. [22] Khan AS, Huang S. Continuum theory of plasticity. New York: John Wiley & Sons; 1995. [23] Brebbia CA, Telles JCF, Wrobel LC. Boundary element techniques. Berlin: Springer-Verlag; 1984. [24] Brebbia CA, Dominguez J. Boundary elements, an introductory course. Southampton: WIT Press; 1998. [25] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [26] Belytschko T, Liu WK, Moran B. Nonlinear finite elements for continua & structures. New York: J. Wiley & Sons; 2000. [27] Chen L, Zhang CS. AFEM@matlab: a Matlab package of adaptive finite element methods. University of Maryland at College Park; 2006.