How to correctly simulate a pre-twisted rod with one end fixed and one end free?

[Qin-Shaowen]

I am trying to simulate a rod that is initially twisted and then released, with one end fixed and the other end free to untwist. The setup:

• Rod length: 0.4 m, radius: 1.2 mm, circular cross-section
• Young's modulus: 0.8 MPa, density: 1200 kg/m³
• One end is fixed, the other end is attached to a rigid mass (modeled as a short stiff rod segment)
• The rod is initially twisted (10 turns) and then released

I am using PyElastica 1.0.0 with Python 3.13.7.

I have tried:

1. Setting rod.rest_sigma[2, :] = initial_twist_rate → all displacement data becomes nan after the first timestep.
2. Manually rotating each element's director_collection to create an initial twisted geometry, while keeping rest_sigma zero → still nan.
3. Writing a custom constraint class based on HelicalBucklingBC to apply twist to one end for a short time, then release → still nan.
4. Reducing twist to only 2 turns, using very small time steps (2e-5), adding strong damping, using fewer elements → still nan.

What is the correct way to initialize a CosseratRod with a uniform pre-twist and simulate its free untwisting?

Thank you.

[skim0119]

@Qin-Shaowen Please share your script.

Just to clarify, you want to simulate a rod that has rest configuration to be straight, and set the simulation that start with twisted configuration?

• rest_sigma is shear/stretch, so it is irrelevant.
• (2) is probably the closest one.

[Qin-Shaowen]

Yes, you are correct. The rod's rest configuration is straight, but I want the simulation to start from a twisted state.

Here is the script I tried for method (2). It manually rotates the director_collection of each element to create the initial twist, while keeping rest_sigma zero. However, all displacement data in the CSV files ends up as nan after the first timestep.

import numpy as np
import matplotlib.pyplot as plt
from elastica import *

rod_length = 0.4
radius = 1.2e-3
density = 1200.0
E = 0.8e6
G = E / 3.0

I = np.pi * radius**4 / 4.0
J_polar = np.pi * radius**4 / 2.0

ball_mass = 20.0e-3
ball_radius = 0.01

twist_turns = 2
total_twist_angle = twist_turns * 2.0 * np.pi

simulation_time = 2.0
time_step = 2.0e-5
n_steps = int(simulation_time / time_step)
output_interval = 1000

marker_spacing = 0.04
marker_positions = np.arange(0, rod_length + 1e-12, marker_spacing)

class RubberBandBall(BaseSystemCollection, Constraints, Forcing, CallBacks, Damping, Connections):
    pass

system = RubberBandBall()

rod = CosseratRod.straight_rod(
    n_elements=200,
    start=np.array([0.0, 0.0, 0.0]),
    direction=np.array([0.0, 0.0, -1.0]),
    normal=np.array([1.0, 0.0, 0.0]),
    base_length=rod_length,
    base_radius=radius,
    density=density,
    youngs_modulus=E,
    shear_modulus=G,
)

n_elems = rod.n_elems
s_vec = np.linspace(0, rod_length, n_elems + 1)
s_elems = (s_vec[:-1] + s_vec[1:]) / 2.0
twist_angles = total_twist_angle * s_elems / rod_length

init_directors = rod.director_collection.copy()
for i in range(n_elems):
    angle = twist_angles[i]
    cos_a = np.cos(angle)
    sin_a = np.sin(angle)
    Rz = np.array([[cos_a, -sin_a, 0],
                   [sin_a,  cos_a, 0],
                   [0,       0,     1]])
    init_directors[:, :, i] = Rz @ init_directors[:, :, i]
rod.director_collection = init_directors
rod.rest_sigma = np.zeros((3, n_elems))

bend_matrix = np.zeros((3, 3, n_elems - 1))
bend_matrix[0, 0, :] = E * I
bend_matrix[1, 1, :] = E * I
bend_matrix[2, 2, :] = G * J_polar
rod.bend_matrix = bend_matrix

