Mooney Rivlin in Shell. Results slightly off and Augmented Lagrangian not working

[MiguelA]

I created a Disk with a Mooney-Rivlin material to test FEBio against the analytical solution of a bulge test as given by this paper [1].

The parameters are:

• Radius = 4mm
• Thickness = 0.05mm
• mu = 4MPa -> c1=2MPa
• Pressure = 28 kPa

Max displacement from analytical solution: 1.5302715

In this figure is the displacement that I obtain as a function of the bulk modulus. I varied the bulk modulus between the most extreme values that still converged:
EffectOfBulkModulus.png

As you can see I can get a result that is close to the analytical solution for k/c1 in [10,10000]. However, the result seems to be consistently 6% above the analytical solution (i.e. more compliant than expected). This didn't change with mesh refinement.

I then tried to see if I could get the augmented Lagrangian enforcement of incompressibility by adding "<laugon>1</laugon>" and "<atol>0.01</atol>" but unfortunately that seemed to have no effect. Am I doing something wrong? Attached is my .feb file.


AugLagTest.feb

PS: If I switch the BCs to FixedDisplacement instead of FixedShellDisplacement the error reduces to 3% instead of 6%.


[1] Sheng, J.-Y., Zhang, L.-Y., Li, B., Wang, G.-F., & Feng, X.-Q. (2017). Bulge test method for measuring the hyperelastic parameters of soft membranes. Acta Mechanica. https://doi.org/10.1007/s00707-017-1945-x

[ateshian]

Hi Miguel,

Thanks for providing the literature reference and a clear presentation of your problem. I was able to replicate the issues you found using the file you attached. Here is how to fix the analysis:

1) For this particular problem, as you noted, you should fix the shell front nodal displacements on the boundary (x,y,z), instead of the back nodal displacements (sx,sy,sz), since the pressure is being applied on the front face of the shell by default.

2) The shell formulation you are using is new to FEBio 2.7. The old shell formulation did not implement the standard three-field formulation needed for uncoupled materials, so the new formulation does not turn on that feature by default. You need to turn it on explicitly using this (unfortunately not yet documented) command:<use_three_field_shell>1</use_three_field_shell> which should be inserted in the <Control> section. This will allow you to do Lagrange augmentations with shells. You don't need to set the tolerance to 1e-08, I think 1e-3 should suffice.

3) For large deformations in problems that use the pressure boundary condition, I recommend setting <symmetric_stiffness>0</symmetric_stiffness> in the <Control> section. I also recommend using full Newton iterations instead of Broyden/BFGS.

I ran the analysis that you attached under three different conditions (I modified the time steps and increments to get better convergence, but that's easy to do by trial and error):

a) k/c1 = 1000, laugon=0.
b) k/c1 = 100, laugon=0.
c) k/c1 = 100, laugon=1.

I analyzed the response for pressure ranging from 0 to 90 kPa and I plotted the FEBio displacement at the center of the circular membrane on top of Figure 3 of the paper you cited.
AugLagTest.pdf
Results show that case (a) agrees well with the ABAQUS simulation in the paper, case (b) deviates a little bit, and case (c) gives results as good as (a) or perhaps slightly better, when compared to the analytical solution.

Best,

Gerard

[MiguelA]

Thank you professor Ateshian,
following your suggestions I was able to get the augmented lagrangian approach to incompressibility working.

The results are very satisfying reaching the best accuracy of all alternatives like your indicated and independent of the bulk modulus (figure attached).
ComparisonOfDefault3FieldAndAugLag.png

Here is the feb file with the changes you mentioned commented if someone would like to try it out:
AugLagTest_fixed.feb

I have two questions:
1) In your reply you mentioned that this new shell formulation by default does not turn on the 3-Field formulation that is needed for uncoupled materials. What happens when we use the default formulation with an uncoupled material? Because in this example it seems to work fine with the Mooney-Rivlin. Did it switch to the coupled formulation or is the stiffening observed for K/C1>10000 the effect of locking?

Also, may I suggest that if a user uses an uncoupled material without turning on the "use_three_field_shell" that FEBio stops with an error message like this: "ERROR: The material model used is uncoupled but the default shell formulation is displacement based, which can lead to locking issues. To avoid this error switch to a compressible material or add the following line to the control section of the feb file <use_three_field_shell>1</use_three_field_shell>". Just a suggestions that would warn users about this issue. I guess that a value of 2 could be used for not using 3-Field but also skipping the error.

2) This is more a curiosity regarding the shell boundary conditions. If I understood correctly from the theory manual, this formulation accounts for the rotation by the relative positions of nodes below and above the mid surface of the shell (front and back). Does this mean that this formulation could model a boundary condition of fixed displacements and fixed rotations by fixing both front and back displacements?

Also, would it make sense to have the fixed displacements BC be that the average displacement between the front and back nodes be zero?

Thanks again for the help (and for FEBio in general)

All the best,
Miguel

PS: I just realized that setting the "<use_three_field_shell>1</use_three_field_shell>" when the Control section is within the Step section gives the error:
FATAL ERROR: unrecognized tag "use_three_field_shell"

[ateshian]

Hi Miguel,

I am glad you were able to get such good results with the 3-field formulation for shells turned on.

1) In your reply you mentioned that this new shell formulation by default does not turn on the 3-Field formulation that is needed for uncoupled materials. What happens when we use the default formulation with an uncoupled material? Because in this example it seems to work fine with the Mooney-Rivlin. Did it switch to the coupled formulation or is the stiffening observed for K/C1>10000 the effect of locking?

When the 3-field formulation is turned off the analysis still uses the uncoupled formulation, but the relative volume J is not assumed to be uniform over the entire element (which is the assumption of the 3-field formulation). The calculation of the residual is the same with/without 3-field, but the stiffness matrix calculation without 3-field is based on the non-uniform J formulation, whereas with 3-field it accounts for the uniform J. So differences in the results with/without 3-field relate to the difference in calculation of J and yes, locking will occur with increasing values of the bulk modulus K when 3-field is off.

Also, may I suggest that if a user uses an uncoupled material without turning on the "use_three_field_shell" that FEBio stops with an error message like this: "ERROR: The material model used is uncoupled but the default shell formulation is displacement based, which can lead to locking issues. To avoid this error switch to a compressible material or add the following line to the control section of the feb file <use_three_field_shell>1</use_three_field_shell>". Just a suggestions that would warn users about this issue. I guess that a value of 2 could be used for not using 3-Field but also skipping the error.

I agree with your concern. In the next release I'll look into modifying the code so that the 3-field formulation is turned on automatically when using certain shell elements, but we need to do more testing to understand this issue. For example, does the 3-field formulation work equally well with tri3 shell elements? Do we need it with quadratic shell elements (tri6, quad8)? Currently, for solid elements, only hex8, penta6, pyra5 and quadratic tet (e.g., tet10, tet15) elements use the 3-field formulation by default. In other words, tet4, hex20, hex27, penta15 elements don't. We would need to ask Steve for his input on those decisions, but my guess would be as follows: One should never use tet4 elements with uncoupled materials since they lock; this is taken to be a dogma in finite element analysis, even though all finite element codes (commercial and free) provide users the option to use tet4 elements with such materials. I think that 3-field is turned off for those elements by default because the analysis is more efficient that way (even though locking can occur). For quadratic elements such as hex20, hex27 and penta15, I am going to guess that they probably work well and don't lock as easily as linear hex and penta elements (but I haven't personally tested this matter extensively), which is why the 3-field option is turned off by default. (BTW, one can use <use_three_field_hex>1</use_three_field_hex> to turn it on.)

2) This is more a curiosity regarding the shell boundary conditions. If I understood correctly from the theory manual, this formulation accounts for the rotation by the relative positions of nodes below and above the mid surface of the shell (front and back). Does this mean that this formulation could model a boundary condition of fixed displacements and fixed rotations by fixing both front and back displacements?

Yes, that's correct. Fixing (x,y,z) and (sx,sy,sz) fixes translations and rotations.

Also, would it make sense to have the fixed displacements BC be that the average displacement between the front and back nodes be zero?

I tried that in FEBio 2.6, by setting the shell degrees of freedom to mid-shell displacements and mid-shell directors (director = displacement of front - displacement of back). The problem with this approach is that one cannot capture correctly the effect of a pressure applied on the front (or back) face, since the pressure is really being applied at the mid-surface. So the shell would not change its thickness in response to such a pressure and I didn't like this non-physical behavior. This became problematic when a shell would be sandwiched between two solid elements, as we may commonly encounter in biomechanics (see http://biomechanical.asmedigitalcoll...icleid=2696682)

Best,

Gerard

[MiguelA]

Dear Professor Ateshian,
thank you for the replies.

I fear that the issues of locking reach a complexity level that far exceeds my capabilities! Based on my small knowledge, I believe that there is always locking in displacement based elements, but higher order elements show significantly less locking so probably the accuracy gains of using 3-field might not justify the costs. If I recall correctly, an alternative to 3-field for lower order quads and hexes would be reduced integration. I don't know if this would make sense for FEBio but there is an interesting technique called F-bar[1]. The same authors also propose a similar technique for tris and tets but it is a little more complicated because it requires the computation of the average J in a patch of triangles around the element[2].

As an aside, while refreshing my memory about locking I saw a paper about "generalized B-bar" [3]. This is interesting because it allows the FEM to avoid locking in any "rigid" direction, not just volumetric. They seem to be concerned with glass fiber composites, but perhaps connective tissue would also benefit from this treatment. Just though it was interesting and connected to this conversation.

Thank you for the clarification regarding the shell boundary conditions. It makes a lot of sense. I was confused with Fixed Displacement vs Fixed Shell Displacement since I didn't know about the front/back thing.

All the best,
Miguel


[1] de Souza Neto, E. A., Perić, D., Dutko, M., & Owen, D. R. J. (1996). Design of simple low order finite elements for large strain analysis of nearly incompressible solids. International Journal of Solids and Structures, 33(20-22), 3277-3296.
[2] Neto, E. D. S., Pires, F. A., & Owen, D. R. J. (2005). F‐bar‐based linear triangles and tetrahedra for finite strain analysis of nearly incompressible solids. Part I: formulation and benchmarking. International Journal for Numerical Methods in Engineering, 62(3), 353-383.
[3] Oberrecht, S. P., & Krysl, P. (2016). Generalized B‐bar: unlocking multiple compliance modes of finite elements with stiff‐fiber anisotropy. International Journal for Numerical Methods in Engineering, 108(6), 568-587.