Proper contact on folding surfaces

[Fluvio]

Hello,

Is there an appropriate way of setting-up a contact BC for a surface that folds onto itself?
state0.PNGstate75.PNG

Would it make sense to use sliding contact with the surface as a master and a slave?


Thank you!

[ateshian]

Hi Fluvio,

Self-contact is allowed in the sliding contact interfaces. Just select the same surface(s) as primary and secondary. When using sliding-elastic contact, make sure to reduce the search_radius to a small enough value.

Best,

Gerard

[Fluvio]

ateshian,

Apologies it took me this long to respond. I have been trying to simplify these meshes so our simulation does not fail before getting to the contact portion of the problem. I am also now using the latest version of FEBioStudio.

My model is a little different now, but the problem is still the same (Figure 1-2). I followed your suggestion and went from the default search_radius = 1 to 0.1 and 0.01. In either case, the simulation ended at 76-77%, with 0.1 getting further.
Here is a summary of the sliding-elastic contact parameters;

<Contact>[INDENT]<contact type="sliding-elastic" name="SlidingElasticContact1" surface_pair="SlidingElasticContact1">
<laugon>0</laugon>
<tolerance>0.2</tolerance>
<gaptol>0</gaptol>
[B]<penalty>10</penalty>[/B]
<auto_penalty>1</auto_penalty>
<two_pass>0</two_pass>
<search_tol>0.01</search_tol>
<symmetric_stiffness>0</symmetric_stiffness>
[B]<search_radius>0.1</search_radius>[/B]
<seg_up>0</seg_up>
<tension>0</tension>
<minaug>0</minaug>
<maxaug>10</maxaug>
<fric_coeff>0</fric_coeff>
<smooth_aug>0</smooth_aug>
<node_reloc>0</node_reloc>
<flip_primary>0</flip_primary>
<flip_secondary>0</flip_secondary>
<knmult>0</knmult>
<update_penalty>0</update_penalty>[/INDENT]
 </contact>

Are there any other parameters that I should be considering?
How does the penalty factor play into this type of contact?

Here are some pictures for references;

state0.PNG
state077.PNG

As always, thank you for your help!

[ateshian]

Hi Fluvio,

The fact that you got so far means that the analysis was set up properly, so that's a good place to be.

Difficulties in contact convergence can occur to a variey of factors, but here are a few more things you can try:

• Make sure that the analysis step uses a non-symmetric stiffness matrix (because this sliding-elastic contact interface has produces a non-symmetric stiffness).
• You can turn laugon=1 and reduce the penalty from 10 to 1. Start with the default tolerance (0.2), but you can reduce it if you find that the contacting surfaces overlap too much.

Best,

Gerard

[Fluvio]

Gerard,

Great news!
I actually got the simulation to converge using the penalty method (penalty = 1, search_radius = 0.1).
First, I confirmed that all of the simulations were non-symmetric and the quasi-newton method used was broyden.

Here is a summary of the simulations I ran;

auto_penalty = 1, penalty = 10, search_radius = 0.1 -----------> 77%
auto_penalty = 1, penalty = 10, search_radius = 0.01 -----------> 76%
auto_penalty = 1, penalty = 100, search_radius = 0.1 -----------> 71%
auto_penalty = 1, penalty = 1, search_radius = 0.1 -----------> 100%

I am never clear of the interaction between the auto_penalty flag and the penalty value itself, when both are active of course. I wonder if this had an effect?

[ateshian]

When auto-penalty=0, the penalty value represents the contact stiffness, you need to provide it in units of force/length.

