8.29. Consider the cantilever beam in Example 8.1 with a concentrated load corresponding to degree of freedom $U_{4}$ . The governing equations are $$ \left[ \begin{array}{r r r r} 7 & - 4 & 1 & 0 \\ - 4 & 6 & - 4 & 1 \\ 1 & - 4 & 5 & - 2 \\ 0 & 1 & - 2 & 1 \end{array} \right] \left[ \begin{array}{l} U _ {1} \\ U _ {2} \\ U _ {3} \\ U _ {4} \end{array} \right] = \left[ \begin{array}{l} 0 \\ 0 \\ 0 \\ 1 \end{array} \right] $$ Solve this set of equations by the conjugate gradient method using two different preconditioners that you shall propose. 8.30. Assume that each eight-node element in Fig. E8.21 is replaced by a 20-node three-dimensional element. Estimate the number of matrix elements under the skyline (to be stored in a direct solution) and the number of actually nonzero matrix elements (to be stored in an iterative solution.) # 8.4 SOLUTION OF NONLINEAR EQUATIONS We discussed in Sections 6.1 and 6.2 that the basic equations to be solved in nonlinear analysis are, at time $t + \Delta t$ , $$ ^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} = \mathbf {0} \tag {8.79} $$ where the vector $^{t+\Delta t}\mathbf{R}$ stores the externally applied nodal loads and $^{t+\Delta t}\mathbf{F}$ is the vector of nodal point forces that are equivalent to the element stresses. Both vectors in (8.79) are evaluated using the principle of virtual displacements. Since the nodal point forces $^{t+\Delta t}\mathbf{F}$ depend nonlinearly on the nodal point displacements, it is necessary to iterate in the solution of (8.79). We introduced in Section 6.1 the Newton–Raphson iteration, in which, assuming that the loads are independent of the deformations, we solve for $i = 1, 2, 3, \ldots$ . $$ \Delta \mathbf {R} ^ {(i - 1)} = ^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i - 1)} \tag {8.80} $$ $$ { } ^ { t + \Delta t } \mathbf { K } ^ { ( i - 1 ) } \Delta \mathbf { U } ^ { ( i ) } = \Delta \mathbf { R } ^ { ( i - 1 ) } \tag {8.81} $$ $$ { } ^ { t + \Delta t } \mathbf { U } ^ { ( i ) } = { } ^ { t + \Delta t } \mathbf { U } ^ { ( i - 1 ) } + \Delta \mathbf { U } ^ { ( i ) } \tag {8.82} $$ with $t + \Delta t\mathbf{U}^{(0)} = {}^t\mathbf{U};\qquad {}^{t + \Delta t}\mathbf{F}^{(0)} = {}^t\mathbf{F}$ (8.83) These equations were obtained by linearizing the response of the finite element system about the conditions at time $t + \Delta t$ , iteration $(i - 1)$ . In each iteration we calculate in (8.80) an out-of-balance load vector that yields an increment in displacements obtained in (8.81), and we continue the iteration until the out-of-balance load vector $\Delta \mathbf{R}^{(i-1)}$ or the displacement increments $\Delta \mathbf{U}^{(i)}$ are sufficiently small. The objective in this section is to discuss the above iterative scheme and others for the solution of (8.81) in more detail. Important ingredients of all solution schemes to be presented are the calculation of the vector $t + \Delta t \mathbf{F}^{(i)}$ and the tangent stiffness matrix $t + \Delta t \mathbf{K}^{(i-1)}$ , and the solution of equations of the form (8.81). The appropriate evaluation of nodal point force vectors and tangent stiffness matrices was discussed in Chapter 6, and the solution of the linearized equations in (8.81) was presented in Sections 8.2 and 8.3; hence, the only, but very important, aspect of concern now is the construction of iterative schemes of the form (8.80) to (8.82) that display good convergence characteristics and can be employed effectively. The methods we now present are basic techniques that in practice would be combined in a self-adaptive procedure that chooses load steps, iterative method, and convergence criteria automatically depending on the problem considered and solution accuracy sought. # 8.4.1 Newton-Raphson Schemes The most frequently used iteration schemes for the solution of nonlinear finite element equations are the Newton-Raphson iteration given in (8.80) to (8.83) and closely related techniques. Because of the importance of the Newton-Raphson method, let us derive the procedure in a more formal manner. The finite element equilibrium requirements amount to finding the solution of the equations $$ \mathbf {f} (\mathbf {U} ^ {*}) = \mathbf {0} \tag {8.84} $$ where $\mathbf{f}(\mathbf{U}^{*}) = {}^{t + \Delta t}\mathbf{R}(\mathbf{U}^{*}) - {}^{t + \Delta t}\mathbf{F}(\mathbf{U}^{*})$ (8.85) We denote here and in the following the complete array of the solution as U\* but realize that this vector may also contain variables other than displacements, for example, pressure variables and rotations (see Sections 6.4 and 6.5). Assume that in the iterative solution we have evaluated $t + \Delta t \mathbf{U}^{(i-1)}$ ; then a Taylor series expansion gives $$ \mathbf {f} (\mathbf {U} ^ {*}) = \mathbf {f} \left(^ {t + \Delta t} \mathbf {U} ^ {(i - 1)}\right) + \left[ \frac {\partial \mathbf {f}}{\partial \mathbf {U}} \right] \Bigg | _ {t + \Delta t \mathbf {U} ^ {(i - 1)}} \left(\mathbf {U} ^ {*} - ^ {t + \Delta t} \mathbf {U} ^ {(i - 1)}\right) + \text { higher - order terms } \tag {8.86} $$ Substituting from (8.85) into (8.86) and using (8.84), we obtain $$ \left[ \frac {\partial \mathbf {F}}{\partial \mathbf {U}} \right] \bigg | _ {t + \Delta t _ {\mathbf {U} (i - 1)}} \left(\mathbf {U} ^ {*} - ^ {t + \Delta t} \mathbf {U} ^ {(i - 1)}\right) + \text { higher - order terms } = ^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i - 1)} \tag {8.87} $$ where we assumed that the externally applied loads are deformation-independent [see (6.83) and (6.84) regarding deformation-dependent loading]. Neglecting the higher-order terms in (8.87), we can calculate an increment in the displacements, $$ { } ^ { t + \Delta t } \mathbf { K } ^ { ( i - 1 ) } \Delta \mathbf { U } ^ { ( i ) } = { } ^ { t + \Delta t } \mathbf { R } - { } ^ { t + \Delta t } \mathbf { F } ^ { ( i - 1 ) } \tag {8.88} $$ where $t+\Delta t\mathbf{K}^{(i-1)}$ is the current tangent stiffness matrix $$ { } ^ { t + \Delta t } \mathbf { K } ^ { ( i - 1 ) } = \left. \left[ \frac { \partial \mathbf { F } } { \partial \mathbf { U } } \right] \right| _ { t + \Delta t _ { \mathbf { U } } ( t - 1 ) } \tag {8.89} $$ and the improved displacement solution is $$ { } ^ { t + \Delta t } \mathbf { U } ^ { ( i ) } = { } ^ { t + \Delta t } \mathbf { U } ^ { ( i - 1 ) } + \Delta \mathbf { U } ^ { ( i ) } \tag {8.90} $$ The relations in (8.88) and (8.90) constitute the Newton-Raphson solution of (8.79). Since an incremental analysis is performed with time (or load) steps of size $\Delta t$ (see Chapter 6), the initial conditions in this iteration are ${}^{t+\Delta t}\mathbf{K}^{(0)} = {}^{t}\mathbf{K}$ , ${}^{t+\Delta t}\mathbf{F}^{(0)} = {}^{t}\mathbf{F}$ and ${}^{t+\Delta t}\mathbf{U}^{(0)} = {}^{t}\mathbf{U}$ . The iteration is continued until appropriate convergence criteria are satisfied as discussed in Section 8.4.4. A characteristic of this iteration is that a new tangent stiffness matrix is calculated in each iteration, which is why this method is also referred to as the full Newton-Raphson method. We shall mention below methods in which a current tangent stiffness matrix is not used, and these techniques are therefore not full Newton-Raphson methods (but related techniques). Figure 8.11 illustrates the process of solution when used for a single degree of freedom system. The nonlinear response characteristics are such that convergence is rapidly obtained. However, we can imagine a more complex response characteristic with a starting 

