Nonlinear finite elements/Timoshenko beams

Timoshenko Beam

File:TimoshenkoBeam.png Timoshenko beam.

Displacements

$\begin{array}{cc}{u}_1& =u_0(x)+z{\phi }_x\\ u_2& =0\\ u_3& =w_0(x)\end{array}$

Strains

$\begin{array}{cc}{\epsilon }_{xx}& ={\epsilon }_{xx}^0+z{\epsilon }_{xx}^1\\ {\gamma }_{xz}& ={\gamma }_{xz}^0\end{array}$

$\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{\phi }_x}{dx}\\ {\gamma }_{xz}^0& ={\phi }_x+\frac{d{w}_0}{dx}\end{array}$

Principle of Virtual Work

$\delta W_{\text{int}}=\delta W_{\text{ext}}$

where

$\begin{array}{cc}\delta W_{\text{int}}& =\underset{x_a}{\overset{x_b}{\int }}\underset{A}{\int }({\sigma }_{xx}\delta {\epsilon }_{xx}+{\sigma }_{xz}\delta {\gamma }_{xz})dAdx\\ & =\underset{x_a}{\overset{x_b}{\int }}(N_{xx}\delta {\epsilon }_{xx}^0+M_{xx}\delta {\epsilon }_{xx}^1+Q_x\delta {\gamma }_{xz}^0)dx\\ \delta W_{\text{ext}}& =\underset{x_a}{\overset{x_b}{\int }}q\delta w_{0}dx+\underset{x_a}{\overset{x_b}{\int }}f\delta u_{0}dx\end{array}$

$Q_x=K_s\underset{A}{\int }{\sigma }_{xz}dA$

$K_s$ = shear correction factor

Taking Variations

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

Take variation

$\delta {\epsilon }_{xx}^0=\frac{d(\delta u_0)}{dx}+\frac{1}{2}\left(2\frac{d{w}_0}{dx}\right)\left(\frac{d(\delta w_0)}{dx}\right)=\frac{d(\delta u_0)}{dx}+\frac{d{w}_0}{dx}\frac{d(\delta w_0)}{dx}.$

${\epsilon }_{xx}^1=\frac{d{\phi }_x}{dx}$

Take variation

$\delta {\epsilon }_{xx}^1=\frac{d\delta {\phi }_x}{dx}$

${\gamma }_{xz}^0={\phi }_x+\frac{d{w}_0}{dx}$

Take variation

$\delta {\gamma }_{xz}^0=\delta {\phi }_x+\frac{d(\delta w_0)}{dx}$

Internal Virtual Work

$\begin{array}{cc}\underset{x_a}{\overset{x_b}{\int }}{N}_{xx}\delta {\epsilon }_{xx}^{0}dx& =\underset{x_a}{\overset{x_b}{\int }}{N}_{xx}[\frac{d(\delta u_0)}{dx}+\frac{d{w}_0}{dx}\frac{d(\delta w_0)}{dx}]dx\\ \underset{x_a}{\overset{x_b}{\int }}{M}_{xx}\delta {\epsilon }_{xx}^{1}dx& =\underset{x_a}{\overset{x_b}{\int }}{M}_{xx}\left[\frac{d\delta {\phi }_x}{dx}\right]dx\\ \underset{x_a}{\overset{x_b}{\int }}{Q}_x\delta {\gamma }_{xz}^{0}dx& =\underset{x_a}{\overset{x_b}{\int }}{Q}_x[\delta {\phi }_x+\frac{d(\delta w_0)}{dx}]dx\end{array}$

$\begin{array}{cc}\delta W_{\text{int}}=\underset{x_a}{\overset{x_b}{\int }}& \{N_{xx}[\frac{d(\delta u_0)}{dx}+\frac{d{w}_0}{dx}\frac{d(\delta w_0)}{dx}]+\\ & M_{xx}\left[\frac{d\delta {\phi }_x}{dx}\right]+Q_x[\delta {\phi }_x+\frac{d(\delta w_0)}{dx}]\}dx\end{array}$