system.append(rod)

ball_length = 0.02
ball_n_elements = 4
ball_radius_equiv = 0.01
ball_area = np.pi * ball_radius_equiv**2
ball_density = ball_mass / (ball_area * ball_length)

ball_rod = CosseratRod.straight_rod(
    n_elements=ball_n_elements,
    start=rod.position_collection[..., -1].copy(),
    direction=np.array([0.0, 0.0, -1.0]),
    normal=np.array([1.0, 0.0, 0.0]),
    base_length=ball_length,
    base_radius=ball_radius_equiv,
    density=ball_density,
    youngs_modulus=E * 1e2,
    shear_modulus=G * 1e2,
)

last_twist_angle = total_twist_angle
cos_a = np.cos(last_twist_angle)
sin_a = np.sin(last_twist_angle)
Rz_end = np.array([[cos_a, -sin_a, 0],
                   [sin_a,  cos_a, 0],
                   [0,       0,     1]])
for i in range(ball_n_elements):
    ball_rod.director_collection[:, :, i] = Rz_end @ np.eye(3)
ball_rod.rest_sigma = np.zeros((3, ball_n_elements))

ball_I = np.pi * ball_radius_equiv**4 / 4.0
ball_J = np.pi * ball_radius_equiv**4 / 2.0
ball_bend = np.zeros((3, 3, ball_n_elements - 1))
ball_bend[0, 0, :] = E * 1e2 * ball_I
ball_bend[1, 1, :] = E * 1e2 * ball_I
ball_bend[2, 2, :] = G * 1e2 * ball_J
ball_rod.bend_matrix = ball_bend

system.append(ball_rod)

system.constrain(rod).using(
    OneEndFixedBC,
    constrained_position_idx=(0,),
    constrained_director_idx=(0,),
)

system.connect(rod, ball_rod, first_connect_idx=-1, second_connect_idx=0).using(
    FixedJoint, k=1e5, nu=2.0, kt=1e3, nut=0.5,
)

gravity_vec = np.array([0.0, 0.0, -9.81])
system.add_forcing_to(rod).using(GravityForces, acc_gravity=gravity_vec)
system.add_forcing_to(ball_rod).using(GravityForces, acc_gravity=gravity_vec)

system.dampen(rod).using(AnalyticalLinearDamper, damping_constant=0.3, time_step=time_step)
system.dampen(ball_rod).using(AnalyticalLinearDamper, damping_constant=0.3, time_step=time_step)

class Recorder(CallBackBaseClass):
    def __init__(self, rod, ball_rod, marker_positions, output_interval):
        self.rod = rod
        self.ball_rod = ball_rod
        self.marker_positions = marker_positions
        self.output_interval = output_interval
        self.time = []
        self.ball_omega_z = []
        self.marker_x = [[] for _ in marker_positions]
        self.marker_y = [[] for _ in marker_positions]
        self._init_csv_files()

    def _init_csv_files(self):
        with open('ball_omega.csv', 'w') as f:
            f.write('Time(s),Omega_z(rad/s)\n')
        for i, s in enumerate(self.marker_positions):
            label = int(s * 100)
            with open(f'marker_{label:02d}cm_x.csv', 'w') as f:
                f.write('Time(s),x(m)\n')
            with open(f'marker_{label:02d}cm_y.csv', 'w') as f:
                f.write('Time(s),y(m)\n')

    def make_callback(self, system, time, current_step):
        if current_step % self.output_interval == 0:
            self.time.append(time)
            omega_z = self.ball_rod.omega_collection[2, -1]
            self.ball_omega_z.append(omega_z)
            pos = self.rod.position_collection
            s_nodes = np.linspace(0, rod_length, pos.shape[1])
            for i, s in enumerate(self.marker_positions):
                x = np.interp(s, s_nodes, pos[0, :])
                y = np.interp(s, s_nodes, pos[1, :])
                self.marker_x[i].append(x)
                self.marker_y[i].append(y)
            with open('ball_omega.csv', 'a') as f:
                f.write(f'{time},{omega_z}\n')
            for i, s in enumerate(self.marker_positions):
                label = int(s * 100)
                with open(f'marker_{label:02d}cm_x.csv', 'a') as f:
                    f.write(f'{time},{self.marker_x[i][-1]}\n')
                with open(f'marker_{label:02d}cm_y.csv', 'a') as f:
                    f.write(f'{time},{self.marker_y[i][-1]}\n')