line

| Displacement | Load (t+ΔtR) | Load (t+ΔtK(0)) | | ------------ | ------------ | --------------- | | t_u | t_R | Δu(1) | | t+Δt_R | t+Δt_R | t+Δt_F(0) | | t+Δt_F(0) | t+Δt_F(0) | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_F(1) | t+Δt_F(1) | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_F(1) | t+Δt_F(1) | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_F(1) | t+Δt_R | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_F(1) | t+Δt_R | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+ Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_R | t+Δt_R | t+Δt_F(1) | | t+Δt_R (t+Δt_R) | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R) | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R) | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (t+Δt_R) | t+ Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t + Δt_R (t+Δt_R)| t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (0) (t+Δt_R) | t + Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δ t_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R (t+Δt_R)| t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (0) (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (0) (t+Δt_R) | t+Δt_F (1) (t+Δt_R) | | t+Δt_R | t+Δt_R (0) (t+Δt_R) | t+Δt_F (1) (t+Δt_R) | | t+Δt_R | t+Δt_R (0) (t+Δt_R) | t+Δt_F (1) (t+Δt_R) | | t+Δt_R | t+Δt _t+Δt_R | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt _t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δ t_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1)(t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δr_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R)| t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R)| t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R)| t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δr_R)| t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R)| t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R)| t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t + Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δr (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+ Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | t+Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) | | Δt_R | t+Δt_R (t+Δt_R) | t+Δt_F(1) (t+Δt_R) |


