# Nonlinear finite elements/Homework 5/Solutions/Problem 1

## Problem 1: Nonlinear Beam Bending

Given:

The differential equations governing the bending of straight beams are

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle v_{1}}$ be the weighting function for equation (1) and let ${\displaystyle v_{2}}$ be the weighting function for equation (2).

Then the weak forms of the two equations are

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle E}$ is the Young's modulus of the material.

Derive expressions for ${\displaystyle N_{xx}}$ and ${\displaystyle M_{xx}}$ in terms of the displacements ${\displaystyle u_{0}}$ and ${\displaystyle 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

${\displaystyle {\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 ${\displaystyle \varepsilon _{xx}}$ are zero, the above equations reduce to

${\displaystyle \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 ${\displaystyle \sigma _{yy}}$ and ${\displaystyle \sigma _{zz}}$, the only allowable value of ${\displaystyle \nu }$ is zero. Then the stress-strain relations become

${\displaystyle \sigma _{xx}=E\varepsilon _{xx}~.}$

Plugging this relation into the stress and moment of stress resultant equations, we get

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle \varepsilon _{xx}^{0}}$ and ${\displaystyle \varepsilon _{xx}^{1}}$ are independent of ${\displaystyle y}$ and ${\displaystyle z}$, we can take these quantities outside the integrals and get

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle u_{0}}$ and ${\displaystyle w_{0}}$ we substitute the expressions for the von Karman strains to get

{\displaystyle {\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

{\displaystyle {\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:

1. The elastic modulus is constant throughout the cross-section.
2. The ${\displaystyle x}$-axis passes through the centroid of the cross-section.

From the first assumption, we have

${\displaystyle B_{xx}=E\int _{A}z~dA~.}$

From the second assumption, we get

${\displaystyle \int _{A}z~dA=0\qquad \implies \qquad B_{xx}=0~.}$

Then the relations for ${\displaystyle N_{xx}}$ and ${\displaystyle M_{xx}}$ reduce to

{\displaystyle {\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 ${\displaystyle N_{xx}}$ we get

${\displaystyle \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 ${\displaystyle \delta u_{0}:=v_{1}}$, ${\displaystyle Q_{1}:=-N_{xx}(x_{a})}$, and ${\displaystyle Q_{4}:=N_{xx}(x_{b})}$. Then we get

${\displaystyle \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 ${\displaystyle N_{xx}}$ and ${\displaystyle M_{xx}}$, we get

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle \theta =-(dw_{0}/dx)}$.

Compute the element stiffness matrix for the element.

The weak forms of the governing equations are

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle d_{j}}$ are generalized displacements and

${\displaystyle \{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 ${\displaystyle u_{0}}$ and ${\displaystyle w_{0}}$ from equations (19) and (20) into the weak form, and substitute the shape functions ${\displaystyle \psi _{i}}$ for ${\displaystyle \delta u_{0}}$, ${\displaystyle \phi _{i}}$ for ${\displaystyle \delta w_{0}}$.

For the first equation (17) we get

${\displaystyle \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

{\displaystyle {\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

${\displaystyle \sum _{j=1}^{2}K_{ij}^{11}u_{j}+\sum _{j=1}^{4}K_{ij}^{12}d_{j}=F_{i}^{1}}$

where ${\displaystyle i=1,2}$ and

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

${\displaystyle \sum _{j=1}^{2}K_{ij}^{21}u_{j}+\sum _{j=1}^{4}K_{ij}^{22}d_{j}=F_{i}^{2}}$

where ${\displaystyle i=1,2,3,4}$ and

{\displaystyle {\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

${\displaystyle \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

${\displaystyle \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

${\displaystyle {\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

${\displaystyle {\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 ${\displaystyle \mathbf {K} ^{12}}$ contains a factor of ${\displaystyle 1/2}$ while ${\displaystyle \mathbf {K} ^{21}}$ does not. The expressions of these terms are

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle \mathbf {K} ^{22}}$ terms.

Repeating the procedure used in the previous question

{\displaystyle {\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

{\displaystyle {\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

${\displaystyle \sum _{j=1}^{2}K_{ij}^{21}u_{j}+\sum _{j=1}^{4}K_{ij}^{22}d_{j}=F_{i}^{2}}$

where ${\displaystyle i=1,2,3,4}$ and

{\displaystyle {\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

${\displaystyle \mathbf {K} (\mathbf {U} )\mathbf {U} =\mathbf {F} }$

where

{\displaystyle {\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

${\displaystyle \mathbf {R} =\mathbf {K} \mathbf {U} -\mathbf {F} ~.}$

For Newton iterations, we use the algorithm

${\displaystyle \mathbf {U} ^{r+1}=\mathbf {U} ^{r}-(\mathbf {T} ^{r})^{-1}\mathbf {R} ^{r}}$

where the tangent stiffness matrix is given by

${\displaystyle \mathbf {T} ^{r}={\frac {\partial \mathbf {R} ^{r}}{\partial \mathbf {U} }}~.}$

The coefficients of the tangent stiffness matrix are given by

${\displaystyle 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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle \mathbf {K} ^{11}}$, ${\displaystyle \mathbf {K} ^{12}}$, and ${\displaystyle \mathbf {K} ^{21}}$ of the symmetric stiffness matrix are independent of ${\displaystyle u_{1}}$ and ${\displaystyle u_{2}}$. Also, the terms of ${\displaystyle \mathbf {K} ^{11}}$ are independent of the all the generalized displacements. Therefore, the above equations reduce to

{\displaystyle {\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 ${\displaystyle \mathbf {T} ^{12}}$:

${\displaystyle 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

${\displaystyle 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,

{\displaystyle {\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

{\displaystyle {\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 ${\displaystyle \mathbf {T} ^{21}}$:

${\displaystyle T_{nk}^{21}=K_{nk}^{21}+\sum _{q=1}^{4}d_{q}{\frac {\partial K_{nq}^{22}}{\partial u_{k}}}~.}$

The coefficients of ${\displaystyle \mathbf {K} ^{22}}$ are

${\displaystyle 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

{\displaystyle {\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 ${\displaystyle \mathbf {T} ^{21}}$ are

{\displaystyle {\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 ${\displaystyle \mathbf {T} ^{22}}$ coefficients, we start with

${\displaystyle 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

{\displaystyle {\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

{\displaystyle {\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

{\displaystyle {\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