Equations of Motion for GaitLowerExtremity

basvdploeg · December 14, 2010, 12:03pm

Hi,

As part of my Master Assignment within the TLEMsave project of Prof. Verdonschot and Prof. Koopman at the university of Twente, I am trying to validate the Equations of Motion for an Anybody model. I am interested in the femur of GaitLowerExtremity model in the ISB coordinate system. Since this produces large output files, I am using Matlab to process the data.

I have written a code which extracts all insertion positions, forces, moments and gravity vectors per timestep from the AnyMechOutputFileForceExport file in a local reference frame (ISB coordinates)

The input of the Newton Euler equations (on top of every timestep in the AnyMechOutputFile) are the Jmatrix, rCoM, Axes and the accelerations. Theoretically, the resulting Forces and Moments should match the summed Forces and Moments of the AnyMechOutputFile.

Forces are: sum of all forces in x, y or z direction, plus the Gravityfactor in this direction times the mass of the segment.

Moments are: sum of all moments in x, y, or z direction, plus the outer product of the forces times their positions, plus the outer product of the gravity factors times their distance to the centre of mass from the point of rotation.

!!! I attached a file with the formulas I used for clarity !!!

My questions is as follows:

I tried to validate the formulas I use with a small model first. I used a single segment model similar to the one which is used in the tutorials. For this model, the forces and the moments were exactly the same as calculated by the Newton Euler equations! For the complete femur model, only the forces are exactly the same. I think this has to do with the fact that the rotation point of the hip (were the local reference frame is {0,0,0}) has relative motion and rotation with respect to the global reference frame. I think I need to use the Axes to rotate the Jmatrix, but everything I tried so far did not work. (I also extensively searched through the forum and the reference manual) Does somebody know what needs to be multiplied with which rotation matrix to achieve the right moments? (forces are already correct). Are the accelerations with respect to the global reference frame or with respect to the pelvis?

Thanks in advance,

Bas van der Ploeg
Twente University

toerholm · December 17, 2010, 10:33am

Hi Bas,

Sorry for the slow feedback you have not been forgotten, these are good questions.

We will look into this and give you some feedback next week but I can give some basic advices right away though.

The accelerations measured on a segment is always wrt global.

If was you i would continue with the simple example and make it more complex to reflect the situation in the leg model.

So move the origin of the segment coordinate system away from the center of mass and make use of the rCom property, also move the center of rotation away from global origin.

Best regards
Søren

toerholm · December 20, 2010, 11:24am

Hi Bas,

I have got some help from our developers, here is a more detailed answer that i hope will help you:

Your equations do indeed look right, but there may be a few of the quantities that are not completely clear (I acknowledge the documentation is not very extensive on this file format).

Jmatrix in the file (AnyMechOutputFileForceExport) is “The actual inertia matrix referring the body-fixed segment coordinate system.”, see the AnyScript Reference Manual for AnySeg because it is direct output of this one.
This implies that it must be what your term: (Jc - m[c][c]) and it is local coordinates (which makes it constant over time. AnyBody internally uses local version of the rotational equation for this reason, but of course it should add up in any system of reference.)
This also implies that for the rotational equations you must either transform all torques into the local frame or the inertia terms out into the global frame for instance JMatrixGlobal = AJMatrixA^T.
rCoM is the actual position of CoM (global or in the reference frame chosen for the file). This is implies it is not your c-vector.
In AnyBody the r-vector are always “global” position whereas ‘s’ is typically used for relative vector so you have relationships like

r_p = r_0 + A_0 * s_p’ (where the mark on s_p indicate a vector of local components)

