Cable Composite Object Using Wrong Area Moment of Inertia

[sriddle97]

Intro

Hello,

I am a grad student using Mujoco to simulate a soft bodied robot made with the cable composite elements. Despite having the proper dimensions and Twist/Bend properties for the material I noticed my deflections were not matching the physical robot we used as our comparison point so I did a few tests and dove into the base code.

My setup

Mujoco, Python

What's happening? What did you expect?

In the plugin file cable.cc (https://github.com/google-deepmind/mujoco/blob/main/plugin/elasticity/cable.cc#L134) lines 180-201 show the area moment of inertia calculation for the cable elements. I have selected capsule in the xml code which should mean Iy, Iz, and J for the circular cross section are used (lines 184-185). However, when I run a simple beam bending test, maintaining small deflection to conserve the small angle approximation, it produces stiffness that lines up almost perfectly with the rectangular cross section equations of the box elements (lines 193-195). Is it possible that the if statement in line 180/181 is not being executed properly and is instead always doing the box geom calculations for some reason?

Steps for reproduction

1. Put the python code (copied into a python file) and the xml model (in the Minimal Model section) in the same folder. Note that the model is attached as a txt file and will need to be copy/pasted into an xml file as xml is apparently not a supported attachment type.
2. Run the code.
3. You will see the rectangular area moment of inertia calculation and the measured area moment of inertia line up while the circular cross section, which should be used for the capsule type, is off by a lot.

Minimal model for reproduction

cable_xml.txt

Code required for reproduction

import numpy as np
import mujoco
import mediapy as media
import matplotlib.pyplot as plt
import os
import math

work_dir = os.getcwd()
xml_path = work_dir + r"\cable.xml"
model = mujoco.MjModel.from_xml_path(xml_path)
data = mujoco.MjData(model)

mujoco.mj_forward(model, data)

muscle_ind = 0      # might be 2??
tmax = 4
dt = 0.1/1000
t = np.arange(0, tmax, dt)

model.opt.timestep = dt

max_act = 0.01
force_app = 0.1

ctrl_ramp = np.ones(len(t))*max_act
ctrl_ramp[0:int(len(t)/8)] = np.linspace(0,max_act,int(len(t)/8))
frames = []
framerate = 60
muscle_force = np.zeros(np.size(t))
x_pos = np.zeros(np.size(t))
tendon_length = np.zeros(np.size(t))
slider_id = mujoco.mj_name2id(model,mujoco.mjtObj.mjOBJ_BODY,"slider")
for i in range(len(t)):
    data.xfrc_applied[slider_id] = [force_app, 0, 0, 0, 0, 0]
    mujoco.mj_step(model, data)
    x_pos[i] = data.xpos[slider_id,0]
# print(ctrl_ramp)

E_mod = 3e8
# E_mod = 10e8
cable_L = 0.3

radius = 0.00299
I_calc = (math.pi/4)*pow(radius,4)
print("Calculated Area Moment of Inertia for Circular C/S: ", I_calc)

I_calc_sq = pow(2*radius,3)*2*radius/12
print("Calculated Area Moment of Inertia for Rectangular C/S: ", I_calc_sq)

print("Force = ", force_app)
Disp = np.max(x_pos)
print("Displacement = ", Disp)

I_meas = (force_app*pow(cable_L,3))/(Disp*3*E_mod)
print("Measured Area Moment of Inertia: ", I_meas)

Confirmations

• I searched the latest documentation thoroughly before posting.
• I searched previous Issues and Discussions, I am certain this has not been raised before.

[quagla]

Hi this is unlikely as we have a test for the correctness of the bending stiffness of the capsule https://github.com/google-deepmind/mujoco/blob/main/test/plugin/elasticity/elasticity_test.cc#L255 Are you sure about your formula for I_meas?

[sriddle97]

Hello, I am not sure how the built-in correctness test works but I am certain my I_meas is correct. The code I included above performs a cantilever, end-loaded beam bending analysis and the I_meas is calculated using max displacement, beam length, elastic modulus, and applied force all of which are prescribed beforehand or measured during the simulation. The equation I used is a re-organized version of a common statics beam bending equation found here https://mechanicalc.com/reference/beam-deflection-tables, where I multiplied both sides by I and divided both sides by disp, ie. swapped the I and disp variables. I made sure to run the simulation long enough that the beam stopped moving to achieve max displacement. The output I get when I run the code is as follows:

Calculated Area Moment of Inertia for Circular C/S: 6.277325295188244e-11
Calculated Area Moment of Inertia for Rectangular C/S: 1.0656718401333334e-10
Force = 0.1
Displacement = 0.028637255408429853
Measured Area Moment of Inertia: 1.0475864244717039e-10

Given that the measured value is off by just 2% for the rectangular I and a significant 49% for the circular I, it's clear to me that the rectangular is being used in the simulation.

If I made a mistake while constructing the model I am happy to fix it but as you can see in the included xml (.txt file) I have the geom type set to capsule. The only other place I could see being an issue is the default section in lines 16-24 where I do not define a type in the geom part, but when I commented that out it produced the same results.

[quagla]

2% seems a bit off to draw a conclusion.. Could you debug it by putting a printf directly here https://github.com/google-deepmind/mujoco/blob/main/plugin/elasticity/cable.cc#L200 ?

[sriddle97]

I ran a further analysis with a variety of different conditions and it seems that my initial method was flawed. By fixing the first element of the cable I inadvertently reduced the bendable length of the beam but used the full length in the calculation. Increasing the number of elements in the cable got me closer and closer to the expected circular area moment of inertia value as the fixed length decreased and the bendable length increased closer to the full beam length. Apparently it was just coincidental that the rectangular calculation was very close.