Nonlinear finite elements/Euler Bernoulli beams

Euler-Bernoulli Beam

Euler-Bernoulli beam

Displacements

$\begin{array}{cc}{u}_1& =u_0(x)-z\frac{d{w}_0}{dx}\\ u_2& =0\\ u_3& =w_0(x)\end{array}$

Strains

${\epsilon }_{11}={\epsilon }_{xx}={\epsilon }_{xx}^0+z{\epsilon }_{xx}^1$

$\begin{array}{cc}{\epsilon }_{xx}^0& =\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2\\ {\epsilon }_{xx}^1& =-\frac{d^2{w}_0}{d{x}^2}\end{array}$

Strain-Displacement Relations

${\epsilon }_{ij}=\frac{1}{2}(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i})+\frac{1}{2}\left(\frac{\partial u_m}{\partial x_i}\frac{\partial u_m}{\partial x_j}\right)$

$\begin{array}{cc}{\epsilon }_{11}& =\frac{1}{2}(\frac{\partial u_1}{\partial x_1}+\frac{\partial u_1}{\partial x_1})+\frac{1}{2}(\frac{\partial u_1}{\partial x_1}\frac{\partial u_1}{\partial x_1}+\frac{\partial u_2}{\partial x_1}\frac{\partial u_2}{\partial x_1}+\frac{\partial u_3}{\partial x_1}\frac{\partial u_3}{\partial x_1})\\ {\epsilon }_{22}& =\frac{1}{2}(\frac{\partial u_2}{\partial x_2}+\frac{\partial u_2}{\partial x_2})+\frac{1}{2}(\frac{\partial u_1}{\partial x_2}\frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_2}\frac{\partial u_2}{\partial x_2}+\frac{\partial u_3}{\partial x_2}\frac{\partial u_3}{\partial x_2})\\ {\epsilon }_{33}& =\frac{1}{2}(\frac{\partial u_3}{\partial x_3}+\frac{\partial u_3}{\partial x_3})+\frac{1}{2}(\frac{\partial u_1}{\partial x_3}\frac{\partial u_1}{\partial x_3}+\frac{\partial u_2}{\partial x_3}\frac{\partial u_2}{\partial x_3}+\frac{\partial u_3}{\partial x_3}\frac{\partial u_3}{\partial x_3})\\ {\epsilon }_{23}& =\frac{1}{2}(\frac{\partial u_2}{\partial x_3}+\frac{\partial u_3}{\partial x_2})+\frac{1}{2}(\frac{\partial u_1}{\partial x_2}\frac{\partial u_1}{\partial x_3}+\frac{\partial u_2}{\partial x_2}\frac{\partial u_2}{\partial x_3}+\frac{\partial u_3}{\partial x_2}\frac{\partial u_3}{\partial x_3})\\ {\epsilon }_{31}& =\frac{1}{2}(\frac{\partial u_3}{\partial x_1}+\frac{\partial u_1}{\partial x_3})+\frac{1}{2}(\frac{\partial u_1}{\partial x_3}\frac{\partial u_1}{\partial x_1}+\frac{\partial u_2}{\partial x_3}\frac{\partial u_2}{\partial x_1}+\frac{\partial u_3}{\partial x_3}\frac{\partial u_3}{\partial x_1})\\ {\epsilon }_{12}& =\frac{1}{2}(\frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1})+\frac{1}{2}(\frac{\partial u_1}{\partial x_1}\frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1}\frac{\partial u_2}{\partial x_2}+\frac{\partial u_3}{\partial x_1}\frac{\partial u_3}{\partial x_2})\end{array}$

The displacements

$u_1=u_0(x_1)-x_3\frac{d{w}_0}{d{x}_1};u_2=0;u_3=w_0(x_1)$

The derivatives

$\begin{array}{cccccc}\frac{\partial u_1}{\partial x_1}& =\frac{d{u}_0}{d{x}_1}-x_3\frac{d^2{w}_0}{d{x}_1^2};& \frac{\partial u_1}{\partial x_2}& =0;& \frac{\partial u_1}{\partial x_3}& =-\frac{d{w}_0}{d{x}_1}\\ \frac{\partial u_2}{\partial x_1}& =0;& \frac{\partial u_2}{\partial x_2}& =0;& \frac{\partial u_2}{\partial x_3}& =0\\ \frac{\partial u_3}{\partial x_1}& =\frac{d{w}_0}{d{x}_1};& \frac{\partial u_3}{\partial x_2}& =0;& \frac{\partial u_3}{\partial x_3}& =0\end{array}$

