User:Bdhasp/sandbox

Nonlinear Axial Bar
We describe the deformation of a nonlinear elastic bar. The non-linearity to be discussed is

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.

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} -           u_j\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}(u_j) = AE_0 \int_{\Omega} \left[\cfrac{dN_i}{dx}\cfrac{dN_j}{dx} - u_j\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}$$. Also, we have not been able to remove the dependence of the stiffness matrix on the displacement.

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.