Compressible Neo-Hookean materials

Author: Alex Sadovsky

Notation

In what follows, the symbols \(\mathbf{u}, \bF, \bB, \bC, \stress\) denote, respectively, the displacement field, the deformation gradient, the left Cauchy-Green strain tensor \(\bB = \bF \bF^T\), the right Cauchy-Green strain tensor \(\bC =\bF^T \bF\), and the Cauchy stress tensor.

We also introduce the symbols \(I_1 := \tr \bC\) and \(J := \det\bF\). Use will be made of the identity:

\[{\p J \over \p \bC} = J \bC^{-1}\]

The symbol \(\Id\) denotes the identity tensor. The symbol \(\Omega_{0}\) denotes the reference configuration of the body to be deformed. The unit volume in the reference (resp., deformed) configuration is denoted \(dV\) (resp., \(dV_{0}\)); these two are related by:

\[dV = J dV_{0},\]

which allows an integral over \(\Omega\) involving the Cauchy stress \(\bT\) to be rewritten as an integral of the Kirchhoff stress \(\kappa = J \bT\) over \(\Omega_{0}\).

Recommended References

For an exposition of nonlinear elasticity and of the underlying linear and tensor algebra, see [OGDEN1984]. For an advanced mathematical analysis of the Finite Element Method, see [RAVIART1998].

A Neo-Hookean Compressible Material

Constitutive Theory and Tangent Stress Measures

The strain energy density function is given by:

\[W = {\mu \over 2}(I_1 - \tr \Id - 2 \ln J)\]

(see [HORGAN2004], formula (12)).

The corresponding 2nd Piola-Kirchoff stress tensor is given by:

\[\bS_{n} := {\p W \over \p\bE} (\bF_{n}) = \mu (\Id - \bC^{-1})\]

The Kirchhoff stress, then, is:

\[\kappa = \bF \bS \bF^{T} = \mu (\bB - \Id)\]

The tangent Kirchhoff stress tensor at \(\bF_{n}\) acting on \(\delta \bF_{n+1}\) is, consequently:

\[{\p \kappa \over \p \bF} (\bF_{n}) \delta \bF_{n+1} = \mu \left[ \bF_{n} (\delta \bF_{n+1})^T + \delta \bF_{n+1} (\bF_{n})^T \right]\]

The Weak Form of the BVP in the Absence of Body (External) Forces

The \(\Omega_0\) we are considering is an elliptical annulus, whose boundary consists of two concentric ellipses (each allowed to be a circle as a special case), with the major axes parallel. Let \(P\) denote the dead stress load (traction) on a portion \(\partial \Omega_0^{t}\) (= the inner ellipse) of the boundary \(\partial \Omega_0\). On the rest of the boundary, we prescribe zero displacement.

The weak formulation of the boundary value problem is:

\[\begin{split}\arr{lll} 0 & = & \int_{\Omega_0} \kappa[\bF] \: : \: \left\{ (\Grad \otimes \mathbf{w}) (\bF)^{-1} \right\}\\ & - & \int_{\p \Omega_0^{t}} P \cdot \hat{N}_0\\ \rra\end{split}\]

For brevity, in the rest of this section we assume \(P = 0\). The provided FreeFEM code, however, does not rely on this assumption and allows for a general value and direction of \(P\).

Given a Newton approximation \(\mathbf{u}_n\) of the displacement field \(\mathbf{u}\) satisfying the BVP, we seek the correction \(\delta \mathbf{u}_{n+1}\) to obtain a better approximation:

\[\mathbf{u}_{n+1} = \mathbf{u}_{n} + \delta \mathbf{u}_{n+1}\]

by solving the weak formulation:

