Stabilized mixed triangular elements with area bubble functions at small and large deformations

Computers and Structures 138 (2014) 172–182

Contents lists available at ScienceDirect

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

Stabilized mixed triangular elements with area bubble functions at small and large deformations Ismail Caylak ⇑, Rolf Mahnken University of Paderborn, Chair of Engineering Mechanics, Warburger Str. 100, D-33098 Paderborn, Germany

a r t i c l e

i n f o

Article history: Received 19 November 2012 Accepted 3 January 2014 Available online 11 February 2014 Keywords: Finite elements Triangular elements Mixed elements Area bubble functions

a b s t r a c t A well known approach for stabilization of mixed finite triangular elements are volume bubble functions. In this paper a new concept of area bubble functions in geometrically linear (elastic and elasto-plastic) and geometrically non-linear elastic regimes is presented. These functions are used in order to enrich the displacement field. In the numerical examples, simulations for Cook’s membrane problem and a block under central pressure demonstrate, that the proposed formulation avoids locking and reduces stress oscillations for incompressible materials. Furthermore, the results reveal that area bubble functions are superior to volume bubble functions. This paper is an extension of the works published by the authors regarding tetrahedral elements. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Triangular elements can be easily meshed and therefore are very attractive for industrial applications. However, low order one-field triangular elements show locking and stress oscillations due to the incompressibility constraint and therefore lead to unstable finite elements. One possibility to prevent locking is a mixed finite element formulation. Standard mixed triangular elements use the same order of spaces for displacement and the pressure, respectively. Since the pressure space is too rich compared to the displacement space, stress oscillations occur which is well known and visible in the numerical examples. According to Rannacher [1] there exist essentially two approaches for stabilization: (1) adding stabilization terms in the continuity equation; and (2) augmenting the displacement field by so-called bubble functions. Methods for approach (1) are least-square stabilization, subgrid scale and orthogonal sub-grid scale are methods which consider additional stabilization terms in the continuity equation, e.g. in fluid mechanics [2] and in non-linear solid mechanics [3,4]. A mixed finite element for small strain elasto-viscoplastic material behavior based on the least-squares method is presented in [5,6]. A difficulty of least-square methods lies in the appropriate choice of the stabilization factor. The starting point of the orthogonal sub-grid scales is the least-square stabilization and was introduced in [7] and followed by the works of [8–10] for stabilized formulations at small and large deformations. The basic idea is ⇑ Corresponding author. Tel.: +49 5251 602285; fax: +49 5251 603483. E-mail addresses: [email protected] [email protected] (R. Mahnken). http://dx.doi.org/10.1016/j.compstruc.2014.01.006 0045-7949/Ó 2014 Elsevier Ltd. All rights reserved.

(I.

Caylak),

rolf.

the split to different scales of space resolution, the resolvable scales and unresolvable scales, or sub-grid scales. The sub-grid scale method circumvents the LBB-condition. It is a locking-free formulation with no pressure oscillations. Augmenting the displacement field by so-called bubble functions in approach (2), arises from the ‘‘MINI’’ element introduced by Arnold et al. in [11] and discussed in [12,13]. The use of bubble functions satisfies the LBB-condition [11], however, a mesh-dependent parameter is required, which tends to instability at large deformations. A shortcoming of additional bubble functions is the presence of small pressure oscillations. Bubble functions are also used in mixed-enhanced elements based on the classical enriched strain method introduced by Simo and Rifai [14]. A mixed-enhanced formulation using volume bubble functions for the enhanced part is developed in [15] for tetrahedral elements at small and large deformations and an optimization has been presented in [16]. Recent works which also consider volume bubble functions are [17–19]. Furthermore, bubble functions could be used for the method of incompatible modes introduced by Wilson et al. [20] and is used in [14,21]. The drawback of the classical formulation of Simo and Rifai [14] is the evaluation of complex non-linear terms which arise due to the enriched part. Apart from stabilization methods of mixed finite elements, which will not be discussed in detail here, further methods like nodal pressures/ strains [22–26,17], or reduced integration [27,28], followed by recent works [29,30] only for quadrilaterals. In this paper we concentrate on the stabilization approach (2), where the displacement field is augmented by area bubble functions and the focus holds on the following aspects:

173

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

 In contrast to [31], where volume bubble functions are used, in this paper new area bubble functions are introduced for stabilization of mixed triangular elements.  In order to avoid the non-linear terms of classical methods, like the method of incompatible modes, constant stabilization terms are considered.  In contrast to the B-bar, F-bar and reduced integration method or other approaches, our mixed formulation is able to solve problems with m ¼ 0:5 and to avoid volume locking for incompressible materials. Additionally, stabilization of the mixed method with constant stabilization terms damps drastically stress oscillatory effects. The structure of this paper is as follows: In the next section new area bubble functions and constant stabilization matrices are presented. Unnecessary repetitions, like the governing equations which do not differ from [18,19,32,31] are not presented here. In the representative examples Cook’s membrane problem and a block under central pressure illustrate the performance of the presented approaches in comparison to existing finite element formulations for linear elastic, elasto-plastic and non-linear elastic materials. 2. Finite element matrix formulations for triangular elements

