Issue with ComputeFiniteStrain - seems unable to solve problem with zero rotation

[richmondodufisan]

Hi, I am solving a simple problem which involves stretching of a membrane in one direction. I get the following error with ADComputeFiniteStrain:

Cannot take square root of a number less than or equal to zero in the calculation of C1 for the Rashid approximation for the rotation tensor. This zero or negative number may occur when elements become heavily distorted.
Nonlinear solve did not converge due to DIVERGED_FNORM_NAN iterations 0

However, if I switch to ComputeIncrementalSmallStrain it works fine. I initially thought the issue was with the stretch ratio of 1.1 (i.e 1.1 times it's original size) being to small to be considered Finite Strain, but then when I increased the ratio to 1.8, I got the same error. I am using an incredibly fine mesh, and since the poisson's ratio is close to 0.5, I also applied a volumetric locking correction.

I am not sure how to solve this issue. There is no rotation in my problem so the rotation tensor should be zero.

[GiudGiud]

@gambka @jiangwen84

[dschwen]

ComputeFiniteStrain - seems unable to solve problem with zero rotation

If "no rotation" really was the issue then we wouldn't even be able to compute the initial residual or a body at rest with no loading. I'm pretty sure that's not the issue here. It would help if you added a minimal non-working example here.

[richmondodufisan]

I see. Here's the input file with the issue:

https://github.com/richmondodufisan/purple/blob/main/Cornea_Linear_Elastic/Cornea_Stretch.i

[richmondodufisan]

Hi @dschwen , were you able to take a look? I haven't been able to figure out the reason for the error

[gambka]

What order elements are you using? The volumetric_locking_correction is not required for higher (second) order elements. What happens if you set the decomposition_method = EigenSolution instead of the default TaylorExpansion?

[richmondodufisan]

I was using using first order elements. When I switched the decomposition_method, it didn't converge, the non-linear residual stayed at the same order of magnitude. I tried switching to 2nd order elements (and I took out the volumetric_locking_correction) but it still did the same thing. Here is the input file:

https://github.com/richmondodufisan/purple/blob/main/Cornea_Linear_Elastic/Cornea_Stretch.i

[richmondodufisan]

Hi @gambka , were you able to take a look? I haven't been able to figure out the issue yet

[hugary1995]

A few suggestions:

1. Ramp up displacement bc with time -- allow the solve to fail and time step to be cut.
2. Set preset=false in the bcs.
3. Turn off line search.

[richmondodufisan]

@hugary1995

Hi, I implemented with the new system. I couldn't get it to converge without restraining the top and bottom boundaries.

-The "base_name" functionality didn't work, I had to hardcode the name of the elasticity tensor in the input file
-I set large_kinematics = true globally
-I also set stabilize_strain = true globally, for volumetric locking correction
-It didn't converge, so I tried switching to second order elements but it also did not converge

In all simulations, the residual would decrease and then remain constant for several iterations until the "solve did not converge" error shows up. Then I tried setting the displacement smaller, and the lower the applied displacement, the lower (and closer to convergence) the residual would get. When the displacement became small enough (1%, i.e stretch_ratio = 1.001), it converged as it was essentially a small strain problem.

Does this mean that I cannot get a convergent result unless I use a hyperelastic model? I know that for large stretches, using a linear elastic model is unphysical but I didn't expect it to not converge at all. Also, looking at the dimensions of modules/solid_mechanics/test/tests/lagrangian/cartesian/total/patch/large_patch.i, the applied displacement seems greater than 1% of the mesh dimensions.

Here is the relevant input file: https://github.com/richmondodufisan/purple/blob/main/Cornea_Linear_Elastic/Cornea_Stretch.i

---

If this is all true I want to try implementing a hyperelastic model. I read through the documentation but want to ask about my approach before trying. The steps I plan to take are:

-Implement a custom strain energy function, W(F), with a custom material
-Implement a material that inherits from ComputeLagrangianStressPK1, and define the stress as the derivative of W wrt F

-To calculate the derivative, I want to follow the steps in https://mooseframework.inl.gov/source/materials/DerivativeMaterialInterface.html#:~:text=The%20derivative%20of%20a%20material,derivatives%20of%20other%20material%20properties.

In this case, I will "get" the custom strain energy material I created, and differentiate wrt "deformation_gradient" (the declared name of F as seen in ComputeLagrangianStrainBase.C)

Am I missing any steps?

[richmondodufisan]

@hugary1995 did you get a chance to give this a look? your feedback will be much appreciated

[richmondodufisan]

Update:

Implemented the hyperelastic (neo Hookean) model as described, but again, cannot get convergence past 1% strain for this incompressible formulation.

Input file: test_3.i
https://github.com/richmondodufisan/purple/tree/main/Cornea_Hyperelastic_NeoHookean

[richmondodufisan]

Hi, I found the source of the problem. To get it to converge to the right solution I:

-Switched to manual numerical differentiation, was able to get at least one converged timestep. Apparently, the getMaterialPropertyDerivative cannot do tensor differentiation.

-Switched to a penalty method. The PC failed due to FACTOR_NUMERIC_ZEROPIVOT was coming from the extra field (pressure), and there was no boundary condition to set to it. After implementing that that went away.

-Switched to PK2 formulation (probably didnt do as much to help)

-Set low floating point values to zero, they were coming up in the deformation gradient and seemed to throw the solution in the opposite direction? (see the last result above). I'm not sure how but I remember this being the last thing I did before it worked.

All this was because I couldn't get a similar problem to converge for MOOSE's linear elastic material for any deformation more than 1%. Maybe using a nearly incompressible poissons ratio on a linear elastic material really is that stiff even with stabilize_strain?