Nonlinear finite elements/Timoshenko beams

Timoshenko Beam

 Timoshenko beam.

Displacements

{\displaystyle {\begin{aligned}u_{1}&=u_{0}(x)+z\varphi _{x}\\u_{2}&=0\\u_{3}&=w_{0}(x)\end{aligned}}}

Strains

{\displaystyle {\begin{aligned}\varepsilon _{xx}&=\varepsilon _{xx}^{0}+z\varepsilon _{xx}^{1}\\\gamma _{xz}&=\gamma _{xz}^{0}\end{aligned}}}
{\displaystyle {\begin{aligned}\varepsilon _{xx}^{0}&={\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\\\varepsilon _{xx}^{1}&={\cfrac {d\varphi _{x}}{dx}}\\\gamma _{xz}^{0}&=\varphi _{x}+{\cfrac {dw_{0}}{dx}}\end{aligned}}}

Principle of Virtual Work

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

where

{\displaystyle {\begin{aligned}\delta W_{\text{int}}&=\int _{x_{a}}^{x_{b}}\int _{A}(\sigma _{xx}\delta \varepsilon _{xx}+\sigma _{xz}\delta \gamma _{xz})dA~dx\\&=\int _{x_{a}}^{x_{b}}(N_{xx}\delta \varepsilon _{xx}^{0}+M_{xx}\delta \varepsilon _{xx}^{1}+Q_{x}\delta \gamma _{xz}^{0})~dx\\\delta W_{\text{ext}}&=\int _{x_{a}}^{x_{b}}q\delta w_{0}~dx+\int _{x_{a}}^{x_{b}}f\delta u_{0}~dx\end{aligned}}}
${\displaystyle Q_{x}={K_{s}}\int _{A}\sigma _{xz}~dA}$

${\displaystyle K_{s}}$ = shear correction factor

Taking Variations

${\displaystyle \varepsilon _{xx}^{0}={\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}}$

Take variation

${\displaystyle \delta \varepsilon _{xx}^{0}={\cfrac {d(\delta u_{0})}{dx}}+{\frac {1}{2}}\left(2{\cfrac {dw_{0}}{dx}}\right)\left({\cfrac {d(\delta w_{0})}{dx}}\right)={\cfrac {d(\delta u_{0})}{dx}}+{\cfrac {dw_{0}}{dx}}{\cfrac {d(\delta w_{0})}{dx}}~.}$
${\displaystyle \varepsilon _{xx}^{1}={\cfrac {d\varphi _{x}}{dx}}}$

Take variation

${\displaystyle \delta \varepsilon _{xx}^{1}={\cfrac {d\delta \varphi _{x}}{dx}}}$
${\displaystyle \gamma _{xz}^{0}=\varphi _{x}+{\cfrac {dw_{0}}{dx}}}$

Take variation

${\displaystyle \delta \gamma _{xz}^{0}=\delta \varphi _{x}+{\cfrac {d(\delta w_{0})}{dx}}}$

Internal Virtual Work

{\displaystyle {\begin{aligned}\int _{x_{a}}^{x_{b}}N_{xx}\delta \varepsilon _{xx}^{0}~dx&=\int _{x_{a}}^{x_{b}}N_{xx}\left[{\cfrac {d(\delta u_{0})}{dx}}+{\cfrac {dw_{0}}{dx}}{\cfrac {d(\delta w_{0})}{dx}}\right]~dx\\\int _{x_{a}}^{x_{b}}M_{xx}\delta \varepsilon _{xx}^{1}~dx&=\int _{x_{a}}^{x_{b}}M_{xx}\left[{\cfrac {d\delta \varphi _{x}}{dx}}\right]~dx\\\int _{x_{a}}^{x_{b}}Q_{x}\delta \gamma _{xz}^{0}~dx&=\int _{x_{a}}^{x_{b}}Q_{x}\left[\delta \varphi _{x}+{\cfrac {d(\delta w_{0})}{dx}}\right]~dx\end{aligned}}}
{\displaystyle {\begin{aligned}\delta W_{\text{int}}=\int _{x_{a}}^{x_{b}}&\left\{N_{xx}\left[{\cfrac {d(\delta u_{0})}{dx}}+{\cfrac {dw_{0}}{dx}}{\cfrac {d(\delta w_{0})}{dx}}\right]\right.+\\&M_{xx}\left[{\cfrac {d\delta \varphi _{x}}{dx}}\right]+\left.Q_{x}\left[\delta \varphi _{x}+{\cfrac {d(\delta w_{0})}{dx}}\right]\right\}~dx\end{aligned}}}

Integrate by Parts

Get rid of derivatives of the variations.

{\displaystyle {\begin{aligned}\int _{x_{a}}^{x_{b}}N_{xx}&\left[{\cfrac {d(\delta u_{0})}{dx}}+{\cfrac {dw_{0}}{dx}}{\cfrac {d(\delta w_{0})}{dx}}\right]~dx=\left[N_{xx}\delta u_{0}\right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}{\cfrac {dN_{xx}}{dx}}\delta u_{0}~dx+\\&\qquad \qquad \qquad \left[N_{xx}{\cfrac {dw_{0}}{dx}}\delta w_{0}\right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}{\cfrac {d}{dx}}\left(N_{xx}{\cfrac {dw_{0}}{dx}}\right)\delta w_{0}~dx\\\\\int _{x_{a}}^{x_{b}}M_{xx}&\left[{\cfrac {d(\delta \varphi _{x})}{dx}}\right]~dx=\left[M_{xx}\delta \varphi _{x}\right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}{\cfrac {dM_{xx}}{dx}}\delta \varphi _{x}~dx\\\\\int _{x_{a}}^{x_{b}}Q_{x}&\left[\delta \varphi _{x}+{\cfrac {d(\delta w_{0})}{dx}}\right]~dx=\int _{x_{a}}^{x_{b}}Q_{x}\delta \varphi _{x}~dx+\left[Q_{x}\delta w_{0}\right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}{\cfrac {dQ_{x}}{dx}}\delta w_{0}~dx\end{aligned}}}

Collect terms

{\displaystyle {\begin{aligned}&\left[N_{xx}{\delta u_{0}}\right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}{\cfrac {dN_{xx}}{dx}}{\delta u_{0}}~dx+\\&\left[N_{xx}{\cfrac {dw_{0}}{dx}}\delta w_{0}+Q_{x}\delta w_{0}\right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}\left[{\cfrac {d}{dx}}\left(N_{xx}{\cfrac {dw_{0}}{dx}}\right)+{\cfrac {dQ_{x}}{dx}}\right]\delta w_{0}~dx+\\&\left[M_{xx}\delta \varphi _{x}\right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}\left({\cfrac {dM_{xx}}{dx}}-Q_{x}\right)\delta \varphi _{x}~dx\\&=\int _{x_{a}}^{x_{b}}q~\delta w_{0}~dx+\int _{x_{a}}^{x_{b}}f{\delta u_{0}}~dx\end{aligned}}}

Euler-Lagrange Equations

{\displaystyle {\begin{aligned}\int _{x_{a}}^{x_{b}}\left({\cfrac {dN_{xx}}{dx}}+f\right){\delta u_{0}}~dx&=\left[N_{xx}{\delta u_{0}}\right]_{x_{a}}^{x_{b}}\\\int _{x_{a}}^{x_{b}}\left[{\cfrac {d}{dx}}\left(N_{xx}{\cfrac {dw_{0}}{dx}}\right)+{\cfrac {dQ_{x}}{dx}}+q\right]~\delta w_{0}~dx&=\left[N_{xx}{\cfrac {dw_{0}}{dx}}~\delta w_{0}+Q_{x}~\delta w_{0}\right]_{x_{a}}^{x_{b}}\\\int _{x_{a}}^{x_{b}}\left({\cfrac {dM_{xx}}{dx}}-Q_{x}\right)~\delta \varphi _{x}~dx&=\left[M_{xx}~\delta \varphi _{x}\right]_{x_{a}}^{x_{b}}\end{aligned}}}
{\displaystyle {\begin{aligned}{\cfrac {dN_{xx}}{dx}}+f&=0\\{\cfrac {d}{dx}}\left(N_{xx}{\cfrac {dw_{0}}{dx}}\right)+{\cfrac {dQ_{x}}{dx}}+q&=0\\{\cfrac {dM_{xx}}{dx}}-Q_{x}&=0\end{aligned}}}

Constitutive Relations

${\displaystyle \sigma _{xx}=E\varepsilon _{xx}~;\qquad \sigma _{xz}=G\gamma _{xz}}$

Then,

{\displaystyle {\begin{aligned}N_{xx}&=A_{xx}~\varepsilon _{xx}^{0}+B_{xx}~\varepsilon _{xx}^{1}\\\\M_{xx}&=B_{xx}~\varepsilon _{xx}^{0}+D_{xx}~\varepsilon _{xx}^{1}\\\\Q_{x}&=S_{xx}~\gamma _{xz}^{0}\end{aligned}}}

where

${\displaystyle S_{xx}=K_{s}\int _{A}G~dA\qquad \leftarrow \qquad {\text{shear stiffness}}}$

Equilibrium Equations

{\displaystyle {\begin{aligned}{\cfrac {d}{dx}}\left\{A_{xx}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]\right\}+f&=0\\{\cfrac {d}{dx}}\left\{S_{xx}\left({\cfrac {dw_{0}}{dx}}+\varphi _{x}\right)+A_{xx}{\cfrac {dw_{0}}{dx}}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]\right\}+q&=0\\{\cfrac {d}{dx}}\left(D_{xx}{\cfrac {d\varphi _{x}}{dx}}\right)+S_{xx}\left({\cfrac {dw_{0}}{dx}}+\varphi _{x}\right)&=0\end{aligned}}}

