05.4 - Kinematic Analysis
Movement Mechanism
Inspiration
Reuleaux Triangle Animation [1]
Concept & Mobility Calculations
Position & Velocity Calculations
The location of Point D was then calculated by translating the location of Point C by length l in the x-direction. Vx and Vy for Point C were calculated by taking the difference of the past x-position/y-position from the current x-position/y-position.
Force Calculation
# Force Calculations
# Motor Torque: 15 kgf·cm
# Need lbf·in
# 1 kilogram-force centimeter = 0.8679616621 in*lbf
# Motor Torque: 15*0.8679616621 lbf·in
# Force = Torque/radius
# Radius = 2.5 in
# Motor Force: 15*0.8679616621/2.5 lbf
F = 15*0.8679616621/2.5
print(F)
# 5.2077699725999995 lbf
The force could not be animated due to the code running endlessly. (The same unedited previous code also did this). However, the animation would look like the force was always constant in the direction of triangle motion.
No acceleration analysis was conducted as the input speed would be constant and the acceleration of Points C and D would not affect the motion of the mechanism.
Plots
Animation
Punching Mechanism
Inspiration
Velocity Analysis Practice Problems - Problem 2
Position & Velocity Calculations
Force Calculations
Plots
Animation
Velocity Animation
Forces Animation
Movement Mechanism Code
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
class Visualize():
# Constructor
def __init__(self):
self.l = 2.5
self.lk = 2
self.ls = 2
self.x1 = -1.25
self.x2 = 1.25
self.y1 = -1.25
self.y2 = -1.25
self.ylim = 6
self.xlim = 6
ds = list(range(361))
self.xAs = []
self.xBs = []
self.ys = []
for d in ds:
if d < 90:
xA = 0
xB = self.l
y = (self.l/90)*d
elif d < 180:
xA = (-self.l/90)*(d-90)
xB = (-self.l/90)*(d-90)+self.l
y = self.l
elif d < 270:
xA = -self.l
xB = 0
y = ((-self.l/90)*(d-180))+self.l
else:
xA = ((self.l/90)*(d-270))-self.l
xB = ((self.l/90)*(d-270))
y = 0
self.xAs.append(xA)
self.xBs.append(xB)
self.ys.append(y)
self.yCs = []
self.xCs = []
for d in ds:
xA = self.xAs[d]
yA = self.ys[d]
lT = np.sqrt(((xA-self.x1)**2)+((yA-self.y1)**2))
C = np.arccos(((lT**2)+(self.ls**2)-(self.lk**2))/(2*lT*self.ls))
base = abs(xA-self.x1)
if base == 0:
al = np.deg2rad(90)
else:
al = np.arctan(abs(yA-self.y1)/base)
if xA < self.x1:
angle = np.deg2rad(180) - al - C
else:
angle = al - C
yC = self.y1 + self.ls*np.sin(angle)
xC = self.x1 + self.ls*np.cos(angle)
tx = np.sqrt(((yC-yA)**2)+((xC-xA)**2))
self.xCs.append(xC)
self.yCs.append(yC)
self.vxAs = []
self.vyAs = []
self.vxCs = []
self.vyCs = []
for d in ds:
vxA = self.xAs[d]-self.xAs[d-1]
self.vxAs.append(vxA)
vyA = self.ys[d]-self.ys[d-1]
self.vyAs.append(vyA)
vxC = self.xCs[d]-self.xCs[d-1]
self.vxCs.append(vxC)
vyC = self.yCs[d]-self.yCs[d-1]
self.vyCs.append(vyC)
plt.figure()
plt.grid()
plt.plot(ds, self.xAs)
plt.plot(ds, self.xBs)
plt.plot(ds, self.xCs)
plt.plot(ds, [x + self.l for x in self.xCs])
plt.title('Input Degree Versus X-Position')
plt.xlabel('Input Degree (°)')
plt.ylabel('X-Position (in)')
plt.legend(['Point A', 'Point B', 'Point C', 'Point D'])
plt.savefig('X-Postions Plot.png')
plt.figure()
plt.grid()
plt.plot(ds, self.ys)
plt.plot(ds, self.yCs)
plt.title('Input Degree Versus Y-Position')
plt.xlabel('Input Degree (°)')
plt.ylabel('Y-Position (in)')
plt.legend(['Points A & B', 'Points C & D'])
plt.savefig('Y-Postions Plot.png')
plt.figure()
plt.grid()
plt.plot(ds, self.vxAs)
plt.plot(ds, self.vxCs)
plt.title('Input Degree Versus X-Velocity')
plt.xlabel('Input Degree (°)')
plt.ylabel('X-Velocity (in/°)')
plt.legend(['Points A & B', 'Points C & D'])
plt.savefig('X-Velocities Plot.png')
plt.figure()
plt.grid()
plt.plot(ds, self.vyAs)
plt.plot(ds, self.vyCs)
plt.title('Input Degree Versus Y-Velocity')
plt.xlabel('Input Degree (°)')
plt.ylabel('Y-Velocity (in/°)')
plt.legend(['Points A & B', 'Points C & D'])
plt.savefig('Y-Velocities Plot.png')
self.fig, self.ax = plt.subplots(1,1,figsize=(6,6)) # New figure
self.ax.set_aspect('equal')
self.fig.set_tight_layout(True)
self.ax.set_xlim([-self.xlim,self.xlim])
self.ax.set_ylim([-self.ylim,self.ylim])
# Centers the axes
self.ax.spines['left'].set_position('center')
self.ax.spines['bottom'].set_position('center')
self.ax.spines['right'].set_color('none')
self.ax.spines['top'].set_color('none')
def plot(self,d):
d = int(d)
self.ax.clear()
self.ax.set_aspect('equal')
self.fig.set_tight_layout(True)
self.ax.set_xlim([-self.xlim,self.xlim])
self.ax.set_ylim([-self.ylim,self.ylim])
# Centers the axes
self.ax.spines['left'].set_position('center')
self.ax.spines['bottom'].set_position('center')
self.ax.spines['right'].set_color('none')
self.ax.spines['top'].set_color('none')
self.ax.plot([self.xAs[d],self.xBs[d]],[self.ys[d],self.ys[d]],color='purple', linestyle='solid')
self.ax.plot([self.x1,self.x2],[self.y1,self.y2],color='black', linestyle='solid')
self.ax.plot([self.xCs[d],self.xCs[d]+self.l],[self.yCs[d],self.yCs[d]],color='magenta', linestyle='solid')
self.ax.plot([self.xAs[d],self.xCs[d]],[self.ys[d],self.yCs[d]],color='pink', linestyle='solid')
self.ax.plot([self.xBs[d],self.xCs[d]+self.l],[self.ys[d],self.yCs[d]],color='blue', linestyle='solid')
self.ax.plot([self.xCs[d],self.x1],[self.yCs[d],self.y1],color='orange', linestyle='solid')
self.ax.plot([self.xCs[d]+self.l,self.x2],[self.yCs[d],self.y2],color='red', linestyle='solid')
self.ax.plot([self.xAs[d],self.xAs[d]],[self.ys[d],self.ys[d]+self.l],color='green', linestyle='solid')
self.ax.plot([self.xBs[d],self.xBs[d]],[self.ys[d],self.ys[d]+self.l],color='green', linestyle='solid')
self.ax.plot([self.xAs[d],self.xBs[d]],[self.ys[d]+self.l,self.ys[d]+self.l],color='green', linestyle='solid')
self.ax.arrow(self.xAs[d],self.ys[d],50*self.vxAs[d],50*self.vyAs[d],head_width=0.1,width=0.001,length_includes_head=True,linestyle='dashed',color='brown',overhang=1.0)
self.ax.arrow(self.xBs[d],self.ys[d],50*self.vxAs[d],50*self.vyAs[d],head_width=0.1,width=0.001,length_includes_head=True,linestyle='dashed',color='brown',overhang=1.0)
self.ax.arrow(self.xCs[d],self.yCs[d],50*self.vxCs[d],50*self.vyCs[d],head_width=0.1,width=0.001,length_includes_head=True,linestyle='dashed',color='cyan',overhang=1.0)
self.ax.arrow(self.xCs[d]+self.l,self.yCs[d],50*self.vxCs[d],50*self.vyCs[d],head_width=0.1,width=0.001,length_includes_head=True,linestyle='dashed',color='cyan',overhang=1.0)
self.ax.arrow(self.xAs[d],self.ys[d]+self.l,50*self.vxAs[d],50*self.vyAs[d],head_width=0.1,width=0.001,length_includes_head=True,linestyle='dashed',color='olive',overhang=1.0)
self.ax.arrow(self.xBs[d],self.ys[d]+self.l,50*self.vxAs[d],50*self.vyAs[d],head_width=0.1,width=0.001,length_includes_head=True,linestyle='dashed',color='olive',overhang=1.0)
self.ax.scatter([self.xAs[d]],[self.ys[d]], s=200, marker='.', color='black')
self.ax.scatter([self.xBs[d]],[self.ys[d]], s=200, marker='.', color='black')
self.ax.scatter([self.xCs[d]],[self.yCs[d]], s=200, marker='.', color='black')
self.ax.scatter([self.xCs[d]+self.l],[self.yCs[d]], s=200, marker='.', color='black')
self.ax.scatter([self.x1],[self.y1], s=200, marker='.', color='black')
self.ax.scatter([self.x2],[self.y2], s=200, marker='.', color='black')
self.ax.scatter([self.xAs[d]],[self.ys[d]+self.l], s=200, marker='.', color='black')
self.ax.scatter([self.xBs[d]],[self.ys[d]+self.l], s=200, marker='.', color='black')
self.ax.set_aspect('equal')
self.fig.set_tight_layout(True)
self.ax.set_xlim([-self.xlim,self.xlim])
self.ax.set_ylim([-self.ylim,self.ylim])
# Centers the axes
self.ax.spines['left'].set_position('center')
self.ax.spines['bottom'].set_position('center')
self.ax.spines['right'].set_color('none')
self.ax.spines['top'].set_color('none')
# Function to animate the crank
def animate(self):
# Runs on the first animation frame
def init():
self.ax.clear()
self.ax.set_aspect('equal')
self.fig.set_tight_layout(True)
self.ax.set_xlim([-self.xlim,self.xlim])
self.ax.set_ylim([-self.ylim,self.ylim])
# Centers the axes
self.ax.spines['left'].set_position('center')
self.ax.spines['bottom'].set_position('center')
self.ax.spines['right'].set_color('none')
self.ax.spines['top'].set_color('none')
# Adds the four bar plotting function to the animation
def getFrames(frames):
self.plot(frames)
# Animation function
delay_interval = 0.25
theta_frames = np.linspace(0,360,361)
anim = FuncAnimation(self.fig, getFrames, init_func=init, frames=theta_frames[::1], interval=delay_interval,repeat_delay=0.25)
# Save as GIF
anim.save('Triangle.gif', dpi=200, writer='pillow')
if __name__ == '__main__':
mechanism = Visualize()
mechanism.animate()
Punching Mechanism Analysis Code
"""
Team 5
Six-Bar Position & Velocity Analysis
(Four-Bar + Slider-Crank with offset)
"""
import numpy as np
import matplotlib.pyplot as plt
# offset of 0.5
# First all the helper functions
def fourbarpos(l, th_1, th_2, delta):
th_1_rad = th_1 * np.pi / 180
th_2_rad = th_2 * np.pi / 180
# Find my K constants for link lengths
K1 = l[0] / l[1]
K2 = l[0] / l[3]
K3 = (l[1] ** 2 - l[2] ** 2 + l[3] ** 2 + l[0] ** 2) / (2 * l[1] * l[3])
K4 = l[0] / l[2]
K5 = (l[3] ** 2 - l[0] ** 2 - l[1] ** 2 - l[2] ** 2) / (2 * l[1] * l[2])
# Quadratic constants with input angle th_2
A = -K1 - K2 * np.cos(th_2_rad) + K3 + np.cos(th_2_rad)
B = -2 * np.sin(th_2_rad)
C = K1 - K2 * np.cos(th_2_rad) + K3 - np.cos(th_2_rad)
D = np.cos(th_2_rad) - K1 + K4 * (np.cos(th_2_rad)) + K5
E = -2 * np.sin(th_2_rad)
F = K1 + (K4 - 1) * np.cos(th_2_rad) + K5
# Check for imaginary roots. We won't handle those here. -- should be eliminated with grashof condition but this can help check
disc_4 = B**2 - 4 * A * C
disc_3 = E**2 - 4 * D * F
if disc_4 < 0 or disc_3 < 0:
print("This function does not handle imaginary roots")
return
th_4_rad = 2 * np.arctan2((-B + delta * np.sqrt(disc_4)), (2 * A))
th_3_rad = 2 * np.arctan2((-E + delta * np.sqrt(disc_3)), (2 * D))
th_vec_rad = [th_1_rad, th_2_rad, th_3_rad, th_4_rad]
th_vec = [(x * 180 / np.pi) % 360 for x in th_vec_rad]
return th_vec
def fourbarvel(l, th_vec, w_2):
th_vec = [x * np.pi / 180 for x in th_vec]
w_3 = (
w_2
* l[1]
/ l[2]
* (np.sin(th_vec[3] - th_vec[1]))
/ (np.sin(th_vec[2] - th_vec[3]))
)
w_4 = (
w_2
* l[1]
/ l[3]
* (np.sin(th_vec[1] - th_vec[2]))
/ (np.sin(th_vec[3] - th_vec[2]))
)
VA = np.multiply([-np.sin(th_vec[1]), np.cos(th_vec[1])], (l[1] * w_2))
VBA = np.multiply([-np.sin(th_vec[2]), np.cos(th_vec[2])], (l[2] * w_3))
VB = np.multiply([-np.sin(th_vec[3]), np.cos(th_vec[3])], (l[3] * w_4))
w_vec_out = [0, w_2, w_3, w_4]
return w_vec_out, VA, VBA, VB
def slidercrankpos(l, th_1, th_2, th_4, delta):
th_1 = th_1 * np.pi / 180
th_2 = th_2 * np.pi / 180
th_4 = th_4 * np.pi / 180
if delta == -1:
# Obtain theta3 for open chain
th_3 = np.arcsin(-(l[0] * np.sin(th_2) - l[2]) / l[1]) + np.pi
else:
# Obtain theta3 for cross chain
th_3 = np.arcsin((l[0] * np.sin(th_2) - l[2]) / l[1])
# Compute l1
l1 = l[0] * np.cos(th_2) - l[1] * np.cos(th_3)
th_vec = [th_1, th_2, th_3, th_4]
return l1, [x * 180 / np.pi for x in th_vec]
def slidercrankvel(l_vec, th_vec, w_2):
th_vec = [x * np.pi / 180 for x in th_vec]
# w_3 for either configuration
w_3 = (l_vec[1]) / l_vec[2] * (np.cos(th_vec[1]) / np.cos(th_vec[2])) * w_2
# Compute l1_dot
l1_dot = -l_vec[1] * np.sin(th_vec[1]) * w_2 + l_vec[2] * w_3 * np.sin(th_vec[2])
return l1_dot, w_3
def compute_slider_force_MA(
theta_2_vec_local,
th_vec_out,
th_vec_slide,
l2,
l4,
l6,
T_in,
alpha_global=90,
units="SI",
):
"""
Computes mechanical advantage and output slider force using
transmission angle formulas.
Returns:
MA_total : array
Total mechanical advantage vs input angle
F_out : array
Output slider force vs input angle
"""
# --- Unit conversion ---
IN_TO_M = 0.0254
if units == "SI":
l2_m = l2 * IN_TO_M
l4_m = l4 * IN_TO_M
l6_m = l6 * IN_TO_M
else:
l2_m = l2
l4_m = l4
l6_m = l6
# --- Extract angles (convert to radians) ---
theta2 = np.deg2rad(theta_2_vec_local)
theta3 = np.deg2rad(th_vec_out[2, :])
# theta4 = np.deg2rad(th_vec_out[3, :]) not working
theta4_local = np.deg2rad(th_vec_out[3, :])
theta4_global = np.deg2rad((th_vec_out[3, :] - alpha_global) % 360)
theta5 = np.deg2rad(th_vec_slide[2, :])
theta6 = np.deg2rad(90.0) # slider constraint
mu1 = np.abs(theta4_local - theta3)
nu1 = np.abs(theta2 - theta3)
nu2 = np.abs(theta4_global - theta5)
mu2 = np.abs(theta5 - theta6)
# --- Mechanical advantage (lecture formulas) including clipped singularities which pop up at extreme angles ---
mu_min = np.deg2rad(10) # 5–10 deg is standard
nu_min = np.deg2rad(10)
mu1_c = np.maximum(mu1, mu_min)
nu1_c = np.maximum(nu1, nu_min)
mu2_c = np.maximum(mu2, mu_min)
nu2_c = np.maximum(nu2, nu_min)
MA_1 = (l4_m * np.sin(mu1_c)) / (l2_m * np.sin(nu1_c))
MA_2 = (l6_m * np.sin(mu2_c)) / (l4_m * np.sin(nu2_c))
MA_total = MA_1 * MA_2
F_out = MA_total * T_in / l2_m
return MA_total, F_out
# %% Problem Code starts here
# Link lengths and input angle
l1 = 2.00 # [inches]
l2 = 0.80 # [inches]
l3 = 1.60 # [inches]
l4 = 2.30 # [inches]
l5 = 4 # [inches]
l6 = 0.689 # [inches] no offset
alpha_global = 90 # deg
w_2 = -30 # rad/sec
# Just some set up of variables
theta_2_vec = np.linspace(0, 360, 1000)
theta_2_vec_local = (theta_2_vec + 102) % 360
l_vec = [l1, l2, l3, l4]
delta1 = -1 # open circuit
delta2 = -1 # open circuit
# Setting up arrays
l1_slide = np.zeros([1, len(theta_2_vec)])
l4_arr = np.full([1, len(theta_2_vec)], l4)
l_vec_slide = np.concatenate(
(
l1_slide,
np.full([1, len(theta_2_vec)], l4),
np.full([1, len(theta_2_vec)], l5),
np.full([1, len(theta_2_vec)], l6),
),
axis=0,
)
th_vec_out = np.zeros([4, len(theta_2_vec)])
th_vec_slide = np.zeros([4, len(theta_2_vec)])
w_vec_out = np.zeros([4, len(theta_2_vec)])
VA = np.zeros([2, len(theta_2_vec)])
VBA = np.zeros([2, len(theta_2_vec)])
VB = np.zeros([2, len(theta_2_vec)])
l_1_dot_b = np.zeros(len(theta_2_vec))
w_5 = np.zeros(len(theta_2_vec))
VC_out = np.zeros(len(theta_2_vec))
# %% Calculations start here
for thIdx in range(0, len(theta_2_vec)):
# Original fourbar position analysis
th_vec_out[:, thIdx] = fourbarpos(l_vec, 0, theta_2_vec_local[thIdx], delta1)
# Put the important outputs into global terms
th_vec_out_global = [
0,
(th_vec_out[3, thIdx] - alpha_global) % 360,
0,
90,
] # th_5 is dummy
# Slider block analysis
l_vec_slide[0, thIdx], th_vec_slide[:, thIdx] = slidercrankpos(
l_vec_slide[1:, thIdx],
th_vec_out_global[0],
th_vec_out_global[1],
th_vec_out_global[3],
delta2,
)
# Fourbar velocity analysis
w_vec_out[:, thIdx], VA[:, thIdx], VBA[:, thIdx], VB[:, thIdx] = fourbarvel(
l_vec, th_vec_out[:, thIdx], w_2
)
# Slider block velocity analysis
l_1_dot_b[thIdx], w_5[thIdx] = slidercrankvel(
l_vec_slide[:, thIdx], th_vec_slide[:, thIdx], w_vec_out[3, thIdx]
)
VC_out[thIdx] = l_1_dot_b[thIdx]
T_in = 0.01 # N*m (small hobby motor) STILL A LARGE ESTIMATE NEAR STALL TORQUE!!
MA_total, F_out = compute_slider_force_MA(
theta_2_vec_local, th_vec_out, th_vec_slide, l2, l4, l6, T_in
)
# %% Plotting
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec, (th_vec_out[2, :] - alpha_global) % 360)
ax.set_title(r"$\theta_3$ vs input angle")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$\theta_3$ (deg)")
ax.grid()
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec, (th_vec_out[3, :] - alpha_global) % 360)
ax.set_title(r"$\theta_4$ vs input angle")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$\theta_4$ (deg)")
ax.grid()
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec - alpha_global, (th_vec_out[3, :]) % 360)
ax.set_title(r"$\theta_4$ vs input angle (local)")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$\theta_4$ (deg)")
ax.grid()
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec, (th_vec_slide[2, :]) % 360)
ax.set_title(r"$\theta_5$ vs input angle (global)")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$\theta_5$ (deg)")
ax.grid()
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec, l_vec_slide[0, :])
ax.set_title(r"$l_{O_{4}O_{6}}$ vs input angle (global)")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$l_{O_{4}O_{6}}$")
ax.grid()
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec, w_vec_out[2, :])
ax.set_title(r"$\omega_3$ vs input angle")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$\omega_3$ (rad/sec)")
ax.grid()
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec, w_vec_out[3, :])
ax.set_title(r"$\omega_4$ vs input angle ($\theta_{2_{global}}=0$))")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$\omega_4$ (rad/sec)")
ax.grid()
th2idx = np.where(abs(theta_2_vec - 258) < 0.05)[0][0]
w_vec_plot = np.concatenate([w_vec_out[3, th2idx:], w_vec_out[3, 0:th2idx]])
fig, ax = plt.subplots()
plt.subplots_adjust(hspace=0.6)
ax.plot(theta_2_vec, w_vec_plot)
ax.set_title(r"$\omega_4$ vs input angle ($\theta_{2_{local}}=0$))")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$\omega_4$ (rad/sec)")
ax.grid()
fig, ax = plt.subplots()
ax.plot(theta_2_vec, VC_out)
ax.set_title(r"Slider Output Velocity $V_S$ vs Input Angle $\theta_2$")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel(r"$V_S$ (in/s)")
ax.grid()
fig, ax = plt.subplots()
ax.plot(theta_2_vec, MA_total)
ax.set_title("Total Mechanical Advantage vs Input Angle")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel("Mechanical Advantage")
ax.grid()
fig, ax = plt.subplots()
ax.plot(theta_2_vec, F_out)
ax.set_title("Output Slider Force vs Input Angle")
ax.set_xlabel(r"$\theta_2$ (deg)")
ax.set_ylabel("Slider Force (N)")
ax.grid()
plt.show()