When auto-penalty=1, the penalty value is a scale factor that multiplies an estimate of the contact stiffness obtained from the material properties (Young's modulus) and local finite element thickness (element volume divided by area of the contact face).

Best,

Gerard

[Fluvio]

Gerard,

After getting this convergence on the folding surface, I went back and added a second contact condition that describes the tissue sliding along the curved chest wall.
The second contact condition is also defined by a sliding-elastic formulation, with tension enabled.
I am having issues with converging using the penalty or the augmented Lagrangian method, so just wanted to get your thoughts...

First, I used the penalty method;

<contact type="sliding-elastic" name="SlidingElasticContact2" surface_pair="SlidingElasticContact2">[INDENT]<laugon>0</laugon>
<tolerance>0.2</tolerance>
<gaptol>0</gaptol>
[B]<penalty>120</penalty>[/B]
[B]<auto_penalty>1</auto_penalty>[/B]
<two_pass>0</two_pass>
<search_tol>0.01</search_tol>
<symmetric_stiffness>0</symmetric_stiffness>
<search_radius>1</search_radius>
<seg_up>0</seg_up>
[B]<tension>1</tension>[/B]
<minaug>0</minaug>
<maxaug>10</maxaug>
<fric_coeff>0</fric_coeff>
<smooth_aug>0</smooth_aug>
<node_reloc>0</node_reloc>
<flip_primary>0</flip_primary>
<flip_secondary>0</flip_secondary>
<knmult>0</knmult>
<update_penalty>0</update_penalty>[/INDENT]
 </contact>

As opposed to the previous contact problem, larger penalty values helped with convergence, but only until 120. With a penalty = 120, the simulation fails at 0.453 (45%).
When penalty was set to >120 (e.g. 130, 150, and 200), I did not get the first step to converge. Note that I am using the auto_penalty too.

Here I figured I could try using the augmented Lagrangian method for once.

<contact type="sliding-elastic" name="SlidingElasticContact2" surface_pair="SlidingElasticContact2">[INDENT]<laugon>1</laugon>
[B]<tolerance>0.01</tolerance>[/B]
[B]<gaptol>1</gaptol>[/B]
<penalty>1</penalty>
<auto_penalty>0</auto_penalty>
<two_pass>0</two_pass>
<search_tol>0.01</search_tol>
<symmetric_stiffness>0</symmetric_stiffness>
<search_radius>1</search_radius>
<seg_up>0</seg_up>
[B]<tension>1</tension>
<minaug>5</minaug>
<maxaug>10</maxaug>[/B]
<fric_coeff>0</fric_coeff>
<smooth_aug>0</smooth_aug>
<node_reloc>0</node_reloc>
<flip_primary>0</flip_primary>
<flip_secondary>0</flip_secondary>
<knmult>0</knmult>
<update_penalty>0</update_penalty>[/INDENT]
 </contact>

With a gap_tol = 0.1 and a tol = 0.01 I get as far as 0.30 (30%). I got the furthest with a gap_tol = 1 (0.36 or 36%). But increasing the gap_tol > 1 (e.g. 10) or increasing the tol > 0.01 (e.g. 0.2), I only got to 0.25 (25%).
I know there is a risk with having both gap_tol and tol on at the same time, according to the manual, but perhaps I am missing something else?

[ateshian]

Hi Fluvio,

In general, when using a stringent value for <gaptol>, you should turn off the <tolerance> (set it to 0). Similarly, if you turn on <tolerance>, I recommend turning off <gaptol>. When either of these tolerances is set to be too tight, you may find that your analysis can never achieve those desired tolerances, especially if the mesh is relatively coarse. In that case however, the code will simply augment until maxaug is reached, then continue with the next time step. So, in principle, the values of the tolerances should not affect the termination time of the analysis. This makes me suspect that you might be using a quasi-Newton method to solve this problem (e.g., BROYDEN or BFGS)? If so, I recommend switching to a full-Newton analysis. The reason for that is that the stiffness matrix in a contact analysis changes as the contact interface evolves, so a quasi-Newton method may lead to poorer convergence and, possibly, premature termination.

For the problem description you have provided, it is unclear to me why the tension flag should be turned on. Have you tried running the analysis with that flag turned off?

Best,

Gerard

[Fluvio]

Gerard,

Let me clarify. I have two contact conditions going on.

The first contact condition (left) is for the skin folding over itself (this was the original problem, which I was able to solve with a sliding-elastic contact(auto_penalty = 1, penalty = 1, search_radius = 0.1) and using BROYDEN
The second contact condition (right), on the same problem, is for the tissue sliding along the chest wall (red arrows). I am also using a sliding sliding-elastic contact and having the issues described above. This also uses BROYDEN
We have applied sliding to previous models. If you recall, my colleague Tres and I worked on a normal displacement constraint this past summer.

folding_s1.PNGslinding_and_folding_s1.PNG

I will try using the full-Newton method. In order to set-up the full method, do I set the number of max_updates = 0?
I will try a few things, but this response should sort of clear the approach.

[ateshian]

Hi Fluvio,

Thanks for the clarification. Yes of course, I remember that you and Tres worked on this issue.

Setting max_updates to 0 indeed implements the full-Newton method.

Best,

Gerard

[Fluvio]

Gerard,

We are getting closer, but we still have some issues...
I went ahead and started using a full-Newton method by setting <max_ups = 0> and also increased <max_refs = 50>, which I found advised on the FEBio user's manual.
Sadly, I still do not get a convergence with both contact conditions together. In the best scenario, my simulation fails between 46% and 48% of the way through.

I went ahead and split the problem in two parts. First, applying only the sliding elastic contact (Contact #1) to avoid surface penetration when the tissues folds over.
I was able to get through the entire simulation with the following parameters;

<Control>[INDENT]<analysis>STATIC</analysis>
<time_steps>10</time_steps>
<step_size>0.1</step_size>
<solver>
[B]<max_refs>50</max_refs>[/B]
[B]<max_ups>0</max_ups>[/B]
<diverge_reform>1</diverge_reform>
<reform_each_time_step>1</reform_each_time_step>
<dtol>0.001</dtol>
<etol>0.01</etol>
<rtol>0</rtol>
<lstol>0.9</lstol>
<min_residual>1e-20</min_residual>
<qnmethod>BROYDEN</qnmethod>
<rhoi>0</rhoi>
<symmetric_stiffness>0</symmetric_stiffness>
</solver>
<time_stepper>
<dtmin>0.01</dtmin>
<dtmax>0.1</dtmax>
<max_retries>5</max_retries>
<opt_iter>10</opt_iter>
<aggressiveness>1</aggressiveness>
</time_stepper>[/INDENT]
 </Control>
...
<Contact>[INDENT]<contact type="sliding-elastic" name="SlidingElasticContact1" surface_pair="SlidingElasticContact1">[/INDENT][INDENT=2]<laugon>0</laugon>
<tolerance>0.2</tolerance>
<gaptol>0</gaptol>
[B]<penalty>1</penalty>[/B]
[B]<auto_penalty>1</auto_penalty>[/B]
<two_pass>1</two_pass>
<search_tol>0.01</search_tol>
<symmetric_stiffness>0</symmetric_stiffness>
[B]<search_radius>0.1</search_radius>[/B]
<seg_up>0</seg_up>
<tension>0</tension>
<minaug>0</minaug>
<maxaug>10</maxaug>
<fric_coeff>0</fric_coeff>
<smooth_aug>0</smooth_aug>
<node_reloc>0</node_reloc>
<flip_primary>0</flip_primary>
<flip_secondary>0</flip_secondary>
<knmult>0</knmult>
<update_penalty>0</update_penalty>[/INDENT][INDENT]</contact>[/INDENT]
 </Contact>

Then, I created a model with only with the sliding elastic contact allowing the sliding of the tissue along the chest wall (Contact #2).
I believe the entire issue is related to this contact condition, but do not know if changing the simulation parameters is enough...
I have tried;
- penalty = 120, auto_penalty = 1 -----------------> 48%
- augmentation_tol = 0.2, gap_tol = 0 -------------> 46%
- augmentation_tol = 0, gap_tol = 1 -------------> 46%

I have tried other combinations, with higher and lower tolerances and penalty factors, but these are the best results I get... always around 46-48%
Is this an indication that the mesh needs refinement?
...maybe re-meshing?

Thank you!

[ateshian]

Hi Fluvio,
Would you be willing to share your model file with me? Perhaps I can provide more insights if I could try it out myself.
Best,
Gerard

[Fluvio]

Gerard,

Of course. I was kind of hoping you would ask. Can I email it to you?

Thank you!

[ateshian]

Yes, ateshian@columbia.edu.

[ateshian]

Hi Fluvio,

Thanks for emailing me the file. I tried a few things to get the model to converge, here I summarize what I did to get it to run to completion:
- Because the model was taking too long to run, I remeshed the "skin" domain using "Tet Remesh" with Element size: 6, Min element size: 3. This brought down the number of equations from > 500 K to < 100 K.
- For the sliding-elastic contact with tension=No, I set laugon=1 and kept tolerance=0.2, penalty=1, auto-penalty=Yes, two-pass=Yes.
- For the sliding-elastic contact with tension=Yes, I set laugon=1, tolerance=0, gap tolerance=0.15, penalty=10, auto-penalty=Yes, two-pass=No.
- I added a boundary condition: Fix X on the surfaces for which you had also fixed Y. This prevented part of the torso from expanding laterally along X as gravity increased, while another part was fixed from the fixed BC you had applied along X.
- I increased the optimal iterations (opt_iter) to 20.

With these changes the model ran to completion in about 11 minutes. I have emailed you a link to my updated model files so that you can view and run this model yourself. Let me know if you are unable to reproduce those results, or if you require clarifications about these choices.

Best,

Gerard