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

### Solution
#### Part 1

Show that the weak forms of these equations can be written as
    ${\begin{aligned}\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^{2}v_{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{aligned}}$
First we get rid of the shear force term by combining the second and third equations to get
    ${\begin{aligned}{\cfrac {dN_{xx}}{dx}}+f(x)&=0{\text{(1)}}\qquad \\{\cfrac {d^{2}M_{xx}}{dx^{2}}}+q(x)+{\cfrac {d}{dx}}\left(N_{xx}{\cfrac {dw_{0}}{dx}}\right)&=0{\text{(2)}}\qquad \end{aligned}}$
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{aligned}\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^{2}M_{xx}}{dx^{2}}}+q(x)+{\cfrac {d}{dx}}\left(N_{xx}{\cfrac {dw_{0}}{dx}}\right)\right]~dx&=0{\text{(4)}}\qquad \end{aligned}}$
To get the symmetric weak forms, we integrate by parts (even though the symmetry is not obvious at this stage) to get
    ${\begin{aligned}\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{aligned}}$
Integrate equation (6) again by parts, and get
    ${\begin{aligned}\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^{2}v_{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{aligned}}$
Collect terms and rearrange equations (5) and (7) to get
    ${\begin{aligned}\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^{2}v_{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{aligned}}$
Rewriting the equations, we get
    ${\begin{aligned}\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^{2}v_{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{aligned}}{\text{(10)}}\qquad$
Hence shown.

#### Part 2

The von Karman strains are related to the displacements by
    ${\begin{aligned}\varepsilon _{xx}&=\varepsilon _{xx}^{0}+z\varepsilon _{xx}^{1}\\\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 stress and moment resultants are defined as
    ${\begin{aligned}N_{xx}&=\int _{A}\sigma _{xx}~dA\\M_{xx}&=\int _{A}z\sigma _{xx}~dA\end{aligned}}$
For a linear elastic material, the stiffnesses of the beam in extension and bending are defined as
    ${\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}}$
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{aligned}N_{xx}&=\int _{A}E\varepsilon _{xx}~dA\\M_{xx}&=\int _{A}zE\varepsilon _{xx}~dA\end{aligned}}$
Plugging in the relations for the strain we get
    ${\begin{aligned}N_{xx}&=\int _{A}E(\varepsilon _{xx}^{0}+z\varepsilon _{xx}^{1})~dA=\int _{A}E\varepsilon _{xx}^{0}~dA+\int _{A}zE\varepsilon _{xx}^{1}~dA\\M_{xx}&=\int _{A}zE(\varepsilon _{xx}^{0}+z\varepsilon _{xx}^{1})~dA=\int _{A}zE\varepsilon _{xx}^{0}~dA+\int _{A}z^{2}E\varepsilon _{xx}^{1}~dA~.\end{aligned}}$
Since both $\varepsilon _{xx}^{0}$ and $\varepsilon _{xx}^{1}$ are independent of $y$ and $z$, we can take these quantities outside the integrals and get
    ${\begin{aligned}N_{xx}&=\varepsilon _{xx}^{0}\int _{A}E~dA+\varepsilon _{xx}^{1}\int _{A}zE~dA\\M_{xx}&=\varepsilon _{xx}^{0}\int _{A}zE~dA+\varepsilon _{xx}^{1}\int _{A}z^{2}E~dA~.\end{aligned}}$
Using the definitions of the extensional, extensional-bending, and bending stiffness, we can then write
    ${\begin{aligned}N_{xx}&=\varepsilon _{xx}^{0}A_{xx}+\varepsilon _{xx}^{1}B_{xx}\\M_{xx}&=\varepsilon _{xx}^{0}B_{xx}+\varepsilon _{xx}^{1}D_{xx}~.\end{aligned}}$
To write these relations in terms of $u_{0}$ and $w_{0}$ we substitute the expressions for the von Karman strains to get
    ${\begin{aligned}N_{xx}&=\left[{\cfrac {du_{0}}{dx}}+{\frac {1}{2}}\left({\cfrac {dw_{0}}{dx}}\right)^{2}\right]A_{xx}-{\cfrac {d^{2}w_{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^{2}w_{0}}{dx^{2}}}D_{xx}~.\end{aligned}}$
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{aligned}\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^{2}v_{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{aligned}}$
At this stage we make two more assumptions:
* The elastic modulus is constant throughout the cross-section.
* The $x$-axis passes through the centroid of the cross-section.
From the first assumption, we have
    $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{aligned}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^{2}w_{0}}{dx^{2}}}D_{xx}{\text{(14)}}\qquad ~.\end{aligned}}$
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{aligned}\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{aligned}}$
We can write the above in terms of the virtual displacements and the generalized forces by defining
    ${\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}}$
to get
    ${\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}}$

