Muscular Snake Example Questions

[jam474cwru]

Hello! I particularly would like to inquire about the muscular snake example. What is the difference between nu and poissons ratio when defining the bodies? what about n_elem_body? How are the force amplitudes determined and for which muscles (I see the equations in theory but wanted answers relating to this example)?

I really appreciate your help through these questions

[jam474cwru]

I also would like to ask effects of changing variables of the snake body. My goal is to observe the effects of changing mechanical variables like density, modulus of elasticity, and poissons ratio to identify the effects on the flexure and speed of flexure of the snake. Is there a proper way to observe this other than just changing the respective inputs in the code?

[armantekinalp]

@xzhan139 could you please take care this?

[jam474cwru]

I have changed those key inputs to the snake body and notice the same result in terms of video. Does something with muscles also need to be changed? I am seeking to see what effects the flexure of the snake in addition to my previous questions. I am eager to hear the responses and discuss. Thank you.

[xzhan139]

Hello - Sorry for the delayed response.

"Hello! I particularly would like to inquire about the muscular snake example. What is the difference between nu and poissons ratio when defining the bodies? what about n_elem_body? How are the force amplitudes determined and for which muscles (I see the equations in theory but wanted answers relating to this example)?"

• nu is a damping term that aims to capture the environmental viscoelasticity and internal damping effect due to all the connectivity. It is also a meaningful term in terms of numerical stability, such that a higher nu usually leads to more stable behavior. In the code, the overall effect of nu is adding a small penalty on the velocity of the filament. Poissons ratio, instead, defines the deformation behavior of the filament, just like how it is classically defined in mechanics. n_elem_body is the number of discretized element of any body (filament) in the system. It may be different for each element, say the n_elem_body of the snake body could be different from the n_elem_body of the muscles. Force amplitude is determined through an optimization course, through which we determine the optimal forces of each muscle so that the snake will move fastest forward. You can refer to both RSOS 2018 and Nat Comm 2019 papers for more detailed information about all these parameter.

"I also would like to ask effects of changing variables of the snake body. My goal is to observe the effects of changing mechanical variables like density, modulus of elasticity, and poissons ratio to identify the effects on the flexure and speed of flexure of the snake. Is there a proper way to observe this other than just changing the respective inputs in the code?"
"I have changed those key inputs to the snake body and notice the same result in terms of video. Does something with muscles also need to be changed? I am seeking to see what affects the flexure of the snake in addition to my previous questions. I am eager to hear the responses and discuss. Thank you."

• The term that will change the flexure of the snake most is the Young's modulus of the snake body. Or, if going a little bit more in detail, is the flexure rigidity (or moment of area) in the horizontal (binormal) direction of the snake. You can hack in all of these terms to change its bending behavior. Just be mindful that if a snake gets too rigid, you may have numerically instability that leads to divergence. In terms of changing the parameters, if you have a certain objective in mind, you can always run an optimization course for parameters selection. Again, there are reference on this in both papers I mentioned above. The alternative is that you can write a small code that runs in parallel a sweep of parameters parsing to your elastica simulation. So that you don't have to manually run each case.

Hope these will help!

[jam474cwru]

Thank you so much for your response! Even hacking in different numbers for density, modulus of elasticity, and poissons ratio, I have not noticed any sort of difference in the final plotting of the video. The speed of movement, how much the snake flexes, and how fast it flexes seems to be the same. Is there a different way to visualize these concepts we are discussing?

How would you advise going about attempting the small code that runs a parallel sweep of parameters parsing the the elastica simulation? I am getting back into coding and am relatively rusty. I appreciate your help in this discussion.

[xzhan139]

Hi -
"Even hacking in different numbers for density, modulus of elasticity, and poissons ratio, I have not noticed any sort of difference in the final plotting of the video. The speed of movement, how much the snake flexes, and how fast it flexes seems to be the same."

Can you send me a couple of examples on this. What parameter did you change, how much did you change, and the corresponding video? Normally, increasing E of the snake body should give you a visible difference in terms of movement. Another general advice is that maybe you can consider whether your problem can be simplified and translated to "continuumsnake" case. The main difference is that in there, you don't have muscles as real cosserat rod, but you can still apply torque to the snake body at any location to micmic the muscular actuation. The benefit is then, your simulation is much simpler, much quicker, and the experimentation on the parameters can be done much easier.