von Karman strains

The von Karman strains

$\begin{array}{cc}{\epsilon }_{11}& =\frac{d{u}_0}{d{x}_1}-x_3\frac{d^2{w}_0}{d{x}_1^2}+\frac{1}{2}[{(\frac{d{u}_0}{d{x}_1}-x_3\frac{d^2{w}_0}{d{x}_1^2})}^2+{\left(\frac{d{w}_0}{d{x}_1}\right)}^2]\\ {\epsilon }_{22}& =0\\ {\epsilon }_{33}& =\frac{1}{2}{\left(\frac{d{w}_0}{d{x}_1}\right)}^2\\ {\epsilon }_{23}& =0\\ {\epsilon }_{31}& =\frac{1}{2}(\frac{d{w}_0}{d{x}_1}-\frac{d{w}_0}{d{x}_1})-\frac{1}{2}\left[(\frac{d{u}_0}{d{x}_1}-x_3\frac{d^2{w}_0}{d{x}_1^2})\left(\frac{d{w}_0}{d{x}_1}\right)\right]\\ {\epsilon }_{12}& =0\end{array}$

Equilibrium Equations

Balance of forces

$\begin{array}{cc}\frac{d{N}_{xx}}{dx}+f(x)& =0\\ \frac{d^2{M}_{xx}}{d{x}^2}+q(x)+\frac{d}{dx}\left(N_{xx}\frac{d{w}_0}{dx}\right)& =0\end{array}$

Stress Resultants

$\begin{array}{cc}{N}_{xx}& =\underset{A}{\int }{\sigma }_{xx}dA\\ M_{xx}& =\underset{A}{\int }z{\sigma }_{xx}dA\end{array}$

Constitutive Relations

Stress-Strain equation

${\sigma }_{xx}=E{\epsilon }_{xx}$

Stress Resultant - Displacement relations

$\begin{array}{cc}{N}_{xx}& =A_{xx}{\epsilon }_{xx}^0+B_{xx}{\epsilon }_{xx}^1=A_{xx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]-B_{xx}\frac{d^2{w}_0}{d{x}^2}\\ M_{xx}& =B_{xx}{\epsilon }_{xx}^0+D_{xx}{\epsilon }_{xx}^1=B_{xx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]-D_{xx}\frac{d^2{w}_0}{d{x}^2}\end{array}$

Extensional/Bending Stiffness

$\begin{array}{cc}{A}_{xx}& =\underset{A}{\int }EdA\leftarrow \text{extensional stiffness}\\ B_{xx}& =\underset{A}{\int }zEdA\leftarrow \text{extensional-bending stiffness}\\ D_{xx}& =\underset{A}{\int }{z}^{2}EdA\leftarrow \text{bending stiffness}\end{array}$

If $E$ is constant, and $x$-axis passes through centroid

$\begin{array}{cc}{A}_{xx}& =E\underset{A}{\int }dA=EA\\ B_{xx}& =E\underset{A}{\int }zdA=0\\ D_{xx}& =E\underset{A}{\int }{z}^{2}dA=EI\end{array}$

Weak Forms

Axial Equation

$\begin{array}{cc}\underset{x_a}{\overset{x_b}{\int }}\frac{d(\delta u_0)}{dx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]A_{xx}dx& =\underset{x_a}{\overset{x_b}{\int }}(\delta u_0)fdx+\\ & \delta u_0(x_a)Q_1+\delta u_0(x_b)Q_4\end{array}$

where

$\begin{array}{cc}\delta u_0& :=v_1\\ Q_1& :=-N_{xx}(x_a)\\ Q_4& :=N_{xx}(x_b)\end{array}$

Bending Equation

