help seeking on VUMAT for Neo-hooke model in ABAQUS
Dr. Bergstrom and other experts,
I programmed a VUMAT of neo-hooke model in ABAQUS, according to the example in ABAQUS example problems manual 1.1.14 (UMAT). I used the same constitutive law and formulation, but in the single elament testing phase, it failed due to not agreeing with ABAQUS built-in neo-hooke model with same parameters on S11, S12 etc. distributions while having rotation on the element. Uniaxial tension is fine, through. I post the subroutine code below, if you can give me any lights on it, that will be greatly appreciated. The code is exactly a VUMAT version of the subroutine in ABAQUS example problems manual 1.1.14 bootseal_umat.f.
Thanks a lot.
subroutine vumat (
C Read only -
* nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
* stepTime, totalTime, dt, cmname, coordMp, charLength,
* props, density, strainInc, relSpinInc,
* tempOld, stretchOld, defgradOld, fieldOld,
* stressOld, stateOld, enerInternOld, enerInelasOld,
* tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
* stressNew, stateNew, enerInternNew, enerInelasNew )
dimension coordMp(nblock,*), charLength(nblock), props(nprops),
1 density(nblock), strainInc(nblock,ndir+nshr),
2 relSpinInc(nblock,nshr), tempOld(nblock),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
1 BBar(nblock,4), defgradNewBar(nblock,5), intv(2)
parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
* four = 4.d0)
intv(1) = ndir
intv(2) = nshr
C if (ndir .ne. 3 .or. nshr .ne. 1) then
C call xplb_abqerr(1,Subroutine VUMAT is implemented //
C * only for plane strain and axisymmetric cases //
C * (ndir=3 and nshr=1),0,zero, )
C call xplb_abqerr(-2,Subroutine VUMAT has been called //
C * with ndir=%I and nshr=%I,intv,zero, )
C call xplb_exit
C end if
C10 = props(1)
D1 = props(2)
C C10=1.11619E6 D1=4.48E-8
C if stepTime equals zero, assume pure elastic material and use initial elastic modulus
if(stepTime .EQ. zero) then
trace1 = strainInc(k,1) + strainInc(k,2) + strainInc(k,3)
stressNew(k,1) = stressOld(k,1)
* + twomu * strainInc(k,1) + alamda * trace1
stressNew(k,2) = stressOld(k,2)
* + twomu * strainInc(k,2) + alamda * trace1
stressNew(k,3) = stressOld(k,3)
* + twomu * strainInc(k,3) + alamda * trace1
stressNew(k,4) = stressOld(k,4)
* + twomu * strainInc(k,4)
C write(6,*) totalTime,stressNew(k,1),
C 1 stressNew(k,2),stressNew(k,3),stressNew(k,4)
C JACOBIAN AND DISTORTION TENSOR (F is non-symmetric)
det=defgradNew(k, 1)*defgradNew(k, 2)*defgradNew(k, 3) -
1 defgradNew(k, 3)*defgradNew(k, 4)*defgradNew(k, 5)
defgradNewBar(k, 1)=defgradNew(k, 1)*scale
defgradNewBar(k, 2)=defgradNew(k, 2)*scale
defgradNewBar(k, 3)=defgradNew(k, 3)*scale
defgradNewBar(k, 4)=defgradNew(k, 4)*scale
defgradNewBar(k, 5)=defgradNew(k, 5)*scale
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)
BBar(k,1)=defgradNewBar(k, 1)**two+defgradNewBar(k, 4)**two
BBar(k,2)=defgradNewBar(k, 2)**two+defgradNewBar(k, 5)**two
BBar(k,4)=defgradNewBar(k, 1)*defgradNewBar(k, 5) +
1 defgradNewBar(k, 2)*defgradNewBar(k, 4)
C CALCULATE STRESS tensor
stressNew(k,1)=EG*(BBar(k,1)-TRBBar/Three) + PR
stressNew(k,2)=EG*(BBar(k,2)-TRBBar/Three) + PR
stressNew(k,3)=EG*(BBar(k,3)-TRBBar/Three) + PR
ABAQUS VUMAT for the NH model
Thanks for your request.
Instead of debugging the code that you submitted I created a completely new VUMAT for the neo-hookean model. You can download and study this file in the [url= https://www.polymerfem.com/modules.php?name=Downloads&d_op=viewdownloaddetails&lid=33&ttitle=vumat_nh.f ] downloads section[/url].
Thank you so much, Dr. Bergstrom
Very few advisors would like to write codes for their students and you are one of them.
I have several questions on the coding after reading and testing it in my ABAQUS testing file.
1) F(i,j) calculated in your subroutine is really deformation gradient, isnt it? So why dont you just use the F1(nblock,ndi+nshr+nshr) in the parameter list directly?
Also note that F is unsymmetric tensor, so if you want to calculate it, do you need to calculate F(2,1) as well since it is not equal to F(1,2)?
Since F=RU, do we need to consider the rotation tensor R if rotation can not be ignored? In my subroutine, I just use defgradNew(k, ndir+nshr+nshr) as F, that would be easier, right?
2)F(2,1) is not defined ahead when calculating B(2,2).
3)F is unsymmetric tensor, but the steps calculating B(i,j) depend on a symmetric F.
4)Stress = mu/J * Dev(Bstar) + kappa*log(J)/J * Eye
Stress = mu/J * Dev(Bstar) + kappa*(J-1) * Eye
I even saw a formula
Stress = mu/J * Dev(Bstar) + kappa*log(J)/J *(J-1)* Eye
Which one is correct?
Thanks again for your advising!
Good questions. Here are my answers:
(1) The variable F used in the subroutine is the really the symmetric tensor U. When writing a VUMAT for ABAQUS you need to make sure that you return the stress tensor in an appropriately rotated coordinate frame. One way to do this is to use the deformation F to calculate the stress, and then rotate the stress to the right coordinate frame. Another method is to simply use the tensor U instead of the deformation gradient F.
(2) Since I am using U instead of F, and since U is symmetric, I can simplify the calcuations and only use the upper-diagonal terms of the U tensor (which is called F in the subroutine). You are right, however, that the term F(2,1) was used inappropriately. My mistake, it should have been F(1,2). I have updated the VUMAT routine in the downloads section.
(3) About the volumetric part. There is not one correct version for the dependence on the term J. All three of your expressions can be considered, which one is best depends on the material that you are intrested in modeling. For most elastomeric materials, the volumetric deformations are small and the all three expressions give very similar results. For compressible materials (like foams) the situation is very different.
Thanks for very clear explaination.
I did tests for the NH VUMAT.
Everything looks fine with 1 element tension, biaxial tension, simple shear tests.
However, for 2 and more elements tests, the results from VUMAT are quite different from those of Neo-hookean model in ABAQUS even with same mu and kappa values.
The images I am sending to your email box include results of uniaxial 2 elements test. Upper part has VUMAT, and lower has neo-hooke from ABAQUS. Both used same parameters. Quite confused with the results shown.
Here are a few suggestions/comments:
[list]:arrow: What value do you use for the bulk modulus? In explicit simulations the bulk modulus cannot be too large[/list]
[list]:arrow: Did you make the time increments too large? As you know, in explicit simulatins the maximum time increment has to be rather small.[/list]
[list]:arrow: When you run a one element simulation, do you get the same answer between the VUMAT and the built-in NH model? Do you get the stress and strain response that you expect?[/list]
[list]:arrow: When you run a multiple element simulation, do both the VUMAT and the built-in NH model work? I.e. Do you get deformation fields that you expect? Do you get similar stress fields?[/list]
[list]:arrow: Is there any particular reason you would like to combine elements with the VUMAT and elements with the build-in NH model in the same simulation?[/list]
Thanks again for your kindly reply.
In my subroutine, I am using C10=1.11619E6 D1=4.48E-8 for NH model, which correspond to a bulk modulus kappa=4.46E7 and shear modulus mu=2.23E6. The ratio for kappa/mu=20 which is the suggested value in Explicit.
In my model, time increments are decided by ABAQUS itself. ABAQUS find the most critical element and then decide the max. time increment. I did not specify a value for that. Below is the related lines in .sta file on the max. time inc used.
Most critical elements :
Element number Rank Time increment Increment ratio
1 1 1.603471E-03 1.000000E+00
In one element simulation, I do get the perfectly same stress and strain history curve between the VUMAT and the built-in NH model and the stress and strain response is as expected.
Then I run a 4-element simulation, the VUMAT and the built-in NH model both give some similar results but the difference is large.
Using the VUMAT and the build-in NH model in one model is for benchmark purpose only, these procedures are just used to make sure my VUMAT is working correctly. Then I can have confidence to use the VUMAT instead of built-in NH model in my complicated simulations with element deletion option.
DO NOT KNOW WHY THE PICTURES ARE NOT SHOWN HERE PROPERLY. YOU CAN COPY THE IMAGE ADDR FROM ITS PROPERTIES AND PASTE IT TO A new IE WINDOWS ADDR BAR.
Mises stress contour comparison for uniaxial tension (upper:VUMAT lower:Built-in NH) difference is apparent.
S11 history plot comparison (Red:NH Green:VUMAT)
They do have the same trends but there is difference. While in one element test, 2 curves are overlapping together.