Newton-Raphson method

In general, the equilibrium equations obtained by discretising the virtual work equations can be written at time \(t\) symbolically as:

\[ \boldsymbol{F}_\alpha^t(\boldsymbol{d}_\beta^t) = 0 \]

Where \(\boldsymbol{F}^t_\alpha\) is the force component conjugate to the \(\alpha\)th variable in the problem and \(\boldsymbol{d}^t_\beta\) is the value of the \(\beta\)th variable. In the following the subscripts \(\alpha\) and \(\beta\) are dropped for the sake of simplicity. Focusing on history-depending problems, the solution must be developed by a series of small increments \(\Delta \sqcup\):

\[ \boldsymbol{d}^t = \boldsymbol{d}^{(t_0+\Delta t)} = \boldsymbol{d}^{t_0} + \Delta \boldsymbol{d} \]

The total solution increment \(\Delta \boldsymbol{d}\) is adapted iteratively by iterative increments \(\boldsymbol{c}^i\) until equilibrium is reached, up to a prescribed tolerance. Indicating the iteration number \(i\) with a right subscript, the incremental solution at iteration \(i + 1\) is calculated from

\[ \Delta \boldsymbol{d}^{i+1} = \Delta \boldsymbol{d}^i + \boldsymbol{c}^{i+1} \]

To determine the iterative increments \(\boldsymbol{c}^{i+1}\) numgeo uses the Newton's method as a numerical technique for solving the equilibrium equations.

The system of equations soe to be solved can be summarized as follows:

\[ \boldsymbol{r}\left(\boldsymbol{d},\dot{\boldsymbol{d}},\ddot{\boldsymbol{d}}\right) = s_D \boldsymbol{M}\ddot{\boldsymbol{d}} - \boldsymbol{F}\left(\boldsymbol{d},\dot{\boldsymbol{d}}\right) \label{eq:solution-residual} \]

Therein \(\boldsymbol{r}=[\boldsymbol{r}^s,\boldsymbol{r}^w,\boldsymbol{r}^a]\) is the residuum, \(\boldsymbol{d}\) are the unknown nodal variables, \(\boldsymbol{M}\ddot{\boldsymbol{d}}\) corresponds to the inertia forces and \(\boldsymbol{F}(\boldsymbol{d},\dot{\boldsymbol{d}})\) gathers the internal and external nodal forces. For dynamic simulations \(s_D=1\) holds, whereas \(s_D=0\) for (quasi-)static simulations.

Time integration of the spatially discretized governing equations is one of the important steps in numerical analysis of dynamic problems for achieving accurate results and simultaneously saving computational time. To integrate the spatially discretized governing equations, the improved Newmark's method 1 by Hilber et al. (1977)2 is used. The so-called Hilber-Hughes-Taylor-Method (HHT) incorporates the parameter \(\alpha\) to damp out any spurious participation of high-frequency modes. The smaller the value of the parameter \(\alpha\), the greater the numerical dissipation. Eq. (\(\ref{eq:solution-residual}\)) is modified as follows:

\[ \boldsymbol{r}\left(\boldsymbol{d},\dot{\boldsymbol{d}},\ddot{\boldsymbol{d}}\right) = s_D \boldsymbol{M}\ddot{\boldsymbol{d}} - (1+\alpha)\boldsymbol{F}\left(\boldsymbol{d},\dot{\boldsymbol{d}}\right) - \alpha \boldsymbol{F}\left(\boldsymbol{d}^{(t)},\ddot{\boldsymbol{d}}^{(t)}\right) \]

with:

\[ \begin{align} \dot{\boldsymbol{d}} = & \frac{\gamma}{\beta \Delta t}\left(\boldsymbol{d}-\boldsymbol{d}^{(t)}\right)+\frac{\beta-\gamma}{\beta} \dot{\boldsymbol{d}}^{(t)} + \frac{\Delta t (\beta - 1/2 \gamma)}{\beta} \ddot{\boldsymbol{d}}^{(t)} \\ \ddot{\boldsymbol{d}} = & \frac{\gamma}{\beta \Delta t^2} \left(\boldsymbol{d} - \boldsymbol{d}^{(t)}\right) + \frac{1}{\beta \Delta t} \dot{\boldsymbol{d}}^{(t)} + \frac{\beta-1/2}{\beta} \ddot{\boldsymbol{d}}^{(t)} \end{align} \]

Variables with the index \((t)\) correspond to the known (converged) state at time \(t\), where the index is omitted the variable belongs to the (unknown) state at the end of the current increment \(t+\Delta t\). The parameters of numerical dissipation are selected such that \(-1/3 \leq \alpha \leq 0\), \(\beta = (1-\alpha^2)/4\) and \(\gamma = 1/2-\alpha\). For \(\alpha = 0\) (\(\beta=1/4\) and \(\gamma=1/2\)) the Newmark's method is recovered. Notice that only then a second order accuracy is achieved in linear problems, otherwise the method has only first order accuracy.

Additionally, the Backward Euler integration scheme is obtained when using Eq. (\(\ref{eq:backward-euler}\)) for the update of \(\dot{\boldsymbol{d}}\) and \(\ddot{\boldsymbol{d}}\). \(\alpha=0\) is implied (solving of Eq. (\(\ref{eq:solution-residual}\)) at the end of the time step \(t+\Delta t\)).

\[ \dot{\boldsymbol{d}} = \frac{\boldsymbol{d}-\boldsymbol{d}^{(t)}}{\Delta t} ~~~ ; ~~~ \ddot{\boldsymbol{d}} = \frac{\dot{\boldsymbol{d}}-\dot{\boldsymbol{d}}^{(t)}}{\Delta t} \label{eq:backward-euler} \]

The nonlinear system of equations is solved using the iterative Newton-Raphson-Method, where the estimated solution \(\boldsymbol{d}^{(i)}\) at the beginning of the iteration \((i)\) is enhanced with each iteration \((i+1)\) by adding a correction \(\boldsymbol{c}^{(i+1)}\). The nonlinear system at time \(t+\Delta t\) is linearised using the first two terms of a Taylor series:

\[ \frac{\partial \boldsymbol{r}\left(\boldsymbol{d}^{(i)}\right)}{\partial \boldsymbol{d}}\boldsymbol{c}^{(i+1)} = - \boldsymbol{r}\left(\boldsymbol{d}^{(i)}\right) \label{eq:eq1} \]

\(\dfrac{\boldsymbol{r}\left(\boldsymbol{d}^{(i)}\right)}{\partial \boldsymbol{d}}\) is the Jacobian of the system of equations. Convergence is reached when Eq. \(\ref{eq:eq1}\) and \(\boldsymbol{d}^{(i+1)} = \boldsymbol{d}^{(i)}\) is satisfied within a given tolerance. Details on the different criteria used in numgeo to judge the convergence of the simulation are presented in the next Section

  1. N. M. Newmark. A method of computation for structural dynamics. Journal of the engineering mechanics division, 85(3):67–94, 1959. ↩

  2. Hans M. Hilber, Thomas J. R. Hughes, and Robert L. Taylor. Improved numerical dissipation for time integration algorithms in structural dynamics. Earthquake Engineering & Structural Dynamics, 5(3):283–292, 7 1977. doi:10.1002/eqe.4290050306. ↩