1: J ¼ const: and j ¼ detðJÞ ¼ const: 2 3 X N i;g yi 6 6 1 6 i¼1 2: J is invertible and J 1 ¼ 3 detðJÞ 6 4 X  Ni;g xi i¼1

3 3 X  Ni;n yi 7 7 i¼1 7: ð3Þ 7 3 X 5 Ni;n xi i¼1

Remarks 1. The standard terms of a triangular element are given in the Pascal triangle of Fig. 5. Volume bubble functions introduce additional terms ng, n2 g and ng2 . 2. The derivatives of the volume bubble functions with respect to the local coordinates n, g are quadratic and obtained from Eq. (2) as

1: N;n ¼ 27ðgu  ngÞ

ð4Þ

2: N;g ¼ 27ðnu  ngÞ: 2.2. Area bubble functions for triangular elements

In addition to the volume bubble shape function in Eq. (2), new area bubble functions are introduced in order to eliminate the locking effects and to reduce stress oscillatory effects as follows: 2

1: N 1 ðnÞ ¼ 4gu þ n2 þ n

2.1. Isoparametric concept for triangular elements h

2: N 2 ðnÞ ¼ 4nu þ

Sne

In Fig. 1a discretization of the domain B  B ¼ e¼1 Be into ne triangular elements with three nodes is considered. Each element occupies a subdomain Be . In the isoparametric domain Be0 , a local coordinate system with coordinates n; g is introduced, with the properties

0 6 n; g 6 1:

ð1Þ

In addition to the standard shape functions Ni ; i ¼ 1; 2; 3 : Be0 # Be of the isoparametric concept, the volume bubble function [11] is expressed as

NðnÞ ¼ 27ngu:

ð2Þ

This function is illustrated in Fig. 2. It gives zero contribution along the three edges and has the value 1 at the center of the triangle, see e.g., [33]. The derivatives are easily converted from the local coordinate system to the global coordinate system by using the Jacobian matrix J. In Fig. 3 the transformation from the local coordinates to the global coordinates is illustrated, which for linear triangular elements has the important properties

3

3: N ðnÞ ¼ 4ng þ

g2 2

u2 2

þg

ð5Þ

þ u:

An illustration of the area bubble function N 3 is given in Fig. 4(b) left. The function value is zero at node 1 and node 2, and one at the center of the considered edge u ¼ 0. At node 3, N3 has a maximum value of 3=2. A quantitative illustration of N 3 along line a and c of Fig. 4(b) left is given in Fig. 4(b) right. Remarks 1. Compared to the volume bubble shape functions in Eq. (2), area bubble functions introduce additional terms n2 , ng and g2 , which are summarized in the Pascal triangle in Fig. 5. 2. The derivatives of the area bubble functions with respect to the local coordinates n, g are obtained from Eq. (5) as

1: N1;n ¼ 4g þ n þ 1; N2;n ¼ 4ðu  nÞ; N3;n ¼ 4g  u  1; 2: N1;g ¼ 4ðu  gÞ; N2;g ¼ 4n þ g þ 1; N3;g ¼ 4n  u  1: ð6Þ In Eq. (6), the derivatives of the area bubble functions are linear. An illustration of the derivative of the area bubble function N3 ; N 3;n is given in Fig. 6. 3. As shown in Eq. (17) of Appendix A the following integrals over the domain Be0 of the reference element in the isoparametric space have the properties

Z

Ni;n dV 0 ¼ 0 and Be0

Z

Ni;g dV 0 ¼ 0

ð7Þ

Be0

R RR and where Be0 dV 0 ¼ dndg. In Fig. 6, therefore, the volume above the triangle D123 is equal to the volume below the triangle. The property in Eq. (7) is exploited for verification of the patch test. For specific explanation please see [18,19,32]. 2.3. Finite element matrix formulation for the mixed method of incompatible modes

Fig. 1. Discretization in triangular elements.

In the mixed method of incompatible modes the trial functions ~ 2 R2 , the pressure for the enriched displacement interpolation u

174

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

(a)

(b) 2

2 N

P 3

3

1

1

Fig. 2. Illustration of volume bubble function: (a) point P with maximum value of the volume bubble function; (b) volume bubble N.

Fig. 3. Isoparametric mapping to the reference configuration.

(a)

η

2 P

(b) N3 d

P

1

3

P

3

ξ

3

N

1

c

3 2

1

2

η

2

3 2

b a

3

ξ

1

c

1

a

3

1

ξ

Fig. 4. Illustration of area bubble functions: (a) points P i ; i ¼ 1; 2; 3 for all three area bubble functions; (b) different views of area bubble function N 3 .

~ ðnÞ ¼ 1: u Standard

ξ ξ ξ

3

2: pðnÞ ¼

nB X Ni ðnÞv^ i ¼ uðnÞ þ v ðnÞ; i¼1

3 X

^i ¼ Np ^; Ni ðnÞp

ð8Þ

i¼1

2

ξη η 2 3 2 ξ η ξη η

area bubble

3:

eðu^; v^ Þ ¼

nB 3 X X ^i þ ^ þ Bv v^ ; Biu ðnÞu Biv ðnÞv^ i ¼ Bu u i¼1

Fig. 5. The Pascal triangle.

field interpolation p 2 R and the strain field interpolation within each element read

i¼1

