Small scale rod simulation

[Junang-Wang]

Hello,

I'm currently working on a project that involves simulating a small-scale rod with specific dimensions: length = 3mm, radius = 0.3mm, Young's modulus E = 14MPa, shear modulus G = E/3, and a density of 2.3 g/cm^3. I've encountered an issue where the time step needs to be set at a very low value (1e-7s) to prevent numerical instability, which makes the simulation extremely time-consuming. Does anyone have any insights or suggestions on how to address this challenge when simulating a small-scale rod?

Thank you in advance for any assistance or advice you can provide!

[Junang-Wang]

The following is my code snippet.

class PendulumSimulator(BaseSystemCollection, Constraints, Forcing, Damping): pass

Pendulum_Sim = PendulumSimulator()

# Create Rod 
n_elem = 40

#--------------mmGS(mm, mg, s) Unit------------------
density = 2.273  # mg/mm^3
gravitational_acc = -9.80665e3  #mm/s^2
base_length = 3 #mm 
base_radius = 0.3 #mm
E = 1.4e9 # mg/(mm*s^2)

dl = base_length / n_elem 
# dt = 0.002 * dl 
dt = 1.0e-7
nu = 1500 #damping constant
shear_modulus = E / (3.0)

start = np.array([0.0,0.0,0.0])
direction = np.array([0.0, 1.0, 0.0])
normal = np.array([1.0, 0.0, 0.0])
base_area = np.pi * base_radius**2
# base_area = base_radius*np.sqrt(3)/3

origin_force = np.array([0.0, 0.0, 0.0])
end_force = np.array([0.0, 0.0, -end_force_zs[i]])
ramp_up_time = 0.1

softer_rod = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    youngs_modulus=E,
    shear_modulus = shear_modulus,
)

stiffer_rod = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    youngs_modulus = E*10,
    shear_modulus = shear_modulus*10,
)

Pendulum_Sim.append(softer_rod)
Pendulum_Sim.append(stiffer_rod)
Pendulum_Sim.dampen(softer_rod).using(
    AnalyticalLinearDamper,
    damping_constant = nu,
    time_step = dt,
)
Pendulum_Sim.dampen(stiffer_rod).using(
    AnalyticalLinearDamper,
    damping_constant = nu,
    time_step = dt,
)

Pendulum_Sim.constrain(softer_rod).using(
    OneEndFixedRod,constrained_position_idx = (0,), constrained_director_idx= (0,)
)
Pendulum_Sim.constrain(stiffer_rod).using(
    OneEndFixedRod,constrained_position_idx = (0,), constrained_director_idx= (0,)
)
print("One end of the rod is now fixed in place")

Pendulum_Sim.add_forcing_to(softer_rod).using(
    EndpointForces, origin_force, end_force, ramp_up_time = ramp_up_time
)
Pendulum_Sim.add_forcing_to(stiffer_rod).using(
    EndpointForces, origin_force, end_force, ramp_up_time=ramp_up_time
)
print("Forces added to the rod")

Pendulum_Sim.finalize()
print("System finalized")

#----------------------------time integration----------------------------------
final_time = 2.0
total_steps = int(final_time / dt)
timestepper = PositionVerlet()
integrate(timestepper, Pendulum_Sim, final_time, total_steps)

[Ali-7800]

Hi @Junang-Wang,

What I think is causing the numerical issues in your simulation is that your Young's modulus value is too big compared to the other values of your parameters causing double-precision float arithmetic issues when solving the cosserat rod equations.

A possible solution is to scale the Young's modulus and the forces in the simulation to avoid the those issues, For example you can apply a force F on a rod with Young's modulus E to it to stretch it by dx, but you can also apply a force F/1000 on a rod with Young's modulus E/1000 to get the same stretch value dx given that everything else is the same. Just keep in mind that if you need to take a look at the internal forces later, that the forces there are also scaled.

One other note, your damping value is quite high which might affect the dynamics of your simulation, I would recommend keeping it low in order to have a more dynamically accurate simulation.

[Junang-Wang]

Hi @Ali-7800,

You made an excellent point regarding the factor causing numerical numerical issue. I think you are kind of corrent. In static state, I believe the result by using your method is as same as the origin one. However, for dynamtic states, the two results may not be the same. As the density remains unchanged, we apply a force F/1000 on a mass m which leads to a acceleration a/1000. So the dynamic path may be different. What do you think? Should I scale the density 1000 times to fix this problem.

Additionally, you mentioned that the damping value in my code appears to quite high. Do you know any guidelines or recommended process I can follow to choose a more suitable value?

Thank you once again for your prompt response.

[Ali-7800]

Yes, you should also adjust the density as well to account for that if you care about the dynamics. For the damping I would recommend starting as close as possible to zero then increasing it very slightly to make sure that the dynamics don't change as the damping increases, once you get a stable simulation that is dynamically accurate I would stop.

[AaronRapopot]

Hellow, I met a same problem like yours, a small scale rod simulation. I check the position of the rod. The shape of the position is (250, 3, 101). There were a total of 250 steps. In the initial dozens of steps, numerical values were still present (although highly unstable), but all subsequent values displayed "nan". Have you met a similar situation? I tried to set a small time step value but it did't works. Thanks for your help. Like this: ''' array([[ 0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan], [ 0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan], [ 0., nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan]]) '''

[Junang-Wang]

Hi @Ali-7800,

I scale down the mass m to m/1000 as we discussed earlier. And we found calculation timesteps have to be extremely small, which will significantly elongate total calculation time. We thought it might be due to the small moment of inertia of elements in small-scale soft robots, which led to enormous velocity increments in a single timestep. Any solutions in your mind can solve this problem?

My system tends to diverge easily when the damping constant approaches zero. Could you please advise if there is anything I am doing wrong?

Thank you again for your help. Your assistance is greatly appreciated.

[Junang-Wang]

Hellow, I met a same problem like yours, a small scale rod simulation. '''

Hi @AaronRapopot, you can check what @Ali-7800 has discussed above. @Ali-7800 has provided a nice approximation method.

[Junang-Wang]

Hi @Ali-7800,

Another quick question. Do I need to scale down torque, when I scale down young's modulus and shear modulus according to our previous discussion?

A possible solution is to scale the Young's modulus and the forces in the simulation to avoid the those issues, For example you can apply a force F on a rod with Young's modulus E to it to stretch it by dx, but you can also apply a force F/1000 on a rod with Young's modulus E/1000 to get the same stretch value dx given that everything else is the same. Just keep in mind that if you need to take a look at the internal forces later, that the forces there are also scaled.

[bhosale2]

@Junang-Wang, similar to what @Ali-7800 mentioned, scaling down torque maybe needed if you need to maintain the same angular deflection, since scaling down E also scales down the rods bending and torsional stiffness. A good way to approach these systems is non dimensionalize the parameters involved for your reference, which gives you the guidance by what factors should you scale your forces and torques. @Ali-7800 feel free to add or correct as needed.

[bhosale2]

@Junang-Wang Has your issue been resolved? Please let us know, or else we shall close this issue due the inactivity in a few days.

[bhosale2]

@Ali-7800 and @armantekinalp closing this issue due to inactivity, as mentioned above. Feel free to open if deemed necessary.