Nonlinear finite elements/Euler Bernoulli beams

Euler-Bernoulli Beam

Euler-Bernoulli beam

Displacements

${\begin{aligned}u_{1}&=u_{0}(x)-z{\cfrac {dw_{0}}{dx}}\\u_{2}&=0\\u_{3}&=w_{0}(x)\end{aligned}}$

Strains

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

${\begin{aligned}\varepsilon _{xx}^{0}&={\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\\\varepsilon _{xx}^{1}&=-{\cfrac {d^{2}w_{0}}{dx^{2}}}\end{aligned}}$

Strain-Displacement Relations

$\varepsilon _{ij}={\frac {1}{2}}\left({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}}\right)+{\frac {1}{2}}\left({\frac {\partial u_{m}}{\partial x_{i}}}{\frac {\partial u_{m}}{\partial x_{j}}}\right)$

${\begin{aligned}\varepsilon _{11}&={\frac {1}{2}}\left({\frac {\partial u_{1}}{\partial x_{1}}}+{\frac {\partial u_{1}}{\partial x_{1}}}\right)+{\frac {1}{2}}\left({\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}}}\right)\\\varepsilon _{22}&={\frac {1}{2}}\left({\frac {\partial u_{2}}{\partial x_{2}}}+{\frac {\partial u_{2}}{\partial x_{2}}}\right)+{\frac {1}{2}}\left({\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}}}\right)\\\varepsilon _{33}&={\frac {1}{2}}\left({\frac {\partial u_{3}}{\partial x_{3}}}+{\frac {\partial u_{3}}{\partial x_{3}}}\right)+{\frac {1}{2}}\left({\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}}}\right)\\\varepsilon _{23}&={\frac {1}{2}}\left({\frac {\partial u_{2}}{\partial x_{3}}}+{\frac {\partial u_{3}}{\partial x_{2}}}\right)+{\frac {1}{2}}\left({\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}}}\right)\\\varepsilon _{31}&={\frac {1}{2}}\left({\frac {\partial u_{3}}{\partial x_{1}}}+{\frac {\partial u_{1}}{\partial x_{3}}}\right)+{\frac {1}{2}}\left({\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}}}\right)\\\varepsilon _{12}&={\frac {1}{2}}\left({\frac {\partial u_{1}}{\partial x_{2}}}+{\frac {\partial u_{2}}{\partial x_{1}}}\right)+{\frac {1}{2}}\left({\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}}}\right)\end{aligned}}$

The displacements

$u_{1}=u_{0}(x_{1})-x_{3}{\cfrac {dw_{0}}{dx_{1}}}~;~~u_{2}=0~;~~u_{3}=w_{0}(x_{1})$

The derivatives

${\begin{aligned}{\frac {\partial u_{1}}{\partial x_{1}}}&={\cfrac {du_{0}}{dx_{1}}}-x_{3}{\cfrac {d^{2}w_{0}}{dx_{1}^{2}}}~;~~&{\frac {\partial u_{1}}{\partial x_{2}}}&=0~;~~&{\frac {\partial u_{1}}{\partial x_{3}}}&=-{\cfrac {dw_{0}}{dx_{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}}}&={\cfrac {dw_{0}}{dx_{1}}}~;~~&{\frac {\partial u_{3}}{\partial x_{2}}}&=0~;~~&{\frac {\partial u_{3}}{\partial x_{3}}}&=0\end{aligned}}$

von Karman strains

The von Karman strains

${\begin{aligned}\varepsilon _{11}&={\cfrac {du_{0}}{dx_{1}}}-x_{3}{\cfrac {d^{2}w_{0}}{dx_{1}^{2}}}+{\frac {1}{2}}\left[\left({\cfrac {du_{0}}{dx_{1}}}-x_{3}{\cfrac {d^{2}w_{0}}{dx_{1}^{2}}}\right)^{2}+\left({\cfrac {dw_{0}}{dx_{1}}}\right)^{2}\right]\\\varepsilon _{22}&=0\\\varepsilon _{33}&={\frac {1}{2}}\left({\cfrac {dw_{0}}{dx_{1}}}\right)^{2}\\\varepsilon _{23}&=0\\\varepsilon _{31}&={\frac {1}{2}}\left({\cfrac {dw_{0}}{dx_{1}}}-{\cfrac {dw_{0}}{dx_{1}}}\right)-{\frac {1}{2}}\left[\left({\cfrac {du_{0}}{dx_{1}}}-x_{3}{\cfrac {d^{2}w_{0}}{dx_{1}^{2}}}\right)\left({\cfrac {dw_{0}}{dx_{1}}}\right)\right]\\\varepsilon _{12}&=0\end{aligned}}$

Equilibrium Equations

Balance of forces

