TensorMechanics: Incorrect stress and strain from ABAQUS UMAT interface

[WulfHans]

Bug Description

I did a simulation of a simple shear test with the UMAT NeoHooke implementation provided by ABAQUS. The expected stresses can be easily computed or just found on wikipedia.

The resulting shear stress was slightly incorrect. The difference also depends on the time step size. The reason is that the interface rotates the stress that comes out of UMAT with the rotation_increment. This is incorrect, as the UMAT returns the cauchy stress and does not require any further modification. However, if this stress is passed back into UMAT during the next time step als stress_old, it must be rotated. Solution: The rotation must be done before calling UMAT.

The normal stress are also incorrect, in fact cauchy_xx and cauchy_yy are swapped. The reason is: C++ uses row-major, while fortran uses column major ordering. Consequence: If a 3x3 matrix is passed from C to Fortran as an array of 9 numbers, the matrix is effectively transposed just by the function call. Inside UMAT, the deformation gradient and all strain tensors computed from it are already wrong. In the case of a simple shear test, this results in swapping the xx and yy stresses. Solution: Transpose all non-symmetric second order tensors before passing them into UMAT.

I also tried to compare with the Moose NeoHooke implementation. Unfortunately, the material law implemented by ComputeNeoHookeanStress is not the same as NeoHooke in commercial programs like ABAQUS. Thats why I have attached our own implementation of NeoHooke, that matches the ABAQUS NeoHooke.

Steps to Reproduce

Compile a test app with the supplied NeoHooke class and the supplied Neo_Hooke.f plugin. Run the attached input file (maybe change the path of the plugin). Compare the stresses. demonstrateUmatBug.zip

Impact

AbaqusUMATStress yields incorrect stresses and should not be used currently.

Fixing the issue

I already have a corrected version of AbaqusUMATStress available. As soon as this is acknowledged as a bug, I would like to supply the fix myself. I just need to setup the workflow for contributing first.

[dschwen]

@WulfHans thanks for you investigation here. @recuero and I had a discussion about this and think there are some things that need clarification.

The resulting shear stress was slightly incorrect. The difference also depends on the time step size.

The shear stress as a function of the deformation state should not depend on the time step size, so your observation is a bit confusing. Note, that MOOSE is using the Rashid formulation for incremental kinematics (International Journal for Numerical Methods in Engineering, 36(23):3937–3956, 1993) rather the Hughes-Winget formulation. That would explain small differences (increasing with rotation and time step) between MOOSE and Abaqus. A user (@jessecarterMOOSE ) recently implemented the Hughes-Winget formulation and reported agreement with Abaqus (maybe he can comment here, too).

The reason is that the interface rotates the stress that comes out of UMAT with the rotation_increment. This is incorrect, as the UMAT returns the cauchy stress and does not require any further modification. However, if this stress is passed back into UMAT during the next time step als stress_old, it must be rotated. Solution: The rotation must be done before calling UMAT.

We can see that this could be a potential issue and would be happy to check such a change agains some internal verification cases.

The normal stress are also incorrect, in fact cauchy_xx and cauchy_yy are swapped. The reason is: C++ uses row-major, while fortran uses column major ordering

This is odd. How would a transpose (or lack thereof) swap component indices this way (and ultimately on-diagonal components)? We do have tests with orthotropic elasticity tensors, that exhibit shear anisotropy. And there a UMAT implementation reproduces the results of the native MOOSE implementation exactly.

