Good mornig everybody ,
I wrote a Umat routine for the Yeoh hyperelastic model,
In order to compare the built in yeoh model and my Umat subroutine , I apply a uniaxial displacement on a single element and compare the stress strain rates
the problem is that I obtain different curves for the strees strain rates , and a different deformation I can't figure out why ,
My question is, do you have any documents or links describing the theory for the Yeoh model, and could you please have a look at my Umat code, attached here.
Thank you in advance for your response.
Best regards
Hamza
c
subroutine umat(stress,statev,ddsdde,sse,spd,scd,rpl,ddsddt,
1 drplde,drpldt,stran,dstran,time,dtime,temp,
2 dtemp,predef,dpred,cmname,ndi,nshr,ntens,nstatv,pr ops,
3 nprops,coords,drot,pnewdt,celent,
4 dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)
include 'aba_param.inc'
character*80 cmname
dimension stress(ntens),statev(nstatv),ddsdde(ntens,ntens),
1 ddsddt(ntens),drplde(ntens),stran(ntens),dstran(nt ens),
2 predef(1),dpred(1),props(nprops),coords(3),drot(3, 3),
3 time(2),dfgrd0(3,3),dfgrd1(3,3)
dimension bbar(6)
dimension y(6,6),bb(6,6),cs(6,6)
dimension distgr(3,3)
dimension xh(6,6)
parameter(zero=0.d0, one=1.d0, two=2.d0, three=3.d0, six=6.d0)
parameter(half = 0.5d0,third = 1.d0/3.d0,twothrd=2.d0/3.d0)
parameter(four=4.d0)
c10=props(1)
c20=props(2)
c30=props(3)
c
det = dfgrd1(1,1) * dfgrd1(2,2) * dfgrd1(3,3)
1 - dfgrd1(1,2) * dfgrd1(2,1) * dfgrd1(3,3)
c
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 = abs(det)**(-one /three)
c
do k1 = 1, 3
do k2 = 1, 3
distgr(k2,k1) = scale * dfgrd1(k2,k1)
end do
end do
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
trbbar = (bbar(1) + bbar(2) + bbar(3)) / three
xinv = bbar(1) + bbar(2) + bbar(3)
der1=c10+two*c20*(xinv-three)+three*c30*(xinv-three)**2
der2=two*c20+six*c30*(xinv-three)
c
eg=two*der1/det
c
do k1=1,ndi
stress(k1)=eg*(bbar(k1)-trbbar)+pr
end do
do k1=ndi+1,ntens
stress(k1)=eg*bbar(k1)
end do
c
a1=two/det*der1
a2=-four/(three*det)*(der1+xinv*der2)
a3=four/det*der2
xh(1,1)=two*bbar(1)*a1
xh(1,2)=zero
xh(1,3)=zero
xh(1,4)=bbar(4)*a1
xh(1,5)=bbar(5)*a1
xh(1,6)=zero
xh(2,2)=two*bbar(2)*a1
xh(2,3)=zero
xh(2,4)=bbar(4)*a1
xh(2,5)=zero
xh(2,6)=bbar(6)*a1
xh(3,3)=two*bbar(3)*a1
xh(3,4)=zero
xh(3,5)=bbar(5)*a1
xh(3,6)=bbar(6)*a1
xh(4,4)=(bbar(1)+bbar(2))/two*a1
xh(4,5)=bbar(6)/two*a1
xh(4,6)=bbar(5)/two*a1
xh(5,5)=(bbar(1)+bbar(3))/two*a1
xh(5,6)=bbar(4)/two*a1
xh(6,6)=(bbar(2)+bbar(3))/two*a1
y(1,1)=two*bbar(1)*a2
y(1,2)=(bbar(1)+bbar(2))*a2
y(1,3)=(bbar(1)+bbar(3))*a2
y(1,4)=bbar(4)*a2
y(1,5)=bbar(5)*a2
y(1,6)=bbar(6)*a2
y(2,2)=two*bbar(2)*a2
y(2,3)=(bbar(2)+bbar(3))*a2
y(2,4)=bbar(4)*a2
y(2,5)=bbar(5)*a2
y(2,6)=bbar(6)*a2
y(3,3)=two*bbar(3)*a2
y(3,4)=bbar(4)*a2
y(3,5)=bbar(5)*a2
y(3,6)=bbar(6)*a2
y(4,4)=zero
y(4,5)=zero
y(4,6)=zero
y(5,5)=zero
y(5,6)=zero
y(6,6)=zero
bb(1,1)=bbar(1)**2*a3
bb(1,2)=bbar(1)*bbar(2)*a3
bb(1,3)=bbar(1)*bbar(3)*a3
bb(1,4)=bbar(1)*bbar(4)*a3
bb(1,5)=bbar(1)*bbar(5)*a3
bb(1,6)=bbar(1)*bbar(6)*a3
bb(2,2)=bbar(2)**2*a3
bb(2,3)=bbar(2)*bbar(3)*a3
bb(2,4)=bbar(2)*bbar(4)*a3
bb(2,5)=bbar(2)*bbar(5)*a3
bb(2,6)=bbar(2)*bbar(6)*a3
bb(3,3)=bbar(3)**2*a3
bb(3,4)=bbar(3)*bbar(4)*a3
bb(3,5)=bbar(3)*bbar(5)*a3
bb(3,6)=bbar(3)*bbar(6)*a3
bb(4,4)=bbar(4)**2*a3
bb(4,5)=bbar(4)*bbar(5)*a3
bb(4,6)=bbar(4)*bbar(6)*a3
bb(5,5)=bbar(5)**2*a3
bb(5,6)=bbar(5)*bbar(6)*a3
bb(6,6)=bbar(6)**2*a3
cs(1,1)=xh(1,1)+y(1,1)+bb(1,1)
cs(1,2)=xh(1,2)+y(1,2)+bb(1,2)
cs(1,3)=xh(1,3)+y(1,3)+bb(1,3)
cs(1,4)=xh(1,4)+y(1,4)+bb(1,4)
cs(1,5)=xh(1,5)+y(1,5)+bb(1,5)
cs(1,6)=xh(1,6)+y(1,6)+bb(1,6)
cs(2,2)=xh(2,2)+y(2,2)+bb(2,2)
cs(2,3)=xh(2,3)+y(2,3)+bb(2,3)
cs(2,4)=xh(2,4)+y(2,4)+bb(2,4)
cs(2,5)=xh(2,5)+y(2,5)+bb(2,5)
cs(2,6)=xh(2,6)+y(2,6)+bb(2,6)
cs(3,3)=xh(3,3)+y(3,3)+bb(3,3)
cs(3,4)=xh(3,4)+y(3,4)+bb(3,4)
cs(3,5)=xh(3,5)+y(3,5)+bb(3,5)
cs(3,6)=xh(3,6)+y(3,6)+bb(3,6)
cs(4,4)=xh(4,4)+y(4,4)+bb(4,4)
cs(4,5)=xh(4,5)+y(4,5)+bb(4,5)
cs(4,6)=xh(4,6)+y(4,6)+bb(4,6)
cs(5,5)=xh(5,5)+y(5,5)+bb(5,5)
cs(5,6)=xh(5,6)+y(5,6)+bb(5,6)
cs(6,6)=xh(6,6)+y(6,6)+bb(6,6)
do k1 = 1, ntens
do k2 = 1, k1 - 1
cs(k1, k2) = cs(k2, k1)
end do
end do
ddsdde(1,1) = cs(1,1) - one/three*(cs(1,1)+cs(1,2)+cs(1,3))
ddsdde(1,2) = cs(1,2) - one/three*(cs(1,1)+cs(1,2)+cs(1,3))
ddsdde(1,3) = cs(1,3) - one/three*(cs(1,1)+cs(1,2)+cs(1,3))
ddsdde(1,4) = cs(1,4)
ddsdde(1,5) = cs(1,5)
ddsdde(1,6) = cs(1,6)
ddsdde(2,2) = cs(2,2) - one/three*(cs(2,1)+cs(2,2)+cs(2,3))
ddsdde(2,3) = cs(2,3) - one/three*(cs(2,1)+cs(2,2)+cs(2,3))
ddsdde(2,4) = cs(2,4)
ddsdde(2,5) = cs(2,5)
ddsdde(2,6) = cs(2,6)
ddsdde(3,3) = cs(3,3) - one/three*(cs(3,1)+cs(3,2)+cs(3,3))
ddsdde(3,4) = cs(3,4)
ddsdde(3,5) = cs(3,5)
ddsdde(3,6) = cs(3,6)
ddsdde(4,4) = cs(4,4)
ddsdde(4,5) = cs(4,5)
ddsdde(4,6) = cs(4,6)
ddsdde(5,5) = cs(5,5)
ddsdde(5,6) = cs(5,6)
ddsdde(6,6) = cs(6,6)
do k1 = 1, ntens
do k2 = 1, k1 - 1
ddsdde(k1, k2) = ddsdde(k2, k1)
end do
end do
return
end
| Attachment | Size |
|---|---|
| stress-strain-rates.JPG | 91.8 KB |
