14.12. Newton-Raphson Procedure

The following Newton-Raphson procedure topics are available:

14.12.1. Overview

The finite element discretization process yields a set of simultaneous equations:

 (14–146)

where:

 [K] = coefficient matrix {u} = vector of unknown DOF (degree of freedom) values {Fa} = vector of applied loads

If the coefficient matrix [K] is itself a function of the unknown DOF values (or their derivatives) then Equation 14–146 is a nonlinear equation. The Newton-Raphson method is an iterative process of solving the nonlinear equations and can be written as (Bathe([2])):

 (14–147)
 (14–148)

where:

 = Jacobian matrix (tangent matrix) i = subscript representing the current equilibrium iteration = vector of restoring loads corresponding to the element internal loads

Both and are evaluated based on the values given by {ui}. The right-hand side of Equation 14–147 is the residual or out-of-balance load vector; i.e., the amount the system is out of equilibrium. A single solution iteration is depicted graphically in Figure 14.9: Newton-Raphson Solution - One Iteration for a one DOF model. In a structural analysis, is the tangent stiffness matrix, {ui} is the displacement vector and is the restoring force vector calculated from the element stresses. In a thermal analysis, is the conductivity matrix, {ui} is the temperature vector and is the resisting load vector calculated from the element heat flows. In an electromagnetic analysis, is the Dirichlet matrix, {ui} is the magnetic potential vector, and is the resisting load vector calculated from element magnetic fluxes. In a transient analysis, is the effective coefficient matrix and is the effective applied load vector which includes the inertia and damping effects.

As seen in the following figures, more than one Newton-Raphson iteration is needed to obtain a converged solution. The general algorithm proceeds as follows:

1. Assume {u0}. {u0} is usually the converged solution from the previous time step. On the first time step, {u0} = {0}.

2. Compute the updated tangent matrix and the restoring load from configuration {ui}.

3. Calculate {Δui} from Equation 14–147.

4. Add {Δui} to {ui} in order to obtain the next approximation {ui + 1} (Equation 14–148).

5. Repeat steps 2 to 4 until convergence is obtained.

Figure 14.10: Newton-Raphson Solution - Next Iteration shows the solution of the next iteration (i + 1) of the example from Figure 14.9: Newton-Raphson Solution - One Iteration. The subsequent iterations would proceed in a similar manner.

The solution obtained at the end of the iteration process would correspond to load level {Fa}. The final converged solution would be in equilibrium, such that the restoring load vector (computed from the current stress state, heat flows, etc.) would equal the applied load vector {Fa} (or at least to within some tolerance). None of the intermediate solutions would be in equilibrium.

If the analysis included path-dependent nonlinearities (such as plasticity), then the solution process requires that some intermediate steps be in equilibrium in order to correctly follow the load path. This is accomplished effectively by specifying a step-by-step incremental analysis; i.e., the final load vector {Fa} is reached by applying the load in increments and performing the Newton-Raphson iterations at each step:

 (14–149)

where:

 [Kn,i] = tangent matrix for time step n, iteration i = total applied force vector at time step n = restoring force vector for time step n, iteration i

This process is the incremental Newton-Raphson procedure and is shown in Figure 14.11: Incremental Newton-Raphson Procedure. The Newton-Raphson procedure guarantees convergence if and only if the solution at any iteration {ui} is “near” the exact solution. Therefore, even without a path-dependent nonlinearity, the incremental approach (i.e., applying the loads in increments) is sometimes required in order to obtain a solution corresponding to the final load level.

When the stiffness matrix is updated every iteration (as indicated in Equation 14–147 and Equation 14–149) the process is termed a full Newton-Raphson solution procedure ( NROPT,FULL or NROPT,UNSYM). Alternatively, the stiffness matrix could be updated less frequently using the modified Newton-Raphson procedure (NROPT,MODI). Specifically, for static or transient analyses, it would be updated only during the first or second iteration of each substep, respectively. Use of the initial-stiffness procedure (NROPT,INIT) prevents any updating of the stiffness matrix, as shown in Figure 14.12: Initial-Stiffness Newton-Raphson. If a multistatus element is in the model, however, it would be updated at iteration in which it changes status, irrespective of the Newton-Raphson option. The modified and initial-stiffness Newton-Raphson procedures converge more slowly than the full Newton-Raphson procedure, but they require fewer matrix reformulations and inversions. A few elements form an approximate tangent matrix so that the convergence characteristics are somewhat different.

14.12.2. Convergence

The iteration process described in the previous section continues until convergence is achieved. The maximum number of allowed equilibrium iterations (input on NEQIT command) are performed in order to obtain convergence.

Convergence is assumed when

 (14–150)

and/or

 (14–151)

where {R} is the residual vector:

 (14–152)

which is the right-hand side of the Newton-Raphson Equation 14–147. {Δui} is the DOF increment vector, εR and εu are tolerances (CNVTOL,,,`TOLER`) and Rref and uref are reference values (CNVTOL,,`VALUE`). |||| is a vector norm; that is, a scalar measure of the magnitude of the vector (defined below).