#### Part 4

Assume that the approximate solutions for the axial displacement and the transverse deflection over a two noded element are given by
    ${\begin{aligned}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{aligned}}$
where $\theta =-(dw_{0}/dx)$.
Compute the element stiffness matrix for the element.
The weak forms of the governing equations are
    ${\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}{\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^{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}~.{\text{(18)}}\qquad \end{aligned}}$
Let us first write the approximate solutions as
    ${\begin{aligned}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{aligned}}$
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{aligned}\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{aligned}}$
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{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\\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{aligned}}$
For the second equation (18) we get
    ${\begin{aligned}\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{aligned}}$
After rearranging we get
    ${\begin{aligned}\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{aligned}}$
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{aligned}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{aligned}}$
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{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,\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{aligned}}$
To get a symmetric stiffness matrix, we write equation (18) as
    ${\begin{aligned}\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^{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}}$
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{aligned}\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{aligned}}$
After rearranging we get
    ${\begin{aligned}\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{aligned}}$
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{aligned}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{aligned}}$
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{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
    $\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{aligned}\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{aligned}}$
where the subscripts have been changed to avoid confusion.
Therefore, the residuals are
    ${\begin{aligned}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{aligned}}$
The derivatives of the residuals with respect to the generalized displacements are
    ${\begin{aligned}T_{mk}^{11}={\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_{ml}^{12}={\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_{nk}^{21}={\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_{nl}^{22}={\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{aligned}}$
Differentiating, we get
    ${\begin{aligned}T_{mk}^{11}&=\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_{ml}^{12}&=\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_{nk}^{21}&=\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_{nl}^{22}&=\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{aligned}}$
These equations can therefore be written as
    ${\begin{aligned}T_{mk}^{11}&=K_{mk}^{11}+\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_{ml}^{12}&=K_{ml}^{12}+\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_{nk}^{21}&=K_{nk}^{21}+\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_{nl}^{22}&=K_{nl}^{22}+\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{aligned}}$
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{aligned}T_{mk}^{11}&=K_{mk}^{11},&&\qquad m=1,2;~~k=1,2\\T_{ml}^{12}&=K_{ml}^{12}+\sum _{q=1}^{4}d_{q}{\frac {\partial K_{mq}^{12}}{\partial d_{l}}},&&\qquad m=1,2;~~l=1,2,3,4\\T_{nk}^{21}&=K_{nk}^{21}+\sum _{q=1}^{4}d_{q}{\frac {\partial K_{nq}^{22}}{\partial u_{k}}},&&\qquad n=1,2,3,4;~~k=1,2\\T_{nl}^{22}&=K_{nl}^{22}+\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{aligned}}$
Consider the coefficients of $\mathbf {T} ^{12}$:
    $T_{ml}^{12}=K_{ml}^{12}+\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{aligned}{\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{aligned}}$
The tangent stiffness matrix coefficients are therefore
    ${\begin{aligned}{T_{ml}^{12}}&=K_{ml}^{12}+\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_{ml}^{12}+{\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_{ml}^{12}+{\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_{ml}^{12}}~.\end{aligned}}$
Next, {consider the coefficients of $\mathbf {T} ^{21}$:
    $T_{nk}^{21}=K_{nk}^{21}+\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{aligned}{\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{aligned}}$
Therefore the coefficients of $\mathbf {T} ^{21}$ are
    ${\begin{aligned}{T_{nk}^{21}}&=K_{nk}^{21}+\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_{nk}^{21}+{\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_{nk}^{21}+{\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_{nk}^{21}}~.\end{aligned}}$
Finally, for the $\mathbf {T} ^{22}$ coefficients, we start with
    $T_{nl}^{22}=K_{nl}^{22}+\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{aligned}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{aligned}}$ \\
The derivatives are
    ${\begin{aligned}{\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{aligned}}$
and
    ${\begin{aligned}{\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{aligned}}$
Therefore, the coefficients of the tangent stiffness matrix can be written as
    ${\begin{aligned}{T_{nl}^{22}}&=K_{nl}^{22}+\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_{nl}^{22}+{\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_{nl}^{22}+{\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_{nl}^{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 _{l}}{dx}}{\cfrac {d\phi _{n}}{dx}}~dx}~.\end{aligned}}$
The final expressions for the tangent stiffness matrix terms are
    ${\begin{aligned}T_{ij}^{11}&=K_{ij}^{11},&&\qquad i=1,2;~~j=1,2\\T_{ij}^{12}&=2K_{ij}^{12},&&\qquad i=1,2;~~j=1,2,3,4\\T_{ij}^{21}&=2K_{ij}^{21},&&\qquad i=1,2,3,4;~~j=1,2\\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&&\qquad i=1,2,3,4;j=1,2,3,4~.\end{aligned}}$