large compression contact problem

[remiverhoeven]

Hi!

I am currently working on a contact problem and am hoping for some feedback to make it more stable.

The model concists of a 1/8th sphere contacting a rectangular box. The materials used in the sphere are all Ogden constrained, whilest the box is neo-Hookean. I am using the Paradiso (default) solver and a 'sliding2' contact type.

The model represents a part of the buttocks (spheres) that indents a cushion (the box) with a specific force (weight). I've tried a default analysis with fixed displacement aswell as a dynamic analysis, but this seemed to make little difference in stability.

The problems I have been having are:

1. the (mesh)size of the box.
I can not increase the mesh size (amount of elements/nodes) of the box without making the situation unstable. Ideally I would like a 10x10x4 element box, but this will not run in my present configuration. I have tried varying the models parameters (penalty, gaptol, dtol, tolerance, time_step)m but they all seem to end up with either high-penetration results or unstable solutions (negative jacobian).

2. Jumping nodes.
In the first and last steps of the solution the sphere 'moves' up nodes to connect them to the box. Basically all parameter variations I have tried that found a solution have ended up with these contact problems.

3. Time consuming.
Although this is less of a problem, the present model - with just 1395 nodes - takes about an hour to run. Most of this time is lost on the first and last iterations on remodeling the 'maximum gap' to obey the gap tolerance. These situations are also the situations that the jumping nodes usually occur.

Are there any ways you might be able to assist me with these problems and/or give me some pointers based on the model? I have attached the model for your convenience.

[ateshian]

Hi,

I looked at your file and here are a few recommendations:

- For your contact interface you selected auto-penalty="Yes" (a good choice), but then you selected a penalty factor of 0.001. It is best to select a penalty factor between 1 and 100. A higher value will produce less penetration (more faithful contact), but might lead to negative jacobians and slow down convergence. A lower value improves convergence but may produce undesirable penetration. You can start with a value of 10 and see how that works. Also, let augmented Lagrangian="No", unless you are not satisfied with the amount of penetration. This will make the analysis complete faster.

- In your contact problem, the bottom rectangular body is meshed with a single element. This means that the contact surface only consists of one face (coarse mesh). In these types of problems you should use two-pass="No" (i.e., a single-pass analysis) and select the coarsely-meshed surface as the master surface. This produces a more accurate outcome, while giving you the benefit of using a coarse mesh where feasible.

- Your contact analysis is purely under load control: You push one body against the next by prescribing a load. Unless the contacting bodies are initially overlapping a little bit, a pure load control analysis will not work well for contact problems. You should consider performing a two-step analysis: In Step 01, prescribe a small displacement to push one body against the other. In Step 02, apply the load, starting from a small (but non-zero) value.

- Your bottom rectangular block appears to be very stiff compared to the "buttock". You could consider modeling it as a rigid body. This would have the benefit of facilitating the application of prescribed displacement and forces. (Remember to use single-pass analysis when performing contact with a rigid body, and pick the rigid contact surface as the master.)

If you implement these changes, I expect that you will see much improved behavior.

Best,

Gerard

[remiverhoeven]

Hi Gerard,

Thank you for your swift reply. I have looked over the points you adressed.

One of my initial issues was the detail in the mesh of the box; I couldn't get it to work with a finer mesh then the 1 element version. Due to your hints I have decreased the stiffness of the box to a more acceptable and realistic stifness and now it seems to work just fine. I also had to use the Langrangian="Yes" after changing all parameters you suggested because of the terrible penetration I got with acceptable penalty factors that don't end up in negative Jacobians.

I did implement the two loadcurve setup you described. This seems to work out alright, although the load application clips the FEBio application when I enter slightly higher loads (in excess of 1000 per node or a total of 100000). Would you recommend ramping up the load to get to the target load? How would I go about this task? Would specifying a more exponential loadpoint distribution (e.g. 1,0 ; 1.1,0.004 ; 1.2,0.012 ; 1.3,0.025 ; 1.4,0.076 ; 1.5,0.130 etc.) work to stabalise the calculation?

Best
Remi

[ateshian]

Hi Remi,

I have made some modifications to your sample file to make it run better. Please see the attached version.

Here is a brief summary of what I did:

