Nonlinear finite elements/Nonlinear axial bar weak form

From Wikiversity

Jump to: navigation, search

Contents

[edit] 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} ~.

[edit] 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.

[edit] 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 Kij 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 ui 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).

[edit] 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 x1 is the coordinate of node 1 of the element, and x2 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 u1 is the displacement of local node 1, and u2 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 u1, u2, and u3.

The boundary condition at x = 0 is u1 = 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 u2 and u3.

How can we find u2 and u3 under these circumstances? We will use the Newton-Raphson numerical technique.

Retrieved from "http://en.wikiversity.org/wiki/Nonlinear_finite_elements/Nonlinear_axial_bar_weak_form"