$\begin{array}{cc}\underset{x_a}{\overset{x_b}{\int }}\{\frac{d(\delta w_0)}{dx}& [\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]\frac{d{w}_0}{dx}{A}_{xx}+\frac{d^2(\delta w_0)}{d{x}^2}\left(\frac{d^2{w}_0}{d{x}^2}\right)D_{xx}\}dx=\\ & \underset{x_a}{\overset{x_b}{\int }}(\delta w_0)qdx+\delta w_0(x_a)Q_2+\delta w_0(x_b)Q_5+\delta \theta (x_a)Q_3+\delta \theta (x_b)Q_6.\end{array}$

where

$\begin{array}{cccc}\delta w_0& :=v_2& \delta \theta & :=\frac{d{v}_2}{dx}\\ Q_2& :=-{[\frac{d{M}_{xx}}{dx}+N_{xx}\frac{d{w}_0}{dx}]}_{x_a}& Q_5& :={[\frac{d{M}_{xx}}{dx}+N_{xx}\frac{d{w}_0}{dx}]}_{x_b}\\ Q_3& :=-M_{xx}(x_a)& Q_6& :=M_{xx}(x_b)\end{array}$

Finite Element Model

File:EulerBeamFE.png Finite element model for Euler Bernoulli beam

$\begin{array}{cc}{u}_0(x)& =u_1{\psi }_1(x)+u_2{\psi }_2(x)\\ w_0(x)& =w_1{\varphi }_1(x)+{\theta }_1{\varphi }_2(x)+w_2{\varphi }_3(x)+{\theta }_2{\varphi }_4(x)\end{array}$

where $\theta =-(d{w}_0/dx)$.

Hermite Cubic Shape Functions

File:HermiteShape.png Hermite shape functions for beam

Finite Element Equations

$\left[\begin{array}{ccc}{𝐊}^{11}& ⋮& {𝐊}^{12}\\ & ⋮& \\ {𝐊}^{21}& ⋮& {𝐊}^{22}\end{array}\right]\left[\begin{array}{c}𝐮\\ \\ 𝐝\end{array}\right]=\left[\begin{array}{c}{𝐅}^1\\ \\ {𝐅}^2\end{array}\right]$

where

$\begin{array}{cc}𝐮& =[u_1{u}_2{]}^T\\ 𝐝& =[w_1{\theta }_1{w}_2{\theta }_2{]}^T\end{array}$

$\begin{array}{cccc}{𝐊}^{11}& =2\times 2;& {𝐊}^{12}& =2\times 4\\ {𝐊}^{21}& =4\times 2;& {𝐊}^{22}& =4\times 4\end{array}$

Symmetric Stiffness Matrix

$\begin{array}{cc}{K}_{ij}^{11}& =\underset{x_a}{\overset{x_b}{\int }}{A}_{xx}\frac{d{\psi }_i}{dx}\frac{d{\psi }_j}{dx}dx\\ K_{ij}^{12}& =\frac{1}{2}\underset{x_a}{\overset{x_b}{\int }}\left(A_{xx}\frac{d{w}_0}{dx}\right)\frac{d{\psi }_i}{dx}\frac{d{\varphi }_j}{dx}dx\\ K_{ij}^{21}& =\frac{1}{2}\underset{x_a}{\overset{x_b}{\int }}\left(A_{xx}\frac{d{w}_0}{dx}\right)\frac{d{\varphi }_i}{dx}\frac{d{\psi }_j}{dx}dx\\ K_{ij}^{22}& =\underset{x_a}{\overset{x_b}{\int }}\{\frac{1}{2}{A}_{xx}[\frac{d{u}_0}{dx}+{\left(\frac{d{w}_0}{dx}\right)}^2]\frac{d{\varphi }_i}{dx}\frac{d{\varphi }_j}{dx}+D_{xx}\frac{d^2{\varphi }_i}{d{x}^2}\frac{d^2{\varphi }_j}{d{x}^2}\}dx\end{array}$

Load Vector

$\begin{array}{cc}{F}_i^1& =\underset{x_a}{\overset{x_b}{\int }}{\psi }_{i}fdx+{\psi }_i(x_a)Q_1+{\psi }_i(x_b)Q_4\\ F_i^2& =\underset{x_a}{\overset{x_b}{\int }}{\varphi }_{i}qdx+{\varphi }_i(x_a)Q_2+{\varphi }_i(x_b)Q_5+\frac{d{\varphi }_i}{dx}(x_a)Q_3+\frac{d{\varphi }_i}{dx}(x_b)Q_6\end{array}$

