Nonlinear finite elements/Homework 5/Solutions/Problem 1

Problem 1: Nonlinear Beam Bending
 Given:

The differential equations governing the bending of straight beams are

\begin{align} \cfrac{dN_{xx}}{dx} + f(x) & = 0 \\ \cfrac{dV}{dx} + q(x) & = 0 \\ \cfrac{dM_{xx}}{dx} -V + N_{xx}\cfrac{dw_0}{dx} & = 0 \end{align} $$

Part 1
Show that the weak forms of these equations can be written as
 * $$\begin{align}

\int_{x_a}^{x_b} \cfrac{dv_1}{dx}N_{xx}~dx & = \int_{x_a}^{x_b} v_1 f~dx + \left.v_1 N_{xx}\right|_{x_a}^{x_b}\\ \int_{x_a}^{x_b} \left[ \cfrac{dv_2}{dx}\left(\cfrac{dw_0}{dx}N_{xx}\right) - \cfrac{d^2v_2}{dx^2}M_{xx}\right]~dx & = \int_{x_a}^{x_b} v_2 q~dx + \left.v_2\left(\cfrac{dw_0}{dx}N_{xx}+\cfrac{dM_{xx}}{dx}\right) \right|_{x_a}^{x_b} - \left.\cfrac{dv_2}{dx}M_{xx} \right|_{x_a}^{x_b} \end{align}$$

First we get rid of the shear force term by combining the second and third equations to get
 * $$\begin{align}

\cfrac{dN_{xx}}{dx} + f(x) & = 0 \text{(1)} \qquad \\ \cfrac{d^2M_{xx}}{dx^2} + q(x) + \cfrac{d}{dx}\left(N_{xx}\cfrac{dw_0}{dx}\right) & = 0 \text{(2)} \qquad \end{align}$$ Let $$v_1$$ be the weighting function for equation (1) and let $$v_2$$ be the weighting function for equation (2).

Then the weak forms of the two equations are
 * $$\begin{align}

\int_{x_a}^{x_b} v_1\left(\cfrac{dN_{xx}}{dx} + f(x)\right)~dx & = 0 \text{(3)} \qquad \\ \int_{x_a}^{x_b} v_2\left[\cfrac{d^2M_{xx}}{dx^2} + q(x) + \cfrac{d}{dx}\left(N_{xx}\cfrac{dw_0}{dx}\right)\right]~dx & = 0 \text{(4)} \qquad \end{align}$$ To get the symmetric weak forms, we integrate by parts (even though the symmetry is not obvious at this stage) to get
 * $$\begin{align}

\left.(v_1~N_{xx})\right|_{x_a}^{x_b} & - \int_{x_a}^{x_b} \cfrac{dv_1}{dx} N_{xx}~dx + \int_{x_a}^{x_b} v_1~f(x)~dx = 0 \text{(5)} \qquad \\ \left.\left(v_2~\cfrac{dM_{xx}}{dx}\right)\right|_{x_a}^{x_b} & - \int_{x_a}^{x_b} \cfrac{dv_2}{dx}\cfrac{dM_{xx}}{dx}~dx + \int_{x_a}^{x_b} v_2~q(x)~dx + \\ & \left.\left(v_2~N_{xx}\cfrac{dw_0}{dx}\right)\right|_{x_a}^{x_b} - \int_{x_a}^{x_b} \cfrac{dv_2}{dx}\left(N_{xx}\cfrac{dw_0}{dx}\right)~dx = 0 \text{(6)} \qquad \end{align}$$ Integrate equation (6) again by parts, and get
 * $$\begin{align}

\left.\left(v_2~\cfrac{dM_{xx}}{dx}\right)\right|_{x_a}^{x_b} & - \left.\left(\cfrac{dv_2}{dx}~M_{xx}\right)\right|_{x_a}^{x_b} + \int_{x_a}^{x_b} \cfrac{d^2v_2}{dx^2}M_{xx}~dx + \int_{x_a}^{x_b} v_2~q(x)~dx + \\ & \left.\left(v_2~N_{xx}\cfrac{dw_0}{dx}\right)\right|_{x_a}^{x_b} - \int_{x_a}^{x_b} \cfrac{dv_2}{dx}\left(N_{xx}\cfrac{dw_0}{dx}\right)~dx = 0 \text{(7)} \qquad \end{align}$$ Collect terms and rearrange equations (5) and (7) to get
 * $$\begin{align}

\int_{x_a}^{x_b} \cfrac{dv_1}{dx} N_{xx}~dx & = \int_{x_a}^{x_b} v_1~f(x)~dx + \left.(v_1~N_{xx})\right|_{x_a}^{x_b} \text{(8)} \qquad \\ \int_{x_a}^{x_b} \left[ \cfrac{dv_2}{dx}\left(N_{xx}\cfrac{dw_0}{dx}\right) - \cfrac{d^2v_2}{dx^2}M_{xx}\right]~dx & = \int_{x_a}^{x_b} v_2~q(x)~dx +\\ &\left[\left(v_2~\cfrac{dM_{xx}}{dx}\right) - \left(\cfrac{dv_2}{dx}~M_{xx}\right) + \left(v_2~N_{xx}\cfrac{dw_0}{dx}\right)\right]_{x_a}^{x_b} \text{(9)} \qquad \end{align}$$ Rewriting the equations, we get

