Fig. 3.1 Simple nonlinear problem
The formulation of a nonlinear problem in FEM leads to a system of nonlinear equations whose unknown values represent the numerical values of the sought functions at the nodes of the computational model. The total number of equations is equal to the total number of degrees of freedom in the model. To solve the system of nonlinear equations FEM programs use algorithms in the form of programmed solvers based on basic mathematical methods for solving such problems. Most often they use an iterative approach based on the Newton-Raphson method which we now remind with a special focus on the equations of structural problems.
Consider a bar of elastic material with Young's modulus $E_0$ and a cross-sectional area $S_0$ fixed and loaded as shown in Fig. 3.1. The bar has only one degree of freedom, that is, vertical displacement $u$, whose dependence on the load force $F$ can be approximately expressed by a simple nonlinear equation [1]
With numerical values $E_0 = 200\,000$ MPa, $S_0 = 800$ mm², $l_0 = 2000$ mm, $h = 50$ mm, equation (3.1) changes to
$$ 0.01\left(u^3 + 150u^2 + 5000u\right) = F \tag{3.2} $$If we want to find the displacement $u$ for a certain force value $F$ iteratively, it is advantageous to write equation (3.2) in a form that makes it easy to control the iterative process
$$ R(u) = 0.01\left(u^3 + 150u^2 + 5000u\right) - F = 0 \tag{3.3} $$where $R(u)$ is called the non-equilibrium force, because the equilibrium equation of the bar (3.2) is fulfilled only when this force equals zero. If, for example, for the calculation $u$ we choose force $F = 4000$ N, the non-equilibrium force will be
$$ R(u) = 0.01\left(u^3 + 150u^2 + 5000u\right) - 4000 \tag{3.4} $$A graphical representation of this function in Fig. 3.2 shows that the displacement of the bar at $F = 4000$ N is approximately $u^* = 35$ mm.
The principle of the Newton-Raphson method can be easily explained by an equation with one unknown. We choose the initial value $u_0$ of the displacement (Fig. 3.2) and express the non-equilibrium force $R(u_0)$. At this point a tangent to the curve is made (this can be called linearization of the function around this point), which marks the displacement $u_1$ on the $u$-axis as the first approximation to the solution. This displacement is easy to determine because, as you can see from the picture, it applies
$$ \tan(\varphi_0) = (dR/du)_{u_0} = \frac{R(u_0)}{u_0 - u_1} \tag{3.5} $$and from that
$$ u_1 = u_0 - \frac{R(u_0)}{(dR/du)_{u_0}} = u_0 - \frac{R(u_0)}{(K_T)_{u_0}} \tag{3.6} $$An analogous repetition of this procedure with values $u_1$, $u_2$ and so on, is the basis for the necessary iterative algorithm. The derivative
$$ K_T(u) = \frac{dR(u)}{du} \tag{3.7} $$is the so-called tangential stiffness which, as can be seen, is not constant in this process.
In the Newton-Raphson method, the root (solution) of a one-dimensional equation is searched as follows:
1. The gradient of the function $R(u)$, the so-called tangential stiffness, is calculated
$$ K_T(u) = \frac{dR(u)}{du} \tag{3.8} $$2. The initial displacement $u_0$ is selected
3. In the cycle from $i=0$ to the selected value of the maximum allowable number of iterations $i_{max}$
a) the following displacement value is calculated
$$ u_{i+1} = u_i - \frac{R(u_i)}{K_T(u_i)} \tag{3.9} $$b) the absolute value of the non-equilibrium force $|R_{i+1}| = |R(u_{i+1})|$ is compared with the prescribed tolerance $R_{tol}$ and the decision-making steps are taken
• if $|R_{i+1}| \ge R_{tol}$, go to 3, if not, finish the solution
• if $i > i_{max}$, finish the solution
For the bar fastened and loaded according to Fig. 3.3, the vertical displacement of the end point of the bar $u_A$ is calculated using the Newton-Raphson method with a permissible solution error of $R_{tol} = 1$ N. Given: $F = 4000$ N, $l_0 = 2000$ mm, $h = 50$ mm, $E_0 = 200\,000$ MPa, and $S_0 = 800$ mm².
The relationship between the displacement of the bar end point and the force $F$ in the form of an non-equilibrium force is given by Equation (3.4)
$$ R(u) = 0.01(u^3+150u^2+5000u)-4000 $$with tangential stiffness $K_T(u)=\frac{dR(u)}{du}=0.01(3u^2+300u+5000)$
The first step (i=0): The initial displacement (selected) $u_0=0$
The first displacement value is $u_1 = u_0 - \frac{R(0)}{K_T(0)} = 0 - \frac{-4000}{50} = 80$ mm
The absolute value of the non-equilibrium force at this offset value was (3.4) $|R_1| = 14720$ N
Second step (i=1): The corrected displacement value is $u_2 = u_1 - \frac{R(80)}{K_T(80)} = 80 - \frac{14720}{482} = 49.46$ mm
The non-equilibrium force is $|R_2| = 3352.5$ N
The third step (i=2): The corrected displacement value is $u_3 = u_2 - \frac{R(49.46)}{K_T(49.46)} = 49.46 - \frac{3352.53}{271.77} = 37.12$ mm
The non-equilibrium force dropped to $|R_3| = 435.28$ N.
To calculate the next values we used the program Mathematica 7 shown in Fig. 3.4, which also shows the calculated values of the example in the next steps. A graphical illustration of the iterative solution process is presented in Fig. 3.5.
From the result of the iteration it can be seen that a sufficiently accurate solution of equation (3.2) is the displacement with the value $u_A = 34.92$ mm.
For the general use of the method it is useful to recall its mathematical principle: the function $r(u)$ can be approximated (linearized) around the point $u_i$ by decomposition into a truncated Taylor series
$$ r(u) \approx r(u_i) + \frac{dr(u_i)}{du}(u-u_i) \tag{3.10} $$The solution of the equation can then be obtained by the successive iterative solution of the equation
$$ r(u_i) + \frac{dr(u_i)}{du}\Delta u = 0 \quad i=0,1,2,... $$ $$ u_{i+1} = u_i + \Delta u \tag{3.11} $$to the extent of accuracy specified by the selected convergence condition (the convergence is not guaranteed for any type of function $r$).
The Newton-Raphson method can also be generalized for a system of nonlinear equations with several unknowns. This modification is normally used in FEM because the application of this method to a nonlinear problem leads to a system of such equations. The generalized displacements of the nodal points arranged in vector of functions $\mathbf{u} = [u_1,u_2,u_3,...,u_n]^T$ are unknown. In matrix notation, the system of equations of body equilibrium has the following shape
$$ \mathbf{r}(\mathbf{u}) = \mathbf{q}(\mathbf{u}) - \mathbf{f} = \mathbf{0} \tag{3.12} $$where $\mathbf{f}$ is the known external nodal force vector of the body and $\mathbf{q}(\mathbf{u})$ is the internal nodal force vector of the FEM model calculated from finite element stresses containing nonlinear functions of unknowns $\mathbf{u}$.
Equations (3.12) express that the body is in equilibrium when the internal nodal forces are equal to the external nodal forces at each node (in the nodes without external forces, the internal nodal forces are in equilibrium). Using FEM, we create a system (3.12) and then, based on the selected initial values, iteratively adjust the displacements so that the non-equilibrium forces in the vector drop to sufficiently small values according to a certain prescribed convergence condition.
The Newton-Raphson method for the system of nonlinear equations $\mathbf{r}(\mathbf{u}) = \mathbf{0}$, analogous to the one-dimensional case, consists of linearizing the components of the vector function $\mathbf{r}(\mathbf{u}) = [f_1(\mathbf{u}),f_2(\mathbf{u}),...,f_n(\mathbf{u})]^T$ using their first differential and replacing the nonlinear system of equations with a linear system to correct the arguments by $\Delta\mathbf{u}$.
By linearizing the system $\mathbf{r}(\mathbf{u}) = \mathbf{0}$ around the initial approximation $\mathbf{u}_0$ we get
$$ \mathbf{r}_L(\mathbf{u}) = \mathbf{r}(\mathbf{u}_0) + \frac{\partial\mathbf{r}(\mathbf{u}_0)}{\partial\mathbf{u}}(\mathbf{u}-\mathbf{u}_0) $$We are looking for the first approximation $\mathbf{u}_1$ such that $\mathbf{r}_L(\mathbf{u}_1) = \mathbf{0}$. Then it applies
$$ \frac{\partial\mathbf{r}}{\partial\mathbf{u}}(\mathbf{u}_0)(\mathbf{u}_1-\mathbf{u}_0) = -\mathbf{r}(\mathbf{u}_0) \tag{3.13} $$After selecting the initial displacement values $\mathbf{u}_0$, (3.13) represents a linear system of equations with an unknown vector $\Delta\mathbf{u}_0 = \mathbf{u}_1 - \mathbf{u}_0$. If the matrix of this system is regular then its solution is
$$ \Delta\mathbf{u}_0 = -[\mathbf{K}_T(\mathbf{u}_0)]^{-1}\mathbf{r}(\mathbf{u}_0) $$with the so-called tangential body stiffness matrix (Jacobi matrix)
$$ \mathbf{K}_T = \frac{\partial\mathbf{r}}{\partial\mathbf{u}} = \begin{bmatrix} \partial R_1/\partial u_1 & \dots & \partial R_1/\partial u_n \\ \vdots & \ddots & \vdots \\ \partial R_n/\partial u_1 & \dots & \partial R_n/\partial u_n \end{bmatrix} \tag{3.14} $$The desired approximation is then given by $\mathbf{u}_1 = \mathbf{u}_0 + \Delta\mathbf{u}_0$. By repeating this procedure, we obtain the relation for the sequence of approximations
$$ \mathbf{u}_{i+1} = \mathbf{u}_i - [\mathbf{K}_T(\mathbf{u}_i)]^{-1}\mathbf{r}(\mathbf{u}_i) \quad i=0,1,2,...,i_{max} \tag{3.15} $$If this sequence converges to the solution of the system of equations (3.12), the calculation is terminated by the condition for $\|\mathbf{u}_{i+1}-\mathbf{u}_i\|$ or $\|\mathbf{r}_{i+1}-\mathbf{r}_i\|$.
In cases where the external nodal forces in (3.12) are not dependent on displacements, the tangential stiffness matrix can be expressed from the internal nodal forces as follows:
$$ \mathbf{K}_T = \frac{\partial\mathbf{r}}{\partial\mathbf{u}} = \frac{\partial\mathbf{q}}{\partial\mathbf{u}} \tag{3.16} $$The procedure, as shown in (3.15), requires inversion of the tangential body stiffness matrix at each iteration step, which can be time-consuming for large problems. However, in this case, the quadratic convergence of the iterative process is guaranteed and we are talking about the full Newton-Raphson method. In FEM programs it is usually possible to choose the modified Newton-Raphson method, where the inversion of the tangential stiffness matrix is performed only once at the beginning of the iteration process and does not change in the next iteration steps. Thus, the same matrix was used in all iteration steps. The iterative process is accelerated, but the number of iterative steps increases as the rate of convergence decreases.
Equation (3.15) can be adjusted to
$$ \Delta\mathbf{u}_{i+1} = \mathbf{u}_{i+1}-\mathbf{u}_i = -[\mathbf{K}_T(\mathbf{u}_i)]^{-1}(\mathbf{q}_i-\mathbf{f}) \tag{3.17} $$and further multiplied from the left by $\mathbf{K}_T(\mathbf{u}_i)$ to the form in which it usually occurs in program manuals
$$ \mathbf{K}_T(\mathbf{u}_i)\Delta\mathbf{u}_{i+1} = \mathbf{f}-\mathbf{q}_i \tag{3.18} $$Suppose that the iterative procedure begins with an unloaded body. Then, the node displacements are zero ($\mathbf{u}=0 \rightarrow \mathbf{K}_T=\mathbf{K}_{lin}$) and since the stresses are zero ($\mathbf{q}=0$) the initial iteration step is performed with a linear system $\mathbf{K}_{lin}\mathbf{u}_1 = \mathbf{f}$. From the calculated linear solution $\mathbf{u}_1$, the deformations and stresses in the elements are determined, from which $\mathbf{q}_1$ and the vector of non-equilibrium nodal forces $\Delta\mathbf{r}_1 = \mathbf{f}-\mathbf{q}_1$ are determined. A tangential stiffness matrix is also created for the next iteration step. According to (3.17), the displacement corrections $\Delta\mathbf{u}_{i+1}$ are then calculated $$ \mathbf{u}_{i+1} = \mathbf{u}_i + \Delta\mathbf{u}_{i+1} $$
With the new displacements, the iteration process continued until the convergence conditions were met. In the case of a linear problem, convergence is achieved in the first iteration.
In the case where the body and its loading form a so-called conservative system (meaning that no energy is dissipated during loading and the forces are conservative), it is also possible to formulate the equilibrium equations based on the principle of minimum total potential energy of the body, for which the following holds
$$ \Pi = U - W = \text{minimum} \tag{3.19} $$where $U$ is the strain energy of the body and $W$ is the potential energy of the external forces, referenced to the state of the body before deformation (which implies its negative sign, and the products of the external forces with displacements are not multiplied by one-half). The vector of unbalanced forces can then be expressed as the gradient of the extremum of this energy with respect to the displacement vector of the body’s nodal points
$$ \mathbf{r} = \frac{\partial\Pi}{\partial\mathbf{u}} = \mathbf{0} \tag{3.20} $$According to (3.20), (3.19), and (3.12) we then obtain
$$ \mathbf{r} = \frac{\partial U}{\partial\mathbf{u}} - \frac{\partial W}{\partial\mathbf{u}} = \mathbf{q}-\mathbf{f}=0 \tag{3.21} $$For the internal nodal forces, it follows that
$$ \mathbf{q} = \frac{\partial U}{\partial\mathbf{u}} \tag{3.22} $$The application of FEM to solve a simple geometrically non-linear truss system with a single free nodal point with three degrees of freedom led to a system of three nodal equations of equilibrium
$$ R_1(u_1,u_2,u_3) = (0,01u_1^2 + 50)u_1 + 1,5u_1u_2^2 - 1,5u_1u_3^2 - 600 = 0 $$ $$ R_2(u_1,u_2,u_3) = 1,5u_1^2u_2 + (0,01u_2^2 + 50)u_2 + 1,5u_2u_3^2 - 800 = 0 $$ $$ R_3(u_1,u_2,u_3) = -1,5u_1^2u_3 + 1,5u_2^2u_3 + (0,01u_3^2 + 50)u_3 - 500 = 0 $$in matrix notation
$$ \mathbf{r}(\mathbf{u}) = \mathbf{q}(\mathbf{u}) - \mathbf{f} = \mathbf{0} $$where vector $\mathbf{f}$ contains the orthogonal components of the external force acting at the node. The nodal displacement components $\mathbf{u} = [u_1\quad u_2\quad u_3]^T$ are determined using the Newton-Raphson method.
The Mathematica 7 solution is shown in Fig. 3.6, where a set of equations is given in part a), the program for solving the problem is given in part b), the program start command with the necessary input parameters is in part c), and part d) contains an iterative sequence of solutions with the resulting non-equilibrium force vector satisfying the convergence condition.
From the results it can be seen that iteration begins with zero displacement components, the vector $\mathbf{u}_1$ contains a linear solution and finally the vector $\mathbf{u}_8$ provides the values of the rectangular displacement components in the node (in mm) that satisfy the convergence condition.
For the bar fixed and loaded as shown in the figure Fig. 3.3, ANSYS Mechanical APDL was used to calculate the vertical displacement of the bar end point $u_A$ when $F = 4000$ N, $l_0 = 2000$ mm, $h = 50$ mm, $E_0 = 200000$ MPa, and $S_0 = 800$ mm$^2$. See attached Video 1.
The result $u_A = 34.952$ mm is in good agreement with the result in Example 3.1.
[1] Benča, Š.: Solution of nonlinear structural problems using FEM. Slovak university of technology, Bratislava, 2009 (in Slovak)