# Nonlinear finite elements/Nonlinear axial bar weak form

## Symmetric weak form for the axially loaded bar

Let us start with the weak form and derive a symmetric stiffness matrix. We have,

${ A~\cfrac{d}{dx}\left[E(u)\cfrac{du}{dx}\right] + ax = 0 ~. }$

Once again, we multiply with a weighting function and integrate over an element to get

$\int_{\Omega} A~\cfrac{d}{dx}\left[E(u)\cfrac{du}{dx}\right]w~dx + \int_{\Omega} axw = 0 ~.$

Integrating the first term by parts, we get

$\left.AE(u)\cfrac{du}{dx}w\right|_{\Gamma_t} - \int_{\Omega} AE(u)\cfrac{du}{dx}\cfrac{dw}{dx}~dx + \int_{\Omega} axw = 0 ~.$

Taking the internal force terms to the left and the external force terms to the right, we get

$\int_{\Omega} AE(u)\cfrac{du}{dx}\cfrac{dw}{dx}~dx = \int_{\Omega} axw + \left.AE(u)\cfrac{du}{dx}w\right|_{\Gamma_t} ~.$

Note that, at $x = L$,

$f = A\sigma = AE(u)\cfrac{du}{dx} = R ~.$

Also, at $x = 0$, we have $w = 0$. The external surface force is zero at all other points. Therefore,

$\int_{\Omega} AE(u)\cfrac{du}{dx}\cfrac{dw}{dx}~dx = \int_{\Omega} axw + Rw|_{x=L} ~.$

## The finite element approximation

Let, the trial and weighting functions be

$u = \sum_j u_j N_j~;~\qquad w = N_i ~.$

Then, we get

$\int_{\Omega} AE(u) \left(\sum_j u_j \cfrac{dN_j}{dx}\right) \cfrac{dN_i}{dx}~dx = \int_{\Omega} axN_i + RN_i|_{x=L} ~.$

Take the constants outside the integral sign as usual, and we can write

$A\sum_j\left(\int_{\Omega} E(u) \cfrac{dN_i}{dx}\cfrac{dN_j}{dx}~dx\right) u_j = \int_{\Omega} axN_i + RN_i|_{x=L} ~.$

The above can be written in the usual form as

$K_{ij} u_j = f_i \,$

where

{ \begin{align} K_{ij} & := A\int_{\Omega} E(u) \cfrac{dN_i}{dx}\cfrac{dN_j}{dx}~dx \\ f_i & := \int_{\Omega} axN_i + RN_i|_{x=L} ~. \end{align} }

This stiffness matrix is symmetric. However, we are still left with the problem of dealing with the $E(u)$ part.

### Dealing with $E(u)$

The expression of $E(u)$ is

$E(u) := E_0\left(1-\cfrac{du}{dx}\right)~.$

In terms of the trial function, we have

$E(u) = E_0\left(1-\sum_k u_k\cfrac{dN_k}{dx}\right)~.$

I have chosen to call the index $k$ in order to distinguish it from the index names that we have already used, that is $i$ and $j$.

Plugging the expansion for $E(u)$ into the expression for $K_{ij}$ we get

$K_{ij} = AE_0\int_{\Omega}\left(1-\sum_k u_k\cfrac{dN_k}{dx}\right) \cfrac{dN_i}{dx}\cfrac{dN_j}{dx}~dx ~.$

Notice that the stiffness matrix continues to remain symmetric when we do this because the sum over $k$ includes all the relevant nodes.

The stiffness matrix is now a function of the unknown displacements $u_i$ and the system of equations can be written as

$\mathbf{K}(\mathbf{u}) \mathbf{u} = \mathbf{f}~.$

This is a nonlinear system of equations and cannot be solved directly (as was done for linear systems).

### Making things a bit more concrete

Let us get a bit more explicit about what all that leads to.

Assume that the bar is divided into two linear elements. Therefore, there are three global nodes: $1$, $2$, and $3$. Each element has two nodes which have local numbers $1$ and $2$.

