Nonlinear finite elements/Newton method for axial bar

Newton's method applied to the axial bar problem
Let us now return to the axial bar problem and see how the Newton-Raphson scheme can be used to solve the system of equations. The reduced system of equations is

\mathbf{K} \mathbf{u} = \mathbf{f} \qquad \rightarrow \qquad \cfrac{AE_0}{h^2} \begin{bmatrix} (2h-u_3) & -(h+u_2-u_3) \\ -(h+u_2-u_3) & (h+u_2-u_3) \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} 0 \\ R \end{bmatrix} ~. $$ Therefore, the residual is

\mathbf{r} = \mathbf{K} \mathbf{u} - \mathbf{f} \qquad \rightarrow \qquad \mathbf{r} = \cfrac{AE_0}{h^2} \begin{bmatrix} (2h-u_3) & -(h+u_2-u_3) \\ -(h+u_2-u_3) & (h+u_2-u_3) \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} - \begin{bmatrix} 0 \\ R \end{bmatrix} ~. $$ The terms of the tangent stiffness matrix are given by

\mathbf{T} := \mathbf{K} + \frac{\partial \mathbf{K}}{\partial \mathbf{u}}~\mathbf{u} \qquad \rightarrow \qquad T_{ij} := K_{ij} + \sum_{m=1}^n \frac{\partial K_{im}}{\partial u_j} u_m ~. $$ For our $$2\times 2$$ reduced stiffness matrix, we have
 * $$\begin{align}

T_{22} & = K_{22} + \frac{\partial K_{22}}{\partial u_2} u_2 + \frac{\partial K_{23}}{\partial u_2} u_3 \\ T_{23} & = K_{23} + \frac{\partial K_{22}}{\partial u_3} u_2 + \frac{\partial K_{23}}{\partial u_3} u_3 \\ T_{32} & = K_{32} + \frac{\partial K_{32}}{\partial u_2} u_2 + \frac{\partial K_{33}}{\partial u_2} u_3 \\ T_{33} & = K_{33} + \frac{\partial K_{32}}{\partial u_3} u_2 + \frac{\partial K_{33}}{\partial u_3} u_3~. \end{align}$$ Now, from the stiffness matrix, we see that
 * $$\begin{align}

K_{22} & = \cfrac{AE_0}{h^2} (2h - u_3) & K_{23} &= -\cfrac{AE_0}{h^2} (h + u_2 - u_3) \\ K_{32} & = -\cfrac{AE_0}{h^2} (h + u_2 - u_3) & K_{33} &= \cfrac{AE_0}{h^2} (h + u_2 - u_3) ~. \end{align}$$ The derivatives are,
 * $$\begin{align}

\frac{\partial K_{22}}{\partial u_2} & = 0 & \frac{\partial K_{23}}{\partial u_2} &= -\cfrac{AE_0}{h^2} & \frac{\partial K_{22}}{\partial u_3} & = -\cfrac{AE_0}{h^2} & \frac{\partial K_{23}}{\partial u_3} &= \cfrac{AE_0}{h^2} \\ \frac{\partial K_{32}}{\partial u_2} &= -\cfrac{AE_0}{h^2} & \frac{\partial K_{33}}{\partial u_2} &= \cfrac{AE_0}{h^2} & \frac{\partial K_{32}}{\partial u_3} &= \cfrac{AE_0}{h^2} & \frac{\partial K_{33}}{\partial u_3} &= -\cfrac{AE_0}{h^2} ~. \end{align}$$ Therefore, the components of $$\mathbf{T}$$ are
 * $$\begin{align}

T_{22} & = \cfrac{AE_0}{h^2} (2h - u_3) + 0 - \cfrac{AE_0}{h^2} u_3 & T_{23} & = -\cfrac{AE_0}{h^2} (h + u_2 - u_3) - \cfrac{AE_0}{h^2} u_2 + \cfrac{AE_0}{h^2} u_3 \\ T_{32} & = -\cfrac{AE_0}{h^2} (h + u_2 - u_3) - \cfrac{AE_0}{h^2} u_2 + \cfrac{AE_0}{h^2} u_3 & T_{33} & = \cfrac{AE_0}{h^2} (h + u_2 - u_3) + \cfrac{AE_0}{h^2} u_2 - \cfrac{AE_0}{h^2} u_3~. \end{align}$$ In matrix form,
 * $$ \text{(8)} \qquad

{   \mathbf{T} = \cfrac{AE_0}{h^2} \begin{bmatrix} 2(h - u_3)        & -h - 2(u_2 - u_3) \\ -h - 2(u_2 - u_3) &  h + 2(u_2 - u_3) \end{bmatrix}~. } $$ The inverse of $$\mathbf{T}$$ (the reduced form shown in equation (8)) is

\mathbf{T}^{-1} = \cfrac{h^2}{AE_0(h - 2u_2)} \begin{bmatrix} 1 & 1 \\          1 & \cfrac{2h - 2u_3}{h + 2u_2 - 2u_3} \end{bmatrix} ~. $$

Maple code
11 A Maple text script for doing these computations symbolically is show below.

Can the tangent stiffness matrix be assembled too?
We have derived the tangent stiffness matrix $$\mathbf{T}$$ from the assembled global stiffness matrix. If the stiffness matrix is large that does not seem to be the way to go. So the question is, can the tangent stiffness matrix be calculated on an element by element basis and assembled just like the global stiffness matrix?

Let us try to figure that out.

The element stiffness matrix has the form

\mathbf{K}^e = \cfrac{AE_0}{h^2}(h + u_1 - u_2)\begin{bmatrix}1 & -1 \\-1 & 1 \end{bmatrix} ~. $$ The tangent stiffness matrix is given by

