Nonlinear finite elements/Euler Bernoulli beams

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

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

${\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^{2}w_{0}}{dx^{2}}}\end{aligned}}}$

The von Karman strains

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

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

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

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

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

${\displaystyle {\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 ${\displaystyle E}$ is constant, and ${\displaystyle x}$-axis passes through centroid

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

${\displaystyle {\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

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

${\displaystyle {\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

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

${\displaystyle {\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

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

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

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

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

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

where

${\displaystyle {\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

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

For Newton iterations, we use the algorithm

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

where the tangent stiffness matrix is given by

${\displaystyle \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~.}$

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

Recall

${\displaystyle 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.

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

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

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

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

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

• Continue until F is reached.

Membrane strain:

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

or

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

Hence, shape functions should be such that

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

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

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

Assume ${\displaystyle A_{xx}}$ = constant.

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

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