DeepAI

# Four-Field Mixed Finite Element Methods for Incompressible Nonlinear Elasticity

We introduce conformal mixed finite element methods for 2D and 3D incompressible nonlinear elasticity in terms of displacement, displacement gradient, the first Piola-Kirchhoff stress tensor, and pressure, where finite elements for the curl and the div operators are used to discretize strain and stress, respectively. These choices of elements follow from the strain compatibility and the momentum balance law. Some inf-sup conditions are derived to study the stability of methods. By considering 96 choices of simplicial finite elements of degree less than or equal to 2 in 2D and 3D, we conclude that 28 choices in 2D and 6 choices in 3D satisfy these inf-sup conditions. The performance of stable finite element choices are numerically studied. Although the proposed methods are computationally more expensive than the standard two-field methods for incompressible elasticity, they are potentially useful for accurate approximations of strain and stress as they are independently computed in the solution process.

10/20/2019

### A Conformal Three-Field Formulation for Nonlinear Elasticity: From Differential Complexes to Mixed Finite Element Methods

We introduce a new class of mixed finite element methods for 2D and 3D c...
06/23/2019

### Compatible-Strain Mixed Finite Element Methods for 3D Compressible and Incompressible Nonlinear Elasticity

A new family of mixed finite element methods-compatible-strain mixed fin...
07/12/2020

### A computationally tractable framework for nonlinear dynamic multiscale modeling of membrane fabric

A general-purpose computational homogenization framework is proposed for...
05/09/2018

### Spatial Poisson Processes for Fatigue Crack Initiation

In this work we propose a stochastic model for estimating the occurrence...
09/08/2020

### Three-field mixed finite element methods for nonlinear elasticity

In this paper, we extend the tangential-displacement normal-normal-stres...
08/31/2021

### A mixed method for 3D nonlinear elasticity using finite element exterior calculus

This article discusses a mixed FE technique for 3D nonlinear elasticity ...
06/04/2021

### Finite element stress analysis of a combined stacker-reclaimer machine: A design audit report

Design audit or design verification is an important step in engineering ...

## 1 Introduction

It is well-known that incompressible nonlinear elasticity is useful to describe the mechanical behavior of various soft tissues [1]. Deriving stable numerical methods for modeling such tissues that may undergo large deformations is a challenging task, for example see [2] and references therein. In this paper, we develop four-field mixed finite element methods for incompressible nonlinear elasticity by extending three-field mixed finite element methods introduced in [3] for compressible nonlinear elasticity. These four-field methods are based on a mixed formulation with displacement, displacement gradient, the first Piola-Kirchhoff stress tensor, and pressure as the independent unknowns. A potential advantage of this mixed formulation is providing accurate approximations of strain and stress, since they are independent unknowns and are explicitly calculated in the solution process, that is, no post-processing is required to indirectly compute strain and stress from displacement.

Similar to [3], finite elements suitable for the and the operators are respectively employed to discretize displacement gradient and stress in the proposed conformal mixed finite element methods. One can show that these choices follow from the strain compatibility and the balance law [4, 5]. To discretize displacement and pressure, - and -conformal finite element spaces are used, respectively.

A stability analysis based on the abstract theory of [6, 7] for Galerkin approximations of nonlinear problems is provided. More specifically, an inf-sup condition is written which is locally sufficient for the convergence of solutions of nonlinear finite element methods to a regular solution of incompressible nonlinear elasticity. This inf-sup condition can be interpreted as the necessary and sufficient condition for the well-posedness of the linear problems associated to Newtons’ iterations for solving nonlinear finite element methods. Based on this interpretation, we also write other inf-sup conditions that are necessary for the stability of Newtons’ iterations. Using these conditions, we derive some relations between the dimension of finite element spaces for different unknowns, which are necessary for the stability of Newtons’ iterations.

By considering choices of simplicial finite elements of degree less than or equal to for obtaining mixed finite element methods, we conclude that choices in D and choices in D do not satisfy all these inf-sup conditions, in general, and lead to unstable methods. The performance of stable choices in D and D are studied by solving numerical examples.

This paper is organized as follows: The four-field mixed formulation and the associated conformal finite element methods are discussed in Sections 2 and 3. Newtons’ iterations for solving nonlinear finite element methods are provided in Section 4. In Section 5, different inf-sup conditions are derived for studying the stability of finite element methods and it is shown that some choices of finite elements lead to unstable methods as they violate the inf-sup conditions. In Section 6, the inf-sup conditions are numerically studied. Moreover, by solving two numerical examples, the performance of stable finite element methods are investigated in D and D. Finally, some concluding remarks are given in Section 7.

