Nonlinear finite elements/Euler Bernoulli beams

Displacements


\begin{align} u_1 & = u_0(x) - z \cfrac{dw_0}{dx} \\ u_2 & = 0 \\ u_3 & = w_0(x) \end{align} $$

Strains


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



\begin{align} \varepsilon_{xx}^0 & = \cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2 \\ \varepsilon_{xx}^1 & = -\cfrac{d^2w_0}{dx^2} \end{align} $$

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

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{align}

\frac{\partial u_1}{\partial x_1} & = \cfrac{du_0}{dx_1} - x_3\cfrac{d^2w_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{align}$$

von Karman strains
The von Karman strains

\begin{align} \varepsilon_{11} & = \cfrac{du_0}{dx_1} - x_3\cfrac{d^2w_0}{dx_1^2} + \frac{1}{2}\left[ \left(\cfrac{du_0}{dx_1}-x_3\cfrac{d^2w_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^2w_0}{dx_1^2}\right) \left(\cfrac{dw_0}{dx_1}\right)\right] \\ \varepsilon_{12} & = 0 \end{align} $$

Balance of forces


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

Stress Resultants

 * $$\begin{align}

N_{xx} & = \int_A \sigma_{xx}~ dA \\ M_{xx} & = \int_A z\sigma_{xx}~ dA \end{align}$$

Stress-Strain equation


\sigma_{xx} = E \varepsilon_{xx} $$

Stress Resultant - Displacement relations


\begin{align} 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^2w_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^2w_0}{dx^2} \end{align} $$

Extensional/Bending Stiffness


\begin{align} 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^2E~dA\qquad \leftarrow \qquad \text{bending stiffness} \end{align} $$

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

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

Axial Equation


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

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

Bending Equation


\begin{align} \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^2w_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{align} $$ where
 * $$\begin{align}

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

Finite Element Model

 * $$\begin{align}

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{align}$$ where $$\theta = -(dw_0/dx)$$.

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{align}

\mathbf{u} & = [u_1 \quad u_2]^T \\ \mathbf{d} & = [w_1 \quad \theta_1 \quad w_2 \quad \theta_2]^T \end{align}$$


 * $$\begin{align}

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

Symmetric Stiffness Matrix


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

Load Vector


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

Newton-Raphson Solution


\mathbf{K}(\mathbf{U}) \mathbf{U} = \mathbf{F} $$ where
 * $$\begin{align}

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^1_2, ~ F_3= F^2_1, ~ F_4= F^2_2, ~ F_5= F^2_3, ~ F_6= F^2_4 \end{align}$$ 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{align} i=1\dots2;~j=1\dots2 &: \\ & { T^{11}_{ij} = K^{11}_{ij}} \\ \\ i=1\dots2;~j=1\dots4 &: \\ & { T^{12}_{ij} = 2 K^{12}_{ij}} \\ \\ i=1\dots4;~j=1\dots2 &: \\ & { T^{21}_{ij} = 2 K^{21}_{ij}} \\ \\ i=1\dots4;~j=1\dots4 &: \\ & { T^{22}_{ij} = K^{22}_{ij} + \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{align} $$

Load Steps
Recall

N_{xx} = A_{xx}\left[\cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2 \right] - \cfrac{d^2w_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 $$


 * 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{align} 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{align} $$

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

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