Table of contents

Truss analysis overview

The simple truss is represented by two rods that are attached to each other in the right corner (node) using a pin connector which is assumed to be moment free. The other two ends of the rods are firmly fixed in place, which is denoted by the two triangles. A point force is placed at the node, denoted by PP. The rods can only carry axial force, meaning that loads cannot be applied anywhere but at the nodes. This problem is a 2D planar problem, so every node has two degrees of freedom i.e., it can be translated in the xx and yy-direction.

We begin our analysis by numbering each rod and each node. Each rod has four quantities of interest: normal forces, NiN_i, normal stresses, σi\sigma_i, strains, εi\varepsilon_i and deformations δi\delta_i where i=1,2i=1,2. Additionally each node has two quantities of interest: nodal forces fjf_j and nodal displacements uju_j where j=1,2,3j=1,2,3.

For a structure like this, the interesting questions are:

Is the structure able to carry the load? The the structural requirements are typically posed on the maximum stress and displacements as well as reaction forces, i.e., nodal forces at the boundary conditions.

An overview of the relations between nodal quantities and rod quantities is given by the figure below

These relations essentially form the equation

δ=FLEA \delta=\dfrac{FL}{EA}

which is the simplest solution (setting f=0f=0) to the boundary value problem:

BVP{ddx(EAdudx)=ffor x[0,L]u(0)=LEAdudx(L)=F \boxed{ \text{BVP}\begin{cases} -\dfrac{d}{dx}\left(EA\dfrac{du}{dx}\right)=f & \text{for }x\in[0,L]\\ u(0)=L\\ EA\dfrac{du}{dx}(L)=F \end{cases} }

This is why the method in this (and the following) section is also known as the direct approach, since we are tackling the differential equation directly, but also with a loss of generality.

In what follows, we shall derive a method for determining the displacements uu of all nodes given a set of external forces ff.

For our problem, we can see that we have four boundary conditions: u2x=u2y=u3x=u3y=0u_2^x =u_2^y =u_3^x =u_3^y =0. The known nodal loads are f1y=Pf_1^y =-P. Furthermore we know the rod quantities E1=E2=EE_1 =E_2 =E, A1=A2=AA_1 =A_2 =A, L1=LL_1 =L and L2=2LL_2 =\sqrt{2}L.

We want to know the displacements in node 1, the reaction forces in nodes, 2 and 3 and all rod normal forces, stresses, strains and deformations.

The node 1 is said to be free, i.e., able to be translated in xx and yy. We say that node 1 has two degrees of freedom.

We now introduce the displacement field (vector) which corresponds to each degree of freedom.

u=[u1xu1y] \def\arraystretch{1.5} \bm u=\left\lbrack \begin{array}{c} u_1^x \\ u_1^y \end{array}\right\rbrack

Additionally, we introduce the nodal load vector

f=[f1xf1y] \def\arraystretch{1.5} \bm f =\left\lbrack \begin{array}{c} f_1^x \\ f_1^y \end{array}\right\rbrack

We will replace the known loads later, i.e., we will evaluate f1y=Pf_1^y =-P at a later time.

The problem is to find a relation between the loads f\bm f and the displacements u\bm u. This is done by computing the stiffness matrix for the structure such that

Su=fu=S1f\boldsymbol S \boldsymbol u = \bm f \Rightarrow \boldsymbol u = \boldsymbol S^{-1} \bm f

Equilibrium of an isolated node leads to the matrix AT{\mathit{\mathbf{A} } }^T directly, this approach is known as method of joints.

Setting up material relations leads to the Se{\mathit{\mathbf{S} } }^e matrix directly.

Thus having AT\boldsymbol A^T and Se\boldsymbol S^e gives

S=ATSeAu=S1f\boldsymbol S = \boldsymbol A^T \boldsymbol S^e \boldsymbol A \Rightarrow \boldsymbol u = \boldsymbol S^{-1} \bm f

Once the displacements are known, we can compute the rod quantities:

δ=Au \boldsymbol \delta = \boldsymbol A \boldsymbol u, N=Seδ\boldsymbol N = \boldsymbol S^e \boldsymbol \delta and finally σi=Ni/Ai\sigma_i =N_i /A_i

Example 1

Step 1: Number all rod and nodes

Step 2: Introduce a coordinate system

Step 3: Introduce the displacement vector and load vector

u=[u1xu1y]\mathit{\mathbf{u} }=\left\lbrack \begin{array}{c} u_1^x \\ u_1^y \end{array}\right\rbrack f=[f1xf1y]\mathit{\mathbf{f} }=\left\lbrack \begin{array}{c} f_1^x \\ f_1^y \end{array}\right\rbrack

Step 4: Rod fields

δ=[δ1δ2δ3],N=[N1N2N3]\boldsymbol \delta =\left\lbrack \begin{array}{c} \delta_1 \\ \delta_2 \\ \delta_3 \end{array}\right\rbrack ,\mathit{\mathbf{N} }=\left\lbrack \begin{array}{c} N_1 \\ N_2 \\ N_3 \end{array}\right\rbrack

Step 5: Determine either AT{\mathit{\mathbf{A} } }^T by equilibrium using the method of joints or determine A\mathit{\mathbf{A} } using the unit displacements method. We shall do both and you can decide what method you prefer.

Displacement method

Use the deformation relation

δ=A  u\boldsymbol \delta =\mathit{\mathbf{A} }\;\mathit{\mathbf{u} }

first, set u=[10]:\mathit{\mathbf{u} }=\left\lbrack \begin{array}{c} 1\\ 0 \end{array}\right\rbrack :

[δ1δ2δ3]=[A11A12A21A22A31A32][10]=[A11A21A31]\left\lbrack \begin{array}{c} \delta_1 \\ \delta_2 \\ \delta_3 \end{array}\right\rbrack =\left\lbrack \begin{array}{cc} A_{11} & A_{12} \\ A_{21} & A_{22} \\ A_{31} & A_{32} \end{array}\right\rbrack \left\lbrack \begin{array}{c} 1\\ 0 \end{array}\right\rbrack =\left\lbrack \begin{array}{c} A_{11} \\ A_{21} \\ A_{31} \end{array}\right\rbrack

i.e., we get the first column of A\mathit{\mathbf{A} }, likewise, setting u=[01]T\mathit{\mathbf{u} }={\left\lbrack \begin{array}{cc} 0 & 1 \end{array}\right\rbrack }^T will give us the second column.

Now we apply this using a deformation figure

Assuming small deformations (small deformation theory) all rods must be only translated such that they stay parallell to each other.

We look at the deformation in the xx-direction:

δ1x=1δ2x=cos45°=12δ3x=0\begin{array}{l} \delta_1^x =1\\ \delta_2^x =\textrm{cos45}\degree =\dfrac{1}{\sqrt{2} }\\ \delta_3^x =0 \end{array}

Now in the yy-direction:

δ1y=0δ2y=sin45°=12δ3y=1\begin{array}{l} \delta_1^y =0\\ \delta_2^y =-\textrm{sin45}\degree =-\dfrac{1}{\sqrt{2} }\\ \delta_3^y =-1 \end{array}

Note, here the deformation are negative, since we are compressing the rods.

These deformations lead to

A=[δ1xδ1yδ2xδ2yδ3xδ3y]=[10121201]\Rightarrow A=\left\lbrack \begin{array}{cc} \delta_1^x & \delta_1^y \\ \delta_2^x & \delta_2^y \\ \delta_3^x & \delta_3^y \end{array}\right\rbrack =\left\lbrack \begin{array}{cc} 1 & 0\\ \dfrac{1}{\sqrt{2} } & -\dfrac{1}{\sqrt{2} }\\ 0 & -1 \end{array}\right\rbrack

This matrix is purely geometrically determined!

Force method (method of joints)

f=ATN\mathit{\mathbf{f} }={\mathit{\mathbf{A} } }^T \mathit{\mathbf{N} }

equilibrium of the node (joint):

Note: Be aware of the directions considering the coordinate system!

:f1xN1N2cos45°=0:f1y+N3+N2sin45°=0\begin{array}{l} \to :f_1^x -N_1 -N_2 \textrm{cos45}\degree =0\\ \downarrow :f_1^y +N_3 +N_2 \textrm{sin45}\degree =0 \end{array} f1x=N1+N212f1y=N212N3\begin{array}{l} \Rightarrow f_1^x =N_1 +N_2 \dfrac{1}{\sqrt{2} }\\ f_1^y =-N_2 \dfrac{1}{\sqrt{2} }-N_3 \end{array}

On matrix form:

[f1xf1y]=[11200121]AT[N1N2N3]\left\lbrack \begin{array}{c} f_1^x \\ f_1^y \end{array}\right\rbrack =\underset{ {\mathit{\mathbf{A} } }^T }{\underbrace{\left\lbrack \begin{array}{ccc} 1 & \dfrac{1}{\sqrt{2} } & 0\\ 0 & -\dfrac{1}{\sqrt{2} } & -1 \end{array}\right\rbrack } } \left\lbrack \begin{array}{c} N_1 \\ N_2 \\ N_3 \end{array}\right\rbrack

Note: Determine either A\mathit{\mathbf{A} } or AT{\mathit{\mathbf{A} } }^T using which ever method you feel is easier and then just transpose!.

Step 6: Material relations, N=Seδ\mathit{\mathbf{N} }={\mathit{\mathbf{S} } }^e \boldsymbol \delta

For each rod we have

δi=NiLiEiAiNi=EiAiLiδi\delta_i =\dfrac{N_i L_i }{E_i A_i }\Rightarrow N_i =\dfrac{E_i A_i }{L_i }\delta_i

or on matrix form

[N1N2N3]=[E1A1L1000E2A2L2000E3A3L3]Se[δ1δ2δ3]\left\lbrack \begin{array}{c} N_1 \\ N_2 \\ N_3 \end{array}\right\rbrack =\underset{ {\mathit{\mathbf{S} } }^e }{\underbrace{\left\lbrack \begin{array}{ccc} \frac{E_1 A_1 }{L_1 } & 0 & 0\\ 0 & \frac{E_2 A_2 }{L_2 } & 0\\ 0 & 0 & \frac{E_3 A_3 }{L_3 } \end{array}\right\rbrack } } \left\lbrack \begin{array}{c} \delta_1 \\ \delta_2 \\ \delta_3 \end{array}\right\rbrack

or

[N1N2N3]=[E1A1L1E2A2L2E3A3L3][100010001]ISe[δ1δ2δ3]\left\lbrack \begin{array}{c} N_1 \\ N_2 \\ N_3 \end{array}\right\rbrack =\underset{ {\mathit{\mathbf{S} } }^e }{\underbrace{\left\lbrack \begin{array}{ccc} \dfrac{E_1 A_1 }{L_1 } & \dfrac{E_2 A_2 }{L_2 } & \dfrac{E_3 A_3 }{L_3 } \end{array}\right\rbrack \odot \underset{\mathit{\mathbf{I} } }{\underbrace{\left\lbrack \begin{array}{ccc} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{array}\right\rbrack } } } } \left\lbrack \begin{array}{c} \delta_1 \\ \delta_2 \\ \delta_3 \end{array}\right\rbrack

where \odot is the Hadamard product, i.e., element-wise multiplication. In Matlab this is written as A .* B.

E1=E2=E3=EE_1 =E_2 =E_3 =E, A1=A2=A3=AA_1 =A_2 =A_3 =A, L1=L3=LL_1 =L_3 =L and L2=2LL_2 =\sqrt{2}L.

Leading to

Se=E  AL[1000120001]{\Rightarrow \mathit{\mathbf{S} } }^e =\dfrac{E\;A}{L}\left\lbrack \begin{array}{ccc} 1 & 0 & 0\\ 0 & \dfrac{1}{\sqrt{2} } & 0\\ 0 & 0 & 1 \end{array}\right\rbrack

Step 7: Put everything together

f=ATN=ATSeδ=ATSeAS  u\mathit{\mathbf{f} }={\mathit{\mathbf{A} } }^T \mathit{\mathbf{N} }={\mathit{\mathbf{A} } }^T {\mathit{\mathbf{S} } }^e \boldsymbol \delta =\underset{\mathit{\mathbf{S} } }{\underbrace{ {\mathit{\mathbf{A} } }^T {\mathit{\mathbf{S} } }^e \mathit{\mathbf{A} } } } \;\mathit{\mathbf{u} } f=S  u\Rightarrow \mathit{\mathbf{f} }=\mathit{\mathbf{S} }\;\mathit{\mathbf{u} } S=E  AL[11200121][1000120001][10121201]=E  AL[1+1221221221+122]\begin{array}{l} S=\dfrac{E\;A}{L}\left\lbrack \begin{array}{ccc} 1 & \dfrac{1}{\sqrt{2} } & 0\\ 0 & -\dfrac{1}{\sqrt{2} } & -1 \end{array}\right\rbrack \left\lbrack \begin{array}{ccc} 1 & 0 & 0\\ 0 & \dfrac{1}{\sqrt{2} } & 0\\ 0 & 0 & 1 \end{array}\right\rbrack \left\lbrack \begin{array}{cc} 1 & 0\\ \dfrac{1}{\sqrt{2} } & -\dfrac{1}{\sqrt{2} }\\ 0 & -1 \end{array}\right\rbrack \\ =\dfrac{E\;A}{L}\left\lbrack \begin{array}{cc} 1+\dfrac{1}{2\sqrt{2} } & -\dfrac{1}{2\sqrt{2} }\\ -\dfrac{1}{2\sqrt{2} } & 1+\dfrac{1}{2\sqrt{2} } \end{array}\right\rbrack \end{array}

Step 8: Apply known values on f\mathit{\mathbf{f} } (and/ or u\mathit{\mathbf{u} })

f=[0P]=P[01]\mathit{\mathbf{f} }=\left\lbrack \begin{array}{c} 0\\ P \end{array}\right\rbrack =P\left\lbrack \begin{array}{c} 0\\ 1 \end{array}\right\rbrack

solve the system: u=S1f={MATLABS\f}=P  LE  A[22(2+2)(2232)]P  LE  A[0.20710.7929]\mathit{\mathbf{u} }={\mathit{\mathbf{S} } }^{-1} \mathit{\mathbf{f} }=\left\lbrace \underset{\mathrm{S}\backslash \mathrm{f} }{\textrm{MATLAB} } \right\rbrace =\frac{P\;L}{E\;A}\left\lbrack \begin{array}{c} \frac{\sqrt{2} }{2\left(\sqrt{2}+2\right)}\\ -\left(\frac{\sqrt{2} }{2}-\frac{3}{2}\right) \end{array}\right\rbrack \approx \frac{P\;L}{E\;A}\left\lbrack \begin{array}{c} 0\ldotp 2071\\ 0\ldotp 7929 \end{array}\right\rbrack

Step 9: The rod forces, strains, deformations and stresses are computed using u\mathit{\mathbf{u} }!

N=SeA  uP[0.20710.29290.7929]N={\mathit{\mathbf{S} } }^e \mathit{\mathbf{A} }\;\mathit{\mathbf{u} }\approx P\left\lbrack \begin{array}{c} 0\ldotp 2071\\ -0\ldotp 2929\\ -0\ldotp 7929 \end{array}\right\rbrack

Note! In cases when AT{\mathit{\mathbf{A} } }^T is square (static determinate structure) we can get the rod forces directly from the loads: N=(AT)1f\mathit{\mathbf{N} }={\left({\mathit{\mathbf{A} } }^T \right)}^{-1} \mathit{\mathbf{f} }.

δ=A  u\boldsymbol \delta =\mathit{\mathbf{A} }\;\mathit{\mathbf{u} } σi=NiAi\sigma_i =\dfrac{N_i }{A_i }

Example 2

The boundary condition at node 1 means u1x=u1y=0u_1^x =u_1^y =0 and at node two, we have the boundary condition u2y=0u_2^y =0.

Solution

We start by posing equilibrium equations on nodes which are partly free, i.e., nodes 2, 3 and 4.

Node 2:

:f2x=N1+1/2  N5:f2y=N21/2  N5\begin{array}{l} \to :f_2^x =N_1 +{1/\sqrt{2}\;N}_5 \\ \uparrow :f_2^y ={-N}_2 -{1/\sqrt{2}\;N}_5 \end{array}

Node 3:

:f3x=N3:f3y=N2\begin{array}{l} \to :f_3^x =N_3 \\ \uparrow :f_3^y =N_2 \end{array}

Node 4:

:f4x=N31/2N5:f4y=N4+1/2N5\begin{array}{l} \to :f_4^x ={-N}_3 -1/\sqrt{2}N_5 \\ \uparrow :f_4^y =N_4 +1/\sqrt{2}N_5 \end{array}

Note how we are skipping node 1, since all degrees of freedom are prescribed to zero. We could add it though, and deal with eliminating equations later.

Then we can combine the node equations

[f2xf2yf3xf3yf4xf4y]=[10001/201001/2001000100000101/200011/2][N1N2N3N4N5]\left\lbrack \begin{array}{c} f_2^x \\ f_2^y \\ f_3^x \\ f_3^y \\ f_4^x \\ f_4^y \end{array}\right\rbrack =\left\lbrack \begin{array}{ccccc} 1 & 0 & 0 & 0 & 1/\sqrt{2}\\ 0 & -1 & 0 & 0 & -1/\sqrt{2}\\ 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & -1/\sqrt{2}\\ 0 & 0 & 0 & 1 & 1/\sqrt{2} \end{array}\right\rbrack \left\lbrack \begin{array}{c} N_1 \\ N_2 \\ N_3 \\ N_4 \\ N_5 \end{array}\right\rbrack

or

f=ATN\mathit{\mathbf{f} }={\mathit{\mathbf{A} } }^T \mathit{\mathbf{N} }

from which we identify AT{\mathit{\mathbf{A} } }^T.

Material relations

Here we formulate the lengths of the bars, each entry along the diagonal corresponds to a bar.

Se=E  A[1L000001L000001  L000001  L0000012L]S^e =E\;A\left\lbrack \begin{array}{ccccc} \dfrac{1}{L} & 0 & 0 & 0 & 0\\ 0 & \dfrac{1}{L} & 0 & 0 & 0\\ 0 & 0 & \dfrac{1}{\;L} & 0 & 0\\ 0 & 0 & 0 & \dfrac{1}{\;L} & 0\\ 0 & 0 & 0 & 0 & \dfrac{1}{\sqrt{2}L} \end{array}\right\rbrack

or

Se=E  A  [1L,1L,1L,1L,12L]I{\mathit{\mathbf{S} } }^e =E\;A\;\left\lbrack \dfrac{1}{L},\dfrac{1}{L},\dfrac{1}{L},\dfrac{1}{L},\dfrac{1}{\sqrt{2}L}\right\rbrack \odot \mathit{\mathbf{I} }

Stiffness matrix

S=ATSeA\mathit{\mathbf{S} }={\mathit{\mathbf{A} } }^T {\mathit{\mathbf{S} } }^e \mathit{\mathbf{A} }

Now we will implement this using Matlab since it's kindof cumbersome to solve a 6×66\times 6 system by hand.

syms E A L P
sq = 1/sqrt(2);
AT = [1   0   0    0   sq
      0  -1   0    0  -sq
      0   0   1    0    0
      0   1   0    0    0
      0   0  -1    0   -sq
      0   0   0    1    sq];
Se = E*A*[1/L 1/L 1/L 1/L sq*1/L] .* eye(5);
S = AT*Se*AT'
S=(σ2(σ1)σ100σ1σ1σ1σ2(σ1)0σ2σ1σ100AEL0σ200σ20AEL00σ1σ1σ20σ2(σ1)σ1σ1σ100σ1σ2(σ1))where    σ1=2AE4L    σ2=AEL \begin{array}{l} \bm S = \left(\begin{array}{cccccc} -\sigma_2 -{\left(-\sigma_1 \right)} & -\sigma_1 & 0 & 0 & -\sigma_1 & \sigma_1 \\ -\sigma_1 & -\sigma_2 -{\left(-\sigma_1 \right)} & 0 & \sigma_2 & \sigma_1 & -\sigma_1 \\ 0 & 0 & \frac{A\,\textrm{E} }{L} & 0 & \sigma_2 & 0\\ 0 & \sigma_2 & 0 & \frac{A\,\textrm{E} }{L} & 0 & 0\\ -\sigma_1 & \sigma_1 & \sigma_2 & 0 & -\sigma_2 -{\left(-\sigma_1 \right)} & -\sigma_1 \\ \sigma_1 & -\sigma_1 & 0 & 0 & -\sigma_1 & -\sigma_2 -{\left(-\sigma_1 \right)} \end{array}\right)\\ \mathrm{}\\ \textrm{where}\\ \mathrm{}\\ \;\;\sigma_1 =\frac{\sqrt{2}\,A\,\textrm{E} }{4\,L}\\ \mathrm{}\\ \;\;\sigma_2 =-\frac{A\,\textrm{E} }{L} \end{array}

Now we insert what we know into the system.

[S11S12S13S14S15S16S21S31S41S51S61S66][u2xu2y=0u3xu3yu4xu4y]=[f2x=0f2y=0f3x=0.7Pf3y=0.3Pf4x=3/5Pf4y=4/5P] \def\arraystretch{1.5} \left\lbrack \begin{array}{cccccc} S_{11} & S_{12} & S_{13} & S_{14} & S_{15} & S_{16} \\ S_{21} & & & & & \\ S_{31} & & \ddots & & & \vdots \\ S_{41} & & & & & \\ S_{51} & & & & & \\ S_{61} & & \ldots & & & S_{66} \end{array}\right\rbrack \left\lbrack \begin{array}{c} u_2^x \\ u_2^y =0\\ u_3^x \\ u_3^y \\ u_4^x \\ u_4^y \end{array}\right\rbrack =\left\lbrack \begin{array}{c} f_2^x =0\\ f_2^y =0\\ f_3^x =0\ldotp 7P\\ f_3^y =-0\ldotp 3P\\ f_4^x =3/5P\\ f_4^y =4/5P \end{array}\right\rbrack

We have to reduce the system, by removing equations where the solution uu is known. This means, deleting row 2 and column 2.

f = [0
     0
     0.7*P
     0.3*P
     3/5*P
     4/5*P];
u = sym(zeros(6,1));
presc = 2;
free = [1,3,4,5,6];
u(free) = S(free,free)\f(free)
(13LP10AE0L(41P+262P)10AE3LP10AE2LP(172+26)10AE21LP10AE) \def\arraystretch{1.5} \left(\begin{array}{c} \frac{13\,L\,P}{10\,A\,\textrm{E} }\\ 0\\ \frac{L\,{\left(41\,P+26\,\sqrt{2}\,P\right)} }{10\,A\,\textrm{E} }\\ \frac{3\,L\,P}{10\,A\,\textrm{E} }\\ \frac{\sqrt{2}\,L\,P\,{\left(17\,\sqrt{2}+26\right)} }{10\,A\,\textrm{E} }\\ \frac{21\,L\,P}{10\,A\,\textrm{E} } \end{array}\right)

Visualization

In order to visualize the displacements, we need numeric data, in this case, we were not provided with data, so we set everything to 1.

E = 1; A = 1; L = 1; P = 1;
u = double(subs(u))
u = 6x1    
    1.3000
         0
    7.7770
    0.3000
    7.0770
    2.1000

We split the x and y components of the displacements into a matrix with the first column being the x-components and the second column being the y-components.

U = [0 0
     u(1:2:end), u(2:2:end)]
U = 4x2    
         0         0
    1.3000         0
    7.7770    0.3000
    7.0770    2.1000

Now, we shall create the topology of the structure, first the coordinates of our structure:

p = [0 0
     L 0
     L L
     0 L];

Then the edge connectivity.

edges = [1 2
         2 3
         3 4
         1 4
         2 4];

Now we can send this off to the patch function, which draws the whole structure in one simple line.

figure
patch('Faces',edges,'Vertices',p,'LineWidth',2)
axis equal
hold on

Now we shall overlay the undeformed structure with our deformed structure. We do this by adding the displacement field to our coordinate field such that

Pdeformed=Pundeformed+U  sP_{\textrm{deformed} } =P_{\textrm{undeformed} } +U\;s

where we need to scale the displacement field such that it shows the trend without being overly deformed.

scale = 0.05;
patch('Faces',edges,'Vertices',p+U*scale,'LineWidth',1,'LineStyle',':')

Lastly, we add node numbering and edge numbering.

plot(p(:,1),p(:,2),'ok','MarkerFaceColor','c','MarkerSize',14)
text(p(:,1)-0.01,p(:,2),cellstr(num2str([1:4]')))
for i=1:5
    inod = edges(i,:);
    xm = mean(p(inod,1));
    ym = mean(p(inod,2));
    text(xm,ym,num2str(i),'backgroundColor','y')
end

Prescribed displacement

We start by converting previously symbolic variables to numeric.

S = double(subs(S));
f = double(subs(f));
[S11S12S13S14S15S16S21S31S41S51S61S66][u2xu2y=15P  LE  Au3xu3yu4xu4y]=[f2x=0f2y=0f3x=0.7Pf3y=0.3Pf4x=3/5Pf4y=4/5P]\left\lbrack \begin{array}{cccccc} S_{11} & S_{12} & S_{13} & S_{14} & S_{15} & S_{16} \\ S_{21} & & & & & \\ S_{31} & & \ddots & & & \vdots \\ S_{41} & & & & & \\ S_{51} & & & & & \\ S_{61} & & \ldots & & & S_{66} \end{array}\right\rbrack \left\lbrack \begin{array}{c} u_2^x \\ u_2^y =-\dfrac{1}{5}\dfrac{P\;L}{E\;A}\\ u_3^x \\ u_3^y \\ u_4^x \\ u_4^y \end{array}\right\rbrack =\left\lbrack \begin{array}{c} f_2^x =0\\ f_2^y =0\\ f_3^x =0\ldotp 7P\\ f_3^y =-0\ldotp 3P\\ f_4^x =3/5P\\ f_4^y =4/5P \end{array}\right\rbrack

In this case, the displacement u2yu_2^y is not zero, and causes a reaction force which affects the whole structure. We deal with this by modifying the right hand side. In other words, we need to subtract the reaction forces which this displacement causes from the load vector. The reaction force is

FR=u2y[S21S22S23S24S25S26]F_R =u_2^y \left\lbrack \begin{array}{c} S_{21} \\ S_{22} \\ S_{23} \\ S_{24} \\ S_{25} \\ S_{26} \end{array}\right\rbrack

We end up with the reduced system of equations.

[S11S13S14S15S16S31S41S51S61S66][u2xu3xu3yu4xu4y]=[f2x=0f3x=0.7Pf3y=0.3Pf4x=3/5Pf4y=4/5P]u2y[S21S23S24S25S26]\left\lbrack \begin{array}{ccccc} S_{11} & S_{13} & S_{14} & S_{15} & S_{16} \\ S_{31} & \ddots & & & \vdots \\ S_{41} & & & & \\ S_{51} & & & & \\ S_{61} & \ldots & & & S_{66} \end{array}\right\rbrack \left\lbrack \begin{array}{c} u_2^x \\ u_3^x \\ u_3^y \\ u_4^x \\ u_4^y \end{array}\right\rbrack =\left\lbrack \begin{array}{c} f_2^x =0\\ f_3^x =0\ldotp 7P\\ f_3^y =-0\ldotp 3P\\ f_4^x =3/5P\\ f_4^y =4/5P \end{array}\right\rbrack -u_2^y \left\lbrack \begin{array}{c} S_{21} \\ S_{23} \\ S_{24} \\ S_{25} \\ S_{26} \end{array}\right\rbrack

We do this in Matlab by defining a variable for the prescribed degree of freedom.

presc = 2;

Then we modify the right-hand side

u = zeros(6,1);
u(presc) = -1/5*P*L/(E*A)
u = 6x1    
    0
   -0.2000
    0
    0
    0
    0
fr = S(:,presc)*u(presc)
fr = 6x1    
    0.0707
   -0.2707
         0
    0.2000
   -0.0707
    0.0707
f = f - fr
f = 6x1    
   -0.0707
    0.2707
    0.7000
    0.1000
    0.6707
    0.7293
free = [1,3,4,5,6]; %Free degrees of freedom
S(free,free)
ans = 5x5    
    1.3536         0         0   -0.3536    0.3536
         0    1.0000         0   -1.0000         0
         0         0    1.0000         0         0
   -0.3536   -1.0000         0    1.3536   -0.3536
    0.3536         0         0   -0.3536    1.3536
f(free)
ans = 5x1    
   -0.0707
    0.7000
    0.1000
    0.6707
    0.7293

Now we solve the for the free DOFs.

u(free) = S(free,free)\f(free)
u = 6x1    
    1.3000
   -0.2000
    7.9770
    0.1000
    7.2770
    2.1000

Now, let's visualize!

U = [0 0
     u(1:2:end), u(2:2:end)]
U = 4x2    
         0         0
    1.3000   -0.2000
    7.9770    0.1000
    7.2770    2.1000
figure
patch('Faces',edges,'Vertices',p,'LineWidth',2)
axis equal
hold on
scale = 0.05;
patch('Faces',edges,'Vertices',p+U*scale,'LineWidth',1,'LineStyle',':')
plot(p(:,1),p(:,2),'ok','MarkerFaceColor','c','MarkerSize',14)
text(p(:,1)-0.01,p(:,2),cellstr(num2str([1:4]')))
for i=1:5
    inod = edges(i,:);
    xm = mean(p(inod,1));
    ym = mean(p(inod,2));
    text(xm,ym,num2str(i),'backgroundColor','y')
end

We can see a rather small displacement in the negative y-direction in node 2. Comparing the numerical values shows that all displacements are displaced slightly downward.