\begin{align} \int_{x_a}^{x_b} \cfrac{dv_1}{dx} N_{xx}~dx & = \int_{x_a}^{x_b} v_1 f~dx + \left.v_1 N_{xx}\right|_{x_a}^{x_b} \\ \int_{x_a}^{x_b} \left[ \cfrac{dv_2}{dx}\left(N_{xx}\cfrac{dw_0}{dx}\right) - \cfrac{d^2v_2}{dx^2}M_{xx}\right]~dx & = \int_{x_a}^{x_b} v_2 q~dx + \left[v_2\left(\cfrac{dM_{xx}}{dx} + N_{xx}\cfrac{dw_0}{dx}\right) \right]_{x_a}^{x_b} - \left[\cfrac{dv_2}{dx}~M_{xx}\right]_{x_a}^{x_b} \end{align} \text{(10)} \qquad $$ Hence shown.

Part 2
The  von Karman strains are related to the displacements by
 * $$\begin{align}

\varepsilon_{xx} & = \varepsilon^0_{xx} + z\varepsilon^1_{xx} \\ \varepsilon^0_{xx} & = \cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2 \\ \varepsilon^1_{xx} & = -\cfrac{d^2w_0}{dx^2} \end{align}$$ The stress and moment resultants are defined as
 * $$\begin{align}

N_{xx} & = \int_A \sigma_{xx}~ dA \\ M_{xx} & = \int_A z\sigma_{xx}~ dA \end{align}$$ For a linear elastic material, the stiffnesses of the beam in extension and bending are defined as
 * $$\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}$$ where $$E$$ is the Young's modulus of the material.

Derive expressions for $$N_{xx}$$ and $$M_{xx}$$ in terms of the displacements $$u_0$$ and $$w_0$$ and the extensional and bending stiffnesses of the beam assuming a linear elastic material.

The stress-strain relations for an isotropic linear elastic material are

\begin{bmatrix} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{zz} \\ \sigma_{yz} \\ \sigma_{zx} \\ \sigma_{xy} \end{bmatrix} = \cfrac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0& 0 & 0 & (1-2\nu) & 0 & 0 \\ 0& 0 & 0 & 0 & (1-2\nu) & 0 \\ 0& 0 & 0 & 0 & 0 & (1-2\nu) \end{bmatrix} \begin{bmatrix} \varepsilon_{xx} \\ \varepsilon_{yy} \\ \varepsilon_{zz} \\ \varepsilon_{yz} \\ \varepsilon_{zx} \\ \varepsilon_{xy} \end{bmatrix} $$ Since all strains other than $$\varepsilon_{xx}$$ are zero, the above equations reduce to

\sigma_{xx} = \cfrac{E(1-\nu)}{(1+\nu)(1-2\nu)}\varepsilon_{xx}; \qquad \sigma_{yy} = \cfrac{E\nu}{(1+\nu)(1-2\nu)}\varepsilon_{xx}; \qquad \sigma_{zz} = \cfrac{E\nu}{(1+\nu)(1-2\nu)}\varepsilon_{xx}~. $$ If we ignore the stresses $$\sigma_{yy}$$ and $$\sigma_{zz}$$, the only allowable value of $$\nu$$ is zero. Then the stress-strain relations become

\sigma_{xx} = E\varepsilon_{xx} ~. $$ Plugging this relation into the stress and moment of stress resultant equations, we get
 * $$\begin{align}

N_{xx} & = \int_A E\varepsilon_{xx}~ dA \\ M_{xx} & = \int_A zE\varepsilon_{xx}~ dA \end{align}$$ Plugging in the relations for the strain we get
 * $$\begin{align}

N_{xx} & = \int_A E(\varepsilon^0_{xx} + z\varepsilon^1_{xx})~ dA   = \int_A E\varepsilon^0_{xx}~dA + \int_A zE\varepsilon^1_{xx}~dA \\ M_{xx} & = \int_A zE(\varepsilon^0_{xx} + z\varepsilon^1_{xx})~ dA   = \int_A zE\varepsilon^0_{xx}~dA + \int_A z^2E\varepsilon^1_{xx}~dA ~. \end{align}$$ Since both $$\varepsilon^0_{xx}$$ and $$\varepsilon^1_{xx}$$ are independent of $$y$$ and $$z$$, we can take these quantities outside the integrals and get
 * $$\begin{align}

N_{xx} & = \varepsilon^0_{xx}\int_A E~dA + \varepsilon^1_{xx}\int_A zE~dA \\ M_{xx} & = \varepsilon^0_{xx}\int_A zE~dA + \varepsilon^1_{xx}\int_A z^2E~dA ~. \end{align}$$ Using the definitions of the extensional, extensional-bending, and bending stiffness, we can then write
 * $$\begin{align}

N_{xx} & = \varepsilon^0_{xx} A_{xx} + \varepsilon^1_{xx} B_{xx} \\ M_{xx} & = \varepsilon^0_{xx} B_{xx} + \varepsilon^1_{xx} D_{xx} ~. \end{align}$$ To write these relations in terms of $$u_0$$ and $$w_0$$ we substitute the expressions for the von Karman strains to get

\begin{align} N_{xx} & = \left[\cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2\right] A_{xx} -\cfrac{d^2w_0}{dx^2} B_{xx} \\ M_{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} D_{xx} ~. \end{align} $$ These are the expressions of the resultants in terms of the displacements.

