1. Introduction

In this blog post, 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 scenario; 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 not orientation preserving, i.e., a non-physical solution, and hence is not the best model to simulate the block compression scenario. Indeed, the St-Venant Kirchhoff model is usually best to simulate scenarios involving small strains ¹. The convergence issue is further aggravated as the model for finer grid sizes, which are needed to avoid the effect of locking when using low order finite elements.

In this post, we are interested to see how far we can push the St-Venant Kirchhoff model, and if there is a way to make modifications which can help us better approximate and capture displacement profiles under large forces and for fine grids. To this end, we focus on building a new St-Venant Kirchhoff type model by adding a linear term which can help us approximate large displacements while avoiding the above challenges. Our motivation comes from an apriori knowledge of the lack of well-posedness of the system and we introduce physical “stiffness” which helps us solve this issue to some extent and also allows the new model to capture displacement profiles under large forces and on finer grids.

2. Details of the Computational Solver

Following the notation we set up in our previous blog 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{eq: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_2} 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 finer grid sizes, or (ii) convergence to different solution profiles for different initial guesses in the nonlinear solver (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 the low order 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. Or simply, one may use Q2 elements which, in our experience as well, reduces locking. In our experiments, however, since we work with finer grids to test the robustness of our model, we do not encounter locking to a high degree when using 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 clamped bar under a constant force scenario from our previous blog post. Extending into 2D, we now consider $\Omega = (0, 1) \times (0, 0.4)$ [m $^2$]. The boundary conditions are

The Young’s modulus is $E_Y = 10^7$ [Pa] and the Poisson ratio is $\nu = 0.48$ [-]. For an external force of $f = 4 \times 10^7$ [N / m $^3$], our computational solver does not converge, and hence we choose a smaller force of $f = 3 \times 10^6$. The results are shown in Fig. 2 for a uniform grid size of $0.025 \times 0.025$ [m $^2$] (i.e. $h = 0.025$).

The profile shows the bar deformed with more mass accumulating towards the right boundary. We also tabulate the performance of the computational solver in Tab. 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$, i.e., for a fine grid. 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.

We note that in both the examples above our model and solver struggles to converge for large $f$ or $t_N$, values which may have led to convergence in 1D. Moreover, even for small $f$ or $t_N$, if we refine the grid sufficiently we reach a stage where we either get a non-physical solution, or our solver does not converge, behaviour that was also observed in 1D in our previous blog post discussion. We now work on addressing this issue by introducing a new parameter in the St-Venant Kirchhoff model.

4. Introducing a New Model

The above examples (and the discussion in our previous blog post) 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 non-physical profile (such as where det$(F) < 0$), issues which we now try to address through the lack of well-posedness of the St-Venant Kirchhoff model. We now work towards improving the model through a stabilizing term to see if we can attempt to have well-posedness in some form.

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 for large displacement valuesr 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.

We now proceed towards obtaining a numerical solution to the system arising from \eqref{eq:modified_law}. 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$.

Existence of solution to Newton’s update. We note 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.

Well-posedness in 1D. We now provide an outline of a proof of well-posedness for the new model. Let us consider the variational form \eqref{eq:modified_variational_form} with homogeneous Dirichlet boundary conditions for simplicity and clarity of exposition. We also now denote by $V_h \subset H_0^1(\Omega)$ the finite dimensional subspace of piecewise-linear functions which vanish on $\partial \Omega$. Rewriting \eqref{eq:modified_variational_form} using a nonlinear operator $a : V_h \rightarrow V_h’$, where $V_h’$ is the dual of $V_h$, the variational form becomes: we seek $u_h \in V_h$ such that

where

and $S \in V_h’$ is defined as

We now investigate 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 the Minty-Browder’s theorem [Ciarlet’ 2013, Theorem 9.14-1] ⁸. That, however, is not a trivial task, and instead we will prove the properties in 1D while also making 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. The cut-off operator is used in literature in other problems as well: see for example ⁹ for an application in thermo-poroelastic settings.

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_var_form} 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 us with well-posedness of our discretized modified problem.

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, and a parameter that we may fine tune if we know apriori the grid scales that we are dealing with.

5. Testing Our New Model: Numerical Results Revisited

We are now ready to test our new model on the physical scenarios discussed above. We begin by returning 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 solver.

5.1. 1D Clamped Bar Revisited