Convergence, therefore, is obtained when size of the residual (disequilibrium) is less than a tolerance times a reference value and/or when the size of the DOF increment is less than a tolerance times a reference value. The default is to use out-of-balance convergence checking only.

There are three available norms (CNVTOL,,,,`NORM`) to choose from:

1. Infinite norm

2. L1 norm

3. L2 norm

For DOF increment convergence, substitute Δu for R in the above equations. The infinite norm is simply the maximum value in the vector (maximum residual or maximum DOF increment), the L1 norm is the sum of the absolute value of the terms, and the L2 norm is the square root of the sum of the squares (SRSS) value of the terms, also called the Euclidean norm. The default is to use the L2 norm.

The default out-of-balance reference value Rref is ||{Fa}||. For DOFs with imposed displacement constraints, the default Rref value is {Fnr} at those DOFs.

In the following cases, the default Rref value is the specified or default minimum reference value set via the CNVTOL,,,,,`MINREF` command:

• For structural DOFs if Rref falls below 1.0E-2 (typically occurring in rigid-body motion analyses, such as those involving stress-free rotation)

• For thermal DOFs if Rref falls below 1.0E-6

• For EMAG DOFs if Rref falls below 1.0E-12

• For all other DOFs if Rref is equal to 0.0

The default reference value uref is ||{u}||.

14.12.3. Predictor

The solution used for the start of each time step n {un,0} is usually equal to the current DOF solution {un -1}. The tangent matrix [Kn,0] and restoring load {Fn,0} are based on this configuration. The predictor option (PRED command) extrapolates the DOF solution using the previous history in order to take a better guess at the next solution.

In static analyses, the prediction is based on the displacement increments accumulated over the previous time step, factored by the time-step size:

 (14–153)

where:

 {Δun} = displacement increment accumulated over the previous time step n = current time step
 (14–154)

and β is defined as:

 (14–155)

where:

 Δtn = current time-step size Δtn-1 = previous time-step size

β is not allowed to be greater than 5.

In transient analyses, the prediction is based on the current velocities and accelerations using the Newmark formulas for structural DOFs:

 (14–156)

where:

 = current displacements, velocities and accelerations Δtn = current time-step size α = Newmark parameter (input on TINTP command)

For thermal, magnetic and other first order systems, the prediction is based on the trapezoidal formula:

 (14–157)

where:

 {un - 1} = current temperatures (or magnetic potentials) = current rates of these quantities θ = trapezoidal time integration parameter (input on TINTP command)

See Transient Analysis for more details on the transient procedures.

The subsequent equilibrium iterations provide DOF increments {Δu} with respect to the predicted DOF value {un,0}, hence this is a predictor-corrector algorithm.

Adaptive descent (Adptky on the NROPT command) is a technique which switches to a “stiffer” matrix if convergence difficulties are encountered, and switches back to the full tangent as the solution convergences, resulting in the desired rapid convergence rate (Eggert([152])).

The matrix used in the Newton-Raphson equation (Equation 14–147) is defined as the sum of two matrices:

 (14–158)

where:

 [KS] = secant (or most stable) matrix [KT] = tangent matrix ξ = descent parameter

The program adaptively adjusts the descent parameter (ξ) during the equilibrium iterations as follows:

1. Start each substep using the tangent matrix (ξ = 0).

2. Monitor the change in the residual ||{R}||2 over the equilibrium iterations:

If it increases (indicating possible divergence):

• remove the current solution if ξ < 1, reset ξ to 1 and redo the iteration using the secant matrix

• if already at ξ = 1, continue iterating

If it decreases (indicating converging solution):

• If ξ = 1 (secant matrix) and the residual has decreased for three iterations in a row (or 2 if ξ was increased to 1 during the equilibrium iteration process by (a.) above), then reduce ξ by a factor of 1/4 (set it to 0.25) and continue iterating.

• If the ξ < 1, decrease it again by a factor of 1/4 and continue iterating. Once ξ is below 0.0156, set it to 0.0 (use the tangent matrix).

3. If a negative pivot message is encountered (indicating an ill-conditioned matrix):

• If ξ < 1, remove the current solution, reset ξ = 1 and redo the iteration using the secant matrix.

• If ξ = 1, bisect the time step if automatic time stepping is active, otherwise terminate the execution.

The nonlinearities which make use of adaptive descent (that is, they form a secant matrix if ξ > 0) include: plasticity, contact, stress stiffness with large strain, nonlinear magnetics using the scalar potential formulation, the concrete element SOLID65 with KEYOPT(7) = 1, and the membrane shell element SHELL41 with KEYOPT(1) = 2. Adaptive descent is used by default in these cases unless the line search or arc-length options are on. It is only available with full Newton-Raphson, where the matrix is updated every iteration. Full Newton-Raphson is also the default for plasticity, contact and large strain nonlinearities.

14.12.5. Line Search

