Force dependent kinematics with STL surfaces

Main Forum
1

Hi guys,

I’m working with the new functionnality of the version 5.1, which puts together the Force dependent kinematics algorithm with contact simulation based on STL surfaces. In the tutorials however, you don’t specify how the direction of the contact forces is defined. In my simulation (humerus on glenoid contact with a stiffness of 5e7 N/m), I can see that some forces are not oriented towards the centre of the joint. Thus, does these forces are still correct ? Can I improve this by, for example, refine the mesh ?

Thank you.

Pierre

2

Hi Pierre,

the forces calculated in the AnyForceSurfaceContact are calculaed as weighted averages of the surface normals. That means that the direction of the forces applied by this object (which are Fmaster and Fslave) (…and also the moments) depend on the STL surfaces you are using.
So yes, you are right: if you would like to have a more spherical behaviour for your joint you could use some external tool to change your mesh to better fit to a sphere.
The refinement options in the AnyForceSurfaceContact class (MeshRefinementMaster and MeshRefinementSlave) are refining the STL meshes by cutting the triangles, but it will not change the normals of the mesh or smoothen the surface.

Best regards,
Daniel

3

Hi Daniel,

Thank you for this answer. I will try to refine the mesh.

I also have another problem. When I do the simulation with a contact between the scapula and the humerus with the AnyForceSurfaceContact, there is a big penetration between the surfaces (see the attached image called penetration.png).

Can you explain it ? Is it caused by the poor quality of the humerus and scapula STL bones, provided by the AnyBody shoulder model ?

Thank you.

Pierre

4

Hi Pierre,
one thing is that the bones in the Repository are rather course surface meshes, so that could cause some problems using them for the contact calculation. Remeshing and smoothing of the surface could help in this case.

Another way to decrease the penetration depth is to increase the PressureModule (which is a value in N / m^3).

Best regards,
Daniel

5

Hi Daniel,

Ok I get it. I refined my mesh in the contact area of both STLs and now it works much better.

Does the pressure modulus is the Bulk modulus ? For me, it looks like it. If it is, I think you should reconsider and increase its value for the prostheses in the femur-tibia tutorial. I found that a value of 160 GPa should be more appropriate for material like steel. I also found a value of 12.5 GPa for bone. Do you know if these values are correct ?

I also noticed a bug when re-playing an analysis with these contact STLs. In fact, plotting the contact forces at each frame does not work. It stays frozen on the last contact frame, while the bones are moving. Is this normal ?

Thank you.

Pierre

6

Hi Pierre,

for calculating our pressure module, we oriented at the publication of Y. Bei and B.J. Fregly: Multibody dynamic simulation of knee contact mechanics. Here is a calculation of a pressure using a simple spring model as a function of the penetration. You can do a similar thing to get an idea of the pressure module.
When you calculate values for your model, you should consider that there is very often a cartilage layer between two bones. It depends a bit on the shape of your surfaces, but I think you might have difficulties with the FDK solver if you choose a too large pressure module. In generel it is a good approach to start with a smaller value and increase it until the kinematics look ok.

The problem with the vector field in the replay is a know issue and it is only a tool that should help to debug a model. If you want to replay it, you have to capture the simulation e.g. using an AnyCamera. The vector field is also not saved in the output of the simulation.

Best regards
Daniel

1 Like
7

Hi Daniel,

The paper from Bei and Fregly should help us to understand and eventually explain this STL contact method if we publish a paper in which we use it.

Another problem that I have happens when I use a mesh that is too much refined, in which case I get forces at nodes that do not touch anything. In this case, the simulation does not converge and I get a lot of “Degenerate triangle #…” errors.

Finally, since prostheses ofter produce friction, it would be nice to include a coulomb-based friction force at each penetrated vertex. To do so, you would rely on a measure such as:

AnyKinRotational GHRot = {
AnyRefNode &ref1 = …Seg.Scapula.gh;
AnyRefNode &ref2 = …Seg.Humerus.gh.;
Type = RotVector;
};

This would allow to determine the instantaneous direction of rotation and the friction direction would simply be the reverse (minus) of this direction. The magnitude of friction would thus be the multiplication of the contact force by a friction coefficient. Thereafter, you would just have to add a friction coefficient input variable in the STL function.

I think that this would be a major improvement, especially when simulating hemi-prosthesis contacts, for which one of the two bones is not replaced.

Best regards.

Pierre

8

Hi Pierre,

We also saw the problem that the amount of messages of degenerate triangels are a little bit annoying. We will do something about it in the next update.

The reason is that by converting STL files to AnySurf files, vertices which are close to each other are collapsed into a single vertex. This has some positive effects on drawing the surfaces, unfortunately it can change the topology of the surface mesh a little bit.
To avoid this, you should coarsen the stl a little bit, maybe it might be enough to coarsen sharp and small details.

