Need help in modeling indentation stress relaxation of cartilage

[vgupta0603]

Hi

I am a graduate student at University of Kansas. I have been trying to model the mechanical testing I did for one of my projects. I am trying to simulate a indentation test, first part of which includes compressing the tissue to 15% strain at a rate of 0.15%/sec and the second part of the test allows the tissue to relax for a period of 1000 seconds.

When I try to model the experiment, the simulation fails to converge at 36 seconds as negative jacobian is detected. I have tried reducing the time step but it doesn't seem to solve the problem.

I would greatly appreciate if you could have a look at my files and give me some tips/pointers to fix the problem.


Regards
Vineet Gupta

[ateshian]

Hi Vineet,

Here are a few suggestions that might help:
- Make sure that the tension-compression contact is set to sustain contact.
- Decrease the penalty factor for the biphasic contact from 100 to 10 (or even 1, in which case you may need to turn on Lagrangian augmentation). A penalty value that is too high may produce convergence problems.
- Change the geometry of your indenter to match the wedge geometry of the biphasic layer. This will enforce axisymmetric conditions better.
- If these changes are not sufficient, consider using mesh refinements to produce a finer mesh near the edge of the indenter.

Let me know if these suggestions work.

Best,

Gerard

[vgupta0603]

Dr. Ateshian,

This is what I tried based on your suggestions: -
- Turned on the tension flag for the TC contact.
- Decreased the penalty factor from 100 to 1 (with Lagrangian augmentation turned on)
- Changed the geometry of the indenter.

Still the model fails to converge. The message displayed at termination is "Maximum no. of reformations reached". Do you think increasing the no. of reformations will do any good?

Lastly, I am not sure how to create a finer mesh at the indenter edges. Is there a way in Preview to create a biased mesh? Or do I have to use another mesh generator software?

I am unable to attach my .feb file here but if you would like to have a look, I can email it to you.

Thanks and Regards
Vineet

[ateshian]

Hi Vineet,

For this indentation problem you can create a biased mesh as follows:

1) Create a wedge with height equal to the tissue layer and a radius equal to the indenter radius. Bias the mesh to make it finer at the radial edge. Also use double bias to refine the mesh near the surface and bottom of the layer. Convert the mesh to editable form.
2) Create a hollow tube whose inner radius is equal to the indenter radius, outer radius equal to the desired value, and height equal to the layer thickness. Bias the mesh to make it finer near the inner edge. Use the same bias through the thickness as for the cylinder. Delete elements to produce a wedge.
3) Attach the two geometries and weld the nodes at the interface where the geometries meet.

This allows you to create a biased mesh with finer elements near the edge of the indenter.

Best,

Gerard

[vgupta0603]

Dear Dr. Ateshian

I tied to incorporate your suggestions into the model but still it doesn't work. I have attached the modified file here. Please let me know where I made mistakes.


Regards
Vineet

[ateshian]

Hi Vineet,

I have attached a working model that reproduces your analysis. Here are a couple of additional remarks that should help explain why this problem is challenging:

1) The sharp edge of the rigid flat indenter will cause very large distortions to the elements that come into contact with it. As a result, I found it necessary to reduce the penalty factor further, down to 0.2. This reduction effectively "softens" the rigid indenter, allowing more overlap with the biphasic layer. To reduce this overlap I decreased the tolerance on the augmented Lagrangian from the default value of 0.20 down to 0.05. Of course, this increases the total number of iterations.

2) I added the lateral (cylindrical) face of the indenter to the master surface of the contact interface, since the biphasic material bulges out significantly during the deformation, coming into contact with that lateral side.

Check out CartIndSR2.feb, once you run the analysis you'll be able to see the large deformation and significant distortion that occurs at the edge of the indenter.

Best,

Gerard

[vgupta0603]

Thanks a lot Dr. Ateshian! The model works and I can see the large deformation occurring near the edges of the indenter.

But I have a few questions regarding the model that I want to clarify: -
1) It looks like the meshing of the indenter different from that of the tissue. Is it?
2) You have used the neo-hookean solid for solid part of the tissue and I used isotropic elastic in my earlier models. How does that impact the overall analysis?
2) You defined the rigid constraints in the second step of the analysis. Is there a reason for that?

Also, I tried to curve fit my experimental data. The optimization runs up to 15 major iterations and 64 minor iterations and takes too much time. Is there a way to minimize the time for optimization? I tried to decrease the <f_diff_scale> to 0.01 but then the optimization terminates prematurely. Also if I use any other initial guesses than what I used in the attached file the optimization terminates prematurely. Lastly, when I plot the data from the .log file the differences between the model load and experimental loads seem quite large. Do you think the model needs further refinement?

I have attached a figure as well as my optimization input files. Please let me know what you think.

Thanks once again for your help. I greatly appreciate it. This is my first time working with FEBio and I don't have much experience in the Bio mechanics modeling as well.

Regards
Vineet

[vgupta0603]