The line search option (accessed with LNSRCH command) attempts to improve a Newton-Raphson solution {Δui} by scaling the solution vector by a scalar value termed the line search parameter.

Consider Equation 14–148 again:

 (14–159)

In some solution situations, the use of the full {Δui} leads to solution instabilities. Hence, if the line search option is used, Equation 14–159 is modified to be:

 (14–160)

where:

 s = line search parameter, 0.05 < s < 1.0

s is automatically determined by minimizing the energy of the system, which reduces to finding the zero of the nonlinear equation:

 (14–161)

where:

 gs = gradient of the potential energy with respect to s

An iterative solution scheme based on regula falsi is used to solve Equation 14–161 (Schweizerhof and Wriggers([153])). Iterations are continued until either:

1. gs is less than 0.5 go, where go is the value of Equation 14–161 at s = 0.0 (that is, using for {Fnr (s{Δu})}).

2. gs is not changing significantly between iterations.

3. Six iterations have been performed.

If go > 0.0, no iterations are performed and s is set to 1.0. s is not allowed below 0.05.

The scaled solution {Δui} is used to update the current DOF values {ui+1} in Equation 14–148 and the next equilibrium iteration is performed.

14.12.6. Arc-Length Method

The arc-length method (accessed with ARCLEN,ON) is suitable for nonlinear static equilibrium solutions of unstable problems.

Application of the arc-length method involves the tracing of a complex path in the load-displacement response into the buckling/post buckling regimes. The arc-length method uses Crisfield’s method as describe in ([174]) to prevent any fluctuation of the step size during equilibrium iterations. It is assumed that all load magnitudes can be controlled by a single scalar parameter (that is, the total load factor).

An unsmooth or discontinuous load-displacement response, which is sometimes seen in contact analyses and elastic-perfectly plastic analyses, cannot be traced effectively by the arc-length solution method.

Mathematically, the arc-length method can be viewed as the trace of a single equilibrium curve in a space spanned by the nodal displacement variables and the total load factor. Therefore, all options of the Newton-Raphson method are still the basic method for the arc-length solution. As the displacement vectors and the scalar load factor are treated as unknowns, the arc-length method itself is an automatic load step method; therefore, AUTOTS,ON is not needed.

For problems with sharp turns in the load-displacement curve or with path-dependent materials, it is necessary to limit the initial arc-length radius (NSUBST) and the arc-length radius augmentation (via the `MAXARC` argument of the ARCLEN command). During the solution, the arc-length method will vary the arc-length radius at each arc-length substep according to the degree of nonlinearity that is involved. The range of variation of the arc-length radius is limited by the maximum and minimum multipliers (`MAXARC` and `MINARC` on the ARCLEN command).

In the arc-length procedure, Equation 14–147 is recast and associated with the total load factor λ:

 (14–162)

Figure 14.13: Arc-Length Approach with Full Newton-Raphson Method is a graphical representation of the arc-length method. Writing the proportional loading factor λ in an incremental form at substep n and iteration i yields:

 (14–163)

Following Equation 14–163, the incremental displacement {Δui} can be written in two parts:

 (14–164)

where:

 (14–165)
 (14–166)

In each arc-length iteration, it is necessary to use Equation 14–165 and Equation 14–166 to solve for and .

We can define one vector between the previous equilibrium point and a point determined in iteration i by:

 (14–167)

where:

 {Δun} = the displacement increment accumulated over the current time step β = the scaling vector for unit displacements

The norm of this vector is:

 (14–168)

At iteration i+1, the vector becomes:

 (14–169)

Crisfield’s method assumes that the norm of the vector is constant along the equilibrium iterations, that is:

 (14–170)

Substituting Equation 14–167 and Equation 14–169 and using Equation 14–164, the following quadratic equation results:

 (14–171)

where:

This system has two real roots, Δλ(1) and Δλ(2), which both satisfy the constant arc-length radius.

For each of these roots, we can define the angle between vector tn-1 in the previously converged substep, and vector ti+1 in the current substep. This angle may be obtained from:

 (14–172)

To move the equilibrium path forward, we choose the route for which the cosine of the associated angle is the closest to 1.

Finally, the solution vectors are updated according to (see Figure 14.13: Arc-Length Approach with Full Newton-Raphson Method):

 (14–173)

and

 (14–174)

Values of λn and Δλ are available in POST26 corresponding to labels ALLF and ALDLF, respectively, on the SOLU command. The normalized arc-length radius label ARCL (also on SOLU) corresponds to value , where is the initial arc-length radius defined through Equation 14–168 and Equation 14–175 (an arc-length radius at the first iteration of the first substep).

Defining K0 as the initial tangent matrix, Fa as the full external load, and λ0 as the initial load factor (specified using the NSUBST command), the initial arc-length radius t0 is determined by:

 (14–175)

where:

The factors `MAXARC` and `MINARC` (input on the ARCLEN command) are used to define the limits of the arc-length radius by the following formulas:

tMAX = `MAXARC` * t0

tMIN = `MINARC` * t0