where uðnÞ and v ðnÞ denote the compatible and incompatible part of the displacement field, respectively. The unknowns in Eq. (8) are: ^ i ¼ ½u ^ xi ; u ^ yi T ; i ¼ 1; 2; 3; p ^i ; i ¼ 1; 2; 3; v ^ i ¼ ½v^ xi ; v^ yi T ; i ¼ 1; . . . ; nB , u where nB ¼ 1 for volume bubble functions and nB ¼ 3 for area bubble functions. Additionally

volume bubble

e 2 R3

 Bu ¼ B1u ^ ¼ ½u ^1 u

B2u ^2 u

N3,˫

ˤ 2

 1 2 3 B3u ; N ¼ ½N ; N ; N ; N ¼ ½N1 ; N2 ; N3 ; ^ 3 T ; v^ ¼ ½ v^ 1 u

v^ 2 v^ 3 T ; p^ ¼ ½ p^1

ˤ 2

3

˫=0 ˤ=0.4

P

P1 3

^i þ Ni ðnÞu

i¼1

η

2

3 X

2

P

1

˫

3

Fig. 6. The derivative of the area bubble function, N 3;n .

˫=0.75 ˤ=0.25

1

˫

^2 p

^ 3 T p

ð9Þ

175

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

with

2

Ni;x 6 Biu ¼ 6 4 0 Ni;y

3

2

Ni;x 7 6 i 6 Ni;y 7 5 and Bv ¼ 4 0 Ni;y N i;x 0

0

3

7 Ni;y 7 5 Ni;x

are given. Inserting the matrix formulations into the weak form, after some rearrangements the finite element residual equations in step 3 of Table 1 are derived. Then in step 4 of Table 1 the non-linear residual equations are solved iteratively with a Newton algorithm. In [32] the corresponding equations are given for large deformation problems. In the following remarks only the general procedure is described. Detailed information is given in [19,32].

4. Regarding the iterative Newton procedure, k^e is kept constant, ^eðkÞ is varying, which is a key idea of this concept. In whereas k T this way quadratic convergence of the Newton algorithm in Table 1 is obtained. 5. In classical methods of augmenting the displacement field by so-called bubble functions one can substitute, e.g., the factor 27 in Eq. (2) by a non-dimensional constant c. For different values of c one would expect slightly different solutions for a given mesh, but converge to exactly the same solution at the same rate. Multiplying Eq. (2) by c modifies the element matrices in Table 1 as 

T

keuv ¼ ckeuv ¼ kev u ;



kevv ¼ c2 kevv ;



T

kev p ¼ ckev p ¼ kepv :

Then, the constant stabilization matrix k^e in Table 1 becomes

Remarks 1. As explained in e.g. [19], for the method of incompatible modes the system consist of three finite element residuals Ru ; Rev and Rp for the standard, incompatible and pressure part. Since the incompatible residual part Rev is valid on the element level, ^ are eliminated by static conthe incompatible displacements v densation. Therefore, the system reduces to the residuals given in Table 1 with the additional term k^e . 2. The matrix k^e in Table 1 is kept constant during the iterative Newton procedure for non-linear problems and could be interpreted as a stabilization of the standard two-field formulation. Furthermore, the constant elasticity matrix C is used in order to determine k^e . ^eðkÞ , namely 3. The reminder terms of the element matrices in k T keuu ; keup ; kepu and kepp , are equivalent to the standard two-field method. Note that for determination of the tangent matrix ^eðkÞ , the consistent tangent matrix Ct is used. k T

" # keuv 1  e 1  e  ^ e k ¼c e k c kv u kpv c2 vv

 kev p ¼ k^e :

This result means, that parameter c in our formulation has no influence, which also has been verified numerically. As explained in Section 2.2 the integral of the derivative of the shape functions must be zero in order to verify the patch test. Consequently, the parameters 1; 12 and 4 of the area bubble functions in Eq. (4) can not been chosen arbitrarily. Alternatively, a parameter c could be introduced in Table 1 by modifying 

kevv ¼ ckevv :

ð10Þ

Then, different parameters for c yield different approximate solutions for a given mesh, but converge to the same solution. This effect is illustrated in the numerical example in Section 3.1. 3. Numerical examples In this section, results for different finite element examples obtained with triangular are shown for linear elastic, elasto-plastic and non-linear elastic problems. The results are also compared to those of different elements familiar from the literature. All finite element formulations are distinguished with the following nomenclature:

Table 1 Newton algorithm for solution of non-linear problems. ^ ðkÞ , p ^ðkÞ , set k ¼ 0 1. Initialize for all elements: u 2. Strains, stresses, consistent tangent, pressure ðkÞ

eðkÞ ¼ Bu u^ðkÞ ; r^ ðkÞ ¼ r^ ðeðkÞ Þ; CðkÞ T ¼

^ dr ; deðkÞ

^ðkÞ pðkÞ ¼ N p

3. Residual including stabilization matrix

"

ðkÞ Ru ðkÞ Rp

#

3 " ðkÞ #   ^ ðkÞ þ pðkÞ 1ÞdV BT ðI r ne ne F ext ^ u Be u dev   5  A k^e ¼ A 4R  A ðkÞ 1 T 1 T ^ ðkÞ e¼1 e¼1 e¼1 0 ^ p N 3 1 r  pðkÞ dV Be K ne

