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.

This article is issued from Wikiversity - version of the Wednesday, March 10, 2010. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.