Skip to content

Newton-Raphson Method

Linearization and Iterative Recurrence

Virtual Work of External Forces and Assembly of the Global Equation gives a nonlinear equation for the nodal displacement at time \(t_{n+1}\), \(\boldsymbol{u}_{n+1}\), which is solved by the Newton-Raphson method. The nodal displacement up to time \(t_n\), \(\boldsymbol{u}_n\), is assumed known, and the displacement increment \(\Delta\boldsymbol{u}\) is taken as the unknown variable to determine

\[ \boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u} \]

Hereafter, dependence of the external-force vector on nodal displacement is neglected, and with \(\boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{F}_{n+1}\),

\[ \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

is solved.

At the current solution \(\Delta\boldsymbol{u}\), define the tangent stiffness

\[ \boldsymbol{K} = \left. \frac{\partial \boldsymbol{Q}}{\partial \boldsymbol{u}} \right|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}} \]

Using this, linearization of the nonlinear equation gives

\[ \boldsymbol{K}\, d\boldsymbol{u} + \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) - \boldsymbol{F}_{n+1} = \boldsymbol{0} \]

Let the correction in the \(i\)-th iteration be \(d\boldsymbol{u}_i\), and let the residual vector at the start of the iteration be

\[ \boldsymbol{R}_{i-1} = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u}) \]

Then the iterative recurrence is

\[ \boldsymbol{K}_i\, d\boldsymbol{u}_i = \boldsymbol{R}_{i-1}, \qquad \Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i \]

Thus, the residual \(\boldsymbol{R}_i\) is a quantity corresponding to the force imbalance from equilibrium.

Construction of the Tangent Stiffness Matrix

The tangent stiffness \(\boldsymbol{K} = \partial\boldsymbol{Q}/\partial\boldsymbol{u}\) is constructed by partially differentiating the element internal-force vector obtained in Discretization of the Virtual Work of Internal Forces with respect to nodal displacement, integrating the resulting element-level integrands over each element domain, and assembling them. Denoting the element-level integrand by \(\boldsymbol{K}^e_X\) (reference-configuration notation, TL formulation) or \(\boldsymbol{K}^e_x\) (current-configuration notation, UL formulation), the element tangent stiffness is

\[ \boldsymbol{K}^e = \int_{\Omega^e_0} \boldsymbol{K}^e_X\, dV \quad (\text{TL}), \qquad \boldsymbol{K}^e = \int_{\Omega^e} \boldsymbol{K}^e_x\, dv \quad (\text{UL}) \]

The following gives the final forms of the TL/UL integrands. In both cases they are decomposed into the sum of a material stiffness term (initial-displacement term) and a geometric stiffness term (initial-stress term).

Total Lagrange Formulation

In the Total Lagrange formulation, a linear relation is assumed between the rate of the second Piola-Kirchhoff stress \(\dot{\boldsymbol{S}}\) and the Green-Lagrange strain rate \(\dot{\boldsymbol{E}}\), namely \(\dot{\boldsymbol{S}} = \boldsymbol{\mathsf{C}}:\dot{\boldsymbol{E}}\). This corresponds to constitutive laws for linear elastic materials (St. Venant-Kirchhoff materials) and hyperelastic materials, and FrontISTR uses the Total Lagrange formulation for these materials. The element tangent-stiffness integrand is then written in tensor form as

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_X\, \dot{\boldsymbol{u}}^e = \dot{\boldsymbol{S}}:\delta\boldsymbol{E} + \boldsymbol{S}:(\delta\boldsymbol{F}^T \dot{\boldsymbol{F}}) \]

The first term on the right-hand side is the material stiffness term (initial-displacement term), and the second term is the geometric stiffness term (initial-stress term).

In the FrontISTR implementation, this integrand is evaluated in matrix form using Voigt notation:

\[ \boldsymbol{K}^e_X = (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T\, \tilde{\boldsymbol{C}}\, (\boldsymbol{B}_L + \boldsymbol{B}_{NL}) + \boldsymbol{F}_9^T\, \boldsymbol{S}_9\, \boldsymbol{F}_9 \]

Each matrix is as follows. \(\boldsymbol{B}_L, \boldsymbol{B}_{NL}\) are the B matrices introduced in Discretization of the Virtual Work of Internal Forces, and \(\tilde{\boldsymbol{C}}\) is the Voigt representation of the constitutive tensor \(\boldsymbol{\mathsf{C}}\), namely a \(6\times 6\) material stiffness matrix (Tensor Notation and Mathematical Fundamentals). \(\boldsymbol{S}_9, \boldsymbol{F}_9\) are the following rearrangement matrices used to express the geometric stiffness term as a matrix product. First, for a \(3\times 3\) second-order tensor \(\boldsymbol{A}\), define the notation \([\,\cdot\,]\) that rearranges it into a 9-component vector as

\[ [\boldsymbol{A}] = (A_{11}, A_{21}, A_{31}, A_{12}, A_{22}, A_{32}, A_{13}, A_{23}, A_{33})^T \]

With this definition, \(\boldsymbol{F}_9\) expresses the variation of the deformation gradient in the form \([\delta\boldsymbol{F}] = \boldsymbol{F}_9\, \delta\boldsymbol{u}^e\) and is a \(9\times d n_e\) matrix. For element node \(\alpha = 1, \ldots, n_e\), the corresponding \(9\times d\) block is

\[ [\boldsymbol{F}_9]_\alpha = \begin{bmatrix} (\partial N_\alpha^e/\partial X_1)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_2)\, \boldsymbol{I} \\ (\partial N_\alpha^e/\partial X_3)\, \boldsymbol{I} \end{bmatrix} \qquad (\boldsymbol{I} \text{ is the } 3\times 3 \text{ identity matrix}) \]

