Close this search box.
Clear all

Plane Stress (Shell) incompressible hyperelasticity VUMAT written in C -- Example & Instability problem

3 Posts
2 Users
0 Reactions
Posts: 2
Topic starter
New Member
Joined: 4 years ago

Based on examples from Dr. Bergstorm, I wrote a VUMAT in C to implement Yeoh hyperelasticity for both 3D and plane stress whereas plane stress assumes incompressible material.  The incompressibility is enforced by F33 = 1./(F11*F22-F12*F21) and plane stress by finding duDi1 so that S33=0. 

The 3D version works well and results match ABQ almost exactly including time increment selected by ABQ for both built-in material model and vumat. 

However, the plane stress version still has problem as below, and was hoping getting some help here. 

When using small fixed time step or use "element by element" option for increment size control, vumat matches ABQ results almost exactly. in this case, vumat is compared with ABQ parallel network model with only the elastic network active as it also assumes incompressible in plane stress.  However, if use automatic time increment control without "element by element" option, vumat goes into instable fairly quickly and almost all of sudden in a few increments after performing well for thousands of increments. I just could not figure out what the problem is in the vumat.    it seems that vumat drives ABQ to choose too large increment resulting in instability. 

files can be downloaded here!AvzF3H-Ljv02g8FsQCGROQuwyuleew?e=E2WWmP

elastic-D.cpp -- vumat  written in c. should work with abaqus 2019 - 2021, double precision

shell-shear-ABQ.inp  pure shear case, ABQ built-in material model

shell-shear-vumat.inp pure shear case, vumat

shell-uniaxial-ABQ.inp uniaxial tensile case, ABQ built-in material model

shell-uniaxial-vumat.inp uniaxial tensile case, vumat



2 Replies
Posts: 3998
Joined: 5 years ago

I have not looked at your files, but I think I have seen that in the past too. I would just make Abaqus pick a smaller dt in order to overcome the issue.


Posts: 2
Topic starter
New Member
Joined: 4 years ago

thanks Jorgen. wonder what might be the cause.  whether it was strainInc(3).  used  -1.*(strainInc(1)+strainInc(2)) vs. (F33-previous increment F33)/F33; results were the same, but stable dt were quite different.  if the incompressible F33 is strictly enforced, should some artificial compressibility be added to stress calculation?