recorder = Recorder(rod, ball_rod, marker_positions, output_interval)
system.collect_diagnostics(rod).using(
    Recorder,
    rod=rod,
    ball_rod=ball_rod,
    marker_positions=marker_positions,
    output_interval=output_interval,
)

print("Simulation starting...")
system.finalize()
integrator = PositionVerlet()
integrate(integrator, system, final_time=simulation_time, n_steps=n_steps)

print("Simulation finished, but all marker data is nan.")

The simulation runs to completion, but all recorded marker displacements in the CSV files are nan after the first output step. The same happens even if I reduce the twist to 2 turns, use smaller time steps, or increase damping. Could you help me identify what is wrong with the initial twisted setup?

Thank you for your help.

To add more context about the original physical problem I'm trying to simulate:

The real-world experiment involves a straight rubber band (natural state) with a metal ball attached to its lower end. The ball is manually twisted many times around the vertical axis, and then released. After release, the rubber band untwists, and transverse standing waves can be observed on the band. The goal is to simulate this entire process — the initial twisted state, the release, and the resulting dynamics of both the ball and the rubber band.

Initially, I tried to use a Sphere as the end mass and connect it to the CosseratRod using a FixedJoint. However, I encountered a dimension mismatch error between the director_collection of the rod (3, 3) and the sphere (3, 3, 1) inside the joint's torque calculation. Since I couldn't resolve this, I temporarily replaced the sphere with a very short, very stiff CosseratRod segment (the "ball_rod" in my script) to act as an equivalent rigid mass.

Is there a proper way to attach a rigid body (Sphere) to the end of a CosseratRod that avoids the dimension mismatch? Any guidance on this would also be greatly appreciated.

[skim0119]

@Qin-Shaowen Here is some changes and comments for the code. I left Note comments where I made changes.
Mainly,

• Showed how to use ea.Sphere with fixed relative orientation (I believe this is what you want: link the angular momentum).
• Tuned some parameters

Most of the code should be the same. Please take a look.

Probably, few physical parameters (contact springs and dissipation) needs to be reconsidered to match the timescale exactly, but I would start from this.

bell_simulation.mp4

import numpy as np
import matplotlib.pyplot as plt
from elastica import *
from elastica._rotations import get_relative_rotation_two_systems

rod_length = 0.4
radius = 1.2e-3
density = 1200.0
E = 0.8e6
G = E / 3.0

I = np.pi * radius**4 / 4.0
J_polar = np.pi * radius**4 / 2.0

ball_mass = 20.0e-3  # seems
ball_radius = 0.01

twist_turns = 10  # 2 NOTE: Testing 10
total_twist_angle = twist_turns * 2.0 * np.pi

simulation_time = 10
time_step = 2.0e-5
output_interval = 1000

marker_spacing = 0.04
marker_positions = np.arange(0, rod_length + 1e-12, marker_spacing)


class RubberBandBall(
    BaseSystemCollection, Constraints, Forcing, CallBacks, Damping, Connections
):
    pass


system = RubberBandBall()

rod = CosseratRod.straight_rod(
    n_elements=200,
    start=np.array([0.0, 0.0, 0.0]),
    direction=np.array([0.0, 0.0, -1.0]),
    normal=np.array([1.0, 0.0, 0.0]),
    base_length=rod_length,
    base_radius=radius,
    density=density,
    youngs_modulus=E,
    shear_modulus=G,
)

