for the thickness problem: I would recommend to look at the results. If the penetration depth looks ok and everything works fine then use it. If not I would try to reduce the number to optain a more stable solution. For me it doesn’t make sense if you choose a very high Pressure Module which gives you the correct penetration but the FDK solver stops with a rather large error (and thus, the position in all other directions may not be solved accurate).
If you think it might help, you can also try to add more constraints. I’m sorry that I can’t help you here since I haven’t looked at the GHreactions.any up to now.
You are right for the lower pressure module. So far, this is the only way we found to converge to a solution.
From your opinion, what would be the maximal acceptable residual force error in Newtons (0.1, 1, 10 N) ? Maybe it would make more sense to compute this maximal residual force as a percentage of the maximum contact force or another relevant force estimated from the model ?
Adding more constraints is not really possible for now. We already use a translational spring which, from our experience, has an optimal stiffness value for the simulated conditions. This value is not far from that used in your knee implant tutorial.
Now the challenge is to be able to justify all these choices. This is always a difficult thing to do during abstract/article writing process and this is the main reason why we communicate with your support group so often.
Hi Daniel,
I’m currently working on the shoulder model with Pierre. I’m looking for validating our model and especially our choice of module pressure. I’m not sure to see exactly the link between the formula described in the paper of [Bei and al, 2004] and the AnyForceSurfaceContact algorithm.
In the tutorials (FDK : lesson 4), you describe the pressure module as ''a constant defining a linear law between penetration volume and force". Does it correspond to the right side of the equation of the paper (parameter “p”), or does it correspond to the material property divided by the thickness of elastic layer, or… ? What is the unit of this pressure module ?
the paper that is mentioned describes how to calculate a pressure p from material properties as the Poisson ratio, the Young’s modulus and a layer thickness h for a deformation d, which corresponds to the penetration in our case. So for a certain area, you’ll get the penetration force by multiplying the pressure by this area. And that’s basically the link to the contact force.
Looking at the units, the pressure module as to be in Pa/m or N/m^3 to give a force in N.
I hope this covers your question.
Best regards
Daniel
Hi Daniel,
Thank you for the explanations, it looks like clearer now.
I just now would like, for a better understand, why you, in your tutorial, chose a pressure module corresponding to a thickness between 1 and 10cm. Does it correspond to any measure of the implant geometry ? If I’m right, th PE layer of a knee implant is thiner than this so why such a choice ?
Moreover, for a shoulder prosthesis, only the glenoid implant is in PE ; does it influence the calculation since elastic layer is not on the both parts, as in knee or in femur implants ? If we refer to [Blankevoort, 1991] the h (layer) in the equation comes from the addition of 2* h/2, from each part of the implant.
The choice of the layer thickness was covered by the idea to get the kinematics right. The example in the tutorial shows an artificial knee joint and thus, the value of the pressure module might also be a little bit artificial.
Which paper are you referring to by [Blankevoort, 1991]? I guess the layer thickness and thus the pressure module might change if you have contact between two elastic layers or one elastic layer and a stiffer surface.
Hi Daniel,
I calculated the elastic layer thickness corresponding to my pressure module and get a value 10* higher than the implant thickness. Do you think that this difference could comes from the material properties hypothesis we did (i.e. PE=linear material but not really the case) ?
Here is the complete citation for the paper : Blankevoort, L., Kuiper, J. H., Huiskes, R., & Grootenboer, H. J. (1991). Articular contact in a three-dimensional model of the knee. Journal of biomechanics, 24(11), 1019-31. Retrieved from http://www.ncbi.nlm.nih.gov/pubmed/1761580
I think this paper corresponds more or less to the origin of the study of [Bei and al., 2004].
As explained in [Li, G., 1997]*, it corresponds to an elastic layer of h/2 thickness on each side of the components in contact (the components being modelised as rigid bodies).
I also think that modeling layer + rigid body could conduct to some differences compared with modeling 2 elastic layers. That the reason of my questions : in the knee, your implant has elastic layers in both parts ; that’s not the case in shoulder prosthesis (humeral head considered rigid since in metal and a young modulus higher than the implant-PE Young modulus). Do you think that this aspect could be takes into account in my model ?
Thank you !
Lauranne
Li, G., Sakamoto, M., & Chao, E. Y. S. (1997). Technical note : A comparison of different methods in predicting static pressure distribution in articulating joints. The Journal of biomechanics, 30(6), 635-638
I think I have to admit that we had the same problem with some of our models were we reduced the pressure module to make the FDK solver work and solve the problem.
From the implementation of your model, I think the variations of modeling the contact are rather limit. Basically the only reasonable parameter you can change to modify the contact force is to use another pressure module and everything else is how you calculate it. At the moment, we are thinking about options to implement different force laws in the next release of AnyBody.
If you want to implement a more sophisticated contact model, my suggestion would be to define different layers of contact surfaces with a some offset and different pressure modules.
For more more detailed models to analyze the results, I would suggest to set up a FE contact model where it is probably easier to change the material properties. Therefor, you should have a look at the Any2FE converter if you want to transfer loads from the AnyBody model into the FE model.
I’m still working on my contact model for shoulder prosthesis. I would be interested in knowing how is computed the joint reaction forces. I understood it is a 3D vector located at the center of pressure. But could you explain to me for example if it is the sum of all reaction forces vectors, calculated at the points of the STL surface where there is penetration ? Or …??
Moreover, regarding the surface area, is there any way to know the shape of this area for both the master and the slaves surfaces ?
I’m not sure if I hit the core of your question, but I’ll just try.
So the joint reaction force is so to say the counterpart of all outer loads, so muscle forces of muscles spanning the joint and forces on the segments on the body. The joint reaction force is then the force needed to get an equilibrium. Thus, it is calculated by solving a system of equations. And as long as all forces have a representation as a 3D vector, the joint reaction force has it as well (compare Damsgaard et al., “Analysis of musculoskeletal systems in the AnyBody Modeling System”)
The AnyForceSurfaceContact applies two forces on the segments containing the surfaces (actio and re-actio so to say). This force is calculated as a using the penetration volume (compare ref manual) and the direction is calculated as weighted average of the vertex normals. This force is then applied to the segments, the same with the moment which then needs some transformation to have the same effect as if the force would be applied in the center of pressure (CoP).
About the surface area: the shape of the surface in contact is the same as the original surface, we are not doing any FE analysis and calculte deformations, so the surfaces are static. Unfortunatly there is no way in outputting the vertices in contact.
Hi Daniel,
thank you for your explanations. It makes my understood clearer.
I also had a look on the ref manual and saw the descriptions to improve/adapt the mesh quality. Did you do any trials to compare results with different mesh qualities (assuming that the study converges…!) in terms of reactions forces, muscle recruitment, etc… ? If yes are they big differences ?
Finally, I did my study converging. My question relates to the choice of the stiffness for the FDK. I began by choosing the same than yours (for the knee) and modifying to perform the analysis. I know this subject has already been discussed but… how to validate the choice of this stiffness ?! In the ref manual the FDK is introduced as a way to better represent the real joint and for example the ligaments. So my guess would be that we could found a relation between the spring of the FDKand the ligaments but… not sure !! Could you developp a little for me your process for choosing your spring stiffness ?
Thank you for your explanations and your time !
Lauranne
we did some comparison on the mesh quality and contact forces, so something like a convergence study, but it was during testing of the AnyForceSurfaceContact class and we haven’t published any results on that. In general I think this kind of study has to be done for each individual model, like it is handles for validation in FE studies.
I’m not sure what you meen by FDK stiffness. Are you talking about the model in the FDK tutorial?
Hi Daniel,
About the stiffness, yes, I’m talking about the tutorial (Lesson 2: Adding FDK Features). The stiffness for the AnyForce class is set at 1000 for x-direction, 5000 for y-direction et none for z-direction since it’s a 2D model. I would guess that these values refer to some previous data related to the measure of the forces of ligaments or such a thing. Consequently I guess you validated your choice as such a manner. Am I right ? Or did you set these values following several trials to make your model converging ?
It is very likely that these numbers are arbitrary since the purpose of this lesson is just the demonstration of the methodology rather than answering any research questions. Thus, in turn, this oversimplified example was not validated.