Part 3
Express the weak forms in terms of the displacements and the extensional and bending stiffnesses.

The weak form equations are
 * $$\begin{align}

\int_{x_a}^{x_b} \cfrac{dv_1}{dx} N_{xx}~dx & = \int_{x_a}^{x_b} v_1 f~dx + \left.v_1 N_{xx}\right|_{x_a}^{x_b} \text{(11)} \qquad \\ \int_{x_a}^{x_b} \left[ \cfrac{dv_2}{dx}\left(N_{xx}\cfrac{dw_0}{dx}\right) - \cfrac{d^2v_2}{dx^2}M_{xx}\right]~dx & = \int_{x_a}^{x_b} v_2 q~dx + \left[v_2\left(\cfrac{dM_{xx}}{dx} + N_{xx}\cfrac{dw_0}{dx}\right) \right]_{x_a}^{x_b} - \left[\cfrac{dv_2}{dx}~M_{xx}\right]_{x_a}^{x_b} \text{(12)} \qquad \end{align}$$ At this stage we make two more assumptions: From the first assumption, we have
 * 1) The elastic modulus is constant throughout the cross-section.
 * 2) The $$x$$-axis passes through the centroid of the cross-section.

B_{xx} = E \int_A z~dA ~. $$ From the second assumption, we get

\int_A z~dA = 0 \qquad \implies \qquad B_{xx} = 0 ~. $$ Then the relations for $$N_{xx}$$ and $$M_{xx}$$ reduce to
 * $$\begin{align}

N_{xx} & = \left[\cfrac{du_0}{dx} + \frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2\right] A_{xx} \text{(13)} \qquad \\ M_{xx} & =- \cfrac{d^2w_0}{dx^2} D_{xx} \text{(14)} \qquad~. \end{align}$$ Let us first consider equation (11). Plugging in the expression for $$N_{xx}$$ we get

\int_{x_a}^{x_b} \cfrac{dv_1}{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} v_1 f~dx + \left.v_1 N_{xx}\right|_{x_a}^{x_b} $$ We can also write the above in terms of virtual displacements by defining $$\delta u_0 := v_1$$, $$Q_1 := -N_{xx}(x_a)$$, and $$Q_4 := N_{xx}(x_b)$$. Then we get

\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 ~. $$ Next, we do the same for equation (12). Plugging in the expressions for $$N_{xx}$$ and $$M_{xx}$$, we get
 * $$\begin{align}

\int_{x_a}^{x_b} \left\{\cfrac{dv_2}{dx}\left(\left[\cfrac{du_0}{dx} + \cfrac{1}{2}~\left(\cfrac{dw_0}{dx}\right)^2\right] A_{xx} \cfrac{dw_0}{dx}\right) \right. & +    \left.\cfrac{d^2 v_2}{dx^2} \left(\cfrac{d^2 w_0}{dx^2}\right) D_{xx} \right\}~dx= \int_{x_a}^{x_b} v_2 q~dx +\\ & \left[v_2\left(\cfrac{dM_{xx}}{dx} + N_{xx}\cfrac{dw_0}{dx}\right)\right]_{x_a}^{x_b} - \left[\cfrac{dv_2}{dx}~M_{xx}\right]_{x_a}^{x_b} \end{align}$$ We can write the above in terms of the virtual displacements and the generalized forces by defining
 * $$\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}$$ to get

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

Part 4
Assume that the approximate solutions for the axial displacement and the transverse deflection over a two noded element are given by
 * $$\begin{align}

u_0(x) & =u_1 \psi_1(x) + u_2 \psi_2(x) \text{(15)} \qquad \\ w_0(x) & =w_1 \phi_1(x) + \theta_1 \phi_2(x) + w_2 \phi_3(x) + \theta_2 \phi_4(x) \text{(16)} \qquad \end{align}$$ where $$\theta = -(dw_0/dx)$$.

Compute the element stiffness matrix for the element.

The weak forms of the governing equations are
 * $$\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 \text{(17)} \qquad \\ \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 ~. \text{(18)} \qquad \end{align}$$

Let us first write the approximate solutions as
 * $$\begin{align}

u_0(x) & = \sum_{j=1}^2 u_j \psi_j \text{(19)} \qquad \\ w_0(x) & = \sum_{j=1}^4 d_j \phi_j \text{(20)} \qquad \end{align}$$ where $$d_j$$ are generalized displacements and

\{d_1, d_2, d_3, d_4\} \equiv \{w_1, \theta_1, w_2, \theta_2\} ~. $$

To formulate the finite element system of equations, we substitute the expressions from $$u_0$$ and $$w_0$$ from equations (19) and (20) into the weak form, and substitute the shape functions $$\psi_i$$ for $$\delta u_0$$, $$\phi_i$$ for $$\delta w_0$$.

For the first equation (17) we get

\int_{x_a}^{x_b} \cfrac{d\psi_i}{dx}\left[ \left(\sum_{j=1}^2 u_j\cfrac{d\psi_j}{dx}\right) + \frac{1}{2}\cfrac{dw_0}{dx} \left(\sum_{j=1}^4 d_j\cfrac{d\phi_j}{dx}\right)\right] A_{xx}~dx = \int_{x_a}^{x_b} \psi_i f~dx + \psi_i(x_a) Q_1 + \psi_i(x_b) Q_4 ~. $$ After reorganizing, we have