n_elems = rod.n_elems
s_vec = np.linspace(0, rod_length, n_elems + 1)
s_elems = (s_vec[:-1] + s_vec[1:]) / 2.0
twist_angles = total_twist_angle * s_elems / rod_length

init_directors = rod.director_collection.copy()
for i in range(n_elems):
    angle = twist_angles[i]
    cos_a = np.cos(angle)
    sin_a = np.sin(angle)
    Rz = np.array([[cos_a, -sin_a, 0], [sin_a, cos_a, 0], [0, 0, 1]])
    init_directors[:, :, i] = Rz @ init_directors[:, :, i]
rod.director_collection = init_directors
# rod.rest_sigma = np.zeros((3, n_elems)) # Note: No need to set rest_sigma

bend_matrix = np.zeros((3, 3, n_elems - 1))
bend_matrix[0, 0, :] = E * I
bend_matrix[1, 1, :] = E * I
bend_matrix[2, 2, :] = G * J_polar
rod.bend_matrix = bend_matrix  # Note: Not sure what the purpose of this line

system.append(rod)

# Note: Replace with ea.Sphere
# ball_length = 0.02
# ball_n_elements = 1
# ball_radius_equiv = 0.01
# ball_area = np.pi * ball_radius_equiv**2
# ball_density = ball_mass / (ball_area * ball_length)

# ball_rod = CosseratRod.straight_rod(
#     n_elements=ball_n_elements,
#     start=rod.position_collection[..., -1].copy(),
#     direction=np.array([0.0, 0.0, -1.0]),
#     normal=np.array([1.0, 0.0, 0.0]),
#     base_length=ball_length,
#     base_radius=ball_radius_equiv,
#     density=ball_density,
#     youngs_modulus=E * 1e2,
#     shear_modulus=G * 1e2,
# )
sphere_density = ball_mass / (4.0 / 3.0 * np.pi * ball_radius**3)
sphere_center = rod.position_collection[..., -1].copy() + np.array(
    [0.0, 0.0, -ball_radius]
)
sphere = Sphere(sphere_center, ball_radius, sphere_density)

# Note: Why ball needs to be twisted?
# last_twist_angle = total_twist_angle
# cos_a = np.cos(last_twist_angle)
# sin_a = np.sin(last_twist_angle)
# Rz_end = np.array([[cos_a, -sin_a, 0], [sin_a, cos_a, 0], [0, 0, 1]])
# for i in range(ball_n_elements):
#     ball_rod.director_collection[:, :, i] = Rz_end @ np.eye(3)
# ball_rod.rest_sigma = np.zeros((3, ball_n_elements))

# ball_I = np.pi * ball_radius_equiv**4 / 4.0
# ball_J = np.pi * ball_radius_equiv**4 / 2.0
# ball_bend = np.zeros((3, 3, ball_n_elements - 1))
# ball_bend[0, 0, :] = E * 1e2 * ball_I
# ball_bend[1, 1, :] = E * 1e2 * ball_I
# ball_bend[2, 2, :] = G * 1e2 * ball_J
# ball_rod.bend_matrix = ball_bend

system.append(sphere)

system.constrain(rod).using(
    OneEndFixedBC,
    constrained_position_idx=(0,),
    constrained_director_idx=(0,),
)

system.connect(rod, sphere, first_connect_idx=-1, second_connect_idx=0).using(
    FixedJoint,
    k=1e2,  # Note: Keep it low.
    nu=0.0,
    kt=1e-3,  # Note: Keep it low.
    nut=0.0,
    rest_rotation_matrix=get_relative_rotation_two_systems(
        rod, -1, sphere, 0
    ),  # Note: Fix relative orientation between tip and sphere director
    # BELOW ARE GIVEN VALUES
    # k=1e5,
    # nu=2.0,
    # kt=1e3,
    # nut=0.5,
)

