Code
import sympy as sp # Symbolical math tools
import numpy as np # Numerical math tools
import matplotlib.pyplot as plt # Visualization tools
Determine how all of the forces in the system varies as the direction of the 4000 N force changes.
This is an example of the minimal amount of “paperwork” necessary before you start working with Python.
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.
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!
Imports
import sympy as sp # Symbolical math tools
import numpy as np # Numerical math tools
import matplotlib.pyplot as plt # Visualization tools
Constants
= 4 # [m] Distance ED
a = 3 # [m] Distance AE
b = 4000 # [N] Force magnitude: B
F_B = 200 # [N] Force magnitude: C F_C
Parameters
= sp.symbols('theta', real=True) # [rad] Angle at B theta
Variables
= sp.symbols('R_Ax, R_Ay, R_Ex', real=True) # [N] Unknowns R_Ax, R_Ay, R_Ex
Force vectors
= sp.Matrix([R_Ax, R_Ay, 0]) # [N] Reaction force A
RR_A = F_B * sp.Matrix([sp.cos(theta), sp.sin(theta), 0])# [N] External force B
FF_B = sp.Matrix([R_Ex, 0, 0]) # [N] Reaction force E
RR_E = F_C * sp.Matrix([0, 1, 0]) # [N] External force C FF_C
Positional vectors
= a * sp.Matrix([1, 0, 0]) # [m] Vector from A-B
rr_AB = sp.Matrix([2*a, -b, 0]) # [m] Vector from A-C
rr_AC = b * sp.Matrix([0, -1, 0]) # [m] Vector from A-E rr_AE
Newton’s second law
= sp.Eq(RR_A + FF_B + RR_E + FF_C,
N2 0, 0, 0]))
sp.Matrix([ N2
Euler’s second law
= sp.Eq(rr_AB.cross(FF_B) + rr_AC.cross(FF_C) + rr_AE.cross(RR_E),
E2 0, 0, 0]))
sp.Matrix([ E2
Solve for unknown variables
= sp.solve([N2, E2], [R_Ax, R_Ay, R_Ex]) sol
Symbolic -> Numerical
= sp.lambdify([theta], sol[R_Ax])
R_Ax_f = sp.lambdify([theta], sol[R_Ay])
R_Ay_f = sp.lambdify([theta], sol[R_Ex]) R_Ex_f
Visualization
# figure settings
'dark_background')
plt.style.use(= plt.figure(figsize=(8.00, 6.00))
fig = fig.add_subplot()
ax set(xlabel='Angle θ [°]', ylabel='Reaction force [N]')
ax.True, which='both', color='w', linestyle='-', alpha=0.1)
plt.grid(
plt.minorticks_on()
# draw values
= np.linspace(0, 2*np.pi, 100)
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')
ax.plot(thetas
# draw legend
'$R_{Ax}$', '$R_{Ay}$', '$R_{Ex}$']) ax.legend([
The beam is supported at A with rotational joints and a wire at B. The beam has a mass of
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.
Constants
### Input
= 9.81 # [m/s^2]
g = 475 # [kg] m
We model the placement of B with a symbolic parameter
The y-coordinate of D can be calculated with
= sp.symbols("r_OBx", real=True, positive=True)
r_OBx
= sp.Matrix([0.12, 0, 0])
rr_OA = sp.Matrix([2.5, 0, 0])
rr_OG = sp.Matrix([3.5, -0.25, 0])
rr_OC = sp.Matrix([r_OBx, 0.25, 0])
rr_OB
= np.tan(25*np.pi/180)*5
r_ODy = sp.Matrix([0, r_ODy, 0]) rr_OD
The direciton for the wire and thus the wire force,
and
# Direction of the wire force
= rr_OD - rr_OB
rr_BD = rr_BD.normalized() ee_BD
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 the force vectors
= sp.symbols("R_Ax, R_Ay, F_B", real=True)
R_Ax, R_Ay, F_B
= sp.Matrix([R_Ax, R_Ay, 0])
RR_A = m*g*sp.Matrix([0, -1, 0])
WW = sp.Matrix([0, -10000, 0])
FF_C = F_B * ee_BD FF_B
We then establish the equations of equilibrium in the following way:
### Establish the equations of equlibrium and solve the system
= sp.Eq(RR_A+WW+FF_C+FF_B,sp.Matrix([0,0,0]))
sumF = 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]))
sumM
= sp.solve([sumF,sumM], [R_Ax, R_Ay, F_B]) sol
Finally, we plot the reaction forces as functions of
### Plotting the reaction forces
# Create a plot and make settings
'seaborn-v0_8-whitegrid')
plt.style.use(= plt.figure(figsize=(7, 7))
fig = fig.add_subplot()
ax set(xlim=[0, 5], ylim=[-20000, 75000], xlabel='Placement of B [m]', ylabel='Force [N]')
ax."Reaction forces as a function to the wire placement")
ax.set_title(True, which='both', color='k', linestyle='-', alpha=0.1)
plt.grid(
plt.minorticks_on()
# Create the numerical data and functions
= np.linspace(0,5,100)
r_OBx_np = sp.lambdify([r_OBx], sol[F_B], 'numpy')
F_B_f = sp.lambdify([r_OBx], sol[R_Ax], 'numpy')
R_Ax_f = sp.lambdify([r_OBx], sol[R_Ay], 'numpy')
R_Ay_f
# Creating each plot
='-', color='black', label="$F_{B}$")
ax.plot(r_OBx_np, F_B_f(r_OBx_np), linestyle='-', color='blue', label="$R_{Ax}$")
ax.plot(r_OBx_np, R_Ax_f(r_OBx_np), linestyle='-', color='red', label="$R_{Ay}$")
ax.plot(r_OBx_np, R_Ay_f(r_OBx_np), linestyle
ax.legend()
It is interesting to note what happens to the reaction forces when the placement is close to the point A. When the
And we can also plot the geometry and force vectors in order to judge if the directions are reasonable given the calculated values for
### Plotting the vectors and points
# Create a plot and make settings
#plt.style.use('seaborn-v0_8-whitegrid')
= plt.figure(figsize=(6, 6))
fig2 = fig2.add_subplot()
ax2 set(xlim=[-1, 6], ylim=[-2, 5], xlabel='X [m]', ylabel='Y [m]')
ax2."Visualization of the beam and forces")
ax2.set_title(True, which='both', color='k', linestyle='-', alpha=0.1)
plt.grid(
plt.minorticks_on()
= rr_OB.subs(r_OBx, 5).evalf()
rr_OB_np = rr_OD.evalf()
rr_OD_np
# Plotting the beam
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)
ax2.plot(
# Plotting the points and labels
= np.array([0, 0, 0])
rr_O = 0.05
offset
*zip(rr_O), 'ok', ms=5)
ax2.plot(0]-5*offset,rr_O[1]+offset,"O")
ax2.text(rr_O[
*rr_OA[:-1], 'ok', ms=5)
ax2.plot(0]+offset, rr_OA[1]+offset,"A")
ax2.text(rr_OA[
*rr_OB_np[:-1], 'ok', ms=5)
ax2.plot(0]+offset, rr_OB_np[1]+offset,"B")
ax2.text(rr_OB_np[
*rr_OG[:-1], 'ok', ms=5)
ax2.plot(0]+offset, rr_OG[1]+offset,"G")
ax2.text(rr_OG[
*rr_OC[:-1], 'ok', ms=5)
ax2.plot(0]+offset, rr_OC[1]+offset,"C")
ax2.text(rr_OC[
*rr_OD[:-1], 'ok', ms=5)
ax2.plot(0]+offset, rr_OD[1]+offset,"D")
ax2.text(rr_OD[
# Create vectors with numeric values
= np.array(rr_OA).astype(np.float64)
rr_OA_np = np.array(rr_OC).astype(np.float64)
rr_OC_np = np.array(rr_OG).astype(np.float64)
rr_OG_np = np.array(WW).astype(np.float64)
WW_np = np.array(FF_C).astype(np.float64)
FF_C_np = np.array(FF_B.subs({r_OBx:5,F_B:F_B_f(5)})).astype(np.float64)
FF_B_np = np.array(RR_A.subs({R_Ax:R_Ax_f(5), R_Ay:R_Ay_f(5)})).astype(np.float64)
RR_A_np
# Plotting the force vectors
*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.quiver(
ax2.legend()
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
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.
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.
Constants
### Input
= 65 * sp.pi/180 #[rad]
theta = 80 #[kg]
m = 9.82 #[m/s^2] g
Based on the free body diagrams and the decisions for what points to calculate moments around, we see the need for the following vectors:
The first elements in
Because the structure is symmetric the first element in
### Define position vectors and vectors for moment calculations
= sp.Matrix([480/sp.tan(theta), 480, 0])
rr_OD = sp.Matrix([905/sp.tan(theta), 905, 0])
rr_OG = sp.Matrix([-480/sp.tan(theta), 480, 0])
rr_CE = sp.Matrix([-905/sp.tan(theta), 905, 0]) rr_CG
To define
and
= rr_OG - (rr_CG-rr_CE)
rr_OE = rr_OE - rr_OD rr_DE
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
= sp.symbols(
R_Ox, R_Oy, R_Dx, R_Dy, R_Gx, R_Gy, R_Ex, R_Ey, R_Cy 'R_Ox, R_Oy, R_Dx, R_Dy, R_Gx, R_Gy, R_Ex, R_Ey, R_Cy', real=True)
= sp.Matrix([R_Ox, R_Oy, 0])
RR_O = sp.Matrix([R_Dx, R_Dy, 0])
RR_D = sp.Matrix([R_Gx, R_Gy, 0])
RR_G = sp.Matrix([R_Ex, R_Ey, 0])
RR_E = sp.Matrix([0, R_Cy, 0])
RR_C = sp.Matrix([0, -m*g, 0]) FF
We then establish the equations of equilibrium in the following way:
Body 1 (ODG)
Body 2 (CEG)
Body 3 (DE)
### Establish the equations of equlibrium and solve the system
# Equilibrium for body 1
= sp.Eq(RR_O + RR_D + FF + RR_G, sp.Matrix([0,0,0]))
sumF1 = sp.Eq(rr_OD.cross(RR_D) + rr_OD.cross(FF) + rr_OG.cross(RR_G), sp.Matrix([0,0,0]))
sumM1
# Equilibrium for body 2
= sp.Eq(RR_C + RR_E + (-RR_G), sp.Matrix([0,0,0]))
sumF2 = sp.Eq(rr_CE.cross(RR_E) + rr_CG.cross(-RR_G), sp.Matrix([0,0,0]))
sumM2
# Equilibrium for body 3
= sp.Eq((-RR_D) + (-RR_E), sp.Matrix([0,0,0]))
sumF3 = sp.Eq(rr_DE.cross(-RR_E), sp.Matrix([0,0,0]))
sumM3
# Solve the system of equations
= sp.solve([sumF1,sumM1,sumF2,sumM2,sumF3,sumM3],
sol [R_Ox, R_Oy, R_Dx, R_Dy, R_Gx, R_Gy, R_Ex, R_Ey, R_Cy])
Our expectation based on intuition is that the reaction forces from the floor should be pointing up. That implies that