\begin{align} \sum_{j=1}^{2} \left[\int_{x_a}^{x_b} A_{xx}\cfrac{d\psi_i}{dx}\cfrac{d\psi_j}{dx}~dx\right]u_j & + \sum_{j=1}^{4} \left[\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\right]d_j = \\ & \int_{x_a}^{x_b} \psi_i f~dx + \psi_i(x_a) Q_1 + \psi_i(x_b) Q_4 ~. \end{align} $$ We can write the above as

\sum_{j=1}^{2} K_{ij}^{11} u_j + \sum_{j=1}^{4} K_{ij}^{12} d_j = F_i^1 $$ where $$i = 1, 2$$ and

\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\\ F_i^1 & = \int_{x_a}^{x_b} \psi_i f~dx + \psi_i(x_a) Q_1 + \psi_i(x_b) Q_4 ~. \end{align} $$

For the second equation (18) we get

\begin{align} \int_{x_a}^{x_b} \left\{\cfrac{d\phi_i}{dx} \right. &    \left[\left(\sum_{j=1}^{2} u_j\cfrac{d\psi_j}{dx}\right) + \frac{1}{2}\cfrac{dw_0}{dx} \left(\sum_{j=1}^{4} d_j\cfrac{d\phi_j}{dx}\right)\right] \cfrac{dw_0}{dx}A_{xx} + \left.\cfrac{d^2\phi_i}{dx^2} \left(\sum_{j=1}^{4} d_j \cfrac{d^2\phi_j}{dx^2}\right)D_{xx} \right\}~dx = \\ & \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} $$ After rearranging we get

\begin{align} \sum_{j=1}^{2} \left[ \int_{x_a}^{x_b} \left(A_{xx} \cfrac{dw_0}{dx}\right) \cfrac{d\phi_i}{dx}\cfrac{d\psi_j}{dx}~dx \right] u_j& + \sum_{j=1}^{4} \left\{\frac{1}{2} \int_{x_a}^{x_b} \left[A_{xx} \left(\cfrac{dw_0}{dx}\right)^2\right] \cfrac{d\phi_i}{dx}\cfrac{d\phi_j}{dx}~dx \right\} d_j\\ & + \sum_{j=1}^{4} \left( \int_{x_a}^{x_b} D_{xx}\cfrac{d^2\phi_i}{dx^2}\cfrac{d^2\phi_j}{dx^2}~dx    \right) d_j = \\ & \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} $$ We can write the above as

\sum_{j=1}^{2} K_{ij}^{21} u_j + \sum_{j=1}^{4} K_{ij}^{22} d_j = F_i^2 $$ where $$i = 1, 2, 3, 4$$ and

\begin{align} K_{ij}^{21} & = \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}\left[A_{xx} \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 \\ 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} $$ In matrix form, we can write

\mathbf{K} = \begin{bmatrix} K_{11}^{11} & K_{12}^{11} & \vdots & K_{11}^{12} & K_{12}^{12} & K_{13}^{12} & K_{14}^{12} \\ K_{21}^{11} & K_{22}^{11} & \vdots & K_{21}^{12} & K_{22}^{12} & K_{23}^{12} & K_{24}^{12} \\ &&\dots &&&&\\ K_{11}^{21} & K_{12}^{21} & \vdots & K_{11}^{22} & K_{12}^{22} & K_{13}^{12} & K_{14}^{12} \\ K_{21}^{21} & K_{22}^{21} & \vdots & K_{21}^{22} & K_{22}^{22} & K_{23}^{12} & K_{24}^{12} \\ K_{31}^{21} & K_{32}^{21} & \vdots & K_{31}^{22} & K_{32}^{22} & K_{33}^{22} & K_{34}^{22} \\ K_{41}^{21} & K_{42}^{21} & \vdots & K_{41}^{22} & K_{42}^{22} & K_{43}^{22} & K_{44}^{22} \\ \end{bmatrix} $$ or

\mathbf{K} = \begin{bmatrix} \mathbf{K}^{11} & \vdots & \mathbf{K}^{12} \\ & \vdots &\\ \mathbf{K}^{21} & \vdots & \mathbf{K}^{22} \end{bmatrix} ~. $$ The finite element system of equations can then be written as
 * $$\text{(21)} \qquad

\begin{bmatrix} K_{11}^{11} & K_{12}^{11} & \vdots & K_{11}^{12} & K_{12}^{12} & K_{13}^{12} & K_{14}^{12} \\ K_{21}^{11} & K_{22}^{11} & \vdots & K_{21}^{12} & K_{22}^{12} & K_{23}^{12} & K_{24}^{12} \\ &&\dots &&&&\\ K_{11}^{21} & K_{12}^{21} & \vdots & K_{11}^{22} & K_{12}^{22} & K_{13}^{12} & K_{14}^{12} \\ K_{21}^{21} & K_{22}^{21} & \vdots & K_{21}^{22} & K_{22}^{22} & K_{23}^{12} & K_{24}^{12} \\ K_{31}^{21} & K_{32}^{21} & \vdots & K_{31}^{22} & K_{32}^{22} & K_{33}^{22} & K_{34}^{22} \\ K_{41}^{21} & K_{42}^{21} & \vdots & K_{41}^{22} & K_{42}^{22} & K_{43}^{22} & K_{44}^{22} \\ \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\\\ d_1 \\ d_2 \\ d_3 \\ d_4 \end{bmatrix} = \begin{bmatrix} F_1^1 \\ F_2^1 \\\\ F_1^2 \\ F_2^2 \\ F_3^2 \\ F_4^2 \end{bmatrix} $$ or

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

