How to simulate the gravity of cable being applied forces?

[woshizh951]

Hi PyElasitica teams,

I want to simulate the process of cable operation. I set one end of the cable fixed, applied an upward force in the middle of the cable, and applied gravity to the whole cable.

The unfixed end of the real cable will sag naturally under the influence of gravity, as shown in above figure, but in my simulation results, the unfixed end is still lifting until the whole cable is straightened, as shown in below figure.

Here is my code:

[woshizh951]

import numpy as np

from elastica.wrappers import BaseSystemCollection, Connections, Constraints, Forcing, CallBacks
from elastica.rod.cosserat_rod import CosseratRod
from elastica.boundary_conditions import FreeRod
from elastica.timestepper.symplectic_steppers import PositionVerlet
from elastica.timestepper import integrate
from elastica import *

# Make simulator
class SutureSimulator(BaseSystemCollection, Connections, Constraints, Forcing, CallBacks):
    pass

# Make constraint to apply velocity to rod midpoint while both ends are held in place
class SutureTest(FreeRod):
    def __init__(self):
        pass

    def constrain_values(self, system, time: np.float64 = 0.0):
        # Anchor both rod ends
        system.position_collection[..., 0] = np.zeros(3)
        # system.position_collection[..., -1] = [0, 0.0, 5.0] # This position is in cm units

    def constrain_rates(self, system, time: np.float64 = 0.0):
        # # Calculate displacement from goal position (offset from midpoint by 8 mm in y direction)
        # displacement = (
        #     np.array([5.0, 0.0, 0.8]) # Goal point, in cm units
        #     - system.position_collection[..., n_elem//2] # Rod midpoint
        # )
        # # Move at low speed in direction of displacement
        # # Stop movement if control node is close enought to goal position
        # if (np.linalg.norm(displacement) > 1e-6):
        #     system.velocity_collection[..., n_elem//2] = displacement / 2 # arbitrarily scale displacement to reduce speed
        # else:
        #     system.velocity_collection[..., n_elem//2] *= 0

        # Anchor both rod ends
        system.velocity_collection[..., 0] *= 0
        # system.velocity_collection[..., -1] *= 0

