13.2 Kinematic Analysis

Kinematic Analysis:

Mobility: M = 3(8 − 1) − 2(10) − 0 = 1 DOF

Grashof condition:
shortest + longest ≤ sum of the other two links

Loop 1: m–k–c–a

m = 90.0 mm
k = 371.4 mm
c = 235.8 mm
a = 228.0 mm

Shortest = m = 90.0 mm
Longest = k = 371.4 mm

Check:

90.0 + 371.4 ≤ 235.8 + 228.0

461.4 ≤ 463.8

Result:
461.4 ≤ 463.8 → Grashof condition satisfied

Loop 2: m–j–b–a

m = 90.0 mm
j = 300.0 mm
b = 249.0 mm
a = 228.0 mm

Shortest = m = 90.0 mm
Longest = j = 300.0 mm

Check:

90.0 + 300.0 ≤ 249.0 + 228.0

390.0 ≤ 477.0

Result:
390.0 ≤ 477.0 → Grashof condition satisfied

Loop 3: d–c–g–f

d = 240.6 mm
c = 235.8 mm
g = 220.2 mm
f = 236.4 mm

Shortest = g = 220.2 mm
Longest = d = 240.6 mm

Check:

220.2 + 240.6 ≤ 235.8 + 236.4

460.8 ≤ 472.2

Result:
460.8 ≤ 472.2 → Grashof condition satisfied

Final summary:
All three four-bar loops satisfy the Grashof condition. Therefore, the mechanism allows continuous crank rotation and should not lock due to non-Grashof geometry.

Our kinematic analysis focuses on characterizing the time-dependent motion of the assembly by evaluating not only the position, but the trace path, velocity, and acceleration of the Toe Point. In this analysis, we make sure that we map out the visual trace path of the toe point from the position vector in order to make sure that the trajectory of the end avoids breaking within our workspace limits. In order to make sure that all of the actuators can actually drive the system, the velocity plots are analyzed in order to make sure the required speeds don’t exceed the motors' peak ratings. We also analyze the acceleration plots to check where exactly we need stronger structural supports to make sure the leg doesn’t vibrate a lot.

Position Plots

Toe Path Plot

Velocity Plots

Acceleration Plots

Below is the code that was used to generate all of the plots used above:

