A stabilized mixed three-field formulation for stress accurate analysis including the incompressible limit in finite strain solid dynamics

A STABILIZED MIXED THREE-FIELD FORMULATION FOR STRESS ACCURATE ANALYSIS INCLUDING THE INCOMPRESSIBLE LIMIT IN FINITE STRAIN SOLID DYNAMICS INOCENCIO CASTAÑAR§, RAMON CODINA§,‡AND JOAN BAIGES§,‡ Abstract. In this work a new methodology for finite strain solid dynamics problems for stress accurate analysis including the incompressible limit is presented. In previous works, the authors have presented the stabilized mixed displacement/pressure formulation to deal with the incompressibility constraint in finite strain solid dynamics. To this end, the momentum equation is complemented with a constitutive law for the pressure which emerges from the deviatoric/volumetric decomposition of the strain energy function for any hyperelastic material model. The incompressible limit is attained automatically depending on the material bulk modulus. This work exploits the concept of mixed methods to formulate stable displacement/pressure/deviatoric stress finite elements. The final goal is to design a finite element technology able to tackle simultaneously problems which may involve incompressible behavior together with a high degree of accuracy of the stress field. The Variational Multi-Scale stabilization technique and, in particular, the Orthogonal Subgrid Scale method allows the use of equal-order interpolations. These stabilization procedures lead to discrete problems which are fully stable, free of volumetric locking, stress oscillations and pressure fluctuations. Numerical benchmarks show that the results obtained compare very favorably with those obtained with the corresponding stabilized mixed displacement/pressure formulation. Keywords: Incompressible hyperelasticity, Solid dynamics, Mixed interpolations, Stabilization methods, Orthogonal subgrid scales. 1. INTRODUCTION Several industrial manufacturing processes such as metal forming, forging, or friction stir welding among many others require, at the same time, stress accuracy and performance in the incompressible limit [1, 2]. It becomes crucial in these cases to use a finite element (FE) technology capable of dealing with complex phenomena such as strain localization [3], the formation of shear bands, the prediction of crack propagation [4] or the isochoric behavior of the inelastic strains [5]. Incompressibility is a widely accepted assumption used in continuum and computational mechanics [6]. In biomechanics, several materials can be modeled as nearly or fully incompressible [7]. Stress accuracy enhancement becomes very useful in many fields, such as cardiac electromechanics [8, 9] in which stress tensor acts as the coupling field with the equations describing electrical propagation in stress-assisted diffusion models [10]. Displacement-based low order FE methods perform poorly in such nearly and fully incompressible scenarios [11]. Volumetric and shear locking, pressure fluctuations and poor performance in bending dominated cases are some of the effects that are often found [12]. Popular solutions to tackle the nearly incompressible limit in the solid mechanics community are reduced and selective integration techniques [13], the B-bar and the F-bar Date: January 12, 2023. §Universitat Politècnica de Catalunya, Barcelona Tech, Jordi Girona 1-3, Edifici C1, 08034 Barcelona, Spain. ‡Centre Internacional de Mètodes Numèrics en Enginyeria (CIMNE), Edifici C1, Campus Nord UPC, Gran Capitán S/N, 08034 Barcelona, Spain. E-mails: [email protected] (IC),[email protected] (RC), [email protected] (JB). 1 I. CASTAÑAR, R. CODINA & J. BAIGES 2 methods [14] or the well-known mean dilatation FE method, which avoid these numerical instabilities by reducing the evaluation of the incompressibility constraints at quadrature points. However, these strategies are only designed to work with structured hexahedral meshes and they are not able to tackle the fully incompressible regime. Mixed formulations are well established and regularly used to avoid these instabilities. The use of different stabilization techniques, and particularly those based on the Variational Multi-Scale (VMS) framework [15], allows for the use of equal-order interpolations for all master fields. A particular formulation of this type, namely, the Orthogonal Subgrid Scale (OSGS) method [16], was used in [17] to design a stabilized FE formulation for the three-field linear Stokes problem, using displacements, pressure and deviatoric stresses as variables. The analysis and FE approximation of Darcy’s problem presented in [18] motivated the introduction of both strains/displacements and stresses/displacements pairs as primary variables in [19] for infinitesimal strain elasticity; in this particular case, one can change the functional framework to increase the accuracy in the calculation of the stresses. To tackle the incompressible limit, the pressure needs to be introduced as a variable [20, 21, 22], although it is also possible to design a formulation using the volumetric strain as unknown [23]. Formulations including stresses as unknowns produce a considerable increase in the number of unknowns per node, but they also increase the accuracy for strains and stresses. Furthermore, in [24] the idea of using a three-field displacement/pressure/deviatoric stress formulation was tested and seen to be very effective when solving incompressible cases in which accurate results for stress and strain fields are required. These FE technologies have demonstrated enhanced stress accuracy as well as the ability to capture stress concentrations and strain localizations guaranteeing stress convergence upon mesh refinement for first order elements. Mixed formulations are also applied under the transient finite strain assumption. In [25, 26] the velocity/pressure pair is taken as unknown of the problem and the displacement field is updated explicitly as a final step. The problem is stabilized with the VMS framework. A family of first-order form of the equations is presented in [27, 28, 29, 30, 31, 32] where the authors propose to use as primary variables the linear momentum p, the deformation gradient F, the cofactor tensor of the deformation gradient Hand the jacobian J; the objective for this choice of variables is to ease dealing with some complex constitutive laws, and in particular with polyconvex hyperelastic potentials. In [33] the incompressibility of the material is treated with the displacement/pressure pair in an updated Lagrangian formulation framework. Another possibility is to consider Finite Volume schemes to present a conservative cell-centered Lagrangian Finite Volume scheme for solving the hyperelasticity equations on unstructured multidimensional grids [34]. In previous works, the authors have applied stabilized mixed formulations for elasticity. Lately, a stabilized mixed displacement/pressure was presented in [35] for both nearly and fully incompressible hyperelastic material models. The system was stabilized by means of the VMS framework. The present work makes a step forward introducing a mixed three-field formulation based on displacement/pressure/deviatoric stress1elements with equal-order interpolations for all master fields. The only requirement is the introduction of the constitutive law for deviatoric stresses in the system of equations to be solved. This technology is expected to enhance stress accuracy as well as to increase the ability to capture stress concentrations with the guarantee of stress convergence upon mesh refinement. This work is organized as follows: In Section 2 the solid dynamics equations in finite strain theory are summarized and a novel mixed three-field formulation is developed. Furthermore, the variational form of the problem, its linearization and several employed time integrators are introduced. In Section 3 some VMS stabilization techniques are presented and the resulting stabilized forms of the three-field formulation are shown. In Section 4 1As it will be shown, the stress used as unknown is deviatoric in the deformed configuration only. I. CASTAÑAR, R. CODINA & J. BAIGES 3 several benchmarks and numerical examples are tested to assess the present formulation and to validate its performance. Also a comparison with its two-field formulation counterpart are highlighted. To end up, in Section 5 some conclusions of the proposed formulation are drawn. 2. SOLID DYNAMICS PROBLEM In this work we employ index notation to identify a vector or tensor with its Cartesian coordinates, either in the reference or the deformed configuration. As usual, repeated indexes imply summation for all space dimensions (see e.g. [6]). To denote scalar, vector and tensor quantities we use uppercase letters when they are evaluated in the reference configuration and lowercase letters if they are reckoned in the deformed one. We employ the index zero for quantities acting in the reference configuration. 2.1. Conservation equations. Let Ω0:= Ω (0) be an open, bounded and polyhedral domain of Rd, where d∈ {2,3}is the number of space dimensions. The initial configuration of the body is Ω0, whereas the current configuration of the body at time tis denoted by Ω (t). The motion is described by a function ψwhich links a material particle X∈Ω0to the spatial configuration x∈Ω (t)according to ψ: Ω0−→ Ω (t),x=ψ(X, t),∀X∈Ω0, t ≥0. The boundary of the reference configuration is denoted as Γ0:=∂Ω0and Γ (t):=∂Ω (t) represents the boundary of the current configuration at time t. We always assume that the mapping between both boundaries is defined through the motion, i.e., ψ(Γ0, t) = Γ (t). We denote as ]0, T [the time interval of analysis. The conservation of linear momentum in finite strain theory in a total Lagrangian formulation framework reads as (1) ρ0 ∂2ua ∂t2−∂ ∂XA {FaBSBA}=ρ0bain Ω0×]0,T[, where ρ0is the initial density, F=∂x ∂Xis the deformation gradient, Sis the second PiolaKirchhoff (PK2) stress tensor and ρ0bare the body forces. Mass conservation implies that (2) ρJ =ρ0, where ρis the density at time tand J=det F>0is the Jacobian of F. With regards to the balance of angular momentum, it implies that the PK2 stress tensor must be symmetric. The objective of this work is to obtain a mixed formulation for stress accurate analysis including the incompressible limit. The volumetric/deviatoric split of the Cauchy stress tensor σis the starting point to develop such formulation: (3) σ=σdev −pI, where σdev is the deviatoric part of σ,pis the pressure and Ithe second-order identity tensor. We can now use the relation between stresses to obtain a proper decomposition for the PK2 stress tensor [36]: (4) SAB =JF −1 Aa F−1 Bb σab (3) =JF −1 Aa F−1 Bb σdev ab −pJC−1 AB :=S0 AB −pJC−1 AB, where we have introduced the ‘deviatoric’ PK2 stresses S0(see Remark 2.1 below) and the right Cauchy-Green tensor C=FTF. Thanks to the decomposition in Eq. (4), the conservation of linear momentum can be reformulated as (5) ρ0 ∂2ua ∂t2−∂ ∂XAFaBS0 BA+∂ ∂XApJF−1 Aa =ρ0bain Ω0×]0,T[. I. CASTAÑAR, R. CODINA & J. BAIGES 4 Remark 2.1. Tensor S0is often referred to as the ‘true’ deviatoric component of S. The trace of σdev is zero by construction. However, it does not imply that the trace of S0also vanishes, and thus S0is not deviatoric in the algebraic sense. In fact, the ‘true’ deviatoric component of Ssatisfies the following equation (see for instance [36]): (6) S0:C= 0, which can be interpreted as the trace with respect to the metric tensor C. The above equation enables the hydrostatic pressure pto be evaluated directly from Sas (7) p=1 3JS:C. 2.2. Constitutive model. Let us restrict ourselves to nonlinear isotropic hyperelastic models (see [6, 11, 36] for further details). These models postulate the existence of a Helmholtz free-energy function (or strain energy function) Ψ. The PK2 stress tensor can be derived by taking derivatives of the Helmholtz free-energy functional with respect to the right Cauchy-Green tensor, namely (8) S= 2∂Ψ (C) ∂C. We want to deal with compressible models that can reach the incompressible limit case. To characterize such models, it is convenient to adopt a decoupled representation of the strain energy function of the specific form (9) Ψ (C) = W¯ C+U(J), where ¯ C=J−2/3Cis the volume-preserving part of C. Let us remark that this decomposition allows one to split the elastic response of the material into the so-called deviatoric and volumetric parts, respectively, measured in the initial configuration. We can now derive the PK2 stress tensor as (10) S= 2∂Ψ ∂C= 2∂W ∂C+ 2∂U ∂C= 2∂W ∂C+dU dJ JC−1. By comparing this definition with Eq. (4) we obtain expressions for both the pressure and the deviatoric PK2 stress tensor (11) S0= 2∂W ∂Cand p=−dU dJ . Several constitutive models for both deviatoric and volumetric components are shown in [35]. Let us just describe the ones which will be applied in this work. Readers are referred to [11, 25] for further details on this kind of models. 2.2.1. Deviatoric models. The strain energy density must be written in terms of the strain invariants, which are defined for the volume-preserving tensor ¯ Cby ¯ I1=trace ¯ C,¯ I2=1 2htrace ¯ C2−trace ¯ C2i.(12) Let us present two suitable functions for the deviatoric component of the strain energy function: •Neo-Hookean model This model results from considering only the first principal invariant: (13) W¯ I1=µ 2¯ I1−3, where µ > 0is the shear modulus. I. CASTAÑAR, R. CODINA & J. BAIGES 5 •Mooney-Rivlin model This model is derived considering the dependance on the second invariant as (14) W¯ I1,¯ I2=α1¯ I1−3+α2¯ I2−3, where α1and α2are material parameters that must satisfy µ= 2 (α1+α2)>0. 2.2.2. Volumetric models. Due to the decoupled form of the strain energy density, compressibility is accounted for by the volumetric strain energy function. Let us now show two models that depend upon the bulk modulus κ=2µ(1+ν) 3(1−2ν), where νis the Poisson ratio. •Quadratic model [37]: (15) U(J) = κ 2(J−1)2;dU dJ =κ(J−1) . •Simo-Taylor model [38]: (16) U(J) = κ 4J2−1−2 log J;dU dJ =κ 2J−1 J. Remark 2.2. The volumetric functions can be written as U(J) = κG(J). Therefore, Eq. (11) can be used to obtain a proper way to impose the incompressibility of an hyperelastic material (17) p=−dU dJ ⇔p=−κdG dJ ⇔p κ+dG dJ = 0. This equation can be applied regardless of the compressibility of the material under study. It is interesting to observe that in the incompressible limit, when Poisson’s ratio ν→0.5 (for isotropic materials) then κ→ ∞ and Eq. (17) reduces automatically to (18) dG dJ = 0. Eq. (18) imposes directly that J= 1, which is in fact the condition that a material must satisfy to be incompressible in finite strain theory. Remark 2.3. It is interesting to show how to impose incompressibility if the real deviatoric/volumetric decomposition of the PK2 stress tensor is considered. The following relation holds: S=Sdev −p∗I=S0−pJC−1, where p∗=1 3traceS. If we take the trace of the PK2 stress tensor, we can obtain an expression for the pseudo-pressure p∗as p∗=−1 3S0 AA +κdG dJ JC−1 AA, which allows us to write the volumetric component of the constitutive equation as 1 κ 3p∗+S0 AA JC−1 AA +dG dJ = 0. Taking into account the widely used decomposition of the strain energy function given by Eq. (9), it seems more natural and effective to consider the classical decomposition, which gives us simpler equations in nearly incompressible scenarios. I. CASTAÑAR, R. CODINA & J. BAIGES 6 2.3. Governing equations. In this section, a novel three-field formulation is introduced. The objective is the definition of a general framework, which includes the mixed two-field formulation presented in [35] to be able to tackle the incompressible limit and introduces S0as primary unknown to obtain a higher accuracy in the computation of stresses in finite strain problems. To this end, let us introduce the three-field mixed upS0formulation. Let D={(X, t)|X∈Ω0,0< t < T}be the space-time domain where the problem is defined. The problem consists of finding a displacement field, u:D−→ Rd, together with a deviatoric component of the PK2 stress tensor, S0:D−→ Rd⊗Rdand a pressure field, p:D−→ Rsuch that ρ0 ∂2ua ∂t2−∂ ∂XAFaBS0 BA+∂ ∂XApJF−1 Aa =ρ0bain Ω0×]0,T[,(19) p κ+dG dJ = 0 in Ω0×]0,T[,(20) S0 AB −2∂W ∂CAB = 0 in Ω0×]0,T[.(21) The governing equations must be supplied with initial conditions of the form u=u0, ∂u ∂t =v0in Ω0at t= 0, with u0and v0given, and a set of boundary conditions which can be split into Dirichlet boundary conditions (22), where the displacement is prescribed, or Neumann boundary conditions (23), where the value of tractions tNare prescribed, i.e.: u=uDon Γ0,D,(22) n0·(F·S)(4) =n0·(F·S0)−pJn0·F−T=tNon Γ0,N ,(23) where n0is the geometric unit outward normal vector on the boundary of the reference configuration Γ0. Remark 2.4. Verifying that this formulation reduces to the three-field mixed formulation in linear elasticty when infinitesimal strains theory is considered is crucial. Regarding to the momentum equation (19), the following simplifications are obtained: FaBS0 BA (4) ≈Jσ0 ab ≈1 + ∂uc ∂xcσ0 ab =σ0 ab + : ≈0 (∇ · u)σ0 ab ≈σ0 ab (24) pJF−1≈p(1 + ∇ · u)I=pI+ : ≈0 p∇ · uI ≈pI(25) and taking into account that both the reference and current configurations match, we obtain the simplified momentum equation for linear elasticity (26) ρ∂2ua ∂t2−∂σ0 ab ∂xa +∂p ∂xa =ρbain Ω0×]0,T[. With respect to the incompressibility equation (20), as it was previously mentioned, dG dJ is a function which imposes J= 1. Therefore, (27) p κ+dG dJ ≈p κ+J−1≈p κ+ (1 + ∇ · u)−1 = p κ+∇ · u= 0 which is exactly the constitutive law of the pressure when considering linear elasticity. Let us recall that in the incompressible limit κ→ ∞, and this equation will reduce automatically to (28) ∇ · u= 0 which is the incompressibility condition for infinitesimal strain theory. Finally, with regards to the deviatoric constitutive law (21) we can observe that S0=JF−1σdevF−T≈(1 + ∇ · u) (I+∇u)−1σdev (I+∇u)−T =σdev +Okuk2≈σdev (29) I. CASTAÑAR, R. CODINA & J. BAIGES 7 and all the deviatoric models of the strain energy presented in this work must satisfy that, in the infinitesimal strain assumption, they recover the deviatoric constitutive law of linear elasticity when second order terms are neglected. Therefore (30) σdev = 2∂W ∂C≈Cdev :ε, where εis the infinitesimal strain field and Cdev is the 4th order deviatoric constitutive tensor for isotropic linear elastic materials and it is defined as (31) Cdev = 2µI−1 3I⊗I, where Iis the 4th order identity tensor. With Eqs. (26, 27 and 30) we recover automatically the three-field formulation for linear elasticity presented in [24]. 2.4. Variational form of the problem. Given a set ω⊂Ω0, we shall use the symbol h·,·iωto refer to the integral of the product of two functions, assuming it well-defined. The subscript is omitted when ω= Ω0. Let V,Qand Tbe, respectively, the proper functional spaces where displacement, pressure and deviatoric PK2 stress solutions are well-defined for each fixed time t∈]0, T [. We denote by V0functions in Vwhich vanish in the Dirichlet boundary Γ0,D. We shall be interested also in the spaces W:=V×Q×Tand W0:=V0×Q×T. The variational statement of the problem is derived by testing system (19)-(21) against arbitrary test functions, V:= [v, q, T]T,v∈V0,q∈Qand T∈T. The weak form of the problem reads: find U:= [u, p, S0]T: ]0, T[→Wsuch that initial conditions and the Dirichlet condition (22) are satisfied and (32) va, ρ0 ∂2ua ∂t2+A(U,V) = F(V)∀V∈W0, where A(U,V)is a semilinear form defined on W×W0as A(U,V):=∂va ∂XA , FaBS0 BA−∂va ∂XA , pJF −1 Aa +q, dG dJ  +Dq, p κE+TAB, S0 AB−TAB,2∂W ∂CAB .(33) In addition, F(V)is a linear form defined on W0as (34) F(V):=hva, ρ0bai+hva, tNaiΓ0,N . Integration by parts has been used in order to decrease the continuity requirements of unknowns pand S0. 2.5. Time discretization. In this work implicit time integrators are considered. When the material under study is either nearly or fully incompressible, the Courant-FriedrichsLewy condition, which involves the bulk modulus, becomes very restrictive. If explicit time integrators are considered, extremely small time steps are needed in order to satisfy it. Furthermore, in the fully incompressible case, as κ→ ∞, solving the problem with explicit time integration is not possible. Although in principle any implicit time discretization method can be applied, it has to be taken into account that a hyperbolic system of equations of second order in time is being solved . It is important to control spurious high-frequency oscillations that might appear in the solution for both nearly and fully incompressible hyperelastic materials. As a consequence, numerical time integrators with high-frequency dissipation will be applied. Let us now consider a partition of the time interval [0, T ]into Ntime steps of size δt, assumed to be constant. I. CASTAÑAR, R. CODINA & J. BAIGES 8 2.5.1. Backward differentiation formula (BDF). Given a generic time dependent function at a time step tn+1 =tn+δt, for n= 0,1,2, ... the approximation of the time derivative of order k= 1,2, ... is written using information from already computed time instants. In our problem, we have to approximate the second time derivative of the displacement, ∂2u ∂t2 n+1 :=an+1. Depending on the accuracy of the method, we can select the specific formulae: BDF1 : an+1 =1 δt2[un+1 −2un+un−1] + O(δt), BDF2 : an+1 =1 δt2[2un+1 −5un+ 4un−1−un−2] + Oδt2. 2.5.2. Newmark-βequations. This is a popular class of time integrators [11]. In this time integration formula, the updated acceleration an+1 and velocity vn+1 are given by: an+1 ≈1 βδt2[un+1 −un−δtvn−δt2 2(1 −2β)an], vn+1 ≈vn+ (1 −γ)δtan+γδtan+1. Here βand γare parameters to be tuned. When β=1 4and γ=1 2the Newmark-βmethod is implicit, unconditionally stable and second-order accurate for linear problems. Remark 2.5. Newmark-βmethod is not unconditionally stable for nonlinear problems. The algorithmic energy conservation was recognized as the key to achieve stability in long term simulations according to the generalized theorem presented in [39]. The stability in time is further studied in [40, 41]. In [42, 43] several energy-momentum consistent timestepping schemes are proposed for some mixed formulations in nonlinear problems. 2.6. Linearization. In order to solve the problem, the system needs to be linearized so that a bilinear operator which allows to compute a correction δUof a given guess for the solution at time tn+1 is obtained, that we denote by Un+1. Iteration counters will be omitted to simplify the notation. After using a Newton-Raphson scheme, we obtain the following linearized form of the problem. Given Un+1 as the solution at time tn+1 and the previous iteration, find a correction δU:= [δu, δp, δS0]T∈W0such that (35) v, ρ0 C δt2δu+B(δU,V) = F(V)− A Un+1,V−v, ρ0an+1∀V∈W0, where B(δU,V)is the bilinear form obtained through the Newton-Raphson linearization and it is defined on W0×W0as B(δU,V) = ∂va ∂XA ,∂δua ∂XB S0 BA+∂va ∂XA , FaBδS0 BA−∂va ∂XA , JpF −1 Bb ∂δub ∂XB F−1 Aa  +∂va ∂XA , JpF −1 Ab ∂δub ∂XB F−1 Ba −∂va ∂XA , JδpF −1 Aa +q, f (J)F−1 Aa ∂δua ∂XA +q, δp κ−TAB,C0 ABCDFaC ∂δua ∂XD+TAB, δS0 AB,(36) where f (J)is a function coming from the linearization of dG dJ and depends upon the volumetric strain energy function into consideration and C0 ABCD = 2 ∂2W ∂CAB∂CCD is the deviatoric constitutive tangent matrix which relates variations of the deviatoric PK2 stress tensor, δS0, with variations of the Right Cauchy tensor, δC. Note that for every implicit time integrator presented here, we can write: ∂2u ∂t2tn+1 ≈C δt2δu+an+1, I. CASTAÑAR, R. CODINA & J. BAIGES 9 where Cis a coefficient depending on the time integration scheme and an+1 is the acceleration computed at the previous iteration, which in the time discretized problem will be given by any of the expressions introduced above. 2.7. Symmetrization. In the way we have written the problem, it is not symmetric. To achieve symmetry2, it is possible to modify Eq. (35) by (37) v, ρ0 C δt2δu+Bmod (δU,V) = F(V)−Amod Un+1,V−v, ρ0an+1∀V∈W0, where Bmod (δU,V)is the bilinear form defined on W0×W0as Bmod (δU,V) = ∂va ∂XA ,∂δua ∂XB S0 BA+∂va ∂XA , FaBδS0 BA−∂va ∂XA , JpF −1 Bb ∂δub ∂XB F−1 Aa  +∂va ∂XA , JpF −1 Ab ∂δub ∂XB F−1 Ba −∂va ∂XA , JδpF −1 Aa +q, JF−1 Aa ∂δua ∂XA +q, J f (J) δp κ−TAB, FaA ∂δua ∂XB+TAB,C0−1 ABCDδSCD,(38) and Amod (U,V)is a semilinear form defined on W×W0as Amod (U,V):=∂va ∂XA , FaBS0 BA−∂va ∂XA , pJF −1 Aa +q, J f (J) dG dJ  +q, J f (J) p κ+TAB,C0−1 ABCDS0 CD−TAB,C0−1 ABCD2∂W ∂CCD ,(39) where we have multiplied the second equation by the linearized term J f(J)and we have introduced C0−1as the inverse deviatoric constitutive tangent matrix. To define such 4th order tensor, it is necessary to obtain the inverse strain energy function, which involves a nonlinear problem. In spite of this difficulty, the symmetric form of the problem can be interesting from both the theoretical and the practical point of views. From the theoretical point of view, the problem to be solved corresponds to the minimization of a certain mechanical energy, whereas from the practical point of view the symmetry can be exploited when solving the linear system. Also from the conceptual standpoint, the test functions for the constitutive equation in the non-symmetric case are in fact strains, whereas in the symmetric case they are stresses. Let us also remark that, while in the infinitesimal strain theory it is equivalent to use stresses or strains as unknowns, in finite strain elasticity the symmetrization of the problem using strains (e.g. the Green-Lagrange strain E) is much more involved than using the PK2 stress that we have presented. In any case, we will not discuss here the introduction of strains as unknowns of the problem. For simplicity, we will employ the non-symmetric form of the problem in what follows, although the use of the symmetric version would be straightforward. 2.8. Galerkin spatial discretization. The standard Galerkin approximation of this abstract variational problem is now straightforward. Let Phdenote a FE partition of the domain Ω0. The diameter of an element domain K∈ Phis denoted by hKand the diameter on the FE partition by h=max{hK|K∈ Ph}. We can now construct conforming FE spaces Vh⊂V,Qh⊂Q,Th⊂Tand Wh=Vh×Qh×Thin the usual manner, as well as the corresponding subspaces Vh,0⊂V0and Wh,0=Vh,0×Qh×Th,Vh,0being made of functions that vanish on the Dirichlet boundary. In principle, functions in Vhare continuous, whereas functions in both Qhand Thnot necessarily. 2In the case of homogeneous boundary conditions, obviously. I. CASTAÑAR, R. CODINA & J. BAIGES 16 For the sake of completeness, Figs. 1b-1d-1f show the nonlinear iteration convergence error for each unknown of the formulation. As it can be seen, a quadratic convergence is attained thanks to the Newton-Raphson linearization of the problem. 10 -3 10 -2 10 -1 10 -8 10 -6 10 -4 10 -2 10 0 (a) Displacement error upon mesh refinement 10 -3 10 -2 10 -1 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 (b) Pressure error upon mesh refinement 10 -3 10 -2 10 -1 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 (c) Deviatoric stress error upon mesh refinement 10 110 210 310 410 510 610 7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 (d) Deviatoric stress error upon number of degrees of freedom 10 010 110 210 310 410 5 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 (e) Deviatoric stress error upon CPU time Figure 2. Manufactured convergence test. Comparison of convergence between the upS0formulation and the upformulation I. CASTAÑAR, R. CODINA & J. BAIGES 17 Mesh Number of elements 16×6×36 ×6 28×8×48 ×6 312 ×12 ×72 ×6 Table 1. Bending Beam. Different mesh and number of elements. More interesting results are obtained when comparing these convergence rates with respect to the ones obtained with the mixed upformulation. Figs 2a-2b show the displacement and pressure convergence rates upon mesh refinement, respectively. Both fields are considered as primary unknowns of both formulations and therefore a similar accuracy for a given mesh size can be expected. Fig. 2c displays the deviatoric PK2 stress convergence rates upon mesh refinement for both formulations. As expected, higher accuracy is achieved for a given mesh size for the mixed upS0formulation. To achieve the same accuracy, e.g. 0.1%of global error, the upformulation requires a mesh size, h, almost 10 times finer (h≈0.003) than the upS0formulation (h≈0.03) as it can be seen in Fig. 2c. Furthermore, to obtain a fairer comparison, the same study is conducted in terms of the number of degrees of freedom (DOFs) in Fig. 2d. To get an error lower than 0.1%, the upformulation requires 6·104DOFs approximately, while the upS0formulation needs less than 2·103DOFs (25 times lesser than the upformulation). Results clearly show that both the upS0and the upformulations deal appropriately with the incompressibility constraint but the three-field formulation exhibits a higher accuracy in the stress field, even for very coarse meshes. For the sake of exhaustiveness, Fig. 2e depicts the total CPU time needed by each FE technology to achieve a given global deviatoric stress accuracy. In particular, to reduce the simulation error below 0.1%, the upS0formulation is more or less 10 times faster compared to the upone. Remark 4.1. Note that the upformulation computes the stresses (locally) at the numerical integration points, while the upS0formulation adopts a continuous stress field. To compare stress accuracy, a local smoothing technique has been applied to the original discontinuous stress fields of the mixed two-field formulation. So, Figs. 2c-2d-2e present the continuous values obtained after the smoothing operation. 4.2. Bending beam. As a second test in finite strain elasticity, we consider a three dimensional beam of square section clamped on its bottom face very similar to the one presented in [25, 35]. The initial geometry is a thick column of dimensions 1×1×6 m as shown in Fig. 3a. We consider stress free conditions in all boundaries except the clamped one in which zero displacement is imposed. An initial linear in space velocity field v0(X, Y, Z) = 5 3Z, 0,0Tm/s is imposed so as to start the column oscillations in time, leading to a large oscillatory bending deformation. A Mooney-Rivlin material with initial density ρ0= 1.1×103kg/m3and material parameters α1= 2.69 MPa and α2= 0.142 MPa is considered. In order to avoid unphysical modes appearing from the time integration scheme, we have selected the mildly-dissipative BDF2 time integrator with time step δt = 0.01 s. The main goal of this example is to show that our stabilized mixed formulation works properly in bending dominated scenarios and in 3D cases. For this reason, we have selected 3 different structured linear tetrahedral meshes (as the one shown in Fig. 3b), specified in Table 1. Let us start showing some interesting properties about the three-field mixed upS0formulation presented here with respect to the classical displacement-based formulation [36] I. CASTAÑAR, R. CODINA & J. BAIGES 18 (a) Geometry (b) Mesh Figure 3. Bending Beam. Geometry (3a) and tetrahedral structured mesh (3b). (from now on named as uformulation). To do so, let us consider the bending beam problem for several compressible regimes. We consider 3 differents scenarios: ν= 0.2, which reproduces a compressible material, ν= 0.49, which mimics a nearly incompressible material and finally, we take ν= 0.5to reach the incompressible limit. All these cases are performed with Mesh 2. Fig. 4a displays the evolution in time up to T= 3 for the first component of the displacement field at point A. As expected, in the compressible regime (ν= 0.2), both formulations exhibit very similar results. More interesting conclusions can be drawn when moving to the nearly incompressible regime (ν= 0.49). In such case, the displacement-based formulation presents volumetric locking, which tends to show smaller displacements than the expected ones. On the contrary, the upS0formulation is able to obtain proper solutions without presenting these instabilities. Furthermore, in the incompressible limit (ν= 0.5), the upS0formulation gives us precise solutions whereas the u formulation fails over the running stage. To end up this study, Figs. 4b-4c-4d show the pressure field and some components of the deviatoric PK2 stress tensor run with the upS0 formulation. As it can be clearly seen, well-defined solutions are obtained regardless the incompressibility of the material and no oscillations can be appreciated even for this coarse mesh. From now on, let us consider a fully incompressible material with ν= 0.5. Fig. 5 presents the evolution of point A along time up to T= 3 s for both upand upS0formulations. Figs. 5a-5b show the L2(Ω0)norm for the displacement field and the pressure field, respectively. As previously commented, both unknowns are considered as main variables of the problem. Very similar results can be observed for the displacement field when comparing both formulations for a fixed mesh. Despite the fact that the pressure field is a master field for both formulations, it turns out that more accurate results are obtained for the upS0formulation for a fixed mesh due to the capability of the method to capture stress concentrations better than the upformulation. It is interesting to remark, that the observed behavior seems to be more dissipative in the three-field formulation. This indicates that including the deviatoric PK2 stress tensor as unknown of the problem in the upS0formulation both enhances the accuracy of the solution and its energy dissipation rate. On the contrary the upformulation shows a less optimal dissipative behavior for the same mesh. Furthermore, Fig. 5c presents the evolution for the L2(Ω0)norm for I. CASTAÑAR, R. CODINA & J. BAIGES 19 (a) X-Displacement component (m) (b) Pressure (Pa) (c) XX-Deviatoric PK2 stress component (Pa) (d) XY-Deviatoric PK2 stress component (Pa) Figure 4. Bending beam. Comparison between uand upS0formulations while increasing the incompressibility of the material at point A. the deviatoric PK2 stress field. As expected, for a fixed mesh, more accurate results are obtained when introducing the deviatoric PK2 stresses S0as an extra unknown of the problem in the three-field formulation. For the sake of thoroughness, we show in Figs. 6-7 the deformed beam at t= 2.25 s and at t= 3 s, respectively run with Mesh 1. First of all, we can observe that very similar deformations and pressure fields are appreciated for both formulations. Finally, regarding the deviatoric PK2 stress tensor, one can see the gain of accuracy in this field for the upS0formulation by including this field as primary unknown of the problem instead of computing it from the displacement derivatives. 4.3. Twisting column. As a final example, we present the twisting column test. This test is widely used to assess the robustness of any formulation in extreme nonlinear deformations [25, 26, 35, 49, 50]. The initial geometry of the column is the same as the one shown in Fig. 3a. We consider stress free conditions and zero displacement initial conditions are applied on the corresponding boundaries. In order to make the column twist, we apply an initial sinusoidal velocity field: (66) v0(X, Y, Z) = ωsin πZ 12 (Y, −X, 0)Tm/s where ω= 100 rad/s. The material is considered to be Neo-Hookean with initial density ρ0= 1.1×103kg/m3, shear modulus µ= 5.7×106Pa and Poisson ratio ν= 0.5, to model a fully incompressible material. To define the deviatoric part of the material, we consider a Simo-Taylor law. Several levels of refinement have been considered to perform the computations. To construct the meshes, we select structured hexahedral elements. So we consider 3 different meshes. Mesh 1, with 6×6×36 trilinear FEs, Mesh 2 with 16 ×16 ×96 FEs and we end up with Mesh 3 with 32 ×32 ×192 FEs. We select a time step δt = 0.002 s. I. CASTAÑAR, R. CODINA & J. BAIGES 20 (a) L2(Ω0)norm displacement (m) (b) Pressure (Pa) (c) L2(Ω0)norm deviatoric PK2 stress (Pa) Figure 5. Bending beam. Evolution of point A along time for both up and upS0formulations. I. CASTAÑAR, R. CODINA & J. BAIGES 21 (a) upS0(b) up(c) upS0(d) up Figure 6. Bending beam. Comparison between upS0and upat t= 2.25 s. Pressure field (Pa) and L2(Ω0)norm of the deviatoric PK2 stress (Pa) (a) upS0(b) up(c) upS0(d) up Figure 7. Bending beam. Comparison between upS0and upat t= 3 s. Pressure field (Pa) and L2(Ω0)norm of the deviatoric PK2 stress (Pa) First of all, let us perform some analysis for the different time integration schemes presented in Section 2. We run the same problem with different schemes and Mesh 2 and the main results can be seen in Fig. 8. On the one hand, left figures display the main unknowns up to T= 0.5s. As it can be seen, both BDF schemes are capable of reproducing the whole event. However, BDF1 is only first-order accurate in time and it is highly dissipative, excessively mitigating physical oscillations. With regards to the BDF2 scheme, it is able to dissipate the nonphysical modes, which helps preventing high frequency oscillations while keeping the second-order accuracy of the method. On the other hand, right figures illustrate the evolution obtained with a Newmark scheme for β=1 4and γ=1 2, which results in a second-order scheme in time. This method does not introduce any numerical dissipation, and therefore, it does not eliminate high frequency nonphysical oscillations. Next, let us fix BDF2 as time integration scheme and perform some comparisons between the upand the upS0formulations. Fig. 9 shows the evolution of point A up to T= 0.5 s. In Fig. 9a we can observe a comparison for the displacement field. As expected, both formulations show very similar results, which become closer upon mesh refinement. Moving to the pressure field in Fig. 9b, one can see that for a fixed mesh, similar evolutions are obtained but the upS0formulation gives more accurate results taking into account the evolution when refining the mesh. More interesting remarks can be made for the deviatoric I. CASTAÑAR, R. CODINA & J. BAIGES 22 BDF1 BDF2 (a) Displacement BDF1 and BDF2 Newmark (b) Displacement Newmark BDF1 BDF2 (c) Pressure BDF1 and BDF2 Newmark (d) Pressure Newmark BDF1 BDF2 (e) Deviatoric PK2 stress BDF1 and BDF2 Newmark (f) Deviatoric PK2 stress Newmark Figure 8. Twisting column. Time integrators comparison. PK2 stress tensor in Fig 9c. As it can be clearly appreciated, for a fixed mesh the threefield formulation attains more accurate results than the two-field version of the problem. In fact, we can observe that we need always an extra level of refinement for the upformulation to be able to achieve the same accuracy as the one given by the upS0formulation. To complete this example, Figs. 10-11 display the evolution of the deformation for the twisting column at different stages of the problem with Mesh 2. As it can be observed, the problem is well-captured, and no numerical oscillations can be seen neither for the pressure field nor for the deviatoric PK2 stress tensor. Let us remark, once more, the capability of the formulation to capture stress concentrations, in this case, placed at the clamped face of the twisting column. 5. CONCLUSIONS In this paper we have described a new stabilized FE method for stress accurate analysis in solid dynamics when considering nearly and fully incompressible materials. The point of departure is the splitting of the Cauchy stress tensor into deviatoric and spherical components, which then translates into a splitting of the second Piola–Kirchhoff stress tensor. The momentum equation is complemented with a constitutive law for the pressure which emerges from the deviatoric/volumetric decomposition of the strain energy function for I. CASTAÑAR, R. CODINA & J. BAIGES 23 (a) L2(Ω0)norm displacement (m) (b) Pressure (c) L2(Ω0)norm deviatoric PK2 stress (Pa) Figure 9. Twisting column. Evolution at point A. I. CASTAÑAR, R. CODINA & J. BAIGES 24 (a) 0.1s(b) 0.2s(c) 0.3s(d) 0.4s(e) 0.5s Figure 10. Twisting column. Deformation and Pressure field (Pa) along time. (a) 0.1s(b) 0.2s(c) 0.3s(d) 0.4s(e) 0.5s Figure 11. Twisting column. Deformation and L2(Ω0)norm deviatoric PK2 stress (Pa) along time. any hyperelastic material model. This law is formulated properly to obtain a simple way to impose the incompressibility of the material automatically. Furthermore, to design a FE technology with a high degree of accuracy of the stress field, the constitutive law for deviatoric stresses is added to the system to obtain a monolithic system of equations for the displacement/pressure/deviatoric stress formulation. The presented three-field approach is able to deal with any hyperelastic material, including fully incompressible cases. We have proposed two residual-based (ASGS and OSGS) and a term-by-term (S-OSGS) type stabilization techniques based on the decomposition of the unknowns into FE scales and SGSs. All stabilization techniques are able to circumvent the compatibility restrictions on the interpolation functions among the primary unknowns of the problem. Furthermore, I. CASTAÑAR, R. CODINA & J. BAIGES 25 the proposed scheme shows the desired rate of convergence upon mesh refinement regardless the stabilization technique. It is interesting to remark that the S-OSGS stabilization technique allows us to obtain well-defined solutions and to neglect terms that do not contribute to stability. This methods turns out to be more robust for solving problems when large stress gradients are present. Likewise, for the examples we have considered, we have been able to assume quasi-static SGSs, although dynamic SGSs might need to be considered if very small time step sizes are required. Concerning the computational cost of the method, we have observed that the proposed methods display quadratic non-linear convergence regardless the stabilization technique, as it is expected from the implementation of a Newton–Raphson iterative procedure. The proposed three-field formulation is convergent upon mesh refinement, virtually free of any volumetric or shear locking. The technology is suitable for engineering applications in which a higher accuracy of stresses is needed. A comparison with the two-field formulation (displacement/pressure) is also carried out. Results clearly show that both the upS0and the upformulations deal appropriately with the incompressibility constraint but the threefield formulation exhibits a higher accuracy in the stress field, even for very coarse meshes. Acknowledgements Inocencio Castañar gratefully acknowledges the support received from the Agència de Gestió d’Ajut i de Recerca through the predoctoral FI grant 2019-FI-B-00649. R. Codina gratefully acknowledges the support received through the ICREA Acadèmia Research Program of the Catalan Government. This work was partially funded through the TOP-FSI: RTI2018-098276-B-I00 project of the Spanish Government. References [1] C.A. Moreira, G.B. Barbat, M. Cervera, and M. Chiumenti. Accurate thermal-induced structural failure analysis under incompressible conditions. Engineering Structures, 261:114213, 2022. [2] N. Dialami, M. Chiumenti, and M. Cervera. Defect formation and material flow in friction stir welding. European Journal of Mechanics A / Solids, 80:103912, 2020. [3] M. Cervera, M. Chiumenti, and R. Codina. Mixed stabilized finite element methods in nonlinear solid mechanics. Part II: Strain localization. Computer Methods in Applied Mechanics and Engineering, 199(37–40):2571–2589, 2010. [4] M. Cervera, G.B. Barbat, M. Chiumenti, and J.Y. Wu. A comparative review of xfem, mixed fem and phase-field models for quasi-brittle cracking. Archives of Computational Methods in Engineering, 29:1009–1083, 2022. [5] M. Chiumenti, M. Cervera, C.A. Moreira, and G.B. Barbat. Stress, strain and dissipation accurate 3-field formulation for inelastic isochoric deformation. Finite Elements in Analysis and Design, 192:103534, 2021. [6] G. A. Holzapfel. Nonlinear solid mechanics. A continuum approach for engineering. Wiley, 2000. [7] G. A. Holzapfel. Biomechanics of soft tissue. The Handbook of Materials Behavior, 3:1049–1063, 2001. [8] P. E. Farrell, L.F. Gatica, B. P. Lamichhane, R. Oyarzúa, and R. Ruiz-Baier. Mixed kirchhoff stressdisplacement-pressure formulations for incompressible hyperelasticity. Computer Methods in Applied Mechanics and Engineering, 374:113562, 2021. [9] R. Ruiz-Baier, A. Gizzi, A. Loppini, C. Cherubini, and S. Filippi. Thermo-electric effects in an anisotropic active-strain electromechanical model. Communications in Computational Physics, 27(1):87–115, 2020. [10] A. Propp, A. Gizzi, F. Levrero-Florencio, and R. Ruiz-Baier. An orthotropic electro-viscoelastic model for the heart with stress-assisted diffusion. Biomechanics and Modeling in Mechanobiology, 19(2):633– 659, 2020. [11] T. Belytschko, W. K. Liu, and B. Moran. Nonlinear Finite Elements for Continua and Structures. Wiley, 2001. [12] T. J. R. Hughes. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Prentice-Hall,Englewood Cliffs, New Jersey, 1987. [13] D. S. Malkus and T. J. R. Hughes. Mixed finite element methods - Reduced and selective integration techniques: A unification of concepts. Computer Methods in Applied Mechanics and Engineering, 15:63–81, 1978.