Integrate by Parts

Get rid of derivatives of the variations.

$\begin{array}{cc}\underset{x_a}{\overset{x_b}{\int }}{N}_{xx}& [\frac{d(\delta u_0)}{dx}+\frac{d{w}_0}{dx}\frac{d(\delta w_0)}{dx}]dx={\left[N_{xx}\delta u_0\right]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}\frac{d{N}_{xx}}{dx}\delta u_{0}dx+\\ & {\left[N_{xx}\frac{d{w}_0}{dx}\delta w_0\right]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}\frac{d}{dx}\left(N_{xx}\frac{d{w}_0}{dx}\right)\delta w_{0}dx\\ \\ \underset{x_a}{\overset{x_b}{\int }}{M}_{xx}& \left[\frac{d(\delta {\phi }_x)}{dx}\right]dx={\left[M_{xx}\delta {\phi }_x\right]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}\frac{d{M}_{xx}}{dx}\delta {\phi }_{x}dx\\ \\ \underset{x_a}{\overset{x_b}{\int }}{Q}_x& [\delta {\phi }_x+\frac{d(\delta w_0)}{dx}]dx=\underset{x_a}{\overset{x_b}{\int }}{Q}_x\delta {\phi }_{x}dx+{\left[Q_x\delta w_0\right]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}\frac{d{Q}_x}{dx}\delta w_{0}dx\end{array}$

Collect terms

$\begin{array}{cc}& {\left[N_{xx}\delta u_0\right]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}\frac{d{N}_{xx}}{dx}\delta u_{0}dx+\\ & {[N_{xx}\frac{d{w}_0}{dx}\delta w_0+Q_x\delta w_0]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}[\frac{d}{dx}\left(N_{xx}\frac{d{w}_0}{dx}\right)+\frac{d{Q}_x}{dx}]\delta w_{0}dx+\\ & {\left[M_{xx}\delta {\phi }_x\right]}_{x_a}^{x_b}-\underset{x_a}{\overset{x_b}{\int }}(\frac{d{M}_{xx}}{dx}-Q_x)\delta {\phi }_{x}dx\\ & =\underset{x_a}{\overset{x_b}{\int }}q\delta w_{0}dx+\underset{x_a}{\overset{x_b}{\int }}f\delta u_{0}dx\end{array}$

Euler-Lagrange Equations

$\begin{array}{cc}\underset{x_a}{\overset{x_b}{\int }}(\frac{d{N}_{xx}}{dx}+f)\delta u_{0}dx& ={\left[N_{xx}\delta u_0\right]}_{x_a}^{x_b}\\ \underset{x_a}{\overset{x_b}{\int }}[\frac{d}{dx}\left(N_{xx}\frac{d{w}_0}{dx}\right)+\frac{d{Q}_x}{dx}+q]\delta w_{0}dx& ={[N_{xx}\frac{d{w}_0}{dx}\delta w_0+Q_x\delta w_0]}_{x_a}^{x_b}\\ \underset{x_a}{\overset{x_b}{\int }}(\frac{d{M}_{xx}}{dx}-Q_x)\delta {\phi }_{x}dx& ={\left[M_{xx}\delta {\phi }_x\right]}_{x_a}^{x_b}\end{array}$

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

Constitutive Relations

${\sigma }_{xx}=E{\epsilon }_{xx};{\sigma }_{xz}=G{\gamma }_{xz}$

Then,

$\begin{array}{cc}{N}_{xx}& =A_{xx}{\epsilon }_{xx}^0+B_{xx}{\epsilon }_{xx}^1\\ \\ M_{xx}& =B_{xx}{\epsilon }_{xx}^0+D_{xx}{\epsilon }_{xx}^1\\ \\ Q_x& =S_{xx}{\gamma }_{xz}^0\end{array}$

where

$S_{xx}=K_s\underset{A}{\int }GdA\leftarrow \text{shear stiffness}$

Equilibrium Equations

