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

(14–164)

where:

If the coefficient matrix [K] is itself a function of the unknown DOF values (or their derivatives) then Equation 14–164 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–165)

(14–166)

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 {u_i}. The right-hand side of Equation 14–165 is the residual or out-of-balance load vector; that is, 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, {u_i} is the displacement vector and is the restoring force vector calculated from the element stresses. In a thermal analysis, is the conductivity matrix, {u_i} is the temperature vector and is the resisting load vector calculated from the element heat flows. In an electromagnetic analysis, is the Dirichlet matrix, {u_i} 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:

Figure 14.9: Newton-Raphson Solution - One Iteration

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 {F^a}. 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 {F^a} (or at least to within some tolerance). None of the intermediate solutions would be in equilibrium.

Figure 14.10: Newton-Raphson Solution - Next Iteration

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; that is, the final load vector {F^a} is reached by applying the load in increments and performing the Newton-Raphson iterations at each step:

(14–167)

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 {u_i} is "near" the exact solution. Therefore, even without a path-dependent nonlinearity, the incremental approach (that is, applying the loads in increments) is sometimes required in order to obtain a solution corresponding to the final load level.

Figure 14.11: Incremental Newton-Raphson Procedure

When the stiffness matrix is updated every iteration (as indicated in Equation 14–165 and Equation 14–167) 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.

Figure 14.12: Initial-Stiffness Newton-Raphson

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–168)

and/or

(14–169)

where {R} is the residual vector:

(14–170)

which is the right-hand side of the Newton-Raphson Equation 14–165. {Δu_i} is the DOF increment vector, ε_R and ε_u are tolerances (CNVTOL,,,TOLER) and R_ref and u_ref 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 the 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.

An alternative convergence criterion is based on the incremental energy calculation:

(14–171)

where is the incremental system energy during a Newton interation :

(14–172)

It can be interpreted as the work done by an out-of-balance force on the incremental displacement. is the tolerance specified by CNVTOL,,,TOLER, and is the reference value specified by CNVTOL,,VALUE.

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

In the above equation, R can be replaced with different convergence labels, such as Δu for DOF increment convergence. 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 R_ref is ||{F^a}||. For DOFs with imposed displacement constraints, the default R_ref value is {F^nr} at those DOFs.

If the reference value R_ref falls below a minimum reference value specified by CNVTOL,,,,,MINREF, then R_ref is replaced by the minimum reference value. The program also internally calculates a reference value and uses the maximum of this internally calculated value and the MINREF default. (See CNVTOL for a list of default minimum reference values for the various convergence criteria.) The internally calculated minimum reference value depends upon the first iteration residual value and contact forces. Furthermore, this minimum reference also depends on the displacement convergence value and R_ref.

Minimum reference is typically useful for problems with almost no external loading. If the program-calculated minimum reference is still too low, you should manually specify a minimum reference that is appropriate for the problem under consideration. Examples where it may be useful to manually specify the minimum reference include, but are not limited to, problems involving contact initial penetration resolution, free thermal expansions, and rigid-body motion analyses with stress-free rotation. When the program-calculated internal reference is used, the program may sometimes tighten the displacement convergence criterion to improve solution accuracy.

Sometimes, the program-calculated minimum reference might be too high. This can lead to premature convergence, and accuracy of the solution may be affected. This would typically be the case for problems where very low levels of strain are expected, and bodies are not deforming much. For such cases, it is recommended you manually define a tight minimum reference (for example, 0.01 or 0.001) to achieve accurate results.

The default reference value u_ref is calculated using the historic infinity norm of incremental displacement.

If you do not specify the reference value for the energy norm, , the program calculates it as:

(14–173)

where is the incremental energy from the 1st Newton-Raphson iteration in substep , and || || norm represents the absolute value of .

If the residual ({F^a} - {F^nr}) in Equation 14–173 is too small, the program uses the applied force ({F^a}) or the internal force ({F^nr}) to calculate the reference.

In theory, the energy convergence check is equivalent to the combination of the out-of-balance force and the displacement check together. However, by default the energy convergence check is more relaxed than this combination. In practice, if only the energy convergence check is invoked with the default tolerance (the out-of-balance force and displacement convergence are off), the solution could converge in a smaller number of cumulative iterations but may not be as accurate as a solution that includes out-of-balance force and displacement convergence. The energy convergence check provides an option to achieve faster runtime along with a meaningful solution.

The energy convergence check is more suited to single-field structural applications (UX, UY, UZ, ROTX, ROTY, ROTZ) or multiphysics applications with structural degrees of freedom because the energy calculated in the check is the real physical energy of the system. This may not be the case for other multiphysics applications.

In static analyses, the solution used for the start of each time step n {u_n,0} is usually equal to the current DOF solution {u_n-1}. The tangent matrix [K_n,0] and restoring load {F^n,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. The prediction is based on the displacement increments accumulated over the previous time step, factored by the time-step size:

(14–174)

where:

(14–175)

and β is defined as:

(14–176)

where:

= current time-step size

= previous time-step size

β is not allowed to be greater than 5.

Equation 14–174 is called the formulation of the linear predictor. It uses two known solution values, and, to extrapolate the solution of the next substep in a static analysis.

The linear predictor normally calculates a good prediction for a solution that is smooth, such as a translatonal displacement (UX, UY, UZ) in a mildly nonlinear problem.

The solution can also be predicted by quadratic extrapolation as described by the formulation shown in Equation 14–177:

(14–177)

where:

= displacement increment accumulated over the time step prior to the previous time step.

c1, c2 are scalable parameters determined by the current, previous, and prior to previous time steps.

The quadradic predictor in Equation 14–177 is normally a better predictor than the linear predictor for solutions with rotational degrees of freedom (ROTX, ROTY, ROTZ).

In transient analyses, the prediction is based on the current velocities and accelerations and the time integration scheme used (Newmark, HHT, Backward Eulter). For example, using the Newmark formulas for structural DOFs:

(14–178)

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–179)

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 {u_n,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 [153]).

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

(14–180)

where:

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

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, and nonlinear magnetics using the scalar potential formulation. 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.

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

Consider Equation 14–166 again:

(14–181)

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

(14–182)

where:

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

(14–183)

where:

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

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

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

The line search algorithm is very susceptible to round-off errors (Crisfield [441]). Truncation of g_s can help solution repeatability (see LStrun on LNSRCH).

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 (Crisfield [176]) 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–165 is recast and associated with the total load factor λ:

(14–184)

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–185)

Figure 14.13: Arc-Length Approach with Full Newton-Raphson Method

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

(14–186)

where:

(14–187)

(14–188)

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

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

(14–189)

where:

The norm of this vector is:

(14–190)

At iteration i+1, the vector becomes:

(14–191)

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

(14–192)

Substituting Equation 14–189 and Equation 14–191 and using Equation 14–186, the following quadratic equation results:

(14–193)

where:

This system has two real roots, Δλ⁽¹⁾ and Δλ⁽²⁾, which both satisfy the constant arc-length radius.

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

(14–194)

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–195)

and

(14–196)

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–190 and Equation 14–197 (an arc-length radius at the first iteration of the first substep).

Defining K₀ as the initial tangent matrix, F^a as the full external load, and λ₀ as the initial load factor (specified using the NSUBST command), the initial arc-length radius t₀ is determined by:

(14–197)

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: