Worm Robot Calculation and Simulations

Here are a list of formulas to create theoretical expectations and to process experimentation data, and the subsequent calculations and simulations.

TODO:

  1. make a measurement block
  2. simulate the movement of the worm
In [1]:
# Resources
%matplotlib inline
import numpy as np
import scipy as sp
import pandas as pd
import seaborn as sns
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
sns.set(style="white")

Mesh manufacturing

In [2]:
# Simulation functions
def mesh_plot_3d (len_mesh,rad_mesh,beta,resolution):
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    theta = np.linspace(0, 3 * np.pi, len_mesh*resolution)
    x = np.linspace(0, len_mesh, len_mesh*resolution)
    r = rad_mesh * np.ones(len_mesh*resolution)
    y1 = r * np.sin(theta)
    y2 = - r * np.sin(theta)
    z = r * np.cos(theta)
    ax.plot(x, y1, z,sns.xkcd_rgb["denim blue"])
    ax.plot(x, y2, z,sns.xkcd_rgb["amber"])

    for _ in range(0,9):
        theta = [x + 0.2*np.pi for x in theta]
        y1 = r * np.sin(theta)
        y2 = - r * np.sin(theta)
        z = r * np.cos(theta)
        ax.plot(x, y1, z, sns.xkcd_rgb["denim blue"])
        ax.plot(x, y2, z, sns.xkcd_rgb["amber"])
    return ax

def mesh_plot_2d (len_mesh,rad_mesh,beta,resolution):
    fig,ax = plt.subplots(figsize=(10,6))
    theta = np.linspace(0, 3 * np.pi, len_mesh*resolution)
    x = np.linspace(0, len_mesh, len_mesh*resolution)
    r = rad_mesh * np.ones(len_mesh*resolution)
    y1 = r * np.sin(theta)
    y2 = - r * np.sin(theta)
    ax.plot(x, y1,sns.xkcd_rgb["denim blue"])
    ax.plot(x, y2,sns.xkcd_rgb["amber"])

    for _ in range(0,9):
        theta = [x + 0.2*np.pi for x in theta]
        y1 = r * np.sin(theta)
        y2 = - r * np.sin(theta)
        ax.plot(x, y1,sns.xkcd_rgb["denim blue"])
        ax.plot(x, y2,sns.xkcd_rgb["amber"])
    ax.set_ylim([-rad_mesh*1.1,rad_mesh*1.1])
    ax.set_title("Mesh body with length "+str(len_mesh)+" mm and radiux "+str(rad_mesh)+" mm")
    ax.set_aspect('equal')
In [3]:
# Measurements (in mm)
len_mesh = 60
rad_mesh = 10
beta = 40 # disected rhombus angle (deg)

# Plot parameter
resolution = 200 # number of points for each mm

mesh_3d = mesh_plot_3d (len_mesh,rad_mesh,beta,resolution)
mesh_2d = mesh_plot_2d (len_mesh,rad_mesh,beta,resolution)

NiTi actuator spring manufacturing

In [4]:
# Calculate the amount of NiTi wire needed to make the coil
def wire_per_coil_length (len_spr, rad_spr, rad_wir):
    return (rad_spr + rad_wir) * np.pi * len_spr / rad_wir
In [5]:
# Measurements (in mm)
len_spr = 30.0 # desired length of spring coil 
rad_spr = 0.2 # inner radius of the spring coil
rad_wir = 0.1 # radius of the NiTi wire

len_wir = wire_per_coil_length (len_spr, rad_spr, rad_wir)
print ('{:.2f} mm'.format(len_wir))

282.74 mm

Actuator Manufacturing

In [6]:
# Calculation the spring length needed to wrap around the mesh
def circular_actuator_length (rad_mesh,rad_spr,rad_wir):
    #r = rad_mesh + rad_spr + 2 * rad_wir
    r = rad_mesh + rad_wir
    return np.pi * r * 2

Notice that we want to know the radius around the mesh when

  1. the mesh is in its unloaded state (no pressure)
  2. the spring is in its martensite (extended form)