## 2 A Mixed Formulation for Incompressible Elasticity

We write a four-field mixed formulation for incompressible nonlinear elasticity in this section. For simplicity, we consider static problems here, however, this formulation can be readily extended to time-dependent problems as well. Let denote the reference configuration of a D or D elastic body and let

be the unit outward normal vector field associated to the boundary

of .

Suppose are respectively displacement, displacement gradient, and the first Piola-Kirchhoff stress tensor. The constitutive relation of incompressible elastic bodies can be expressed as , where is a given function of , the real-valued function is pressure, and is the identity tensor [8, Section ]. The boundary value problem of incompressible nonlinear elastostatics can be stated as:

Given a body force , a displacement of , and a traction vector field on , find such that
(2.1) (2.2)

To write a week form based on the above equations, we consider the following spaces: Let be the standard space of square-integrable functions and let be the Sobolev space of -functions with first-order derivatives of -class. The space of vector fields in , , with components of -class is denoted by . The space of -vector fields that vanish on is written as .

We also consider the spaces and , which are the spaces of second-order tensor fields with -components that their and are respectively of -class. Consider the matrix representing the Cartesian components of a second-order tensor field in . Since and are respectively obtained by applying the standard and of vector fields to each row of , the spaces and can be considered as -copies of the corresponding spaces for vector fields.

For obtaining a concise form for the weak formulation, we employ the following notation: The standard inner product of is denoted by “” and the -inner products of real-valued functions, vector fields, and tensor fields are denoted by , that is, , , and , where we consider the summation convention on repeated indices. By taking the -inner product of the equations (2.1) with suitable test functions and using Green’s formula, one can write the following four-field mixed formulation for incompressible nonlinear elasticity:

Given a body force , a boundary displacement , and a surface traction vector field on , find such that , on , and
(2.3)

The above weak form is an extension of the formulation of [3] for compressible elasticity to incompressible elasticity. The choices of the spaces and guarantee that physically relevant jump conditions, namely, the Hadamard jump condition for strain and the continuity of traction vector fields at interfaces for stress, will be hold in the weak sense and also on the discrete-level. Following the discussion of [3, Remark 2], it is not hard to show that even for hyperelastic materials, the above weak formulation is not associated to a stationary point of any functional.

## 3 Conformal Finite Element Methods

Suppose , , , and are finite element spaces such that , , , and . Also let and let

be the interpolation operator associated to

. Then, one can write the following conformal mixed finite element method for (2.3):

Given a body force , a boundary displacement , and a surface traction vector field on , find such that , on , and
(3.1)

We employ simplicial finite elements shown in Figure 1 to generate the above finite element spaces in D and D. More specifically, we use the standard simplicial Lagrange elements and discontinuous elements respectively for and . As mentioned earlier, the definitions of and for second-order tensors imply that in , the tensorial spaces and are equivalent to -copies of the associated vectorial spaces in the sense that each row of second-order tensors of the - and -classes can be considered as a vector field beloning to the vectorial spaces and corresponding to the standard and operators of vector fields, respectively. We use the Nédélec elements of the first and the second kinds to generate and the Raviart-Thomas elements and the Brezzi-Douglas-Marini elements to generate .

## 4 Newtons’ Iterations

We employ Newton’s method to solve (3.1). Here we mention this approach in some details since, as will be discussed in the following section, it is also useful for stability analysis. Let and . Also let and , with , , and . The nonlinear problem (2.3) can be written as: Find such that

The -th Newton’s iteration for solving the above nonlinear equation reads: Given , find by solving a linear problem and let . This linear problem is obtained by linearizing the above equation and can be stated as: Find such that

 b(z,y;ui)=−⟨H(ui),y⟩,∀y∈Z, (4.2)

where the bilinear form is given by

where is the elasticity tensor in terms of the displacement gradient and . By discretizing the linear problem (4.2) using the finite element spaces of the previous section, one obtains the following Newton’s iteration for the finite element method (3.1):

Let , and . Given , find and let , where is obtained by solving the linear finite element method

Suppose

 (4.5)

and let , denote the total number of degrees of freedom. Then, it is straightforward to see that the stiffness matrix of the above linear problem is of the form

 Snt×nt=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣0\tabularcell@hbox0\tabularcell@hboxS1dn1×nd0[2pt/2pt]Sc1nc×n1\tabularcell@hboxSccnc×nc\tabularcell@hbox00[2pt/2pt]0\tabularcell@hboxSdcnd×nc\tabularcell@hboxSddnd×ndSdDnd×nD[2pt/2pt]0\tabularcell@hboxSDcnD×nc\tabularcell@hbox00⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦. (4.6)

## 5 Stability Analysis

The stability of the nonlinear finite element method (3.1) can be studied by using the general theory for the Galerkin approximation of nonlinear problems discussed in [6, 7]. Let be a regular solution of (2.3), or equivalently (4.1), in the sense that the derivative of the nonlinear mapping defined in (4.1) is nonsingular at . This derivative can be expressed in terms of the bilinear form introduced in (4.3). The theory of [6, 7] states that in a neighborhood of and for sufficiently refined meshes with the maximum element diameter , under some additional mild conditions, the nonlinear finite element method (3.1) has a unique solution that converges to as if there exists a mesh-independent number such that

 infyh∈Zhsupzh∈Zhb(zh,yh;u)∥zh∥Z∥yh∥Z≥α>0, (5.1)

where is the mixed finite element space for the Newton’s iteration (4.4) and is the norm of . Thus, the inf-sup condition (5.1) is a sufficient stability condition.

The condition (5.1) can be interpreted as a stability condition for the linear finite element method of the Newton’s iteration (4.4) as well. More specifically, suppose that at the -th Newton’s iteration, we have . Then, the condition (5.1) implies that the linear problem (4.4) at the the -th iteration has a unique solution and the stiffness matrix introduced in (4.6) is non-singular. One can also approximate the lower bound of (5.1) by using a matrix associated to [9, Proposition 3.4.5]. In particular,

can be approximated by the smallest singular value of the matrix

, where is the symmetric, positive definite matrix associated to the norm .

The non-singularity of the stiffness matrix implies that the submatrices

 S1dn1×nd,SDcnD×nc, B=⎡⎢⎣0\tabularcell@hbox0\tabularcell@hboxS1dn1×nd[2pt/2pt]Sc1nc×n1\tabularcell@hboxSccnc×nc\tabularcell@hbox0⎤⎥⎦, C=⎡⎢⎣0\tabularcell@hboxS1dn1×nd\tabularcell@hbox0[2pt/2pt]Sdcnd×nc\tabularcell@hboxSddnd×nd\tabularcell@hboxSdDnd×nD⎤⎥⎦, (5.2) (5.3)

must be full rank. This result can be stated by using some inf-sup conditions associated to suitable bilinear forms induced by (4.4). The upshot can be stated as follows:

• is full rank if and only if there exists such that

• is full rank if and only if there exists such that

 infqh∈VDhsupMh∈Vch\llangledet(I+Kih)tr[(I+Kih)−1Mh],qh\rrangle∥Mh∥∥qh∥≥αh; (5.5)
• is full rank if and only if there exists such that

 (5.6)

where

 (5.7)
• is full rank if and only if there exists such that

 inf(Υh,πh)∈[V1h,1]n×Vdhsup(Mh,Qh,rh)∈Vch×Vdh×VDhbC((Mh,Qh,rh),(Υh,πh);Kih,pih)∥(Mh,Qh,rh)∥∥(Υh,πh)∥≥αh, (5.8)

where

• is full rank if and only if there exists such that

 inf(λh,qh)∈Vch×VDhsup(Vh,Mh)∈[V1h,1]n×VchbD((Vh,Mh),(λh,qh);Kih)∥(Vh,Mh)∥∥(λh,qh)∥≥αh, (5.11)

where

• is full rank if and only if there exists such that

 inf(πh,qh)∈Vdh×VDhsup(Mh,Qh,rh)∈Vch×Vdh×VDhbE((Mh,Qh,rh),(πh,qh);Kih,pih)∥(Mh,Qh,rh)∥∥(πh,qh)∥≥αh, (5.14)