2 R

"

# Z keuv  e 1  e  T kv u kev p ; keuv ¼ kvv BTu D Bv dV ¼ kev u e kpv Be Z T BTv D Bv dV; kev p ¼ BTv 1 N dV ¼ kepv ; D ¼ Idev C Idev :

where k^e ¼ kevv ¼

Z Be

Be

4. Solve finite element matrix equation ne

A

e¼1



^eðkÞ  k^e k T

keup ¼

Z Be

"

ðkÞ

^ Du ðkÞ Dp^

#

ne

¼ A

e¼1

"

ðkÞ

Ru

#

ðkÞ

Rp

# Z keup ðkÞ ðkÞ ðkÞ ðkÞ BTu DT Bu dV; DT ¼ Idev CT Idev ; keuu ¼ e kpp Be Z 1 T T BTu 1N dV ¼ kepu ; kepp ¼ N N dV: Be K

eðkÞ

^ where k T

"

¼

keuu kepu

ðkÞ

5. Update for all elements

^ ðkþ1Þ ¼ u ^ ðkÞ  Du ^ ðkÞ ! k ¼ k þ 1; u

GOTO 2

Linear displacements Quadratic displacements

T1 T2

Quadratic displacements; linear pressure Linear displacements; reduced integration Linear displacements; linear pressure Linear displacements; linear pressure; volume bubble for IM Linear displacements; linear pressure; area bubble for IM

T2P1 Q 1R T1P1 T1P1IM2ST T1P1IM6ST

Here T1 and T2 are standard triangular with linear and quadratic interpolations for the displacements, respectively. T1P1 includes linear interpolations for the displacements and linear interpolations for the pressure. Note that this element does not satisfy the inf–sup condition. T2P1 includes quadratic interpolations for the displacements and linear interpolations for the pressure. Q 1R is a standard quadrilateral element with linear interpolation for the displacement and reduced integration. T1, T2, T2P1 and Q 1R are available in the commercial finite element program Abaqus[34]. The mixed finite element formulations for triangular elements using volume and area bubble functions are T1P1IM2ST and T1P1IM6ST and have been implemented in a UEL user subroutine of the finite element program Abaqus [34]. Furthermore the abbreviation IM means incompatible mode and ST means stabilization. For linear and non-linear elasticity the very popular Cook’s membrane problem is chosen. A compressible and nearly incompressible material, with m ¼ 0:33 and m ¼ 0:4999, is considered

176

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

48

16

C

t

A

44

t=1 [mm] 2

(a)

3

(b)

1

B Fig. 7. Cook’s membrane problem (linear elastic): Geometry and discretization with 6 elements along the edge A–B.

(a)

(b)

3.7

Displacement at C [mm]

Displacement at C [mm]

3.6 3.5 3.4

T1 T1P1

3.3

T2 T2P1

3.2

3.7

3.6

3.5

Q1R

3.4

T1P1IM2ST

3.3

Q1R

T1P1IM6ST

Q2P1

3.1

3.2

0

22

44

66

88

0

22

No. of elements along A-B

44

66

88

No. of elements along A-B

Fig. 8. Cook’s membrane problem (linear elastic): Convergence study for the displacement at point C for T1, T1P1, T2, T2P1, Q 1R, Q 2P1, T1P1IM2ST and T1P1IM6ST with m ¼ 0:33.

(a)

(b)

3.3

3.15 3.10

2.9

Displacement at C [mm]

Displacement at C [mm]

3.1

2.7 2.5 2.3

T1

2.1

T1P1 T2

1.9

T2P1

0

22

44

66

2.95

Q1R

2.90

T1P1IM2ST T1P1IM6ST

Q2P1

1.5

3.00

2.85

Q1R

1.7

3.05

2.80 88

No. of elements along A-B

0

22

44

66

88

No. of elements along A-B

Fig. 9. Cook’s membrane problem (linear elastic): Convergence study for the displacement at point C for T1, T1P1, T2, T2P1, Q 1R, Q 2P1, T1P1IM2ST and T1P1IM6ST with m ¼ 0:4999.

for the linear elastic case. Firstly, a comparison between our element and linear, quadratic and mixed elements of the finite element program Abaqus [34] is presented. In addition Cook’s membrane is simulated with same conditions as presented in

Chiumenti et al. [35] to compare our element with the stabilized formulation in [35]. For non-linear elasticity Cook’s membrane problem is investigated with the same conditions as presented in de Souza Neto et al. [25], where the F-bar quadrilateral formulation

177

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

Fig. 10. Cook’s membrane problem (linear elastic): Convergence study for the displacement at point C for T1P1IM6ST with m ¼ 0:4999 considering different stabilization parameters.