Newton-Raphson Solution

$𝐊(𝐔)𝐔=𝐅$

where

$\begin{array}{cc}{U}_1& =u_1,U_2=u_2,U_3=d_1,U_4=d_2,U_5=d_3,U_6=d_4\\ F_1& =F_1^1,F_2=F_2^1,F_3=F_1^2,F_4=F_2^2,F_5=F_3^2,F_6=F_4^2\end{array}$

The residual is

$𝐑=𝐊𝐔-𝐅.$

For Newton iterations, we use the algorithm

${𝐔}^{r+1}={𝐔}^r-({𝐓}^r{)}^{-1}{𝐑}^r$

where the tangent stiffness matrix is given by

${𝐓}^r=\frac{\partial {𝐑}^r}{\partial 𝐔};\text{or}{T}_{ij}=\frac{\partial R_i}{\partial U_j},i=1\dots 6,j=1\dots 6.$

Tangent Stiffness Matrix

$\begin{array}{cc}i=1\dots 2;j=1\dots 2& :\\ & T_{ij}^{11}=K_{ij}^{11}\\ \\ i=1\dots 2;j=1\dots 4& :\\ & T_{ij}^{12}=2K_{ij}^{12}\\ \\ i=1\dots 4;j=1\dots 2& :\\ & T_{ij}^{21}=2K_{ij}^{21}\\ \\ i=1\dots 4;j=1\dots 4& :\\ & T_{ij}^{22}=K_{ij}^{22}+\frac{1}{2}\underset{x_a}{\overset{x_b}{\int }}{A}_{xx}[\frac{d{u}_0}{dx}+2{\left(\frac{d{w}_0}{dx}\right)}^2]\frac{d{\varphi }_i}{dx}\frac{d{\varphi }_j}{dx}dx\end{array}$

Load Steps

Recall

$N_{xx}=A_{xx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]-B_{xx}\frac{d^2{w}_0}{d{x}^2}$

• Divide load into small increments.

$F=\sum _{i=1}^N\mathrm{\Delta }{F}_i$

• Compute $𝐮$ and $𝐝$ for first load step,

$𝐊({𝐔}_0){𝐔}_1=\mathrm{\Delta }{F}_1$

File:EulerBeamStiff.png Stiffness of Euler-Bernoulli beam.

• Compute $𝐮$ and $𝐝$ for second load step,

$𝐊({𝐔}_1){𝐔}_2=\mathrm{\Delta }{F}_1+\mathrm{\Delta }{F}_2$

• Continue until F is reached.

Membrane Locking

Recall

$\left[\begin{array}{ccc}{𝐊}^{11}& ⋮& {𝐊}^{12}\\ & ⋮& \\ {𝐊}^{21}& ⋮& {𝐊}^{22}\end{array}\right]\left[\begin{array}{c}𝐮\\ \\ 𝐝\end{array}\right]=\left[\begin{array}{c}{𝐅}^1\\ \\ {𝐅}^2\end{array}\right]$

where

$\begin{array}{cc}{K}_{ij}^{11}& =\underset{x_a}{\overset{x_b}{\int }}{A}_{xx}\frac{d{\psi }_i}{dx}\frac{d{\psi }_j}{dx}dx\\ K_{ij}^{12}& =\frac{1}{2}\underset{x_a}{\overset{x_b}{\int }}\left(A_{xx}\frac{d{w}_0}{dx}\right)\frac{d{\psi }_i}{dx}\frac{d{\varphi }_j}{dx}dx\\ K_{ij}^{21}& =\frac{1}{2}\underset{x_a}{\overset{x_b}{\int }}\left(A_{xx}\frac{d{w}_0}{dx}\right)\frac{d{\varphi }_i}{dx}\frac{d{\psi }_j}{dx}dx\\ K_{ij}^{22}& =\underset{x_a}{\overset{x_b}{\int }}\{\frac{1}{2}{A}_{xx}[\frac{d{u}_0}{dx}+{\left(\frac{d{w}_0}{dx}\right)}^2]\frac{d{\varphi }_i}{dx}\frac{d{\varphi }_j}{dx}+D_{xx}\frac{d^2{\varphi }_i}{d{x}^2}\frac{d^2{\varphi }_j}{d{x}^2}\}dx\end{array}$