"How would you advise going about attempting the small code that runs a parallel sweep of parameters parsing the the elastica simulation? I am getting back into coding and am relatively rusty. I appreciate your help in this discussion."

The Elastica simulation can take in arguments of your choice. For example:

if name == "main":
parser = argparse.ArgumentParser()
parser.add_argument(
"--wave_length", type=float, default=0.5,
)
parser.add_argument(
"--timestep", type=float, default=5e-4,
)
parser.add_argument(
"--damping", type=float, default=1e-1,
)
args = parser.parse_args()

here I have wave_length, timestep and damping as parameters to run sweep on. Then you just need to create your own sweep runner code, using any multiprocessing tool you like to run through all the parameters. One example could be using p_map in p_tqdm to run each of your Elastica simulations as subprocesses. Here are some references:
https://github.com/swansonk14/p_tqdm
https://docs.python.org/3/library/subprocess.html

[jam474cwru]

https://vscode.dev/github/GazzolaLab/PyElastica/blob/master/examples/MuscularSnake/3D_muscular_snake.mp4

Hello, I guess with drastic changes in E, I am kind of able to see it. Here I used the E for TPU as 162E6. It shows the rod for a second and disappears. I am assuming it is too rigid and it diverges? I also used Poissons ratio as 0.37 and density as 1235 kg/m^3. These numbers gave me an awkward result.

On the other hand, I also used the values for PLA as E = 13E6, Poisson's = 0.35, and density as 1240 kg/m^3. These numbers gave me a very similar video to the original numbers in the code and was hard to tell any difference. Maybe because there were small differences compared to TPUs E value?

These are the only things I changed.

[xzhan139]

Hi - Yes, it is very likely that the simulation diverged when running those stiffness. In general, the stiffer the rod gets, the smaller tilmestep is required to keep the simulation running. The default timestep we put in the code are tuned so that it is almost the largest possible for it to run properly, also meaning that the simulation is close to fastest given its stiffness. In your case, E=162E6 is an order of magnitude larger than the original value (1e7), so that some serious parameter tuning is required. Some general advice is: tuning down the tilmestep, make the rod thicker (if possible), make the rod heavier (if possible) and increase the external damping (nu). Another advice is, maybe you are already doing this, other than visualizing the video, plotting the velocity/displacement of the snake will be useful to judge the difference when using closer values, say when you use E = 13E6.

[jam474cwru]

Further, I used these same values in the continuum snake example and did not receive any visual results using my values of PLA and received very weird results using the values for TPU. What could cause the lack of results when running this?

[xzhan139]

I think is again the same numerical diverging issue here. If I remember correctly, in continuum snake we use default E=1e6. So those numbers you were trying are too rigid for the default timestep to run. You can try similar tricks to tune it. Given the continuum snake has much simpler structure, it may be much easier job than tuning the muscular snake case. Another thing to try if you just want to confirm your simulation is running with different parameters, is to make it softer by tuning down the E. This is normally more straightforward.

[jam474cwru]

Thank you for your help. After much time of playing around, I am starting to understand more. I do have some follow up questions:

1. Is there a specific formula/equation for velocity that both the muscular and continuum snake example use? I have read the literature and wanted to specifically ask which equation provides this. I am in search of what characteristics of deformation determine velocity.

2. How can I label the x and y axes in the results of running these simulations? Where are the axis titles and names in the snake examples?

3. Under the continuum snake example, I am noticing it doesn't give the velocity vs time graph under its respective png file but shows it clearly in the normal continuum snake.png file. Are these the same?

I appreciate your help as always.

[xzhan139]

hi

