1. Introduction

In this blog, we continue towards building a robust solver for the St-Venant Kirchhoff hyperelastic model. In our previous blog post, we delved into the challenges that arise when solving for

in a 1D setting, and now we investigate the challenges that arise in a 2D setting. In particular, we are interested in investigating how the St-Venant Kirchhoff model behaves in a well-known benchmark case: the incompressible block compression benchmark; see Fig. 1 for an illustration of the domain and problem.

As we have already observed in our previous blog post, for large compressive forces the St-Venant Kirchhoff model struggles to converge to a solution, or it may converge to a solution that is orientation preserving, i.e., a non-physical solution, and hence is not the best model to simulate the block compression scenario. However, we are interested to see how far we can push the model, and if there is a way to make a slight modification which can help us better approximate and capture the physical displacement profile.

To this end, we analyze the computational solver and …

We also work on modifying the St-Venant Kirchhoff model by looking at the problem from a point of view of convergence and well-posedness, and analyze the results for the new model. In other words, we explore how a simple stabalization of the model can be achieved, and discuss the improvements…

Mention locking…

2. Details of the Computational Solver

Following the notation we set up in our previous post, we recall the St-Venant Kirchhoff model

where $\lambda$ and $\mu$ are Lamé parameters. We concern ourselves with only the reference configuration $\Omega$ in 2D, i.e., $\Omega \subset \mathbb{R}^2$. We consider the system \eqref{eq:governing_system} with mixed boundary conditions

where $t_N$ [Pa] is the given traction which is assumed to be smooth enough (we are mostly concerned with constant values here), and the Dirichlet boundary condition on $u = (u_1, u_2)$ is equivalent to

2.1. Numerical Discretization

The variational problem associated with the above system \eqref{eq:governing_system}-\eqref{eq:mixed_bc} is given as follows: we seek $u_h \in V_h$ such that ¹

where

is the finite dimensional subspace of Q1 bilinear elements on $\Omega_h$. To avoid ambiguity due to notation, we highlight that in \eqref{eq:variational_form}, the basis elements $\phi_i : \mathbb{\Omega}_h \rightarrow \mathbb{R}^2$ are vector-valued functions, and

For any $\phi_h \in V_h$, the $H^1$ semi-norm $\vert \cdot \vert_2$ is defined as

Using the Poincaré inequality, $\vert \cdot \vert_2$ becomes a norm on $(H^1_0(\Omega))^2$.

2.2. Nonlinear Solver: Newton’s Method

We now recall some details of Newton’s method, which we use to solve the nonlinear system of equations resulting from \eqref{eq:variational_form}. Letting $\{ \phi_i \}$ denote the $N$ basis elements for $V_h$, we consider the residuals defined as

and thus we seek a solution $U = \left[U_1 \; U_2 \; \dots \; U_N \right]^T \in \mathbb{R}^{N}$ to

where $u_h = \sum_{i=1}^N U_i \phi_i$ and $\mathcal{T} = [\mathcal{T}_1 \; \mathcal{T}_2 \; \dots \mathcal{T}_N]^T$.

The Newton’s method to solve for \eqref{eq:nonlinear_eq} gives the following algorithm: starting from $u^{(0)}_h \in V_h$, perform the update ²

where $\mathcal{T}^{(m-1)} = \mathcal{T}(U^{(m-1)})$ is the residual at $U^{(m-1)}$ and $\mathcal{J}^{(m-1)} = \left(\nabla_U \mathcal{T}\right) (U^{(m-1)})$ is the Jacobian of $\mathcal{T}$ (here $\nabla_U$ denotes the gradient with respect to $U$).

We now investigate the properties of the Jacobian $\mathcal{J}$ which will help us select an appropriate linear solver to compute its inverse $\mathcal{J}^{-1}$ when performing the update \eqref{eq:newton_update}.

2.2.1. Symmetric Property of the Jacobian