\[\begin{split}\arr{lll} 0 &=& \int_{\Omega_0}\kappa[\bF_{n} + \delta \bF_{n+1}]\: :\: \left\{(\Grad \otimes \mathbf{w}) (\bF_{n} + \delta\bF_{n+1})^{-1}\right\}- \int_{\p \Omega_0} P \cdot \hat{N}_0\\ &=& \int_{\Omega_0}\left\{\kappa[\bF_{n}] + {\p \kappa \over \p \bF}[\bF_{n}]\delta \bF_{n+1}\right\}\: :\: \left\{(\Grad \otimes \mathbf{w})(\bF_{n} + \delta \bF_{n+1})^{-1}\right\}\\ &=& \int_{\Omega_0}\left\{\kappa[\bF_{n}] + {\p \kappa \over \p \bF}[\bF_{n}]\delta \bF_{n+1}\right\}\: :\: \left\{(\Grad \otimes \mathbf{w}) (\bF_{n}^{-1} + \bF_{n}^{-2} \delta \bF_{n+1})\right\}\\ \\ &=& \int_{\Omega_0}\kappa[\bF_{n}]\: :\: \left\{(\Grad \otimes \mathbf{w})\bF_{n}^{-1}\right\}\\ &-& \int_{\Omega_0}\kappa[\bF_{n}]\: :\: \left\{(\Grad \otimes \mathbf{w})(\bF_{n}^{-2} \delta \bF_{n+1})\right\}\\ &+& \int_{\Omega_0}\left\{{\p \kappa \over \p \bF}[\bF_{n}]\delta \bF_{n+1}\right\}\: :\: \left\{(\Grad \otimes \mathbf{w}) \bF_{n}^{-1} \right\} \\ \rra \quad \mbox{for all test functions} \mathbf{w},\end{split}\]

where we have taken:

\[\delta \bF_{n+1} = \Grad \otimes \delta \mathbf{u}_{n+1}\]

Note

Contrary to standard notational use, the symbol \(\delta\) here bears no variational context. By \(\delta\) we mean simply an increment in the sense of Newton’s Method. The role of a variational virtual displacement here is played by \(\mathbf{w}\).

An Approach to Implementation in FreeFEM

Introducing the code-like notation, where a string in \(< >\)’s is to be read as one symbol, the individual components of the tensor:

\[<TanK> := {\p \kappa \over \p \bF}[\bF_{n}] \delta \bF_{n+1}\]

will be implemented as the macros \(<TanK11>\), \(<TanK12>\), …

The individual components of the tensor quantities:

\[\bD_{1} := \bF_{n} (\delta \bF_{n+1})^T + \delta \bF_{n+1} (\bF_{n})^T,\]

\[\bD_{2} := \bF_{n}^{-T} \delta \bF_{n+1},\]

\[\bD_{3} := (\Grad \otimes \mathbf{w}) \bF_{n}^{-2} \delta \bF_{n+1},\]

and

\[\bD_{4} := (\Grad \otimes \mathbf{w}) \bF_{n}^{-1},\]

will be implemented as the macros:

\[\begin{split}\arr{l} <d1Aux11>, <d1Aux12>, \quad \ldots \quad, <d1Aux22>,\\ <d2Aux11>, <d2Aux12>, \quad \ldots \quad, <d2Aux22>\\ <d3Aux11>, <d3Aux12>, \quad \ldots \quad, <d3Aux22>\\ <d4Aux11>, <d4Aux12>, \quad \ldots \quad, <d4Aux22>\\ \rra,\end{split}\]

respectively.

In the above notation, the tangent Kirchhoff stress term becomes

\[{\p \kappa \over \p \bF} (\bF_{n}) \: \delta \bF_{n+1} = \mu \: \bD_{1}\]

while the weak BVP formulation acquires the form:

\[\begin{split}\arr{lll} 0 & = & \int_{\Omega_0} \kappa[\bF_{n}] \: : \: \bD_{4} \\ &-& \int_{\Omega_0} \kappa[\bF_{n}] \: : \: \bD_{3} \\ &+& \int_{\Omega_0} \left\{ {\p \kappa \over \p \bF}[\bF_{n}] \delta \bF_{n+1} \right\} \: : \: \bD_{4} \\ \rra \quad \mbox{for all test functions} \mathbf{w}\end{split}\]

// Macro
//Macros for the gradient of a vector field (u1, u2)
macro grad11(u1, u2) (dx(u1)) //
macro grad21(u1, u2) (dy(u1)) //
macro grad12(u1, u2) (dx(u2)) //
macro grad22(u1, u2) (dy(u2)) //

//Macros for the deformation gradient
macro F11(u1, u2) (1.0 + grad11(u1, u2)) //
macro F12(u1, u2) (0.0 + grad12(u1, u2)) //
macro F21(u1, u2) (0.0 + grad21(u1, u2)) //
macro F22(u1, u2) (1.0 + grad22(u1, u2)) //

//Macros for the incremental deformation gradient
macro dF11(varu1, varu2) (grad11(varu1, varu2)) //
macro dF12(varu1, varu2) (grad12(varu1, varu2)) //
macro dF21(varu1, varu2) (grad21(varu1, varu2)) //
macro dF22(varu1, varu2) (grad22(varu1, varu2)) //

