12  Equilibrium Examples

12.1 Example 1

12.1.1 Problem description

Determine how all of the forces in the system varies as the direction of the 4000 N force changes.

12.1.2 Paperwork

This is an example of the minimal amount of “paperwork” necessary before you start working with Python.

12.1.3 Free body diagram and vector definitions

A free body diagram consists of several definitions: a coordinate system, points of interest, forces, moments and distances. Additionally, clearly define any known constants or other parameters that are allowed to vary. Finally, define all force- and positional vectors in terms of both knowns and unknowns.

12.1.4 Analysis of unknowns and equations of motion

Compare the number of unknown variables to the number of equations you can establish for the given problem. If they are the same, you can proceed to coding. If not, more thinking and/or paperwork is required!

12.1.5 Code

Imports

Code
import sympy as sp                                       # Symbolical math tools
import numpy as np                                       # Numerical math tools
import matplotlib.pyplot as plt                          # Visualization tools

Constants

Code
a = 4                                                    # [m] Distance ED
b = 3                                                    # [m] Distance AE
F_B = 4000                                               # [N] Force magnitude: B
F_C = 200                                                # [N] Force magnitude: C

Parameters

Code
theta = sp.symbols('theta', real=True)                   # [rad] Angle at B

Variables

Code
R_Ax, R_Ay, R_Ex = sp.symbols('R_Ax, R_Ay, R_Ex', real=True) # [N] Unknowns

Force vectors

Code
RR_A = sp.Matrix([R_Ax, R_Ay, 0])                        # [N] Reaction force A
FF_B = F_B * sp.Matrix([sp.cos(theta), sp.sin(theta), 0])# [N] External force B
RR_E = sp.Matrix([R_Ex, 0, 0])                           # [N] Reaction force E
FF_C = F_C * sp.Matrix([0, 1, 0])                        # [N] External force C

Positional vectors

Code
rr_AB = a * sp.Matrix([1, 0, 0])                         # [m] Vector from A-B
rr_AC = sp.Matrix([2*a, -b, 0])                          # [m] Vector from A-C
rr_AE = b * sp.Matrix([0, -1, 0])                        # [m] Vector from A-E

Newton’s second law

Code
N2 = sp.Eq(RR_A + FF_B + RR_E + FF_C, 
           sp.Matrix([0, 0, 0]))
N2

[RAx+REx+4000cos(θ)RAy+4000sin(θ)+2000]=[000]

Euler’s second law

Code
E2 = sp.Eq(rr_AB.cross(FF_B) + rr_AC.cross(FF_C) + rr_AE.cross(RR_E), 
           sp.Matrix([0, 0, 0]))
E2

[003REx+16000sin(θ)+1600]=[000]

Solve for unknown variables

Code
sol = sp.solve([N2, E2], [R_Ax, R_Ay, R_Ex])

Symbolic -> Numerical

Code
R_Ax_f = sp.lambdify([theta], sol[R_Ax])
R_Ay_f = sp.lambdify([theta], sol[R_Ay])
R_Ex_f = sp.lambdify([theta], sol[R_Ex])

Visualization

Code
# figure settings
plt.style.use('dark_background')
fig = plt.figure(figsize=(8.00, 6.00))
ax = fig.add_subplot()
ax.set(xlabel='Angle θ [°]', ylabel='Reaction force [N]')
plt.grid(True, which='both', color='w', linestyle='-', alpha=0.1)
plt.minorticks_on()

# draw values
thetas = np.linspace(0, 2*np.pi, 100)
ax.plot(thetas*180/np.pi, R_Ax_f(thetas), 'red')
ax.plot(thetas*180/np.pi, R_Ay_f(thetas), 'green')
ax.plot(thetas*180/np.pi, R_Ex_f(thetas), 'blue')

# draw legend
ax.legend(['$R_{Ax}$', '$R_{Ay}$', '$R_{Ex}$'])

12.2 Equilibrium: Example 1

12.2.1 Problem formulation

The beam is supported at A with rotational joints and a wire at B. The beam has a mass of m=475 kg. There is also a force acting on the hook as shown in the figure.

Figure of the problem.

12.2.2 Do the following:

  • Create a free body diagram of the beam shown in the figure, with all applied forces and reaction forces.
  • Define the force vectors, select a point to define moments around, and define the moments.
  • Calculate the reaction forces at A and the magnitude of the force in the wire connected at B.
  • Create a plot:
    • draw the origin, the beam, all points of interest with both markers and text.
    • plot the force vectors
  • Let the connection of the wire vary along the length of the beam between A and B and plot the magnitude of the wire force and reaction forces for all wire positions.

