Hello,
I am new with VUMAT. Recently I have finished reading VUMAT example for Neo-Hookean material from the Mechanics of Solid Polymers textbook. Could anyone share more VUMAT examples of Hyperelastic materials for learning purposes? Thank you very much.
Regards,
Willsen
Hello,
I have written a hyperelastic VUMAT that can handle strain-energy functions of the form:
W = K1(I1-3)^m + K2(I1-3)^p + K3(I1-3)^q
which covers the Neo-Hookean, Yeoh, and (almost) Davies-De-Thomas (1994) models. I think it will be published in an appendix of a related work soon. If you want to take a look before then, let me know and I may be allowed to share it with you before it hits the press.
Regards,
Travis
Hi Somaye,
Recalling my experience on how to learn VUMAT, I did a lot of in depth study on the basic of Continuum Mechanics. Through a proper understanding of deformation gradient, different strain and stress measures in large deformation, it becomes really straight forward to understand the different physical quantities available in VUMAT. Jorgens VUMAT for Neo Hookean material was super valuable to me.
My suggestion will be to understand his VUMAT in depth and apply it in a single element test. Do an identical test using ABAQUS in-built Neo-Hookean material model and compare the results.
Build upon that. Try to write a more complicated VUMAT to describe Mooney-Rivlin material and compare it with the equivalent mateiral model from ABAQUS, and so on.
Attached is an example of VUMAT for Mooney Rivlin model:
! This is a VUMAT subroutine for Mooney-Rivlin
subroutine vumat(nblock, ndi, nshr, nstatev, nfieldv, nprops, &
lanneal, stepTime, totTime, dt, cmname, coordMp, &
charLen, props, density, strainInc, relSpinInc, temp0, &
U0, F0, field0, stress0, state0, &
intEne0, inelaEn0, temp1, U1, &
F1, field1, stress1, state1, intEne1, inelaEn1)
implicit none
integer :: nblock, ndi, nshr, nstatev, nfieldv, nprops, lanneal
real :: stepTime, totTime, dt
character*80 :: cmname
real :: coordMp(nblock,*)
real :: charLen, props(nprops), density(nblock), &
strainInc(nblock,ndi+nshr), relSpinInc(nblock,nshr), &
temp0(nblock), U0(nblock,ndi+nshr), &
F0(nblock,ndi+nshr+nshr), field0(nblock,nfieldv), &
stress0(nblock,ndi+nshr), state0(nblock,nstatev), &
intEne0(nblock), inelaEn0(nblock), temp1(nblock), &
U1(nblock,ndi+nshr), F1(nblock,ndi+nshr+nshr), &
field1(nblock,nfieldv), stress1(nblock,ndi+nshr), &
state1(nblock,nstatev), intEne1(nblock), inelaEn1(nblock)
! local variables
real :: J, U(3,3), U2dev(3,3), U4dev(3,3), C10, C01, D1 ,traceU2dev, traceU4dev , constant1, constant2
integer :: i
C10 = props(1)
C01 = props(2)
D1 = props(3)
! loop through all blocks
do i = 1, nblock
! setup F (upper diagonal part)
U(1,1) = U1(i,1)
U(2,2) = U1(i,2)
U(3,3) = U1(i,3)
U(1,2) = U1(i,4)
U(2,1) = U(1,2)
if (nshr .eq. 1) then
U(2,3) = 0.0
U(1,3) = 0.0
else
U(2,3) = U1(i,5)
U(1,3) = U1(i,6)
end if
U(3,2)=U(2,3)
U(3,1)=U(1,3)
! calculate J
J = (U(1,1)*(U(2,2)*U(3,3)-U(2,3)**2)) + (U(1,2)*(U(2,3)*U(1,3)-U(1,2)*U(3,3))) + (U(1,3)*(U(1,2)*U(2,3)-U(2,2)*U(1,3)))
! Define the square of the deviatoric stretch tensor
U2dev=(J**(-2.0/3.0))*matmul(U,U)
! Define the forth power of the deviatoric stretch tensor
U4dev=matmul(U2dev,U2dev)
! Define some constants to simplify the stress expression
traceU2dev=(U2dev(1,1)+U2dev(2,2)+U2dev(3,3))
traceU4dev=(U4dev(1,1)+U4dev(2,2)+U4dev(3,3))
constant1=2.0*C10/J
constant2=2.0*C01/J
!The corotational stress
stress1(i,1) = (constant1+constant2*traceU2dev)*(U2dev(1,1)-1.0/3.0*traceU2dev) - constant2*(U4dev(1,1)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
stress1(i,2) = (constant1+constant2*traceU2dev)*(U2dev(2,2)-1.0/3.0*traceU2dev) - constant2*(U4dev(2,2)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
stress1(i,3) = (constant1+constant2*traceU2dev)*(U2dev(3,3)-1.0/3.0*traceU2dev) - constant2*(U4dev(3,3)-1.0/3.0*traceU4dev) + 2.0/D1*(J-1.0)
stress1(i,4) = (constant1+constant2*traceU2dev)*U2dev(1,2) - constant2*U4dev(1,2)
if (nshr .eq. 3) then
stress1(i,5) = (constant1+constant2*traceU2dev)*U2dev(2,3) - constant2*U4dev(2,3)
stress1(i,6) = (constant1+constant2*traceU2dev)*U2dev(3,1) - constant2*U4dev(3,1)
end if
end do
end subroutine vumat
Thank you so much Willsen!
Can the Jorgens VUMAT for Neo Hookean material or yours be run as .f or .for file in ABAQUS? It seems those formats are kind of different from the umat file I have :
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS, 3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC) C INCLUDE ABA_PARAM.INC C CHARACTER*8 MATERL DIMENSION STRESS(NTENS),STATEV(NSTATV), 1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS), 2 STRAN(NTENS),DSTRAN(NTENS),DFGRD0(3,3),DFGRD1(3,3), 3 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3) C C LOCAL ARRAYS C ---------------------------------------------------------------- C BBAR - DEVIATORIC RIGHT CAUCHY-GREEN TENSOR C DISTGR - DEVIATORIC DEFORMATION GRADIENT (DISTORTION TENSOR) C ---------------------------------------------------------------- C DIMENSION BBAR(6),DISTGR(3,3) C PARAMETER(ZERO=0.D0, ONE=1.D0, TWO=2.D0, THREE=3.D0, FOUR=4.D0) C C ---------------------------------------------------------------- C UMAT FOR COMPRESSIBLE NEO-HOOKEAN HYPERELASTICITY C CANNOT BE USED FOR PLANE STRESS C ---------------------------------------------------------------- C PROPS(1) - C10 C PROPS(2) - C01 C PROPS(3) - D1 C ---------------------------------------------------------------- C C ELASTIC PROPERTIES C C10=PROPS(1) C01=PROPS(2) D1 =PROPS(3) C C JACOBIAN AND DISTORTION TENSOR C DET=DFGRD1(1, 1)*DFGRD1(2, 2)*DFGRD1(3, 3) 1 -DFGRD1(1, 2)*DFGRD1(2, 1)*DFGRD1(3, 3) IF(NSHR.EQ.3) THEN DET=DET+DFGRD1(1, 2)*DFGRD1(2, 3)*DFGRD1(3, 1) 1 +DFGRD1(1, 3)*DFGRD1(3, 2)*DFGRD1(2, 1) 2 -DFGRD1(1, 3)*DFGRD1(3,1)*DFGRD1(2, 2) 3 -DFGRD1(2, 3)*DFGRD1(3, 2)*DFGRD1(1, 1) END IF SCALE=DET**(-ONE/THREE) DO K1=1, 3 DO K2=1, 3 DISTGR(K2, K1)=SCALE*DFGRD1(K2, K1) END DO END DO C C CALCULATE LEFT CAUCHY-GREEN TENSOR C BBAR(1)=DISTGR(1, 1)**2+DISTGR(1, 2)**2+DISTGR(1, 3)**2 BBAR(2)=DISTGR(2, 1)**2+DISTGR(2, 2)**2+DISTGR(2, 3)**2 BBAR(3)=DISTGR(3, 3)**2+DISTGR(3, 1)**2+DISTGR(3, 2)**2 BBAR(4)=DISTGR(1, 1)*DISTGR(2, 1)+DISTGR(1, 2)*DISTGR(2, 2) 1 +DISTGR(1, 3)*DISTGR(2, 3) IF(NSHR.EQ.3) THEN BBAR(5)=DISTGR(1, 1)*DISTGR(3, 1)+DISTGR(1, 2)*DISTGR(3, 2) 1 +DISTGR(1, 3)*DISTGR(3, 3) BBAR(6)=DISTGR(2, 1)*DISTGR(3, 1)+DISTGR(2, 2)*DISTGR(3, 2) 1 +DISTGR(2, 3)*DISTGR(3, 3) END IF C C CALCULATE THE STRESS C TRBBAR=(BBAR(1)+BBAR(2)+BBAR(3))/THREE EG=TWO*C10/DET EK=TWO/D1*(TWO*DET-ONE) PR=TWO/D1*(DET-ONE) DO K1=1,NDI STRESS(K1)=EG*(BBAR(K1)-TRBBAR)+PR END DO DO K1=NDI+1,NDI+NSHR STRESS(K1)=EG*BBAR(K1) END DO C C CALCULATE THE STIFFNESS C EG23=EG*TWO/THREE DDSDDE(1, 1)= EG23*(BBAR(1)+TRBBAR)+EK DDSDDE(2, 2)= EG23*(BBAR(2)+TRBBAR)+EK DDSDDE(3, 3)= EG23*(BBAR(3)+TRBBAR)+EK DDSDDE(1, 2)=-EG23*(BBAR(1)+BBAR(2)-TRBBAR)+EK DDSDDE(1, 3)=-EG23*(BBAR(1)+BBAR(3)-TRBBAR)+EK DDSDDE(2, 3)=-EG23*(BBAR(2)+BBAR(3)-TRBBAR)+EK DDSDDE(1, 4)= EG23*BBAR(4)/TWO DDSDDE(2, 4)= EG23*BBAR(4)/TWO DDSDDE(3, 4)=-EG23*BBAR(4) DDSDDE(4, 4)= EG*(BBAR(1)+BBAR(2))/TWO IF(NSHR.EQ.3) THEN DDSDDE(1, 5)= EG23*BBAR(5)/TWO DDSDDE(2, 5)=-EG23*BBAR(5) DDSDDE(3, 5)= EG23*BBAR(5)/TWO DDSDDE(1, 6)=-EG23*BBAR(6) DDSDDE(2, 6)= EG23*BBAR(6)/TWO DDSDDE(3, 6)= EG23*BBAR(6)/TWO DDSDDE(5, 5)= EG*(BBAR(1)+BBAR(3))/TWO DDSDDE(6, 6)= EG*(BBAR(2)+BBAR(3))/TWO DDSDDE(4,5)= EG*BBAR(6)/TWO DDSDDE(4,6)= EG*BBAR(5)/TWO DDSDDE(5,6)= EG*BBAR(4)/TWO END IF DO K1=1, NTENS DO K2=1, K1-1 DDSDDE(K1, K2)=DDSDDE(K2, K1) END DO END DO C RETURN END