${\begin{aligned}{\cfrac {dN_{xx}}{dx}}+f(x)&=0\\{\cfrac {d^{2}M_{xx}}{dx^{2}}}+q(x)+{\cfrac {d}{dx}}\left(N_{xx}{\cfrac {dw_{0}}{dx}}\right)&=0\end{aligned}}$

Stress Resultants

${\begin{aligned}N_{xx}&=\int _{A}\sigma _{xx}~dA\\M_{xx}&=\int _{A}z\sigma _{xx}~dA\end{aligned}}$

Constitutive Relations

Stress-Strain equation

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

Stress Resultant - Displacement relations

${\begin{aligned}N_{xx}&=A_{xx}\varepsilon _{xx}^{0}+B_{xx}\varepsilon _{xx}^{1}=A_{xx}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]-B_{xx}{\cfrac {d^{2}w_{0}}{dx^{2}}}\\M_{xx}&=B_{xx}\varepsilon _{xx}^{0}+D_{xx}\varepsilon _{xx}^{1}=B_{xx}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]-D_{xx}{\cfrac {d^{2}w_{0}}{dx^{2}}}\end{aligned}}$

Extensional/Bending Stiffness

${\begin{aligned}A_{xx}&=\int _{A}E~dA\qquad \leftarrow \qquad {\text{extensional stiffness}}\\B_{xx}&=\int _{A}zE~dA\qquad \leftarrow \qquad {\text{extensional-bending stiffness}}\\D_{xx}&=\int _{A}z^{2}E~dA\qquad \leftarrow \qquad {\text{bending stiffness}}\end{aligned}}$

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

${\begin{aligned}A_{xx}&=E\int _{A}~dA=EA\\B_{xx}&=E\int _{A}z~dA=0\\D_{xx}&=E\int _{A}z^{2}~dA=EI\end{aligned}}$

Weak Forms

Axial Equation

${\begin{aligned}\int _{x_{a}}^{x_{b}}{\cfrac {d(\delta u_{0})}{dx}}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]A_{xx}~dx&=\int _{x_{a}}^{x_{b}}(\delta u_{0})f~dx+\\&\delta u_{0}(x_{a})Q_{1}+\delta u_{0}(x_{b})Q_{4}\end{aligned}}$

where

${\begin{aligned}\delta u_{0}&:=v_{1}\\Q_{1}&:=-N_{xx}(x_{a})\\Q_{4}&:=N_{xx}(x_{b})\end{aligned}}$

Bending Equation

${\begin{aligned}\int _{x_{a}}^{x_{b}}\left\{{\cfrac {d(\delta w_{0})}{dx}}\right.&\left[{\cfrac {du_{0}}{dx}}+{\cfrac {1}{2}}~\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]{\cfrac {dw_{0}}{dx}}A_{xx}+\left.{\cfrac {d^{2}(\delta w_{0})}{dx^{2}}}\left({\cfrac {d^{2}w_{0}}{dx^{2}}}\right)D_{xx}\right\}~dx=\\&\int _{x_{a}}^{x_{b}}(\delta w_{0})q~dx+\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{aligned}}$

where

${\begin{aligned}\delta w_{0}&:=v_{2}&\delta \theta &:={\cfrac {dv_{2}}{dx}}\\Q_{2}&:=-\left[{\cfrac {dM_{xx}}{dx}}+N_{xx}{\cfrac {dw_{0}}{dx}}\right]_{x_{a}}&Q_{5}&:=\left[{\cfrac {dM_{xx}}{dx}}+N_{xx}{\cfrac {dw_{0}}{dx}}\right]_{x_{b}}\\Q_{3}&:=-M_{xx}(x_{a})&Q_{6}&:=M_{xx}(x_{b})\end{aligned}}$

Finite Element Model

Finite element model for Euler Bernoulli beam

${\begin{aligned}u_{0}(x)&=u_{1}\psi _{1}(x)+u_{2}\psi _{2}(x)\\w_{0}(x)&=w_{1}\phi _{1}(x)+\theta _{1}\phi _{2}(x)+w_{2}\phi _{3}(x)+\theta _{2}\phi _{4}(x)\end{aligned}}$

where $\theta =-(dw_{0}/dx)$.

Hermite Cubic Shape Functions

Hermite shape functions for beam

Finite Element Equations

${\begin{bmatrix}\mathbf {K} ^{11}&\vdots &\mathbf {K} ^{12}\\&\vdots &\\\mathbf {K} ^{21}&\vdots &\mathbf {K} ^{22}\end{bmatrix}}{\begin{bmatrix}\mathbf {u} \\\\\mathbf {d} \end{bmatrix}}={\begin{bmatrix}\mathbf {F} ^{1}\\\\\mathbf {F} ^{2}\end{bmatrix}}$