12.2.3 Solution steps

We start by selecting a point for our moments to be calculated around. In this case we select the center of the leftmost end of the beam and call that point O.

Then we create the necessary vectors between different points and how they relate to each other. We also name the points according the the figure below.

We can then define our force vectors and formulate our equations of equilibrium.

12.2.4 Paperwork

Free body diagram with names of all points and vectors.

Equilibrium equations.

12.2.5 Define the force vectors

Constants

### Input
g = 9.81 # [m/s^2]
m = 475 # [kg]

We model the placement of B with a symbolic parameter rOBx and then procede to create the vectors from O to the points of interest according to the figure.

rOB=[rOBx0.250]

The y-coordinate of D can be calculated with

rODy=5tan25

Code
r_OBx = sp.symbols("r_OBx", real=True, positive=True)

rr_OA = sp.Matrix([0.12, 0, 0])
rr_OG = sp.Matrix([2.5, 0, 0])
rr_OC = sp.Matrix([3.5, -0.25, 0])
rr_OB = sp.Matrix([r_OBx, 0.25, 0])

r_ODy = np.tan(25*np.pi/180)*5
rr_OD = sp.Matrix([0, r_ODy, 0])

The direciton for the wire and thus the wire force, FB, can be calculated in the following way: rBD=rODrOB

and

eBD=rBD||rBD||

Code
# Direction of the wire force
rr_BD = rr_OD - rr_OB
ee_BD = rr_BD.normalized()

When it comes to the forces, we define the force vectors using the definitions from the free body diagram created earlier. We create symbolic variables for the unknown scalar values RAx, RAy, and FB.

Code
### Define the force vectors

R_Ax, R_Ay, F_B = sp.symbols("R_Ax, R_Ay, F_B", real=True)

RR_A = sp.Matrix([R_Ax, R_Ay, 0])
WW = m*g*sp.Matrix([0, -1, 0])
FF_C = sp.Matrix([0, -10000, 0])
FF_B = F_B * ee_BD

We then establish the equations of equilibrium in the following way:

ΣF=0:RA+W+FC+FB=0

ΣM=0:(rOA×RA)+(rOG×W)+(rOC×FC)+(rOB×FB)=0

Code
### Establish the equations of equlibrium and solve the system

sumF = sp.Eq(RR_A+WW+FF_C+FF_B,sp.Matrix([0,0,0]))
sumM = sp.Eq(rr_OA.cross(RR_A)+rr_OG.cross(WW)+rr_OC.cross(FF_C)+rr_OB.cross(FF_B),sp.Matrix([0,0,0]))

sol = sp.solve([sumF,sumM], [R_Ax, R_Ay, F_B])
FB=5.61127562510170.230797548423172rOBx2+114001293545187.0rOBx1500000000000.0
RAx=1.078294000137911021rOBx5.60051741807481016rOBx6.01015
RAy=8.210218521962211017rOBx56005174180748.0rOBx6000000000000.02.33246875101856005174180748.0rOBx6000000000000.0

12.3 Plot the forces, points and bases

Finally, we plot the reaction forces as functions of rOBx.

Show the code
### Plotting the reaction forces

# Create a plot and make settings
plt.style.use('seaborn-v0_8-whitegrid')
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot()
ax.set(xlim=[0, 5], ylim=[-20000, 75000], xlabel='Placement of B [m]', ylabel='Force [N]')
ax.set_title("Reaction forces as a function to the wire placement")
plt.grid(True, which='both', color='k', linestyle='-', alpha=0.1)
plt.minorticks_on()

# Create the numerical data and functions
r_OBx_np = np.linspace(0,5,100)
F_B_f = sp.lambdify([r_OBx], sol[F_B], 'numpy')
R_Ax_f = sp.lambdify([r_OBx], sol[R_Ax], 'numpy')
R_Ay_f = sp.lambdify([r_OBx], sol[R_Ay], 'numpy')

# Creating each plot
ax.plot(r_OBx_np, F_B_f(r_OBx_np), linestyle='-', color='black', label="$F_{B}$")
ax.plot(r_OBx_np, R_Ax_f(r_OBx_np), linestyle='-', color='blue', label="$R_{Ax}$")
ax.plot(r_OBx_np, R_Ay_f(r_OBx_np), linestyle='-', color='red', label="$R_{Ay}$")


ax.legend()

Figure of the reaction forces as a function of the wire placement.

It is interesting to note what happens to the reaction forces when the placement is close to the point A. When the FB vector intercects A the reaction forces are infinite. Also FB is negative when FB intercects between O and A. That is of course not physical behaviour of a wire.