from [25] is taken as the reference solution. We also give results for the case of exact incompressibility, with m ¼ 0:5, where different elements from literature fail to give solutions. The last example, block under central pressure, is used in the literature to test the performance of elements in highly compressed regions where extreme mesh distortion occurs, see e.g. the work of Reese [36] for an elasto-plastic material. In all numerical examples our elements produce good results for convergence and damp the stress oscillations. 3.1. Cook’s membrane problem -linear elasticIn the first example, Cook’s membrane problem is investigated. This example is one of the most popular benchmarks for investigating the finite element performance in case of locking effects. The geometry is shown in Fig. 7(a). The panel is fully constrained on the left side, and is loaded with a surface traction t ¼ 10 MPa on the right side. We consider plane strain conditions for the compressible behavior with Poisson’s ratio m ¼ 0:33 and for the incompressible behavior with m ¼ 0:4999, with Young’s modulus E ¼ 1000 MPa. For the convergence study, different discretizations with 6, 11, 22, 44, and 88 elements along the edge A–B of Fig. 7(a) are examined. The finite element discretization with six elements along the edge A–B is illustrated in Fig. 7(b). Subsequently, we compare in Fig. 8 and 9(a) the convergence behavior of literature elements with m ¼ 0:33 and m ¼ 0:4999, where Q 1R is regarded as the reference solution due to the fact that it has the fastest convergence. As expected, the linear interpolated element T1 renders a large deviation from the quadrilateral

(a)

100

element Q 1R. T1P1 also reveals inadequate performance for coarser meshes. In Figs. 8 and 9(b) the proposed area bubble element T1P1IM6ST is compared with volume bubble element T1P1IM2ST and the reference element Q 1R. All of the elements reveal reasonable results. Furthermore, Figs. 8 and 9(b) shows that our new area bubble element T1P1IM6ST converges faster than the volume bubble element T1P1IM2ST. Regarding area bubble functions, the influence of the stabilization parameter c of Eq. (10) for three different values is given in Fig. 10. As expected, different parameters for c yield different approximate solutions for a given mesh, but converge to the same solution. Additionally, the stress distribution S11 along the clamped edge A-B for the discretization with 22 elements along A–B is considered. In Fig. 11(a) the stress distributions are compared for element T1, T1P1, Q 1R. It is obvious that the curve T1 shows strong oscillations, with large deviation in the stress value at point A of Fig. 11(a) compared with the reference solution Q 1R. T1P1 shows evident oscillations as well. Moreover, in Fig. 11(b), the results of T1P1IM2ST and T1P1IM6ST are compared against the reference solution Q 1R. For both, volume bubbles and new area bubbles, we observe no oscillations in stress distribution. Fig. 12 shows the contour plots for the stress S11 . The contour plot of T1 shows strong oscillations and non-smooth stress distributions, while the reference solution Q 1R exhibits excellent results. T1P1IM2ST and T1P1IM6ST reveal nearly the same stress distributions as Q 1R. In addition Cook’s membrane for linear elasticity is simulated with same conditions as presented in Chiumenti et al. [35] to compare our element with the stabilized formulation in [35]. For nearly incompressible limit (m ¼ 0:4999) Chiumenti et al. [35] demonstrate a faster convergence of their element formulation to the exact solution compared with standard elements Q 1P0, Q 1 and T1. Therefore, in Fig. 13 a comparison between the element formulation of Chiumenti et al. [35] and the stabilized mixed formulation with area bubble functions proposed in this work is presented. The abbreviation T1P1C is used for Chiumenti et al. [35] elements. Fig. 13 shows no significant differences between T1P1IM6ST and T1P1C with good results also for coarser discretization. Several methods like the B-bar, F-bar and reduced integration from the literature fail to give solutions for exact incompressibility with m ¼ 0:5. However, our mixed formulation of Table 1 is not effected by this deficiency which could be observed in Fig. 13. The extension ‘‘Inc.’’ in Fig. 13 means incompressible. As expected, nearly the same results are obtained.

(b)

80

0 0

40 20 0 -20

0

11

22

33

-40 -60 -80 -100

T1 T1P1

44

Stress-S11 [Mpa]

Stress-S11 [MPa]

60

10

11

22

33

44

-10 -20 -30 Q1R -40 T1P1IM2ST -50 T1P1IM6ST

Q1R -60

Fig. 11. Cook’s membrane problem (linear elastic): Stress S11 along the clamped edge side A-B for T1, T1P1, T2, T2P1, Q 1R, Q 2P1, T1P1IM2ST and T1P1IM6ST with m ¼ 0:4999.

178

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

(a)

(b)

2

2 1

3

3

(c)

(d)

2

2

3

1

1

3

1

Fig. 12. Cook’s membrane problem (linear elastic): Stress field S11 for T1, Q 1R, T1P1IM2ST, T1P1IM6ST with

m ¼ 0:4999.

as the reference solution. In Fig. 14 the convergence behavior of different element formulations is considered, where the linear interpolated element T1 renders a large deviation from the reference quadrilateral element F  bar. T1P1 also reveals inadequate performance for coarser meshes. The proposed new area bubble elements show the best results. Fig. 15 shows the stress contour plots S11 . The contour plot of T1P1 shows oscillations in the stress in 1-direction, while the new area bubble functions T1P1IM6ST exhibit excellent results. In Table 2 the development of the largest residual force during the Newton Raphson equilibrium iteration is presented for the new area bubble elements T1P1IM6ST. 3.3. Block under central pressure -elasto-plastic-

Fig. 13. Cook’s membrane problem (linear elastic): Vertical displacement at point C versus number of elements per side for area bubble elements compared with [35].