$\begin{array}{cc}\frac{d}{dx}\left\{A_{xx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]\right\}+f& =0\\ \frac{d}{dx}\{S_{xx}(\frac{d{w}_0}{dx}+{\phi }_x)+A_{xx}\frac{d{w}_0}{dx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]\}+q& =0\\ \frac{d}{dx}\left(D_{xx}\frac{d{\phi }_x}{dx}\right)+S_{xx}(\frac{d{w}_0}{dx}+{\phi }_x)& =0\end{array}$

Weak Form

$\begin{array}{cc}\underset{x_a}{\overset{x_b}{\int }}{A}_{xx}\frac{d(\delta u_0)}{dx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]dx& =\underset{x_a}{\overset{x_b}{\int }}f\delta u_{0}dx+{\left[N_{xx}\delta u_0\right]}_{x_a}^{x_b}\\ \\ \underset{x_a}{\overset{x_b}{\int }}{A}_{xx}\frac{d(\delta w_0)}{dx}\frac{d{w}_0}{dx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]dx& +\underset{x_a}{\overset{x_b}{\int }}{S}_{xx}\frac{d(\delta w_0)}{dx}(\frac{d{w}_0}{dx}+{\phi }_x)dx=\\ & \underset{x_a}{\overset{x_b}{\int }}q\delta u_{0}dx+{\left[(N_{xx}\frac{d{w}_0}{dx}+Q_x)\delta w_0\right]}_{x_a}^{x_b}\\ \\ -\underset{x_a}{\overset{x_b}{\int }}{S}_{xx}\delta {\phi }_x(\frac{d{w}_0}{dx}+{\phi }_x)dx& +\underset{x_a}{\overset{x_b}{\int }}{D}_{xx}\frac{d(\delta {\phi }_x)}{dx}\frac{d{\phi }_x}{dx}dx={\left[M_{xx}\delta {\phi }_x\right]}_{x_a}^{x_b}\end{array}$

Finite element model

Trial Solution

$\begin{array}{cc}{u}_0(x)& =\sum _{j=1}^m{u}_j{\psi }_j^{(1)};\delta u_0={\psi }_i^{(1)}\\ \\ w_0(x)& =\sum _{j=1}^n{w}_j{\psi }_j^{(2)};\delta w_0={\psi }_i^{(2)}\\ \\ {\phi }_x(x)& =\sum _{j=1}^p{s}_j{\psi }_j^{(3)};\delta {\phi }_x={\psi }_i^{(3)}\end{array}$

Element Stiffness Matrix

$\left[\begin{array}{ccccc}{𝐊}^{11}& ⋮& {𝐊}^{12}& ⋮& {𝐊}^{13}\\ & ⋮& & ⋮& \\ {𝐊}^{21}& ⋮& {𝐊}^{22}& ⋮& {𝐊}^{23}\\ & ⋮& & ⋮& \\ {𝐊}^{31}& ⋮& {𝐊}^{32}& ⋮& {𝐊}^{33}\end{array}\right]\left[\begin{array}{c}𝐮\\ \\ 𝐰\\ \\ 𝐬\end{array}\right]=\left[\begin{array}{c}{𝐅}^1\\ \\ {𝐅}^2\\ \\ {𝐅}^3\end{array}\right]$

Choice of Approximate Solutions

Choice 1

${\psi }^{(1)}$ = linear ($m=2$)

${\psi }^{(2)}$ = linear ($n=2$)

${\psi }^{(3)}$ = linear ($p=2$).

Nearly singular stiffness matrix ($6\times 6$).

Choice 2

${\psi }^{(1)}$ = linear ($m=2$)

${\psi }^{(2)}$ = quadratic ($n=3$)

${\psi }^{(3)}$ = linear ($p=2$).

The stiffness matrix is ($7\times 7$). We can statically condense out the interior degree of freedom and get a ($6\times 6$) matrix. The element behaves well.

Choice 3

${\psi }^{(1)}$ = linear ($m=2$)

${\psi }^{(2)}$ = cubic ($n=4$)

${\psi }^{(3)}$ = quadratic ($p=3$)

