Nonlinear finite elements/Nonlinear axially loaded bar

An Example: Axial Bar with Distributed Load

Recall the bar with a distributed axial body load from Handout 8 (see Figure 1).

Figure 1. Distributed axial loading of a bar.

In our previous discussion we had assumed that the bar had a constant Young's modulus of $E$ , that is, the constitutive model of the bar was linear elastic.

Let us now assume that the stress-strain relationship is nonlinear, and of the form

{ \sigma = E(u)~\varepsilon, ~\qquad \varepsilon := \cfrac{du}{dx}, ~\qquad~ E(u) := E_0\left(1-\cfrac{du}{dx}\right) }

where $E_0$ is a material constant which can be interpreted as the initial Young's modulus. Figure 2 shows how the stress-strain relationship for this material looks.

Figure 2. Stress-strain relationship for the nonlinear bar.

The governing differential equation for the bar is

{ A~\cfrac{d\sigma}{dx} + ax = 0 ~. }

If we plug in the stress-strain relation into the governing equation, we get

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

Notice that we have not included the expression for $E(u)$ in the above equation. The reason is that we want to maintain the symmetry of the stiffness matrix.

Effect of including expression for $E(u)$

You can see what happens when we include the expression for $E(u)$ in the steps below. We get,

A~\cfrac{d}{dx}\left[E_0\left(1-\cfrac{du}{dx}\right) \cfrac{du}{dx}\right] + ax = 0 ~.

or,

AE_0~\cfrac{d^2u}{dx^2} - 2AE_0~\cfrac{d^2u}{dx^2}~\cfrac{du}{dx} + ax = 0~.

Now, when we try to derive the weak form, we get

\int_{\Omega} AE_0~\cfrac{d^2u}{dx^2}~w~dx - \int_{\Omega} 2AE_0~\cfrac{d^2u}{dx^2}~\cfrac{du}{dx}~w~dx + \int_{\Omega} ax~w~dx = 0~.

The first integral is the same as before, and we get

\int_{\Omega} AE_0~\cfrac{d^2u}{dx^2}~w~dx = \left. AE_0\cfrac{du}{dx}~w\right|_{\Gamma_t} - \int_{\Omega} AE_0~\cfrac{du}{dx}~\cfrac{dw}{dx}~dx ~.

The third integral remains the same. However, the second integral becomes

\begin{align} \int_{\Omega} 2AE_0~\cfrac{d^2u}{dx^2}~\cfrac{du}{dx}~w~dx & = \left. 2AE_0\cfrac{du}{dx}~w~\cfrac{du}{dx}\right|_{\Gamma_t} - \int_{\Omega} 2AE_0~\cfrac{du}{dx}~\cfrac{d}{dx}\left(w~\cfrac{du}{dx}\right)~dx \\ & = \left. 2AE_0\left(\cfrac{du}{dx}\right)^2w\right|_{\Gamma_t} - \int_{\Omega} 2AE_0~\cfrac{du}{dx}~\left(\cfrac{dw}{dx}\cfrac{du}{dx} + w~\cfrac{d^2u}{dx^2}\right)~dx \\ & = \left. 2AE_0\left(\cfrac{du}{dx}\right)^2w\right|_{\Gamma_t} - \int_{\Omega} 2AE_0~\left(\cfrac{du}{dx}\right)^2\cfrac{dw}{dx}~dx - \int_{\Omega} 2AE_0~\cfrac{d^2u}{dx^2}~\cfrac{du}{dx}~w~dx ~. \end{align}

Rearranging, we get

\int_{\Omega} 2AE_0~\cfrac{d^2u}{dx^2}~\cfrac{du}{dx}~w~dx = \left. AE_0\left(\cfrac{du}{dx}\right)^2w\right|_{\Gamma_t} - \int_{\Omega} AE_0~\left(\cfrac{du}{dx}\right)^2\cfrac{dw}{dx}~dx~.

After collecting all the terms, the weak form becomes

\left. AE_0\cfrac{du}{dx}~w\right|_{\Gamma_t} - \int_{\Omega} AE_0~\cfrac{du}{dx}~\cfrac{dw}{dx}~dx - \left. AE_0\left(\cfrac{du}{dx}\right)^2w\right|_{\Gamma_t} + \int_{\Omega} AE_0~\left(\cfrac{du}{dx}\right)^2\cfrac{dw}{dx}~dx + \int_{\Omega} ax~w~dx = 0 ~.

Separating the terms corresponding to the internal and external forces, we get

\int_{\Omega} AE_0~\cfrac{du}{dx}~\cfrac{dw}{dx}~dx - \int_{\Omega} AE_0~\left(\cfrac{du}{dx}\right)^2\cfrac{dw}{dx}~dx = \int_{\Omega} ax~w~dx + \left. AE_0\left[\cfrac{du}{dx} - \left(\cfrac{du}{dx}\right)^2\right]w\right|_{\Gamma_t}~.

If we choose trial and weighting functions

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

and plug them into the weak form, we get

\int_{\Omega} AE_0\left(\sum_j u_j\cfrac{dN_j}{dx}\right)\cfrac{dN_i}{dx}~dx - \int_{\Omega} AE_0~\left(\sum_j u_j\cfrac{dN_j}{dx}\right)^2 \cfrac{dN_i}{dx}~dx = \int_{\Omega} ax~N_i~dx + \left. AE_0\left[\cfrac{du}{dx} - \left(\cfrac{du}{dx}\right)^2\right]w\right|_{\Gamma_t}~.

Taking the constants out of the integrals, we get (regarding this step, read the Discuss page of this article)

AE_0 \sum_j \left[\int_{\Omega} \left(\cfrac{dN_j}{dx}\cfrac{dN_i}{dx} - \left(\cfrac{dN_j}{dx}\right)^2\cfrac{dN_i}{dx}\right)~dx \right]u_j = \int_{\Omega} ax~N_i~dx + \left. AE_0\left[\cfrac{du}{dx} - \left(\cfrac{du}{dx}\right)^2\right]w\right|_{\Gamma_t}~.

Therefore, the coefficients of the stiffness matrix are given by

K_{ij} = AE_0 \int_{\Omega} \left[\cfrac{dN_i}{dx}\cfrac{dN_j}{dx} - \cfrac{dN_i}{dx}\left(\cfrac{dN_j}{dx}\right)^2\right]~dx ~.

This stiffness matrix is not symmetric because of the second term inside the integral. That is, $K_{ij} \ne K_{ji}$ .

We would like to work with a symmetric stiffness matrix for computational reasons, i.e., the amount of storage needed is smaller and there are a number of very fast solvers for symmetric matrices.

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