\mathbf{T}^e = \mathbf{K}^e + \frac{\partial \mathbf{K}^e}{\partial \mathbf{u}}~\mathbf{u}~. $$ After going through the algebra, we get

\mathbf{T}^e = \cfrac{AE_0}{h^2}(h + 2u_1 - 2u_2)\begin{bmatrix}1 & -1 \\-1 & 1 \end{bmatrix} ~. $$ For the two element mesh, the assembled global tangent stiffness matrix is

\mathbf{T} = \cfrac{AE_0}{h^2} \begin{bmatrix} h + 2(u_1 - u_2) & -h - 2(u_1 - u_2) & 0 \\ & 2h + 2(u_1 - u_3) & -h - 2(u_2 - u_3)\\ \text{Symm.}    &                   & h + 2(u_2 - u_3) \\ \end{bmatrix} $$ Applying the boundary condition $$u_1 = 0$$, we get

\mathbf{T} = \cfrac{AE_0}{h^2} \begin{bmatrix} h - 2u_2    & -h + 2u_2  & 0 \\ & 2(h - u_3) & -h - 2(u_2 - u_3)\\ \text{Symm.} &           & h + 2(u_2 - u_3) \end{bmatrix} $$ This matrix is the same as that shown in equation (8). Therefore, we can actually assemble the global tangent stiffness matrix. You could also show that in more general terms. But we avoid that proof here.

Doing the actual computation
The Newton iteration scheme is

\mathbf{u}_{r+1} = \mathbf{u}_r - \mathbf{T}^{-1}(\mathbf{u}_r)~\mathbf{r}(\mathbf{u}_r) $$

Let us go through the process. We will have to start with some assumed values of the parameters.

Let $$A = 0.01$$ m$$^2$$, $$L = 2$$ m, $$E_0 = 100$$ MPa, $$R = 100$$ kN. Since there are two elements, $$h = 1$$ m.

For the first Newton iteration, let the initial guess be

\mathbf{u}_0 = \begin{bmatrix} u_2^0 \\ u_3^0 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix} $$ Then, the first Newton iteration gives

\mathbf{K}_0 = 10^6\begin{bmatrix} 2 & -1 \\ -1 & 1 \end{bmatrix}, \qquad \mathbf{T}_0 = 10^6\begin{bmatrix} 2 & -1 \\ -1 & 1 \end{bmatrix}, \qquad \mathbf{u}_1 = \begin{bmatrix} 0.1 \\ 0.2 \end{bmatrix} ~. $$ The second Newton iteration gives

\mathbf{K}_1 = 10^6\begin{bmatrix} 1.8 & -0.9 \\ -0.9 & 0.9 \end{bmatrix}, \qquad \mathbf{T}_1 = 10^6\begin{bmatrix} 1.6 & -0.8 \\ -0.8 & 0.8 \end{bmatrix}, \qquad \mathbf{u}_2 = \begin{bmatrix} 0.1125 \\ 0.2250 \end{bmatrix} ~. $$ The third Newton iteration gives

\mathbf{K}_2 = 10^6\begin{bmatrix}1.775 & -0.8875\\-0.8875 & 0.8875 \end{bmatrix}, \qquad \mathbf{T}_2 = 10^6\begin{bmatrix}1.55 & -0.775\\-0.775 & 0.775 \end{bmatrix}, \qquad \mathbf{u}_3 = \begin{bmatrix} 0.1127 \\ 0.2254 \end{bmatrix} ~. $$ The fourth Newton iteration gives

\mathbf{K}_3 = 10^6\begin{bmatrix}1.7746&-0.8873\\-0.8873 & 0.8873 \end{bmatrix}, \qquad \mathbf{T}_3 = 10^6\begin{bmatrix}1.5492&-0.7746\\-0.7746 & 0.7746 \end{bmatrix}, \qquad \mathbf{u}_4 = \begin{bmatrix} 0.1127 \\ 0.2254 \end{bmatrix} ~. $$ At this stage we stop the computation because the solution has converged. A plot of the convergence of the Newton iterations in shown in Figure 1. Once we have computed the nodal displacements, we can compute the element stresses in the usual manner. We write the strains as

\varepsilon = \cfrac{du}{dx} = u_1\cfrac{dN_1}{dx} + u_2 \cfrac{dN_2}{dx}~. $$ The stresses are computed using

\sigma = E_0(1-\varepsilon)\varepsilon ~. $$ For the axially loaded bar with no body forces, we have a uniform state of strain (and therefore stress). For the numbers that we have used, these stresses are 10 MPa in each element.

A slightly more complicated situation
At this stage, let us add in the body force terms and work out the solution for another test case.

Let us assume that the input parameters are $$E_0 = 10$$, $$A = 1$$, $$L = 1$$, $$a = 1$$, and $$R = 2$$. The stress-strain curve for this material is shown in Figure 2. We use an error tolerance of $$10^{-6}$$ as a flag to stop the Newton iterations.

Mesh with 10 elements
If we use a mesh containing $$10$$ equally sized elements, we get the solutions shown in Figure 3.

Mesh with 100 elements
If we use a mesh containing $$100$$ equally sized elements, we get the solutions shown in Figure 4.

Even though the displacements and strains are quite different in the linear and the nonlinear cases, the stresses are remarkably exactly the same for the linear and nonlinear cases. This is a result of the form of the modulus that we have chosen to use for this problem. Notice that the error norms show nearly quadratic convergence.

There is however one issue. Due to the shape of the stress-strain curve of the material, Newton iterations do not converge when the applied external load exceeds the amount needed for the strains to exceed 0.5. A method other than Newton-Raphson is needed to solve the problem for such situations.

The Matlab script used for the above calculations is shown below.