We now study the simple 1D clamped bar example to understand the performance of our new model compared with the St-Venant Kirchhoff model. We follow the same parameters and tolerances as mentioned in our previous blog post, and consider different grid sizes and initial guesses.

Scenario results with new model. We first simulate the scenario using grid size $h = 0.05$ [m] and by choosing $\gamma = 10^7$ [Pa]. The results are shown in Fig. 4.

Robustness testing. As discussed in our previous blog 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}. The results are shown in Fig. 5 and Fig. 6 below.

Performance summary. Let us now discuss the performance of the model and the computational solver. On the physical side, as expected, the material response is much more stiffer, resulting in smaller displacements for the new model. On the computational side, there is good improvement, and the results are tabulated in Tab. 2. In particular, we observe convergence of the solver for finer grids and different initial guesses, behaviour that was not observed when using the St-Venant Kirchhoff model.

We also note that for smaller $\gamma$, the difference between the St-Venant Kirchhoff model and our new model is less (as expected), and consequently convergence issues also arise.

2D results. We now revisit the example above in 2D with $f = 4 \times 10^7$ and $\gamma = 1.55 \times 10^7$. The displacment profile is shown in Fig. 7 for a grid size of $0.0125 \times 0.0125$ [m $^2$] (which led to no-convergence for the St-Venant Kirchhoff model in our discussion above).

Unlike before, the solver now converges in NS iter. 11 and LS iter. 4653, and the maximum displacement magnitude is more than twice what it was in Fig. 2.

5.2. Incompressible Block With New Model

We are now in a position to provide the numerical results for the incompressible block scenario as above. We use $h = 0.05$ [mm] and $\gamma = 4 \times 10^8$ [Pa] (after some trial and error) for $t_N = 4 \times 10^8$. The results are shown in Fig. 8 below.

Fig. 8 (left) shows a much better displacement profile with a maximum value of 4.46 [mm], better than the values reported previously in Fig. 3, although still short of the benchmark values reported in ⁷. The NS iter. in this case are 13 and the LS iter. are 78411.

Sensitivity to $\gamma$. The model is expectedly sensitive to $\gamma$. A small $\gamma$ means we tend towards the St-Venant Kirchhoff model, and a large $\gamma$ means the material responds very stiffly. We attempt to increase the maximum displacement in the incompressible block scenario by varying $\gamma$. Indeed, after some more trial and error, it was observed that when $\gamma = 3.819 \times 10^8$, the maximum displacement improves to 4.89 [mm]; the profile is shown in Fig. 7 (right). In this case, the NS iter. taken are 21 and the LS iter. is 143151.

To conclude the above numerical results, we note that we are able to obtain physical displacement profiles for grid sizes and forces that the St-Venant Kirchhoff model could provide for. This is consistent with the well-posedness that we can hope to achieve with our new model for a large enough $\gamma$.

Further Reading and Thoughts

The above discussion provides an attempt at improving the physical and computational aspect of the St-Venant Kirchhoff model under large compressive forces. Although the new model introduced above is not found on any experimental data, we successfully demonstrate how the modification can help us achieve better results with improved performance, particularly by obtaining agreeable physical displacement profiles under large forces as opposed to the St-Venant Kirchhoff model which does not converge under the same forces. As highlighted, however, it should be noted that the new model does not inherently address the fundamental issue of the St-Venant Kirchhoff model which is that the stored energy function $W(F) \rightarrow 0$ as det$(F) \rightarrow 0$. Indeed, the stored energy of the new model is

which, too, does not provide any resistance to det$(F) \rightarrow 0$ and will face the same issues as the St-Kirchhoff Venant model for large forces.

Moreover, a major drawback is that the parameter $\gamma$ depends on the space $V_h$, i.e., on the grid size $h$ for model and performance improvements, and we also cannot hope that $\gamma$ is dimension invariant: for example, if a particular $\gamma$ value gives good results in 2D, it may still require fine-tuning in 3D. Another issue is that for large $\gamma$ and for small forces or traction, the material response is too stiff which would result in very small and inaccurate displacements. To address this, a more versatile parameter could be a monotnone $\gamma = \gamma(u)$ which may still be able to provide coercivity and monotonicty of the operator associated with $\widetilde{T}$ (provided we may impose something like $\gamma(0) = 0$ and Lipschitz continuity), and would slowly introduce stiffness as $u$ increases with increasing $f$ or $t_N$. That, however, would be a topic for a future blog post.

References