line

| Point Label | Description | |-------------|-------------| | t+ΔtR | t+ΔtR (t+ΔtR) | | t+ΔtF(0) | t+ΔtF(0) | | t+ΔtF(1) | t+ΔtF(1) | | Δu(1) | Δu(1) | | Δu(2) | Δu(2) | | t+ΔtU | t+ΔtU |

Figure 8.11 Illustration of Newton-Raphson iteration in solution of a (generic) single degree of freedom system. Top shows load-displacement relation, bottom shows iteration for zero of function f used in (8.84). Here $f = r^{+ \Delta t} R - r^{+ \Delta t} F(u)$ . point of iteration for which the procedure does not converge (see Exercise 8.31). Hence, the representation in Fig. 8.11 is rather simplistic because a very special case is considered—that of a well-behaved single degree of freedom system. In the solution of systems with many degrees of freedom, the response curves will in general be rather nonsmooth and complicated. The Newton-Raphson iteration is demonstrated in the solution of a simple problem as follows. EXAMPLE 8.23: For a single degree of freedom system we have $$ ^ {t + \Delta t} R = 1 0; \quad^ {t + \Delta t} F = 4 + 2 \left| (^ {t + \Delta t} U) ^ {1 / 2} \right| $$ and $^{\prime}U = 1$ . Use the Newton-Raphson iteration to calculate $^{t + \Delta t}U$ . In this case we have, corresponding to (8.88) as the governing equation, $$ \left(\frac {1}{\left| \left(t ^ {+ \Delta t} U ^ {(i - 1)}\right) ^ {1 / 2} \right|}\right) \Delta U ^ {(i)} = 6 - 2 \left| \sqrt {t ^ {+ \Delta t} U ^ {(i - 1)}} \right| \tag {a} $$ Using (a) with $r + \Delta r U^{(0)} = 1$ , we obtain $$ { } ^ { t + \Delta t } U ^ { ( 1 ) } = 5 . 0 0 0 0 ; \quad { } ^ { t + \Delta t } U ^ { ( 2 ) } = 8 . 4 1 6 4 $$ $$ ^ {t + \Delta t} U ^ {(3)} = 8. 9 9 0 2; \quad^ {t + \Delta t} U ^ {(4)} = 9. 0 0 0 0 $$ and convergence is achieved in four iterations. Since the Newton-Raphson iteration is so widely used in finite element analysis and indeed represents the primary solution scheme for nonlinear finite element equations, it is appropriate that we summarize some major properties of the method (see, for example, D. P. Bertsekas [A]). # The first property is If the tangent stiffness matrix $^{t + \Delta t}\mathbf{K}^{(i - 1)}$ is nonsingular, if $\mathbf{f}$ and its first derivatives with respect to the solution variables (i.e. the elements of the tangent stiffness matrix) are continuous in a neighborhood of $\mathbf{U}^*$ , and if $^{t + \Delta t}\mathbf{U}^{(i - 1)}$ lies in that neighborhood, then $^{t + \Delta t}\mathbf{U}^{(i)}$ will be closer to $\mathbf{U}^*$ than $^{t + \Delta t}\mathbf{U}^{(i - 1)}$ and the sequence of iterative solutions generated by the algorithm (8.88) to (8.90) converges to $\mathbf{U}^*$ . # The second property is If the tangent stiffness matrix also satisfies $$ \left\| ^ {t + \Delta t} \mathbf {K} \left| \mathbf {u} _ {1} - ^ {t + \Delta t} \mathbf {K} \right| \mathbf {u} _ {2} \right\| \leq L \left\| \mathbf {U} _ {1} - \mathbf {U} _ {2} \right\| \tag {8.91} $$ for all $U_{1}$ and $U_{2}$ in the neighborhood of $U^{*}$ and L > 0, then convergence is quadratic. This means that if the error after iteration i is of order $\epsilon$ , then the error after iteration $i + 1$ will be of the order $\epsilon^{2}$ . The condition in (8.91) is referred to as Lipschitz continuity; it is stronger than mere continuity in the stiffness matrix but weaker than differentiability of the matrix. The practical consequence of these properties is that if the current solution iterate is sufficiently close to the solution $U^{*}$ and if the tangent stiffness matrix does not change abruptly, we can expect rapid (i.e., quadratic) convergence. The assumption is of course that the exact tangent stiffness matrix is used in the iteration; that is, (8.89) must be satisfied, meaning that $t^{+\Delta t}\mathbf{K}^{(i-1)}$ must be evaluated consistent with the evaluation of $t^{+\Delta t}\mathbf{F}^{(i-1)}$ (see Chapter 6 and in particular Sections 6.3.1 and 6.6.3). On the other hand, if the current solution iterate is not sufficiently close to $U^{*}$ and/or the stiffness matrix used is not the exact tangent matrix and/or changes abruptly, then the iteration may diverge. In an effective finite element program, the exact tangent stiffness matrix will be used, if possible, and hence the primary procedure for reaching convergence (if convergence difficulties are encountered) is to decrease the magnitude of the load step. Considering the Newton-Raphson iteration it is recognized that in general the major computational cost per iteration lies in the calculation and factorization of the tangent stiffness matrix. Since these calculations can be quite expensive when large-order systems are considered, the use of a modification of the full Newton-Raphson algorithm can be effective. 