//Macro for the determinant of the deformation gradient
macro J(u1, u2) (
      F11(u1, u2)*F22(u1, u2)
    - F12(u1, u2)*F21(u1, u2)
) //

//Macros for the inverse of the deformation gradient
macro Finv11 (u1, u2) (
      F22(u1, u2) / J(u1, u2)
) //
macro Finv22 (u1, u2) (
      F11(u1, u2) / J(u1, u2)
) //
macro Finv12 (u1, u2) (
    - F12(u1, u2) / J(u1, u2)
) //
macro Finv21 (u1, u2) (
    - F21(u1, u2) / J(u1, u2)
) //

//Macros for the square of the inverse of the deformation gradient
macro FFinv11 (u1, u2) (
      Finv11(u1, u2)^2
    + Finv12(u1, u2)*Finv21(u1, u2)
) //

macro FFinv12 (u1, u2) (
      Finv12(u1, u2)*(Finv11(u1, u2)
    + Finv22(u1, u2))
) //

macro FFinv21 (u1, u2) (
      Finv21(u1, u2)*(Finv11(u1, u2)
    + Finv22(u1, u2))
) //

macro FFinv22 (u1, u2) (
      Finv12(u1, u2)*Finv21(u1, u2)
    + Finv22(u1, u2)^2
) //

//Macros for the inverse of the transpose of the deformation gradient
macro FinvT11(u1, u2) (Finv11(u1, u2)) //
macro FinvT12(u1, u2) (Finv21(u1, u2)) //
macro FinvT21(u1, u2) (Finv12(u1, u2)) //
macro FinvT22(u1, u2) (Finv22(u1, u2)) //

//The left Cauchy-Green strain tensor
macro B11(u1, u2) (
      F11(u1, u2)^2 + F12(u1, u2)^2
)//

macro B12(u1, u2) (
      F11(u1, u2)*F21(u1, u2)
    + F12(u1, u2)*F22(u1, u2)
)//

macro B21(u1, u2) (
      F11(u1, u2)*F21(u1, u2)
    + F12(u1, u2)*F22(u1, u2)
)//

macro B22(u1, u2)(
      F21(u1, u2)^2 + F22(u1, u2)^2
)//

//The macros for the auxiliary tensors (D0, D1, D2, ...): Begin
//The tensor quantity D0 = F{n} (delta F{n+1})^T
macro d0Aux11 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * F11(u1, u2)
    + dF12(varu1, varu2) * F12(u1, u2)
)//

macro d0Aux12 (u1, u2, varu1, varu2) (
      dF21(varu1, varu2) * F11(u1, u2)
    + dF22(varu1, varu2) * F12(u1, u2)
)//

macro d0Aux21 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * F21(u1, u2)
    + dF12(varu1, varu2) * F22(u1, u2)
)//

macro d0Aux22 (u1, u2, varu1, varu2) (
      dF21(varu1, varu2) * F21(u1, u2)
    + dF22(varu1, varu2) * F22(u1, u2)
)//

///The tensor quantity D1 = D0 + D0^T
macro d1Aux11 (u1, u2, varu1, varu2) (
      2.0 * d0Aux11 (u1, u2, varu1, varu2)
)//

macro d1Aux12 (u1, u2, varu1, varu2) (
      d0Aux12 (u1, u2, varu1, varu2)
    + d0Aux21 (u1, u2, varu1, varu2)
)//

macro d1Aux21 (u1, u2, varu1, varu2) (
      d1Aux12 (u1, u2, varu1, varu2)
)//

macro d1Aux22 (u1, u2, varu1, varu2)(
      2.0 * d0Aux22 (u1, u2, varu1, varu2)
)//

///The tensor quantity D2 = F^{-T}_{n} dF_{n+1}
macro d2Aux11 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * FinvT11(u1, u2)
    + dF21(varu1, varu2) * FinvT12(u1, u2)
)//

macro d2Aux12 (u1, u2, varu1, varu2) (
      dF12(varu1, varu2) * FinvT11(u1, u2)
    + dF22(varu1, varu2) * FinvT12(u1, u2)
)//

macro d2Aux21 (u1, u2, varu1, varu2) (
      dF11(varu1, varu2) * FinvT21(u1, u2)
    + dF21(varu1, varu2) * FinvT22(u1, u2)
)//

macro d2Aux22 (u1, u2, varu1, varu2) (
      dF12(varu1, varu2) * FinvT21(u1, u2)
    + dF22(varu1, varu2) * FinvT22(u1, u2)
)//

///The tensor quantity D3 = F^{-2}_{n} dF_{n+1}
macro d3Aux11 (u1, u2, varu1, varu2, w1, w2) (
      dF11(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
    + dF21(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
)//

macro d3Aux12 (u1, u2, varu1, varu2, w1, w2) (
      dF12(varu1, varu2) * FFinv11(u1, u2) * grad11(w1, w2)
    + dF22(varu1, varu2) * FFinv12(u1, u2) * grad11(w1, w2)
    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad12(w1, w2)
    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad12(w1, w2)
)//

macro d3Aux21 (u1, u2, varu1, varu2, w1, w2) (
      dF11(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
    + dF21(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
    + dF11(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
    + dF21(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
)//

macro d3Aux22 (u1, u2, varu1, varu2, w1, w2) (
      dF12(varu1, varu2) * FFinv11(u1, u2) * grad21(w1, w2)
    + dF22(varu1, varu2) * FFinv12(u1, u2) * grad21(w1, w2)
    + dF12(varu1, varu2) * FFinv21(u1, u2) * grad22(w1, w2)
    + dF22(varu1, varu2) * FFinv22(u1, u2) * grad22(w1, w2)
)//

///The tensor quantity D4 = (grad w) * Finv
macro d4Aux11 (w1, w2, u1, u2) (
      Finv11(u1, u2)*grad11(w1, w2)
    + Finv21(u1, u2)*grad12(w1, w2)
)//

macro d4Aux12 (w1, w2, u1, u2) (
      Finv12(u1, u2)*grad11(w1, w2)
    + Finv22(u1, u2)*grad12(w1, w2)
)//

macro d4Aux21 (w1, w2, u1, u2) (
      Finv11(u1, u2)*grad21(w1, w2)
    + Finv21(u1, u2)*grad22(w1, w2)
)//

macro d4Aux22 (w1, w2, u1, u2) (
      Finv12(u1, u2)*grad21(w1, w2)
    + Finv22(u1, u2)*grad22(w1, w2)
)//
//The macros for the auxiliary tensors (D0, D1, D2, ...): End

//The macros for the various stress measures: BEGIN
//The Kirchhoff stress tensor
macro StressK11(u1, u2) (
      mu * (B11(u1, u2) - 1.0)
)//

//The Kirchhoff stress tensor
macro StressK12(u1, u2) (
      mu * B12(u1, u2)
)//

//The Kirchhoff stress tensor
macro StressK21(u1, u2) (
      mu * B21(u1, u2)
)//

//The Kirchhoff stress tensor
macro StressK22(u1, u2) (
      mu * (B22(u1, u2) - 1.0)
)//

//The tangent Kirchhoff stress tensor
macro TanK11(u1, u2, varu1, varu2) (
      mu * d1Aux11(u1, u2, varu1, varu2)
)//

macro TanK12(u1, u2, varu1, varu2) (
      mu * d1Aux12(u1, u2, varu1, varu2)
)//

macro TanK21(u1, u2, varu1, varu2) (
      mu * d1Aux21(u1, u2, varu1, varu2)
)//

macro TanK22(u1, u2, varu1, varu2) (
      mu * d1Aux22(u1, u2, varu1, varu2)
)//
//The macros for the stress tensor components: END

// Parameters
real mu = 5.e2; //Elastic coefficients (kg/cm^2)
real D = 1.e3; //(1 / compressibility)
real Pa = -3.e2; //Stress loads

real InnerRadius = 1.e0; //The wound radius
real OuterRadius = 4.e0; //The outer (truncated) radius
real tol = 1.e-4; //Tolerance (L^2)
real InnerEllipseExtension = 1.e0; //Extension of the inner ellipse ((major axis) - (minor axis))

int m = 40, n = 20;

// Mesh
border InnerEdge(t=0, 2.*pi){x=(1.0 + InnerEllipseExtension)*InnerRadius*cos(t); y=InnerRadius*sin(t); label=1;}
border OuterEdge(t=0, 2.*pi){x=(1.0 + 0.0*InnerEllipseExtension)*OuterRadius*cos(t); y=OuterRadius*sin(t); label=2;}
mesh Th = buildmesh(InnerEdge(-m) + OuterEdge(n));
int bottom = 1, right = 2, upper = 3, left = 4;
plot(Th);

// Fespace
fespace Wh(Th, P1dc);
fespace Vh(Th, [P1, P1]);
Vh [w1, w2], [u1n,u2n], [varu1, varu2];
Vh [ehat1x, ehat1y], [ehat2x, ehat2y];
Vh [auxVec1, auxVec2]; //The individual elements of the total 1st Piola-Kirchoff stress
Vh [ef1, ef2];

fespace Sh(Th, P1);
Sh p, ppp;
Sh StrK11, StrK12, StrK21, StrK22;
Sh u1, u2;

// Problem
varf vfMass1D(p, q) = int2d(Th)(p*q);
matrix Mass1D = vfMass1D(Sh, Sh, solver=CG);

p[] = 1;
ppp[] = Mass1D * p[];

real DomainMass = ppp[].sum;
cout << "DomainMass = " << DomainMass << endl;

varf vmass ([u1, u2], [v1, v2], solver=CG)
    = int2d(Th)( (u1*v1 + u2*v2) / DomainMass );

matrix Mass = vmass(Vh, Vh);

matrix Id = vmass(Vh, Vh);

//Define the standard Euclidean basis functions
[ehat1x, ehat1y] = [1.0, 0.0];
[ehat2x, ehat2y] = [0.0, 1.0];

real ContParam, dContParam;

problem neoHookeanInc ([varu1, varu2], [w1, w2], solver=LU)
    = int2d(Th, qforder=1)(
        - (
              StressK11 (u1n, u2n) * d3Aux11(u1n, u2n, varu1, varu2, w1, w2)
            + StressK12 (u1n, u2n) * d3Aux12(u1n, u2n, varu1, varu2, w1, w2)
            + StressK21 (u1n, u2n) * d3Aux21(u1n, u2n, varu1, varu2, w1, w2)
            + StressK22 (u1n, u2n) * d3Aux22(u1n, u2n, varu1, varu2, w1, w2)
        )
        + TanK11 (u1n, u2n, varu1, varu2) * d4Aux11(w1, w2, u1n, u2n)
        + TanK12 (u1n, u2n, varu1, varu2) * d4Aux12(w1, w2, u1n, u2n)
        + TanK21 (u1n, u2n, varu1, varu2) * d4Aux21(w1, w2, u1n, u2n)
        + TanK22 (u1n, u2n, varu1, varu2) * d4Aux22(w1, w2, u1n, u2n)
    )
    + int2d(Th, qforder=1)(
          StressK11 (u1n, u2n) * d4Aux11(w1, w2, u1n, u2n)
        + StressK12 (u1n, u2n) * d4Aux12(w1, w2, u1n, u2n)
        + StressK21 (u1n, u2n) * d4Aux21(w1, w2, u1n, u2n)
        + StressK22 (u1n, u2n) * d4Aux22(w1, w2, u1n, u2n)
    )
    //Choose one of the following two boundary conditions involving Pa:
    // Load vectors normal to the boundary:
    - int1d(Th, 1)(
          Pa * (w1*N.x + w2*N.y)
    )
    //Load vectors tangential to the boundary:
    //- int1d(Th, 1)(
    //    Pa * (w1*N.y - w2*N.x)
    //)
    + on(2, varu1=0, varu2=0)
    ;

//Auxiliary variables
matrix auxMat;

// Newton's method
ContParam = 0.;
dContParam = 0.01;

//Initialization:
[varu1, varu2] = [0., 0.];
[u1n, u2n] = [0., 0.];
real res = 2. * tol;
real eforceres;
int loopcount = 0;
int loopmax = 45;

// Iterations
while (loopcount <= loopmax && res >= tol){
    loopcount ++;
    cout << "Loop " << loopcount << endl;

    // Solve
    neoHookeanInc;

    // Update
    u1 = varu1;
    u2 = varu2;

    // Residual
    w1[] = Mass*varu1[];
    res = sqrt(w1[]' * varu1[]); //L^2 norm of [varu1, varu2]
    cout << " L^2 residual = " << res << endl;

    // Newton
    u1n[] += varu1[];

    // Plot
    plot([u1n,u2n], cmm="displacement");
}

// Plot
plot(Th, wait=true);

// Movemesh
mesh Th1 = movemesh(Th, [x+u1n, y+u2n]);

// Plot
plot(Th1, wait=true);
plot([u1n,u2n]);

Here you can use the search bar to search in all FreeFEM examples instead

The list of all matching examples on the right is updated for each new query

Clicking on a search hit or on an example on the right brings up the corresponding script from GitHub

You can also filter the list on the right by selecting one or more intersecting categories on the left tree