Help debugging finite element solution in nonlinear elasticity

Question

I'm writing some code to solve problems in nonlinear elasticity using finite element methods. I have been following Bathe's book but I am having trouble with some nagging details.

My question is related to which configuration variables to use for calculating different quantities. I'm working on static problems in the updated Lagrangian formulation, so for each "time" $t$ I have a loading coefficient $\lambda_t$, and a displacement $\sideset{^t}{}{U}$. In chapter 6, Bathe derives a linearization of the principle of virtual work, resulting in the equation (6.103): \begin{equation} \left(\sideset{^t_{t\,}}{_\text{L}}{K}+\sideset{^t_{t\,}}{_\text{NL}}{K}\right)U=\sideset{^{t+\Delta t}}{}R-\sideset{^t_t}{}F \end{equation} where $\sideset{^t_{t\,}}{_\text{L}}{K}$ is the "linear" part of the incremental stiffness matrix in the configuration with displacement $\sideset{^t}{}{U}$, $\sideset{^t_{t\,}}{_\text{NL}}{K}$ is the "nonlinear" part, $U$ is the displacement increment, $\sideset{^{t+\Delta t}}{}R$ is the externally applied load at "time" $t+\Delta t$, and $\sideset{^t_t}{}F$ is the nodal point force resulting from the stress, at "time" $t$.

My question is about actually implementing this as a system of nonlinear equations to be solved via Newton's method. As I understand it, the idea is to consider a function $H(\sideset{^{t}}{}U)=\sideset{^{t+\Delta t}}{}R-\sideset{^t_t}{}F$, where $\sideset{^t_t}{}F$ is a function of $\sideset{^{t}}{}U$ via the stress.

We want to find a displacement $\sideset{^{t+\Delta t}}{}U$ s.t. the unbalanced force vector is zero: \begin{equation} 0=H(\sideset{^{t+\Delta t}}{}U)\approx H(\sideset{^t}{}U)+\nabla H(\sideset{^t}{}U)\left(\sideset{^{t+\Delta t}}{}U-\sideset{^t}{}U\right)=\sideset{^{t+\Delta t}}{}R-\sideset{^t_t}{}F - \left(\sideset{^t_{t\,}}{_\text{L}}{K}+\sideset{^t_{t\,}}{_\text{NL}}{K}\right)U \end{equation} so that $\left(\sideset{^t_{t\,}}{_\text{L}}{K}+\sideset{^t_{t\,}}{_\text{NL}}{K}\right)$ is minus the Jacobian of the unbalanced force vector $\sideset{^{t+\Delta t}}{}R-\sideset{^t_t}{}F$. Is this interpretation correct? Is the sum of the incremental stiffness matrices derived by Bathe exactly equal to minus the Jacobian of the unbalanced force vector, when all are evaluated in the same configuration? In implementing Newton's method, an iteration is introduced, so that we can write $U^{(k)}$ for the incremental displacement at the $k$th iteration. How does this affect the terms in the system of equations? At iteration $k$, are the integrals still calculated over the volume $\sideset{^t}{}V$, or over $\sideset{^{t+\Delta t}}{}V^{(k-1)}$? What about the stress, and the stiffness matrices?

It seems to me that a full implementation of Newton's method, with no modifications, would involve calculating the unbalanced force vector AND the stiffness matrices in the most recently calculated configuration $\sideset{^{t+\Delta t}}{^{(k-1)}}U$: \begin{equation} 0=H(\sideset{^{t+\Delta t}}{^{(k)}}U)\approx H(\sideset{^{t+\Delta t}}{^{(k-1)}}U)+\nabla H(\sideset{^{t+\Delta t}}{^{(k-1)}}U)\left(\sideset{^{t+\Delta t}}{^{(k)}}U-\sideset{^{t+\Delta t}}{^{(k-1)}}U\right)=\sideset{^{t+\Delta t}}{}R-\sideset{^{t+\Delta t}_{t+\Delta t}\,}{^{(k-1)}}F - \left(\sideset{^{t+\Delta t}_{t+\Delta t\,}}{^{(k-1)}_\text{L}}{K}+\sideset{^{t+\Delta t}_{t+\Delta t\,}}{^{(k-1)}_\text{NL}}{K}\right)U \end{equation} and solving the linear system for the incremental displacement $U$. However, when I set things up this way, the Jacobian of $H$ at displacement $\sideset{^{t+\Delta t}}{^{(k-1)}}U$, as calculated by finite difference, is not equal to $-\left(\sideset{^{t+\Delta t}_{t+\Delta t\,}}{^{(k-1)}_\text{L}}{K}+\sideset{^{t+\Delta t}_{t+\Delta t\,}}{^{(k-1)}_\text{NL}}{K}\right)$.

One clue is that the coded Jacobian matches the finite difference Jacobian almost exactly, for all iterations at $t=0$. My first instinct was to look for errors coming from the "nonlinear" stiffness matrix $K_{\text{NL}}$, but that matrix has entries which are far too small to explain the disparity.

I am fairly confident that I've implemented the calculation properly, provided my understanding is correct. I can provide details if anyone is willing to help me check them, and help me debug the many steps in this calculation - Green-Lagrange strain, deformation gradient, second Piola-Kirchhoff stress, Cauchy stress, etc. Suggestions for sanity checks are much appreciated as well. Thanks!

Answer 1

I find Bathe's notation and presentation confusing. First, it is hard to keep track of all the superscripts and subscripts. And second, writing the equations in an incremental form introduces a confusion between load increments and Newton iteration increments. Here is the way I like to think about the nonlinear FE equations that seems simpler to me.

Assuming we want to solve for the displacements, $U$, of a body under some set of loads we start with Bathe's equation 6.2

$$ ^tR - {^tF} = 0 $$ I'll drop the superscript $t$ for simplicity. Note, also, that I'm defining $U$ as the total displacements of the body, NOT a displacement increment. We assume that $F$, but not $R$, is a function of the displacements so the equation is nonlinear.

We want to solve this equation by Newton iteration as follows

$$ K^{(i-1)}\Delta U^{(i)} = R - F^{(i-1)} $$ $$ U^{(i)} = U^{(i-1)} + \Delta U^{(i)} $$ where $i$ is the Newton iteration and $$ K = \frac{\partial F}{\partial U} $$ This is Bathe's equation 6.8.

One more point concerns the "updated" versus total Lagrangian formulations. Bathe presents the equations using both formulations but states that, theoretically they are equivalent but that the updated Lagrangian may, in some cases, be more "numerically effective." But, as far as I am aware, he gives no specific examples to show this. For simplicity, I suggest using the total Lagrangian formulation where the reference configuration is always the initial configuration. That is, integrals over the element volume required for evaluating the internal forces and tangent stiffness matrices are always the initial, undeformed volume.

Your test of making sure the tangent stiffness matrix calculated in closed form is equal to the one calculated using finite differences of the internal force vector is a good one.