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

(14–152)

where:

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

(14–154)

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

(14–155)

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 (i.e., 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–153 and Equation 14–155) 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–156)

and/or

(14–157)

where {R} is the residual vector:

(14–158)

which is the right-hand side of the Newton-Raphson Equation 14–153. {Δ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 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:

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 R_ref is ||{F^a}||. For DOFs with imposed displacement constraints, the default R_ref value is {F^nr} at those DOFs.

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

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

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.

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

(14–159)

where:

(14–160)

and β is defined as:

(14–161)

where:

β 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–162)

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

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([152])).

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

(14–164)

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–154 again:

(14–165)

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–165 is modified to be:

(14–166)

where:

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

(14–167)

where:

An iterative solution scheme based on regula falsi is used to solve Equation 14–167 (Schweizerhof and Wriggers([153])). 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–154 and the next equilibrium iteration is performed.

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

(14–168)

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

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

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

(14–170)

where:

(14–171)

(14–172)

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

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

(14–173)

where:

The norm of this vector is:

(14–174)

At iteration i+1, the vector becomes:

(14–175)

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

(14–176)

Substituting Equation 14–173 and Equation 14–175 and using Equation 14–170, the following quadratic equation results:

(14–177)

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

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

and

(14–180)

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–174 and Equation 14–181 (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–181)

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: