Skip to article frontmatterSkip to article content

Lesson Monday September 9th

During today’s lesson it’s demonstrated how you to use differential equations to solve structural problems.

Demonstration

Given a structure as shown below:

We can calculate the force distribution and displacements of this structure using differential equations.

First, let’s define the segments with continuous loads. This requires splitting the top element DB\text{DB} into two because the distributed load causes a discontinuity. This leads to four segments:

  • AC\text{AC}
  • AD\text{AD}
  • DE\text{DE}
  • EB\text{EB}

Now, each segment can be drawn separately, indicating the loads, coordinate systems, dimensions and section forces at the edges of each segment.

Element AC\text{AC}

As element is hinged at both ends and is only loaded at its end, it’s a zero-force member. The free-body diagram therefore only includes normal forces:

Free body diagram element \text{AC}

Free body diagram element AC\text{AC}

As this element is unloaded, the loading equation for this element gives:

qx,AC=0q_{x,\text{AC}} = 0

leading to the equations:

NAC(x)=C1N_\text{AC}\left( x \right) = C_1
uAC(x)=C1x+C220000u_\text{AC}\left( x \right) = \cfrac{C_1 \cdot x + C_2}{20000}

Because of the support at C\text{C} the first boundary condition is as follows:

uAC(10)=0u_\text{AC} \left(10\right) = 0

At C\text{C} the normal force is unknown (equilibrium with support reactions).

At A\text{A} the extension is unknown as A\text{A} can move horizontally this leads to a potential displacement of A\text{A} in the longitudinal direction of the element.

The normal force at A\text{A} follows from the equilibrium of node A\text{A}. It’s free-body-diagram is:

Free body diagram node \text{A}

Figure 3:Free body diagram node A\text{A}

Horizontal equilibrium gives the second boundary condition:

FhA=045NAC+150=0\sum {{{\left. {{F_{\rm{h}}}} \right|}^{\rm{A}}} = 0} \to \frac{4}{5} \cdot {N_{{\rm{AC}}}} + 150 = 0

This element could be solved independently as there are two unknowns (two integration constants) and two equations ((4) and (5)).

Element AD\text{AD}

Again, this element is a zero-force member. The free-body diagram therefore only includes normal forces:

Free body diagram element \text{AD}

Free body diagram element AD\text{AD}

As this element is unloaded too, the load equation for this element gives:

qx,AD=0q_{x,\text{AD}}=0

leading to the equations:

NAD(x)=C3N_\text{AD}\left( x \right) = C_3
uAD(x)=C3x+C420000u_\text{AD}\left( x \right) = \cfrac{C_3 \cdot x + C_4}{20000}

Because of the support at A\text{A} the first boundary condition is as follows:

uAD(9)=0u_\text{AD} \left(9\right) = 0

At D\text{D} another boundary condition can be identified from the consistency of displacements: the downwards displacement of the top of element AD\text{AD} must be the same as the downwards displacement of the left of element DE\text{DE}:

uAD(0)=wDE(0){u_{{\rm{AD}}}}\left( 0 \right) = {w_{{\rm{DE}}}}\left( 0 \right)

And finally, a boundary condition can be formulated as the normal force in D\text{D} is in equilibrium with the forces from element DE\text{DE}. The free-body-diagram of D\text{D} is:

Free body diagram node \text{D}

Free body diagram node D\text{D}

Vertical equilibrium gives:

FvD=0NAD+VDDE=0\sum {{{\left. {{F_{\rm{v}}}} \right|}^{\rm{D}}} = 0} \to {N_{{\rm{AD}}}} + V_{\rm{D}}^{{\rm{DE}}} = 0

Element DE\text{DE}

The free-body diagram of this element is:

Free body diagram element \text{ED}

Free body diagram element ED\text{ED}

As this element is unloaded, the load equation for this element gives:

qz,DE=0q_{z,\text{DE}}=0

leading to the equations:

VDE(x)=C5V_\text{DE} \left(x\right) = C_5
MDE(x)=C5x+C6M_\text{DE} \left(x\right) = C_5 \cdot x + C_6
κDE(x)=C5x+C65000\kappa_\text{DE} \left(x\right) = \cfrac{C_5 \cdot x + C_6}{5000}
φDE(x)=12C5x2+C6x5000+C7\varphi_\text{DE} \left(x\right) = \cfrac{\cfrac{1}{2}\cdot C_5 \cdot x^2 + C_6 \cdot x}{5000} + C_7
wDE(x)=16C5x312C6x25000C7x+C8w_\text{DE} \left(x\right) = \cfrac{-\cfrac{1}{6}\cdot C_5 \cdot x^3 - \cfrac{1}{2} C_6 \cdot x^2}{5000} - C_7 \cdot x + C_8

Because of the hinge at D\text{D}, a boundary condition can be formulated as:

MDE(0)=0M_\text{DE} \left(0\right) = 0

At E\text{E} another boundary condition can be identified from the consistency of displacements: both the downwards and rotational displacement of the right of element DE\text{DE} must be the same as the downwards and rotational displacement of the left of element BE\text{BE}:

wDE(4)=wBE(0)w_\text{DE} \left(4\right) = w_\text{BE} \left(0\right)
φDE(4)=φEB(0)\varphi_\text{DE} \left(4\right) = \varphi_\text{EB} \left(0 \right)

Finally, there should be consistency with the section forces too, as shown in the free body diagram:

Free body diagram node \text{E}

Free body diagram node E\text{E}

Vertical equilibrium gives:

FvE=0VEDE+VEBE=0\sum {{{\left. {{F_{\rm{v}}}} \right|}^{\rm{E}}} = 0} \to -{V_{{\rm{E}}^\text{DE}}} + V_{\rm{E}}^{{\rm{BE}}} = 0

Moment equilibrium gives:

TEE=0MEDE+MEBE=0\sum {{{\left. {{T}} \right|}_\text{E}^{\rm{E}}} = 0} \to -{M_{{\rm{E}}^\text{DE}}} + M_{\rm{E}}^{{\rm{BE}}} = 0

Element BE\text{BE}

The free-body diagram of this element is:

Free body diagram element \text{BE}

Free body diagram element BE\text{BE}

As this element is loaded, the load equation for this element gives:

qz,BE=10q_{z,\text{BE}}=10

leading to the equations:

VBE(x)=10x+C9V_\text{BE} \left(x\right) = -10 \cdot x + C_9
MBE(x)=5x2+C9x+C10M_\text{BE} \left(x\right) = -5 \cdot x^2 + C_9 \cdot x + C_10
κBE(x)=5x2+C9x+C105000\kappa_\text{BE} \left(x\right) = \cfrac{-5 \cdot x^2 + C_9 \cdot x + C_10}{5000}
φBE(x)=53x3+12C9x2+C10x5000+C11\varphi_\text{BE} \left(x\right) = \cfrac{-\cfrac{5}{3} \cdot x^3 + \cfrac{1}{2}\cdot C_9 \cdot x^2 + C_10 \cdot x}{5000} + C_11
wBE(x)=512x416C9x312C10x25000C11x+C12w_\text{BE} \left(x\right) = \cfrac{\cfrac{5}{12} \cdot x^4 -\cfrac{1}{6}\cdot C_9 \cdot x^3 - \cfrac{1}{2} C_10 \cdot x^2}{5000} - C_11 \cdot x + C_12

Because of the support at D, two boundary conditions can be formulated directly:

φBE(4)=0\varphi_\text{BE} \left(4\right) = 0
wBE(4)=0w_\text{BE} \left(4\right) = 0

Solve system of equations

12 boundary conditions have been formulated with 12 integration constants. These can now be solved. You’re advised to use your calculator or a symbolic programming tool to solve this.

Filling in the boundary conditions gives:

C12000+C2=04C15+150=09C320000+C4=0C4+C8=0C3+C5=0C6=0C12+4C51875+C6625+4C7C8=0C11C5625C61250C7=0C5+C9=0C104C5C6=0C101250+C11+C96258375=0C106254C11+C124C91875+8375=0\begin{align} \frac{C_{1}}{2000} + C_{2} = 0 \\ \frac{4 C_{1}}{5} + 150 = 0 \\ \frac{9 C_{3}}{20000} + C_{4} = 0 \\ - C_{4} + C_{8} = 0 \\ C_{3} + C_{5} = 0 \\ C_{6} = 0 \\ C_{12} + \frac{4 C_{5}}{1875} + \frac{C_{6}}{625} + 4 C_{7} - C_{8} = 0 \\ C_{11} - \frac{C_{5}}{625} - \frac{C_{6}}{1250} - C_{7} = 0 \\ - C_{5} + C_{9} = 0 \\ C_{10} - 4 C_{5} - C_{6} = 0 \\ \frac{C_{10}}{1250} + C_{11} + \frac{C_{9}}{625} - \frac{8}{375} = 0 \\ - \frac{C_{10}}{625} - 4 C_{11} + C_{12} - \frac{4 C_{9}}{1875} + \frac{8}{375} = 0 \\ \end{align}

Solving this gives:

C1=3752C2=332C3=1792415C4=504259375C5=1792415C6=0C7=4904778125C8=504259375C9=1792415C10=7168415C11=472778125C12=2792155625\begin{align} C_{1} = - \frac{375}{2} \\ C_{2} = \frac{3}{32}\\ C_{3} = - \frac{1792}{415}\\ C_{4} = \frac{504}{259375}\\ C_{5} = \frac{1792}{415} \\ C_{6} = 0\\ C_{7} = - \frac{4904}{778125} \\ C_{8} = \frac{504}{259375}\\ C_{9} = \frac{1792}{415} \\ C_{10} = \frac{7168}{415}\\ C_{11} = \frac{472}{778125}\\ C_{12} = \frac{2792}{155625}\\ \end{align}