and is given by \(\boldsymbol{F}_9 = [[\boldsymbol{F}_9]_1, \ldots, [\boldsymbol{F}_9]_{n_e}]\), with the blocks arranged horizontally in element-node order. \(\boldsymbol{S}_9\) is selected so that, when combined with this matrix, the geometric stiffness term is expressed as \(\delta\boldsymbol{u}^{eT}\, \boldsymbol{F}_9^T \boldsymbol{S}_9 \boldsymbol{F}_9\, \dot{\boldsymbol{u}}^e\); it is the following \(9\times 9\) matrix

\[ \boldsymbol{S}_9 = \begin{bmatrix} S_{11} \boldsymbol{I} & S_{12} \boldsymbol{I} & S_{13} \boldsymbol{I} \\ S_{21} \boldsymbol{I} & S_{22} \boldsymbol{I} & S_{23} \boldsymbol{I} \\ S_{31} \boldsymbol{I} & S_{32} \boldsymbol{I} & S_{33} \boldsymbol{I} \end{bmatrix} \]

This is the resulting matrix.

Updated Lagrange Formulation

In the Updated Lagrange formulation, a linear relation is assumed between the Jaumann rate of the relative Kirchhoff stress tensor \(\hat{\boldsymbol{\sigma}}^{\nabla J}\) and the rate-of-deformation tensor \(\boldsymbol{D}\), namely \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\). This is the form of a hypoelastic constitutive law common to linear elastic, elastoplastic, and creep materials, and FrontISTR uses the Updated Lagrange formulation for these materials. The element tangent-stiffness integrand expressed in the current configuration is then written in tensor form as

\[ \delta\boldsymbol{u}^{eT}\, \boldsymbol{K}^e_x\, \dot{\boldsymbol{u}}^e = \boldsymbol{\sigma}^{\nabla T}:\delta\boldsymbol{A}_{(L)} + \boldsymbol{\sigma}:(\delta\boldsymbol{F}_t^T\, \boldsymbol{L}) \]

where \(\boldsymbol{\sigma}^{\nabla T}\) is the Truesdell rate, \(\boldsymbol{A}_{(L)}\) is the linear part of the Almansi strain, \(\boldsymbol{F}_t = \partial\boldsymbol{u}/\partial\boldsymbol{x}\) is the displacement gradient with respect to the current configuration, and \(\boldsymbol{L}\) is the velocity-gradient tensor. The first term on the right-hand side is the material stiffness term, and the second term is the geometric stiffness term.

In the FrontISTR implementation, this integrand is evaluated in matrix form using Voigt notation:

\[ \boldsymbol{K}^e_x = \boldsymbol{b}^T\, (\tilde{\boldsymbol{C}} - \boldsymbol{G})\, \boldsymbol{b} + \boldsymbol{f}_9^T\, \boldsymbol{\sigma}_9\, \boldsymbol{f}_9 \]

Here, \(\boldsymbol{b}\) is the B matrix constructed in the current configuration (Discretization of the Virtual Work of Internal Forces). \(\boldsymbol{\sigma}_9, \boldsymbol{f}_9\) are obtained from \(\boldsymbol{S}_9, \boldsymbol{F}_9\) defined for the TL formulation by replacing the second PK stress \(\boldsymbol{S}\) with the Cauchy stress \(\boldsymbol{\sigma}\) and the reference-configuration gradient \(\partial N_\alpha^e/\partial X_i\) with the current-configuration gradient \(\partial N_\alpha^e/\partial x_i\).

\(\boldsymbol{G}\) is a Cauchy-stress-dependent correction matrix required to make the hypoelastic constitutive law \(\hat{\boldsymbol{\sigma}}^{\nabla J} = \boldsymbol{\mathsf{C}}:\boldsymbol{D}\) consistent with the tangent-stiffness framework as a Truesdell-rate-based constitutive law. It is obtained by arranging the fourth-order tensor components \(G_{ijkl} = \delta_{il}\sigma_{kj} + \delta_{jl}\sigma_{ik}\) in \(6\times 6\) Voigt form as

\[ \boldsymbol{G} = \begin{bmatrix} 2\sigma_{11} & 0 & 0 & \sigma_{12} & 0 & \sigma_{31} \\ 0 & 2\sigma_{22} & 0 & \sigma_{12} & \sigma_{23} & 0 \\ 0 & 0 & 2\sigma_{33} & 0 & \sigma_{23} & \sigma_{31} \\ \sigma_{12} & \sigma_{12} & 0 & \tfrac{\sigma_{11}+\sigma_{22}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{23}}{2} \\ 0 & \sigma_{23} & \sigma_{23} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{22}+\sigma_{33}}{2} & \tfrac{\sigma_{12}}{2} \\ \sigma_{31} & 0 & \sigma_{31} & \tfrac{\sigma_{23}}{2} & \tfrac{\sigma_{12}}{2} & \tfrac{\sigma_{33}+\sigma_{11}}{2} \end{bmatrix} \]

This is the resulting matrix.

Assembly of the Global Stiffness Matrix

The global tangent stiffness \(\boldsymbol{K}\) is obtained by dividing each element stiffness \(\boldsymbol{K}^e\) into \(d\times d\) blocks \(\boldsymbol{K}^e_{\alpha\beta}\) for each pair of nodes and using the second-order-tensor assembly set \(\mathcal{E}^2(i_g, i_h)\) introduced in Assembly of Element-Nodal Physical Quantities:

\[ \boldsymbol{K}_{i_gi_h} = \sum_{(e,\alpha,\beta) \in \mathcal{E}^2(i_g, i_h)} \boldsymbol{K}^e_{\alpha\beta} \]

The resulting values are arranged as a matrix with row \(i_g\) and column \(i_h\). In the implementation, the set \(\mathcal{E}^2\) is not constructed explicitly; instead, the corresponding blocks are added directly within the element loop. The matrix is square with dimension equal to degrees of freedom per node \(\times\) total number of nodes \(n_g\), but because components other than those between nodes connected through elements are \(0\), it is stored in sparse-matrix form.

The element stiffness matrices for the TL and UL formulations have the same form except for switching the reference configuration (nodal coordinates and the source used to construct the B matrix) and the presence or absence of the \(\boldsymbol{G}\) matrix. FrontISTR therefore implements both formulations in a common subroutine.

Iteration Algorithm

Summarizing the above, at the start of the iteration set \(\Delta\boldsymbol{u} = \boldsymbol{0}\) and compute the initial residual \(\boldsymbol{R}_0 = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n)\). Then, in the \(i\)-th iteration, perform the following procedure.

  1. At the current displacement \(\boldsymbol{u}_n + \Delta\boldsymbol{u}\), compute the tangent stiffness \(\boldsymbol{K}_i = \partial\boldsymbol{Q}/\partial\boldsymbol{u}|_{\boldsymbol{u}_n + \Delta\boldsymbol{u}}\) using the procedure in Construction of the Tangent Stiffness Matrix.
  2. To impose geometric boundary conditions, modify the tangent stiffness matrix and residual vector for the degrees of freedom subject to displacement constraints, obtaining \(\tilde{\boldsymbol{K}}_i, \tilde{\boldsymbol{R}}_{i-1}\) (Treatment of Geometric Boundary Conditions).
  3. Solve the linear equation \(\tilde{\boldsymbol{K}}_i\, d\boldsymbol{u}_i = \tilde{\boldsymbol{R}}_{i-1}\) to obtain the correction \(d\boldsymbol{u}_i\). This procedure often accounts for most of the computational cost of the iterative calculation.
  4. Update the displacement increment as \(\Delta\boldsymbol{u} \leftarrow \Delta\boldsymbol{u} + d\boldsymbol{u}_i\), and accordingly compute the internal-force vector \(\boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\) and residual \(\boldsymbol{R}_i = \boldsymbol{F}_{n+1} - \boldsymbol{Q}(\boldsymbol{u}_n + \Delta\boldsymbol{u})\).
  5. Check convergence and terminate the iteration if convergence is achieved. Components corresponding to constrained reactions appear in the residual \(\boldsymbol{R}_i\) at degrees of freedom subject to geometric boundary conditions, so the convergence indicator is constructed from \(\tilde{\boldsymbol{R}}_i\) after excluding these components. Specific convergence indicators and thresholds are described in Convergence Criteria. If convergence is not achieved and the iteration limit is reached, the iteration is treated as failed.

When the iteration converges, add the converged \(\Delta\boldsymbol{u}\) to \(\boldsymbol{u}_n\) to obtain the accumulated displacement at time \(t_{n+1}\), \(\boldsymbol{u}_{n+1} = \boldsymbol{u}_n + \Delta\boldsymbol{u}\), and proceed to the next time step.