...unless our testing UMAT was accidentally written with a row-major ordering... (we'll double check that)

Anyway, always great to see users providing feedback. Would you mind sharing which institution/company you work at?

[WulfHans]

I realized, my input file already uses my local version of the UMAT interface (AbaqusUMATStressCustom). Just delete the Custom and it should work. If you want to have a look at it: Here is my modified version that I believe fixes the issues. AbaqusUMATStressCustom.zip

The shear stress as a function of the deformation state should not depend on the time step size, so your observation is a bit confusing.

I have seen that most testing Umats use the incremental formulation, whereas my umat directly computes the stress from the deformation gradient. Perhaps this contributes to hiding the problem. For the direct stress computation, the reason is straightforward: The stress tensor gets an additional multiplication with rotation_increment. The smaller the time step, the smaller this additional rotation and the smaller the error.

The normal stress are also incorrect, in fact cauchy_xx and cauchy_yy are swapped. The reason is: C++ uses row-major, while fortran uses column major ordering

This is odd. How would a transpose (or lack thereof) swap component indices this way (and ultimately on-diagonal components)?

Its easiest to just compute: If F is [1,1,0; 0,1,0; 0,0,1] (shear F_xy= 1), then C++ will store it in exactly this order in memory, but Fortran will read it as [1,0,0; 1,1,0; 0,0,1] (shear F_yx= 1). The left Cauchy-Green tensor b is defined by b = F * F.transposed(). So on C++ side we get [b] = [2,1,0; 1,1,0; 0,0,1] whereas Fortran computes [b] = [1,1,0; 1,2,0; 0,0,1]. The shear strain component b_xy is actually correct (which could be a reason why some tests pass). But the components b_xx and b_yy are switched and so will be the stress components. This switching of xx and yy is a specialty of the simple shear test. If you pass in an arbitrary unsymmetric F, the strain will be wrong in a more complex way. The direct method to confirm the issue is: Debug into the UMAT or just add some write statements in there. You will see, that the F(1,2) component passed from moose will end up in the F(2,1) spot in UMAT.

We do have tests with orthotropic elasticity tensors, that exhibit shear anisotropy. And there a UMAT implementation reproduces the results of the native MOOSE implementation exactly.

The problem is: If the deformation gradient F is symmetric, you will not see any of the problems. There will be not rotation_increment, and transposing F before passing to UMAT does not change anything obviously. This includes uniaxial, biaxial and volumetric tests. In other words: If you have some orthotropic material, where the orthotropic directions are aligned with x,y,z and you pull it in x,y,z you will nicely see the orthotropy and all stresses will be correct.

...unless our testing UMAT was accidentally written with a row-major ordering... (we'll double check that)

I was also unsure while reporting, if I am really right on this, because its difficult to believe that such an error slips trough the net of tests. I have started to investigate, which tests actually create a unsymmetric F. What I have found so far:

print_shear.i somehow investigates the shear component of the strain. I admit, I have not really understood, what it does so far and why it passes. Perhaps due to the fact, that the shear strain actually comes out correctly.

The combination of the elastic_print.f and elastic_print_c.C plugins definitely has the potential to highlight the problem. But I have not found an actual test, that compares these two under shear loading. I will see, if I can put together such a test.

shear_order.i uses shear, but the magnitude of deformation is very small (1e-5). Recall: The small strain tensor can be defined by symmetrizising (F-I), which does not care about transposing F or the rotations. Hence, as long as the nonlinear strain tensors are very close to the small strain tensor, you will hardly see any of the issues. The difference between small and large strain is a second order effect, something in the magnitude of 1e-10 in this case. This is easy to miss.

Anyway, always great to see users providing feedback. Would you mind sharing which institution/company you work at?

I work as postdoc at TU Chemnitz (Germany). We are just starting using Moose, but I hope it will be the basis for our future FEM software projects also hope to contribute here and there.

[jessecarterMOOSE]

FWIW, there's a commit in my fork (jessecarterMOOSE/moose@f3ef7260cafd189935266becfd29f7bf6c65521c) that takes a stab at doing Hughes-Winget kinematics for comparisons with Abaqus (I piggy-backed on the existing Rashid decomposiiton method enum, which probably isn't the cleanest thing to do). Can't say my analysis was extremely thorough, but I did use a few different UMAT's to test. The comparison between MOOSE and Abaqus was done by a rather crude method of printing out the return values from inside the UMAT before exiting the UMAT call (so I didn't have to postprocess with Abaqus CAE which only uses single precision). For a single-element model undergoing simple shear, the two codes agreed to within about the precision of the solver (e.g. 1e-10). This was true for all quadrature points on the element and for all time steps.

[dschwen]

Hans, could you take a look at https://github.com/jessecarterMOOSE/moose/commit/f3ef7260cafd189935266becfd29f7bf6c65521c . This changes the sequence of rotations. In your proposed modification I see that you commented out the stress rotation after the UMAT call, but didn't add a rotation of the old stress going into the UMAT (which Jesse has in his Hughes-Winget patch).

Anyways, could you go ahead and open a pull request for your change? Let me know if you need any help with that. We have a bit of a time difference, but if you want to have a video call to get a contribution workflow set up let me know. Anything after your 4:30pm could work for me. Und wir koennen uns dann auch auf Deutsch unterhalten wenn das einfacher ist.

[jessecarterMOOSE]

After talking with my colleague a bit on this, seems like for my test cases, I'm using the Hughes-Winget strain increment directly, which is symmetric, meaning my test cases wouldn't be affected by any alleged symmetry bug. With a Neo-Hookean model, one must compute the strain in the UMAT based on the deformation gradient, so I can understand why it must be right.

Regarding the tensor ordering issue, UMAT specification requires tensors to be passed in as arrays, which the AbaqusUMATStress class does, right? And then the user is expected to re-build the tensors themselves in the UMAT from those arrays? The order of an array is not ambiguous, so can the Abaqus manual be consulted to check for the proper ordering?

[dschwen]

The order of an array is not ambiguous, so can the Abaqus manual be consulted to check for the proper ordering?

The Abaqus manual is not explicit here, but all examples are fortran code and the fortran ordering is column major. Hans is correct there.

Our RankTwoTensor is derived form libMesh::TypeTensor, and that class is indeed row major. Ouch! That's quite an oversight...

[WulfHans]

Minimal proof for the the issue (btw: It would be really useful, if you could attach .i files directly): UMAT_test_min.zip

The prescribed boundary condition is shear F_xy=0.5 Your can switch around the plugin (elastic_print and elastic_print_c) and see the difference.

The C plugin gives:

dfgrd1 0 = 1.0000000
dfgrd1 1 = 0.5000000
dfgrd1 2 = 0.0000000
dfgrd1 3 = 0.0000000
dfgrd1 4 = 1.0000000
dfgrd1 5 = 0.0000000
dfgrd1 6 = 0.0000000
dfgrd1 7 = 0.0000000
dfgrd1 8 = 1.0000000

The Fortran plugin gives:

 DFGRD1_ 1 1    1.0000000
 DFGRD1_ 1 2    0.0000000
 DFGRD1_ 1 3    0.0000000
 DFGRD1_ 2 1    0.5000000
 DFGRD1_ 2 2    1.0000000
 DFGRD1_ 2 3    0.0000000
 DFGRD1_ 3 1    0.0000000
 DFGRD1_ 3 2    0.0000000
 DFGRD1_ 3 3    1.0000000

I have forked the repo now and started working on the fix - can you assign it to me please?

The first thing I tried was just copying my local fix in and check if any tests now fail. And indeed one does:

test:umat/print.print_shear_print ........................................... FAILED (EXPECTED OUTPUT MISSING)

I believe this is one of the rare cases, where the test is wrong and not the program. Perhaps it was designed with some wishful thinking. This test is somehow strange anyways - it uses some strain-dependent material parameters and claims to test proper passing of the parameters as a side effect. Question: Is it my job to fix this test? Because this could be rather annoying.

I would much prefer to create a fresh test, based on my investigations. Another question: Should I try to resuse the existing .f plugin, or can I add my own NeoHooke.f ? I believe this example is a classic and should be available for tests anyways.

[dschwen]

Minimal proof for the the issue (btw: It would be really useful, if you could attach .i files directly):

The workflow here would be to start a pull request that adds the minimal non-working example as a new test that will right away fail in our CI system. And as you provide the fix will then start passing.

I have forked the repo now and started working on the fix - can you assign it to me please?

Done

test:umat/print.print_shear_print ........................................... FAILED (EXPECTED OUTPUT MISSING)

This is expected. That test just checks the value of all components to guard against a future regression. I can help you regold it once your PR is up.

I would much prefer to create a fresh test, based on my investigations. Another question: Should I try to resuse the existing .f plugin, or can I add my own NeoHooke.f ? I believe this example is a classic and should be available for tests anyways.

Feel free to add more tests and if you can share the NeoHooke.f plugin under LGPL that would be a welcome addition as well.

[WulfHans]

The workflow here would be to start a pull request that adds the minimal non-working example as a new test that will right away fail in our CI system. And as you provide the fix will then start passing.

Done. Should I already PR the fix as well?

This is expected. That test just checks the value of all components to guard against a future regression. I can help you regold it once your PR is up.

I have a fix for this test. Simply transposing the offending tensors in the tests file works.

Feel free to add more tests and if you can share the NeoHooke.f plugin under LGPL that would be a welcome addition as well.

This is tricky - this UMAT is one of the standard examples and more or less everybody has it. If you just google "how to write UMAT", you will find this (or the incremental version already in git under the name elastic.f). But I have no idea if it has been officially released with any license. Its not my personal code (except for some small changes). But there is also not much you could do differently, as it straightforward computes the outputs from the inputs. Is that a problem? In my opinion, if the current elastic.f was admitted this should as well.

[WulfHans]

PR contains the fix for the transposing issue.

Concerning the rotation problems: I will put this in a new branch and PR, so you can review accept it separately. Good idea? For the test, I have decided to use the Neo_Hooke UMAT, that computes the stress from the strain directly. Reason: The incremental formulation (elastic.f) is dependent time step size in either case. For this model, you can also find a series of deformations which end in the undeformed state, but it still computes a significant non-zero stress. In other words: This is perfect for highlighting the design flaws in Abaqus, but it does not behave how Hyperelasticity should behave. This is not a proper basis for this test, thats why the other NeoHooke implementation is used.

The test is already up, but the fix breaks several existing tests again. This is not very surprising, but needs more work.

[WulfHans]

PR with fix for the rotation problem is up.

The success with fixing the other tests (which now partially fail, no surprise) is mixed: umat/print.print_shear_print is easy to fix, just provide the correct outputs umat/print_c is also easy, C and Fortran still yield the same output, you just need to update the gold file

umat/print.print_shear and umat/predef.predef_multiple are nasty. They are both based on the same strange test, where the stiffness is scaled with some arbitrary factor, that depends on some strain components. What happens in the simulation, is that first strain_xx becomes inhomogeneous because of the fixed support at the bottom, the slightly different stiffness creates a gradient in strain_yy and then the element starts to shear. These make problems, the MOOSE computation and UMAT dont do exactly the same.

I tried to investigate, the source of the differences. It roughly goes:

• even before the first linear step, there is already a very small difference in the passed DROT. For some reason, it has off-diagonal elements different from zero (magnitude 1e-9) and differences in the magnitude 1e-16
• then, the first linear step has a difference in |R| of about 1e-10 (|R| about 1e-7). Seems that the first step is already slightly different
• then, slightly different strain increments are passed, yielding different stress. The relative difference is only about 1e-8 (so the first 7-8 digits are the same), but its enough
• from here, everything slowly but surely diverges, at some point one computation does an additional linear step, the other one does not and it all falls apart

Overall: The correct way I see to debug this is to compare both computations down to the quadrature point and find the difference, but this is very tedious and frustrating. In addition, I do not know if the native MOOSE computation is really guaranteed to do exactly the same (the reference) as the UMAT construct in this scenario. I also do not trust this test in general, as it managed to pass serious errors as correct in areas it claims to cover. Thats why I stopped here for now.

Final remark: I will be on vacation for two weeks and I will not answer or provide anything here. But I have not abandoned the issue and will be back afterwards.

[jessecarterMOOSE]

Rotating the old stress before the call to the UMAT does get you closer to Abaqus, but the rotation increment which you are rotating by is calculated differently for Rashid (MOOSE) and Hughes-Winget (Abaqus). That's why I had to touch on other files in my branch mentioned above to get Abaqus and MOOSE to match.

[jessecarterMOOSE]

This was fixed #21872

[recuero]

Thanks, @jessecarterMOOSE. Closing it.