We begin by simplifying the Jacobian. Note that from \eqref{eq:nonlinear_residual}, we have for any given $U$

Let us simplify the expression \eqref{eq:jacobian_def}. First, note that by definition we have

and thus

Now, since $E = \frac{1}{2}\left(F^T F - I \right)$, we have from \eqref{eq:F_gradient}

Using \eqref{F_gradient} and \eqref{eq:E_gradient} in \eqref{eq:st_venant_kirchhoff_def} we get

Thus, we have the expression

We now prove an that the Jacobian is symmetric.

Lemma 2.2.1.1. The Jacobian $\mathcal{J}$ defined by \eqref{eq:jacobian_def} is symmetric.

Proof. We consider individual terms in the integrand of \eqref{eq:jacobian_def_2}. Note that the term

is trivially symmetric (by interchanging $i$ and $j$). We now consider the other terms.

Note that for any two square matrices $A$ and $B$, the double contraction operator is defined as

and $\text{tr}(AB) = \text{tr}(BA)$.

We now consider the term

First note that using the above identities $\text{tr}\left(\nabla \phi_j^T F \right) = \text{tr}\left( F^T \nabla \phi_j \right)$. Also, we can rewrite \eqref{eq:second_term} as

which is symmetric in $i$ and $j$. Thirdly, we consider the term

which we rewrite as

Further, using commutativity, we can rewrite \eqref{eq:third_term_2} as

which is symmetric in $i$ and $j$. Fourthly, and finally, we consider the term

which we rewrite as

Using the definition of $:$ we can rewrite \eqref{eq:fourth_term} as

\(\text{tr} \left(F^T \nabla \phi_j F^T \nabla \phi_i \right) + F^T \nabla phi_j : F^T \nabla \phi_i.\) By the commutativity of trace and $:$, the symmetry in $i$ and $j$ follows. This proves that $\mathcal{J}$ is symmetric.

□

The above lemma guides us in choosing a linear solver to compute the inverse $\mathcal{J}^{-1}$. We choose a solver that leverages the symmetric property of $\mathcal{J}$, but we cannot naively expect too much out of, say, Conjugate Gradient since $\mathcal{J}$ need not be positive definite. Thus, we proceed with minimal residual (MINRES) method to handle our symmetric indefinite Jacobian ³. However, note that the above lemma does not prove that $\mathcal{J}$ is invertible, and proving that would require investigating the bilinear operator associated with each Newton update, which we do not discuss here.

2.2. Anticipated Challenges

In our previous post, we exemplified the issues of existence and uniqueness for simple scenarios. This led to, for example, (i) lack of convergence of computational solver for refined grids, or (ii) convergence to different solution profiles for different initial guesses in the Newton’s method. Indeed, there is no reason for us to not expect such challenges in 2D (or 3D). Moreover, a very well known challenge in mechanics is expected to arise given our choice of Q1 finite elements: locking ⁴. Locking occurs for incompressible materials, and is defined as the inability of the computational solver to provide an accurate solution profile (displacements) for coarse meshes. For linear elasticity, locking can be explained by loss of coercivity of the associated bilinear form when $\frac{\lambda}{\mu} \rightarrow \infty$, which leads to large error for a given grid size $h$ ⁵.

Locking is well-studied in literature, and there exists multiple formulations (eg. mixed, weak Galerkin etc.) that are known to reduce the effect of locking. We will explore the extent to which locking affects our solution with the (naive) Q1 elements.

3. Numerical Results

We now investigate the performance of our computational solver in physical scenarios. We implement the computational solver using the finite element library deal.II ⁶. We consider the MINRES solver as the linear solver when computing the solution $\delta U = \mathcal{J}^{-1} \mathcal{T}$. The absolute and relative tolerances to which we solve the system are $\epsilon_{ns, abs} = 10^{-8}$ and $\epsilon_{ns, rel} = 10^{-10}$. The MINRES algorithm stops when a tolerance of $\epsilon_{ls} = 10^{-8} \times \lVert \mathcal{T} \rVert_2$ is achieved.

