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!