Weak Form

{\displaystyle {\begin{aligned}\int _{x_{a}}^{x_{b}}A_{xx}{\cfrac {d(\delta u_{0})}{dx}}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]~dx&=\int _{x_{a}}^{x_{b}}f\delta u_{0}~dx+\left[N_{xx}\delta u_{0}\right]_{x_{a}}^{x_{b}}\\\\\int _{x_{a}}^{x_{b}}A_{xx}{\cfrac {d(\delta w_{0})}{dx}}{\cfrac {dw_{0}}{dx}}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]~dx&+\int _{x_{a}}^{x_{b}}S_{xx}{\cfrac {d(\delta w_{0})}{dx}}\left({\cfrac {dw_{0}}{dx}}+\varphi _{x}\right)~dx=\\&\int _{x_{a}}^{x_{b}}q\delta u_{0}~dx+\left[\left(N_{xx}{\cfrac {dw_{0}}{dx}}+Q_{x}\right)\delta w_{0}\right]_{x_{a}}^{x_{b}}\\\\-\int _{x_{a}}^{x_{b}}S_{xx}\delta \varphi _{x}\left({\cfrac {dw_{0}}{dx}}+\varphi _{x}\right)~dx&+\int _{x_{a}}^{x_{b}}D_{xx}{\cfrac {d(\delta \varphi _{x})}{dx}}{\cfrac {d\varphi _{x}}{dx}}~dx=\left[M_{xx}\delta \varphi _{x}\right]_{x_{a}}^{x_{b}}\end{aligned}}}