In [7]:
# Calculate effective displacement (length difference during transformation) of the coil
def effective_displacement (len_spr,rad_spr,rad_wir,shear_modulus,load_force):
    n_coil = len_spr / rad_wir # number of coils in the spring
    D = (rad_spr*2+rad_wir*4) # over all diameter
    C = D / (rad_wir*2) # spring index
    kappa = 0.615/C # stress correction factor
    tau = (8*load_force*D*kappa)/(np.pi*(rad_wir*2)**3)
    gamma = tau/shear_modulus
    free_length = (np.pi*gamma*n_coil*D**2)/((rad_wir*2)*kappa)
    # assume that detwinning cancels out
    detwinning_mart = 0 # estimate later
    detwinning_aust = 0 # estimate later
    return free_length + detwinning_mart - detwinning_aust

The following calculation is using the austenite instead until the compression calculation is in place

In [8]:
len_spr = circular_actuator_length (rad_mesh,rad_spr,rad_wir)
len_wir = wire_per_coil_length (len_spr, rad_spr, rad_wir)
print ('{:.2f} mm'.format(len_wir))

598.10 mm

The following calculation is estimated until springs are manufactured

In [9]:
# Measurements (in m,N,Pa)
rad_mesh = 5.0*10**-3 # desired radius of the mesh when compressed
rad_spr = 0.1*10**-3 # inner radius of the spring coil
rad_wir = 0.05*10**-3 # radius of the NiTi wire
load_force = 432.5*10**-3 # force from the mesh (estimate from paper)
shear_modulus = 31.125*10**9 # shear modulus of the wire in austenite state (estimate)

len_spr = circular_actuator_length (rad_mesh,rad_spr,rad_wir) # already in meters
len_spr *= 5
ed = effective_displacement (len_spr,rad_spr,rad_wir,shear_modulus,load_force)
len_wir = wire_per_coil_length (len_spr, rad_spr, rad_wir)
print ('length of the spring with the given spec: {:.2f} mm'.format(len_spr*10**3))
print ('effective displacement with the given spec: {:.2f} mm'.format(ed*10**3))
print ('NiTi needed per actuator with the given spec: {:.2f} mm'.format(len_wir*10**3))

length of the spring with the given spec: 158.65 mm
effective displacement with the given spec: 225.74 mm
NiTi needed per actuator with the given spec: 1495.25 mm

Actuation Effect

In [10]:
def custom_mesh_plot_2d (len_mesh,rad_mesh,beta,resolution):
    fig,ax = plt.subplots(figsize=(10,6))
    theta1 = np.linspace(0, 3 * np.pi, len_mesh*resolution)
    theta2 = [x+0.2*np.pi for x in theta1]
    x = np.linspace(0, len_mesh, len_mesh*resolution)
    
    mu = 10
    variance = 20
    sigma = np.sqrt(variance)
    
    r = rad_mesh * np.ones(len_mesh*resolution) - 30 * mpl.mlab.normpdf(x, mu, sigma)
    y1 = r * np.sin(theta1)
    z1 = r * np.cos(theta1)
    y2 = r * np.sin(theta2)
    z2 = r * np.cos(theta2)
    ax.plot(x, y1,sns.xkcd_rgb["denim blue"])
    ax.plot(x, y2,sns.xkcd_rgb["amber"])

    for _ in range(0,5):
        theta1 = [x+0.4*np.pi for x in theta1]
        theta2 = [x+0.4*np.pi for x in theta2]
        y1 = r * np.sin(theta1)
        z1 = r * np.cos(theta1)
        y2 = r * np.sin(theta2)
        z2 = r * np.cos(theta2)
        ax.plot(x, y1,sns.xkcd_rgb["denim blue"])
        ax.plot(x, y2,sns.xkcd_rgb["amber"])
    ax.set_ylim([-rad_mesh*1.1,rad_mesh*1.1])
    ax.set_aspect('equal')
In [11]:
# Measurements (in mm)
len_mesh = 60
rad_mesh = 10
beta = 40 # disected rhombus angle (deg)

# Plot parameter
resolution = 200 # number of points for each mm

mesh_2d = mesh_plot_2d (len_mesh,rad_mesh,beta,resolution)
custom_mesh_2d = custom_mesh_plot_2d (len_mesh,rad_mesh,beta,resolution)
In [12]:
# single rhombus extension due to compression
l = 1.0 
beta = np.linspace(0, np.pi/2, 100)
t = l * np.sin(beta)
delta = 0.2 * np.ones(100)
tn = t - delta
ext = np.sqrt(l**2 - tn ** 2) - np.cos(beta)

fig,ax = plt.subplots(figsize=(10,6))
ax.plot(beta, ext,sns.xkcd_rgb["denim blue"])

plt.show()
In [ ]:
In [ ]: