Nonlinear finite elements/Euler Bernoulli beams

Euler-Bernoulli Beam

Euler-Bernoulli beam

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}$

Equilibrium Equations

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}$

Constitutive Relations

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}$

Weak Forms

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

Finite element model for Euler Bernoulli beam

$\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)$.

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{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] - {{B_{xx}}}~~\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$

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

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