3.1. Clamped Bar Under Dead Load

We begin with the example visited in our (previous post)[(https://nvohra0016.github.io/talks/hyperelasticity1/)], i.e., of a clamped bar under a constant force. We consider $\Omega = (0, 1) \times (0, 0.4)$ [m $^2$], and the external force $f = 3 \times 10^6$ [N / m $^3$], along with the following boundary conditions

The Young’s modulus is $E_Y = 10^7$ [Pa] and the Poisson ratio is $\nu = 0.48$ [-]. The results are shown in Fig. 2 for a uniform grid size of $0.025 \times 0.025$ [m $^2$].

The profile shows the bar deformed… We also tabulate the performance of the computational solver in Table 1. It can be observed that the number of nonlinear and linear iterations increase as the mesh is refined, and we quickly reach a stage of no convergence for a small enough grid size $h$. No locking was observed in this example, which is unsurprising since the material is not highly incompressible (owing to $\nu = 0.48$).

It is also noted that when the force $f$ is increased further, the computational solver does not converge (or exhibits non-physical scenarios associated with det$(F) < 0$). Moreover, the LS iter. increase proportionally to the size of the system, i.e., the number of DoFs N, which is not ideal as we will encounter the “curse of dimensionality” as we move into 3D. Although, this we can hope to resolve to some extent with a suitable preconditioner.

3.2. Incompressible Block Compression

We now proceed with investigating the performance of the computational solver on an incompressible block compression scenario as depicted in Fig. 1. The benchmark is well-known and is taken from ⁷ among other references, and is usually simulated using a Neo-Hookean model. The material parameters are chosen from ⁸. In particular, the Young’s modulus is taken to be $E_Y = 240.565 \times 10^6$ [Pa] and the Poisson ratio is taken to be $\nu = 0.4999$ [-].

Re-selecting $t_N$ to obtain convergence. We consider traction vector as $t_N = (-4 \times 10^8, 0)$ [Pa] as in ⁸, however, our computational solver does not converge. To obtain convergence, we reduce the traction magnitude to $t_N = -7 \times 10^7$ (after trial and error with different grid sizes). The results are shown in Fig. 3 for different grid sizes.

We can now observe some severe locking taking place in the solution. Indeed, the deformed configurations look very different and the maximum displacement magnitude of the two profiles differs by a factor of ~$1.7$. Moreover, the maximum vertical displacement we are able to achieve is $2.95$ [mm] at point $(10, 10)$, whereas the results reported in ⁸ are close to $6$ [mm], a big difference.

On the solver side, for grid size $h = 1$ and $h = 0.02$ convergence is achieved in NS iter. 6 and LS iter. 3922, and NS iter. 14 and LS iter. 138170, numbers which can definitely take major improvements.

4. Improving the Performance of the Solver: Adding Stabilization

The above examples show how poorly the St-Venant Kirchhoff model performs for typical physical scenarios, on the physical and computational side. In particular, for large forces or traction, the solver does not converge, or the solution converges to an un-physical profile (such as where det$(F) < 0$), all of which is exacerbated by the issue of existence and uniqueness. Let us now now work towards improving the model through a stabilizing term.

Note. The discussion below is not found on any experimental data, and is just to explore how a material responds to a “modified” St-Venant Kirchhoff type model. We try to provide mathematical rigour on some clear assumptions, which, however, are not guaranteed to hold true for all physical scenarios, although we have seen agreeable improvement in some examples.

Let $\gamma > 0$ be a parameter (we list the exact dependence of $\gamma$ later, and here we note that the units of $\gamma$ are [Pa]). Consider now the modified constitutive law

where we have now introduced an diffusion term with the material (when we take the divergence of \eqref{eq:modified_law}). The reader may link this to the concept of artificial diffusion to improve the stability of numerical methods, but physically this terms adds extra stiffness at large displacement values. That is, we now expect a stiffer response from the material, but only enough so that we obtain physically sound solution profiles when external forces are large. In terms of the strain energy function, this is equivalent to

Let us now work towards finding a solution…

The variational form now becomes: find $u_h \in V_h$ such that

where now it is understood that the traction boundary condition is applied as $\widetilde{T} n = t_N$.

Now we explore the implications of adding this additional term. On the computational side, one thing becomes clear that if $\gamma$ is large enough, then the Jacobian associated with $\widetilde{T}$ becomes symmetric positive definite (SPD) since

where

is an SPD matrix. Thus for large $\gamma$, the eigenvalues of $\widetilde{\mathcal{J}}$ can be made positive since $\mathcal{A}$ has positive eigenvalues (it is SPD) and $\mathcal{J}$ has real eigenvalues (it is symmetric). This further guarantees existence of a unique update $\delta U$ in each Newton iteration.

Now let us consider the variational form \eqref{eq:modified_variational_form} with homogeneous Dirichlet boundary conditions for simplicity and clarity of exposition. For the exposition below, we use $V_h \subset H_0^1(\Omega)$ the finite dimensional subspace of piecewise-linear functions which vanish on $\partial \Omega$. We first rewrite \eqref{eq:modified_variational_form} using a nonlinear operator $a : V_h \rightarrow V_h’$, where $V_h’$ is the dual of $V_h$: we seek $u_h \in V_h$ such that

where

and $S \in V_h’$ is defined as

We now investigate the properties of the operator $a$. In particular, we are interested in establishing some form of coercivity and monotonicity which will guarantee the existence of a unique solution using Minty-Browder’s theorem (Ciarlet ⁹, Theorem 9.14-1). That, however, is not a trivial task, and instead we will prove the properties in 1D. We also make use of a cut-off operator to compute the tensors $F$ and $E$ as

where

and the constant $C_0 > 0$ is taken to be large enough. That is, the cut-off operator $C : \mathbb{R} \rightarrow \mathbb{R}$ is defined as

such that

Note that $C$ is also Lipschitz continuous, and we denote the Lipschitz constant by $L_C$.

Now, define

and consider the following variational formulation: we seek $u_h \in V_h$ such that

where

and $S$ is defined as in \eqref{eq:S_def}, i.e.,

Note on the cut-off operator. The idea of the cut-off operator is that for a large enough $C_0$, the system \eqref{eq:modified_aC} is close to \eqref{eq:modified_problem}, and hence a numerical implementation of the cut-off operator is unnecessary for large $C_0$. However, we rely on the cut-off operator for formal analysis and an existence result as proved below. For an application of the cut-off operator in thermo-poroelastic setting, see ¹⁰.

We now state an existence result making use of the concepts of monotonicity, coercivity, and hemicontinuity of operators. For a brief introduction and definitions, see also the monograph ¹⁰. Note that since $V_h$ is a finite dimensional Hilbert space, it is reflexive and separable.

Theorem 4.1. Fix $h > 0$ and $C_0 > 0$. Then, the operator $a_C : V_h \rightarrow V_h’$ is strictly monotone, hemicontinuous, and coercive for a large enough $\gamma > 0$. Thus $\exists$ a unique solution to \eqref{eq:modified_aC} for a smooth enough $f$.

Proof. First note that since $\frac{du_h}{dX}$ is a piece-wise constant function, we can rewrite

where the polynomial $g_1 : \mathbb{R} \rightarrow \mathbb{R}$ is defined as $g_1(\alpha) = \alpha(\alpha + 1)(\alpha + 2)$. Now, let $u_h, v_h \in V_h$. To prove that $a_C$ is hemicontinuous, we need to show that

is continuous on an interval $(-t_0, t_0), t_0 > 0$ for any $u_h, v_h, \phi_h \in V_h$. This follows easily since if $\{ t_n \}$ is a sequence such that $t_n \rightarrow 0$ then

The second term in \eqref{eq:step-1_proof} can be bounded by Hölder’s inquality since $v_h, \phi_h \in H^1_0(\Omega_h)$, and hence $\rightarrow 0$ as $t \rightarrow 0$. For the first term, note that since $g$ is a polynomial, and $(x - y)$ divivides $(x^m - y^m)$ for any $x, y \in \mathbb{R}$ and $m \in \mathbb{Z}$, $m > 0$, we can write $g(x) - g(y) = (x - y) g_2(x, y)$, for some polynomial $g_2 : \mathbb{R} \times \mathbb{R} \rightarrow \mathbb{R}$. Thus, we can rewrite the first term in \eqref{eq:step-1_proof} as

By the Lipschitz continuity of $C$, and using the boundedness of $C$, we have that $A_0 \rightarrow 0$ as $t \rightarrow 0$. Thus, $B(t)$ is hemicontinuous at $t = 0$.

Now we work towards the montononicity of $a_C$. We have

Using the factorization above, we can rewrite \eqref{eq:step0_proof} as

Now, by using the first inequality in \eqref{eq:C_bounds} and its Lipschitz continuity we can bound the first term in \eqref{eq:step1_proof}

by

where $C_1 = C_1(C_0) = O(C_0^2)$ is some constant. Thus, we have from \eqref{eq:step2_proof} and using the triangle inequality

It can be observed that the RHS of \eqref{eq:step4_proof} is positive for a large enough $\gamma$. This proves that $a_c$ is stricitly montotone.

Similarly, we can prove the coercivity of $a_C$ by observing that since

we can bound the first term in \eqref{eq:step4_proof}

by using the inequalities in \eqref{eq:C_bounds} as

Thus we have

Since the RHS of \eqref{eq:step5_proof} $\rightarrow \infty$ as $\vert u_h \vert_2 \rightarrow \infty$ for large enough $\gamma$, the coercivity of $a_C$ is also established. The result follows with the application of the Minty-Browder theorem.

□

The above theorem can be extended to higher dimensional settings as well, and provides with with well-posedness of our discretized modified problem. We are now ready to test our new model on the physical scenarios discussed above. First, let us return to the 1D clamped bar scenario as in our previous blog post to highlight the importance of well-posedness as proved above and see if we indeed have any improvement in the robustness of the computational solver.

Note. Perhaps the biggest (if not one of the biggest) problems with the above approach is that the parameter $\gamma$ depends on the grid size $h$. This means that for a fixed $\gamma$, if we keep refining the mesh, we will eventually land an ill-posed problem, where the computational solver will not converge and we can expect the same challenges as the St-Venant Kirchhoff model. This is duly noted, and is an interesting hard issue to solve, but in our regime of numerical results we don’t face much of an issue. We may think of $\gamma$ as a parameter that helps obtain a physically sound displacement profile for large compressive forces.

5. Testing Our New Model: Numerical Results Revisited

5.1. Robustness Testing: 1D Clamped Bar Revisited

As discussed in our previous post, the 1D clamped bar example was not robust with the choice of a different initial guess $U^{(0)}$ or if we refine the grid too much when $f$ is large. We now test these aspects of robustness with our new modified model as in \eqref{eq:modified_law}. We follow the same parameters and tolerances as mentioned in our previous blog post, and test the robustness of the solver with respect to different grid sizes and different initial guesses. The results are shown in Fig. 4 and Fig. 5 below.

Let us now discuss the performance of the model. On the physical side, as expected, the material response is much more stiffer, resulting in smaller displacements. On the computational side, there is good improvement: convergence …

5.2. Incompressible Block With New Model

We are now in a state to provide the numerical results for the incompressible block scenario as above. After some trial and error, we choose the value $\gamma = 4 \times 10^8$ [Pa] for $t_N = 4 \times 10^8$. The results are shown in Fig. XXX below.

Further Reading and Thoughts

The above exposition demonstrates th…

References