where

 bE((Mh,Qh,rh), (5.15) +\llanglepih(I+Kih)−TMTh(I+Kih)−T,πh\rrangle−\llanglerh(I+Kih)−T,πh\rrangle (5.16) +\llangledet(I+Kih)tr[(I+Kih)−1Mh],qh\rrangle. (5.17)

A simple inspection of the number of rows and columns of the above matrices together with the rank-nullity theorem yield the following necessary condition for the validity of (5.1) and the well-posedness of (4.4), where we use the notation of (4.5).

###### Theorem 1.

The inf-sup condition (5.1) does not hold and the Newton’s iteration (4.4) does not admit a unique solution if at least one of the following inequalities holds: (i) ; (ii) ; (iii) ; (iv) .

By using this result, one can specify some unstable combinations of finite elements. For example, consider the D elements of Figure 1 and suppose , , and are respectively the number of vertices, edges, and elements of a D simplicial mesh. One can write the relations

 Nel−Ned+Nv=1−I, and 2Ned−N∂ed=3Nel, (5.18)

where is the number of holes and is the number boundary edges [10]. By using the notation of Figure 1, we obtain the following relations on a D simplicial mesh:

 nD={3Nel, D1 element,6Nel, D2 element, n1={2Nv, L1 % element,2(Nv+Ned), L2 element, (5.19) nc=2Ned, N11 element, nd=2Ned, R1 element. (5.20)

A simple varification of the inequalities of Theorem 1 by using the relations (5.18) and (5.19) yields the following result.

###### Corollary 2.

Let , , , and be respectively arbitrary finite element spaces for , , , and spaces. In D, the choices , , , and for the finite element methods (4.4) violate at least one of the inequalities of Theorem 1 and therefore, lead to unstable finite element methods.

## 6 Numerical Results

Some numerical examples are discussed in this section to study the stability and the convergence of the mixed finite element methods (3.1). We consider the simplicial finite elements of Figure 1 and use FEniCS [11] to implement these numerical examples. The underlying domains are assumed to be a unit square in D and a unit cube in D with simplicial meshes shown in Figure 2. We consider incompressible Neo-Hookean materials with the stored energy function , , and the constitutive equation

 P(F)=μF−pF−T. (6.1)

The term of the bilinear form (4.3) then simply reads .

#### Notation.

To refer to any choice of the elements of Figure 1 for the discretization of (2.3), we use the abbreviated names shown in Figure 1. For example, indicates the choice of the Lagrange element of degree , the degree Nédélec element of the first kind, the Raviart-Thomas element of degree , and the discontinuous element of degree respectively for the discretization of .

### 6.1 Stability Study

The inf-sup conditions introduced earlier hold only if the corresponding matrices are full rank. Therefore, a simple approach for computationally studying these inf-sup conditions is to compute the rank deficiency of the matrices , , , , , . Some choices of finite elements that violate these inf-sup conditions are shown in Figure 3. In this figure, the full-rankness of the above matrices are plotted for several choices of finite elements by using D and D meshes of Figure 2, where the full-rankness of a matrix is the rank of divided by its maximum possible rank. Thus, , if is full rank and , otherwise.

The results of Figure 3 are computed using the incompressible Neo-Hookean constitutive equation (6.1) with near the reference configuration, that is, the matrices are associated to the first Newton’s iteration (4.4) starting at the reference configuration. The rank deficiency of some finite elements such as for , for , for , and for , is a simple consequence of the number of columns of the associated submatrix being smaller than the number of its rows as mentioned in Theorem 1. For some finite elements such as for , the associated submatrix is nearly square and as shown in Figure 3, it is also nearly full rank. In this situation, the so-called locking phenomena may occur in the sense that the associated submatrix may represent an injective operator, and thus, it only admits a unique solution.

Our computations suggest that out of possible choices of simplicial elements of Figure 1 for displacement, displacement gradient, stress, and pressure, choices in D and choices in D satisfy all inf-sup conditions. More specifically, the choices that violate each inf-sup condition are as follow:

• In D and D, the inf-sup condition (5.4) does not hold for the choices of elements and for displacement and stress.

• The inf-sup condition (5.5) does not hold with the following choices of elements for displacement gradient and pressure:

• 2D: and ;

• 3D: and , .

• The inf-sup condition (5.6) does not hold with the following choices of elements for displacement, displacement gradient, and stress:

• 2D: , , , and , ;

• 3D: and , .

• In D and D, the inf-sup condition (5.8) does not hold with the following choices of elements for displacement, displacement gradient, stress, and pressure: , , , and , , .

• The inf-sup condition (5.11) does not hold with the following choices of elements for displacement, displacement gradient, and pressure:

• 2D: and , ;

• 3D: and , , .

• The inf-sup condition (5.14) does not hold with the following choices of elements for displacement gradient, stress, and pressure:

• 2D: and , ;

• 3D: , , , and , .

The inf-sup conditions of Section 5 are not independent. For example, the inf-sup conditions (5.4) and (5.6) are closely related. Such dependencies can be read off from the above numerical results. Common strategies may be employed to stabilize choices that violate the above inf-sup conditions such as using different underlying meshes for different unknowns, enriching trial and test spaces using bubble functions, and using perturbed formulations [9].