Part 5
Show the alternate procedure by which the element stiffness matrix can be made symmetric.

The stiffness matrix is unsymmetric because $$\mathbf{K}^{12}$$ contains a factor of $$1/2$$ while $$\mathbf{K}^{21}$$ does not. The expressions of these terms are
 * $$\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, \qquad i=1,2, ~\text{and}~ j=1,2,3,4 \\ K_{ij}^{21} & = \int_{x_a}^{x_b} \left(A_{xx} \cfrac{dw_0}{dx}\right) \cfrac{d\phi_i}{dx}\cfrac{d\psi_j}{dx}~dx, \qquad i=1,2,3,4, ~\text{and}~ j=1,2 ~. \end{align}$$

To get a symmetric stiffness matrix, we write equation (18) as

\begin{align} \int_{x_a}^{x_b} \left\{\cfrac{d(\delta w_0)}{dx} \right. & \left[\frac{1}{2}\cfrac{du_0}{dx}+\frac{1}{2}\left(\cfrac{dw_0}{dx}\right)^2+ {\frac{1}{2}\cfrac{du_0}{dx}}\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} $$ The quantity is green is assumed to be known from a previous iteration and adds to the $$\mathbf{K}^{22}$$ terms.

Repeating the procedure used in the previous question

\begin{align} \int_{x_a}^{x_b} \left\{\cfrac{d\phi_i}{dx} \right. & \left[\frac{1}{2}\left(\sum_{j=1}^{2} u_j\cfrac{d\psi_j}{dx}\right) + \frac{1}{2}\cfrac{dw_0}{dx} \left(\sum_{j=1}^{4} d_j\cfrac{d\phi_j}{dx}\right) \right] \cfrac{dw_0}{dx}A_{xx} + { \frac{1}{2}\cfrac{d\phi_i}{dx} \cfrac{du_0}{dx}\left(\sum_{j=1}^{4} d_j\cfrac{d\phi_j}{dx}\right)A_{xx} + }\\ &\left.\cfrac{d^2\phi_i}{dx^2} \left(\sum_{j=1}^{4} d_j \cfrac{d^2\phi_j}{dx^2}\right)D_{xx} \right\}~dx = \\ & \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} $$ After rearranging we get

\begin{align} \sum_{j=1}^{2} & \left[\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 \right] u_j + \sum_{j=1}^{4} \left\{\frac{1}{2} \int_{x_a}^{x_b} \left[A_{xx} \left(\cfrac{dw_0}{dx}\right)^2\right] \cfrac{d\phi_i}{dx}\cfrac{d\phi_j}{dx}~dx \right\} d_j + \\ & { \sum_{j=1}^{4} \left[\frac{1}{2}\int_{x_a}^{x_b} A_{xx}\cfrac{du_0}{dx}\cfrac{d\phi_i}{dx} \cfrac{d\phi_j}{dx}~dx\right] d_j } + \sum_{j=1}^{4} \left( \int_{x_a}^{x_b} D_{xx}\cfrac{d^2\phi_i}{dx^2}\cfrac{d^2\phi_j}{dx^2}~dx \right) d_j = \\ & \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} $$

We can write the above as

\sum_{j=1}^{2} K_{ij}^{21} u_j + \sum_{j=1}^{4} K_{ij}^{22} d_j = F_i^2 $$ where $$i = 1, 2, 3, 4$$ and

\begin{align} 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 \\ 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} $$ This gives us a symmetric stiffness matrix.

Part 6
Derive the element tangent stiffness matrix for the element.

Equation (21) can be written as

\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}} ~. $$ The coefficients of the tangent stiffness matrix are given by

T_{ij} = \frac{\partial R_i}{\partial U_j}, \qquad i=1 \dots 6, j=1 \dots 6~. $$ Recall that the finite element system of equations can be written as

\begin{align} \sum_{p=1}^{2} K_{mp}^{11} u_p & + \sum_{q=1}^{4} K_{mq}^{12} d_q = F_m^1, \qquad m=1,2 \\ \sum_{p=1}^{2} K_{np}^{21} u_p & + \sum_{q=1}^{4} K_{nq}^{22} d_q = F_n^2, \qquad n=1,2,3,4 \end{align} $$ where the subscripts have been changed to avoid confusion.

Therefore, the residuals are

\begin{align} R_m^1 & = \sum_{p=1}^{2} K_{mp}^{11} u_p + \sum_{q=1}^{4} K_{mq}^{12} d_q - F_m^1, \qquad m=1,2 \\ R_n^2 & = \sum_{p=1}^{2} K_{np}^{21} u_p + \sum_{q=1}^{4} K_{nq}^{22} d_q - F_n^2, \qquad n=1,2,3,4 ~. \end{align} $$ The derivatives of the residuals with respect to the generalized displacements are

