home

PDF (letter size)

Chapter 1
Introduction

A short review for solving the beam problem in 2D is given. The deflection curve, bending moment and shear force diagrams are calculated for a beam subject to bending moment and shear force using direct stiffness method and then using finite elements method by adding more elements. The problem is solved first by finding the stiffness matrix using the direct method and then using the virtual work method.

Chapter 2
Direct method

Starting with only one element beam which is subject to bending and shear forces. There are 4 nodal degrees of freedom. Rotation at the left and right nodes of the beam and transverse displacements at the left and right nodes. The following diagram shows the sign convention used for external forces. Moments are always positive when anti-clockwise direction and vertical forces are positive when in the positive $y$ direction.

The two nodes are numbered $1$ and $2$ from left to right. $M_1$ is the moment at the left node (node 1), $M_2$ is the moment at the right node (node 2). $V_1$ is the vertical force at the left node and $V_2$ is the vertical force at the right node.

The above diagram shows the signs used for the applied forces direction when acting in the positive sense. Since this is a one dimensional problem, the displacement field (the unknown being solved for) will be a function of one independent variable which is the $x$ coordinate. The displacement field in the vertical direction is called $v\left(x\right)$. This is the vertical displacement of point $x$ on the beam from the original $x-axis$. The following diagram shows the notation used for the coordinates.

Angular displacement at distance $x$ on the beam is found using $\theta \left(x\right)=\frac{dv\left(x\right)}{dx}.$ At the left node, the degrees of freedom or the displacements, are called $v_1,{\theta }_1$ and at the right node they are called $v_2,{\theta }_2$. At an arbitrary location $x$ in the beam, the vertical displacement is $v\left(x\right)$ and the rotation at that location is $\theta \left(x\right)$.

The following diagram shows the displacement field $v\left(x\right)$

In the direct method of finding the stiffness matrix, the forces at the ends of the beam are found directly by the use of beam theory. In beam theory the signs are different from what is given in the first diagram above. Therefore, the moment and shear forces obtained using beam theory ($M_B$ and $V_B$ in the diagram below) will have different signs when compared to the external forces. The signs are then adjusted to reflect the convention as shown in the diagram above using $M$ and $V$.

For an example, the external moment $M_1$ is opposite in sign to $M_{B1}$ and the reaction $V_2$ is opposite to $V_{B2}$. To illustrate this more, a diagram with both sign conventions is given below.

The goal now is to obtain expressions for external loads $M_i$ and $R_i$ in the above diagram as function of the displacements at the nodes $\left\{d\right\}=\{v_1,{\theta }_1,v_2,{\theta }_2{\}}^T$.

In other words, the goal is to obtain an expression of the form $\left\{p\right\}=\left[K\right]\left\{d\right\}$ where $\left[K\right]$ is the stiffness matrix where $\left\{p\right\}={\{V_1,M_1,R_2,M_2\}}^T$ is the nodal forces or load vector, and $\left\{d\right\}$ is the nodal displacement vector.

In this case $\left[K\right]$ will be a $4\times 4$ matrix and $\left\{p\right\}$ is a $4\times 1$ vector and $\left\{d\right\}$ is a $4\times 1$ vector.

Starting with $V_1$. It is in the same direction as the shear force $V_{B1}$. Since $V_{B1}=\frac{d{M}_{B1}}{dx}$ then$$V_1=\frac{d{M}_{B1}}{dx}$$ Since from beam theory $M_{B1}=-\sigma \left(x\right)\frac{I}{y}$, the above becomes$$V_1=-\frac{I}{y}\frac{d\sigma \left(x\right)}{dx}$$ But $\sigma \left(x\right)=E\epsilon \left(x\right)$ and $\epsilon \left(x\right)=\frac{-y}{\rho }$ where $\rho$ is radius of curvature, therefore the above becomes$$V_1=EI\frac{d}{dx}\left(\frac{1}{\rho }\right)$$ Since $\frac{1}{\rho }=\frac{\left(\frac{d^{2}u}{d{x}^2}\right)}{{(1+{\left(\frac{du}{dx}\right)}^2)}^{3/2}}$ and for a small angle of deflection $\frac{du}{dx}\ll 1$ then $\frac{1}{\rho }=\left(\frac{d^{2}u}{d{x}^2}\right)$, and the above now becomes$$V_1=EI\frac{d^{3}u\left(x\right)}{d{x}^3}$$ Before continuing, the following diagram illustrates the above derivation. This comes from beam theory.

Now $M_1$ is obtained. $M_1$ is in the opposite sense of the bending moment $M_{B1}$ hence a negative sign is added giving $M_1=-M_{B1}$. But $M_{B1}=-\sigma \left(x\right)\frac{I}{y}$ therefore$$\begin{array}{rl}{M}_1& =\sigma \left(x\right)\frac{I}{y}\\ & =E\epsilon \left(x\right)\frac{I}{y}\\ & =E\left(\frac{-y}{\rho }\right)\frac{I}{y}\\ & =-EI\left(\frac{1}{\rho }\right)\\ & =-EI\frac{d^{2}w}{d{x}^2}\end{array}$$

$V_2$ is now found. It is in the opposite sense of the shear force $V_{B2}$, hence a negative sign is added giving

$$\begin{array}{rl}{V}_2& =-V_{B2}\\ & =-\frac{d{M}_{B2}}{dx}\end{array}$$

Since $M_{B2}=-\sigma \left(x\right)\frac{I}{y}$, the above becomes$$V_2=\frac{I}{y}\frac{d\sigma \left(x\right)}{dx}$$ But $\sigma \left(x\right)=E\epsilon \left(x\right)$ and $\epsilon \left(x\right)=\frac{-y}{\rho }$ where $\rho$ is radius of curvature. The above becomes$$V_2=-EI\frac{d}{dx}\left(\frac{1}{\rho }\right)$$ But $\frac{1}{\rho }=\frac{\left(\frac{d^{2}w}{d{x}^2}\right)}{{(1+{\left(\frac{dw}{dx}\right)}^2)}^{3/2}}$ and for small angle of deflection $\frac{dw}{dx}\ll 1$ hence $\frac{1}{\rho }=\left(\frac{d^{2}u}{d{x}^2}\right)$ then the above becomes$$V_2=-EI\frac{d^{3}u\left(x\right)}{d{x}^3}$$ Finally $M_2$ is in the same direction as $M_{B2}$ so no sign change is needed. $M_{B2}=-\sigma \left(x\right)\frac{I}{y}$, therefore$$\begin{array}{rl}{M}_2& =-\sigma \left(x\right)\frac{I}{y}\\ & =-E\epsilon \left(x\right)\frac{I}{y}\\ & =-E\left(\frac{-y}{\rho }\right)\frac{I}{y}\\ & =EI\left(\frac{1}{\rho }\right)\\ & =EI\frac{d^{2}u}{d{x}^2}\end{array}$$

The following is a summary of what was found so far. Notice that the above expressions are evaluates at $x=0$ and at $x=L$. Accordingly to obtain the nodal end forces vector $\left\{p\right\}$$$\begin{array}{}\text{(1)}& \left\{p\right\}=\left\{\begin{array}{c}{V}_1\\ M_1\\ V_2\\ M_2\end{array}\right\}=\left\{\begin{array}{c}EI{\frac{d^{3}u\left(x\right)}{d{x}^3}|}_{x=0}\\ -EI{\frac{d^{2}u}{d{x}^2}|}_{x=0}\\ -EI{\frac{d^{3}u\left(x\right)}{d{x}^3}|}_{x=L}\\ EI{\frac{d^{2}u}{d{x}^2}|}_{x=L}\end{array}\right\}\end{array}$$ The RHS of the above is now expressed as function of the nodal displacements $v_1,{\theta }_1,v_2,{\theta }_2$. To do that, the field displacement $v\left(x\right)$ which is the transverse displacement of the beam is assumed to be a polynomial in $x$ of degree $3$ or$$\begin{array}{rl}v\left(x\right)& =a_0+a_{1}x+a_2{x}^2+a_3{x}^3\\ \text{(A)}& \theta \left(x\right)& =\frac{dv\left(x\right)}{dx}=a_1+2a_{2}x+3a_3{x}^2\end{array}$$

Polynomial of degree $3$ is used since there are $4$ degrees of freedom, and a minimum of $4$ free parameters is needed. Hence$$\begin{array}{}\text{(2)}& v_1={v\left(x\right)|}_{x=0}=a_0\end{array}$$ And $$\begin{array}{}\text{(3)}& {\theta }_1={\theta \left(x\right)|}_{x=0}=a_1\end{array}$$ Assuming the length of the beam is $L$, then$$\begin{array}{}\text{(4)}& v_2={v\left(x\right)|}_{x=L}=a_0+a_{1}L+a_2{L}^2+a_3{L}^3\end{array}$$ And$$\begin{array}{}\text{(5)}& {\theta }_2={\theta \left(x\right)|}_{x=L}=a_1+2a_{2}L+3a_3{L}^2\end{array}$$ Equations (2-5) gives$$\left\{d\right\}=\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}=\left\{\begin{array}{c}{a}_0\\ a_1\\ a_0+a_{1}L+a_2{L}^2+a_3{L}^3\\ a_1+2a_{2}L+3a_3{L}^2\end{array}\right\}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 1& L& L^2& L^3\\ 0& 1& 2L& 3L^2\end{array}\right)\left\{\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\end{array}\right\}$$ Solving for $a_i$ gives$$\begin{array}{rl}\left\{\begin{array}{c}{a}_0\\ a_1\\ a_2\\ a_3\end{array}\right\}& =\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ -\frac{3}{L^2}& -\frac{2}{L}& \frac{3}{L^2}& -\frac{1}{L}\\ \frac{2}{L^3}& \frac{1}{L^2}& -\frac{2}{L^3}& \frac{1}{L^2}\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\\ & =\left(\begin{array}{c}{v}_1\\ {\theta }_1\\ \frac{3}{L^2}{v}_2-\frac{1}{L}{\theta }_2-\frac{3}{L^2}{v}_1-\frac{2}{L}{\theta }_1\\ \frac{1}{L^2}{\theta }_1+\frac{1}{L^2}{\theta }_2+\frac{2}{L^3}{v}_1-\frac{2}{L^3}{v}_2\end{array}\right)\end{array}$$

$v\left(x\right)$, the field displacement function from Eq. (A) can now be written as a function of the nodal displacements$$\begin{array}{rl}v\left(x\right)& =a_0+a_{1}x+a_2{x}^2+a_3{x}^3\\ & =v_1+{\theta }_{1}x+(\frac{3}{L^2}{v}_2-\frac{1}{L}{\theta }_2-\frac{3}{L^2}{v}_1-\frac{2}{L}{\theta }_1)x^2+(\frac{1}{L^2}{\theta }_1+\frac{1}{L^2}{\theta }_2+\frac{2}{L^3}{v}_1-\frac{2}{L^3}{v}_2)x^3\\ & =v_1+x{\theta }_1-2\frac{x^2}{L}{\theta }_1+\frac{x^3}{L^2}{\theta }_1-\frac{x^2}{L}{\theta }_2+\frac{x^3}{L^2}{\theta }_2-3\frac{x^2}{L^2}{v}_1+2\frac{x^3}{L^3}{v}_1+3\frac{x^2}{L^2}{v}_2-2\frac{x^3}{L^3}{v}_2\end{array}$$

Or in matrix form$$\begin{array}{rl}v\left(x\right)& =\left(\begin{array}{cccc}1-3\frac{x^2}{L^2}+2\frac{x^3}{L^3}& x-2\frac{x^2}{L}+\frac{x^3}{L^2}& 3\frac{x^2}{L^2}-2\frac{x^3}{L^3}& -\frac{x^2}{L}+\frac{x^3}{L^2}\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\\ & =\left(\begin{array}{cccc}\frac{1}{L^3}(L^3-3L{x}^2+2x^3)& \frac{1}{L^2}(L^{2}x-2L{x}^2+x^3)& \frac{1}{L^3}(3L{x}^2-2x^3)& \frac{1}{L^2}(-L{x}^2+x^3)\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\end{array}$$

The above can be written as$$\begin{array}{rl}\text{(5A)}& v\left(x\right)& =\left(\begin{array}{cccc}{N}_1\left(x\right)& N_2\left(x\right)& N_3\left(x\right)& N_4\left(x\right)\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}v\left(x\right)& =\left[N\right]\left\{d\right\}\end{array}$$

Where $N_i$ are called the shape functions. The shape functions are$$\left(\begin{array}{c}{N}_1\left(x\right)\\ N_2\left(x\right)\\ N_3\left(x\right)\\ N_4\left(x\right)\end{array}\right)=\left(\begin{array}{c}\frac{1}{L^3}(L^3-3L{x}^2+2x^3)\\ \frac{1}{L^2}(L^{2}x-2L{x}^2+x^3)\\ \frac{1}{L^3}(3L{x}^2-2x^3)\\ \frac{1}{L^2}(-L{x}^2+x^3)\end{array}\right)$$ We notice that $N_1\left(0\right)=1$ and $N_1\left(L\right)=0$ as expected. Also $${\frac{d{N}_2\left(x\right)}{dx}|}_{x=0}={\frac{1}{L^2}(L^2-4Lx+3x^2)|}_{x=0}=1$$ And $${\frac{d{N}_2\left(x\right)}{dx}|}_{x=L}={\frac{1}{L^2}(L^2-4Lx+3x^2)|}_{x=L}=0$$ Also $N_3\left(0\right)=0$ and $N_3\left(L\right)=1$ and $${\frac{d{N}_4\left(x\right)}{dx}|}_{x=0}={\frac{1}{L^2}(-2Lx+3x^2)|}_{x=0}=0$$ and $${\frac{d{N}_4\left(x\right)}{dx}|}_{x=L}={\frac{1}{L^2}(-2Lx+3x^2)|}_{x=L}=1$$ The shape functions have thus been verified. The stiffness matrix is now found by substituting Eq. (5A) into Eq. (1), repeated below$$\left\{p\right\}=\left\{\begin{array}{c}{V}_1\\ M_1\\ V_2\\ M_2\end{array}\right\}=\left(\begin{array}{c}EI{\frac{d^{3}v\left(x\right)}{d{x}^3}|}_{x=0}\\ -EI{\frac{d^{2}v}{d{x}^2}|}_{x=0}\\ -EI{\frac{d^{3}v\left(x\right)}{d{x}^3}|}_{x=L}\\ EI{\frac{d^{2}v}{d{x}^2}|}_{x=L}\end{array}\right)$$ Hence$$\begin{array}{}\text{(6)}& \left\{p\right\}=\left\{\begin{array}{c}{V}_1\\ M_1\\ V_2\\ M_2\end{array}\right\}=\left(\begin{array}{c}EI\frac{d^3}{d{x}^3}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)\\ -EI\frac{d^2}{d{x}^2}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)\\ -EI\frac{d^3}{d{x}^3}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)\\ EI\frac{d^2}{d{x}^2}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)\end{array}\right)\end{array}$$ But $$\begin{array}{rl}\frac{d^3}{d{x}^3}{N}_1\left(x\right)& =\frac{1}{L^3}\frac{d^3}{d{x}^3}(L^3-3L{x}^2+2x^3)\\ & =\frac{12}{L^3}\end{array}$$

And$$\begin{array}{rl}\frac{d^3}{d{x}^3}{N}_2\left(x\right)& =\frac{1}{L^2}\frac{d^3}{d{x}^3}(L^{2}x-2L{x}^2+x^3)\\ & =\frac{6}{L^2}\end{array}$$

And$$\begin{array}{rl}\frac{d^3}{d{x}^3}{N}_3\left(x\right)& =\frac{1}{L^3}\frac{d^3}{d{x}^3}(3L{x}^2-2x^3)\\ & =\frac{-12}{L^3}\end{array}$$

And$$\begin{array}{rl}\frac{d^3}{d{x}^3}{N}_4\left(x\right)& =\frac{1}{L^2}\frac{d^3}{d{x}^3}(-L{x}^2+x^3)\\ & =\frac{6}{L^2}\end{array}$$

For the second derivatives$$\begin{array}{rl}\frac{d^2}{d{x}^2}{N}_1\left(x\right)& =\frac{1}{L^3}\frac{d^2}{d{x}^2}(L^3-3L{x}^2+2x^3)\\ & =\frac{1}{L^3}(12x-6L)\end{array}$$

And$$\begin{array}{rl}\frac{d^2}{d{x}^2}{N}_2\left(x\right)& =\frac{1}{L^2}\frac{d^2}{d{x}^2}(L^{2}x-2L{x}^2+x^3)\\ & =\frac{1}{L^2}(6x-4L)\end{array}$$

And$$\begin{array}{rl}\frac{d^2}{d{x}^2}{N}_3\left(x\right)& =\frac{1}{L^3}\frac{d^2}{d{x}^2}(3L{x}^2-2x^3)\\ & =\frac{1}{L^3}(6L-12x)\end{array}$$

And$$\begin{array}{rl}\frac{d^2}{d{x}^2}{N}_4\left(x\right)& =\frac{1}{L^2}\frac{d^2}{d{x}^2}(-L{x}^2+x^3)\\ & =\frac{1}{L^2}(6x-2L)\end{array}$$

Hence Eq. (6) becomes$$\begin{array}{rl}\left\{p\right\}& =\left\{\begin{array}{c}{V}_1\\ M_1\\ V_2\\ M_2\end{array}\right\}=\left(\begin{array}{c}EI{\frac{d^3}{d{x}^3}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)|}_{x=0}\\ -EI{\frac{d^2}{d{x}^2}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)|}_{x=0}\\ -EI{\frac{d^3}{d{x}^3}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)|}_{x=L}\\ EI{\frac{d^2}{d{x}^2}(N_1{v}_1+N_2{\theta }_1+N_3{v}_2+N_4{\theta }_2)|}_{x=L}\end{array}\right)\\ & =\left(\begin{array}{c}EI{(\frac{12}{L^3}{v}_1+\frac{6}{L^2}{\theta }_1-\frac{12}{L^3}{v}_2+\frac{6}{L^2}{\theta }_2)}_{x=0}\\ -EI{(\frac{1}{L^3}(12x-6L)v_1+\frac{1}{L^2}(6x-4L){\theta }_1+\frac{1}{L^3}(6L-12x)v_2+\frac{1}{L^2}(6x-2L){\theta }_2)}_{x=0}\\ -EI{(\frac{12}{L^3}{v}_1+\frac{6}{L^2}{\theta }_1-\frac{12}{L^3}{v}_2+\frac{6}{L^2}{\theta }_2)}_{x=L}\\ EI{(\frac{1}{L^3}(12x-6L)v_1+\frac{1}{L^2}(6x-4L){\theta }_1+\frac{1}{L^3}(6L-12x)v_2+\frac{1}{L^2}(6x-2L){\theta }_2)}_{x=L}\end{array}\right)\\ \text{(7)}& & =\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ -{(12x-6L)}_{x=0}& -L{(6x-4L)}_{x=0}& -{(6L-12x)}_{x=0}& -L{(6x-2L)}_{x=0}\\ -12& -6L& 12& -6L\\ {(12x-6L)}_{x=L}& L{(6x-4L)}_{x=L}& {(6L-12x)}_{x=L}& L{(6x-2L)}_{x=L}\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\end{array}$$

Or in matrix form, after evaluating the expressions above for $x=L$ and $x=0$ as$$\begin{array}{rl}\left\{\begin{array}{c}{V}_1\\ M_1\\ V_2\\ M_2\end{array}\right\}& =\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -\frac{12}{L^3}& -\frac{6}{L^2}& \frac{12}{L^3}& -\frac{6}{L^2}\\ (12L-6L)& L(6L-4L)& (6L-12L)& L(6L-2L)\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\\ & =\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\end{array}$$

The above now is in the form$$\left\{p\right\}=\left[K\right]\left\{d\right\}$$ Hence the stiffness matrix is$$\left[K\right]=\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right)$$ Knowing the stiffness matrix means knowing the nodal displacements $\left\{d\right\}$ when given the forces at the nodes. The power of the finite element method now comes after all the nodal displacements $v_1,{\theta }_1,$ $v_2,{\theta }_2$ are calculated by solving $\left\{p\right\}=\left[K\right]\left\{d\right\}$. This is because the polynomial $v\left(x\right)$ is now completely determined and hence $v\left(x\right)$ and $\theta \left(x\right)$ can now be evaluated for any $x$ along the beam and not just at its end nodes as the case with finite difference method. Eq. 5A above can now be used to find the displacement $v\left(x\right)$ and $\theta \left(x\right)$ everywhere.$$\begin{array}{rl}v\left(x\right)& =\left[N\right]\left\{d\right\}\\ v\left(x\right)& =\left(\begin{array}{cccc}{N}_1\left(x\right)& N_2\left(x\right)& N_3\left(x\right)& N_4\left(x\right)\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\end{array}$$

To summarise, these are the steps to obtain $v\left(x\right)$

1. An expression for $\left[K\right]$ is found.
2. $\left\{p\right\}=\left[K\right]\left\{d\right\}$ is solved for $\left\{d\right\}$
3. $v\left(x\right)=\left[N\right]\left\{d\right\}$ is calculated by assuming $v\left(x\right)$ is a polynomial. This gives the displacement $v\left(x\right)$ to use to evaluate the transverse displacement anywhere on the beam and not just at the end nodes.
4. $\theta \left(x\right)=\frac{dv\left(x\right)}{dx}=\frac{d}{dx}\left[N\right]\left\{d\right\}$ is obtained to evaluate the rotation of the beam any where and not just at the end nodes.
5. The strain $ϵ\left(x\right)=-y\left[B\right]\left\{d\right\}$ is found, where $\left[B\right]$ is the gradient matrix $\left[B\right]=\frac{d^2}{d{x}^2}\left[N\right]$.
6. The stress from $\sigma =Eϵ=-Ey\left[B\right]\left\{d\right\}$ is found.
7. The bending moment diagram from $M\left(x\right)=EI\left[B\right]\left\{d\right\}$ is found.
8. The shear force diagram from $V\left(x\right)=\frac{d}{dx}M\left(x\right)$ is found.

2.1 Examples using the direct beam stiffness matrix

The beam stiffness matrix is now used to solve few beam problems. Starting with simple one span beam

2.1.1 Example 1

A one span beam, a cantilever beam of length $L$, with point load $P$ at the free end

The first step is to make the free body diagram and show all moments and forces at the nodes

$P$ is the given force. $M_2=0$ since there is no external moment at the right end. Hence $\left\{p\right\}=\left[K\right]\left\{d\right\}$ for this system is$$\left\{\begin{array}{c}R\\ M_1\\ -P\\ 0\end{array}\right\}=\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}$$ Now is an important step. The known end displacements from boundary conditions is substituted into $\left\{d\right\}$, and the corresponding row and columns from the above system of equations are removed¹. Boundary conditions indicates that there is no rotation on the left end (since fixed). Hence $v_1=0$ and ${\theta }_1=0$. Hence the only unknowns are $v_2$ and ${\theta }_2$. Therefore the first and the second rows and columns are removed, giving$$\left\{\begin{array}{c}-P\\ 0\end{array}\right\}=\frac{EI}{L^3}\left(\begin{array}{cc}12& -6L\\ -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}{v}_2\\ {\theta }_2\end{array}\right\}$$ Now the above is solved for $\left\{\begin{array}{c}{v}_2\\ {\theta }_2\end{array}\right\}$. Let $E=30\times {10}^6$ psi (steel), $I=57$ $i{n}^4$, $L=144in$, and $P=400$ lb, hence

P=400; 
L=144; 
E=30*10^6; 
I=57.1; 
A=(E*I/L^3)*[12       6*L            -12      6*L; 
6*L      4*L^2          -6*L     2*L^2; 
-12     -6*L            12      -6*L; 
6*L      2*L^2          -6*L     4*L^2 
]; 
load=[P;P*L;-P;0]; 
x=A(3:end,3:end)\load(3:end)
 

which gives

x = 
-0.2324 
-0.0024
 

Therefore the vertical displacement at the right end is $v_2=0.2324$ inches (downwards) and ${\theta }_2=-0.0024$ radians. Now that all nodal displacements are found, the field displacement function is completely determined.$$\left\{d\right\}=\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ -0.2324\\ -0.0024\end{array}\right\}$$ From Eq. 5A$$\begin{array}{rl}v\left(x\right)& =\left[N\right]\left\{d\right\}\\ & =\left(\begin{array}{cccc}{N}_1\left(x\right)& N_2\left(x\right)& N_3\left(x\right)& N_4\left(x\right)\end{array}\right)\left\{\begin{array}{c}0\\ 0\\ -0.2324\\ -0.0024\end{array}\right\}\\ & =\left(\begin{array}{cccc}\frac{1}{L^3}(L^3-3L{x}^2+2x^3)& \frac{1}{L^2}(L^{2}x-2L{x}^2+x^3)& \frac{1}{L^3}(3L{x}^2-2x^3)& \frac{1}{L^2}(-L{x}^2+x^3)\end{array}\right)\left\{\begin{array}{c}0\\ 0\\ -0.2324\\ -0.0024\end{array}\right\}\\ & =\frac{0.0024}{L^2}(L{x}^2-x^3)+\frac{0.2324}{L^3}(2x^3-3L{x}^2)\end{array}$$

But $L=144$ inches, and the above becomes$$v\left(x\right)=3.9920\times {10}^{-8}{x}^3-1.6956\times {10}^{-5}{x}^2$$ To verify, let $x=144$ in the above$$\begin{array}{rl}v(x=144)& =3.9920\times {10}^{-8}{144}^3-1.6956\times {10}^{-5}{144}^2\\ & =-0.23240\end{array}$$

The following is a plot of the deflection curve for the beam

v=@(x) 3.992*10^-8*x.^3-1.6956*10^-5*x.^2 
x=0:0.1:144; 
plot(x,v(x),'r-','LineWidth',2); 
ylim([-0.8 0.3]); 
title('beam deflection curve'); 
xlabel('x inch'); ylabel('deflection inch'); 
grid
 

Now instead of removing rows/columns for the known boundary conditions, a $1$ is put on the diagonal. Starting again with$$\left\{\begin{array}{c}R\\ M_1\\ -P\\ 0\end{array}\right\}=\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}$$ Since $v_1=0$ and ${\theta }_1=0$, then$$\left\{\begin{array}{c}0\\ 0\\ -P\\ 0\end{array}\right\}=\frac{EI}{L^3}\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 12& -6L\\ 0& 0& -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}0\\ 0\\ v_2\\ {\theta }_2\end{array}\right\}$$ The above system is now solved as before. $E=30\times {10}^6$ psi, $I=57$ $i{n}^4$, $L=144$in, $P=400$ lb

P=400; 
L=144; 
E=30*10^6; 
I=57.1; 
A=(E*I/L^3)*[12       6*L            -12      6*L; 
6*L      4*L^2          -6*L     2*L^2; 
-12     -6*L            12      -6*L; 
6*L      2*L^2          -6*L     4*L^2 
]; 
load=[0;0;-P;0];  %put zeros for known B.C. 
A(:,1)=0; A(1,:)=0; A(1,1)=1;  %put 1 on diagonal 
A(:,2)=0; A(2,:)=0; A(2,2)=1;  %put 1 on diagonal 
A
 

Gives

A = 
1.0e+07 * 
0.0000         0         0         0 
0    0.0000         0         0 
0         0    0.0007   -0.0496 
0         0   -0.0496    4.7583
 

Then

sol=A\load  %SOLVE
                                                                                    
                                                                                    
 

Gives

sol = 
0 
0 
-0.2324 
-0.0024
 

The same solution is obtained as before, but without the need to remove rows/column from the stiffness matrix. This method might be easier for programming than the first method of removing rows/columns.

The rest now is the same as was done earlier and will not be repeated.

2.1.2 Example 2

This is the same example as above, but the vertical load $P$ is now placed in the middle of the beam

In using stiffness method, all loads must be on the nodes. The vector $\left\{p\right\}$ is the nodal forces vector. Hence equivalent nodal loads are found for the load in the middle of the beam. The equivalent loading is the following

Therefore, the problem is as if it was the following problem

Now that equivalent loading is in place, we continue as before. Making a free body diagram showing all loads (including reaction forces)

The stiffness equation is now written as$$\begin{array}{rl}\left\{p\right\}& =\left[K\right]\left\{d\right\}\\ \left\{\begin{array}{c}R-P/2\\ M_1-PL/8\\ -P/2\\ PL/8\end{array}\right\}& =\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}\end{array}$$

There is no need to determine $R$ and $M_1$ at this point since these rows will be removed due to boundary conditions $v_1=0$ and ${\theta }_1=0$ and hence those quantities are not needed to solve the equations. Note that the rows and columns are removed for the known boundary displacements before solving $\left\{p\right\}=\left[K\right]\left\{d\right\}$. Hence, after removing the first two rows and columns, the above system simplifies to$$\left\{\begin{array}{c}-P/2\\ PL/8\end{array}\right\}=\frac{EI}{L^3}\left(\begin{array}{cc}12& -6L\\ -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}{v}_2\\ {\theta }_2\end{array}\right\}$$ The above is now solved for $\left\{\begin{array}{c}{v}_2\\ {\theta }_2\end{array}\right\}$ using the same numerical values for $P,E,I,L$ as in the first example

P=400; 
L=144; 
E=30*10^6; 
I=57.1; 
A=(E*I/L^3)*[12       6*L            -12      6*L; 
6*L      4*L^2          -6*L     2*L^2; 
-12     -6*L            12      -6*L; 
6*L      2*L^2          -6*L     4*L^2 
]; 
load=[-P/2;P*L/8]; 
x=A(3:end,3:end)\load
 

Gives

x = 
-0.072630472854641 
-0.000605253940455
 

Therefore$$\left\{d\right\}=\left\{\begin{array}{c}{v}_1\\ {\theta }_1\\ v_2\\ {\theta }_2\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ -0.072630472854641\\ -0.000605253940455\end{array}\right\}$$ This is enough to obtain $v\left(x\right)$ as before. Now the reactions $R$ and $M_1$ can be determined if needed. Going back to the full $\left\{p\right\}=\left[K\right]\left\{d\right\}$, results in$$\begin{array}{rl}\left\{\begin{array}{c}R-P/2\\ M_1-PL/8\\ -P/2\\ PL/8\end{array}\right\}& =\frac{EI}{L^3}\left(\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right)\left\{\begin{array}{c}0\\ 0\\ -0.072630472854641\\ -0.000605253940455\end{array}\right\}\\ & =\frac{EI}{L^3}\left(\begin{array}{c}0.87157-3.6315\times {10}^{-3}L\\ 0.43578L-1.2105\times {10}^{-3}{L}^2\\ 3.6315\times {10}^{-3}L-0.87157\\ 0.43578L-2.421\times {10}^{-3}{L}^2\end{array}\right)\end{array}$$