About Problem in UMAT of Drucker Prager

[willness]

Hi, dear Jorgen and all

I wrote a subroutine of Drucker Prager constitutive model (four material parameters, E, v, fai, fi) and want to incorporate it into ABAQUS. The subroutine was written by following the Back-Euler return integration shown in the Book "Nonlinear Finite Element Analysis of Solids and Struction" (by Crisfield).

It works well when I verify it by simulating trixial compression test (all the stresses are negative) with one element. However, when I verify it by simulating trixial extension test, it doesn't look well. The analysis cannot converge. Is there any difference in the Back-Euler procedure when the stresses turn from negative to positive? Because when we look at the yield function (alfa)*(I)+sqrt(J)=0, as stresses turns positive, this equation may not be satisfied. Or could you please suggest a few books in this topic.

If you are willing, I am happy to post the subroutine and the *.inp file so that you can have a quick check if you can use ABAQUS.

Thanks a lot for your concern and help!

willness

[Jorgen]

Hi Willness,

How do you like Crisfield's books? I have not read them yet.

Have you tried to use the built-in Drucker-Prager model with your input file? That could be a good first check. If you want, I can take a quick lock at your code.

- Jorgen

[willness]

Many many many thanks to you, Jorgen!

I do like Crisfield's book, possiblely because there are some numerical examples. For those who are not very strong in mechanics and mathmatics, like me, these examples provide an effective way to help understand the theory.

OK, may I turn to the problem I encounted. Actually, I have done what you suggested to do, although not used Drucker-Prager, but Mohr-Coulomb model built in Abaqus. I did found something meanful to me. When I used built-in Mohr-Coulomb model, generally, the S22 will not larger than a certain value (for the example in the attached *,inP file, S22 will not larger than around -33.***kPa). However, for my own subroutine, this value is larger (around -13.**kPa), and there is a warning "THE SYSTEM MATRIX HAS 1 NEGATIVE EIGENVALUES". I think the problem occours here, although analysis may go a little further, and finially the force cannot converge.

I suspect that in the built-in model, there is some kind of trick (like tension cuf-off?) I don't know.

Attached are the Fortran code and Back.inp. I am very gratful if you are willing to spent some time to have a quick look. I wrote many "write" in the code, so that after each step, there is some information printed out to *.msg file.

Looking forward to your commments!

By the way, the *inp file is from Abaqus standard example problems manual Vol.1 2.2.7. I changed the CamClay model to DP. The code attached is just for verification the Back-Euler procedure. The plastic strain is not stored. Thethe fifth parameter and statev is not meanful. No hardening considered.

[willness]

I found that when the traction displacement is not very large, although there is warning message of NEGATIVE EIGENVALUES(this message emerges at where plastic strain happens), the analysis still can be finished if the load step is not very small. It looks like that in the computation, there is error, but the error is not significant sufficiently to affect the analysis if the load step and the magnintute of the load is controlled. Otherwise, the error will accumulate significantly to affect the convergance of the force or displacement.

However, this problem doen't occour when the load is a contraction displacement, no matter how much the displacement is and how small the load step is.

I am wondering what is the direction of a way to solve the problem. Is it possible there is mistake in the direction of the tensor in calculation?

[Jorgen]

I am afraid that I will not have enough time to carefully study your user-subroutine :(

My recommendation is that you keep printing out intermediate results from within your code and then validate those results with hand calculations...

Best of luck,
Jorgen

[willness]

Hi, Jorgen, thanks for your kind suggestion.

Frankly, I am doing what you suggested. I used MATLAB to help check the stiffness matrix and other vectors printed out. When plasticity happens, the stiffness matrix do have a negative eigenvalue (very close to zero). But I don't know exactly what cause the problem (I know roughly that non-associated flow rule leads to unsymmetric stiffness matrix and probabily causes the problem) and how to solve it.

I am afraid I have to read some papers on return algorithm for non-associated plasticity, becasue I am afraid the theory that the code based on is not suitable for the non-associated plasticity.

Any further recommendation, please do leave a message here.

Thanks again.

Best regards,
Willness

[willness]

I doubt something not correct in the updating of Jacobian matrix. I found a big surprise. For the simple loading condition such as trixial compression or extension, the code works if I use the elastic matrix directly even if the plastic deformation is very large. For complex loading (say the shear is significant), it doen't work.

When I use build-in Mohr-Coulomb model, there is warning message of "Negative Eigenvalue" too. However, the analysis can be finished.

What do you think of my problem based on the phonemon described above? Do you think it is caused by the non-associated flow rule (so that the [K] matrix (negative eigenvalue))? Could you please suggest some papers or books on this topic if you know?

Thanks a lot!

Best regards,
willness