\begin{align} T^{11}_{mk} = \frac{\partial R_m^1}{\partial u_k} & = \sum_{p=1}^{2} \frac{\partial }{\partial u_k}(K_{mp}^{11} u_p) + \sum_{q=1}^{4} \frac{\partial }{\partial u_k}(K_{mq}^{12} d_q) , & & \qquad m=1,2; k=1,2 \\ T^{12}_{ml} = \frac{\partial R_m^1}{\partial d_l} & = \sum_{p=1}^{2} \frac{\partial }{\partial d_l}(K_{mp}^{11} u_p) + \sum_{q=1}^{4} \frac{\partial }{\partial d_l}(K_{mq}^{12} d_q) , & & \qquad m=1,2; l=1,2,3,4 \\ T^{21}_{nk} = \frac{\partial R_n^2}{\partial u_k} & = \sum_{p=1}^{2} \frac{\partial }{\partial u_k}(K_{np}^{21} u_p) + \sum_{q=1}^{4} \frac{\partial }{\partial u_k}(K_{nq}^{22} d_q) , & & \qquad n=1,2,3,4; k=1,2 \\ T^{22}_{nl} = \frac{\partial R_n^2}{\partial d_l} & = \sum_{p=1}^{2} \frac{\partial }{\partial d_l}(K_{np}^{21} u_p) + \sum_{q=1}^{4} \frac{\partial }{\partial d_l}(K_{nq}^{22} d_q) , & & \qquad n=1,2,3,4; l=1,2,3,4 ~. \end{align} $$ Differentiating, we get

\begin{align} T^{11}_{mk} & = \sum_{p=1}^{2} \left(u_p\frac{\partial K_{mp}^{11}}{\partial u_k} + K_{mp}^{11}\frac{\partial u_p}{\partial u_k}\right) + \sum_{q=1}^{4} \left(d_q\frac{\partial K_{mq}^{12}}{\partial u_k} + K_{mq}^{12}\frac{\partial d_q}{\partial u_k}\right), & &\qquad m=1,2; k=1,2 \\ T^{12}_{ml} & = \sum_{p=1}^{2} \left(u_p\frac{\partial K_{mp}^{11}}{\partial d_l} + K_{mp}^{11}\frac{\partial u_p}{\partial d_l}\right) + \sum_{q=1}^{4} \left(d_q\frac{\partial K_{mq}^{12}}{\partial d_l} + K_{mq}^{12}\frac{\partial d_q}{\partial d_l}\right), & &\qquad m=1,2; l=1,2,3,4 \\ T^{21}_{nk} & = \sum_{p=1}^{2} \left(u_p\frac{\partial K_{np}^{21}}{\partial u_k} + K_{np}^{21}\frac{\partial u_p}{\partial u_k}\right) + \sum_{q=1}^{4} \left(d_q\frac{\partial K_{nq}^{22}}{\partial u_k} + K_{nq}^{22}\frac{\partial d_q}{\partial u_k}\right), & &\qquad n=1,2,3,4; k=1,2 \\ T^{22}_{nl} & = \sum_{p=1}^{2} \left(u_p\frac{\partial K_{np}^{21}}{\partial d_l} + K_{np}^{21}\frac{\partial u_p}{\partial d_l}\right) + \sum_{q=1}^{4} \left(d_q\frac{\partial K_{nq}^{22}}{\partial d_l} + K_{nq}^{22}\frac{\partial d_q}{\partial d_l}\right), & &\qquad n=1,2,3,4; l=1,2,3,4 ~. \end{align} $$ These equations can therefore be written as

\begin{align} T^{11}_{mk} & = K^{11}_{mk} + \sum_{p=1}^{2} u_p\frac{\partial K_{mp}^{11}}{\partial u_k} + \sum_{q=1}^{4} d_q\frac{\partial K_{mq}^{12}}{\partial u_k}, & &\qquad m=1,2; k=1,2 \\ T^{12}_{ml} & = K^{12}_{ml} + \sum_{p=1}^{2} u_p\frac{\partial K_{mp}^{11}}{\partial d_l} + \sum_{q=1}^{4} d_q\frac{\partial K_{mq}^{12}}{\partial d_l}, & &\qquad m=1,2; l=1,2,3,4 \\ T^{21}_{nk} & = K^{21}_{nk} + \sum_{p=1}^{2} u_p\frac{\partial K_{np}^{21}}{\partial u_k} + \sum_{q=1}^{4} d_q\frac{\partial K_{nq}^{22}}{\partial u_k}, & &\qquad n=1,2,3,4; k=1,2 \\ T^{22}_{nl} & = K^{22}_{nl} + \sum_{p=1}^{2} u_p\frac{\partial K_{np}^{21}}{\partial d_l} + \sum_{q=1}^{4} d_q\frac{\partial K_{nq}^{22}}{\partial d_l}, & &\qquad n=1,2,3,4; l=1,2,3,4 ~. \end{align} $$ Now, the coefficients of $$\mathbf{K}^{11}$$, $$\mathbf{K}^{12}$$, and $$\mathbf{K}^{21}$$ of the symmetric stiffness matrix are independent of $$u_1$$ and $$u_2$$. Also, the terms of $$\mathbf{K}^{11}$$ are independent of the all the generalized displacements. Therefore, the above equations reduce to

\begin{align} T^{11}_{mk} & = K^{11}_{mk}, & &\qquad m=1,2; k=1,2 \\ T^{12}_{ml} & = K^{12}_{ml} + \sum_{q=1}^{4} d_q\frac{\partial K_{mq}^{12}}{\partial d_l}, & & \qquad m=1,2; l=1,2,3,4 \\ T^{21}_{nk} & = K^{21}_{nk} + \sum_{q=1}^{4} d_q\frac{\partial K_{nq}^{22}}{\partial u_k}, & & \qquad n=1,2,3,4; k=1,2 \\ T^{22}_{nl} & = K^{22}_{nl} + \sum_{p=1}^{2} u_p\frac{\partial K_{np}^{21}}{\partial d_l} + \sum_{q=1}^{4} d_q\frac{\partial K_{nq}^{22}}{\partial d_l}, & & \qquad n=1,2,3,4; l=1,2,3,4 ~. \end{align} $$

