Nonlinear finite elements/Timoshenko beams

Displacements


\begin{align} u_1 & = u_0(x) + z \varphi_x \\ u_2 & = 0 \\ u_3 & = w_0(x) \end{align} $$

Strains


\begin{align} \varepsilon_{xx}& = \varepsilon_{xx}^0 + z \varepsilon_{xx}^1 \\ \gamma_{xz}& = \gamma_{xz}^0 \end{align} $$



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

Principle of Virtual Work


\delta W_{\text{int}} = \delta W_{\text{ext}} $$ where
 * $$\begin{align}

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

Q_x = { K_s} \int_A \sigma_{xz}~dA $$ $$K_s$$ = shear correction factor

Taking Variations


\varepsilon_{xx}^0= \cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2 $$ Take variation

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



\varepsilon_{xx}^1= \cfrac{d\varphi_x}{dx} $$ Take variation

\delta \varepsilon_{xx}^1= \cfrac{d\delta \varphi_x}{dx} $$



\gamma_{xz}^0= \varphi_x + \cfrac{dw_0}{dx} $$ Take variation

\delta \gamma_{xz}^0= \delta \varphi_x + \cfrac{d(\delta w_0)}{dx} $$

Internal Virtual Work

 * $$\begin{align}

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

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

Integrate by Parts
Get rid of derivatives of the variations.
 * $$\begin{align}

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

Collect terms

 * $$\begin{align}

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

Euler-Lagrange Equations

 * $$\begin{align}

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



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

Constitutive Relations


\sigma_{xx} = E\varepsilon_{xx} ~; \qquad \sigma_{xz} = G\gamma_{xz} $$ Then,

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

S_{xx} = K_s\int_A G~dA \qquad \leftarrow \qquad \text{shear stiffness} $$

Equilibrium Equations

 * $$\begin{align}

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

Weak Form

 * $$\begin{align}

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

Trial Solution

 * $$\begin{align}

u_0(x) & = \sum_{j=1}^m u_j~\psi^{(1)}_j~; \qquad \qquad \delta u_0 = \psi^{(1)}_i \\ \\ w_0(x) & = \sum_{j=1}^n w_j~\psi^{(2)}_j~; \qquad \qquad \delta w_0 = \psi^{(2)}_i \\ \\ \varphi_x(x) & = \sum_{j=1}^p s_j~\psi^{(3)}_j~; \qquad \qquad \delta \varphi_x = \psi^{(3)}_i \end{align}$$

Element Stiffness Matrix


\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 1

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

Nearly singular stiffness matrix ($$6 \times 6$$).

Choice 2

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

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

Choice 3

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

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

Example Case
Linear $$u_0$$, Linear $$w_0$$, Linear $$\varphi_x$$.



\cfrac{dw_0}{dx} = \text{constant}. $$ But, for thin beams,

\cfrac{dw_0}{dx} = \text{slope} = -\varphi_x \leftarrow (\text{linear!}) $$ If constant $$\varphi_x$$

\cfrac{d\varphi_x}{dx} = 0 $$

Also
 * 1) $$Q_x = S_{xx}\varphi_x \ne 0 \implies$$ Non-zero transverse shear.
 * 2) $$M_{xx} = D_{xx} \cfrac{d\varphi_x}{dx} = 0 \implies $$ Zero bending energy.

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

Recall

\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,

\cfrac{d}{dx}\left\{ A_{xx} \varepsilon_{xx}^0 \right\} + f= 0 $$ If $$f = 0$$ and $$A_{xx} = $$ constant

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,

\varepsilon_{xx}^0 = 0 = \cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2 $$ Hence,

\cfrac{du_0}{dx} \approx \left(\cfrac{dw_0}{dx}\right)^2 $$

Also recall:

\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,

\cfrac{d}{dx}\left\{ S_{xx} \gamma_{xz}^0 + A_{xx} \cfrac{dw_0}{dx} \varepsilon_{xx}^0 \right\} + q = 0 $$ If $$q = 0$$ and $$S_{xx} = $$ constant, and no membrane strains

S_{xx}\cfrac{d}{dx} (\gamma_{xz}^0)= 0 \qquad \implies \qquad {\gamma_{xz} = ~\text{constant}}~ = \cfrac{dw_0}{dx} + \varphi_x $$ Hence,

\varphi_x \approx \cfrac{dw_0}{dx} $$

Shape functions need to satisfy: $$ \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 $$u_0$$, Linear $$w_0$$, Linear $$\varphi_x$$.


 * First condition $$\implies$$ constant $$=$$ constant. Passes! No Membrane Locking.
 * Second condition $$\implies$$ linear $$=$$ constant. Fails! Shear Locking.

Example Case 2
Linear $$u_0$$, Quadratic $$w_0$$, Linear $$\varphi_x$$.


 * First condition $$\implies$$ constant $$=$$ quadratic. Fails!  Membrane Locking.
 * Second condition $$\implies$$ linear $$=$$ linear. Passes! No Shear Locking.

Example Case 3
Quadratic $$u_0$$, Quadratic $$w_0$$, Linear $$\varphi_x$$.


 * First condition $$\implies$$ linear $$=$$ quadratic. Fails! Membrane Locking.
 * Second condition $$\implies$$ linear $$=$$ linear. Passes! No Shear Locking.

Example Case 4
Cubic $$u_0$$, Quadratic $$w_0$$, Linear $$\varphi_x$$.


 * First condition $$\implies$$ quadratic $$=$$ quadratic. Passes!  No Membrane Locking.
 * Second condition $$\implies$$ linear $$=$$ linear. Passes!  No Shear Locking.

Option 1

 * Linear $$u_0$$, linear $$w_0$$, linear $$\varphi_x$$.
 * Equal interpolation for both $$w_0$$ and $$\phi_x$$.
 * Reduced integration for terms containing $$\phi_x$$ - treat as constant.

Option 2

 * Cubic $$u_0$$, quadratic $$w_0$$, linear $$\varphi_x$$.
 * Stiffness matrix is $$9\times9$$.
 * Hard to implement.