自己编辑个JC本构，不考虑温度和损伤的，但是调试没有数值结果求大佬看看，指导一下
subroutine vumat (C Read only -
1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
2 stepTime, totalTime, dt, cmname, coordMp, charLength,
3 props, density, strainInc, relSpinInc,
4 tempOld, stretchOld, defgradOld, fieldOld,
5 stressOld, stateOld, enerInternOld, enerInelasOld,
6 tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
7 stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
dimension coordMp(nblock,*), charLength(nblock), props(nprops),
1 density(nblock), strainInc(nblock,ndir+nshr),
2 relSpinInc(nblock,nshr), tempOld(nblock),
3 stretchOld(nblock,ndir+nshr),
4 defgradOld(nblock,ndir+nshr+nshr),
5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(nblock),
8 stretchNew(nblock,ndir+nshr),
9 defgradNew(nblock,ndir+nshr+nshr),
1 fieldNew(nblock,nfieldv),
2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
3 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
parameter(zero=0.d0,one=1.d0,two=2.d0,three=3.d0,
* third=one/three,half=0.5d0,threeHalfs=1.5d0,rstrain=1e0)
C moxingcanshu
e=props(1)
xnu=props(2)
A=props(3)
B=props(4)
C=props(5)
n=props(6)
C Lchangshu
twomu=e/(one+xnu)
alamda=xnu*twomu/(one-two*xnu)
thremu=threeHalfs*twomu
C Trial stress
if(stepTime.eq.zero)then
do k=1,nblock
trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
stressNew(k,1)=stressOld(k,1)
1 +twomu*strainInc(k,1)+alamda*trace
stressNew(k,2)=stressOld(k,2)
1 +twomu*strainInc(k,2)+alamda*trace
stressNew(k,3)=stressOld(k,3)
1 +twomu*strainInc(k,3)+alamda*trace
stressNew(k,4)=stressOld(k,4)+twomu*strainInc(k,4)
if(nshr.gt.1)then
stressNew(k,5)=stressOld(k,5)+twomu*strainInc(k,5)
stressNew(k,6)=stressOld(k,6)+twomu*strainInc(k,6)
end if
end do
else
do k=1,nlock
st1=stateOld(k,1) !塑性应变
st2=stateOld(k,2) !应变率
C diaoyongjc
call xvhard(yieldOld,hard,st1,A,B,
1 C,n,st2)
C traceyingli
trace=strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
s11=stressOld(k,1)
1 +twomu*strainInc(k,1)+alamda*trace
s22=stressOld(k,2)
1 +twomu*strainInc(k,2)+alamda*trace
s33=stressOld(k,3)
1 +twomu*strainInc(k,3)+alamda*trace
s12=stressOld(k,4)+twomu*strainInc(k,4)
if ( nshr .gt. 1 ) then
s13=stressOld(k,5) + twomu * strainInc(k,5)
s23=stressOld(k,6) + twomu * strainInc(k,6)
end if
C Mieseyingli
smean=third*(s11+s22+s33)
s11=s11-smean
s22=s22-smean
s33=s33-smean
if ( nshr .eq. 1 ) then
vmises=sqrt(threeHalfs*
1 (s11*s11+s22*s22+s33*s33+two*s12*s12))
else
vmises = sqrt(op5*(s11*s11+s22*s22+s33*s33+
1 two*s12*s12+two*s13*s13+two*s23*s23))
end if
write(*,*) 'vmises',vmises
C suxingpanding
sigdif=vmises-yieldOld
facyld=zero
if(sigdif.gt.zero)then
facyld=one
end if
C dxsuxingyb
deqps=facyld*sigdif/(thremu+hard)
st2=deqps/dt
Cgengxinyingli
yieldNew=yieldOld+hard*deqps
factor=yieldNew/(yieldNew+thremu*deqps)
stressNew(k,1)=s11*factor+smean
stressNew(k,2)=s22*factor+smean
stressNew(k,3)=s33*factor+smean
stressNew(k,4)=s12*factor
if ( nshr .gt. 1 ) then
stressNew(k,5) = s13 * factor
stressNew(k,6) = s23 * factor
end if
Cgxzt
stateNew(k,1)=st1+deqps
stateNew(k,2)=st2
Cgxnengliang
if ( nshr .eq. 1 ) then
stressPower=half*(
1 (stressOld(k,1)+stressNew(k,1))*strainInc(k,1)+
1 (stressOld(k,2)+stressNew(k,2))*strainInc(k,2)+
1 (stressOld(k,3)+stressNew(k,3))*strainInc(k,3)+
1 (stressOld(k,4)+stressNew(k,4))*strainInc(k,4))
else
stressPower = half * (
* ( stressOld(k,1) + stressNew(k,1) ) * strainInc(k,1) +
* ( stressOld(k,2) + stressNew(k,2) ) * strainInc(k,2) +
* ( stressOld(k,3) + stressNew(k,3) ) * strainInc(k,3))+
* ( stressOld(k,4) + stressNew(k,4) ) * strainInc(k,4) +
* ( stressOld(k,5) + stressNew(k,5) ) * strainInc(k,5) +
* ( stressOld(k,6) + stressNew(k,6) ) * strainInc(k,6)
end if
enerInternNew(k)=enerInternOld(k)
1 +stressPower/density(k)
Cgengxfeitan
plasticWorkInc=half*(yieldOld+yieldNew)*deqps
enerInelasNew(k)=enerInelasOld(k)
1 +plasticWorkInc/density(k)
end do
end if
return
end
C
subroutine xvhard(syield,hard,eqpals,table3,table4,
1 table5,table6,sr)
include 'vaba_param.inc'
C
clnrstrain=log(sr/restrainr)
SRc=1.d0+table5*clnrstrain
syield=(table3+table4*eqpals**table6)*SRc
hard=table6*eqpals**(table6-1)*SRc
end subroutine
求大佬指点，大佬们求求了，大神们，我还是子程序的萌新，希望指点一二
页:
[1]