- 积分
- 0
- 注册时间
- 2013-4-2
- 仿真币
-
- 最后登录
- 1970-1-1
|
RT,我是参照庄茁老师的umat,然后改编的vumat程序,但是通不过,不知道是为何,大侠能否出现,帮忙看看啊
subroutine vumat(
C Read only (unmodifiable)variables -
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 (modifiable) variables -
7 stressNew, stateNew, enerInternNew, enerInelasNew )
C
include 'vaba_param.inc'
C
dimension props(nprops), density(nblock), coordMp(nblock,*),
1 charLength(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),
8 defgradNew(nblock,ndir+nshr+nshr),
9 fieldNew(nblock,nfieldv),
1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
2 enerInternNew(nblock), enerInelasNew(nblock)
C
character*80 cmname
C
dimension eelas(6),eplas(6),flow(6),strain(6)
parameter (one=1.0D0,two=2.0D0,three=3.0D0,six=6.0D0, half =0.5d0,
1 third=1.d0/3.d0)
data newton,toler/40,1.D-6/
con=sqrt(2.d0/3.d0)
C props(1) - yang's Modulus
C props(2) - possion radio
C props(3) - inelastic heat fraction
C parameters of johnson cook model
C props(4) - A
C props(5) - B
C props(6) - n
C props(7) - C
C props(8) - m
C
C state(*,1)=elastic strain component 11
C state(*,2)=elastic strain component 22
C state(*,3)=elastic strain component 33
C state(*,4)=elastic strain component 12
C state(*,5)=elastic strain component 23
C state(*,6)=elastic strain component 31
C
C state(*,7)=plastic strain component 11
C state(*,8)=plastic strain component 22
C state(*,9)=plastic strain component 33
C state(*,10)=plastic strain component 12
C state(*,11)=plastic strain component 23
C state(*,12)=plastic strain component 31
C
C state(*,13)=equivalent plastic strain
C
C state(*,14)=total strain component 11
C state(*,15)=total strain component 22
C state(*,16)=total strain component 33
C state(*,17)=total strain component 12
C state(*,18)=total strain component 23
C state(*,19)=total strain component 31
C
C state(*,20)=equivalent strain
C
C state(*,21)=total stress component 11
C state(*,22)=total stress component 22
C state(*,23)=total stress component 33
C state(*,24)=total stress component 12
C state(*,25)=total stress component 23
C state(*,26)=total stress component 31
C
C state(*,27)==smises stress
emod=props(1)
enu=props(2)
eg=emod/(one+enu)/two
eg2=eg*two
eg6=eg*six
elam=eg2*(emod-eg2)/(eg6-two*emod)
! if ( stepTime .eq. zero ) then
! do i = 1, nblock
! trace = strainInc(i,1) + strainInc(i,2) + strainInc(i,3)
! stressNew(i,1) = stressOld(i,1)
! * + eg2 * strainInc(i,1) + elam * trace
! stressNew(i,2) = stressOld(i,2)
! * + eg2 * strainInc(i,2) + elam * trace
! stressNew(i,3) = stressOld(i,3)
! * + eg2 * strainInc(i,3) + elam * trace
! stressNew(i,4)=stressOld(i,4) + eg2 * strainInc(i,4)
! stressNew(i,5)=stressOld(i,5) + eg2 * strainInc(i,5)
! stressNew(i,6)=stressOld(i,6) + eg2 * strainInc(i,6)
!
! end do
! else
do i=1,nblock
trace=strainInc(i,1)+strainInc(i,2)+strainInc(i,3)
sig1=stressOld(i,1)+elam*trace+eg2*strainInc(i,1)
sig2=stressOld(i,2)+elam*trace+eg2*strainInc(i,2)
sig3=stressOld(i,3)+elam*trace+eg2*strainInc(i,3)
sig4=stressOld(i,4) +eg2*strainInc(i,4)
sig5=stressOld(i,5) +eg2*strainInc(i,5)
sig6=stressOld(i,6) +eg2*strainInc(i,6)
smean=third*(sig1+sig2+sig3)
sig11=sig1-smean
sig22=sig2-smean
sig33=sig3-smean
smisesold=(one/con)*sqrt(sig11**2+sig22**2+sig33**2+two*sig4**2
1 +two*sig5**2+two*sig6**2)
do k1=1,6
eelas(k1)=stateOld(i,k1)+strainInc(i,k1)
eplas(k1)=stateOld(i,k1+6)
stateNew(i,k1+13)=stateOld(i,k1+13)+strainInc(i,k1)
enddo
eqplas=stateNew(i,13)
str=(stateNew(i,14)+stateNew(i,15)+stateNew(i,16))/3
str1=stateNew(i,14)-str
str2=stateNew(i,15)-str
str3=stateNew(i,16)-str
str4=stateNew(i,17)
str5=stateNew(i,18)
str6=stateNew(i,19)
c equstrain:equivalent strain
equstrain=con*(str1**2+str2**2+str3**2+two*str4**2+
1 two*str5**2+two*str6**2)
stateNew(i,20)=equstrain
C
C calculate miser stress
C
C call userhard subroutine, get harding rate and yield stress
call userhard(syiel0,hard,eqplas,props(4))
C determine if yielding is actived
if (smises > (1+toler)*syiel0) then
C if plastic deformation has occured, determine the flow direction
shydr0=(stressNew(i,1)+stressNew(i,2)+stressNew(i,3))/3.d0
onesy=one/smises
do k1=1,3
flow(k1)=onesy*(stressNew(i,k1)-shydr0)
enddo
do k1=4,6
flow(k1)=onesy*stressNew(i,k1)
enddo
C read the johnson cook parameters
A=props(4)
B=props(5)
EN=props(6)
C=props(7)
EM=props(8)
C begin newton iteration
syield = syiel0
deqpl=(smises-syield)/eg3
dstres=toler*syiel0/eg3
deqmin=half*dt*exp(1.d0-4/c)
do knewton=1,newton
deqpl=max(deqpl,deqmin)
call userhard(syield,hard,eqplas+deqpl,props(4))
tvp=c*log(deqpl/dt)
tvp1=tvp+one
hard1=hard*tvp1+syield*c/deqpl
syield=syield*tvp1
rhs=smises-eg3*deqpl-syield
deqpl=deqpl+rhs/(eg3+hard1)
if(abs(rhs/eg3) < dstres) goto 140
enddo
140 continue
C update stress and strain
C whether the following codes is true?divited by 2?
do k1=1,3
stressNew(i,k1) = flow(k1)*syield+shydr0
eelas(k1) = eelas(k1)-three * flow(k1)*deqpl/two
eplas(k1) = eplas(k1)+three * flow(k1)*deqpl/two
enddo
do k1=4,6
stressNew(i,k1) = flow(k1)*syield
eelas(k1) = eelas(k1)-three * flow(k1)*deqpl
eplas(k1) = eplas(k1)+three * flow(k1)*deqpl
enddo
eqplas=eqplas+deqpl
spd=edqpl*(syiel0+syield)/two
rpl=props(3)*spd/dt
endif
smises=(stressNew(i,1)-stressNew(i,2)) *
1 (stressNew(i,1)-stressNew(i,2)) +
2 (stressNew(i,2)-stressNew(i,3)) *
3 (stressNew(i,2)-stressNew(i,3)) +
4 (stressNew(i,3)-stressNew(i,1)) *
5 (stressNew(i,3)-stressNew(i,1))
do k1=4,6
smises=smises+six*stressNew(i,k1)*stressNew(i,k1)
enddo
smises=sqrt(smises/two)
do k1=1,6
stateNew(i,k1) = eelas(k1)
stateNew(i,k1+6) = eplas(k1)
enddo
do k1=1,6
stateNew(i,20+k1) = stressNew(i,k1)
enddo
stateNew(i,13)=eqplas
stateNew(i,27)=smises
stressPower = half * (
* ( stressOld(i,1) + stressNew(i,1) ) * strainInc(i,1) +
* ( stressOld(i,2) + stressNew(i,2) ) * strainInc(i,2) +
* ( stressOld(i,3) + stressNew(i,3) ) * strainInc(i,3) ) +
* ( stressOld(i,4) + stressNew(i,4) ) * strainInc(i,4) +
* ( stressOld(i,5) + stressNew(i,5) ) * strainInc(i,5) +
* ( stressOld(i,6) + stressNew(i,6) ) * strainInc(i,6)
enerInternNew(i) = enerInternOld(i) + stressPower / density(i)
C
C Update the dissipated inelastic specific energy -
C
plasticWorkInc = half * ( smises + smisesold ) * deqps
enerInelasNew(i) = enerInelasOld(i)
* + plasticWorkInc / density(i)
enddo
! endif
return
end
C
C
C
subroutine userhard(syield,hard,eqplas,table)
include 'vaba_param.inc'
dimension table(3)
A=table(1)
B=table(2)
EN=table(3)
hard=0.0
C
C calcular current yield stress and hardening rate
C
if (eqplas .eq. 0.0) then
syield =A
else
hard = EN*B*eqplas**(EN-1)
syield=A+B*eqplas**EN
endif
return
end
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
|