Consider the coefficients of $$\mathbf{T}^{12}$$:

T^{12}_{ml} = K^{12}_{ml} + \sum_{q=1}^{4} d_q\frac{\partial K_{mq}^{12}}{\partial d_l}, \qquad m=1,2; l=1,2,3,4 ~. $$ From our previous derivation, we have

K_{mq}^{12} =\frac{1}{2} \int_{x_a}^{x_b} \left(A_{xx}\cfrac{dw_0}{dx}\right) \cfrac{d\psi_m}{dx} \cfrac{d\phi_q}{dx}~dx~. $$ Therefore,
 * $$\begin{align}

\frac{\partial K_{mq}^{12}}{\partial d_l} &= \frac{1}{2}\int_{x_a}^{x_b}\left[A_{xx}\frac{\partial }{\partial d_l}\left(\cfrac{dw_0}{dx}\right)\right] \cfrac{d\psi_m}{dx} \cfrac{d\phi_q}{dx}~dx \\ & = \frac{1}{2}\int_{x_a}^{x_b}\left[A_{xx}\left(\sum_{T=1}^{4}\frac{\partial d_T}{\partial d_l} \cfrac{d\phi_T}{dx}\right)\right] \cfrac{d\psi_m}{dx} \cfrac{d\phi_q}{dx}~dx \\ & = \frac{1}{2}\int_{x_a}^{x_b}\left[A_{xx}\cfrac{d\phi_l}{dx}\right] \cfrac{d\psi_m}{dx} \cfrac{d\phi_q}{dx}~dx \end{align}$$ The tangent stiffness matrix coefficients are therefore
 * $$\begin{align}

{ T^{12}_{ml}} & = K^{12}_{ml} + \sum_{q=1}^{4} d_q\left\{ \frac{1}{2}\int_{x_a}^{x_b}\left[A_{xx}\cfrac{d\phi_l}{dx}\right] \cfrac{d\psi_m}{dx} \cfrac{d\phi_q}{dx}~dx \right\} \qquad m=1,2; l=1,2,3,4 ~. \\ & = K^{12}_{ml} + \frac{1}{2}\int_{x_a}^{x_b} A_{xx}\cfrac{d\phi_l}{dx} \cfrac{d\psi_m}{dx}\left(\sum_{q=1}^{4}d_q\cfrac{d\phi_q}{dx}\right)~dx \\ & = K^{12}_{ml} + \frac{1}{2}\int_{x_a}^{x_b} A_{xx}\cfrac{d\phi_l}{dx} \cfrac{d\psi_m}{dx}\cfrac{dw_0}{dx}~dx\\ & = { 2K^{12}_{ml}} ~. \end{align}$$ Next, {consider the coefficients of $$\mathbf{T}^{21}$$:

T^{21}_{nk} = K^{21}_{nk} + \sum_{q=1}^{4} d_q\frac{\partial K_{nq}^{22}}{\partial u_k} ~. $$ The coefficients of $$\mathbf{K}^{22}$$ are

K_{nq}^{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_n}{dx}\cfrac{d\phi_q}{dx} + D_{xx}\cfrac{d^2\phi_n}{dx^2}\cfrac{d^2\phi_q}{dx^2}\right\}~dx ~. $$ Therefore, the derivatives are
 * $$\begin{align}

\frac{\partial K_{nq}^{22}}{\partial u_k} & = \int_{x_a}^{x_b} \frac{1}{2} A_{xx} \left[\sum_{T=1}^{2}\frac{\partial u_T}{\partial u_k}\cfrac{d\psi_T}{dx}+ 2\cfrac{dw_0}{dx}\left(\sum_{T=1}^{4}\frac{\partial d_T}{\partial u_k}\cfrac{d\phi_T}{dx} \right)\right] \cfrac{d\phi_n}{dx}\cfrac{d\phi_q}{dx}~dx \\ &= \frac{1}{2}\int_{x_a}^{x_b} A_{xx} \cfrac{d\psi_k}{dx} \cfrac{d\phi_n}{dx}\cfrac{d\phi_q}{dx}~dx \\ \end{align}$$ Therefore the coefficients of $$\mathbf{T}^{21}$$ are
 * $$\begin{align}

{ T^{21}_{nk}} & = K^{21}_{nk} + \sum_{q=1}^{4} d_q\left(\frac{1}{2}\int_{x_a}^{x_b} A_{xx} \cfrac{d\psi_k}{dx} \cfrac{d\phi_n}{dx}\cfrac{d\phi_q}{dx}~dx\right) \\ & = K^{21}_{nk} + \frac{1}{2}\int_{x_a}^{x_b} A_{xx}\cfrac{d\psi_k}{dx}\cfrac{d\phi_n}{dx} \left(\sum_{q=1}^{4} d_q\cfrac{d\phi_q}{dx}\right)~dx \\ & = K^{21}_{nk} + \frac{1}{2}\int_{x_a}^{x_b} A_{xx}\cfrac{d\psi_k}{dx}\cfrac{d\phi_n}{dx} \cfrac{dw_0}{dx}~dx \\ & = { 2K^{21}_{nk}} ~. \end{align}$$ Finally, for the $$\mathbf{T}^{22}$$ coefficients, we start with

T^{22}_{nl} = K^{22}_{nl} + \sum_{p=1}^{2} u_p\frac{\partial K_{np}^{21}}{\partial d_l} + \sum_{q=1}^{4} d_q\frac{\partial K_{nq}^{22}}{\partial d_l}, \qquad n=1,2,3,4; l=1,2,3,4 $$ and plug in the derivatives of the stiffness matrix coefficients
 * $$\begin{align}

K_{np}^{21} & = \frac{1}{2} \int_{x_a}^{x_b} \left(A_{xx} \cfrac{dw_0}{dx}\right) \cfrac{d\phi_n}{dx}\cfrac{d\psi_p}{dx}~dx\\ K_{nq}^{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_n}{dx}\cfrac{d\phi_q}{dx} + D_{xx}\cfrac{d^2\phi_n}{dx^2}\cfrac{d^2\phi_q}{dx^2}\right\}~dx ~. \end{align}$$ \\ The derivatives are
 * $$\begin{align}

\frac{\partial K_{np}^{21}}{\partial d_l} & = \frac{1}{2} \int_{x_a}^{x_b} \left[A_{xx} \frac{\partial }{\partial d_l}\left(\cfrac{dw_0}{dx}\right) \right] \cfrac{d\phi_n}{dx}\cfrac{d\psi_p}{dx}~dx\\ & = \frac{1}{2} \int_{x_a}^{x_b} A_{xx}\left[\sum_{T=1}^{4}\frac{\partial d_T}{\partial d_l}\cfrac{d\phi_T}{dx} \right] \cfrac{d\phi_n}{dx}\cfrac{d\psi_p}{dx}~dx \\ & = \frac{1}{2} \int_{x_a}^{x_b} A_{xx}\cfrac{d\phi_l}{dx}\cfrac{d\phi_n}{dx} \cfrac{d\psi_p}{dx}~dx \end{align}$$ and
 * $$\begin{align}

\frac{\partial K_{nq}^{22}}{\partial d_l} & = \int_{x_a}^{x_b} \frac{1}{2} A_{xx} \left[\sum_{T=1}^{2}\frac{\partial u_T}{\partial d_l}\cfrac{d\psi_T}{dx}+ 2\cfrac{dw_0}{dx}\left(\sum_{T=1}^{4} \frac{\partial d_T}{\partial d_l}\cfrac{d\phi_T}{dx}\right)\right] \cfrac{d\phi_n}{dx}\cfrac{d\phi_q}{dx}~dx \\ & = \int_{x_a}^{x_b} A_{xx} \cfrac{dw_0}{dx} \cfrac{d\phi_l}{dx} \cfrac{d\phi_n}{dx}\cfrac{d\phi_q}{dx}~dx \end{align}$$ Therefore, the coefficients of the tangent stiffness matrix can be written as
 * $$\begin{align}

{ T^{22}_{nl}} & = K^{22}_{nl} + \sum_{p=1}^{2} u_p \left[ \frac{1}{2} \int_{x_a}^{x_b} A_{xx}\cfrac{d\phi_l}{dx}\cfrac{d\phi_n}{dx} \cfrac{d\psi_p}{dx}~dx \right] + \sum_{q=1}^{4} d_q \left[ \int_{x_a}^{x_b} A_{xx} \cfrac{dw_0}{dx} \cfrac{d\phi_l}{dx} \cfrac{d\phi_n}{dx}\cfrac{d\phi_q}{dx}~dx \right]\\ & = K^{22}_{nl} + \frac{1}{2} \int_{x_a}^{x_b} A_{xx}\cfrac{d\phi_l}{dx}\cfrac{d\phi_n}{dx} \left(\sum_{p=1}^{2} u_p \cfrac{d\psi_p}{dx}\right)~dx+ \int_{x_a}^{x_b} A_{xx} \cfrac{dw_0}{dx} \cfrac{d\phi_l}{dx} \cfrac{d\phi_n}{dx} \left(\sum_{q=1}^{4} d_q \cfrac{d\phi_q}{dx}\right)~dx\\ & = K^{22}_{nl} + \frac{1}{2} \int_{x_a}^{x_b} A_{xx}\cfrac{d\phi_l}{dx}\cfrac{d\phi_n}{dx} \cfrac{du_0}{dx}~dx + \int_{x_a}^{x_b} A_{xx} \cfrac{dw_0}{dx} \cfrac{d\phi_l}{dx} \cfrac{d\phi_n}{dx} \cfrac{dw_0}{dx}~dx \\ & = { K^{22}_{nl} + \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_l}{dx}\cfrac{d\phi_n}{dx}~dx }~. \end{align}$$

The final expressions for the tangent stiffness matrix terms are

\begin{align} T^{11}_{ij} & = K^{11}_{ij}, & &\qquad i=1,2; j=1,2 \\ T^{12}_{ij} & = 2 K^{12}_{ij}, & & \qquad i=1,2; j=1,2,3,4 \\ T^{21}_{ij} & = 2 K^{21}_{ij}, & & \qquad i=1,2,3,4; j=1,2 \\ 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 & & \qquad i=1,2,3,4; j=1,2,3,4 ~. \end{align} $$