vumat-hyperelastic-results mismatch with inbuilt abaqus model

[nepolean.r]

Dear Professor Jorgen,

I am Nepolean. Doing my Masters in Univ of Cincinnati. (Mech - CAE)
This is the first time I am writing to the group, but constantly go thru all the emails which are very useful to learn many new things.

My research area has a major part in material model development.
Recently I started to learn writing VUMAT .

To start with , I wrote VUMAT for hyperelastic mooney rivlin model (compressible) and tested on single C3D8R element using ABAQUS 6.9.1

VUMAT is based on U tensor - since manual says better to use U instead of F. I tried using F also. The results were almost same - either using F or U.

I compared the results of VUMAT (using U) with the results of inbuilt abaqus mooney rivlin model.

- for uniaxial tension and compression loads, results were closely matching.

- for shear stress load (element bottom constrained and top right nodes has load in x direction)
VUMAT underpredicts von mises stress by almost 20 times and underpredicts other directional stresses by 10 to 15 times

- for combined load cases (tension + shear , compression + shear)
Only the vonmises stress results are closely matching, all other individual direction stresses are out of match.

- combined loading in 2D elements, overpredicts stress (mises + all directional) by 25% to 40 % and displacement results are closely matching

I have also attached the .for file for your reference

possible mistake I thought
-corotational system : does the final stress vector needs to be multiplied by R tensor ? from where can we get the R tensor
but the manual says, only if F is used, then stress needs to be rotated back

kindly give me any possible feedback

thanks a lot for your time and hoping for your possible help.

Nepolean
Grad Student
Univ of Cincinnati

[nepolean.r]

c
subroutine vumat(nblock, ndi, nshr, nstatev, nfieldv, nprops,
. lanneal, stepTime, totTime, dt, cmname, coordMp,
. charLen, props, density, Dstrain, rSpinInc, temp0,
. U0, F0, field0, stressVec0, state0,
. intEne0, inelaEn0, temp1, U1,
. F1, field1, stressVec1, state1, intEne1, inelaEn1)
implicit none
integer nblock, ndi, nshr, nstatev, nfieldv, nprops, lanneal
real stepTime, totTime, dt
character*80 cmname
real coordMp(nblock,*)
integer charLen
real props(nprops), density(nblock),
. Dstrain(nblock,ndi+nshr), rSpinInc(nblock,nshr),
. temp0(nblock), U0(nblock,ndi+nshr),
. F0(nblock,ndi+nshr+nshr), field0(nblock,nfieldv),
. stressVec0(nblock,ndi+nshr), state0(nblock,nstatev),
. intEne0(nblock), inelaEn0(nblock), temp1(nblock),
. U1(nblock,ndi+nshr), F1(nblock,ndi+nshr+nshr),
. field1(nblock,nfieldv), stressVec1(nblock,ndi+nshr),
. state1(nblock,nstatev), intEne1(nblock), inelaEn1(nblock)

c local variables
real F(3,3), B(3,3), Bsqr(3,3),J, t1,C10,C01,D1
real p1,p2,p3,s1,s2,s3,s4,s5,I1bar
integer i

C10 = props(1)
C01 = props(2)
D1 = props(3)

c loop through all blocks
do i = 1, nblock
c setup F (upper diagonal part)
F(1,1) = U1(i,1)
F(2,2) = U1(i,2)
F(3,3) = U1(i,3)
F(1,2) = U1(i,4)
if (nshr .eq. 1) then
F(2,3) = 0.0
F(1,3) = 0.0
else
F(2,3) = U1(i,5)
F(1,3) = U1(i,6)
end if

c calculate J
p1 = F(1,1) * (F(2,2)*F(3,3) - F(2,3)**2)
p2 = F(1,2) * (F(2,3)*F(1,3) - F(1,2)*F(3,3))
p3 = F(1,3) * (F(1,2)*F(2,3) - F(2,2)*F(1,3))
J = p1 + p2 + p3
c from Abaqus Theory Manual : F = F * J^-1/3
c so for B = F*Ft, J^(-2/3) will appear since we are using F tensor directly

t1 = J**(-2.0/3.0)

c B = J^(-2/3) F Ft (upper diagonal part)
B(1,1) = t1*(F(1,1)**2 + F(1,2)**2 + F(1,3)**2)
B(2,2) = t1*(F(1,2)**2 + F(2,2)**2 + F(2,3)**2)
B(3,3) = t1*(F(1,3)**2 + F(2,3)**2 + F(3,3)**2)
B(1,2) = t1*(F(1,1)*F(1,2)+F(1,2)*F(2,2)+F(1,3)*F(2,3))
B(2,3) = t1*(F(1,2)*F(1,3)+F(2,2)*F(2,3)+F(2,3)*F(3,3))
B(1,3) = t1*(F(1,1)*F(1,3)+F(1,2)*F(2,3)+F(1,3)*F(3,3))

c Bsqr = BxB (upper diagonal part)
Bsqr(1,1) = B(1,1)**2 + B(1,2)**2 + B(1,3)**2
Bsqr(2,2) = B(1,2)**2 + B(2,2)**2 + B(2,3)**2
Bsqr(3,3) = B(1,3)**2 + B(2,3)**2 + B(3,3)**2
Bsqr(1,2) = B(1,1)*B(1,2)+B(1,2)*B(2,2)+B(1,3)*B(2,3)
Bsqr(2,3) = B(1,2)*B(1,3)+B(2,2)*B(2,3)+B(2,3)*B(3,3)
Bsqr(1,3) = B(1,1)*B(1,3)+B(1,2)*B(2,3)+B(1,3)*B(3,3)

c Stress = 2*C10/J * Dev(B) + 2/J*C01*I1bar*Dev(B)- 2/J*C01*Dev(Bsqr) + 2*(J-1)/D1 * Eye
s1 = (B(1,1) + B(2,2) + B(3,3)) / 3.0
s2 = 2*C10/J
s3 = 2/D1*(J-1)
s4 = 2/J*C01
I1bar=B(1,1)+B(2,2)+B(3,3)
s5=(Bsqr(1,1) + Bsqr(2,2) + Bsqr(3,3)) / 3.0
c
StressVec1(i,1) = s2 * (B(1,1) - s1) + s4*I1bar*(B(1,1) - s1)
1 -s4*(Bsqr(1,1)-s5) + s3
StressVec1(i,2) = s2 * (B(2,2) - s1) + s4*I1bar*(B(2,2) - s1)
1 -s4*(Bsqr(2,2)-s5) + s3
StressVec1(i,3) = s2 * (B(3,3) - s1) + s4*I1bar*(B(3,3) - s1)
1 -s4*(Bsqr(3,3)-s5) + s3
StressVec1(i,4) = s2 * (B(1,2)) + s4*I1bar*(B(1,2))
1 -s4*(Bsqr(1,2))
if (nshr .eq. 3) then
StressVec1(i,5) = s2 * (B(2,3)) + s4*I1bar*(B(2,3))
1 -s4*(Bsqr(2,3))
StressVec1(i,6) = s2 * (B(1,3)) + s4*I1bar*(B(1,3))
1 -s4*(Bsqr(1,3))
end if

end do
return
end

[Jorgen]

At first sight it looks like you are doing everything right. I recommend that you add a few print statements to your code in order to verify that the calculated F, B, Stress, etc are what you expect. There might be some small bug in there :confused:

-Jorgen

[nepolean.r]

Dear Professor,
thanks for the feedback.
I will do that and get back to you

Nepolean