1) I created a rigid body that I attached to the bottom of the rectangular slab using a rigid interface.

2) I created a two-step analysis, in the first step (from t=0 to t=1) I prescribed a displacement to the rigid body along z (from 0 to -1). In the second step (from t=1 to t=2), I prescribed a load to the rigid body along z (from -2e4 to -2e6). I saved the Fz component of the rigid body force in the log file and I first ran the analysis from t=0 to t=1 to figure out the value of Fz at t=1.

3) For the contact interface I used a single-pass analysis, with a auto-penalty on and a penalty factor of 100.

On my computer this analysis runs in 38 s and it completes satisfactorily. Try it out and see if that works for you.

Best,

Gerard

[remiverhoeven]

Thank you so much Gerard!

[remiverhoeven]

Dear Gerard,


I have now managed to make the model run with a not-rigid cushion with a fixed displacement in step 1. I have been trying to get the model running with a force in step 2, just like you managed to do in the rigid body version of the model. However, I have not succeded in getting succesful iterations in the second step.

In step 1 I'm using:

<time_steps>20</time_steps>
<step_size>0.05</step_size>
...
<dtol>0.01</dtol>
<etol>0.01</etol>
<rtol>0</rtol>
...

with a compression of about 25% of the model.

In step 2 I'm using:

<time_steps>80</time_steps>
<step_size>0.05</step_size>
...
<dtol>0</dtol>
<etol>0.01</etol>
<rtol>0</rtol>
...
<analysis type="dynamic"/>
...
<loadcurve id="2" type="smooth">[INDENT]<loadpoint>0,0</loadpoint>
<loadpoint>1,-0.1</loadpoint>
<loadpoint>4,-1</loadpoint>[/INDENT]
...

I have tried increasing the mesh-density of the cushion and the model, but neither seem to work.

I am also bound to the more general parameters for the sliding interface:

<contact type="sliding2">
<laugon>1</laugon>
<penalty>0.6</penalty>
<auto_penalty>1</auto_penalty>
<two_pass>0</two_pass>
<gaptol>0.001</gaptol>

I have tried a higher penalty but it wont run on my system. Using Langrangian="No" will result in pretty high penetration (makes sence with a small penalty factor...)

Would you know what I could improve to get this setup to work in both the displacement and the force step?
120309_RV_public_v3.feb

[ateshian]

Hi Remi,

When you transition from displacement control to load control, you need to know how much load was being applied at the end of the displacement control step. This value should then serve as your starting point for the force load curve in the load control analysis.

Since you are no longer using a rigid interface as the substrate of the rectangular part, you need to evaluate the load in PostView. There are different ways of doing this, for example you can integrate the normal stress on the bottom face of the rectangular part (select the faces, select the Z-Stress, then Post->Integrate).

Best,

Gerard

[remiverhoeven]

Dear Gerard,

Thanks for the information.

I am currently working on an individualised model to calculate the strains in the layers of the buttocks and comparing them between patients. Ideally I would like to automate this process for every patient. Currently I am working on running the model twice in a row: first to obtain the right force after the displacement step, second to actually apply this force to the model.

Is it possible to 'pause' or 'interupt' the model, change the initial load in the file and then continue with the second step? If it is not I dont have a problem running the model twice, but if there is a nice solution to this problem, I'd like to use it

Best,
Remi

[ateshian]

Hi Remi,

My understanding is that there is a restart option in FEBio that should allow you to do precisely that. However I haven't used it, so I can only refer you to the user's manual: Chapter 5 Restart Input file

Also check out 2.8. Advanced Options for a description of how to run a restart file.

More entries can be found on the online help manual if you search for "restart".

Let me know if that works for you.

Best,

Gerard

[remiverhoeven]

Dear Gerard,

Thanks for the information.
I have looked through the restart input file information, but I am afraid that will complicate the input required from the research assistants. I have now implemented a routine that runs febio with the displacement only. It calculates the force based on its output and then creates a new feb-file with displacement and force. Since the displacement step doesn't take that much time that wouldn't be to much of a time-sink.

Thank you very much though, for all the help you have provided


Remi

[ateshian]

Hi Remi,

I am glad you were able to figure out a workable approach.

Best,

Gerard