Finite element model

Trial Solution

{\displaystyle {\begin{aligned}u_{0}(x)&=\sum _{j=1}^{m}u_{j}~\psi _{j}^{(1)}~;\qquad \qquad \delta u_{0}=\psi _{i}^{(1)}\\\\w_{0}(x)&=\sum _{j=1}^{n}w_{j}~\psi _{j}^{(2)}~;\qquad \qquad \delta w_{0}=\psi _{i}^{(2)}\\\\\varphi _{x}(x)&=\sum _{j=1}^{p}s_{j}~\psi _{j}^{(3)}~;\qquad \qquad \delta \varphi _{x}=\psi _{i}^{(3)}\end{aligned}}}

Element Stiffness Matrix

${\displaystyle {\begin{bmatrix}\mathbf {K} ^{11}&\vdots &\mathbf {K} ^{12}&\vdots &\mathbf {K} ^{13}\\&\vdots &&\vdots &\\\mathbf {K} ^{21}&\vdots &\mathbf {K} ^{22}&\vdots &\mathbf {K} ^{23}\\&\vdots &&\vdots &\\\mathbf {K} ^{31}&\vdots &\mathbf {K} ^{32}&\vdots &\mathbf {K} ^{33}\\\end{bmatrix}}{\begin{bmatrix}\mathbf {u} \\\\\mathbf {w} \\\\\mathbf {s} \end{bmatrix}}={\begin{bmatrix}\mathbf {F} ^{1}\\\\\mathbf {F} ^{2}\\\\\mathbf {F} ^{3}\end{bmatrix}}}$

Choice of Approximate Solutions

Choice 1

${\displaystyle \psi ^{(1)}}$ = linear (${\displaystyle m=2}$)
${\displaystyle \psi ^{(2)}}$ = linear (${\displaystyle n=2}$)
${\displaystyle \psi ^{(3)}}$ = linear (${\displaystyle p=2}$).

Nearly singular stiffness matrix (${\displaystyle 6\times 6}$).

Choice 2

${\displaystyle \psi ^{(1)}}$ = linear (${\displaystyle m=2}$)
${\displaystyle \psi ^{(2)}}$ = quadratic (${\displaystyle n=3}$)
${\displaystyle \psi ^{(3)}}$ = linear (${\displaystyle p=2}$).

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

Choice 3

${\displaystyle \psi ^{(1)}}$ = linear (${\displaystyle m=2}$)
${\displaystyle \psi ^{(2)}}$ = cubic (${\displaystyle n=4}$)
${\displaystyle \psi ^{(3)}}$ = quadratic (${\displaystyle p=3}$)

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

Shear Locking

Example Case

Linear ${\displaystyle u_{0}}$, Linear ${\displaystyle w_{0}}$, Linear ${\displaystyle \varphi _{x}}$.

${\displaystyle {\cfrac {dw_{0}}{dx}}=~~{\text{constant}}.}$

But, for thin beams,

${\displaystyle {\cfrac {dw_{0}}{dx}}=~~{\text{slope}}~~=-\varphi _{x}~~\leftarrow ~~({\text{linear!}})}$

If constant ${\displaystyle \varphi _{x}}$

${\displaystyle {\cfrac {d\varphi _{x}}{dx}}=0}$

Also

1. ${\displaystyle Q_{x}=S_{xx}\varphi _{x}\neq 0\implies }$ Non-zero transverse shear.
2. ${\displaystyle M_{xx}=D_{xx}{\cfrac {d\varphi _{x}}{dx}}=0\implies }$ Zero bending energy.

Result: Zero displacements and rotations ${\displaystyle \implies }$ Shear Locking!

Recall

${\displaystyle {\cfrac {d}{dx}}\left\{A_{xx}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]\right\}+f=0}$

or,

${\displaystyle {\cfrac {d}{dx}}\left\{A_{xx}\varepsilon _{xx}^{0}\right\}+f=0}$

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

${\displaystyle A_{xx}{\cfrac {d}{dx}}(\varepsilon _{xx}^{0})=0\qquad \implies \qquad \varepsilon _{xx}^{0}=~{\text{constant}}.}$