class MoveMidPointByApplyingForces(NoForces):

    def __init__(self, k, ramp_up_time, target_mid_point_position):
        self.k = k
        self.ramp_up_time = ramp_up_time
        self.target_mid_point_position = target_mid_point_position

    def apply_forces(self, system, time: np.float = 0.0):
        factor = min(time/self.ramp_up_time, 1.0)

        displacement = self.target_mid_point_position -system.position_collection[...,n_elem//2]

        force = self.k * factor * displacement

        system.external_forces[...,n_elem//2] += force

class ContinuumSnakeCallBack(CallBackBaseClass):
    """
    Call back function for continuum snake
    """

    def __init__(self, step_skip: int, callback_params: dict):
        CallBackBaseClass.__init__(self)
        self.every = step_skip
        self.callback_params = callback_params

    def make_callback(self, system, time, current_step: int):

        if current_step % self.every == 0:

            self.callback_params["time"].append(time)
            self.callback_params["position"].append(
                system.position_collection.copy()
            )

            return

suture_sim = SutureSimulator()

#initialize all relevant rod and time variables
n_elem = 50
start = np.array([0.0, 0.0, 0.0])
direction = np.array([0.0, 0.0,1.0])
normal = np.array([0.0, 1.0, 0.0])
# Below values are in units of centimeters (cm) and grams (g)
base_length = 5 # cm
base_radius = 0.2 # cm
density = 3e-3 # g/cm^3
nu = 0.1 # g/s
E = 7.5e3 # g/(cm * s^2)
youngs_module = 0.03

final_time = 1
dl = base_length / n_elem
dt = 1e-4 * dl
total_steps = int(final_time / dt)

timestepper = PositionVerlet()

suture_thread = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E,
    youngs_module
)
suture_sim.append(suture_thread)

suture_sim.constrain(suture_thread).using(
    SutureTest
)

suture_sim.add_forcing_to(suture_thread).using(MoveMidPointByApplyingForces, k=8000, ramp_up_time=1.0, target_mid_point_position=np.array([0, 0.8, 2.5]))
gravitational_acc = -9.80665
suture_sim.add_forcing_to(suture_thread).using(
        GravityForces, acc_gravity=np.array([0.0, gravitational_acc, 0.0])
    )

pp_list = defaultdict(list)
suture_sim.collect_diagnostics(suture_thread).using(ContinuumSnakeCallBack, step_skip=200, callback_params=pp_list
    )

suture_sim.finalize()
integrate(timestepper, suture_sim, final_time, total_steps)

SAVE_RESULTS = True

if SAVE_RESULTS:
    import pickle

    filename = "suture_thread.dat"
    file = open(filename, "wb")
    pickle.dump(pp_list, file)
    file.close()`

Could you please kindly give me some advices to simulate the influence of gravity to the cable corrcetly?

Best Regards.

[nmnaughton]

Hi @woshizh951,

It looks like you are operating in grams/mm/sec units? Then I think you need to rescale your gravitational force. You have gravitational_acc = -9.80665, but in g/mm/s I think it should be scaled by 1e6. That would probably explain the issue then as your gravity is essentially not there so elasticity is dominating.

[woshizh951]

Hi @nmnaughton,

Thank you so much for you advice, i adjust the value of g. When i multiple g with 1e6, the position of elements all become to Nan. So i change the scale num to 1e3. And the result looks like what i want.

In addition, i have one more question searching for your help. Could Pyelstica planning the force applied on the rod with time step going on, like a robot hold the rod and do some actions. In example cases , i can only find that the force was designed like a Sinusoidal with time step process. Can force be applied more freely？

I can show an example of force i want to apply on end of rods.: 1-100step first stage a force vector [ 0 , 0 , 1] 100- 200 step second stage change the force vector to [1, 0, 0], and the rod can be pull away to another place. .... ....

Best regards.

[skim0119]

@woshizh951 Hi

I see you defined the custom force MoveMidPointByApplyingForces. When you define apply_forces method, you can use time parameter to custom define the force during the simulation.

class CustomForce(NoForces):

    def __init__(self, k, ramp_up_time, node):
        self.k = k
        self.node = node
        self.ramp_up_time = ramp_up_time

    def apply_forces(self, system, time: np.float = 0.0):
        factor = min(time/self.ramp_up_time, 1.0)

        if time > 0.5:
            force = self.k * factor * np.array([0,0,1])
        else:
            force = self.k * factor * np.array([1,0,0])
        system.external_forces[...,self.node] += force

In my experience, giving sudden changes in force can be unstable, so I would recommend using a smooth transition of force.

Let me know if you have any question.

[nmnaughton]

@woshizh951 One other quick comment. Regarding the stability when you properly scale gravity, you might want to check your material properties. I noticed you have E = 7.5e3 Pa. That is very low. Are you sure that is the value you want? I messed around with your code sample and found E = 7.5e6 Pa along with a dt 10x smaller was stable for a proper gravity value (note that I also scaled the k value in your applied force by 1e3 too since gravity went up). If you really are interested in modeling a suture, I think your young's modulus will probably be in the 1e6-1e9 range.

[woshizh951]

@nmnaughton
Hi, Thanks for your comment. I'm trying to simualte the process of manipulation normal cable such as network cable, the parameters was initialized from another examples. I would modify the E and k follow you advices. And i will adjust the physical parameters follow a tested cable in one paper.

I noticed that the shear_module, which can not find more information in code of pyelastic, can be set when generating a straight rod, but how can i set the bending young's module? I found only bend_matrix in the build-in parameter of the function, could you please tell me how can i set all the parameter i wanted or where can i find the document/reference about bend_module?

Best regards.

[armantekinalp]

Hi @woshizh951

You can set stretching Young's modulus and shear modulus while initializing the rod, as you can see from the below code snippet. But we don't have functionality to set Bending Young's modulus separately. However you can modify bend_matrix after initializing the rod. For example.

suture_thread = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E=E_stretch,
    shear_modulus=shear_modulus,
)
I = np.pi/4 * base_radius**4
suture_thread.bend_matrix[0,0,:] = I * E_bending
suture_thread.bend_matrix[1,1,:] = I * E_bending
suture_sim.append(suture_thread)

[woshizh951]

Hi @armantekinalp,

Thank you. I have tested the bend matrix and ensured that it worked. However, when I set the parameter of straight_rod like your example, an error occurred and said the youngs module is missing. So I just set the parameter like this:

suture_thread = CosseratRod.straight_rod(
    n_elem,
    start,
    direction,
    normal,
    base_length,
    base_radius,
    density,
    nu,
    E_st,
    shear_module = shear
)
I = np.pi/4 * base_radius**4
suture_thread.bend_matrix[0,0,:] = I * E_bending
suture_thread.bend_matrix[1,1,:] = I * E_bending
`suture_sim.append(suture_thread)

As I know, the basic definition of young's module is the same as the stretching young's module, so the way I set is also right?

[armantekinalp]

Hi @woshizh951

Yes how you set the rod is correct.

If you don't have any other questions I will close the issue.

[woshizh951]

Hi @armantekinalp

I have no question right now. Thanks to all members of Pyelastic team help to me with patience.

Best regrads,

[skim0119]

I'm closing the issue. @woshizh951 Feel free to make another issue if you have any other question.