line

| Displacement | Load (t+Δt_R) | Load (t+Δt_R - t+Δt_F^(0)) | Load (t+Δt_R - t+Δt_F^(1)) | Δu^(1) | Δu^(2) | | ------------ | ------------- | --------------------------- | --------------------------- | ------ | ------ | | t_u | t_R | t+Δt_R | t+Δt_R | t+Δt_R | t+Δt_R | | t+Δt_u | t+Δt_R | t+Δt_R | t+Δt_R | t+Δt_R | t+Δt_R |


line

| Displacement | Load (t+ΔtR) | Load (t+ΔtF(0)) | Load (t+ΔtF(1)) | Slope (t+ΔtK(0)) = tK | | ------------ | ------------ | --------------- | --------------- | ---------------------- | | t_u | t_R | t_R | t_R | t_R | | t+Δt_u | t+Δt_R | t+Δt_R | t+Δt_R | t+Δt_R | | t+Δt_F(0) | t+Δt_F(0) | t+Δt_F(0) | t+Δt_F(0) | t+Δt_F(0) | | t+Δt_F(1) | t+Δt_F(1) | t+Δt_F(1) | t+Δt_F(1) | t+Δt_F(1) |

Figure 8.12 Illustration of initial stress and modified Newton-Raphson methods. The problem in Fig. 8.11 is considered. One such modification is to use the initial stiffness matrix ${}^{0}K$ in (8.88) and thus operate on the equations: $$ { } ^ { 0 } \mathbf { K } \Delta \mathbf { U } ^ { ( i ) } = { } ^ { t + \Delta t } \mathbf { R } - { } ^ { t + \Delta t } \mathbf { F } ^ { ( i - 1 ) } \tag {8.92} $$ with the initial conditions ${}^{t+\Delta t}\mathbf{F}^{(0)} = {}^{t}\mathbf{F}, {}^{t+\Delta t}\mathbf{U}^{(0)} = {}^{t}\mathbf{U}$ . In this case only the matrix ${}^{0}K$ needs to be factorized, thus avoiding the expense of recalculating and factorizing many times the coefficient matrix in (8.88). This “initial stress” method corresponds to a linearization of the response about the initial configuration of the finite element system and may converge very slowly or even diverge. In the modified Newton-Raphson iteration an approach somewhat in between the full Newton-Raphson iteration and the initial stress method is employed. In this method we use $$ { } ^ { \tau } \mathbf { K } \Delta \mathbf { U } ^ { ( i ) } = { } ^ { t + \Delta t } \mathbf { R } - { } ^ { t + \Delta t } \mathbf { F } ^ { ( i - 1 ) } \tag {8.93} $$ with the initial conditions ${}^{t+\Delta t}\mathbf{F}^{(0)} = {}^{t}\mathbf{F}, {}^{t+\Delta t}\mathbf{U}^{(0)} = {}^{t}\mathbf{U}$ , and $\tau$ corresponds to one of the accepted equilibrium configurations at times 0, $\Delta t$ , $2 \Delta t$ , $\ldots$ , or t (see Example 6.4). The modified Newton-Raphson iteration involves fewer stiffness reformations than the full Newton-Raphson iteration and bases the stiffness matrix update on an accepted equilibrium configuration. The choice of time steps when the stiffness matrix should be updated depends on the degree of nonlinearity in the system response; i.e. the more nonlinear the response, the more often the updating should be performed. Figure 8.12 illustrates the performance of the initial stress and the modified Newton-Raphson methods for the single degree of freedom system already considered in Fig. 8.11. With the very large range of system properties and nonlinearities that may be encountered in engineering analysis, we find that the effectiveness of the above solution approaches depends on the specific problem considered. The most powerful procedure for reaching convergence is the full Newton-Raphson iteration in (8.88) to (8.90), but if the initial stress or modified Newton-Raphson method can be employed, the solution cost may be reduced significantly. Hence, in practice, these solution options can also be very valuable, and an automatic procedure that self-adaptively chooses an effective technique is most attractive. # 8.4.2 The BFGS Method As an alternative to forms of Newton-Raphson iteration, a class of methods known as matrix update methods or quasi-Newton methods has been developed for iteration on nonlinear systems of equations (see J. E. Dennis, Jr. [A]). These methods involve updating the coefficient matrix (or rather its inverse) to provide a secant approximation to the matrix from iteration $(i - 1)$ to $(i)$ . That is, defining a displacement increment $$ \boldsymbol {\delta} ^ {(i)} = ^ {t + \Delta t} \mathbf {U} ^ {(i)} - ^ {t + \Delta t} \mathbf {U} ^ {(i - 1)} \tag {8.94} $$ and an increment in the out-of-balance loads, using (8.80), $$ \boldsymbol {\gamma} ^ {(i)} = \Delta \mathbf {R} ^ {(i - 1)} - \Delta \mathbf {R} ^ {(i)} \tag {8.95} $$ the updated matrix $^{t+\Delta t}\mathbf{K}^{(i)}$ should satisfy the quasi-Newton equation $$ { } ^ { t + \Delta t } \mathbf { K } ^ { ( i ) } \delta ^ { ( i ) } = \gamma ^ { ( i ) } \tag {8.96} $$ These quasi-Newton methods provide a compromise between the full re-formation of the stiffness matrix performed in the full Newton-Raphson method and the use of a stiffness matrix from a previous configuration as is done in the modified Newton-Raphson method. Among the quasi-Newton methods available, the BFGS (Broyden-Fletcher-Goldfarb-Shanno) method appears to be most effective. In the BFGS method, the following procedure is employed in iteration (i) to evaluate $t^{+\Delta t}U^{(i)}$ and $t^{+\Delta t}K^{(i)}$ , where $t^{+\Delta t}K^{(0)} = {}^{T}K$ (see H. Matthies and G. Strang [A] and K. J. Bathe and A. P. Cimento [A]). Step 1: Evaluate a displacement vector increment: $$ \Delta \overline {{{\mathbf {U}}}} = (t ^ {+ \Delta t} \mathbf {K} ^ {- 1}) ^ {(i - 1)} (t ^ {+ \Delta t} \mathbf {R} - t ^ {+ \Delta t} \mathbf {F} ^ {(i - 1)}) \tag {8.97} $$ This displacement vector defines a “direction” for the actual displacement increment. Step 2: Perform a line search in the direction $\Delta\overline{U}$ to satisfy “equilibrium” in this direction. In this line search we evaluate the displacement vector $$ { } ^ { t + \Delta t } \mathbf { U } ^ { ( i ) } = { } ^ { t + \Delta t } \mathbf { U } ^ { ( i - 1 ) } + \beta \Delta \overline { { \mathbf { U } } } \tag {8.98} $$ where $\beta$ is a scalar multiplier, and we calculate the out-of-balance loads corresponding to these displacements $(t^{+\Delta t}\mathbf{R}-t^{+\Delta t}\mathbf{F}^{(i)})$ . The parameter $\beta$ is varied until the component of the out-of-balance loads in the direction $\Delta\overline{U}$ , as defined by the inner product $\Delta\overline{\mathbf{U}}^{T}(t^{+\Delta t}\mathbf{R}-t^{+\Delta t}\mathbf{F}^{(i)})$ , is small. This condition is satisfied when, for a convergence tolerance STOL, the following equation is satisfied: $$ \Delta \overline {{{\mathbf {U}}}} ^ {T} \left(^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i)}\right) \leq \text { STOL } \Delta \overline {{{\mathbf {U}}}} ^ {T} \left(^ {t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i - 1)}\right) \tag {8.99} $$ The final value of $\beta$ for which (8.99) is satisfied determines ${}^{t+\Delta t}\mathbf{U}^{(i)}$ in (8.98). We can now calculate $\mathbf{\delta}^{(i)}$ and $\mathbf{\gamma}^{(i)}$ using (8.94) and (8.95) and proceed with the evaluation of the matrix update that satisfies (8.96). Step 3: Evaluate the correction to the coefficient matrix. In the BFGS method the updated matrix can be expressed in product form: $$ (^ {i + \Delta t} \mathbf {K} ^ {- 1}) ^ {(i)} = \mathbf {A} ^ {(i) ^ {T}} (^ {i + \Delta t} \mathbf {K} ^ {- 1}) ^ {(i - 1)} \mathbf {A} ^ {(i)} \tag {8.100} $$ where the matrix $\mathbf{A}^{(i)}$ is an $n \times n$ matrix of the simple form $$ \mathbf {A} ^ {(i)} = \mathbf {I} + \mathbf {v} ^ {(i)} \mathbf {w} ^ {(i) T} \tag {8.101} $$ The vectors $\mathbf{v}^{(i)}$ and $\mathbf{w}^{(i)}$ are calculated from the known nodal point forces and displacements using $$ \mathbf {v} ^ {(i)} = - \left(\frac {\boldsymbol {\delta} ^ {(i) ^ {T}} \boldsymbol {\gamma} ^ {(i)}}{\boldsymbol {\delta} ^ {(i) ^ {T}} \iota + \Delta t \mathbf {K} ^ {(i - 1)} \boldsymbol {\delta} ^ {(i)}}\right) ^ {1 / 2} \iota + \Delta t \mathbf {K} ^ {(i - 1)} \boldsymbol {\delta} ^ {(i)} - \boldsymbol {\gamma} ^ {(i)} \tag {8.102} $$ and $$ \mathbf {w} ^ {(i)} = \frac {\delta^ {(i)}}{\delta^ {(i) T} \gamma^ {(i)}} \tag {8.103} $$ The vector $^{t + \Delta t}\mathbf{K}^{(i - 1)}\delta^{(i)}$ in (8.102) is equal to $\beta(^{t + \Delta t}\mathbf{R} - ^{t + \Delta t}\mathbf{F}^{(i - 1)})$ and has already been computed. Since the product defined in (8.100) is positive definite and symmetric, to avoid numerically dangerous updates, the condition number $c^{(i)}$ of the updating matrix $\mathbf{A}^{(i)}$ is calculated: $$ c ^ {(i)} = \left(\frac {\boldsymbol {\delta} ^ {(i) ^ {T}} \boldsymbol {\gamma} ^ {(i)}}{\boldsymbol {\delta} ^ {(i) ^ {T}} {} _ {t + \Delta t} \mathbf {K} ^ {(i - 1)} \boldsymbol {\delta} ^ {(i)}}\right) ^ {1 / 2} \tag {8.104} $$ This condition number is then compared with some preset tolerance, a large number, and the updating is not performed if the condition number exceeds this tolerance. Considering the actual computations involved, it should be recognized that using the matrix updates defined above, the calculation of the search direction in (8.97) can be rewritten as $$ \begin{array}{l} \Delta \overline {{{\mathbf {U}}}} = (\mathbf {I} + \mathbf {w} ^ {(i - 1)} \mathbf {v} ^ {(i - 1) T}) \dots (\mathbf {I} + \mathbf {w} ^ {(1)} \mathbf {v} ^ {(1) T}) ^ {\tau} \mathbf {K} ^ {- 1} (\mathbf {I} + \mathbf {v} ^ {(1)} \mathbf {w} ^ {(1) T}) \dots \\ (\mathbf {I} + \mathbf {v} ^ {(i - 1)} \mathbf {w} ^ {(i - 1) T}) ^ {(t + \Delta t} \mathbf {R} - ^ {t + \Delta t} \mathbf {F} ^ {(i - 1)}) \tag {8.105} \\ \end{array} $$ Hence, the search direction can be computed without explicitly calculating the updated matrices or performing any additional costly matrix factorizations as required in the full Newton-Raphson method. As pointed out above, the line search is an integral part of the solution method. Of course, such line searches as performed in (8.98) and (8.99) can also be used in the Newton-Raphson methods presented in Section 8.4.1. With the line search performed within an iteration (i), the expense of the iteration increases, but fewer iterations may be needed for convergence. Also, the line search may prevent divergence of the iterations, and in practice, this increased robustness is the major reason why a line search can in general be effective. We demonstrate the BFGS iteration in the following simple example. EXAMPLE 8.24: Use the BFGS iteration method to solve for $t+\Delta tU$ of the system considered in Example 8.23. Omit the line searches in the solution. Since this is a single degree of freedom system, the solution for $^{t + \Delta t}U$ could be evaluated using only line searches, i.e., by applying (8.99), provided STOL is a tight enough convergence tolerance. However, to demonstrate in this example the basic steps of the BFGS method more clearly [the use of relations (8.94) to (8.96)], we do not include line searches in the iterative solution. In this analysis (8.97) reduces to $$ \Delta \overline {{{U}}} = (^ {t + \Delta t} K ^ {- 1}) ^ {(i - 1)} (6 - 2 | \sqrt {^ {t + \Delta t} U ^ {(i - 1)}} |) $$ with $(^{t + \Delta t}K^{-1})^{(0)} = 1$ , $^{t + \Delta t}U^{(0)} = 1$ and using $\beta = 1.0$ , we obtain the following values:

i	$^{t+\Delta t}U^{(i-1)}$	$\Delta\overline{U}=\delta^{(i)}$	$^{t+\Delta t}U^{(i)}$	$\gamma^{(i)}$	$(^{t+\Delta t}K^{-1})^{(i)}$
1	1.000	4.000	5.000	2.472	1.618
2	5.000	2.472	7.472	0.995	2.485
3	7.472	1.324	8.796	0.465	2.850
4	8.796	0.194	8.991	0.065	2.982
5	8.991	0.009	9.000	0.003	2.999

and convergence is achieved after five iterations. # 8.4.3 Load-Displacement-Constraint Methods An important requirement of a nonlinear analysis is frequently the calculation of the collapse load of a structure. Figure 8.13 illustrates schematically the response of a structural model that we might seek. For very small loads the load-displacement response is linear. 

text_image

Load Postcollapse response A Load decreases Load increments should be smaller Load increments can be large Displacement

Figure 8.13 Collapse response of a structural model (the graph depicts schematically the load that is carried by a structure) Then, as the load increases, the structural response becomes increasingly nonlinear and at point A the collapse load is reached. The response beyond point A is referred to as postcollapse or postbuckling response. We note that in Fig. 8.13 the load first decreases in this regime and then increases again as the displacement increases. Of course, the response depicted in Fig. 8.13 is a simplistic and generic representation because in the analysis of a multiple degree of freedom system a multidimensional “response surface” must be imagined, but Fig. 8.13 shows the essence of our requirements. In order to calculate the response in Fig. 8.13 initially relatively large load increments can be employed, but as the collapse of the structural model is approached, the load increments must become smaller and there is also the difficulty of traversing the collapse point. At that point, the stiffness matrix is singular (the slope of the load-displacement response curve is zero), and beyond that point a special solution procedure that allows for a decrease in load and an increase in displacement must be used to calculate the ensuing response. To solve for the response shown schematically in Fig. 8.13, a load-displacement-constraint method can be used, as in essence proposed by E. Riks [A]. The basic idea of such methods is to introduce a load multiplier that increases or decreases the intensity of the applied loads, so as to obtain fast convergence in each load step, to be able to traverse the collapse point and evaluate the postcollapse response. Various efficient schemes have been proposed in which some numerical details can be important (see M. A. Crisfield [A], E. Ramm [A], and K. J. Bathe and E. N. Dvorkin [C]). However, we shall present in the following only the general approach of these methods and will omit some details that can be found in the references. A basic assumption in the analysis is that the load vector varies proportionally during the response calculation. The governing finite element equations at time $t + \Delta t$ are $$ ^ {t + \Delta t} \lambda \mathbf {R} - ^ {t + \Delta t} \mathbf {F} = \mathbf {0} \tag {8.106} $$ where $t^{+\Delta t}\lambda$ is a (scalar) load multiplier, unknown and to be determined, and R is the reference load vector for the n degrees of freedom of the finite element model. This vector can contain any loading on the structure but is constant throughout the response calculation. The vector $t^{+\Delta t}F$ is our usual vector of n nodal point forces corresponding to the element stresses at time $t+\Delta t$ [see (8.79)]. The value of the load multiplier can increase or decrease, and the increment per step should in general also change, depending on the structural response characteristics. Since (8.106) represents n equations in $n + 1$ unknowns, we need an additional equation that is used to determine the load multiplier. If we apply one of the previously presented methods to solve (8.106), we obtain $$ { } ^ { \tau } \mathbf { K } \Delta \mathbf { U } ^ { ( i ) } = ( { } ^ { t + \Delta t } \lambda ^ { ( i - 1 ) } + \Delta \lambda ^ { ( i ) } ) \mathbf { R } - { } ^ { t + \Delta t } \mathbf { F } ^ { ( i - 1 ) } \tag {8.107} $$ where the coefficient matrix $^{7}K$ corresponds to the solution schemes discussed in the previous sections. The unknowns in the n equations (8.107) are the vector of displacement increments $^{2}$ $\Delta U^{(i)}$ and the load multiplier increment $\Delta\lambda^{(i)}$ . The additional equation required for solution is a constraint equation between $\Delta\lambda^{(i)}$ and $\Delta U^{(i)}$ , of the form, $$ f (\Delta \lambda^ {(i)}, \Delta \mathbf {U} ^ {(i)}) = 0 \tag {8.108} $$ Let us define within a load step $$ \mathbf {U} ^ {(i)} = ^ {t + \Delta t} \mathbf {U} ^ {(i)} - ^ {t} \mathbf {U} \tag {8.109} $$ and $\lambda^{(i)} = {}^{t + \Delta t}\lambda^{(i)} - {}^{t}\lambda$ (8.110) Hence, $\mathbf{U}^{(i)}$ represents the total increment in displacements within the load step [up to iteration (i)] and $\lambda^{(i)}$ represents the corresponding total increment in load multiplier. An effective constraint equation is given by the spherical constant arc length criterion (see M. A. Crisfield [A] and E. Ramm [A]), $$ (\lambda^ {(i)}) ^ {2} + \frac {\mathbf {U} ^ {(i) ^ {T}} \mathbf {U} ^ {(i)}}{\beta} = (\Delta l) ^ {2} \tag {8.111} $$ where $\Delta l$ is the arc length for the step and $\beta$ is a normalizing factor (to render the terms dimensionless). Figure 8.14(a) illustrates this criterion. In practice, the magnitude of $\Delta l$ is selected based on the history of iterations in the previous steps and is reduced in the current step if convergence difficulties are encountered. Typically, $\Delta l$ should be large when the response is almost linear and should be small when the response is highly nonlinear. 

line

| Displacement | Load | | ------------ | ---- | | t+ΔtλR | tλR | | ΔI | ΔI |

(a) Spherical constant arc length criterion 

line

| Displacement | Load | |-------------|------| | t+Δtλ(1)R | t+Δtλ(1)R | | tλR | tλR | | tU | tU | | t+ΔtU(1) | t+ΔtU(1) |

(b) Constant increment of external work criterion Figure 8.14 Load-displacement-constraint criteria (single degree of freedom simplification) Another effective constraint equation (see K. J. Bathe and E. N. Dvorkin [C]), is the scheme of constant increment of external work given by