If there is only bending but no stretching,

${\displaystyle \varepsilon _{xx}^{0}=0={\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}}$

Hence,

${\displaystyle {\cfrac {du_{0}}{dx}}\approx \left({\cfrac {dw_{0}}{dx}}\right)^{2}}$

Also recall:

${\displaystyle {\cfrac {d}{dx}}\left\{S_{xx}\left({\cfrac {dw_{0}}{dx}}+\varphi _{x}\right)+A_{xx}{\cfrac {dw_{0}}{dx}}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]\right\}+q=0}$

or,

${\displaystyle {\cfrac {d}{dx}}\left\{S_{xx}\gamma _{xz}^{0}+A_{xx}{\cfrac {dw_{0}}{dx}}\varepsilon _{xx}^{0}\right\}+q=0}$

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

${\displaystyle S_{xx}{\cfrac {d}{dx}}(\gamma _{xz}^{0})=0\qquad \implies \qquad {\gamma _{xz}=~{\text{constant}}}~={\cfrac {dw_{0}}{dx}}+\varphi _{x}}$

Hence,

${\displaystyle \varphi _{x}\approx {\cfrac {dw_{0}}{dx}}}$

Shape functions need to satisfy:

${\displaystyle {\cfrac {du_{0}}{dx}}\approx \left({\cfrac {dw_{0}}{dx}}\right)^{2}~;\qquad {\text{and}}\qquad \varphi _{x}\approx {\cfrac {dw_{0}}{dx}}}$


Example Case 1

Linear ${\displaystyle u_{0}}$, Linear ${\displaystyle w_{0}}$, Linear ${\displaystyle \varphi _{x}}$.

• First condition ${\displaystyle \implies }$ constant ${\displaystyle =}$ constant. Passes! No Membrane Locking.
• Second condition ${\displaystyle \implies }$ linear ${\displaystyle =}$ constant. Fails! Shear Locking.

Example Case 2

Linear ${\displaystyle u_{0}}$, Quadratic ${\displaystyle w_{0}}$, Linear ${\displaystyle \varphi _{x}}$.

• First condition ${\displaystyle \implies }$ constant ${\displaystyle =}$ quadratic. Fails! Membrane Locking.
• Second condition ${\displaystyle \implies }$ linear ${\displaystyle =}$ linear. Passes! No Shear Locking.

Example Case 3

Quadratic ${\displaystyle u_{0}}$, Quadratic ${\displaystyle w_{0}}$, Linear ${\displaystyle \varphi _{x}}$.

• First condition ${\displaystyle \implies }$ linear ${\displaystyle =}$ quadratic. Fails! Membrane Locking.
• Second condition ${\displaystyle \implies }$ linear ${\displaystyle =}$ linear. Passes! No Shear Locking.

Example Case 4

Cubic ${\displaystyle u_{0}}$, Quadratic ${\displaystyle w_{0}}$, Linear ${\displaystyle \varphi _{x}}$.

• First condition ${\displaystyle \implies }$ quadratic ${\displaystyle =}$ quadratic. Passes! No Membrane Locking.
• Second condition ${\displaystyle \implies }$ linear ${\displaystyle =}$ linear. Passes! No Shear Locking.

Overcoming Shear Locking

Option 1

• Linear ${\displaystyle u_{0}}$, linear ${\displaystyle w_{0}}$, linear ${\displaystyle \varphi _{x}}$.
• Equal interpolation for both ${\displaystyle w_{0}}$ and ${\displaystyle \phi _{x}}$.
• Reduced integration for terms containing ${\displaystyle \phi _{x}}$ - treat as constant.

Option 2

• Cubic ${\displaystyle u_{0}}$, quadratic ${\displaystyle w_{0}}$, linear ${\displaystyle \varphi _{x}}$.
• Stiffness matrix is ${\displaystyle 9\times 9}$.
• Hard to implement.