where

${\begin{aligned}\mathbf {u} &=[u_{1}\quad u_{2}]^{T}\\\mathbf {d} &=[w_{1}\quad \theta _{1}\quad w_{2}\quad \theta _{2}]^{T}\end{aligned}}$

${\begin{aligned}\mathbf {K} ^{11}&=2\times 2;\qquad &\mathbf {K} ^{12}&=2\times 4\\\mathbf {K} ^{21}&=4\times 2;\qquad &\mathbf {K} ^{22}&=4\times 4\end{aligned}}$

Symmetric Stiffness Matrix

${\begin{aligned}K_{ij}^{11}&=\int _{x_{a}}^{x_{b}}A_{xx}{\cfrac {d\psi _{i}}{dx}}{\cfrac {d\psi _{j}}{dx}}~dx\\K_{ij}^{12}&={\frac {1}{2}}\int _{x_{a}}^{x_{b}}\left(A_{xx}{\cfrac {dw_{0}}{dx}}\right){\cfrac {d\psi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}~dx\\K_{ij}^{21}&={\frac {1}{2}}\int _{x_{a}}^{x_{b}}\left(A_{xx}{\cfrac {dw_{0}}{dx}}\right){\cfrac {d\phi _{i}}{dx}}{\cfrac {d\psi _{j}}{dx}}~dx\\K_{ij}^{22}&=\int _{x_{a}}^{x_{b}}\left\{{\frac {1}{2}}A_{xx}\left[{\cfrac {du_{0}}{dx}}+\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]{\cfrac {d\phi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}+D_{xx}{\cfrac {d^{2}\phi _{i}}{dx^{2}}}{\cfrac {d^{2}\phi _{j}}{dx^{2}}}\right\}~dx\end{aligned}}$

Load Vector

${\begin{aligned}F_{i}^{1}&=\int _{x_{a}}^{x_{b}}\psi _{i}f~dx+\psi _{i}(x_{a})Q_{1}+\psi _{i}(x_{b})Q_{4}\\F_{i}^{2}&=\int _{x_{a}}^{x_{b}}\phi _{i}q~dx+\phi _{i}(x_{a})Q_{2}+\phi _{i}(x_{b})Q_{5}+{\cfrac {d\phi _{i}}{dx}}(x_{a})Q_{3}+{\cfrac {d\phi _{i}}{dx}}(x_{b})Q_{6}\end{aligned}}$

Newton-Raphson Solution

$\mathbf {K} (\mathbf {U} )\mathbf {U} =\mathbf {F}$

where

${\begin{aligned}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{aligned}}$

The residual is

$\mathbf {R} =\mathbf {K} \mathbf {U} -\mathbf {F} ~.$

For Newton iterations, we use the algorithm

$\mathbf {U} ^{r+1}=\mathbf {U} ^{r}-(\mathbf {T} ^{r})^{-1}\mathbf {R} ^{r}$

where the tangent stiffness matrix is given by

$\mathbf {T} ^{r}={\frac {\partial \mathbf {R} ^{r}}{\partial \mathbf {U} }};\quad {\text{or}}\quad T_{ij}={\frac {\partial R_{i}}{\partial U_{j}}},\qquad i=1\dots 6,j=1\dots 6~.$

Tangent Stiffness Matrix

${\begin{aligned}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}}\int _{x_{a}}^{x_{b}}A_{xx}\left[{\cfrac {du_{0}}{dx}}+2\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]{\cfrac {d\phi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}~dx}\end{aligned}}$

Load Steps

Recall

$N_{xx}=A_{xx}\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]-{B_{xx}}~~{\cfrac {d^{2}w_{0}}{dx^{2}}}$

• Divide load into small increments.

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

• Compute $\mathbf {u}$ and $\mathbf {d}$ for first load step,

$\mathbf {K} (\mathbf {U} _{0})\mathbf {U} _{1}=\Delta F_{1}$

Stiffness of Euler-Bernoulli beam.

• Compute $\mathbf {u}$ and $\mathbf {d}$ for second load step,

$\mathbf {K} (\mathbf {U} _{1})\mathbf {U} _{2}=\Delta F_{1}+\Delta F_{2}$

• Continue until F is reached.

Membrane Locking

Recall

${\begin{bmatrix}\mathbf {K} ^{11}&\vdots &\mathbf {K} ^{12}\\&\vdots &\\\mathbf {K} ^{21}&\vdots &\mathbf {K} ^{22}\end{bmatrix}}{\begin{bmatrix}\mathbf {u} \\\\\mathbf {d} \end{bmatrix}}={\begin{bmatrix}\mathbf {F} ^{1}\\\\\mathbf {F} ^{2}\end{bmatrix}}$

where

