6.1 Newmark Integration
FEBio can perform a nonlinear dynamic analysis by iteratively solving the following nonlinear semi-discrete finite element equations .
\[ \begin{equation} \begin{aligned} & \mathbf{M\ddot{d}}_{n+1}^{k}+\mathbf{K}\Delta\mathbf{d}^{k}=\mathbf{T}_{n+1}^{k}-\mathbf{F}_{n+1}\\ & \mathbf{d}_{n+1}^{k}=\mathbf{d}_{n+1}^{k-1}+\Delta\mathbf{d}^{k} \end{aligned} \,.\label{eq717} \end{equation} \]
Here, \(\mathbf{M}\) is the mass matrix, \(\mathbf{K}\) the stiffness matrix, \(\mathbf{T}\) the internal force (stress) vector and \(\mathbf{F}\) the externally applied loads. The upperscript index \(k\) refers to the iteration number, the subscript \(n\)refers to the time increment. The trapezoidal (or midpoint) rule is used to perform the time integration. This results in the following approximations for the displacement and velocity updates.
\[ \begin{equation} \begin{aligned}\mathbf{d}_{n+1} & =\mathbf{d}_{n}+\frac{h}{2}\left(\mathbf{\dot{d}}_{n}+\mathbf{\dot{d}}_{n+1}\right)\\ \mathbf{\dot{d}}_{n+1} & =\mathbf{\dot{d}}_{n}+\frac{h}{2}\left(\mathbf{\ddot{d}}_{n}+\mathbf{\ddot{d}}_{n+1}\right) \end{aligned} \,.\label{eq718} \end{equation} \]
Using \eqref{eq718} we can solve for \(\mathbf{\ddot{d}}_{n+1}\),
\[ \begin{equation} \mathbf{\ddot{d}}_{n+1}^{k}=\frac{4}{h^{2}}\left(\mathbf{d}_{n+1}^{k-1}-\mathbf{d}_{n}+\Delta\mathbf{d}^{k}\right)-\frac{4}{h}\mathbf{\dot{d}}_{n}-\mathbf{\ddot{d}}_{n}\,.\label{eq719} \end{equation} \]
Substituting this into equation \eqref{eq717} results in the following linear system of equations.
\[ \begin{equation} \left(\frac{4}{h^{2}}\mathbf{M}+\mathbf{K}\right)\Delta\mathbf{d}^{k}=\mathbf{T}_{n+1}^{k}-\mathbf{F}_{n+1}-\mathbf{M}\left(\frac{4}{h^{2}}\left(\mathbf{d}_{n+1}^{k-1}-\mathbf{d}_{n}\right)-\frac{4}{h}\mathbf{\dot{d}}_{n}-\mathbf{\ddot{d}}_{n}\right)\,.\label{eq720} \end{equation} \]
Solving this equation for \(\Delta\mathbf{d}^{k}\) and using \eqref{eq717} gives the new displacement vector \(\mathbf{d}_{n+1}^{k}\). The acceleration vector \(\mathbf{\ddot{d}}_{n+1}^{k}\) can then be found from \eqref{eq719} and the velocity vector \(\mathbf{\dot{d}}_{n+1}^{k}\) from\eqref{eq719}. This algorithm is repeated until convergence is reached.
To solve a differential equation which is second-order in time, we need to perform a numerical integration in the time domain. Let \(\theta\left(t\right)\) denote the function of interest and let \(t_{n}\) and \(t_{n+1}\) represent consecutive time steps such that \(\Delta t=t_{n+1}-t_{n}\). The function \(\theta\left(t\right)\) may be represented at each time point as \(\theta_{n}=\theta\left(t_{n}\right)\) and \(\theta_{n+1}=\theta\left(t_{n+1}\right)\). The Newmark integration formulas are used to evaluate \(\theta_{n+1}\) and \(\dot{\theta}_{n+1}\) at time \(\theta_{n+1}\), assuming that they can be integrated from a judiciously selected \(\ddot{\theta}\left(t_{n+\gamma}\right)\) in the time interval \(\left[t_{n},t_{n+1}\right]\). Using the mean value theorem for definite integrals, we know that an exact solution may be found for the integral according to
\[ \begin{equation} \int_{t_{n}}^{t_{n+1}}\ddot{\theta}\left(t\right)\,dt=\dot{\theta}_{n+1}-\dot{\theta}_{n}\equiv\ddot{\theta}\left(t_{n+\gamma}\right)\Delta t\,,\label{eq:Newmark-mean-value-theorem} \end{equation} \]
where \(\gamma\) is generally unknown a priori. In the Newmark integration scheme we let
\[ \begin{equation} \ddot{\theta}\left(t_{n+\gamma}\right)=\gamma\ddot{\theta}_{n+1}+\left(1-\gamma\right)\ddot{\theta}_{n}\,,\label{eq:Newmark-gamma-deriv} \end{equation} \]
where \(\gamma\) is a user-selected parameter in the range \(0\) to \(1\). It follows that
\[ \begin{equation} \dot{\theta}_{n+1}=\dot{\theta}_{n}+\Delta t\left[\gamma\ddot{\theta}_{n+1}+\left(1-\gamma\right)\ddot{\theta}_{n}\right]\,.\label{eq:Newmark-velocity} \end{equation} \]
We can similarly integrate this function twice to obtain \(\theta_{n+1}\),
\[ \int_{t_{n}}^{t_{n+1}}\int_{t_{n}}^{t}\ddot{\theta}\left(\tau\right)\,d\tau\,dt=\int_{t_{n}}^{t_{n+1}}\left(\dot{\theta}\left(t\right)-\dot{\theta}_{n}\right)\,dt=\theta_{n+1}-\theta_{n}-\dot{\theta}_{n}\Delta t\equiv\ddot{\theta}\left(t_{n+2\beta}\right)\frac{\Delta t^{2}}{2} \]
where we let
\[ \ddot{\theta}\left(t_{n+2\beta}\right)=2\beta\ddot{\theta}_{n+1}+\left(1-2\beta\right)\ddot{\theta}_{n}\,. \]
Here, \(\beta\) represents a parameter that varies from \(0\) to \(\frac{1}{2}\). It follows that
\[ \begin{equation} \theta_{n+1}=\theta_{n}+\dot{\theta}_{n}\Delta t+\frac{\Delta t^{2}}{2}\left[2\beta\ddot{\theta}_{n+1}+\left(1-2\beta\right)\ddot{\theta}_{n}\right]\,,\label{eq:Newmark-displacement} \end{equation} \]
or alternatively,
\[ \begin{equation} \ddot{\theta}_{n+1}=\frac{1}{\beta\Delta t^{2}}\left(\theta_{n+1}-\theta_{n}-\dot{\theta}_{n}\Delta t\right)-\left(\frac{1}{2\beta}-1\right)\ddot{\theta}_{n}\,,\label{eq:Newmark-acceleration} \end{equation} \]
from which we may re-evaluate \eqref{eq:Newmark-velocity} as
\[ \begin{equation} \dot{\theta}_{n+1}=\dot{\theta}_{n}+\Delta t\left[\frac{\gamma}{\beta\Delta t^{2}}\left(\theta_{n+1}-\theta_{n}-\dot{\theta}_{n}\Delta t\right)+\left(1-\frac{\gamma}{2\beta}\right)\ddot{\theta}_{n}\right]\,.\label{eq:Newmark-velocity-redux} \end{equation} \]
Stability of this integration scheme is guaranteed when
\[ \begin{equation} \gamma\ge\frac{1}{2}\,,\quad\beta\ge\frac{\left(\gamma+\frac{1}{2}\right)^{2}}{4}\,.\label{eq:Newmark-stability} \end{equation} \]