Follower Forces and Torques

[SvRichter]

Hi,

I am trying to implement the helical motion benchmark from this paper:

https://www.researchgate.net/publication/329252597_A_Geometrically_Exact_Model_for_Soft_Continuum_Robots_The_Finite_Element_Deformation_Space_Formulation

To do this I need to define a follower force and follower torque at the tip of the rod. The force needs to remain perpendicular to the rod under deformation. This is my current attempt at defining a follower force.

class endpointTorque(NoForces):
    def __init__(
        self,
        torque=np.array([0.0, 0.0, 0.0]),
        rampUpTime = 0.0
    ):  
        self.torque=torque
        self.rampUpTime=rampUpTime

    def apply_torques(self, system, time: np.float64 = 0.0):
        if self.rampUpTime != 0:
            factor = min(1, time / self.rampUpTime)
            system.external_torques[..., -1] += system.director_collection[...,-1] @ (self.torque*factor)
        else:
            system.external_torques[..., -1] += system.director_collection[...,-1] @ self.torque



class endpointForce(NoForces):
    def __init__(
        self,
        force=np.array([0.0, 0.0, 0.0]),
        rampUpTime = 0.0
    ):  
        self.force=force
        self.rampUpTime = rampUpTime

    def apply_forces(self, system, time: np.float64 = 0.0):
        if self.rampUpTime != 0:
            factor = min(1, time / self.rampUpTime)
            system.external_forces[..., -1] += system.director_collection[...,-1] @ (self.force*factor)
        else:
            system.external_forces[..., -1] += system.director_collection[...,-1] @ self.force

Unfortunately, this doesn't result in the desired helical twisting of the rod. Instead, the rod is stretched in one direction.

Instead, I would expect something like this:

Do you know what I am doing wrong? Thank you very much.

[armantekinalp]

Hi @SvRichter

This is very interesting benchmark, I am not sure the forcing function you wrote is the same as the one given in paper. In the paper they are saying S-shaped profile. I haven't read paper in detail, but do they mention the profile somewhere else?

In you implementation the force seems constant and not changing, I think force direction has to be changing.

[SvRichter]

Hi @armantekinalp sorry for the late reply. I would need a force that remains constant relative to the rod tip no matter how the rod is deformed.

[fepauly]

Hey @SvRichter ,
have you found a way to implement the follower force?

best regards

[skim0119]

Hi @armantekinalp sorry for the late reply. I would need a force that remains constant relative to the rod tip no matter how the rod is deformed.

@SvRichter @fepauly
In PyElastica, force is defined in Langrangian coordinates, and torque is defined in Eulerian coordinates. Your implementation for force seems correct that you transformed coordinate by multiplying director system.director_collection[...,-1] @ (self.force*factor). For torque, this transformation is not needed.

Please try as below:

class EndpointTorque(NoForces):
    def __init__(
        self,
        torque=np.array([0.0, 0.0, 0.0]),
        rampUpTime = 0.0
    ):  
        self.torque=torque
        self.rampUpTime=rampUpTime

    def apply_torques(self, system, time: np.float64 = 0.0):
        if self.rampUpTime != 0:
            factor = min(1, time / self.rampUpTime)
-            system.external_torques[..., -1] += system.director_collection[...,-1] @ (self.torque*factor)
+.           system.external_torques[..., -1] += self.torque * factor
        else:
-            system.external_torques[..., -1] += system.director_collection[...,-1] @ self.torque
+.           system.external_torques[..., -1] += self.torque



class EndpointForce(NoForces):
    def __init__(
        self,
        force=np.array([0.0, 0.0, 0.0]),
        rampUpTime = 0.0
    ):  
        self.force=force
        self.rampUpTime = rampUpTime

    def apply_forces(self, system, time: np.float64 = 0.0):
        if self.rampUpTime != 0:
            factor = min(1, time / self.rampUpTime)
            system.external_forces[..., -1] += system.director_collection[...,-1] @ (self.force*factor)
        else:
            system.external_forces[..., -1] += system.director_collection[...,-1] @ self.force

[skim0119]

Hi -- I'm closing this issue. If you have further question, feel free to re-open.

[gabotuzl]

Hello,
I have a question regarding :

For torque, this transformation is not needed.

I am confused as to why all the external_forces classes that apply torque, use _batch_matvec(system.director_collection, torque) instead of directly broadcasting the batch torque vector into system.external_forces. Trying out both of these brings about different results, so I would like to ask which one is the correct one and why.