The shape functions within each element are

$N_1 = \cfrac{x_2 - x}{h}, \qquad N_2 = \cfrac{x - x_1}{h}$

where $x_1$ is the coordinate of node $1$ of the element, and $x_2$ is the coordinate of node $2$ of the element. The length of the element is $h$.

The derivatives of the shape functions are

$\cfrac{dN_1}{dx} = -\cfrac{1}{h}, \qquad \cfrac{dN_2}{dx} = \cfrac{1}{h}~.$

The element stiffness matrix terms can be written as

$K_{ij} = AE_0\int_{x_1}^{x_2} \left(1-u_1\cfrac{dN_1}{dx}-u_2\cfrac{dN_2}{dx}\right) \cfrac{dN_i}{dx}\cfrac{dN_j}{dx}~dx$

where $u_1$ is the displacement of local node $1$, and $u_2$ is the displacement of local node $2$.

Plugging in the values of the derivatives of the shape functions, we get

\begin{align} K_{11} & = \cfrac{AE_0}{h^2} \int_{x_1}^{x_2}\left(1+\cfrac{u_1}{h}-\cfrac{u_2}{h}\right)~dx \\ K_{12} & = -\cfrac{AE_0}{h^2} \int_{x_1}^{x_2}\left(1+\cfrac{u_1}{h}-\cfrac{u_2}{h}\right)~dx \\ K_{21} & = -\cfrac{AE_0}{h^2} \int_{x_1}^{x_2}\left(1+\cfrac{u_1}{h}-\cfrac{u_2}{h}\right)~dx \\ K_{22} & = \cfrac{AE_0}{h^2} \int_{x_1}^{x_2}\left(1+\cfrac{u_1}{h}-\cfrac{u_2}{h}\right)~dx \end{align}

The terms inside the integrals are the same for all the coefficients of the stiffness matrix (this is not true in general), and we have

$\int_{x_1}^{x_2}\left(1+\cfrac{u_1}{h}-\cfrac{u_2}{h}\right)~dx = h + u_1 - u_2 ~.$

Therefore,

\begin{align} K_{11} & = \cfrac{AE_0}{h^2}(h + u_1 - u_2) \\ K_{12} & = -\cfrac{AE_0}{h^2}(h + u_1 - u_2) \\ K_{21} & = -\cfrac{AE_0}{h^2}(h + u_1 - u_2) \\ K_{22} & = \cfrac{AE_0}{h^2}(h + u_1 - u_2) \end{align}

In matrix form, the element stiffness matrix is

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

For simplicity, let us set the body force term $a$ to zero. Let us also assume that both elements have equal lengths. Then the assembled global system of equations (for the two element mesh) is

$\cfrac{AE_0}{h^2} \begin{bmatrix} h+u_1-u_2 & -(h+u_1-u_2) & 0 \\ -(h+u_1-u_2) & (h+u_1-u_2)+(h+u_2-u_3) & -(h+u_2-u_3) \\ 0 & -(h+u_2-u_3) & (h+u_2-u_3) \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} f_1 \\ 0 \\ R \end{bmatrix} ~.$

After simplification, we have

$\cfrac{AE_0}{h^2} \begin{bmatrix} h+u_1-u_2 & -(h+u_1-u_2) & 0 \\ -(h+u_1-u_2) & (2h+u_1-u_3) & -(h+u_2-u_3) \\ 0 & -(h+u_2-u_3) & (h+u_2-u_3) \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} f_1 \\ 0 \\ R \end{bmatrix} ~.$

The stiffness matrix is a function of $u_1$, $u_2$, and $u_3$.

The boundary condition at $x = 0$ is $u_1 = 0$. Therefore, the reduced system of equations can be written as

$\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 reduced stiffness matrix is a function of $u_2$ and $u_3$.

How can we find $u_2$ and $u_3$ under these circumstances? We will use the Newton-Raphson numerical technique.