1. Just to clarify a bit, both examples don't use any specific formula as input. The velocity is only the result of the physical settings. We haven't done too much of systematic studies on different snake characteristics vs snake velocity. But you can check this paper Zhang, X., Naughton, N., Parthasarathy, T. and Gazzola, M., 2021. Friction modulation in limbless, three-dimensional gaits and heterogeneous terrains. Nature communications, 12(1), p.6076.
   where you can find an online sandbox that allows you to set different input parameters and visualize the snake behavior. It is based on a simplified mathematical model, but should be a good reference for elastica simulations.
   2 and 3. All the plotting functions are in a separate file "continuum_snake_postprocessing.py". You can go there set all the information you need.

[jam474cwru]

Thank you so much for your timely response! It really helps me move forward in my research. I just need a little help understanding further. How then are the velocity graphs obtained with the forward, normal, and lateral graphs in elastica? I have changed modulus of elasticity, density, and poissons ratio in the muscular and continuum snake examples and received different visualizations of the snake motion, but do not understand how it changes velocity.

In the article, I am seeing pose angle and steering rate effect undulatory motion (Figure 3) of the snake and the effective speed, but I want to understand exactly how (with the goal being to see how geometry of the snake does effect and maximize snake velocity as this is my research). What formula gives those values? I don't understand what is meant by physical settings.

Thank you always for your help.

[jam474cwru]

Hello. I would like to check on the above. Thanks!

[xzhan139]

Hi:

1. Like I wrote in the response for your previous question, all functions related to plotting is in "continuum_snake_postprocessing.py". There you have functions for retrieving all the information saved at the end of the snake simulation, and then process it for plotting. Say"plot_snake_velocity" is one example. It should be quite straightforward to figure out what these functions do.
2. "with the goal being to see how geometry of the snake does effect and maximize snake velocity as this is my research" -
   Sorry for the confusion, what I meant is that the paper itself may not directly relate to your question, but the sandbox attached to that paper may be doing exactly what you need. I attach the link here: https://gazzolalab.github.io/kinematic_snake_sandbox/snake_sandbox.html
   You can define your activation function and other snake properties, and visualize the results. Again, be mindful that this is based on a simplified math model, not the full elastica simulation. It works as a quick, user-friendly GUI for parameter search for whatever snake simulation you want to run, and then you can confirm using elastica.

[jam474cwru]

Thank you for the help. I see the the program function for compute_projected_velocity. I am referring to the exact equation to receive that function. Maybe there is a misunderstanding somewhere. I just am not understanding how that function provides a result with no mathematical equation. The green notes in post processing say that rod velocity needs computed in rod direction with respect to the new tangents and roll directions. How is it computing any velocity? All I am looking for are any equations in regards to velocity for the elastica simulations.

[jam474cwru]

Upon another look in the article, I believe it is equations 2.12-2.15 that help define velocity as the slithering snake example also points to. Can you clarify this? I believe this may be the answer. Maybe just more explanation regarding it? Thank you so much.

https://www.cosseratrods.org/publications/pubs/2018_RSOS.pdf

[xzhan139]

Hi,
If you go line-by-line in compute_projected_velocity:

Compute the projected rod velocity in the direction of the rod

velocity_mag_in_direction_of_rod = np.einsum(
    "ji,i->j", avg_velocity, direction_of_rod
)
velocity_in_direction_of_rod = np.einsum(
    "j,i->ji", velocity_mag_in_direction_of_rod, direction_of_rod
)

these einstein sum functions are for computing the project.
Reference if you are not familiar with einstein sum:
https://numpy.org/doc/stable/reference/generated/numpy.einsum.html

[jam474cwru]

I see. The explanation of Einstein functions states it calculates the array of operations. Which operations is it referring to? What exactly are these einstein functions pulling to calculate?

[jam474cwru]

I would assume its the system of equations to get velocity from the article you sent (Equations 2.12-2.15) right?:

https://www.cosseratrods.org/publications/pubs/2018_RSOS.pdf

[jam474cwru]

Good Afternoon, I wanted to follow up on the answer for the above. Can I please have a concise answer as to the exact equations that give the velocity results for the simulations? I want to confirm they are equations 2.12-2.15 in the link that you gave. I just need the actual equations. Which operations are the einstein sums referring to or what are they pulling to calculate? This is the last piece of information I am in need of. Please advise. Thank you for the help!

https://www.cosseratrods.org/publications/pubs/2018_RSOS.pdf