Now, the full force distribution and displacements of the structure has been solved. For example the moment distribution in BE\text{BE} is 5x2+1792415x+7168415- 5 \cdot x^{2} + \cfrac{1792}{415} \cdot x + \cfrac{7168}{415} which can be plotted as:

Loading...

Solve using digital tools

To solve the system of equations using numerical tools, you can write down the equations as a matrix formulation Ax=bAx=b which can easily be solved by ie. your graphical calculator or python.

[00120001000000000004500000000009200001000000000001035000000000101000001000001000100000000000001000000000004187516254100010000016251125010001000000100010000000041000100000000000162511250100][C1C2C3C4C5C6C7C8C10C11C12]=[01500000000008375]\left[\begin{array}{ccccccccccccc}0 & 0 & \frac{1}{2000} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \frac{4}{5} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\frac{9}{20000} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\1 & 0 & \frac{3}{5} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1\\0 & -1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \frac{4}{1875} & \frac{1}{625} & 4 & -1 & 0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & - \frac{1}{625} & - \frac{1}{1250} & -1 & 0 & 0 & 0 & 1 & 0 & 0\\0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & -4 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{625} & \frac{1}{1250} & 1 & 0 & 0\end{array}\right] \left[ \begin{array} {ccccccccccccc}C_1\\C_2\\C_3\\C_4\\C_5\\C_6\\C_7\\C_8\\C_{10}\\C_{11}\\C_{12}\end{array}\right] = \left[\begin{matrix}0\\-150\\0\\0\\0\\0\\0\\0\\0\\0\\0\\\frac{8}{375}\end{matrix}\right]
np.linalg.solve(A,b)
array([[-1.87500000e+02], [ 9.37500000e-02], [-4.31807229e+00], [ 1.94313253e-03], [ 4.31807229e+00], [ 0.00000000e+00], [-6.30232932e-03], [ 1.94313253e-03], [ 4.31807229e+00], [ 1.72722892e+01], [ 6.06586345e-04], [ 1.79405622e-02]])

Alternatively, a symbolic programming language can be used to solve the equations exactly. As an example the sympy package in Python is used:

import sympy as sym
sym.init_printing()

x = sym.symbols('x')

q = 10 
EI = 5000 
EA = 20000 
F = 150

q_AC = 0
q_AD = 0
q_CD = 0
q_DE = 0
q_EB = q

C_1, C_2 = sym.symbols('C_1 C_2')
C_3, C_4 = sym.symbols('C_3 C_4')
C_5, C_6, C_7, C_8 = sym.symbols('C_5 C_6 C_7 C_8')
C_9, C_10, C_11, C_12 = sym.symbols('C_9 C_10 C_11 C_12')

N_AC = -sym.integrate(q_AC, x) + C_1
eps_AC = N_AC/EA
u_AC = sym.integrate(eps_AC, x) + C_2

N_AD = -sym.integrate(q_AD, x) + C_3
eps_AD  = N_AD/EA
u_AD = sym.integrate(eps_AD, x) + C_4

V_DE = -sym.integrate(q_DE, x) + C_5
M_DE = sym.integrate(V_DE, x) + C_6
kappa_DE = M_DE/EI
phi_DE = sym.integrate(kappa_DE, x) + C_7
w_DE = -sym.integrate(phi_DE, x) + C_8

V_EB = -sym.integrate(q_EB, x) + C_9
M_EB = sym.integrate(V_EB, x) + C_10
kappa_EB = M_EB/EI
phi_EB = sym.integrate(kappa_EB, x) + C_11
w_EB = -sym.integrate(phi_EB, x) + C_12

eq1 = sym.Eq(u_AC.subs(x,10),0)
eq2 = sym.Eq(F + N_AC.subs(x,0)/5*4,0)
eq3 = sym.Eq(u_AD.subs(x,9),0)
eq4 = sym.Eq(w_DE.subs(x,0)-u_AD.subs(x,0),0)
eq5 = sym.Eq(V_DE.subs(x,0) + N_AD.subs(x,0),0)
eq6 = sym.Eq(M_DE.subs(x,0),0)
eq7 = sym.Eq(w_EB.subs(x,0)-w_DE.subs(x,4),0)
eq8 = sym.Eq(phi_EB.subs(x,0)-phi_DE.subs(x,4),0)
eq9 = sym.Eq(V_EB.subs(x,0) - V_DE.subs(x,4),0)
eq10 = sym.Eq(M_EB.subs(x,0) - M_DE.subs(x,4),0)
eq11 = sym.Eq(phi_EB.subs(x,4),0)
eq12 = sym.Eq(w_EB.subs(x,4),0)

sol = sym.solve([eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9,eq10,eq11,eq12],(C_1,C_2,C_3,C_4,C_5,C_6,C_7,C_8,C_9,C_10,C_11,C_12))
display(sol)
Loading...