And we can also plot the geometry and force vectors in order to judge if the directions are reasonable given the calculated values for rOBx=5.

Show the code
### Plotting the vectors and points

# Create a plot and make settings
#plt.style.use('seaborn-v0_8-whitegrid')
fig2 = plt.figure(figsize=(6, 6))
ax2 = fig2.add_subplot()
ax2.set(xlim=[-1, 6], ylim=[-2, 5], xlabel='X [m]', ylabel='Y [m]')
ax2.set_title("Visualization of the beam and forces")
plt.grid(True, which='both', color='k', linestyle='-', alpha=0.1)
plt.minorticks_on()

rr_OB_np = rr_OB.subs(r_OBx, 5).evalf()
rr_OD_np = rr_OD.evalf()

# Plotting the beam
ax2.plot([0,0,5,5,0], [-0.25,0.25,0.25,-0.25,-0.25], linestyle='-', color='black')
ax2.plot(*zip(rr_OB.subs(r_OBx, 5)[:-1],rr_OD[:-1]), linestyle='--', color='black', zorder=0)


# Plotting the points and labels
rr_O = np.array([0, 0, 0])
offset = 0.05

ax2.plot(*zip(rr_O), 'ok', ms=5)
ax2.text(rr_O[0]-5*offset,rr_O[1]+offset,"O")

ax2.plot(*rr_OA[:-1], 'ok', ms=5)
ax2.text(rr_OA[0]+offset, rr_OA[1]+offset,"A")

ax2.plot(*rr_OB_np[:-1], 'ok', ms=5)
ax2.text(rr_OB_np[0]+offset, rr_OB_np[1]+offset,"B")

ax2.plot(*rr_OG[:-1], 'ok', ms=5)
ax2.text(rr_OG[0]+offset, rr_OG[1]+offset,"G")

ax2.plot(*rr_OC[:-1], 'ok', ms=5)
ax2.text(rr_OC[0]+offset, rr_OC[1]+offset,"C")

ax2.plot(*rr_OD[:-1], 'ok', ms=5)
ax2.text(rr_OD[0]+offset, rr_OD[1]+offset,"D")


# Create vectors with numeric values
rr_OA_np = np.array(rr_OA).astype(np.float64)
rr_OC_np = np.array(rr_OC).astype(np.float64)
rr_OG_np = np.array(rr_OG).astype(np.float64)
WW_np = np.array(WW).astype(np.float64)
FF_C_np = np.array(FF_C).astype(np.float64)
FF_B_np = np.array(FF_B.subs({r_OBx:5,F_B:F_B_f(5)})).astype(np.float64)
RR_A_np = np.array(RR_A.subs({R_Ax:R_Ax_f(5), R_Ay:R_Ay_f(5)})).astype(np.float64)

# Plotting the force vectors
ax2.quiver(*rr_OC_np[:-1], *FF_C_np[:-1], color='orange', scale=100000, label="$\\mathbb{F}_C$")
ax2.quiver(*rr_OB_np[:-1], *FF_B_np[:-1], color='red', scale=100000, label="$\\mathbb{F}_B$")
ax2.quiver(*rr_OA_np[:-1], *RR_A_np[:-1], color='red', scale=100000, label="$\\mathbb{R}_A$")
ax2.quiver(*rr_OG_np[:-1], *WW_np[:-1], color='blue', scale=50000, label="$\\mathbb{W}$")


ax2.legend()

Figure of the beam with the points and force vectors.

12.4 Equilibrium: Example 3

12.4.1 Problem formulation

The structure shown in the figure is supported at A, B and C from the floor. The force acting on the point D represents a human with mass m=80 kg.

Figure of the problem.

12.4.2 Do the following:

  • Create the necessary free body diagrams of the structure shown in the figure, with all applied forces and reaction forces.
  • Define the force vectors, select points to define moments around, and define the moments.
  • Calculate the reaction forces at A, B, and C, as well as the joint forces in D, E, and G.

12.4.3 Simplification

Before we begin the formal solution, we can reason about the problem and see if its possible to simplify it. In this case we note that the points D and G are placed in the center of the structure relative to A and B. Based on our intuition and previous knowledge of forces we see that it is reasonable to simplify this problem to 2D as the reaction force in A and B will be equal. The reaction force we would calculate for the point representing A and B will be the sum of the two.

2D-simplification of the problem.

12.4.4 Solution steps

Based of this, we start by determining we need to create free body diagrams of the bodies ODG, CEG, and DE.

We select points for our moments to be calculated around. In this case we select the point O in the body ODG, point C in CEG, and point D in DE.