gravity_vec = np.array([0.0, 0.0, -9.81])
system.add_forcing_to(rod).using(GravityForces, acc_gravity=gravity_vec)
system.add_forcing_to(sphere).using(GravityForces, acc_gravity=gravity_vec)

# Note: Using uniform_damping, or separately treat translational and rotational damping
system.dampen(rod).using(
    AnalyticalLinearDamper,
    translational_damping_constant=1e-3,
    rotational_damping_constant=1e-8,
    time_step=time_step,
)
system.dampen(sphere).using(
    AnalyticalLinearDamper,
    uniform_damping_constant=1e0,
    time_step=time_step,
)


class Recorder(CallBackBaseClass):
    def __init__(self, rod, ball_rod, marker_positions, output_interval):
        self.rod = rod
        self.ball_rod = ball_rod
        self.marker_positions = marker_positions
        self.output_interval = output_interval
        self.time = []
        self.ball_omega_z = []
        self.marker_x = [[] for _ in marker_positions]
        self.marker_y = [[] for _ in marker_positions]
        self._init_csv_files()

    def _init_csv_files(self):
        with open("ball_omega.csv", "w") as f:
            f.write("Time(s),Omega_z(rad/s)\n")
        for i, s in enumerate(self.marker_positions):
            label = int(s * 100)
            with open(f"marker_{label:02d}cm_x.csv", "w") as f:
                f.write("Time(s),x(m)\n")
            with open(f"marker_{label:02d}cm_y.csv", "w") as f:
                f.write("Time(s),y(m)\n")

    def make_callback(self, system, time, current_step):
        if current_step % self.output_interval == 0:
            self.time.append(time)
            omega_z = self.ball_rod.omega_collection[2, -1]
            self.ball_omega_z.append(omega_z)
            pos = self.rod.position_collection
            s_nodes = np.linspace(0, rod_length, pos.shape[1])
            for i, s in enumerate(self.marker_positions):
                x = np.interp(s, s_nodes, pos[0, :])
                y = np.interp(s, s_nodes, pos[1, :])
                self.marker_x[i].append(x)
                self.marker_y[i].append(y)
            with open("ball_omega.csv", "a") as f:
                f.write(f"{time},{omega_z}\n")
            for i, s in enumerate(self.marker_positions):
                label = int(s * 100)
                with open(f"marker_{label:02d}cm_x.csv", "a") as f:
                    f.write(f"{time},{self.marker_x[i][-1]}\n")
                with open(f"marker_{label:02d}cm_y.csv", "a") as f:
                    f.write(f"{time},{self.marker_y[i][-1]}\n")


recorder = Recorder(rod, ball_rod, marker_positions, output_interval)
system.collect_diagnostics(rod).using(
    Recorder,
    rod=rod,
    ball_rod=ball_rod,
    marker_positions=marker_positions,
    output_interval=output_interval,
)

print("Simulation starting...")
system.finalize()
timestepper = PositionVerlet()

# Note: This is same as `integrate` call, just explicitly looping.
total_steps = int(simulation_time / time_step)
print("Total steps", total_steps)
time_step = simulation_time / total_steps
time = 0.0
for _ in range(total_steps):
    time = timestepper.step(system, time, time_step)

print("Simulation finished.")

Hope this helps.
Best,

[Qin-Shaowen]

Thank you for the previous help. The simulation works with k=1e2, kt=1e-3. However, in reality, the ball is tightly tied to the rubber band, meaning the connection should be very stiff, with stiffness comparable to the rubber band itself. When I increase k and kt to physically realistic values (e.g., k = E*A/L_connect), the simulation immediately produces NaN. What is the correct way in PyElastica to model a rigid attachment of a sphere to the end of a CosseratRod? Is there a way to directly attach the sphere without a flexible joint, or should I use a very high stiffness with some special damping settings to avoid numerical instability?