import numpy as np
import matplotlib.pyplot as plt
""" --- POSITION ANALYSIS --- """
# Initialize variables
# Actual CAD measurements of the linkages:
# b = 124.50
# c = 117.9
# d = 120.3
# e = 167.4
# f = 118.2
# g = 110.1
# h = 197.1
# i = 147
# j = 150
# k = 185.7
# m = 45
# x = 114
# y = 23
a = 228 # x
b = 249
c = 235.8
d = 240.6
e = 334.8
f = 236.4
h = 394.2
i = 294
j = 300
k = 372
l = 46.8 #y
m = 90
g = 110.1*2
# #Jansen Holy numbers:
# a = 38  # x-component of ground link
# y = 7.8   # y-component of ground link
# b = 41.5  # Upper 4-bar Rocker
# c = 39.3  # Lower 4-bar Rocker
# d = 40.1
# e = 55.8  # Rigid body "bde"
# f = 39.4  # Parallel-like mechanism
# g = 36.7 
# h = 65.7
# i = 49  # Rigid body "ghi"
# j = 50  # Upper 4-bar Coupler
# k = 61.9  # Lower 4-bar Coupler
# m = 15  # Crank
theta_m_degrees = np.arange(360)
theta_j = []
theta_b = []
theta_k = []
theta_c = []
theta_d = []
theta_f = []
theta_g = []
theta_i = []
theta_p = []
p = []
P_p_x = []
P_p_y = []
X3 = [] 
Y3 = []  
n = np.sqrt(a**2 + l**2)  # Ground link
"""
In this mechanism we refer to FIGURE 0 for the links and linkages of the Theo Jansen Mechanism.
The Jansen mechanism is composed of
- crank (m)
- 2 oscillating rockers (b,c)
- 2 couplers (j, k)
- 2 three-bar linkages (b,d,e) and (g,h,i)
- upper four bar linkages (nm-bj)
- lower four bar linkages (nm-ck)
- open four-bar linkage/parallel lik linkage (dc-fg)
- rigid three-bar (g,h,i)
- Ground (02-04) 
- Toe (P)
We let O2 be the origin of the General Coordinate System.
Note that in general four bar notiation we use FIGURE 1 as a reference for our analysis for
the four bar linkages.
For Figure 1, we can write the vector loop equation as 
d - c - a- b = 0
d exp(jtheta_d) - c exp(jtheta_c) - a exp(jtheta_a) - b exp(jtheta_b) = 0
Using these we can analyze the four linkages (nm-bj), (nm-ck).
"""
for angle in theta_m_degrees:
    # Compute the position of point B
    X_m = a + m * np.cos(np.deg2rad(angle)) #x comp of point B from origin
    Y_m = l + m * np.sin(np.deg2rad(angle)) #y comp of point B from origin
    R = np.sqrt(X_m**2 + Y_m**2) #vector of Point B from origin
    alpha = np.rad2deg(np.arctan2(Y_m, X_m)) #angle of point 1 wrt O2 ground/origin
    """
    ---- VECTOR LOOP (1) UPPER 4-BAR LINKAGE (nm-bj) ---
    Referring to FIGURE 2, we can write the equation of the components of the linkages as:
    n cos(theta_n) + m cos(theta_m) - b cos(theta_b) - j cos(theta_j) = 0
    n sin(theta_n) + m sin(theta_m) - b sin(theta_b) - j sin(theta_j) = 0
    Although a four bar linkage we use R and cos law to get theta_j and theta_b.
    This is done by making a triangle for links {b,j,R}.
    Let A_j be the projecton of j on R.
    Let A_b be the projection of b on R.
    See FIGURE 3.
    """
    A_j = (j**2 - b**2 + R**2) / (2 * j) #projecton of j on R
    A_b = (b**2 - j**2 + R**2) / (2 * b) #projection of b on R
    theta_j.append(np.rad2deg(np.arccos(A_j / R)) + alpha) #theta_j as function of m
    theta_b.append(np.rad2deg(np.arccos(A_b / R)) + alpha) #tehta_b as function of m
    """
    --- VECTOR LOOP (2) LOWER 4-BAR LINAKG (nm-ck) ---
    Again, we use R and cos law to get theta_k and theta_c.
    Note that we are using the same R here since it is the same points for both (nm-bj) and (nm-ck).
    This is done by making a triangle for links {c,k,R}.
    Let A_k be the projecton of k on R.
    Let A_c be the projection of c on R.
    See FIGURE 4 and FIGURE 5.
    """
    A_k = (k**2 - c**2 + R**2) / (2 * k) #projecton of k on R
    A_c = (c**2 - k**2 + R**2) / (2 * c) #projecton of c on R
    theta_k.append(np.rad2deg(np.arccos(A_k / R)) + alpha) #theta_k as function of m
    theta_c.append(np.rad2deg(np.arccos(A_c / R)) + alpha) #theta_c as function of m
    """
    --- PARALLEL LINKAGE (dc-fg) --
    Need to find theta_d since it is the common link between triangle bde and parallel dc-fg.
    Usining cosine law with bde FIGURE 6:
    """
    theta_d_value = np.degrees(np.arccos((b**2 + d**2 - e**2) / (2 * b * d))) + theta_b[-1]
    theta_d.append(theta_d_value)
    """
    Now looking at the parallelogram (FIGURE 7), we take theta d as the new "moving ground".
    We take point 4 as the input since the rotation from m is transmitted to that point.
    Let Point 5 be the output since that point helps determine point P.
    Just like in our analysis with the 2 four bar linkages where we used R to extend a line
    from the origin/input to the output point to help with the analysis, see FIGURE 8. 
    We use R3 to to extend a line from point 4 to point 5. 
    We let X3 and Y3 be the components of R3, and alpha3 be its angle from horizontal.
    """
    X3_value = d * np.cos(np.radians(theta_d_value)) + c * np.cos(np.radians(theta_c[-1]))
    Y3_value = d * np.sin(np.radians(theta_d_value)) + c * np.sin(np.radians(theta_c[-1]))
    X3.append(X3_value)
    Y3.append(Y3_value)
    R3 = np.sqrt(X3[-1]**2 + Y3[-1]**2)
    alpha3 = np.rad2deg(np.arctan(Y3_value/X3_value)) 
    """
    Now that we have alpha3 of point 5, we can use it to find theta_f and theta_g in the parallelogram.
    Just as before, we use projections to help us.
    Let A_f be the projection of f on R3.
    Let A_g be the projection of g on R3.
    We use R3 since this is the third four-bar we are analyzing, 
    and the first two four-bars used the same R.
    """
    A_f = (f**2 - g**2 + R3**2) / (2 * f) #projection of f on R3
    A_g = (g**2 - f**2 + R3**2) / (2 * g) #projection of g on R3
    A_f_clipped = np.clip(A_f / R3, -1, 1) #clipped to take away rouding errors 
    A_g_clipped = np.clip(A_g / R3, -1, 1) #clipped to take away rouding errors 
    theta_f_value = np.rad2deg(np.arccos(A_f_clipped)) + alpha3
    theta_g_value = np.rad2deg(np.arccos(A_g_clipped)) + alpha3
    theta_f.append(theta_f_value) 
    theta_g.append(theta_g_value)
    """
    --- (ghi) TRIANLGE + (cip) TRIANGLE ---
    Referring to Figure 9:
    Now we need find find point p as a function of the crank angle theta_m
    First we find theta i from cosine law of triangle {ghi}
    """
    theta_i_value = theta_g[-1] - np.rad2deg(np.arccos((g**2 + i**2 - h**2) / (2 * g * i)))
    theta_i.append(theta_i_value)
    """
    From the diagram, P is just the sum of the x and y components of link c and link i.
    Therefore, since we already knowt theta_c from the second four-bar analysis and we obtained
    theta i from the ghi triangle, we can find theta p.
    """
    # Angular Position
    theta_p_value = np.rad2deg(np.arctan(
    (c * np.sin(np.deg2rad(theta_c[-1])) + i * np.sin(np.deg2rad(theta_i[-1]))) /
    (c * np.cos(np.deg2rad(theta_c[-1])) + i * np.cos(np.deg2rad(theta_i[-1])))
    ))
    theta_p.append(theta_p_value)
    p_x = c * np.cos(np.deg2rad(theta_c[-1])) + i * np.cos(np.deg2rad(theta_i[-1]))
    p_y = c * np.sin(np.deg2rad(theta_c[-1])) + i * np.sin(np.deg2rad(theta_i[-1]))
    p_value = np.sqrt(p_x**2 + p_y**2)
    p.append(p_value)
    # Linear position
    P_p_y.append(p_y)  # y-component of P
    P_p_x.append(p_x)  # x-component of P
# Position Plot
plt.figure(1)
# Angular Position of Toe Point (P) of Walking Mechanism
plt.subplot(3, 1, 1)
plt.plot(theta_m_degrees, theta_p)  # Plot against theta_m_degrees
plt.title("Angular Position of Toe Point (P) of Walking Mechanism")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel(r'Angular Position (deg) or $\theta_{P}(°)$')
plt.grid(True)
# Position of Toe Point (P) of Walking Mechanism (x-direction)
plt.subplot(3, 1, 2)
plt.plot(theta_m_degrees, P_p_x)  # Plot against theta_m_degrees
plt.title("Position of Toe Point (P) of Walking Mechanism (x-direction)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Toe Point Position (mm)')
plt.grid(True)
# Position of Toe Point (P) of Walking Mechanism (y-direction)
plt.subplot(3, 1, 3)
plt.plot(theta_m_degrees, P_p_y)  # Plot against theta_m_degrees
plt.title("Position of Toe Point (P) of Walking Mechanism (y-direction)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Toe Point Position (mm)')
plt.grid(True)
# Save the figure
plt.tight_layout()  # Adjust subplots to fit into figure area.
# For toe path
plt.figure(2)
plt.plot(P_p_x, P_p_y)
plt.title("Trace Path of Toe Point (P) of Walking Mechanism")
plt.xlabel('Toe Point Position (x-direction)(mm)')
plt.ylabel('Toe Point Position (y-direction)(mm)')
plt.grid(True)
# Velocity and Acceleration Calculations
omega_m = 10.0
d_theta_m = np.deg2rad(1)
V_p_x = np.gradient(P_p_x, d_theta_m) * omega_m
V_p_y = np.gradient(P_p_y, d_theta_m) * omega_m
V_p_mag = np.sqrt(V_p_x**2 + V_p_y**2)
A_p_x = np.gradient(V_p_x, d_theta_m) * omega_m
A_p_y = np.gradient(V_p_y, d_theta_m) * omega_m
A_p_mag = np.sqrt(A_p_x**2 + A_p_y**2)
# Velocity Plots
plt.figure(3)
plt.subplot(3, 1, 1)
plt.plot(theta_m_degrees, V_p_x)
plt.title("Velocity of Toe Point (P) of Walking Mechanism (x-direction)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Velocity (mm/s)')
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(theta_m_degrees, V_p_y)
plt.title("Velocity of Toe Point (P) of Walking Mechanism (y-direction)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Velocity (mm/s)')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(theta_m_degrees, V_p_mag)
plt.title("Magnitude of Velocity of Toe Point (P)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Velocity (mm/s)')
plt.grid(True)
plt.tight_layout()
# Acceleration Plots
plt.figure(4)
plt.subplot(3, 1, 1)
plt.plot(theta_m_degrees, A_p_x)
plt.title("Acceleration of Toe Point (P) of Walking Mechanism (x-direction)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Acceleration ($mm/s^2$)')
plt.grid(True)
plt.subplot(3, 1, 2)
plt.plot(theta_m_degrees, A_p_y)
plt.title("Acceleration of Toe Point (P) of Walking Mechanism (y-direction)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Acceleration ($mm/s^2$)')
plt.grid(True)
plt.subplot(3, 1, 3)
plt.plot(theta_m_degrees, A_p_mag)
plt.title("Magnitude of Acceleration of Toe Point (P)")
plt.xlabel(r'Crank angle(deg) or $\theta_{m}(°)$')
plt.xlim([0, 360])
plt.ylabel('Acceleration ($mm/s^2$)')
plt.grid(True)
plt.tight_layout()
plt.show()

Torque Analysis:

In order to determine what kind of motor we needed to propel the skateboard, we analyzed the torque using an impulse of 50 N·s, representing the standard push of a skateboard. Using the dimensions already used to do the Kinematic Analysis, the script below uses circle-intersection algorithms in order to map both the trajectory and velocity of the foot across an entire 360-degree rotation of the crank.

import numpy as np

# All dimensions already used in the Kinematic Analysis
links = {
'a': 228, # x
'b': 249,
'c': 235.8,
'd': 240.6,
'e': 334.8,
'f': 236.4,
'h': 394.2,
'i': 294,
'j': 300,
'k': 372,
'l': 46.8, #y
'm': 90,
'g': 110.1*2
}

# Kinematic Geometry Solver - Calculating the coordinate or every joint from the crank to the foot.
def intersect(p1, r1, p2, r2, d_dir=1):
    d = np.linalg.norm(p2-p1)
    if d > r1+r2 or d < abs(r1-r2) or d==0: return None
    a = (r1**2 - r2**2 + d**2)/(2*d)
    h = np.sqrt(max(0, r1**2 - a**2))
    p3 = p1+a * (p2-p1)/d
    offset = np.array([-h*(p2[1] - p1[1])/d, h*(p2[0]-p1[0])/d])
    return p3+offset if d_dir == 1 else p3-offset

def get_foot(theta):
    hip = np.array([-links['l'], -links['m']])
    crank = np.array([links['a']*np.cos(theta), links['a']*np.sin(theta)])

    UR = intersect(hip, links['b'], crank, links['j'], 1)
    if UR is None: return None
    LR = intersect(hip, links['c'], crank, links['k'], -1)
    if LR is None: return None
    Knee = intersect(UR, links['e'], LR, links['d'], 1)
    if Knee is None: return None
    Ankle = intersect(LR, links['f'], Knee, links['i'], -1)
    if Ankle is None: return None
    return intersect(Knee, links['h'], Ankle, links['g'], 1)

torques = []
y = []
d_theta = 0.001

for angle in np.linspace(0, 360, 360):
    theta = np.radians(angle)
    f1 = get_foot(theta)
    f2 = get_foot(theta + d_theta)
    if f1 is not None and f2 is not None:
        vx = (f2[0] - f1[0])/d_theta
        torques.append(abs(vx))
        y.append(f1[1])

# Obtaining the peak torque
y_min = min(y)
ground_limit = y_min + (max(y)-y_min)*0.05
peak = max(t for t, z in zip(torques, y) if z<=ground_limit)
meterpeak = peak/1000
print(f"Peak Torque Multiplier: {meterpeak:.2f}")

The code was able to produce the peak normalized torque multiplier of 0.857 N·m of required crank torque per 1 Newton of a horizontal thrust. Since this multiplier is less than one, this means that the machine has a bit of a mechanical advantage in comparison to a 1:1 direct-drive system. Though this means that to have an actually meaningful thrust in the horizontal direction, we will need a very high-torque! To do this, we have two methods: either to use an extremely high-torque gear motor, or to pair a high-torque Brushless DC motor with a planetary gearbox in order to withstand the movement done by the machine.