Then we create the necessary vectors between different points and how they relate to each other.

We can then define our force vectors and formulate our equations of equilibrium.

12.4.5 Paperwork

Free body diagram with names of all points and vectors.

Position vector definitions.

Force vector definitions.

Equilibrium equations.

12.4.6 Define the vectors between points

Constants

### Input
theta = 65 * sp.pi/180 #[rad]
m = 80 #[kg]
g = 9.82 #[m/s^2]

Based on the free body diagrams and the decisions for what points to calculate moments around, we see the need for the following vectors: rOD, rOG, rCE, rCG, and rDE.

The first elements in rOD, rOG, rCE, and rCG can be calculated with for example:

rODx=480tan65

Because the structure is symmetric the first element in rCE will be, rCEx=480tan65

Code
### Define position vectors and vectors for moment calculations

rr_OD = sp.Matrix([480/sp.tan(theta), 480, 0])
rr_OG = sp.Matrix([905/sp.tan(theta), 905, 0])
rr_CE = sp.Matrix([-480/sp.tan(theta), 480, 0])
rr_CG = sp.Matrix([-905/sp.tan(theta), 905, 0])

To define rDE we can use the vectors already defined in the following way:

rOE=rOG(rCGrCE)

and

rDE=rOErOD

Code
rr_OE = rr_OG - (rr_CG-rr_CE)
rr_DE = rr_OE - rr_OD

12.4.7 Define the force vectors and equations of equilibrium

When it comes to the forces, we define the force vectors using the definitions from the free body diagram created earlier. We create symbolic variables for the unknown scalar values.

### Define force vectors

R_Ox, R_Oy, R_Dx, R_Dy, R_Gx, R_Gy, R_Ex, R_Ey, R_Cy = sp.symbols(
  'R_Ox, R_Oy, R_Dx, R_Dy, R_Gx, R_Gy, R_Ex, R_Ey, R_Cy', real=True)

RR_O = sp.Matrix([R_Ox, R_Oy, 0])
RR_D = sp.Matrix([R_Dx, R_Dy, 0])
RR_G = sp.Matrix([R_Gx, R_Gy, 0])
RR_E = sp.Matrix([R_Ex, R_Ey, 0])
RR_C = sp.Matrix([0, R_Cy, 0])
FF = sp.Matrix([0, -m*g, 0])

We then establish the equations of equilibrium in the following way:

Body 1 (ODG)

ΣF:RO+RD+F+RG=0

ΣMO:(rOD×RD)+(rOD×F)+(rOG×RG)=0

Body 2 (CEG)

ΣF:RC+RE+(RG)=0

ΣMC:(rCE×RE)+(rCG×(RG))=0

Body 3 (DE)

ΣF:(RD)+(RE)=0

ΣMD:rDE×(RE)=0

Code
### Establish the equations of equlibrium and solve the system

# Equilibrium for body 1
sumF1 = sp.Eq(RR_O + RR_D + FF + RR_G, sp.Matrix([0,0,0]))
sumM1 = sp.Eq(rr_OD.cross(RR_D) + rr_OD.cross(FF) + rr_OG.cross(RR_G), sp.Matrix([0,0,0]))

# Equilibrium for body 2
sumF2 = sp.Eq(RR_C + RR_E + (-RR_G), sp.Matrix([0,0,0]))
sumM2 = sp.Eq(rr_CE.cross(RR_E) + rr_CG.cross(-RR_G), sp.Matrix([0,0,0]))

# Equilibrium for body 3
sumF3 = sp.Eq((-RR_D) + (-RR_E), sp.Matrix([0,0,0]))
sumM3 = sp.Eq(rr_DE.cross(-RR_E), sp.Matrix([0,0,0]))

# Solve the system of equations
sol = sp.solve([sumF1,sumM1,sumF2,sumM2,sumF3,sumM3], 
               [R_Ox, R_Oy, R_Dx, R_Dy, R_Gx, R_Gy, R_Ex, R_Ey, R_Cy])
ROx=0.0
ROy=577.26408839779
RDx=206.869437880414
RDy=0.0
RGx=206.869437880414
RGy=208.33591160221
REx=206.869437880414
REy=0.0
RCy=208.33591160221

12.4.8 Conclusion

Our expectation based on intuition is that the reaction forces from the floor should be pointing up. That implies that ROy and RCy should be positive, which is what we observe. The larger part of the applied force F is taken up at O, since D is closer to O than to C in the X-direction. Since we have no external force in the X-direction, we expect ROx=0, which we also get in our calculations. From the negative value of REx and positive value of RDx we can conclude from free body diagram 3 that the body DE is in tension.