UMAT and viscohyperelastic constitutive model implementation
There are many questions with regards to the UMAT and the DDSDDE and STRESS computation in ABAQUS. Even though my question might be similar to some of the others, I wanted to go ahead and start a new thread.
I am currently following these steps in UMAT:
1- Computing the STRESS at the end of the current increment (at time t) using deformation gradient, i.e. DFGRD1(i,j).
2- Computing the STRESSINT at an intermediate time, tint, very very close to t.
DDSDDE(1,1) = (STRESS(1) - STRESSINT(1)) / (e1(t)-e1(tint))
DDSDDE(1,2) = (STRESS(1) - STRESSINT(1)) / (e2(t)-e2(tint))
DDSDDE(6,6) = (STRESS(6) - STRESSINT(6)) / (e6(t)-e6(tint)) .
All for a 3D problem.
Unfortunately, things aren't going too well at the moment! From everything I am reading I think my problem is in the second step: computation of the STRESSINT.
I need a few clarifications on:
1- What is the best way of getting the stress at the intermediate step? Better asked, how do I figure out the deformation gradient at this intermediate time?
2- What is the best way to get the strains at this intermediate time as well?
3- If I correctly pick a tint to be very very near t, can I assume linear relationship? Can I do
DFGRD(tint) = DFGRD1(T) - [(DFGRD1-DFGRD0)/(appropriately large const.)] ? If yes, is this valid as well for the strains?
When I do print statements, one issue I am finding is that the DFGRD1 and DFGRD0 are the equal to each other, resulting in equal STRESS and STRESSINT. Any comments on why DFGRD1=DFGRD0.
I am still trying to troubleshoot, in the mean time I really look forward to your responses.
It is OK to assume that the deformation gradient is varying linearly over the time increment. That is, you can estimate DFGRD(tint) by using a linear interpolation. If you get convergence problems then you can try to reduce the max allowed time increment size in order to make the linearization assumption more accurate.
Your scheme to estimate the Jacobian is not a good general purpose approach. The term DDSDDE(1,1) is (basically) the ratio of the change in stress term 1 (sigma11) due to a small change in strain term 1 (strain11). Note that when performing this calculation only one strain term should be varied (in this example strain11). If you simultaneously change all strain terms then you will not be able to estimate the Jacobian terms due to the coupling between the terms.
Thank you very much for the information. So, if I vary e11, e22, e33, ..., one at a time, I can see that the code can become very long. Do you have any suggestions to compute the jacobian?
You are right. Numerical calculation of the Jacobian can make the UMAT run slower. You basically have two options: (1) use a numerical Jacobian; or (2) derive a closed-form approximation of the Jacobian.
Personally, I most often use an approximate closed-form Jacobian when I write UMAT subroutines.
Follow-up with another question
I think that I know have a good handle on the problem. And with that, I feel as though I have a good code. But, I struggling with the following:
I am applying a uniaxial displacement load on a cube. However, the values of the deformation gradient goes much beyond what I am prescribing. I am expecting the stretch ration to be somewhere between 0.8-1.0 and I get 8, 10, or much higher. What would be overriding my prescribed displacements? I check the DTIME parameter, and it seems to be OK.
I look forward to your responses.
That sounds strange.
I recommend the old brute-force method in which you print out a bunch of stuff from within the UMAT until you figure out what is going...
Simple additional question
When computing DDSDDE, for the 1,2,3, I am prescribing uniaxial deformation gradient. For example, for perturbation in E1, I have
F = [L 0 0; 0 1/sqrt(L) 0; 0 0 1/sqrt(L)] where L = 1+E1, the stretch ratio.
Then I compute the stress using this deformation gradient, and
DDSDDE(1,1) = STRESS(1) / E1
DDSDDE(1,2) = STRESS(2) / E1
I repeat this process 6 times for each perturbation to get each row of the Jacobian.
My question is, should I be prescribing a pure or simple shear condition for the 4, 5, 6 perturbations?
You should use a pure shear perturbation.
Another follow-up. I simplified and simplified. Now, I only have a hyperelastic const. model that I am implementing with UMAT. I have already used this model before with UHYPER. So, I am comparing the results I get from UMAT and UHYPER. The ones I get from UHYPER are equal to the theoretical calculation, whereas the UMAT is not so.
This is a 3D incompressible problem with c3d8h elements. The load is uniaxial and the resulting stress field for both UMAT and UHYPER results is uniform and in the direction of the displacement load.
The only issue is that in the UMAT results, the displacements in the 11 and 33 directions are 0. The deformation gradient has [1 lambda 1] as the diagonal terms and not [1/sqrt(lambda) lambda 1/sqrt(lambda)].
I think this is related to the incompressibility constraint somehow, however I thought using hybrid elements was all I needed to account for incompressibility.