View Full Version : 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.
Jun
===========
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 )
C
include 'vaba_param.inc'
C
dimension coordMp(nblock,*), charLength(nblock), props(nprops),
1 density(nblock), strainInc(nblock,ndir+nshr),
2 relSpinInc(nblock,nshr), tempOld(nblock),
3 stretchOld(nblock,ndir+nshr),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
9 defgradNew(nblock,ndir+nshr+nshr),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
dimension devia(nblock,ndir+nshr),
1 BBar(nblock,4), defgradNewBar(nblock,5), intv(2)
C
character*80 cmname
parameter (zero = 0.d0, one = 1.d0, two = 2.d0, three = 3.d0,
* four = 4.d0)
C
intv(1) = ndir
intv(2) = nshr
C
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
C
C10 = props(1)
D1 = props(2)
C C10=1.11619E6 D1=4.48E-8
C
ak=two/D1
twomu=four*C10
amu=two*C10
alamda=(three*ak-twomu)/three
C
C if stepTime equals zero, assume pure elastic material and use initial elastic modulus
C
if(stepTime .EQ. zero) then
do k=1,nblock
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)
end do
else
C
C JACOBIAN AND DISTORTION TENSOR (F is non-symmetric)
C
do k=1,nblock
det=defgradNew(k, 1)*defgradNew(k, 2)*defgradNew(k, 3) -
1 defgradNew(k, 3)*defgradNew(k, 4)*defgradNew(k, 5)
scale=det**(-ONE/THREE)
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
C CALCULATE LEFT CAUCHY-GREEN TENSOR (B is symmetric)
C
BBar(k,1)=defgradNewBar(k, 1)**two+defgradNewBar(k, 4)**two
BBar(k,2)=defgradNewBar(k, 2)**two+defgradNewBar(k, 5)**two
BBar(k,3)=defgradNewBar(k, 3)**two
BBar(k,4)=defgradNewBar(k, 1)*defgradNewBar(k, 5) +
1 defgradNewBar(k, 2)*defgradNewBar(k, 4)
C
C CALCULATE STRESS tensor
C
TRBBar=BBar(k,1)+BBar(k,2)+BBar(k,3)
EG=two*C10/det
PR=TWO/D1*(det-one)
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
stressNew(k,4)=EG* BBar(k,4)
C
write(6,*) totalTime,stressNew(k,1),
1 stressNew(k,2),stressNew(k,3),stressNew(k,4)
end do
end if
return
end
Hello Jun,
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 downloads section (http://www.polymerfem.com/modules.php?name=Downloads&d_op=viewdownloaddetails&lid=33&ttitle=vumat_nh.f).
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, isn't it? So why don't 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.
http://www.imagestation.com/picture/sraid104/p26a2a9915ad131a3e4c0da25843f2ad6/f98ad621.jpg
Here are a few suggestions/comments:
:arrow: What value do you use for the bulk modulus? In explicit simulations the bulk modulus cannot be too large
:arrow: Did you make the time increments too large? As you know, in explicit simulatins the maximum time increment has to be rather small.
: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?
: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?
: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?
Dr. Bergstrom,
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 WINDOW'S ADDR BAR.
Mises stress contour comparison for uniaxial tension (upper:VUMAT lower:Built-in NH) difference is apparent.
http://www.imagestation.com/picture/sraid104/pc2f907d9610db04296eec4eac713533d/f98499df.jpg
S11 history plot comparison (Red:NH Green:VUMAT)
http://www.imagestation.com/picture/sraid104/p4a21c0b8f6da1322e2a7fdd5fca68a6f/f98499e6.jpg
Strain tensors
http://www.imagestation.com/picture/sraid104/p8ee5be7725e8254a6ba4a54ea7bdbdd1/f98499e9.jpg
They do have the same trends but there is difference. While in one element test, 2 curves are overlapping together.
Very interesting comments. I have not been able to look at your images, for some reason I don't have permission to reach the jpg-files.
Here are a few more questions/comments:
(1) What density did you use?
(2) You said that you got perfect agreement for the 1 element case. What loading conditons did you try? Did you try both uniaxial and shear?
(3) I am surprised that you did not get the same answer for the 4-element case. Did you try uniaxial loading for the 4-element case? Do you get the same answer as for the 1-element case?
It sounds like you pretty close to getting it right :)
The density I am using is 1104.43 kg/m3
I run the 1 element case again, and found consistent results for uniaxial tension. For some reason I thought I've got fit results for simple shear, but this run shows there still is problem with simple shear of 1 element.
So I should not conclude "1 element test is good". Sorry for that.
I am sending bunch of plots to your email box with the cae file. NH.for is the NH model VUMAT subroutine downloaded from this website.
Thank you again for taking a look!
Jun
Based on your comments I wanted to make sure that the vumat worked as I though that it should. To verify that this was the case I ran a number of validation experiments with the VUMAT and I got the following results:
:arrow: For uniaxial tension, I get very similar results using these methods:
* Closed-form analytical solution,
* ABAQUS/Explicit with one 3D element using the vumat,
* ABAQUS/Explicit with one 3D element using the built-in NH model,
* ABAQUS/Explicit with one axisymmetric element using the vumat,
* ABAQUS/Explicit with one axisymmetric element using the built-in NH model,
* ABAQUS/Explicit with 4x4x4 3D elements using the vumat
See the following figure
http://polymerFEM.com/polymer_files/vumat_nh/tension.png
:arrow: For simple shear, I get very similar results using these methods:
* Closed-form analytical solution,
* ABAQUS/Explicit with one 3D element using the vumat,
* ABAQUS/Explicit with one 3D element using the built-in NH model,
Due to the deformation state, I get a slightly different stress response for the following two methods:
* ABAQUS/Explicit with 4x4x4 3D elements using the vumat,
* ABAQUS/Explicit with 4x4x4 3D elements using the built-in NH model
Although these two methods are different than the one element case, they are both consistent.
See the following figure
http://polymerFEM.com/polymer_files/vumat_nh/shear.png
These results seem to indicate to me that the vumat is working correctly :D
Three questions here:
1) What are the results for plane strain elements e.g. CPE4R as I used?
2) How are the results compared if strain goes up to, say 500%? Rubber normally has very high stretch ratio before breaking.
3) Do you need to add on a pure elastic code just for the first time increment as suggested for VUMAT in ABAQUS manual?
Answers:
1) The results should be qualitatively similar for plane strain elements. The call structure for axisymmetric and plane strain elements are the same with regard to the vumat (ndi=3, nshr=1).
2) There should be no problem to go to 500% strain. The vumat will behave the same.
3) As far as I know, it it not necessary to make the vumat behave in a certain way for the first time increment. The given vumat works as it is.
If you want, you can try out the response at larger deformations, etc.
Best of luck.
Thank you, Dr. Bergstrom, for helping me on this issue these days.
Btw, could you send me the cae file or inp file for the testing you did on this VUMAT?
You can download the input files here (http://www.polymerfem.com/modules.php?name=Downloads&d_op=viewdownloaddetails&lid=34&ttitle=vumat_nh_inputfiles.zip).
Note that the input files are simple examples aimed at quickly testing some aspects of the accuracy and validity of the VUMAT subroutine.
Could you post the closed-form formulations as well? Sorry for asking too much.
Just run the case01_ABAx.inp, why the volume of the 3D element is increasing? Elastomers are imcompressible in most time, arn't they?
Should the element contracts in 2 & 3 directions under uniaxial tension loading along 1 direction?
The closed-form solution was obtained from a specialized computer program that I am working on. This program is still being developed and I cannot disclose more information about it at this time. For the NH-model, however, it should be easy enough to derive the response in uniaxial tension and simple shear.
For test case case01_ABAx.inp the volume is slightly increasing due to the chosen material parameters. This volume increase is cased by the hydrostatic tension component of the stress state. The amount of volumne change is exactly as expected.
I have a question here, is that Neo-Hookean model is suitable in predicting rubber under compression?
I was given a set of data which provide me force-deflectin data. The rubber is cylindlicar shape, height=0.925 inch, Diameter = 0.984 inc, shape facter =0.2659. Initial young modulus under compression about 2050 ibs, poision ratio take 0.4997 (by guess). The rubber be compressed 0.118 inch (3 mm), 12.77% strain. When I using Ansys, neo-hookean model to compare the test data i given, that is totally unmatch with the force-displacement curve unless first 0.1-0.3 % strain. Afther that the force shoot up without stoping as the stiffness almost not softern but the given test data actually sorften after 2% of compression.
PS: I am not given any other data for using orther material model but just Neo-Hookean. And is that anybody can tell me why this model so unmatch with this case?
The Neo-Hookean material model should work reasonable well for strains up to 12%.
How did you decide which material parameters to use?
Which parameter to use?
If you mean the material parameter need to use in Neo-hookean in ansys it was already set by the program. For neo-hookean it need intial shear modulus and d, which equal to 2/K, where K is bulk modulus.
Up to 12%.... is that apply to compressin also or just tension?
With the rubber disk data given by vender, it is pretty high in initial young modulus under compression. I convert it into initial modulus with
G=E/2(1+2r) and bulk modulus K=E/3(1-2r) and r is potion ratio be set to 0.4997. It end up with the curve almost a straight line with only a bit curving and only match in begining about 0.01% strain only. The rest it just keep on going with almost linear but the graph given by vender actually soften in stiffness and two curve become further and further away.
Thank you for reply
The Neo-Hookean material model should work relatively well for the small strains that you are interested in, and the approach that you are taking seems correct. I cannot say why you don't get a better aggrement :?
http://polymerFEM.com/polymer_files/vumat_nh/shear.png
Dear Jorgen,
Do you know why the stress-strain curves are not smooth for vumat 4x4x4 elements shear tests? is it related to the explicit solver?
Thanks,
Wei
I have not really though about that. Perhaps it is the explicit solver, perhaps there is an instability in the incorporation in the boundary conditions and the out-of-plane displacement. Do you have any other ideas?
I ran two simulations, one with umat, one with vumat for the same material model, with same boundary conditions. The results from the umat one showed a smooth line and the vumat one followed the same curve, however, not smooth. I am not sure what is the reason.
I have been able to reduce the oscillations of the VUMAT output by avoiding the *variable mass scaling, and by changing the *density value only. I would guess it might provide uniform kinetic energy for all elements. It is for quasi-static simulations only.
I ran two simulations, one with umat, one with vumat for the same material model, with same boundary conditions. The results from the umat one showed a smooth line and the vumat one followed the same curve, however, not smooth. I am not sure what is the reason.
vBulletin® v3.7.3, Copyright ©2000-2008, Jelsoft Enterprises Ltd.