A mixed three-field FE formulation for stress accurate analysis including the incompressible limit

A mixed three-…eld FE formulation for stress accurate analysis including the incompressible limit M. Chiumenti, M. Cervera and R. Codina International Center for Numerical Methods in Engineering (CIMNE) Universidad Politécnica de Cataluña (UPC) Edi…cio C1, Campus Norte, Gran Capitán s/n, 08034 Barcelona Spain. e-mail: [email protected], web page: http://www.cimne.com Keywords: Stress accurate, incompressible limit, Mixed three-…eld …nite element technology, Variational Multi Scale (VMS) stabilization. Abstract In previous works, the authors have presented the stabilized mixed displacement/pressure formulation to deal with the incompressibility constraint. More recently, the authors have derived stable mixed stress/displ-acement formulations using linear/linear interpolations to enhance stress accuracy in both linear and non-linear problems. In both cases, the Variational Multi Scale (VMS) stabilization technique and, in particular, the Orthogonal Subgrid Scale (OSS) method allows the use of linear/linear interpolations for triangular and tetrahedral elements bypassing the strictness of the Inf-Sup condition on the choice of the interpolation spaces. These stabilization procedures lead to discrete problems which are fully stable, free of volumetric locking or stress oscillations. This work exploits the concept of mixed …nite element methods to formulate stable displacement/stress/pressure …nite elements aimed for the solution of nonlinear problems for both solid and ‡uid …nite element (FE) analyses. The …nal goal is to design a …nite element technology able to tackle simultaneously problems which may involve isochoric behaviour (preserve the original volume) of the strain …eld together with high degree of accuracy of the stress …eld. These two features are crucial in nonlinear solid and ‡uid mechanics, as used in most numerical simulations of industrial manufacturing processes. Numerical benchmarks show that the results obtained compare very favourably with those obtained with the corresponding mixed displacement/pressure formulation. 1 1 Introduction This work presents a novel …nite element technology with enhanced stress accuracy and, at the same time, able to deal with the fully incompressible behavior. Stress accuracy and performance in the incompressible limit are two requirements which often coexist when addressing the numerical simulation of di¤erent industrial manufacturing processes such as metal forming, forging, extrusion, friction stir welding, cutting or machining operations, among many others. The accuracy of the solution obtained by the numerical simulation of such industrial processes is directly related to the ability of the …nite element technology adopted to deal with complex phenomena such as strain localization, the formation of shear bands, the prediction of crack propagation or the isochoric behavior of the inelastic (plastic) strains (fully deviatoric response of metallic materials under large deformations). In the literature, the incompressible limit case and stress accuracy enhancement are generally treated separately. The problem posed by the incompressibility is to avoid the so called volumetric locking, an undesirable e¤ect exhibited by …nite element approximations based on the standard Galerkin approach. Successful strategies to avoid volumetric locking based on mixed formulations ([29], [2]), enhanced assumed strain methods ([33], [34], [7], [32], [28]) and nodal pressure and strain averaging ([24], [8], [9], [10], [35]) can be found in the literature. Techniques based on the Variational Multi Scale (VMS) approach proposed by Hughes [26] have been applied in the context of solid mechanics in strain localization problems (see [27]). More recently, a method based on the mixed displacement/pressure formulation with Orthogonal Sub Scales (OSS) stabilization technique, introduced by Codina [22], has been applied to incompressible elasticity by the authors (see [20]). E¤ectiveness and robustness of this mixed displacement/pressure technique encouraged the authors to extend this approach to non-linear problems (see [13], [21] and [1]) and to strain localization analysis using both J2 plasticity and J2 damage constitutive models ([14], [15] and [16]). This FE technology has shown a good performance in the case of softening behavior of the material once the elastic limit is reached. In this case, strains concentrate into slip-lines allowing for sliding movements and the mixed formulation is able to capture such strain localization with practically mesh independent solutions. However, even if the incompressible limit has been successfully tackled, the accuracy of the stress …eld approximation was still an open issue solved by local mesh re…nement (see [31]). Lately, the authors have proposed a mixed stress/displacement formulation which uses linear/linear interpolations for both master …elds. Also in this case, the strictness of the inf-sup condition [11], when the standard Galerkin method is applied to mixed elements, imposes severe restrictions on the compatibility of the interpolations used for the displacement and the stress …elds (see [30] and [3, 4] for the analysis of admissible elements in linear elasticity). This di¢ culty can be circumvented again adopting a VMS approach which adds a consistent (reduced upon mesh re…nement) residual-based stabilization to the original problem. This mixed stress/displacement …nite element technology has demonstrated enhanced stress accuracy in both linear and non-linear analysis as well as the ability to capture stress concentrations and strain localizations with the guarantee of stress convergence upon mesh re…nement. This is an essential requirement which cannot be accomplished in a point-wise manner using the standard Galerkin displacement-based formulation. An accurate estimation of the stress …eld even in the strain localization zone drives the crack propagation without the help of ad-hoc tracking algorithms. 2 The present work makes a step forward introducing a mixed three-…eld formulation based on displacement/stress/pressure elements with linear (or, in general, equal) interpolations for all master …elds. The only requirement is the split of the constitutive equation into volumetric and deviatoric parts (more details in Section 3). Once more, the stabilization technique to overcome the Inf-sup condition is presented in terms of the VMS method. The di¤erent assumptions and approximations used to derive this novel …nite element technology are proposed in a very general format, applicable to 2D and 3D problems. Section 4 deals with the implementation and computational aspects. Finally, Section 5 shows the performance of the proposed formulation comparing with the well established mixed displacement-pressure formulation. 2 The continuum problem statement Let us denote by an open and bounded domain in Rndim where ndim is the number of dimensions of the space, and @its boundary. The boundary @is split into @uand @t, being @ = @u[@t such that the prescribed displacements, u, are speci…ed on @u(Dirichlet boundary conditions) and the prescribed tractions,  t, are applied on @t(Neumann boundary conditions). The continuum mechanical problem of linear elasticity to be considered is de…ned by the following system of equations: r  +b=0(1) C:"=0(2) " rsu=0(3) These are 3equations with 3unknown …elds: the displacements, u(x)and both the stress and the strain …elds, (x)and "(x), respectively, de…ned at each point, x, of the continuum. Eq. (1) is the balance of momentum equation, where brepresents the external load per unit of volume and r  ()is the divergence operator. Eq. (2) is the constitutive equation for linear elasticity. Note that Eq. (2) may also correspond to the secant form of a non-linear constitutive equation, where Cis the 4th order (secant) constitutive tensor. Finally, Eq. (3) is the kinematic equation (in the hypothesis of in…nitesimal strains), where rs() = 1 2hr() + r()Tidenotes the symmetric gradient operator and r()is the gradient operator. There are di¤erent alternatives to solve problem (1 3) with the corresponding appropriate boundary conditions described. The classical displacement-based formulation is obtained by substituting Eqs. (2) and (3) into Eq. (1). The result is Navier’s equation r  (C:rsu) + b=0(4) which is written in terms of the displacement …eld only. Alternatively, the mixed u=formulation uses both stresses and displacements as master …elds as: r  +b=0(5) C:rsu=0(6) obtained by substituting Eq.(3) into Eq.(2). 3 3 The volumetric/deviatoric split The objective of this section is the split of the constitutive equation (2) into its volumetric and deviatoric parts. This is possible in the cases of linear elasticity (compressible or incompressible), J2-plasticity (both small and large strain hypotheses), isotropic damage, Newtonian and non-Newtonian ‡ows as Norton-Ho¤, Sheppard-Wright, Bingham visco-plastic ‡ow, among others. The volumetric/deviatoric split is the starting point to develop a formulation able to tackle the incompressible limit. 3.1 Volumetric and deviatoric operators Let us de…ne the 4th rank volumetric and deviatoric tensors, Vand P, as follows: V=1 3(II)(7) P=I1 3(II)(8) P+V=I(9) where I= [ijkl]and I= [ij ]are the 4th rank and the 2nd rank identity tensors, respectively (ij is Kronecker’s delta). Vand Pcan be also thought as operators acting on second order tensors by taking double contraction with them. In this case, Vand Pare orthogonal projectors, that is: V:P=P:V= 0 (10) P:P=P(11) V:V=V(12) 3.2 Split of stress and strain tensors Using Vand Poperators, it is possible to extract the spheric and the deviatoric parts of a generic 2nd order tensor. Particularly, when applied to the stress tensor, , the result is: V:=1 3(II) : =pI(13) P:=I1 3(II):=pI=s(14) where p() = 1 3(:I) = 1 3tr ()is the pressure and s() = P:=dev ()are the deviatoric stresses. On one hand, the stress tensor, , is symmetric and consists of 6independent components. On the other hand, the deviatoric stress tensor, s, is also de…ned by 6components but only 5of them are independent, as the deviatoric stresses must respect the constraint: tr (s) = 0. This poses the question of how to select a frame independent (not unique) basis for the deviatoric stress tensor. This given, the stress tensor can be rebuilt adding both components of the split as: =pI+s(15) In a similar way, it is possible to split the strain tensor, ", as: 4 V:"=1 3(II) : "=1 3evol I(16) P:"=I1 3(II):"="1 3evol I=e(17) where evol =tr (")is the volumetric deformation and e=dev (")accounts for the distortions. The resulting split format of the strain tensor is: "=1 3evol I+e(18) 3.3 Split of the kinematic equation Within the hypothesis of in…nitesimal strains the kinematic equation is expressed as: "=rsu(19) Applying the volumetric/deviatoric operators, equation (19) is split as follows: evol =r  u(20) e=P:rsu=rsu1 3(r  u)I(21) Adding the volumetric and the deviatoric components, the kinematic equation is rebuilt as: "=1 3(r  u)I+P:rsu(22) 3.4 Split of the constitutive equation Let us assume that the constitutive relationship between stresses and strains can be expressed, in secant form, as: =C:"(23) Then the spheric and the deviatoric parts of the constitutive tensor, Cvol and Cdev, respectively, are obtained as: Cvol =V:C(24) Cdev =P:C(25) C=Cvol +Cdev (26) Introducing the split of stresses and strains ((15) and (18), respectively), the constitutive relationship in (23) can be written as: (pI+s) = Cvol +Cdev:1 3evol I+e(27) 5 that is: p=Cvol evol (28) s=Cdev :e(29) which are the volumetric and the deviatoric counterparts of the original constitutive equation, being Cvol =1 9I:Cvol :Ithe bulk modulus of the material. Particularizing to linear isotropic elasticity, the constitutive tensor is given by: C=K(II)+2GI1 3(II)(30) where Kis the (elastic) bulk modulus and Gis the (elastic) shear modulus. Therefore, the spheric and the deviatoric parts of the elastic constitutive tensor (30) are: Cvol =V:C= 3KV=K(II)(31) Cdev =P:C= 2GP= 2GI1 3(II)(32) and Cvol =K. Making use of the compliance (‡exibility) tensor, D=C1, it is possible to express the strains in terms of the stresses as: "=D:(33) Therefore, the spheric and the deviatoric components of the ‡exibility tensor (particularized to linear elasticity) are: Dvol =V:D=1 3KV(34) Ddev =P:D=1 2GP(35) D=Dvol +Ddev (36) Introducing the split format of the stresses (Eq. (15)), the strains (Eq. (18)) and the split of the ‡exibility tensor (Eq. (36)) into Eq. (33), the result is: 1 3evol I+e=Dvol +Ddev: (pI+s)(37) where it is possible to identify the volumetric and the deviatoric expressions of the original constitutive equation (33): evol =Dvol p(38) e=Ddev :s(39) being Dvol =I:Dvol :Ithe compressibility modulus. 6 Finally, making use of the kinematic equations (20) and (21), the previous relations translate into: r  u=Dvol p(40) P:rsu=Ddev :s(41) and, in the incompressible limit Dvol !0, they reduce to: r  u= 0 (42) P:rsu=Ddev :s(43) Observe that the constitutive relationship in the split format above consists of 6equations (the same as for equation (27)). In fact (28) is a single equation while (29) develops into 5equations. The particularization to linear elasticity is: r  u=p K(44) P:rsu=1 2Gs(45) and, in the incompressible limit, K! 1, they reduce to: r  u= 0 (46) P:rsu=1 2Gs(47) 4 The u=s=p three-…eld formulation In this section, a novel three-…eld formulation is introduced. The objective is the de…nition of a general framework, which includes the well-known mixed u=p formulation and the mixed u= formulation as particular cases. To this end, let us de…ne the mixed u=s=p formulation, which uses the displacement …eld, u, together with the deviatoric component of the stresses, s, and the pressure …eld, p, as independent variables. Hence, the governing equations of the problem are rewritten as: r  s+rp+b=0(48) P:rsuDdev :s=0(49) r  uDvol p= 0 (50) where Eq.(48) is the balance of momentum equation in mixed form. The 5equations in (49) together with the scalar equation in (50) are the deviatoric and the volumetric counterparts of a generic constitutive equation as presented in (41 40). The weak form of the mixed u=s=p formulation reads: (v;r  s)+(v;rp)+(v;b)=0in (51) (;P:rsu);Ddev :s= 0 in (52) (q; r  u)Dvol (q; p)=0in (53) 7 where v(vector), q(scalar) and (a tensorial …eld of 5independent components) are the variations of the displacement, the pressure …eld and the deviatoric stresses, respectively, and (;)Ddenotes the integral of the product of two functions in a domain D, which is omitted when D. Integrating Eq.(51) by parts and taking v=0on @u, the problem can be written as: 0 + (rsv;s) + (r  v; p) = F(v) (;P:rsu);Ddev :s+ 0 = 0 (q; r  u) + 0 Dvol (q; p) = 0 (54) where F(v) = (v;b)+(v; t)@t(55) is the work of the external loads. Problem (54) involves the …rst derivatives of u(x). Hence, the natural space for the continuum displacements …eld is: u(x)2V=H1()ndim .Hm() denotes the space of functions whose derivatives (up to order m0) belong to L2(). The corresponding variations are de…ned in: v(x)2V0=fv(x)2Vjv=0for 8x2@ug. The pressure …eld, p, and its variation, q, belong to space Q=L2(), while the natural space for deviatoric part of the stress …eld, s, and its variations, , is b S=fs(x) = [sij (x)] ; sij =sji 2L2() jtr (s) = 0for x2. Other functional settings could be considered by changing the terms integrated by parts. In fact, the formulation that yields optimal stress convergence for equal interpolation for all the unknowns requires more regularity on the stresses. We will not treat this issue in this work (see [5] for similar ideas in the context of Darcy’s problem). Problem (54) is complemented by the Dirichlet boundary conditions in terms of the prescribed displacements, u. Remark 1 To achieve symmetry, it is possible to substitute the second term in the …rst equation, (rsv;s), by the equivalent term: (P:rsv;s). Then, problem (54) reads: 0 + (P:rsv;s) + (r  v; p) = F(v) (;P:rsu);Ddev :s+ 0 = 0 (q; r  u) + 0 Dvol (q; p) = 0 (56) Remark 2 Alternatively, the second equation in (56) can be tested using the test functions, = qI+with q2Qand 2b S, in the form: (;P:rsu);Ddev :s= 0 (57) This equation is equivalent to the original one because the extra terms, (qI;P:rsu)and qI;Ddev :s, are null by construction. However, Eq.(57) develops into 6scalar equations being only 5of them linearly independent. Eq.(57) admits as a solution a stress …eld, ^s, where tr (^s)6= 0 is not de…ned. 8 To overcome this inconvenient, it is necessary to prescribe the volumetric part of the stress …eld ^s. One possibility is to add the term ;Dvol :^sto Eq. (57), which corresponds to the weak form of the constraint: tr (^s) = 0. The result is the following equation: (;P:rsu)(;D:^s) = 0 (58) and problem (56) can be reformulated as: 0 + (P:rsv;^s) + (r  v; p) = F(v) (;P:rsu)(;D:^s) + 0 = 0 (q; r  u) + 0 Dvol (q; p) = 0 (59) where the stress …eld, ^s, is deviatoric in a weak sense only. 5 Discrete approximation for the three-…eld formulation 5.1 Galerkin approach Let us now de…ne the discrete Galerkin …nite element counterpart of problem (56) as: 0 + (P:rsvh;sh) + (r  vh; ph) = F(vh) (h;P:rsuh)h;Ddev :s+ 0 = 0 (qh;r  uh) + 0 Dvol (qh; ph) = 0 (60) where the discrete displacements, uh, and the corresponding test functions, vh, are de…ned in the …nite-dimensional subspaces VhVand V0;h V0, respectively. The approximate counterpart of the stress …eld, shand the pressure …eld, ph, together with their variations, hand qh, belong to the …nite element spaces b Shb Sand QhQ, respectively. Standard conforming approximations are considered. In the following, we will be interested in continuous …nite element spaces Vh,b Shand Qhand, more speci…cally, in equal interpolation for displacements, stresses and pressures. From the computational point of view, it is interesting to adopt Voigt’s notation, which transforms the tensorial format of a generic symmetric tensor into a 6dimensional vector. In the Cartesian system, the stress and strain tensors are expressed using Voigt’s notation as: s=sxx syy szz sxy sxz syz = [si](61) "="xx "yy "zz xy xz yz = ["i](62) where ij = 2"ij are the so-called engineering strains. Let nnod be the number of nodes per element of the …nite element partition. Then, Uh=UA i is an array of dimension 3nnod with 3components (lower subindex, i) of displacements at each node (upper subindex A) of the …nite element. Similarly, Sh=SA iis an array of dimension 6nnod 9 Tensor e C=e Cvol +Cdev (only used for stabilization purposes) is de…ned as: e C=e Cvol (II) + Cdev =8 < : CCompressible case e Cvol =Cvol =K 2GIIncompressible limit e Cvol =2 3G(106) This given, the solution of the problem is approximated as: u'uh+e u=uh+u[r  h+b](107) 'h+e =(1 )h+[C:rsuh]Compressible case (1 )h+[2Grsuh]Incompressible limit (108) and the corresponding stabilized problem for the compressible case is: (rsvh;C:rsuh) + (1 ) (rsvh;h) = F(vh) (1 ) (h;rsuh)(1 ) (h;D:h) + u(r  h;r  h)= 0 (109) The elemental sti¤ness matrix can be expressed as: K(e)=K(e) h(e) uK(e) u(e) K(e)  (110) which shows that there exist two di¤erent contributions adding stability to the original Galerkin’s problem. These contributions are: K(e) u=Z (e) [0] [0] [0]BBTd(111) K(e) =Z (e) BTCB BTN NT BNT DNd(112) In the incompressible limit, the stabilized problem reads: 2G(rsvh;rsuh) + (1 ) (rsvh;h) = F(vh) (1 ) (h;rsuh)(1 )h;Ddev :h+ u(r  h;r  h)= 0 (113) and the corresponding stabilization matrices are: K(e) u=Z (e) [0] [0] [0]BBTd(114) K(e) =Z (e) 2GBTB BTN NTBNTDdevNd(115) 16 Remark 13 The stabilization matrices for both the compressible and incompressible cases can be (formally) obtained summing the second and third rows and columns of the corresponding matrices de…ned for the three-…eld formulation, and assuming the same (discrete) interpolation functions for both the pressure and the deviatoric stress …elds as well as the same stabilization parameters: =s=p. Thus, the u=formulation can be considered a particular case of the u=s=p formulation. The latter allows one to approximate independently the deviatoric and volumetric parts of the stress, whereas in the former they are limited by the stress approximation chosen. 7 Numerical results In this section, the mixed three-…eld formulation presented in this work is tested in a number of numerical benchmarks. The objective is to show the performance of the proposed …nite element technology in terms of both displacement and stress …eld accuracy and its rate of convergence upon mesh re…nement. For the sake of brevity, the incompressible linear elasticity case is studied, although the results can be extended to more complex non-linear constitutive behaviors (allowing for the volumetric/deviator split). Full incompressibility (Poisson’s ratio: = 0:5) is assumed for the di¤erent numerical tests. The performance of the proposed method is compared with the behavior of the mixed displacement/pressure formulation (see [20]). The same displacement sub-grid scale stabilization (u= cuh2 2Gwith cu= 1:0) is adopted for both the u=p and u=s=p formulations, while the stress sub-grid has been introduced in the u=s=p formulation assuming: s=csh Lwith cs= 1 and p= 0. Calculations have been performed with an enhanced version of the …nite element code COMET (see [12]) developed by the authors at the International Center for Numerical Methods in Engineering (CIMNE). 7.1 Plane strain Cook’s membrane problem The Cook’s membrane problem is a bending dominated benchmark used by many authors to test their element formulations (see [34], [19], among others). The problem consists of a tapered panel, clamped on the left hand side and subjected to a shearing load (f= 1) at the right end. Initial geometry of this plane strain problem is shown in Figure 1. Material data have been assigned in terms of Young’s modulus: E= 200 and Poisson’s ratio: = 0:5. For the evaluation of the stabilization parameters in the mixed u=s=p formulation the characteristic length is taken equal to 50. In order to test the convergence behavior of the di¤erent formulations, the problem has been discretized into a series of regular mesh re…nements with 22,44,88,1616,3232,6464, 128 128 and 256 256 elements per side. Both structured quadrilateral and triangular meshes have been tested. In Figure 2 the 16 16 quadrilateral and triangular meshes are shown. Figure 3 compares the performance of the u=s=p vs. u=p formulations. It shows how they both converge to the reference values (best value obtained with the 256256 quadrilateral mesh using the u=s=p formulation) on mesh re…nement. The performance is monitored reporting the displacement of the top right corner (point Ain Figure 1) and both the J2 deviatoric stress and pressure value at the mid point of the bottom side of the Cook’s membrane (point Bin Figure 1). Note that the 17 Figure 1: Cook’s membrane problem: geometry. u=p formulation computes the stresses (locally) at the Gauss points (Figure 4 (a)) while the u=s=p formulation adopts a continuous stress …eld (Figure 4 (b)). To compare stress accuracy, a local smoothing technique has been applied to the original discontinuous stress …elds of the mixed u=p formulation. So, Figure 3 presents the continuous values obtained after the smoothing operation. Results in Figure 3 clearly show that both the u=s=p and the u=p formulations deal appropriately with the incompressibility constraint but the three-…eld formulation exhibits a higher convergence rate in the stress …elds. This translates into an enhanced stress accuracy of the proposed formulation even for very coarse meshes. The same considerations apply for both triangular and quadrilateral mesh discretizations. For the sake of completeness, Figure 3 also shows results obtained for the bilinear displacement/constant pressure Q1P0 element ([33]) and the bilinear displacement with enhanced strains Q1E4 element ([34]). These two elements can only approach the incompressible limit, so a value of Poisson’s ratio = 0:4999 has been used. The same local smoothing technique has been applied to the discontinuous stress …elds that these elements produce. Both quadrilaterals perform satisfactorily on these regular structured meshes, showing an accuracy comparable to that of the mixed u=s=p formulation. However. the …nite element formulations behind these two elements cannot be extended to triangular and tetrahedral meshes. 7.2 Cantilever beam Let us now consider a beam of unit thickness (plane-strain analysis) subjected to a bending moment imposed at right end side by means of a linear normal traction distribution (the maximum traction values is f= 10), as shown in Fig. 5. Both the horizontal and vertical displacements are prescribed at the bottom corner of the left hand side. The horizontal displacement is also …xed at the upper corner. The exact solution of this 18 (a) Structured quadrilateral mesh (b) Structured triangular mesh Figure 2: Cook’s membrane problem: (16 x 16) FE meshes. pure bending problem is given by: u(x; y) = 2f12 EL xL 2y(116) v(x; y) = f12 EL x2+ 1y(yL)(117) where u(x; y)and v(x; y)are the horizontal and vertical displacements (0 xLand 0yH). The beam length is L= 10 and its height is H= 2. Young’s modulus is set to E= 200 and Poisson’s ratio to = 0:5. For the evaluation of the stabilization parameters in the mixed u=s=p formulation the characteristic length is taken equal to H. To assess the accuracy of the proposed formulation, the vertical displacement at the top corner (A)of the right side of the beam is monitored (the analytical value is v(10;2) = 0:375), as well as the maximum horizontal normal stress and the mean stress (pressure) at the mid-point (B)of the bottom side (the analytical values are xx = 2 and p= 1, respectively). The computational domain has been discretized into a series of orthogonal re…ned meshes with 15,210,420,10 50 and 100 500 quadrilateral elements (see Fig. 6). The performance of the proposed three-…eld formulation is shown in Figure 7 and compared to that of the displacement/pressure formulation. The enhanced accuracy of the proposed formulation is clearly demonstrated. Using very coarse meshes such as the 210 mesh (only 2elements in the thickness), the error is less than 5% for the top corner displacement and less that 1% for both horizontal stress and pressure values. Figure 8 depicts the vertical displacement along the top side of the beam using di¤erent mesh resolutions for both the u=s=p and u=p formulations. The analytical solution is given by the parabolic function: v(x; 2) = f12 EL x2= 0:375 102x2(118) 19 (a) QUAD: displacement (b) TRI: displacement (c) QUAD: J2-stress (d) TRI: J2-stress (e) QUAD: pressure (f) TRI: pressure Figure 3: Cook’s membrane problem: performance of the proposed mixed three-…eld formulation compared with the mixed displacement/pressure formulation. 20 (a) u/p formulation (b) u/s/p formulation Figure 4: Cook’s membrane problem (16x16 triangular mesh): J2 deviatoric stress value. Figure 5: Cantilever beam: geometry. Once again, the behavior of the proposed formulation shows a great degree of accuracy even for very coarse meshes. Sensibility of the proposed formulation against element distortion is tested by comparing the above results on orthogonal meshes with those obtained on slanted meshes (see Fig. 6). Table 1 compares errors on maximum vertical displacements, horizontal normal stress and pressure obtained with the u=p and u=s=p formulations on the 10 50 orthogonal and slanted meshes. The accuracy of the proposed u=s=p formulation is almost insensitive to distortion. This is not the case for u=p formulation, where accuracy on stress and pressure deteriorates with slanting. 7.3 Sharp V-notched specimen under tension As a last example, the vertical stretching of a square V-notched specimen, shown in Figure 9, is considered. Dimensions of the sample are 22(width height) and the V-shaped notch has a 21 (a) Orthogonal mesh (b) Slanted mesh Figure 6: Cantilever beam: 10x50 meshes. Max vertical displ. Max horiz. stress Max pressure u=p orth. +10:6 % 0:30 % 8:95 % u=p slanted +2:03 % 13:0 % 12:43 % u=s=p orth. 0:26 % +0:55 % +3:14 % u=s=p slanted 0:53 % +0:35 % +1:74 % Table 1: Cantilever beam. Accuracy on orthogonal and slanted 10x50 meshes. 22 (a) Displacement (b) Max. hor. normal stress (c) Pressure Figure 7: Cantilever beam: performance of the proposed mixed three-…eld formulation compared with the mixed displacement/pressure formulation. 23 (a) u/p formulation (b) u/s/p formulation Figure 8: Cantilever beam problem: vertical displacement along the top side of the beam using di¤erent mesh resolutions. length of 1and a maximum width at the boundary of 0:02. For the evaluation of the stabilization parameters in the mixed formulation, L= 2 is taken as representative length of the problem. Uniform vertical displacements of opposed sign are imposed at the top and bottom boundaries. Figure 9: Geometry for the sharp V-notched specimen under tension In the continuous elastic problem associated to this situation, the strain and stress …elds are singular at the tip of the sharp notch. The discrete model corresponding to the u=p …nite element formulation performs satisfactorily in terms of a global error norm, but approximates very poorly the actual behavior near the singular points. To show this, a coarse structured mesh consisting of 882u=p triangles with a 45obias is 24 (a) u/p formulation: max. pr. stress contour-…ll. (b) u/p formulation: principal stresses. (c) u/s/p formulation: max. pr. stress contour-…ll. (d) u/s/p formulation: principal stresses. Figure 10: Sharp V-notch specimen under tension.Principal stresses. constructed. Figure 10(a) depicts contour-…ll of the maximum principal stresses computed on this mesh. Note the strong mesh bias dependence that is observed in front of and behind the notch tip. In fact, the largest values of the stresses occur behind the tip (left of the tip in the Figures), rather than in front of it (right of the tip in the Figures). Computed stress directions near the tip of the crack also show strong mesh bias dependence (see Figure 10(b)). These severe local errors caused by the mesh alignment are not alleviated by mesh re…nement. Figures 10(c)-(d) show corresponding results obtained using the stabilized mixed u=s=p formulation on the same mesh. The improved accuracy with respect to the u=p formulation is clear. In particular, the maximum principal stress value is detected exactly at the tip of the notch; computed stresses directions are also noticeably improved. The importance of these two features in nonlinear solid and ‡uid mechanics is evident. As it is shown in reference [18], they are crucial in strain 25