A study on finite elements for capturing strong discontinuities

A study on finite elements for capturing strong discontinuities J. Oliver ∗A. E. Huespe † E. Samaniego ‡ E.T.S. Enginyers de Camins, Canals i Ports, Technical University of Catalonia Campus Nord UPC, M`odul C1, Gran Capit´an s/n 08034 Barcelona, Spain March 2, 2002 Abstract The work focuses on the presently existing families of finite elements with embedded discontinuities and explores the possibilities to obtain symmetric statically consistent finite elements that alleviate the stress locking problem. For this purpose, mixed (reduced integration) and assumed enhanced strain techniques are applied to the basic symmetric four-noded element. Numerical simulations show the effectiveness of the proposed measures. keywords: finite elements , embedded discontinuities, strong discontinuities, mixed elements, assumed enhanced strains 1 Motivation In recent years, the so called finite elements with embedded discontinuities have been object of increasing study and development [5], [8], [10], [24], [21], [13], [1], [9], [3] ,[2], [27], [23]. Its ∗Professor. e-mail: oliv[email protected]c.es †Researcher from: CIMEC/ CONICET-UNL, Argentina. e-mail: [email protected] ‡Research assistant 1 rising popularity comes from the fact that, by using these elements, a displacement discontinuity model can be introduced in the bulk of the element (and therefore regardless the orientation of its sides) in combination with appropriate propagation mechanisms. Moreover, it has been shown [13] that in strain localization scenarios, some of those elements completely overcome the well known problem of the spurious mesh size and mesh orientation dependence of the results [4]. Although there are quite different families of such elements, apparently, not all them behave the same way. A fairly comprehensive study of the different families can be found in [6]. They can be classified in three different groups: I) Symmetric statically consistent elements. Traction continuity across the discontinuity interface is introduced in a fully consistent variational environment that results in a symmetric formulation. However the way of introducing the discontinuous kinematics does not guarantee free rigid body relative motions of the two portions of the element split up into by the discontinuity1. A description of such formulation can be found in [9]. II) Symmetric kinematically consistent elements. The kinematics is introduced in a way that does not restrict the rigid body relative motions of the two portions of the element. However, traction continuity across the discontinuous interface is not guaranteed at elemental level. A typical triangular element based on this formulation can be found in [10]. III) Non-symmetrical statically and kinematically consistent elements. Both the rigid body relative motions and the traction continuity are introduced at the elemental level, but the second is introduced in a strong form that makes the resulting formulation nonsymmetrical. Typically, these are the formulations used in [13] and [1]. From the authors experiences the last family of elements iii) is the one that provides more robust an reliable results and they have used it successfully for numerical simulations in many different settings [12], [16], [17], [14], [18] ,[19], [20]. Family ii) exhibits also a robust behaviour although a much slower convergence with mesh refinement whereas family i), as it will be shown below, exhibits in many cases the well known stress locking phenomenon. However, although the fact that formulation iii) is not symmetric is neither a fundamental drawback nor a source of unbearable computational costs, the variational consistency of formulation i) agrees with the traditional finite element technology in Computational Solid Mechanics 1As it will be shown in subsequent sections this can result in stress locking behaviour. 2 and confers additional appeal upon it. More important, the symmetric elements of family i) can be theoretically derived without necessity of the tracking algorithm concept which appears to be a fundamental issue on the other two families and sets practical difficulties to generalize its use for capturing multiple cracking and branching phenomena. This is why this paper is devoted to explore and develop several new types of finite elements with embedded discontinuities belonging to this family. However, and since the aforementioned stress locking phenomena makes the basic formulation unsuitable, specific treatments are explored to face that problem. Mixed and assumed enhanced strain techniques are envisaged as appropriated remedies, so that they are introduced in the basic element and their effects on the stress locking are analyzed. 1.1 The Boundary Value problem Let us consider the body Ω of figure 1-a, undergoing a (rate of) displacement discontinuity [[ ˙u]](x,t) across the material (fixed) surface Sthat splits the body into Ω+(pointed by the unit normal nto S) and Ω−such that Ω+∪Ω−= Ω\S2. The resulting velocity, ˙u(x, t), and strain rate, ˙ε(x, t) , fields3can be expressed as in [19]: ˙u(x, t) = . ¯u(x, t) + MS[[ ˙u]](x, t); [[ ˙u]] = ˙u|x∈∂Ω+∩S −˙u|x∈∂Ω−∩S ;. ¯u =˙u∗in ∂uΩ ˙ε(x, t) = ∇S˙u(x, t) = . ¯ε |{z} regular (bounded) +δS([[ ˙u]] ⊗n)S | {z } singular (unbounded) (1) MS(x) = HS(x)−ϕ(x) ; ϕ(x)∈H1(Ω) =    1∀x∈Ω+\Ωh 0∀x∈Ω−\Ωh (2) where [[ ˙u]] stands for the velocity jump, ∂uΩ is the part of the external boundary of Ω (with outward normal ν) where displacements are prescribed to u∗(x, t), HSis the Heaviside’s jump function placed on S(HS(x) = 1 ∀x∈Ω+and HS(x)=0 ∀x∈Ω−), MSis a unit jump function, whose support is a certain domain Ωhcontaining S(see figure 1-(b)) and constructed 2Notation A\Bstands for the result of substraction of domain Bfrom domain A. 3The mathematical expressions of the resulting continuum format kinematics is referred to as strong discontinuity kinematics [19]. 3 as it is indicated in equation (2), and δSis the Dirac’s delta function placed on Sthat arises from the derivation of the discontinuous function MS(x). The corresponding quasistatic Boundary Value Problem (B.V.P.) can be described, in rate form, as the following three field, u−ε−σ, problem: FIND :         ˙u(x, t) ˙ε(x, t) ˙σ(x, t) satisfying: ∇ · ˙σ−˙ b=0in Ω\S (internal equilibrium) (a) ˙ε− ∇S˙u =0in Ω (kinematical compatibility) (b) ˙σ−˙ Σ(ε) = 0in Ω (constitutive compatibility) (c) ˙σ·ν=˙ text on ∂σΩ (external equilibrim) (d) ˙σΩ+·n−˙σΩ−·n | {z } def = [[ ˙σ]]Ω\S ·n =0on S(outer traction continuity (e) ˙σΩ+·n−˙σS·n | {z } def = [[ ˙σ]]S·n =0on S(inner traction continuity) (f) (3) where σ(x, t) stands for the stresses, b(x, t) are the specific body forces, Σ(ε) stands for the constitutive function returning the stresses in terms of the strains εand ∂σΩ is the part of the boundary of Ω where tractions text are prescribed. Equations (3)-(e) and (3)-(f) state the continuity of the traction vector T=σ·nthrough the discontinuity interface S. After explicit imposition of condition (3)-(c) B.V.P. (3) can be rewritten as a two field, u−ε, problem: : FIND :   ˙u(x, t) ˙ε(x, t) satisfying: ∇ · ˙ Σ−˙ b=0in Ω\S (internal equilibrium) (a) ˙ε−∇S˙u =0in Ω (kinematical compatibility) (b) ˙ Σ·ν=˙ text on ∂σΩ (external equilibrim) (c) ˙ ΣΩ+·n−˙ ΣΩ−·n | {z } def = [[ ˙ Σ]]Ω\S ·n =0on S(outer traction continuity (d) ˙ ΣΩ+·n−˙ ΣS·n | {z } def = [[ ˙ Σ]]S·n =0on S(inner traction continuity) (e) (4) 4 2 Non-symmetric (Petrov-Galerkin) approach A weak form of the B.V.P. (4) can be devised as follows. In view of the velocity field (1) let us consider the functional spaces of the velocities, Vu, and virtual (kinematically admissible) velocities, ¯ V0 u: Vu def ≡©η(x) = ¯η+MSα;η∈[H1(Ω)]ndim ;α∈L2(S)ª ¯ V0 u def ≡©¯η0(x)∈[H1(Ω)]ndim ;¯η0|∂uΩ=0ª(5) where ndim stands for the number of dimensions of the problem, H1(Ω) is the functional space of functions defined in Ω with square integrable first derivatives and L2(S) the functional space of square integrable functions defined in S4. Remark 1 Notice that spaces Vuand ¯ V0 udiffer not only for the homogeneous boundary condition imposed to the elements of ¯ V0 u, but also for the intrinsic nature of their elements η(x) = ¯η+MSαand ¯η0(x)respectively. This fact entails the nonsymmetric character of the resulting formulation. In view of the preceding notation the weak form of the problem can be written: PROBLEM 1 CONTINUOUS NON-SYMMETRIC PROBLEM FIND: ˙u =. ¯u +MS[[ ˙u]] . ;˙u ∈ Vu ˙ε=∇S˙u =. ¯ε+δs([[ ˙u]]⊗n)S(6) SUCH THAT: δΠu(˙ Σ;η) = ZΩ\S ˙ Σ(∇S˙u) : ∇SηdΩ−ZΩ\S ˙ b·ηdΩ + Z∂σΩ ˙ t·ηdΓ | {z } Gext = 0 ∀η∈¯ V0 u(7) Some standard calculations show that the strong form of Problem 1 is: δΠu(u;η) = 0 ⇒ ∇ · ˙ Σ−˙ b=0in Ω\S ˙ Σ·ν=˙ text on ∂σΩ [[ ˙ Σ]]Ω\S ·n= 0 on S (8) 4Roughly H1(·) contains continuous function defined in (·) with discontinuous first derivatives and L2(·) contains piecewise discontinuous bounded functions defined in (·). 5 and, therefore, only condition (4)-(e) remains to be fulfilled in the B.V.P. defined by equations (4). This condition is going to be imposed via a different procedure in the finite element formulation, essentially in strong form. In summary the B.V.P. (4) is approached as: δΠu(˙ Σ(. ¯u,[[ ˙u]]); η) = 0 ∀η∈¯ V0 u(variational/weak form) [[ ˙ ΣS]] ·n=0on S(strong form) (9) 2.1 Finite element discretization (standard non-symmetric element: U4n) Let us consider the material domain Ω discretized in a four-noded5finite element mesh with nelem elements and nnode nodes crossed by the discontinuity interface S(see figure 2-(a)). Let us assume that an available discontinuity tracking algorithm [13] determines the subset Jof the nJelements that are crossed by Sat the considered time t: J:= {e|Ωe∩ S 6=∅} ={ei,...., em,....ep, ...}(10) For every element of J, the tracking algorithm6also provides the position of the elemental discontinuity interface Se(see figure 2-(b) of length lewhich defines the domains Ω+ eand Ω− e and the nodes i+∈ {i+ 1, .., i+ n+ e}and i−∈ {i− 1, .., i− n− e}. Now, the following interpolation of the velocity field ˙u(e)inside a given element eis considered[13]: ˙u(e)(x, t) = Σi=4 i=1 N(e) i(x)˙ di(t) | {z } . ¯u(e) +M(e) S(x) [[ ˙u]]e(t) | {z } . ˜u(e) (11) where . ¯u(e)is the standard C0velocity field, interpolated by the shape functions {N(e) 1, N(e) 2, N(e) 3, N(e) 4}of the linear isoparametric quadrilateral element [28], in terms of the nodal velocities di(t) at node i. The term . ˜u(e), in equation (11), captures the singular (discontinuous) part of the velocity field (1) in terms of the elemental velocity jump [[ ˙u]]eand Me S(x) is the discrete counterpart of the unit jump function in equation (2) defined as follows: M(e) S(x) =          0∀e6∈ J H(e) S(x)−ϕ(e) (ϕ(e)= Σn+ e i+=1Ni+)   ∀e∈ J (12) 5For the sake of simplicity, from now on only two-dimensional problems will be considered. 6This tracking algorithm is a crucial ingredient of the non-symmetric formulations and constitute one of their most typical servitudes. 6 where H(e) Sis the step function. Figure 2-(c) shows the M(e) Sfunction and emphasizes its elemental support. From equations (11) and (12), the discrete (rate of) strain field reads: ˙ε(e)=∇S˙u(e)= Σi=4 i=1 (∇N(e) i⊗˙ di)S−(∇ϕ(e)⊗[[ ˙u]]e)S+δS([[ ˙u]]e⊗n)S(13) Notice that equation (13) matches the strong discontinuity kinematics (1). In order to overcome the numerical difficulties of treating with the Dirac’s delta function δSis replaced by a regularized7function δe Sdefined within the element eas: δ(e) S=µ(e) S 1 k(14) where µ(e) Sis a collocation function whose support is the domain Sk ein figure 2-(d) defined in terms of the regularization parameter k: µ(e) S(x) = 1 ∀x∈ Sk e µ(e) S(x) = 0 ∀x/∈ Sk e (15) By considering equations (14) and (15) the regularized form of the strain rate field reads: ˙ε(e)=∇S˙u(e)= Σi=4 i=1 (∇N(e) i⊗˙ di)S−(∇ϕ(e)⊗[[ ˙u]]e)S+µ(e) S 1 k(n⊗[[ ˙u]]e)S(16) In order to integrate the discontinuous terms emerging from the second term of the righthand-side of equation (16), in addition to the standard sampling points of the linear quadrilateral (PG1 to PG4 in figure 2-(e)), the element is equipped with another integration point (SSP in figure 2-(d)) placed at the center of the element and whose associated area is (see figure 2-(d)) : measure(Sk e) = kle(17) The regularization parameter khas an arbitrary small value (as small as permitted by the machine precision). In this context the inner traction continuity condition in equations (4)-(e) and (9) can be imposed on an element basis in terms of averages: ˙ ΣΩ+·n= ( ˙ ΣΩ−·n) = ˙ ΣS·n→1 ΩeZΩe ˙ Σ·ndΩ | {z } mean value on Ω\S =1 leZSe ˙ Σ·ndΩ | {z } mean value on S ⇒(18) 7This procedure and the resulting formulation have been sometimes termed in the literature the regularized strong discontinuity approach. 7 ⇒ZΩe (µ(e) S 1 k−le Ωe )n·˙ ΣdΩ = 0∀e∈ J (19) In view of the previous finite element discretization, and by resorting to the classical finite element procedures, the discrete counterpart of the B.V.P. in equations (9) can be written as follows: PROBLEM 2 DISCRETE NON-SYMMETRIC PROBLEM GIVEN: Vh u def ≡nηh(x) = Σi=nnode i=1 Ni(x)ηi+ Σe∈J M(e) S(x)αeo ¯ Vh0 u def ≡nηh0(x) = Σi=nnode i=1 Ni(x)η0 i;η0 i|∂uΩ=0o(20) FIND: ˙uh= Σi=nnode i=1 Ni˙ di+ Σe∈J M(e) S(x) [[ ˙u]]e;˙u h∈ Vh u ˙εh=∇S˙uh= Σi=nnode i=1 (∇Ni⊗˙ di)S+ Σe∈J ³[µ(e) S1 kn− ∇ϕ(e)]⊗[[ ˙u]]e´S(21) SUCH THAT: δΠu(˙ Σ;ηh) = Pe=nelem e=1 RΩe∇Sηh:˙ Σ(εh)dΩ−Gext = 0 ∀ηh∈¯ Vh0 u [[ ˙ Σ]]S·n=0→RΩe(µ(e) S1 k−le Ωe)n·˙ ΣdΩ = 0∀e∈ J (22) The structure of equations (22) corresponds to a typical Petrov-Galerkin residual weighting procedure [28] of the original B.V.P. (4). 2.1.1 2D implementation For the two-dimensional case, in a cartesian coordinate system (x, y), using the vector format for the strains {ε}= [εxx, εyy,2εxy]Tand the stresses {σ}= [σxx, σyy,σxy]T(where (·)Tstands for the transpose of (·)), considering the four noded quadrilateral as underlaying element and using the standard finite element B-format [28] equations (22) yield: δΠu(˙ Σ;ηh) = 0 [[ ˙ Σ]]S·n=0   → e=nelem [ e=1 ·ZΩe B(e)T·˙ Σ({ε}(e))dΩ−˙ Fext(e)¸=0(23) where ˙ Fext(e)stands for the classical elemental external forces vector and: {˙ε}(e)=B∗(e)·˙ d(e) ˙ d(e)=h˙ d1,˙ d2,˙ d3,˙ d4,[[ ˙u]]eiT(24) 8 B(e)=hB(e) 1,B(e) 2,B(e) 3,B(e) 4,G(e)i B∗(e)=hB(e) 1,B(e) 2,B(e) 3,B(e) 4,G∗(e)i B(e) i=     ∂xN(e) i0 0∂yN(e) i ∂yN(e) i∂xN(e) i      G(e)= (µ(e) S1 k−le Ωe)     nx0 0ny nynx      G∗(e)=µ(e) S1 k      nx0 0ny nynx      −     ∂xϕ(e)0 0∂yϕ(e) ∂yϕ(e)∂xϕ(e)      (25) where Sstands for the classical assembling operator and n= [nx, ny]T. Remark 2 Notice that matrices B(e)and B∗(e)in equation (25) differ in the terms G(e)6= G∗(e).This fact makes the resulting tangent stiffness matrix non-symmetrical as could be expected from the original non-symmetric character of the approach stated in Remark 1. Remark 3 The structure of equations (23) to (25) suggests the introduction of an internal additional fifth node for each element e, that is activated only for the elements crossed by the discontinuity interface (e∈ J ) and whose corresponding degrees of freedom and associated shape function are, respectively, the displacement jumps [[u]]eand M(e) Sin equations (11) and (12). Since the support of M(e) Sis only Ωe, those internal degrees of freedom can be eventually condensed at the elemental level and removed from the global system of equations. 3 Symmetric assumed enhanced strain approach The assumed enhanced strain methods, which can be considered a particular case of the more general assumed strain methods or mixed methods [25], provide a different setting to approach displacement discontinuities. In next sections the corresponding formulation is presented following the guidelines in ([26]). 9 a) Decrease the polynomial degree of the term ∇S˙u(e)to zero (constant). b) Increase the polynomial degree of the enhancement term to one (linear polynomial). Strategy a) is followed in the mixed approach presented in Section 5 whereas strategy b) leads to the strain re-enhancement methodology presented in Section 612. 5 Mixed approach 5.1 Assumed strain and stress fields Let us consider the following functional spaces: Vu def ≡©η(x)∈[H1(Ω)]2ª(a) V0 u def ≡©η0(x)∈[H1(Ω)]2;η0|∂uΩ=0ª(b) Vε def ≡ {ξ(x) = ¯ ξ+δS(α⊗n)S;¯ ξij ∈L2(Ω) ; αi∈L2(S)}(c) Vσ def ≡nτ(x) ; τij ∈L2(Ω ) ; [[τ]]Ω\S ·n= [[τ]]S·n=0o(d) (52) The variational three field (u−ε−σ) mixed problem can we written as: PROBLEM 5 CONTINUOUS MIXED PROBLEM FIND: ˙u(x, t)˙u ∈ Vu ˙ε(x, t) = . ¯ε+δs(˙ β⊗n)S. ε∈ Vε ˙σ(x,t)˙σ∈ Vσ (53) SUCH THAT δΠu(˙σ;η) = RΩ\S ˙σ:∇SηdΩ−Gext = 0 ∀η∈ V0 u δΠε(˙u,˙ε;τ) = RΩ¡˙ε− ∇S˙u¢:τdΩ = 0 ∀τ∈ Vσ δΠσ(˙σ,˙ Σ;ξ) = RΩ³˙σ−˙ Σ´:ξdΩ = 0 ∀ξ∈ Vε (54) Standard calculations on equations (54) lead to the corresponding strong forms: i) δΠu(˙σ;η) = 0 ⇒ ∇ · ˙σ−˙ b=0in Ω\S ˙σ·ν=˙ ton ∂σΩ [[ ˙σ]]Ω\S ·n=0on S (55) 12Indeed, these strategies for alleviating the stress-locking problem can only be applied to finite elements sensitive to mixed /assumed strain or assumed enhanced strain technologies, like the four-noded quadrilateral. For instance, strategies a) and b) cannot be applied, at elemental level, to linear triangles. 16 ii) δΠε(˙u,˙ε;τ) =0 ⇒RΩ¡˙ε− ∇S˙u¢:τdΩ = 0 ⇒˙ε=∇S˙u in Ω(56) iii) δΠσ(˙σ,˙ Σ;ξ) = 0 ⇒RΩ³˙σ−˙ Σ´:ξdΩ⇒˙σ=˙ Σin Ω(57) Therefore, equations (3)-(a) to (3)-(e) of the original B.V.P. are fulfilled in weak form by the variational equations (54) whereas equation (3)-(f) is fulfilled from the particular choice made for the space Vσin equation (52)-(d). 5.2 Finite element discretization (constant stress/strain mixed element: M4n) For the 4-noded element discretization of figure 2-(a) let us consider the following discrete counterpart of the spaces in equation (52): Vh u def ≡©ηh(x) = Σi=nnode i=1 Ni(x)ηiª(a) V0 u def ≡©ηh(x) =Σi=nnode i=1 Ni(x)ηi;ηi|∂uΩ=0ª(b) Vh ε def ≡ {ξh(x) =Σe=nelem e=1 χe(x)(¯ ξe+µ(e) S1 k(αe⊗n)S)};χe(x) =    1 for x∈Ωe 0 otherwise (c) Vh σ def ≡ {τh(x) =Σe=nelem e=1 χe(x)τe}(d) (58) and the corresponding discrete problem: PROBLEM 6 DISCRETE MIXED PROBLEM FIND: ˙uh= Σi=nnode i=1 Ni˙ di;˙uh∈Vh u(a) ˙εh= Σe=nelem e=1 (χe(x)˙ ¯εe+µ(e) S1 kχe(x)(˙ βe⊗n)S). εh∈Vε(b) ˙σh=Σe=nelem e=1 χe(x). σe˙σh∈Vσ(c) (59) SUCH THAT: δΠu(˙σ;ηh=RΩ˙σh:∇SηhdΩ−Gext = 0 ∀ηh∈Vh0 u(a) δΠε(˙uh,˙εh;τh=RΩ³˙εh− ∇S˙uh´:τhdΩ = 0 ∀τh∈Vh σ(b) δΠσ(˙σh,˙ Σh;ξh) = RΩ³˙σh−˙ Σh´:ξhdΩ = 0 ∀ξh∈Vh ε(c) (60) 17 By inserting the strain rate field ˙εhof equation (59)-(b) into equation (60)-(b) one gets: δΠε(˙uh,˙εh) = 0 ⇒ e=nelem X e=1 ZΩe³˙εh− ∇S˙uh´:τedΩ = (61) = e=nelem X e=1 µΩe . ¯εe+le(˙ βe⊗n)S−ZΩe ∇S˙u hdΩ¶:τe= 0 ∀τe⇒ ⇒µΩe . ¯εe+le(˙ βe⊗n)S−ZΩe ∇S˙u hdΩ¶= 0 e∈ {1...nelem} ⇒ and solving for . ¯εeand, then, for . ¯εhin equation (59)-(b): . ¯εe=1 ΩeZΩe ∇S˙uhdΩ−le Ωe (˙ βe⊗n)Se∈ {1...nelem} ⇒ (62) . ¯εh= e=nelem X e=1 χe(x)³1 ΩeZΩe ∇S˙uhdΩ | {z } def =∇S˙u(e) −(le Ωe −µ(e) S 1 k)( ˙ βe⊗n)S´ where ∇S˙u(e)stands for the mean value of ∇S˙uh(x) inside the element e. In summary: δΠε(˙uh,˙εh:τh) =0 →˙εh=Pe=nelem e=1 χe(x)[∇S˙u(e)+ (µ(e) S1 k−le Ωe)( ˙ βe⊗n)S](63) Now, from equation (60)-(c) and the expansion of ξh(x) in equation (58)-(c): δΠσ(˙σh,˙ Σ;ξh) = 0 ⇒ e=nelem X e=1 ZΩe³˙σe−˙ Σ´:ξhdΩ = (64) = e=nelem X e=1 ZΩe³˙σe−˙ Σ´: (¯ ξe+µ(e) S 1 k(αe⊗n)S)dΩ = = e=nelem X e=1 µ·Ωe˙σe−ZΩe ˙ ΣdΩ¸:¯ ξe+·le˙σe·n−ZSe ˙ Σ·ndS¸·αe¶= 0 ∀¯ ξe∀αe ⇒hΩe˙σe−RΩe ˙ ΣdΩi= 0 hle˙σe·n−RSe ˙ ΣdS · ni   e∈ {1...nelem}(65) and solving for ˙σein equation (65): ˙σe=1 ΩeRΩe ˙ ΣdΩ = ˙ ΣΩe(a) ˙σe·n=1 leRSe ˙ ΣdS · n=˙ ΣS·n(b) (66) 18 where ˙ ΣΩe=1 ΩeRΩe ˙ ΣdΩ and ˙ ΣS=1 leRSe˙ ΣdSare, respectively, the mean values of ˙ Σ(ε(x)) in Ωeand Se.From equation (66) it trivially follows: ˙ ΣΩe·n=˙ ΣS·n(67) which states the inner traction continuity condition (4)-(e) in terms of the mean values of ˙ Σ. Equation (66)-(b) can be rewritten in a more convenient format: 1 leZSe ˙ ΣdS · n−1 ΩeZΩe ˙ ΣdΩ·n=ZΩe (µ(e) S 1 k−le Ωe )˙ Σ·ndΩ = 0(68) In summary, from equations (64) and (68): δΠσ(˙σh,˙ Σ;ξh) = 0 →RΩe(µ(e) S1 k−le Ωe)˙ Σ·ndΩ = 0∀e∈ {1...nelem}(69) Finally, standard algebraic operations on equation (60)-(a), taking into account equations (66) lead to: δΠu(˙σh;ηh) = 0 → e=nelem [ e=1 ZΩe . σe:∇SηhdΩ−Gext = = e=nelem [ e=1 ZΩe ˙ ΣΩe:∇SηhdΩ−Gext = = e=nelem [ e=1 ˙ ΣΩe:ZΩe ∇SηhdΩ | {z } = Ωe∇Sη(e) −Gext = (70) = e=nelem [ e=1 Ωe˙ ΣΩe | {z } RΩe ˙ ΣdΩ :∇Sη(e)−Gext = 0 ⇒ ⇒ e=nelem S e=1 RΩe∇Sηh:˙ Σ(εh)dΩ−Gext = 0 ∀ηh∈ V0 u(71) Remark 6 Notice from equation (63) that the mixed approach recovers the same enhanced strain counterpart . ˜ε= (µ(e) S1 k−le Ωe)( ˙ βe⊗n)Spostulated in equation (37) for the assumed enhanced strain approach. However in the mixed approach the term ∇S˙u(e)in equation (63) is constant and, therefore, so is the component ∇S˙u(e)≈ ∇S˙u(e)in equation (51). This fact is expected to facilitate the cancellation of the elemental bulk strain ˙ε(e) Ω\S =∇S˙u(e)−le Ωe(˙ βe⊗n)S and contribute to alleviating the stress locking phenomenon. 19 5.2.1 2D implementation Again, for the two-dimensional case, and the four noded quadrilateral element, the B-format formulation of the finite element method, from equations (63), (69) and (71) reads: δΠu(˙ Σ;ηh) = 0 δΠσ(˙σh,˙ Σ;ξh)=0 ⇒ e=nelem [ e=1 ·ZΩe ¯ B(e)T·n˙ ΣodΩ−˙ Fext(e)¸=0(72) {˙ε}(e)=¯ B(e)·˙ d(e) ˙ d(e)=h˙ d1,˙ d2,˙ d3,˙ d4,˙ βeiT(73) ¯ B(e)=h¯ B(e) 1,¯ B(e) 2,¯ B(e) 3,¯ B(e) 4,G(e)i ¯ B(e) i=     ∂xN(e) i0 0∂yN(e) i ∂yN(e) i∂xN(e) i      ; (·)(e)def =1 ΩeZΩe (·)dΩ | {z } mean value of (·) on Ωe G(e)= (µ(e) S1 k−le Ωe)     nx0 0ny nynx      (74) Remark 7 For practical purposes the computation of the mean elemental values (·)(e)referred to in equation (74) can be approximately (or even exactly, for undistorted elements) computed, for the quadrilateral element, through sampling at the center of the element. Consequently by comparison of equations (47) and (74) the mixed M4n element can be retrieved from the symmetric S4n element by reducing the four regular sampling points (PG1 to PG4 in figure 2-(e)) to one placed at the center of the element ( (RSP in that figure). Therefore, the classical link between mixed elements and reduced integration ([11]) is also recovered for finite elements with embedded discontinuities. 6 Strain re-enhancement: element E4n Next strategy is based on providing new terms to the enhanced strain component ( . ˜εin equation (28)) of the symmetric element S4n. Two conditions are required on this enhancement : i) Fulfill the orthogonality condition (27) : 20 ZΩ . ˜ε:τdΩ=0 ∀. ˜ε∈ V˜ε∀τ∈Vσ(75) ii) Include linear polynomial terms to contribute to alleviate the stress locking phenomenon. With these conditions in mind, the following strain enhancement is proposed: . ˜εh(e) = (µ(e) S 1 k−le Ωe )(˙ βe⊗n)S | {z } . ˜ε(e) 1 +1 Js˙ Se+1 Jt˙ Te | {z } . ˜ε(e) 2 (76) where sand tstand for the isoparametric coordinates of the standard 4-noded quadrilateral and Jis the jacobian of the isoparametric transformation relating differential areas in the regular and isoparametric spaces through: dΩ = J ds dt (77) In equation (76) . ˜ε(e) 1is the basic strain enhancement, already present in the basic S4n element, and . ˜ε(e) 2is a re-enhancement of the strain field that supplies the required linear polynomial components to the elemental bulk strain. The values n˙ Seo= [ ˙ Sxx,˙ Syy,˙ Sxy]T eand n˙ Teo= [ ˙ Txx,˙ Tyy,˙ Txy]T eare (constant) intensity factors that constitute six (for the 2D problem) additional degrees of freedom of the element. From expression (76) it is clear that the orthogonality condition (75) is fulfilled for the element-wise constant assumed stress field in equation (35) since: ZΩe . ˜ε(e) 1:τedΩ = ZΩe (µ(e) S 1 k−le Ωe )(˙ βe⊗n)S:τedΩ = =ZΩe (µ(e) S 1 k−le Ωe )dΩ | {z } le−le= 0 (˙ βe⊗n)S:τe= 0 (78) ZΩe . ˜ε(e) 2:τedΩ = ZΩe 1 J(s˙ Se+t˙ Te) : τedΩ = (79) =Z+1 −1Z+1 −1 s dsdt | {z } = 0 ˙ Se:τe+Z+1 −1Z+1 −1 t dsdt | {z } = 0 ˙ Te:τe= 0 21 6.1 2D implementation In view of the preceding formulation, the 2D implementation of element E4n becomes similar to that of the element S4n in equations (45) and (47) but including the additional enhanced strain terms in equation (76). That is: δΠu(˙ Σ;ηh) = 0 δΠσ(˙σh,˙ Σ;ξh) = 0 ⇒ e=nelem [ e=1 ·ZΩe B(e)T·n˙ ΣodΩ−˙ Fext(e)¸=0(80) {˙ε}(e)=B(e)·˙ d(e) ˙ d(e)=h˙ d1,˙ d2,˙ d3,˙ d4,[˙ βe,n˙ Seo,n˙ Teo]iT(81) B(e)=hB(e) 1,B(e) 2,B(e) 3,B(e) 4,G(e)i B(e) i=     ∂xN(e) i0 0∂yN(e) i ∂yN(e) i∂xN(e) i      G(e)=     γnx0 0γny γnyγnx | {z } . ˜ε1 s0 0 0s0 00s t0 0 0t0 00t      | {z } . ˜ε2 γ= (µ(e) S1 k−le Ωe) (82) where the eight internal degrees of freedom [ ˙ βe,n˙ Seo,n˙ Teo] can be condensed at elemental level. 7 Numerical tests. 7.1 Homogeneous plate The basic test of Section 4 and figure 3 is now repeated using the modified elements M4n and E4n. The results, again in terms of σxx −δcurves, are presented in figure 4 together with the ones obtained with the original element S4n and the exact result (U4n). It can be observed the dramatic reduction in the stress-locking effect obtained with the new elements which prove the effectiveness of the adopted strategies. Besides, one can observe that the results obtained with element M4n (mixed strategy) and element E4n (re-enhancement 22 strategy) are very similar to each other, which shows that the improvement reasoning based on the cancellation of the bulk strain, stated in Section 4 and implemented in two different strategies, was essentially correct. 7.2 Notched specimen In order to assess the performance of the proposed strategies in more complex problems, the test of figure 5-(a) is also considered. It is a notched specimen whose experimental testing is reported in [7]. The couples of forces F1and F2are progressively applied, as it is indicated, in figure 5-(b) in order to induce a mixed mode crack which develops from the notch tip with an inclination angle of 71o(see figure 5-(c). For the numerical simulation the same material model (isotropic continuum damage with linear softening) as in the previous case has been adopted, with the following parameters : the elastic modulus and Poisson’s ratio are E= 30580.[Mpa] and ν= 0.2 respectively, the fracture energy Gf= 100.[N/m] and the tensile strength is assumed to be σu= 3.[MPa].The width of the specimen is t= 50.8[mm]. The results, in terms of the applied force versus the CMOD, obtained with the U4n element, the original S4n element and the modified symmetric M4n and E4n elements, as well as the experimental values from the above reference, are presented in figure 5-(d). Again the dramatic improvement, in terms of the numerical stress locking effect, obtained with the devised strategies in comparison with the original symmetric element S4n can be observed. Also the iso-displacement curves shown in figure 5-(e) display a modelled discontinuity that matches very well the experimental crack. 7.3 Three point notched beam We consider now the three point notched concrete beam, in plane strain, shown in figure 6-(a). Experimental data for this test have been presented by Peterson [22]. The concrete post-critical behavior is modeled by the same isotropic continuum damage model than in ([15]). The material parameters are taken: E= 30580.[Mpa], ν= 0.2 and the fracture energy Gf= 126.[N/m]. The tensile strength is taken σu= 3.4[MPa].The width of the specimen is t= 100.[mm].The mesh consists of 615 quadrilaterals refined near the notch tip. A crack opening in mode I is expected to develop from the notch tip, propagating in the vertical direction through the beam, as shown in figure 6-(b). The plots in figure 6-(c)(d) 23 correspond, respectively, to the load Fvs. the vertical displacement δcurves using linear and exponential softening laws. Once again it can be noticed the striking behavior already observed in the previous simulations: the M4n and E4n elements dramatically reduce the stress locking phenomenon exhibited by the original S4n element, and provide a response that is very close to the one obtained with the non-symmetric (U4n) element. 8 Concluding remarks Throughout this work we have explored the possibility of using statically consistent symmetric elements with embedded discontinuities to capture strong discontinuities. As a matter of fact it has been shown that the use of mixed (M4n element) or assumed enhanced strain (E4n element) techniques contribute to substantially alleviate the strain locking phenomenon appearing in the original S4n element. In essence these techniques are intended for recovering the capacity of the non-symmetric element U4n to reproduce rigid body motions of the portions of the elements split up into by the discontinuity, while keeping the symmetric character and variational consistency of the original symmetric S4n element. An important issue, that has not been addressed so far, is the stability properties of the derived elements. It is well known in the literature that the use of mixed element techniques [28] and assumed enhanced strain techniques [26] may introduce spurious singular modes whose propagation through the finite element mesh can destroy both the reliability of the solution and the convergence of the non-linear problem. For instance: it is known that the four-noded element with reduced integration (one sampling point) considered here in element M4n, leads to the development and propagation of the so called hour-glass spurious modes [28]. Also the linear re-enhancement modes, . ˜ε(e) 2in equation (76), considered for the E4n element do not fulfill the stability condition (V˜ε∩ Vε={0}) [26]. However, there is a crucial aspect in the way that the elements derived in this work are implemented, that makes their stability behaviour very different from the classical ones. Throughout this work, those finite element formulations have been derived and presented for the B.V.P. in rate (incremental) form (see equations (4),(45),(72) or (80)). This allows to consider the problem as a sequence, along time, of incremental problems, each one having its own finite element formulation. As a matter of fact the implementation is done in such a way 24 that the modifications on the basic element S4n, that lead to the derived elements M4n and E4n, are effective only for the band of elements that captures the discontinuity and only beyond the time that this discontinuity appears . The elements outside this band behave as the original S4n since they do not exhibit the stress locking problem that the modified elements try to overcome. This decreases the number of elements affected by the reduced integration or reenhancement techniques to a single band of essentially one element bandwidth. Consequently the development and propagation of spurious instability modes are extremely restrained. In addition, the computational costs associated to the presence of new degrees of freedom (for the re-enhancement strategy) are then very small. Although specific stability analysis have not been conducted in this paper, and they are left for subsequent works, during the numerical simulations presented above, and others carried out during this study, the possible instability modes have not been observed and have not placed special difficulties on the numerical procedure. However, the authors are aware that this can not be generalized to any type and size of the problems, and that specific studies on the stability issue should be carried out in the future. Also the possibility of using some of the developed elements to dispense with the discontinuity tracking algorithm should be explored in subsequent works. References [1] F. Armero and K. Garikipati. An analysis of strong discontinuities in multiplicative finite strain plasticity and their relation with the numerical simulation of strain localization in solids. Int.J. Solids and Structures, 33(20–22):2863–2885, 1996. [2] T. Belytschko, N. Moes, S. Usui, and C. Parimi. Arbitrary discontinuities in finite elements. Int. J. Numer. Meth. Engng., (50):993–1013, 2001. [3] C.C.Celigoj. On strong discontinuities in anelastic solids. a finite element approach taking a frame indifferent gradient of the discontinuous displacements. Int. J. Numer. Meth. Engng., (49):769–796, 2000. [4] R. de Borst, L. J. Sluys, H. B. Muhlhaus, and J. Pamin. Fundamental issues in finite element analyses of localization of deformation. Engineering Computations, 10:99–121, 1993. 25 0 0.4 0.8 1.2 U4n(non-symmetric) U4n(non-symmetric) M4n(mixed symmetric) M4n(mixed symmetric) E4n (enhanced symmetric) E4n (enhanced symmetric) S4n(symmetric) S4n(symmetric) E= 30580.MPa = 0.20 =0.1mThickness n F[kN] F[kN] 0.1m 1. m1. m 0.2 m F d (a) (b) (d) Experiment Experiment 0 0.4 0.8 1.2 00.2 0.4 0.6 0.8 1. 00.2 0.4 0.6 0.8 1. (c) d[x10 m] -3 d[x10 m] -3 G = 126.Nw/m = 3.4 MPas f u Figure 6: Three points notched beam: a) geometrical model; b) discontinuity path; c) results for linear softening; d) results for exponential softening. 32