The stiffness matrix is ($9\times 9$). We can statically condense out the interior degrees of freedom and get a ($6\times 6$) matrix. If the shear and bending stiffnesses are element-wise constant, this element gives exact results.

Shear Locking

Example Case

Linear $u_0$, Linear $w_0$, Linear ${\phi }_x$.

$\frac{d{w}_0}{dx}=\text{constant}.$

But, for thin beams,

$\frac{d{w}_0}{dx}=\text{slope}=-{\phi }_x\leftarrow (\text{linear!})$

If constant ${\phi }_x$

$\frac{d{\phi }_x}{dx}=0$

Also

1. $Q_x=S_{xx}{\phi }_x\ne 0⟹$ Non-zero transverse shear.
2. $M_{xx}=D_{xx}\frac{d{\phi }_x}{dx}=0⟹$ Zero bending energy.

Result: Zero displacements and rotations $⟹$ Shear Locking!

Recall

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

or,

$\frac{d}{dx}\left\{A_{xx}{\epsilon }_{xx}^0\right\}+f=0$

If $f=0$ and $A_{xx}=$ constant

$A_{xx}\frac{d}{dx}({\epsilon }_{xx}^0)=0⟹{\epsilon }_{xx}^0=\text{constant}.$

If there is only bending but no stretching,

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

Hence,

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

Also recall:

$\frac{d}{dx}\{S_{xx}(\frac{d{w}_0}{dx}+{\phi }_x)+A_{xx}\frac{d{w}_0}{dx}[\frac{d{u}_0}{dx}+\frac{1}{2}{\left(\frac{d{w}_0}{dx}\right)}^2]\}+q=0$

or,

$\frac{d}{dx}\{S_{xx}{\gamma }_{xz}^0+A_{xx}\frac{d{w}_0}{dx}{\epsilon }_{xx}^0\}+q=0$

If $q=0$ and $S_{xx}=$ constant, and no membrane strains

$S_{xx}\frac{d}{dx}({\gamma }_{xz}^0)=0⟹{\gamma }_{xz}=\text{constant}=\frac{d{w}_0}{dx}+{\phi }_x$

Hence,

${\phi }_x\approx \frac{d{w}_0}{dx}$

Shape functions need to satisfy:

Template:Center top
$\frac{d{u}_0}{dx}\approx {\left(\frac{d{w}_0}{dx}\right)}^2;\text{and}{\phi }_x\approx \frac{d{w}_0}{dx}$
Template:Center bottom

Example Case 1

Linear $u_0$, Linear $w_0$, Linear ${\phi }_x$.

• First condition $⟹$ constant $=$ constant. Passes! No Membrane Locking.
• Second condition $⟹$ linear $=$ constant. Fails! Shear Locking.

Example Case 2

Linear $u_0$, Quadratic $w_0$, Linear ${\phi }_x$.

• First condition $⟹$ constant $=$ quadratic. Fails! Membrane Locking.
• Second condition $⟹$ linear $=$ linear. Passes! No Shear Locking.

Example Case 3

Quadratic $u_0$, Quadratic $w_0$, Linear ${\phi }_x$.

• First condition $⟹$ linear $=$ quadratic. Fails! Membrane Locking.
• Second condition $⟹$ linear $=$ linear. Passes! No Shear Locking.

Example Case 4

Cubic $u_0$, Quadratic $w_0$, Linear ${\phi }_x$.

• First condition $⟹$ quadratic $=$ quadratic. Passes! No Membrane Locking.
• Second condition $⟹$ linear $=$ linear. Passes! No Shear Locking.

Overcoming Shear Locking

Option 1

• Linear $u_0$, linear $w_0$, linear ${\phi }_x$.
• Equal interpolation for both $w_0$ and ${\varphi }_x$.
• Reduced integration for terms containing ${\varphi }_x$ - treat as constant.

Option 2

• Cubic $u_0$, quadratic $w_0$, linear ${\phi }_x$.
• Stiffness matrix is $9\times 9$.
• Hard to implement.