3.2. Cook’s membrane problem -non-linear elasticFor non-linear elasticity Cook’s membrane problem is investigated with the same conditions as presented in de Souza Neto et al. [25]. The panel is loaded with a surface traction t ¼ 100 MPa on the right side. An incompressible Neo-Hookean material behavior with compression modulus j ¼ 400942 MPa and shear modulus l ¼ 80:1938 MPa is considered. As in the work of de Souza Neto et al. [25], the F  bar quadrilateral elements are regarded

The last example, block under central pressure, is used in the literature to test the performance of elements in highly compressed regions where extreme mesh distortion occurs, see e.g., the work of Reese [36] for an elasto-plastic material, and is therefore ideally suited to show the difference between locking and non-locking finite elements in J 2 -plasticity. Fig. 16(a) illustrates the geometry of the block. For the convergence study, discretizations with 5, 10, 20, 40 and 80 elements along the edge A–C are examined, see Fig. 16(b). The load is limited to t ¼ 0:025 MPa and the block is fully constrained on the bottom. In order to simulate a nonlinear material behavior of von Mises plasticity with isotropic hardening, the material parameters are chosen as follows: E ¼ 10 MPa, m ¼ 0:333, Y 0 ¼ 0:02 MPa and ET ¼ 0:1 MPa. ET is the elasto-plastic stiffness.

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

179

Fig. 14. Cook’s membrane problem (non-linear elastic): Convergence study for the displacement at point C for T1; T1P1; F  bar, T1P1IM2ST and T1P1IM6ST.

Fig. 15. Cook’s membrane problem (non-linear elastic): Stress field S11 for T1P1 and T1P1IM6ST.

Table 2 Cook’s membrane problem (non-linear elastic): Development of the largest residual force during equilibrium iteration for the area bubble elements T1P1IM6ST. Iteration

1

2

3

4

5

10.7

1.1

5.9E  3

4.4E  07

1.6E  11

In Fig. 17(a), the equivalent plastic strain along section A  C is given. All quadratic triangular elements and the reduced quadrilateral element Q 1R show no significant differences. Due to the incompressible material behavior caused by the plastic

10

(a)

(b)

t

2 3

t=1

deformation, the approach of T1 elements leads to volume locking. The stabilized mixed triangular elements T1P1IM2ST and T1P1IM6ST in Fig. 17(b) show no significant differences compared with the reference solution Q 1R. In Fig. 18 the stress distribution S22 along the clamped edge C–D for various elements is compared. The element T1 shows large deviations from the reference solution Q 1R, which leads to the conclusion that T1 is effected by locking. All other elements from literature, except Q 1R elements, show stress oscillations close to point where the stress gradient changes significantly increases. A better agreement could be observed for volume bubble elements T1P1IM2ST and the new area bubble elements T1P1IM6ST, which

A

B [mm]

1

C

D

20 Fig. 16. Block under central pressure (elasto-plastic): Geometry and discretization with five elements along the edge A–C.

180

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

(a)

(b)

Fig. 17. Block under central pressure (elasto-plastic): Equivalent plastic strain epl v along A  C for T1, T2, T2P1, Q 1R, Q 2P1, T1P1IM2ST and T1P1IM6ST.

(a)

0.03

(b)

T1 T2

0.02

0.03 Q1R 0.02

T1P1IM2ST

Q1R Q2P1

0.00 0

Stress-S22 [MPa]

4

6

8

10

T1P1IM6ST 0.00 0

-0.03

-0.03 -0.04

Section along C-D [mm]

-0.02 3

4

5

6

-0.03 T1

(d)

6

8

10

Section along C-D [mm]

2

3

4

5

6

-0.03 Q1R T1P1IM2ST

T2P1

T1P1ES2ST

Q1R

T1P1IM6ST

Q2P1 Section along C-D [mm]

4

-0.02

T2

-0.04

2

-0.01 -0.02

2

T1P1ES2ST

0.01

-0.02

-0.04

(c)

2

-0.01

Stress-S22 [MPa]

0.01

Stress-S22 [MPa]

Stress-S22 [MPa]

T2P1

-0.04

Section along C-D [mm]

Fig. 18. Block under central pressure (elasto-plastic): Stress S22 along C–D for T1, T2, T2P1, Q 1R, Q 2P1, T1P1IM2ST and T1P1IM6ST.

is illustrated in Fig. 18(b). The zoom in Fig. 18(d) shows low oscillations for volume bubble functions T1P1IM2ST, whilst no oscillations occur for the proposed new area bubble functions T1P1IM6ST.

4. Conclusions In this paper a stabilization of triangular elements using area bubble functions is presented. In addition to volume bubble functions, new area bubble functions are used. This is gained by using

constant stabilization matrices. Numerical examples like Cook’s membrane and block under central pressure are investigated. This stabilization scheme is compared to different element formulations known from the literature and the poor performance of low order triangular elements with respect to locking and stress oscillations is improved. Furthermore, the stabilized mixed finite element formulation with area bubble functions show better results than stabilized mixed finite elements with volume bubble functions. Triangular elements facilitate more convenient manipulations in adaptive mesh-refinement of h-type and should be of low order

181

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

to reduce the computation time. Future developments will be also directed to adaptive mesh refinement of h-type on the basis of error estimators, which plays a prominent role in finite element programs.