${\begin{aligned}K_{ij}^{11}&=\int _{x_{a}}^{x_{b}}A_{xx}{\cfrac {d\psi _{i}}{dx}}{\cfrac {d\psi _{j}}{dx}}~dx\\K_{ij}^{12}&={\frac {1}{2}}\int _{x_{a}}^{x_{b}}\left(A_{xx}{\cfrac {dw_{0}}{dx}}\right){\cfrac {d\psi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}~dx\\K_{ij}^{21}&={\frac {1}{2}}\int _{x_{a}}^{x_{b}}\left(A_{xx}{\cfrac {dw_{0}}{dx}}\right){\cfrac {d\phi _{i}}{dx}}{\cfrac {d\psi _{j}}{dx}}~dx\\K_{ij}^{22}&=\int _{x_{a}}^{x_{b}}\left\{{\frac {1}{2}}A_{xx}\left[{\cfrac {du_{0}}{dx}}+\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]{\cfrac {d\phi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}+D_{xx}{\cfrac {d^{2}\phi _{i}}{dx^{2}}}{\cfrac {d^{2}\phi _{j}}{dx^{2}}}\right\}~dx\end{aligned}}$

Mebrane locking in Euler-Bernoulli beam

For Hinged-Hinged

Membrane strain:

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

or

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

Hence, shape functions should be such that

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

$u_{0}$ linear, $w_{0}$ cubic $\implies$ Element Locks! Too stiff.

Selective Reduced Integration

• Assume $u_{0}$ is linear ;~~ $w_{0}$ is cubic.
• Then ${\cfrac {du_{0}}{dx}}\equiv {\cfrac {d\psi _{i}}{dx}}$ is constant, and ${\cfrac {dw_{0}}{dx}}\equiv {\cfrac {d\phi _{i}}{dx}}$ is quadratic.
• Try to keep $\varepsilon _{xx}^{0}=$ constant.

${\begin{aligned}K_{ij}^{11}&=\int _{x_{a}}^{x_{b}}A_{xx}{\cfrac {d\psi _{i}}{dx}}{\cfrac {d\psi _{j}}{dx}}~dx\\K_{ij}^{12}&={\frac {1}{2}}\int _{x_{a}}^{x_{b}}\left(A_{xx}{\cfrac {dw_{0}}{dx}}\right){\cfrac {d\psi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}~dx\\K_{ij}^{22}&=\int _{x_{a}}^{x_{b}}\left\{{\frac {1}{2}}A_{xx}\left[{\cfrac {du_{0}}{dx}}+\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]{\cfrac {d\phi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}+D_{xx}{\cfrac {d^{2}\phi _{i}}{dx^{2}}}{\cfrac {d^{2}\phi _{j}}{dx^{2}}}\right\}~dx\end{aligned}}$

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

Full integration

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

Assume $A_{xx}$ = constant.

${\begin{aligned}K_{ij}^{11}&=A_{xx}\int _{x_{a}}^{x_{b}}{\cfrac {d\psi _{i}}{dx}}{\cfrac {d\psi _{j}}{dx}}~dx\\&=A_{xx}\int _{-1}^{1}\left(J^{-1}{\cfrac {d\psi _{i}(\xi )}{d\xi }}\right)\left(J^{-1}{\cfrac {d\psi _{j}(\xi )}{d\xi }}\right)J~d\xi =A_{xx}\int _{-1}^{1}F(\xi )~d\xi \\&\approx A_{xx}W_{1}F(\xi _{1})\leftarrow {\text{one-point integration}}\end{aligned}}$

${\begin{aligned}K_{ij}^{12}&={\frac {1}{2}}\int _{x_{a}}^{x_{b}}\left(A_{xx}{\cfrac {dw_{0}}{dx}}\right){\cfrac {d\psi _{i}}{dx}}{\cfrac {d\phi _{j}}{dx}}~dx\\&={\cfrac {A_{xx}}{2}}\int _{-1}^{1}\left(\sum _{i=1}^{4}w_{i}J^{-1}{\cfrac {d\phi _{i}(\xi )}{d\xi }}\right)\left(J^{-1}{\cfrac {d\psi _{i}(\xi )}{d\xi }}\right)\left(J^{-1}{\cfrac {d\phi _{j}(\xi )}{d\xi }}\right)~dx\\&\approx A_{xx}\left[W_{1}F(\xi _{1})+W_{2}F(\xi _{2})+W_{3}F(\xi _{3})\right]\leftarrow {\text{full integration}}\\&\approx A_{xx}\left[W_{1}F(\xi _{1})+W_{2}F(\xi _{2})\right]\leftarrow {\text{reduced integration}}\end{aligned}}$

← Natural vibration

Nonlinear finite elements

Timoshenko beams →