So what you need for c is what we could call sCoM (in AnySeg you find this, with the description “Relative location of the center of mass (CoM).”

sCoM is not in the file, but I guess you could find it as sCoM = (rCoM – r) and to get the local components you need to transform it into the local frame

I hope, I captured everything that could be the course of your difficulties.

Best regards
Søren

basvdploeg · December 21, 2010, 11:35am

Hi Søren,

Thanks for the reply. I think it is really helpful.

I also extended my simple model like you advised earlier on friday. I attached a picture of it. There are two segments, with a local reference frame on the joint between the two segments. In this case, the forces AND the moments are correct! So, I think the troubles with the full model should have to do with the rotation and/or the translation of the hipjoint from the global reference frame.

Thanks again for taking time for me. I will let you know whether I managed to solve the problem in the coming weeks.

Happy Christmas,

Bas

basvdploeg · January 10, 2011, 11:58am

Hi Søren,

Your help provided some new insight in the way the Axes work in Anybody. But unfortunately, the results were the same: The summed forces matched the Newton Euler equations, but the torques did not.

Since my simple leg model did provide the right answers, I started to make my model a bit more complicated to replicate the error.

First I made the model overdetermined by adding a couple of muscles to the system and I also added more degrees of freedom to the hipjoint. Errors only occurred when a wrapping surface was added to the segment. Suddenly there was a residual torque in the outputfile. Closely examining the outputfile resulted in the following: The effect of the wrapping surface is modelled by Forces AND Moments in the Cylinder Centre, which I think is correct. But the exact same Moments are given at the output of the Insertion Node of the muscle, which I think is not correct. This means that when you sum all the moments acting on a segment, the effect of the wrapping surface is added twice. Removing the moments at the insertion node of a muscle with a wrapping surface provided the exact answers corresponding to the Newton Euler equations for the Forces AND the moments.

Force[13]:

Name=Main.Studies.HumanModel.BodyModel.Right.Leg.Mus.GastrocnemiusMedialis1;
Class=AnyMainFolder.AnyFolder.AnyFolder.AnyFolder.AnyFolder.AnyFolder.AnyFolder.AnyShortestPathMuscle;
SegmentName=Main.Studies.HumanModel.BodyModel.Right.Leg.Seg.Thigh;
SegmentID=0;
Pos=1.516149572829147e-002;-2.236731107064947e-001;-8.450775980200284e-003;
RefFrame=Main.Studies.HumanModel.BodyModel.Right.Leg.Seg.Thigh.GastrocnemiusMedialis1Node;
Components(1):
F[0]=-2.877295183599379e+001;-5.332407340307649e+001;-9.087997427162817e+000;
M[0]=4.390214157017099e-002;1.764555138437333e-001;2.506708675195163e-002;

Force[14]:
Name=Main.Studies.HumanModel.BodyModel.Right.Leg.Mus.GastrocnemiusMedialis1;
Class=AnyMainFolder.AnyFolder.AnyFolder.AnyFolder.AnyFolder.AnyFolder.AnyFolder.AnyShortestPathMuscle;
SegmentName=Main.Studies.HumanModel.BodyModel.Right.Leg.Seg.Thigh;
SegmentID=0;
Pos=2.497514693192334e-002;-2.610040272586361e-001;-1.073897890306629e-002;
RefFrame=Main.Studies.HumanModel.BodyModel.Right.Leg.Seg.Thigh.CylCenter;
Components(1):
F[0]=2.451557698298741e+001;-7.610008071970359e+000;9.437994667041821e+000;
M[0]=4.390214157017099e-002;1.764555138437333e-001;2.506708675195163e-002;

The right answers were found when I replaced the red line with: M[0]=0.00000000000e+000;0.00000000000e+000;0.00000000000e+000;

This was exactly true for both my simple model and the GaitLowerExtremityModel (2 wrapping surface muscles: Gastocnemius Lateralis and Medialis)

This means that I think there is an error in the way the system produces the AnyMechanicalOutputFile. Please correct me if I am wrong. I can provide more details about the way I calculated the forces and moments if needed, although the formulas I used are the same as I attached to this thread before.

Best regards,

Bas

toerholm · January 11, 2011, 1:42pm

Hi Bas,

Thanks a lot for making us aware of this, this indeed looks like a bug in the output

We have been able to reproduce the problem and will try to look into it as soon as we can.

Thanks again.

Best regards
Søren

basvdploeg · February 10, 2011, 1:54pm

Hi Soeren,

If it indeed is a bug, can you please let me know? For my research it is not really a problem if the software provides different answers, I can easily remove the redundent moments since I am using Matlab for the output files. But it is important to know whether I am allowed to remove them, for validation purposes. In other words: Can you perhaps confirm it is wrong?

Thanks in advance,

Best Regards,

Bas

toerholm · February 11, 2011, 10:27am

Hi Bas,

It is a bug sorry for not being more clear on this, previously.

We have fix it and it will be part of the next release of AMS.

Best regards
Søren

basvdploeg · May 2, 2011, 1:33pm

Hi Soeren,

I may have found a different error in the outputfile. I may also be that I am not looking careful enough. I hope someone can help.

I am using the ISBCoordinateSystem as a local reference frame for the output in the AnyMechOutputFileForceExport.

The data is used for FEA. I am using Anybody v5.0 and I imported c3d data from S.Kolk. I am interested in the femur,
the initial position of the mannequin is seated and the gravity is pointing in the -z direction. I have watched the webcast
“Streamlining gait Analysis with the Anybody Modeling System v.5”.

In the TrialSpecificData.any folder I set the Gravity to: AnyVec3 Gravity= {0.0, 0.0, -9.81};
I follow all project steps to obtain the outputfile w.r.t the local reference frame. The output file produces the Gravity vector
in the local reference frame and all forces and moments. All insertionpoint positions have the correct coordinates, so I am pretty sure
I use the local reference frame correclty in the outputfile.

If I apply the given Gravity vector times the mass of the segment in the FE model, the model is in “dynamic equilibrium” as expected
(i.e. The reaction forces in the FE software match the reaction forces supplied by AnyMechOutputFileForceExport).
However, If I closely inspect the gravity vector, it is not pointing in the correct direction. If I make a small calculation, it is
more visible.

[Axes][Arel]{Gravity}={0, -9.81 ,0}, whereas you should expect it to be {0, 0, -9,81} in the global reference frame.

In this calculation for the first time increment:

Axes = found at the top of each time sample of the outputfile
= -5.030132066517148e-001,-2.518582194008671e-001,8.267678944263521e-001;3.070927961375282e-001,-9.462599016190475e-001,-1.014209699635392e-001;8.078810113465929e-001,2.028783771330381e-001,5.533250722653647e-001;

Arel = …Seg.Thigh.ISBCoordinateSystem.Arel
= {{0.8836248, -0.1270202, -0.4506364}, {0.1262953, 0.9914821, -0.03182299}, {0.4508401, -0.02879368, 0.8921402}}

Gravity = found at the top of each time sample of the outputfile
=-1.041056459517467e+000;9.557750366246992e+000;1.949797267564597e+000;

I think, the software does not rotate the gravity vector correclty, and subsenquently may provide wrong forces and moments. Or do I have to use another rotation matrix to obtain the Gravity vector in the local reference frame? I hope that it is solved easily.

Best Regards,

Bas van der Ploeg

basvdploeg · May 4, 2011, 12:42pm

Hi Soeren,

Together with the guys from the TLEMsave project, we were able to solve the problem. The gravity vector is defined in the TrialSpecificFolder (LoadC3DMarker.Main.Any). However, there exists no reference between the TrialSpecific Gravity vector and the one that is used in the InverseDynamicsStudy folder. This means that, without adjusting the gravity vector in the InverseDynamicsStudy, the vector is re-defined. To avoid using different vectors, the manual gravity statement in the InverseDynamicsStudy (InverseModel.Main.any) should be changed:

  InverseDynamics.Criterion.UpperBoundOnOff=Off;
  nStep=Main.ModelSetup.nStep;
  Gravity ={0.0, -9.81, 0};
  tStart =Main.TrialSpecificData.tStart+2*Kinematics.ApproxVelAccPerturb; 
  tEnd = Main.TrialSpecificData.tEnd-2*Kinematics.ApproxVelAccPerturb;
  //#include "GCDOutput.any"
};


  InverseDynamics.Criterion.UpperBoundOnOff=Off;
  nStep=Main.ModelSetup.nStep;
  Gravity = Main.TrialSpecificData.Gravity;
  tStart = Main.TrialSpecificData.tStart+2*Kinematics.ApproxVelAccPerturb; 
  tEnd = Main.TrialSpecificData.tEnd-2*Kinematics.ApproxVelAccPerturb;
  //#include "GCDOutput.any"
};

Best Regards,

Bas van der Ploeg

toerholm · May 6, 2011, 1:01pm

Hi Bas,

Great the problem is resolved, you are right about the gravity definition in this model, this has been changed now. Thanks for pointing this out!

Best regards,
Søren