simtk.org

Hello List,

I am looking for some guidance in creating a new friction-force model. I intend to create a variant
of the ElasticFoundationForce that does what I need.

Friction occurs between two bodies. In my usecase one of the bodies is special in that the friction depends upon
the direction of the force on that body.

Let me describe it in more detail. I try to simulate a rower in a rowing boat using OpenSim. The goal is to find the optimal way
to rigger the boat (place of lock, lenght of oars) and move the rower (leg, back, arm coordination) in order to get
the maximum speed, given the capabilities of the rower (size, strength of muscles).
Eventually the rower should be a full body musculoskeletal model, but for now it is kept simple.

The interaction with the water is greatly simplified, the boat is just a big box lying on the ground, and friction
is used to model the resistance to moving forward through the water.
The blade is modelled similar, but this is more complicated. And therein lies my ploblem. This is my current model.
Note that normally a blade is vertical to work properly in the water, but here it is a flat plate with friction to the ground.
The next picture is a top view of the situation. The boat is moving forward in the water, see arrow 1.
The rower is exerting a force (2) so, via the lock, a force (3) is exerted on the blade.
Because the boat is moving there also is a force (4). A moment later the boat and lock have moved a bit so the oar is in
a new position, see the red line (5). The important thing to note is that the blade is resisted very much to move in direction 3, but very little in direction 4.
This of course because of the distinctive shape of the blade.
So if we want to model this by using friction in the plate from our model we need high friction coefficients
in one direction and low friction coefficient in the other.

I hope it is feasable to create a friction in Simbody that implements this, using ElasticFoundationForce as a starting point.

There, forces are calculated in functions calcForce and processContact.
To calculate the resultant force at a certain moment, I need to know the direction of the resultant force
on the blade (arrow 3 and 4 added) with repect to the orientation of the blade.

Is the resultant force one of the forces from bodyForces? Which one? What is the index in this vector?
How can I find the orientation of the body "blade" that has this special behavior?

Hopefully someone can give me a pointer in the right direction. Any help would be much appreciated!

Sietse

To answer my own question: it was not too difficult in the end and I did the following.

• I created a clone of ElasticFoundationFource.cpp and *Impl.h.
  Then changed the function processContact as follows.

• First calculated the angle in the XZ plane between de X-axis and the direction of the speed:

  Code: Select all

     const Rotation body1_rot = body1.getBodyRotation(state);
     const Vec3 xdir = body1_rot*Vec3(1,0,0);
     Vec3 vel = body1.getBodyOriginVelocity(state);
     Real angle;
      if (vel.norm() < 0.0001) {
        angle = 0;                                                                                                                                     
      } else {
        // dot product in the XZ plane                                                                                                                                                                  
        const Real dotprod = vel[0]*xdir[0] + vel[2]*xdir[2];
        const Real prod = std::sqrt(vel[0]*vel[0]+vel[2]*vel[2]) * std::sqrt(xdir[0]*xdir[0]+xdir[2]*xdir[2]);
        angle = std::acos(dotprod/prod);          
        if (angle > Pi/2) angle = Pi - angle;
      }
  

• And finally in the loop over all the springs I use angle to select between the parameters given to the Force-object and an internal set of values that represent the friction in the direction of the blade. Here for the stiffness.

  Code: Select all

      const Real stblade = 1e9;
      Real stiffness = stblade*std::cos(angle) + param.stiffness*std::sin(angle);
           ......
  

  And use there values in calculating the forces.

This works fine when the otherObject is a halfSpace "on the ground". But there is where my "water" is anyway. Also the Mesh has to have a flat surface in contact with the otherObject.

Nice! Congratulations.