File:EulerBeamMembraneLock.png Mebrane locking in Euler-Bernoulli beam

For Hinged-Hinged

Membrane strain:

${\epsilon }_{xx}^0=0$

or

$\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2=0$

Hence, shape functions should be such that

$\frac{d{u}_0}{dx}\approx {\left(\frac{d{w}_0}{dx}\right)}^2$

$u_0$ linear, $w_0$ cubic $⟹$ Element Locks! Too stiff.

Selective Reduced Integration

• Assume $u_0$ is linear ;~~ $w_0$ is cubic.
• Then $\frac{d{u}_0}{dx}\equiv \frac{d{\psi }_i}{dx}$ is constant, and $\frac{d{w}_0}{dx}\equiv \frac{d{\varphi }_i}{dx}$ is quadratic.
• Try to keep ${\epsilon }_{xx}^0=$ constant.

$\begin{array}{cc}{K}_{ij}^{11}& =\underset{x_a}{\overset{x_b}{\int }}{A}_{xx}\frac{d{\psi }_i}{dx}\frac{d{\psi }_j}{dx}dx\\ K_{ij}^{12}& =\frac{1}{2}\underset{x_a}{\overset{x_b}{\int }}\left(A_{xx}\frac{d{w}_0}{dx}\right)\frac{d{\psi }_i}{dx}\frac{d{\varphi }_j}{dx}dx\\ K_{ij}^{22}& =\underset{x_a}{\overset{x_b}{\int }}\{\frac{1}{2}{A}_{xx}[\frac{d{u}_0}{dx}+{\left(\frac{d{w}_0}{dx}\right)}^2]\frac{d{\varphi }_i}{dx}\frac{d{\varphi }_j}{dx}+D_{xx}\frac{d^2{\varphi }_i}{d{x}^2}\frac{d^2{\varphi }_j}{d{x}^2}\}dx\end{array}$

• ${𝐊}^{11}$ integrand is constant, ${𝐊}^{12}$ integrand is fourth-order , ${𝐊}^{22}$ integrand is eighth-order

Full integration

$n_{\text{gauss pt}}=\text{int}[(p+1)/2]+1$

Assume $A_{xx}$ = constant.

$\begin{array}{cc}{K}_{ij}^{11}& =A_{xx}\underset{x_a}{\overset{x_b}{\int }}\frac{d{\psi }_i}{dx}\frac{d{\psi }_j}{dx}dx\\ & =A_{xx}\underset{-1}{\overset{1}{\int }}\left(J^{-1}\frac{d{\psi }_i(\xi )}{d\xi }\right)\left(J^{-1}\frac{d{\psi }_j(\xi )}{d\xi }\right)Jd\xi =A_{xx}\underset{-1}{\overset{1}{\int }}F(\xi )d\xi \\ & \approx A_{xx}{W}_{1}F({\xi }_1)\leftarrow \text{one-point integration}\end{array}$

$\begin{array}{cc}{K}_{ij}^{12}& =\frac{1}{2}\underset{x_a}{\overset{x_b}{\int }}\left(A_{xx}\frac{d{w}_0}{dx}\right)\frac{d{\psi }_i}{dx}\frac{d{\varphi }_j}{dx}dx\\ & =\frac{A_{xx}}{2}\underset{-1}{\overset{1}{\int }}\left(\sum _{i=1}^4{w}_i{J}^{-1}\frac{d{\varphi }_i(\xi )}{d\xi }\right)\left(J^{-1}\frac{d{\psi }_i(\xi )}{d\xi }\right)\left(J^{-1}\frac{d{\varphi }_j(\xi )}{d\xi }\right)dx\\ & \approx A_{xx}[W_{1}F({\xi }_1)+W_{2}F({\xi }_2)+W_{3}F({\xi }_3)]\leftarrow \text{full integration}\\ & \approx A_{xx}[W_{1}F({\xi }_1)+W_{2}F({\xi }_2)]\leftarrow \text{reduced integration}\end{array}$