The authors are grateful to some helpful comments of the unknown reviewers. Appendix A. Integration of the reference element in the isoparametric space For the 2D formulation, the proof is detailed for all integrals as follows: Z 1 Z 1n Z 1 Z 1n Z 1 1n N1;n dgdn ¼ ð4g þ n þ 1Þdgdn ¼ ½2g2 þ ng þ g0 dn 0 0 0 0 0 Z 1 ð2ð1  nÞ2 þ nð1  nÞ þ ð1  nÞÞdn ð11Þ ¼ 0 Z 1 1 ¼ ð1 þ 4n  3n2 Þdn ¼ ½n þ 2n2  n3 0 ¼ 1 þ 2  1 ¼ 0 0 Z 1 Z 1n Z 1 Z 1n N2;n dgdn ¼ 4ðu  nÞdgdn 0 0 0 0 Z 1 Z 1n ¼4 ð1  2n  gÞdgdn 0

¼4

Z

1

0

g  2ng 

0

"

#1

g2

1n

2

dn ¼ 4

Z

1

0

0

 1 3  2n þ n2 dn 2 2

ð12Þ

 n n3 1 1  n2 þ ¼4 1þ ¼0 2 2 2 2 0 Z 1 Z 1n Z 1 Z 1n N3;n dgdn ¼ ð4g  u  1Þdgdn ¼4

0

¼

0 1 Z 1n

Z 0

0

ð2 þ n þ 5gÞdgdn

0

1n  Z 1

5 2 1 3 g dn ¼  2n þ n2 dn 2 2 2 0 0 0 " #1

 3 n n 1 1 ¼  n2 þ ¼0 ¼ 1þ 2 2 2 2 0 Z 1 Z 1n Z 1 Z 1n N1;g dgdn ¼ 4ðu  gÞdgdn ¼

Z

0

1

0

0

¼4

Z

1

2g þ ng þ

0

¼4

Z

¼4 Z

Z

1

0

0

ð1  n  2gÞdgdn

0 1n

½g  ng  g2 0 dn

ð14Þ

1

ð1  n  n þ n2  1 þ 2n  n2 Þdn ¼ 4

0 1n

N2;g dgdn ¼

0

1 

Z

4ng þ

Z 0

g2

1

Z

1

0dn ¼ 0

0

Z

1n

ð4n þ g þ 1Þdgdn

0

1n

þg dn 2 0

  1 1 3 9 3 3 ¼  6n þ n2 dn ¼ n  3n2 þ n3 2 2 2 2 0 0

 3 3 ¼ 3þ ¼0 2 2 Z 1 Z 1n Z Z 1 Z 1n N3;g dgdn ¼ ð4n  u  1Þdgdn ¼ ¼

ð13Þ

1n

1

0

Z

0

Z

0

Z

0

0

1

0

0

1n

Z

1

0

ð15Þ

1

Z

1n

ð2 þ 5n þ gÞdgdn

0

 1 3 9 2g þ 5ng þ g2 dn ¼  þ 6n  n2 dn ¼ 2 2 2 0 0 0  1

 3 3 3 3 ¼0 ¼  n þ 3n2  n3 ¼  þ 3  2 2 2 2 0 Z

V under ¼

Acknowledgement

ð16Þ

In Fig. 6, the volume above the triangle D123 is equal to the volume under the triangle, which is proved as follows:

1 9  ð0:6  0:75  0:5Þ  3 ¼ 3 40 Z 0:75 Z 10:2nþ0:4 Z 1 Z ¼ N3;n dgdn þ

V above ¼

Z Z

0 0:75

0 1

Z

Z

0 10:2nþ0:4

0:75

0

ð17Þ 1n

N 3;n dgdn

ðn þ 5g  2Þdgdn

0 1n

ðn þ 5g  2Þdgdn   Z 1

1 2 2 3 2 1 dn þ dn  n2 þ n  n  2n þ ¼ 10 5 5 2 0 0:75 2 129 3 9  ¼ ¼ 640 128 40

þ

0:75 0 0:75

Z

ð18Þ