About your idea of friction I discussed it with Michael:

AnyBody’s contact does currently not include a friction component and we do not think it will be easy to implement it yourself in AnyScript since it does not offer vertex based output to work with.
One option would be for you to do your own contact forces from scratch using the Python or C++ hook (see e.g. the example of AnyFunEx in the reference manual). But maybe this will be very challenging.

There is a reason why AnyBody’s contact does not have friction. It could quite easily be included in the forces, but AnyBody does not have a solver that deals well with friction whether Coulomb-based or more viscous. History/path dependence of kinematics, which arises with friction, is not handled by AnyBody’s current solvers. The standard AnyBody solver is pure inverse dynamics requiring you to specify all motions and the force dependent kinematics (FDK) solver handles indeterminacy of selected degrees of freedom by assuming a kind of quasi static equilibrium.
To handle the real history/path-dependence introduced by friction would require a more forward dynamics like solution, which would be much more difficult to work with.

We could recommend you to consider transferring the AnyBody predicted forces into an FE model of the implant where you would be able to make a more detailed contact model capturing local phenomena.
Unless you of course believe that friction will influence the kinematics and thereby muscle moment arms significantly. In that case we will have problems.

Regards,
Michael and Daniel

9

I have to short questions, which are related to a similar problem. So I thought I post it under the same thread.

I am trying to setup a simplified model of a TKA using Surface-to-surface contact and FDK.
When it comes to modeling Metal on PE contact, I would like to model a Penetration-Pressure function that is not linear. However, in the current version of the AnyForceSurfaceContact it is only possible to specify a linear PressureModule. Is there a way to specify another non-linear function (such as an Interpolation function)?

I also wanted to play around with external python functions in order to specify a function that has several inputs and outputs. To get an idea how to implement the function I have loaded the demo file AnyFunEx.any, however I am always getting the following error:

ERROR(SCR.EXP.FUN.EXT4) : C:/P…)/A…y/A…1/D…n/A…e/demo/C…s/A…x/AnyFunEx.any : doit : AnyBody.Python extention module is not installed or does not match installed version of Python.

Any idea how to solve that?

Many thanks for you help

Christine

10

Hi Christine,

unfortunately, there is no possibility to use another function for the force law of the contact up to now. But we already thought about having it in the modeling system and there is a good chance that it will come with one of the next updates.

About the Python hook AnyFunEx:
It sounds like there is a Python installation missing on your computer. To use the Python hook, Python version 2.6 or 2.7 which you can find e.g. on python.org needs to be installed on the computer. Please make sure to install the 32 bit version since the 64 bit Version is not supported in the current version of AnyBody.

Best regards
Daniel

11

Hi Daniel

thanks for the quick reply.
I will look out for new releases and will work with the linear relationship in the meantime.

Your solution with Python has worked out. I had Python 3.1 installed on my computer.

Thanks
Christine

12

Hi Daniel,

My colleague Lauranne (laul.sins) and me are working on a shoulder model involving a total anatomical prosthesis (glenoid and humerus components). We are using the FDK algorithm to control the translation of the humeral head with the following functions:

AnyKinLinear GHLin = {
AnyRefNode &scapula_gh = …Seg.Scapula.gh.ghProth;
AnyRefNode &humerus_gh = …Seg.Humerus.gh;
};

AnyKinEqSimpleDriver GHLinCon = {
AnyKinMeasure &ref = .GHLin;
DriverPos = {0.0015,0.004,0.00};
DriverVel = {0,0,0};
Reaction.Type={Off,Off,Off};
CType = {ForceDep,ForceDep,ForceDep};
};

AnyForce TranslationalGHStiffness = {
AnyKinLinear &ref = .GHLin;
F=-7000*{0,(ref.Pos[1]-.GHLinCon.DriverPos[1]),(ref.Pos[2]-.GHLinCon.DriverPos[2])};
};

AnyForceSurfaceContact ProtheseContact = {
AnySurfSTL &master = …Seg.Humerus.HumProthNode.ProtHumSTL;
AnySurfSTL &slave = …Seg.Scapula.GlenPos.GleneSTL;
PressureModule = 5.5e8;
};

The “ghProth” node on the scapula has the orientation of the glenoid component (major axis = y; minor axis = x).

We chose the humeral head component as the master surface, since it is made of an harder material (titanium) compared to the glenoid component (plastic).

The pressure module and the translational stiffness were found by iteration, by looking at the level of penetration and the converge of the solution.

Moreover, in the AnyBodyStudy function, we put:

...

tEnd = 8.0;
Gravity = {0.01, -9.81, 0.01};
nStep = 100;

InverseDynamics.ForceDepKinOnOff=On;
InverseDynamics.Criterion.Type=MR_MinMaxAux;
InverseDynamics.Criterion.AuxQuadraticTerm.Weight = 0.05;    
InverseDynamics.ForceDepKin.MaxIteration = 50;

…

The STL surfaces (in mm) are enclosed with the thread. Is the mesh pattern and size seem correct ?

Thank you.

Pierre

13

Hi,

What means the following error when using the FDK:

Too small steps for GSS search causes looping of Newton method

GSS is an acronym for which term ?

Thank you.

Pierre

14

Hi Daniel,

Does the pressure module is a well-known parameter ? I didn’t find any information on the internet regarding this parameter.

Also, in the paper that you mention, there is nothing related to this pressure module in N/ m^3. Now I must find a value that is representative of UHMWPE, which is the plastic used for the glenoid prosthesis. I already tried a pressure module of 5e8 and it seems to produce an unrealistic high penetration. I also tried a value of 5e9, which produce less penetration but the model becomes highly unstable, which indicates that the joint surfaces are too hard and do not penetrate enough in each others.

Thank you.

Pierre

Also, does the convergence of the FDK solution is sensitive to the approximation of the initial joint position ?

15

Hi Pierre,

GSS stands for golden section search, a well know simple optimization algorithm based on the idea of cutting intervals into pieces.

I think the name pressure module is not too commen. But the given paper describes how to get a value for this parameter.

Since the FDK solver is an iterative solver the convergence is sensitive to the initial guess. So if your joint positions are part of the system to solve with FDK, the solution is sensitve to the initial joint position.

Best regards
Daniel

16

Hi Daniel,

It is important to know how the STL contact force is computed along with the FDK algorithm, in order to explain it well in our fortcoming publications.

In the paper that you mention (Bei et al. 2004), there are different equations for different types of contact (linear vs. non-linear). Which of the four equations is used ? You said that the pressure module is in N/m^3. From equation 1 of Bei et al. 2004, the only way to obtain this unit is to divide the aggregate modulus (N/m^2) by the layer thickness “h” (m) ? However, it does not make sense because the thickness is computed on a element-by-element basis.

I searched for the GSS and found some information on the internet. I guess that the value to minimize is the residual force of the DOFs controlled by the FDK ? How the FDK handles the various degrees of freedom that it controls ? Does it handles them separately ? With this method, what does means the too small steps that cause looping of the Newton method ? Does it means that the algorithm reaches the “MaxNewtonStep” of 0.001 without finding the residual force below the “ForceTol” of 0.001 N ? If yes, does the residual force given in the error message is the force computed in the last step ? Also with this method, what happens if the maximum number of iteration defined by the “MaxIteration” parameter is reached without finding a solution ? Does the position of the lowest residual force is chosen and the corresponding residual force is given in the error message ? Finally, what does means the perturbation parameter (1e-5) ?

17

Hi Pierre,

which equation you use is up to you and depends on the material you want to simulate. For everything I did I assumed a linear elastic behaviour of the material, so I used the first equation with sum assumption of the layer thickness.
But I would like to remind again that the contact algorithm is not meant to replace any FE analysis. It’s main purpose is to get the kinematics right to do calculations with rigid bodies.

Regarding the GSS: as stated in the error message a Newton method is used to solve the FDK. The GSS is just used to calculate a step size for the Newton algorithm to increase the area of convergence and get a higher stability.

Best regards
Daniel

18

Hi Daniel,

What do you mean by “sum assumptions of the layer thickness” ? Did you consider a constant thickness rather than an element-by-element thickness ? If yes, could you tell me which value you used so I would be able to justify the pressure module that I used.

Thank you.

Pierre

19

Hi, Pierre,

sorry, the sum was a type, it should of course be some. What do you mean by element-to-element thickness? There are no elements in AnyBody as long as you don’t use the Any3FE converter.
As far as I remember we choose values between 1 and 10 cm for the calculation of the pressure module to get resonable results, regarding the penetration.

Best regards
Daniel

20

Hi Daniel,

Our glenoid implant is about 5 mm thick. With UHMWPE values (E = 0.8 - 1.6, v = 0.46) of the literature, I found a minimum aggregate modulus of about 3.7e9 N/m^2. If I want to find the pressure module from that value, do I have to divide it by the implant thickness ? If yes, it gives a really high pressure module of 7.4e11 N/m^3 (3.7e9 N/m^2 / 0.005m).

If this is the right value, I do not obtain a stable joint, since this pressure module produces a really low penetration in the glenoid. That brings me to my next question: Do you think that it would be possible to re-create a constraint that would allow a connection between the humeral head translation and the activity of the rotator cuff such as the “GHReactions.any” ? I think that with such constraint, the additional rotator cuff activation would allow a slower translation of the humeral head on the glenoid surface. From my opinion, it will be challenging if not impossible to recreate such constraint, since the contact force is computed with the FDK, which itself depends on the muscle recruitment.

What do you think ?

Pierre