VUMAT状态变量
为什么我的子程序不能运行,求帮助!!!! VUMAT头文件
subroutine vumat(
! Read only (unmodifiable)variables -
1nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2stepTime, totalTime, dt, cmname, coordMp, charLength,
3props, density, strainInc, relSpinInc,
4tempOld, stretchOld, defgradOld, fieldOld,
5stressOld, stateOld, enerInternOld, enerInelasOld,
6tempNew, stretchNew, defgradNew, fieldNew,
! Write only (modifiable) variables -
7stressNew, stateNew, enerInternNew, enerInelasNew )
!
include 'vaba_param.inc'
dimension props(nprops), density(nblock), coordMp(nblock,*),
1charLength(nblock), strainInc(nblock,ndir+nshr),
2relSpinInc(nblock,nshr), tempOld(nblock),
3stretchOld(nblock,ndir+nshr),
4defgradOld(nblock,ndir+nshr+nshr),
5fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6stateOld(nblock,nstatev), enerInternOld(nblock),
7enerInelasOld(nblock), tempNew(nblock),
8stretchNew(nblock,ndir+nshr),
8defgradNew(nblock,ndir+nshr+nshr),
9fieldNew(nblock,nfieldv),
1stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
2enerInternNew(nblock), enerInelasNew(nblock)
!
character*80 cmname
parameter(zero=0.d0,half=0.5,two=2.d0,one=1.d0)
dimension CFULL(6,6)
dimension ss(6)
! 输入材料强度参数,用于失效判断
xt=props(1)
xc=props(2)
yt=props(3)
yc=props(4)
zt=props(5)
zc=props(6)
s12=props(7)
s23=props(8)
s13=props(9)
! 输入材料的弹性常数(下标顺序均按照沈观林的教材规定)
! props(13)=emu21,为主泊松比
! props(14)=emu32,props(15)=emu31为次泊松比
e1=props(10)
e2=props(11)
e3=props(12)
emu12=props(13)
emu23=props(14)
emu13=props(15)
g12=props(16)
g23=props(17)
g13=props(18)
! ******Standard strain****
sn1t=xt/e1
sn1c=xc/e1
sn2t=yt/e2
sn2c=yc/e2
sn3t=zt/e3
sn12=s12/g12
sn23=s23/g23
sn13=s13/g13
E11=e1
E22=e2
E33=e3
GG12=g12
GG23=g23
GG13=g13
! 在初始步给状态参数赋初值,状态参数为0说明材料完好
NTENS=NDIR+NSHR
do 100 i = 1,nblock
if(totalTime .eq.0.02)then
do j=1,6
stateOld(i,j)=0.0
strainInc(i,j)=0.0
end do
endif
do j=1,6
stateNew(i,j)=stateOld(i,j)+strainInc(i,j)
ss(j)=stateOld(i,j)+strainInc(i,j)
end do
damL=0.0
damT=0.0
damZ=0.0
damTZ=0.0
damLZ=0.0
damLT=0.0
term1=(1/sn1t-1/sn1c)*ss(1)+1/sn1t/sn1c*ss(1)*ss(1)
term2=(1/sn2t-1/sn2c)*ss(2)+1/sn2t/sn2c*ss(2)*ss(2)
term3=(1/sn3t-1/sn3c)*ss(3)+1/sn3t/sn3c*ss(3)*ss(3)
term4=(2*ss(5)/sn23)**2
term5=(2*ss(6)/sn13)**2
term6=(2*ss(4)/sn12)**2
term7=1/(sqrt(sn1t*sn1c*sn2t*sn2c))*ss(1)*ss(2)
* +1/(sqrt(sn2t*sn2c*sn3t*sn3c))*ss(2)*ss(3)
* +1/(sqrt(sn3t*sn3c*sn1t*sn1c))*ss(1)*ss(3)
TWD=term1+term2+term3+term4+term5+term6-term7
if(TWD .gt.one) then
if(term1 .ge. term2) then
term8=term1
damL=1.0
damT=0.0
damZ=0.0
damTZ=0.0
damLZ=0.0
damLT=0.0
else
term8=term2
damL=0.0
damT=1.0
damZ=0.0
damTZ=0.0
damLZ=0.0
damLT=0.0
end if
if(term8 .ge. term3) then
else
term8=term3
damL=0.0
damT=0.0
damZ=1.0
damTZ=0.0
damLZ=0.0
damLT=0.0
end if
if(term8 .ge. term4) then
else
term8=term4
damL=0.0
damT=0.0
damZ=0.0
damTZ=1.0
damLZ=0.0
damLT=0.0
end if
if(term8 .GE. term5) then
else
term8=term5
damL=0.0
damT=0.0
damZ=0.0
damTZ=0.0
damLZ=1.0
damLT=0.0
end if
if(term8 .GE. term6) then
else
term8=term6
damL=0.0
damT=0.0
damZ=0.0
damTZ=0.0
damLZ=0.0
damLT=1.0
end if
else
end if
stateNew(i,7)=max(stateOld(i,7),damL)
stateNew(i,8)=max(stateOld(i,8),damT)
stateNew(i,9)=max(stateOld(i,9),damZ)
stateNew(i,10)=max(stateOld(i,10),damLT)
stateNew(i,11)=max(stateOld(i,11),damTZ)
stateNew(i,12)=max(stateOld(i,12),damLZ)
e11=e1
e22=e2
e33=e3
xemu12=emu12
xemu23=emu23
xemu13= emu13
gg12=g12
gg23=g23
gg13=g13
if(stateNew(i,7).EQ.1.0)then
e11=0.01*e1
gg13=0.2*g13
gg12=0.2*g12
end if
if(stateNew(i,8).EQ.1.0)then
e22=0.01*e2
gg23=0.2*g23
gg12=0.2*g12
end if
if(stateNew(i,9).EQ.1.0)then
e33=0.01*e3
gg23=0.2*g23
gg13=0.2*g13
end if
if(stateNew(i,10).EQ.1.0)then
e22=0.01*e2
e33=0.01*e3
gg23=0.2*g23
gg12=0.2*g12
gg13=0.2*g13
end if
if(stateNew(i,11).EQ.1.0)then
e11=0.8*e1
e33=0.01*e3
gg23=0.2*g23
gg13=0.2*g13
end if
if(stateNew(i,12).EQ.1.0)then
e11=0.8*e1
e22=0.01*e2
gg23=0.2*g23
gg12=0.2*g12
end if
!Calculate stress
xemu21 = xemu12 * e22 / e11
xemu31 = xemu13 * e33 / e11
xemu32 = xemu23 * e33 / e22
! Initial analysis,assumely pure elastic behaviour
gg = one / ( one - xemu12*xemu21 - xemu23*xemu32 - xemu31*xemu13
* - two*xemu21*xemu32*xemu13 )
CFULL(1,1) = e11 * ( one - xemu23*xemu32 ) * gg
CFULL(2,2) = e22 * ( one - xemu13*xemu31 ) * gg
CFULL(3,3) = e33 * ( one - xemu12*xemu21 ) * gg
CFULL(1,2) = e11 * ( xemu21 + xemu31*xemu23 ) * gg
CFULL(1,3) = e11 * ( xemu31 + xemu21*xemu32 ) * gg
CFULL(2,3) = e22 * ( xemu32 + xemu12*xemu31 ) * gg
CFULL(4,4) = gg12
CFULL(5,5) = gg23
CFULL(6,6) = gg13
CFULL(2,1) = CFULL(1,2)
CFULL(3,2) = CFULL(2,3)
CFULL(3,1) = CFULL(1,3)
stressNew(i,1)=CFULL(1,1)*stateNew(i,1)+
& CFULL(1,2)*stateNew(i,2)+
& CFULL(1,3)*stateNew(i,3)
stressNew(i,2)=CFULL(2,1)*stateNew(i,1)+
& CFULL(2,2)*stateNew(i,2)+
& CFULL(2,3)*stateNew(i,3)
stressNew(i,3)=CFULL(3,1)*stateNew(i,1)+
& CFULL(3,2)*stateNew(i,2)+
& CFULL(3,3)*stateNew(i,3)
stressNew(i,4)=CFULL(4,4)*2*stateNew(i,4)
stressNew(i,5)=CFULL(5,5)*2*stateNew(i,5)
stressNew(i,6)=CFULL(6,6)*2*stateNew(i,6)
100 continue
return
end
!
为什么删去标红的那几行就能运行了,求帮助 楼主在设置ss(j)的时候有错误,交流qq569250598 你的问题搞定没?我也是做类似的,共同交流,qq569250598 顶一下,楼主问题解决没? 楼主问题解决没?共同探讨,QQ569250598 跟着学习 学习一下
学习一下
页:
[1]