References [1] Rannacher R. A posteriori error estimation in least-squares stabilized finite element schemes. Comput Meth Appl Mech Eng 1998;166:99–114. [2] Hughes TJR, Franca LP. A new finite element formulation for computational fluid dynamics: Vii. the stokes problem with various well-posed boundary conditions: Symmetric formulations that converge for all velocity/pressure spaces. Comput Meth Appl Mech Eng 1987;65(1):85–96. [3] Klaas O, Maniatty A, Shaphard MS. A stabilized mixed finite element method for finite elasticity: Formulation for linear displacement and pressure interpolation. Comput Meth Appl Mech Eng 1999;180(1-2):65–79. [4] Ramesh B, Maniatty A. Stabilized finite element formulation for elastic-plastic finite deformations. Comput Meth Appl Mech Eng 2005;194(1):775–800. [5] Schwarz A, Schröder J, Starke G. Least-squares mixed finite elements for small strain elasto-viscoplasticity. Int J Numer Meth Eng 2009;77:1351–70. [6] Schwarz A, Schröder J, Starke G. A modified least-squares mixed finite element with improved momentum balance. Int J Numer Meth Eng 2010;81:286–306. [7] Codina R. Stabilization of incompressibility and convection through orthogonal subscales in finite element methods. Comput Meth Appl Mech Eng 2000;190(1):1579–99. [8] Cervera M, Chiumenti M, Valverde Q, Agelet de Saracibar C. Mixed linear/linear simplicial elements for incompressible elasticity and plasticity. Comput Meth Appl Mech Eng 2003;192(49-50):5249–63. [9] Agelet de Saracibar C, Chiumenti M, Valverde Q, Cervera M. On the orthogonal subgrid scale pressure stabilization of finite deformation J2 plasticity. Comput Meth Appl Mech Eng 2006;195(1):1224–51. [10] Cervera M, Chiumenti M, Codina R. Mixed stabilized finite element methods in nonlinear solid mechanics. Part I: Formulation. Comput Meth Appl Mech Eng 2010;199(1):2559–70. [11] Arnold DN, Brezzi F, Fortin M. A stable finite element for the Stokes equation. Calcolo 1984;21:337–44. [12] Baiocchi C, Brezzi F, Franca LP. Virtual bubbles and Galerkin-least-squares type methods (Ga.L.S.). Comput Meth Appl Mech Eng 1993;105:125–41. [13] Brezzi F, Fortin M. A minimal stabilisation procedure for mixed finite element methods. Numer Math 2001;89:457–92. [14] Simo JC, Rifai S. A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer Meth Eng 1990;29:1595–638. [15] Taylor RL. A mixed-enhanced formulation tetrahedral finite elements. Int J Numer Meth Eng 2000;47(1-3):205–27. [16] Cisloiu R, Lovell M, Wang J. A stabilized mixed formulation for finite strain deformation for low-order tetrahedral solid elements. Finite Elem Anal Des 2008;44(1):472–82. [17] Lamichhane BP. From the Hu-Washizu formulation to the average nodal strain formulation. Comput Meth Appl Mech Eng 2009;198:3957–61. [18] Mahnken R, Caylak I, Laschet G. Two mixed finite element formulations with area bubble functions for tetrahedral elements. Comput Meth Appl Mech Eng 2007;197:1147–65. [19] Mahnken R, Caylak I. Stabilization of bi-linear finite elements for tetrahedra with enhanced interpolation using volume and area bubble functions. Int J Numer Meth Eng 2008;75:377–413. [20] Wilson EL, Taylor RL, Doherty WP, Ghaboussi J. Incompatible displacement models. In: Fenves SJ et al., editors. Numerical and Computer Models in Structural Mechanics. New York: Academic Press; 1973. p. 43. [21] Simo JC, Armero F. Geometrically non-linear enhanced strain methods and the method of incompatible modes. Int J Numer Meth Eng 1992;33:1413–49. [22] Bonet J, Burton AJ. A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications. Commun Numer Meth Eng 1998;14:437–49. [23] Dohrmann CR, Heinstein MW, Jung J, Key SW, Witkowsky WR. Node-based uniform strain elements for three-node triangular and four-node tetrahedral meshes. Int J Numer Meth Eng 2000;47:1549–68. [24] Bonet J, Marriott H, Hassan O. An averaged nodal deformation gradient linear tetrahedral element for large strain explicit dynamic applications. Commun Numer Meth Eng 2001;17:551–61. [25] Pires FMA, de Souza Neto EA, de la Cuesta Padilla JL. An assessment of the average nodal volume formulation for the analysis of nearly incompressible solids under finite strains. Commun Numer Meth Eng 2004;20(7):569–83.

182

I. Caylak, R. Mahnken / Computers and Structures 138 (2014) 172–182

[26] Puso MA, Solberg J. A stabilized nodally integrated tetrahedral. Int J Numer Meth Eng 2006;67:841–67. [27] Hughes TJR. Equivalence of finite elements for nearly incompressible elasticity. Comput Meth Appl Mech Eng 1977;44:181–3. [28] Hughes TJR. Generalization of selective integration procedures to anisotropic and nonlinear media. Int J Numer Meth Eng 1980;15:1413–8. [29] Reese S. On a physically stabilized one point finite element formulation for three-dimensional finite elasto-plasticity. Comput Meth Appl Mech Eng 2005;194:4685–715. [30] Reese S. A large deformation solid-shell concept based on reduced integration with hourglass stabilization. Int J Numer Meth Eng 2007;69:1671–716. [31] Caylak I, Mahnken R. Mixed finite element formulations with volume bubble functions for triangular elements. Comput Struct 2011;89:1844–51.

[32] Caylak I, Mahnken R. Stabilization of mixed tetrahedral elements at large deformations. Int J Numer Meth Eng 2012;90:218–42. [33] Zienkiewicz OC, Taylor RL. The finite element method, 6th ed., vol. 1. London: McGraw-Hill; 2005. [34] Abaqus. Version 6.9-EF, 2010. [35] Chiumenti M, Valverde Q, Agelet de Saracibar C, Cervera M. A stabilized formulation for incompressible elasticity using linear displacement and pressure interpolations. Comput Meth Appl Mech Eng 2002;191:5253–64. [36] Reese S. On a consistent hourglass stabilization technique to treat large inelastic deformations and thermo-mechanical coupling in plane strain problems. Int J Numer Meth Eng 2003;57:1095–127.