Hyperelastic models are the widely used constitutive equations for rubber-like polymers and soft tissues, and these models are the building blocks for the constitutive relations of fast-growing soft and smart materials such as electroactive and magnetoactive polymers. The deformation behaviour of materials modelled with hyperelastic constitutive equations ranges from the compressible regime in which materials undergo significant volumetric changes to the incompressible regime where the volumetric change is zero.
Hyperelastic models in the incompressible regime require sophisticated finite element formulations for computing accurate numerical results. In the literature, there exist various techniques for addressing the issue of incompressibility. For the truly incompressible case, the mixed displacement-pressure formulation is the only computationally-efficient option, which, however, results in saddle-point problems book-fem-BrezziFortin ; book-fem-ZienkiewiczVol2 ; KadapaPhDThesis ; KadapaCMAME2016elast ; KadapaIJNME2019mixed . To overcome the numerical issues in solving saddle-point problems, the widely followed approach in computational solid mechanics is to impose the incompressibility constraint weakly using quasi-incompressible approximations. Some popular approaches for simulating quasi-incompressible hyperelastic materials are -bar formulation NetoIJSS1996 ; NetoIJNME2005Fbarpatch , enhanced-strain method SimoJAM1986 , average nodal strain formulation BonetCNME1998 ; PiresCNME2004 , reduced integration method with hourglass control FlanaganIJNME1981 ; BelytschkoCMAME1984hglass , enhanced assumed strain (EAS) methods WriggersCM1996 ; KorelcCM2010 , the two-field displacement-pressure formulation ChenCMAME1997 ; KadapaPhDThesis ; KadapaCMAME2016elast , the three-field formulation SimoCMAME1985 , mixed stabilised formulations ChiumentiCMAME2002 ; CerveraCMAME2010a ; CerveraCMAME2010b ; ScovazziIJNME2016 ; ScovazziCMAME2017velocity ; AbboudIJNME2018 ; FrancaNM1988 ; KlassCMAME1999 ; MasudJAM2005 , energy-sampling stabilisation PakravanIJNME2017b ; BijalwanIJMSI2019 , least-squares formulations KadapaIJNME2015 ; KadapaPhDThesis ; ManteuffelSIAMJNA2006 ; MajidiSIAMJNA2001 ; MajidiSIAMJNA2002 , -bar projection methods ElguedjCMAME2008 . Some noteworthy recent contributions to incompressible and quasi-incompressible computational solid mechanics are in LeiEC2016 ; MehnertMMS2017 ; WulfinghoffCMAME2017 ; WriggersCM2017 ; JiangIJCM2018 ; BayatAMSES2018 ; BayatCM2018 ; CoombsCMAME2018 ; SevillaIJNME2018 ; ToghipourCMA2018 ; BijalwanIJMSI2019 ; MoutsanidisCPM2019 ; DalIJNME2019 ; OnishiIJNME2017 ; OnishiIJCM2019 ; SevillaCS2019 ; connolly2018 ; connolly2019 ; ViebahnAMSES2019 ; RajagopalMAMS2018 .
While the two-field mixed displacement-pressure formulation based on the perturbed Lagrangian approach, cf. ChenCMAME1997 ; KadapaPhDThesis ; KadapaCMAME2016elast ; book-fem-BonetWood ; SchroderCM2017 , yields accurate results for the nearly and truly incompressible models, the accuracy of results deteriorates for problems that undergo significant stretches in the compressible regime, i.e., for Poisson’s ratio, . This deterioration of accuracy is due to the absence of the relation between the pressure and the volumetric energy function in the perturbed Lagrangian approach. Although such an issue can be overcome by employing the weak statement for the two-field mixed formulation KadapaIJNME2019mixed ; AbboudIJNME2018 ; ScovazziIJNME2016 , this approach results in unsymmetric matrix systems for the majority of volumetric energy functions. Note that any solver for unsymmetric matrices is expensive when compared to the one for symmetric matrices. Although a mixed displacement-pressure formulation using complementary energy functions yields symmetric matrix systems, the difficulties associated with computing the complementary energy functions limit the applicability of such a formulation to a few simple volumetric energy functions OrtigosaCMAME2016 .
The three-field mixed displacement-pressure-Jacobian formulation is one of the widely-used numerical approaches in recently emerged computational electromechanics and magnetomechanics for rubber-like materials BisharaMMS2018 ; ParkIJSS2012 ; ParkSS2013 ; SeifiCMAME2018 ; AskIJNLM2012 ; PelteretIJNME2016 ; JabareenPIUTAM2015 ; mehnert2018 ; mehnert2019 . The formulation takes into account the relation between the pressure and the volumetric energy functional and also results in symmetric matrix systems. However, since the mechanical part of the constitutive models for rubber-like polymers is decomposed into deviatoric and volumetric parts, the three-field formulation gives no particular computational advantage for such materials modelling because of the fact that the coupling term between the displacement and Jacobian variables vanishes. Besides, the three-field formulation is not applicable for simulating the truly incompressible materials.
Recently, Schröder et al. SchroderCM2017 proposed a novel consistent two-field mixed displacement-pressure formulation that is independent of the penalty function. This approach, however, is limited to incompressible and quasi-incompressible cases as it generalises the classical perturbed Lagrangian approach for different penalty functions in the incompressibility limit. Motivated by the work of Schröder et al. SchroderCM2017 , in this paper, we present a new generalised mixed displacement-pressure formulation that (i) is applicable for compressible, quasi-incompressible as well as truly incompressible cases, (ii) takes into account the relation between the pressure and the volumetric energy function, (iii) yields symmetric matrix systems irrespective of the volumetric energy function, and (iv) proves to be an efficient alternative for the three-field formulation. Such a consistent formulation is necessary for simplified implementations and for viscoelastic and elastoplastic material models in which the deformation behaviour can vary between compressible and incompressible regimes during the course of a simulation.
This paper is organised as follows. The basics of finite strain elasticity are introduced in Section 2. The mixed displacement-pressure formulation based on the perturbed Lagrangian approach and the issued associated with it are discussed in Section 3, followed by the displacement-pressure formulation based on the weak Galerkin statement in Section 4. The proposed mixed displacement-pressure formulation is derived in Section 5, and is compared against the well-established three-field mixed displacement-pressure-Jacobian formulation in Section 6. The accuracy of the proposed formulation is demonstrated with two numerical examples in Section 7. The paper is concluded with Section. 8 with a summary of observations made and conclusions drawn.
2 Governing equations for finite strain elasticity
2.1 Deformation, strain and stress
With as the reference configuration of the solid body under consideration, and as its new configuration under the action of external forces, the nonlinear deformation map that takes a point to a point is denoted as . Now, the displacement field () and the deformation gradient () are defined as
is the second-order identity tensor. Two important strain measures that are used in this work are the right Cauchy-Green deformation tensor () and the determinant of deformation gradient (), and they are defined as,
For truly incompressible materials, the deformation of the solid is such that the total volume change at any point in the domain is zero. This can be represented mathematically as the incompressibility constraint in the finite strain regime as,
For modelling hyperelastic materials in the incompressible finite strain regime, the deformation gradient, , is decomposed into deviatoric and volumetric components as
Using the above definitions, the modified versions of the deformation gradient and the right Cauchy-Green tensor are related as
In this work, the strain energy function for the hyperelastic materials is assumed to be decomposed into a deviatoric part, , and a volumetric part, , as
Here, the deviatoric part of the function is often represented as a function of the invariants of and the volumetric part is only a function of the Jacobian . For comprehensive details on the widely-used strain energy functions for hyperelastic materials, we refer the reader to the works of Steinmann et al. SteinmannAAM2012 , Hossain and Steinmann HossainJMBM2013 , Hossain et al. HossainJMBM2015 , Marckmann and Verron MarckmannRCT2006 , Doll and Schweizorhof DollJAM2000 and Moerman et al. Moerman2019 , and references cited therein.
Now, two important stress measures of interest in the present work are defined as
2.2 Equilibrium equations
In the absence of dynamic effects, the balance of linear momentum and the corresponding boundary conditions in the reference configuration are given as
where, is the body force per unit undeformed volume, is the unit outward normal on the boundary , is the prescribed value of displacement on the Dirichlet boundary and is the specified traction force per unit undeformed area on the Neumann boundary . Here, and are such that and .
where, is deviatoric component of the first Piola Kirchhoff stress which is computed solely from the deviatoric part of the strain energy function, and is an independent approximation for the hydrostatic pressure. For the compressible and the nearly incompressible cases, the pressure is related to the volumetric energy function via the relation
3 Perturbed Lagrangian formulation
For the truly incompressible materials, the volumetric energy function vanishes. Therefore, the incompressibility constraint given by equation (4
) needs to be imposed either strongly by using the Lagrange multiplier approach or weakly by using the penalty approach. The penalty approach results in a formulation that consists of displacement degrees of freedom (DOFs) only; therefore, such an approach suffers from the well-known issues of volumetric locking and spurious pressure fields observed with the pure displacement formulation. Hence, for the truly incompressible case, the Lagrange multiplier approach is the only viable choice. In this approach, the Lagrange multiplier becomes an independent variable for the pressure field.
With as the Lagrangian multiplier, the potential energy functional is given as
where, is the potential energy corresponding to the applied body and surface loads. It is given as
To overcome the numerical difficulties in solving the saddle-point problems resulting from the Lagrangian multiplier approach, the term corresponding to the constraint equation is modified slightly, and such an approach is known as the perturbed Lagrangian approach ChenCMAME1997 ; book-fem-BonetWood ; KadapaCMAME2016elast ; KadapaIJNME2019mixed ; KadapaPhDThesis ; Bercovier ; Wriggers ; Zienkiewicz ; Simo ; Tur . The energy functional for the perturbed Lagrangian approach is given by
where, is the bulk modulus. Note that such a modification introduces slight compressibility into the model, thus, facilitating a numerical framework for modelling nearly incompressible cases, i.e., for or equivalently, . For such cases, the difference in the values of computed from different volumetric energy functions, as shown in Fig. 1, can be considered to be negligible. Note also that the energy functional for the perturbed Lagrangian approach (19) reduces to that of the Lagrange multiplier approach (17) for the truly incompressible case, i.e., for for which .
However, despite its suitability for both nearly and truly incompressible materials modelling, the perturbed Lagrangian approach fails to represent the physical behaviour of the material in the compressible regime. This is due to the lack of accountability for the relation between the pressure and the volumetric energy function, i.e., , in the perturbed Lagrangian approach (19). Because of the significant differences in the pressure values computed using different volumetric energy functions in the compressible regime, as illustrated in Fig. 2 for , the perturbed Lagrangian approach produces incorrect results for problems simulated with compressible material models that experience severe stretches. To overcome this issue, we propose a new generalised mixed displacement-pressure formulation that is applicable for compressible, nearly incompressible as well as truly incompressible cases.
4 Weak Galerkin formulation
The mixed displacement-pressure Galerkin formulation for the finite strain hyperelasticity in the original configuration can be stated as: find the displacement and pressure such that for all the test functions and ,
where, and are the function spaces for the displacement and pressure, and and are the corresponding test function spaces. is the generic constraint that is chosen depending upon whether the material is truly or nearly incompressible. is defined as
For the finite element analysis, the displacement and pressure fields, and , and their test functions, and , are approximated as
are the vectors of displacement degrees of freedom (DOFs), andand are the vectors of the pressure degrees of freedom; and are the matrices of basis functions, respectively, for the displacement and pressure fields. Using the above definitions, the gradient and the divergence of displacement () and its test function () can be written as,
where, , and are the components of displacement, respectively, in X, Y, and Z directions, and and are the gradient-displacement and the divergence-displacement matrices, respectively. and for a single basis function (), are given as
Analogous to the gradient of displacement in Eq. (4), the Cauchy stress tensor is represented in the vector form as
Remark I: The formulation can also be written in terms of the symmetric part of the strain and stress measures by following the Voigt notation, see KadapaPhDThesis ; KadapaIJNME2019bbar ; KadapaIJNME2019mixed ; KadapaCMAME2016elast ; book-fem-ZienkiewiczVol2 . Such an approach, however, leads to complications in transforming the fourth- and third-order tensors resulting from the electrical, magnetic and coupled energy functions in electro- and magnetomechanics problems. For the ease of computer implementation of the finite element discretization, we present the two-field formulations in terms of full strain and stress tensors as vectors. The proposed approach avoids ambiguities in dealing with the derivatives of symmetric and unsymmetric second-order tensors.
where, is the effective Cauchy stress, and it is defined as
Assuming the displacement and pressure DOFs at the th iteration of the Newton-Raphson scheme are and , respectively, the corresponding DOFs at the th iteration become
where and are the incremental displacement and pressure DOFs, respectively. The corresponding increments for the continuum variables are denoted as and .
For the application of Newton-Raphson scheme, the linearisation of the constraint variable yields,
where, the variable is defined as
Following equation (21), we obtain,
and for :
Now, the resulting coupled matrix system for the mixed formulation can be written as
In Eqn (42), is the elasticity tensor of order four, and in the index notation, it is given as,
For the truly incompressible case and for , the coupled matrix system (41) is symmetric; for the other cases, it is unsymmetric. Thus, requiring a solver for unsymmetric sparse matrix systems which is expensive when compared to that of symmetric sparse matrix systems, see Duff et al. book-Matrices-Duff , Gould et al. GouldACMTOMS2007 , Dongarra et al. book-linalgebra-Dongarra , Pissanetzky book-Matrices-Pissanetzky , Saad book-Matrices-Saad , Barrett et al. book-Matrices-Barrett , Nachtigal et al. NachtigalSIAMJMA1992 and Watkins book-Matrices-Watkins , and references therein. This is another motivation behind the novel mixed formulation discussed in the following section.
5 Proposed mixed formulation
The definition of energy functional similar to the ones for the Lagrangian multiplier (17) and the perturbed Lagrangian (19) approaches is not straightforward because of the nonlinear nature of the relation (16) between the pressure and the volumetric energy function. Although we can devise the following energy functional for the mixed formulation,
such a functional results in an unconventional form of the effective Cauchy stress as
Therefore, for the variational approach, the energy functional that yields the conventional form of the Cauchy stress (33) is sought. The energy functional for the proposed mixed formulation should be of the form
with being the corresponding energy function that takes the relation (16) into account. We obtain by linearising . With the linearisation of at the ()th load step as
where . The linearised form of can be written as
Since and the corresponding quantities are evaluated from the displacement from the previously converged load step, and are constants for the current load step. Therefore, there is no need for their linearisations when applying the Newton-Raphson scheme.
Now, the first variation and the associated energy functional for the relation (16) can be written as
Note that the mixed formulation for the truly incompressible case can be recovered from the above equations by simply setting and . Moreover, by setting and , the energy functional for the perturbed Lagrangian formulation (19) can also be recovered. Thus, the proposed formulation not only takes into account the relation between the pressure and the volumetric energy function for the compressible materials but also is applicable for the truly incompressible case.
The total energy functional for the proposed formulation can now be written as
The first () and the second () variations for the energy functional become
Using following identities,
and considering the approximations for the solution variables and their variations as
the coupled matrix system can be written as
It is evident from the above equations that the coupled matrix system (67) of the proposed two-field formulation is always symmetric irrespective of the volumetric energy function, unlike the coupled system (41) of the mixed Galerkin approach which is symmetric only for .
Remark II: A mixed displacement-pressure formulation that yields symmetric matrix system irrespective of the volumetric energy function can also be derived by using complementary potentials to volumetric energy functions, as discussed in Appendix A. However, such an approach requires apriori knowledge of complementary functions, which are uncommon in the constitutive modelling for hyperelastic materials. Alternative techniques based on numerical methods for the computation of complementary potentials and their derivatives pose additional difficulties for the cases with multiple roots, see Appendix A for a detailed discussion. The disadvantages of such a formulation based on complementary energy functions are also motivating factors for the proposed formulation.
Remark III: Note that the Lagrange multiplier () in the displacement-pressure (hybrid) formulation in ABAQUS abaqus-manual is not the same as the independent approximation for the hydrostatic pressure and it is chosen intrusively to get the symmetric expression for the rate of virtual work (equation 3.2.3-11 in ABAQUS’s theory manual abaqus-manual ). This makes ABAQUS’s displacement-pressure formulation much more difficult to implement than the proposed formulation. Moreover, ABAQUS’s formulation ignores changes in the volume for the evaluation of the second variation (, in equation 3.2.3-13 in ABAQUS’s theory manual) which has the potential to deteriorate rate of convergence of Newton-Raphson iterations, while the proposed formulation poses no such issues since it uses consistent linearisation. Furthermore, ABAQUS uses an artificial numerical constant, (refer to equation 3.2.3-4 in ABAQUS’s theory manual), to overcome solver difficulties while the proposed formulation is free from such ad-hoc parameters.
6 Comparison with the three-field formulation
In this section, a comparison between the proposed two-field displacement-pressure mixed formulation and the well-established three-field displacement-pressure-Jacobian formulation SimoCMAME1985 ; book-fem-ZienkiewiczVol2 is presented. For this purpose, we consider the widely-used Q1/P0 finite element discretisation, in which it is possible to condense the pressure DOFs out. Other discretisations, for example, P2/P0 book-fem-BrezziFortin or the recently proposed BT2/BT0 KadapaIJNME2019mixed , are equally applicable.
6.1 Proposed two-field formulation : A condensed form
For the Q1/P0 element, the coupled matrix of the proposed two-field formulation (67) can be condensed as
Using , and simplify to
Since the pressure is assumed to be discontinuous, the effective stiffness matrix, and the effective residual can be computed element-wise before assembling into the global stiffness matrix.
6.2 Three-field mixed formulation
For the three-field mixed displacement-pressure-Jacobian formulation, only those quantities that are relevant to the discussion in the present work are presented here. For comprehensive details on the three-field formulation, we refer the reader to the works of Zienkiewicz and Taylor book-fem-ZienkiewiczVol2 or Kadapa KadapaPhDThesis .
where, represents the change in the volume, and by taking approximations for variables , and as,
the resulting coupled matrix system for the three-field formulation with the energy functional given by Eq. (A.1) can be written as,
For the strain energy functions decomposed into deviatoric and volumetric parts, it can be shown analytically that
See Chapter 4 in Kadapa KadapaPhDThesis for the comprehensive details on the derivation of equations (96). Using the above equations and invoking the constant discontinuous approximation for the Jacobian variable and pressure, we obtain the element-wise stiffness matrices as,