Dear Dr. Ateshian,

Please find attached the image file from the optimization run.



Regards
Vineet

[ateshian]

Hi Vineet,

1) It looks like the meshing of the indenter different from that of the tissue. Is it?

Since the indenter is flat, one element suffices to describe its geometry exactly. Since the indenter is rigid, the contact analysis must be a single-pass which treats the indenter surface as the master surface. This means that the contact analysis only samples points on the biphasic surface (the slave surface) to detect contact with the master surface. That's why it only matters to have a fine mesh on the biphasic surface.

2) You have used the neo-hookean solid for solid part of the tissue and I used isotropic elastic in my earlier models. How does that impact the overall analysis?

As you can see, indentation with a flat punch causes very high strains and rotations of the elements near the indenter edge. The isotropic elastic material is not formulated for large deformations and does not behave well in the limit of large compressions (it does not produce an infinite stress value when compressing the material to zero volume). So it is better to use a neo-Hookean material to minimize the change of getting element inversions during the analysis.

2) You defined the rigid constraints in the second step of the analysis. Is there a reason for that?

No particular reason, this was just the default setting on PreView when I created these constraints. Since there is only one step in this analysis, the results would be the same either way.

Also, I tried to curve fit my experimental data. The optimization runs up to 15 major iterations and 64 minor iterations and takes too much time. Is there a way to minimize the time for optimization? I tried to decrease the <f_diff_scale> to 0.01 but then the optimization terminates prematurely. Also if I use any other initial guesses than what I used in the attached file the optimization terminates prematurely.

I don't know what you mean exactly by too much time. 20 to 30 minutes would be considered normal for this type of problem. If it is taking several hours then I agree there is a problem. What you need to understand is that the conversion characteristics do depend on the values of the material properties, so there may be certain combinations of properties that require many more iterations than the analysis I uploaded for you. Part of this problem is that indentation with a flat frictionless indenter is numerically challenging as discussed above. One way to reduce the complexity of the analysis is to assume that the contact interface between indenter and biphasic layer is adhesive, which fixes the elements under the indenter and overcomes the challenge of elements sliding back and forth across the indenter edge while distorting and rotating. You can create a rigid interface instead of a biphasic contact interface. Keep in mind that the nature of the problem also changes in this case.

Lastly, when I plot the data from the .log file the differences between the model load and experimental loads seem quite large. Do you think the model needs further refinement?

I am not sure why this is the case, but it is possibly related to the implementation of a new optimization algorithm in FEBio 1.8 which apparently is not behaving as we would expect. Try running the optimization analysis with FEBio 1.7 or 1.6 and see if you get the same results. Let me know what you find.

Best,

Gerard

[vgupta0603]

Dear Dr. Ateshian

I was using FEBio 1.5 earlier but now I switched to version 1.7. The optimization results look better but still the experimental data and the FEBio data don't seem to match very well (files attached). What would be your suggestion further to overcome the problem?

Thanks and Regards
Vineet

[ateshian]

Hi Vineet,

The FEBio analysis and curve-fitting algorithm appear to be working fine. The issue you are confronting is more likely related to the choice of constitutive relation for the tissue. I assume you are modeling articular cartilage (based on your filename). In this case please note that the large peak to equilibrium stress response you observed in your experiment is due to the large disparity between the tensile and compressive moduli of the tissue. This disparity results from the fibrillar nature of cartilage whose collagen matrix offers much greater resistance in tension than compression. To improve your experimental fit, you could choose to model the solid matrix as a solid mixture of a neo-Hookean ground matrix and a spherical fiber distribution. This works well for the unconfined compression stress-relaxation experiments we conduct in my lab.

Best,

Gerard

[vgupta0603]

Dear Dr Ateshian,

I changed my material to solid mixture in the .feb file but I am unable to run it in the optimization mode. However, I can run it in the input mode. Can you please have a look and let me know what I did wrong?


Thanks and Regards
Vineet

[ateshian]

Hi Vineet,

I posted a new example for parameter optimization of a biphasic material whose solid matrix is a mixture, see here. This should help you figure out the naming convention for this type of material.

Best,

Gerard

[vgupta0603]

Dear Dr. Ateshian,

Thanks for posting the example. I think I got the naming convention right this time but now the optimization starts but terminates prematurely with the following error message: ************************************************** ***********************
* F A T A L E R R O R *
* *
* FEBio error terminated. Parameter optimization cannot continue. *
************************************************** ***********************

I am sending you my optimization files. I would appreciate any suggestions that you may have.

Thanks and Regards
Vineet

[ateshian]

Hi Vineet,

Parameter optimizations may terminate with an error if the properties used in a trial produce an error during an analysis, such as a negative jacobian. This means that the model is not sufficiently robust to handle the full range of properties specified in an optimization. To understand which error has occurred during the optimization you can choose to print the log of each optimization iteration (either to the screen or the log file). Check out the log_level keyword in the Help manual.

Best,

Gerard