急急急！求大师关于混凝土hjc模型的ABAQUS的vumat二次开发程...
最近在做混凝土的损伤研究，用到了hjc模型，现有以下两个问题向大神请教（1）在调用子程序时，提交成功，但出现这样的错误，该怎么处理？
The job input file "static-123.inp" has been submitted for analysis.
Job static-123: Analysis Input File Processor completed successfully.
ERROR in job messaging system: Error in connection to analysis
Error in job static-123: The executable D:\SIMULIA\Abaqus\6.10-1\exec\package.exe aborted with system error code 5. Please check the .dat, .msg, and .sta files for error messages if the files exist.If there are no error messages and you cannot resolve the problem, please run the command "abaqus job=support information=support" to report and save your system information.Use the same command to run Abaqus that you used when the problem occurred.Please contact your local Abaqus support office and send them the input file, the file support.log which you just created, the executable name, and the error code.
Job static-123 aborted due to errors.
（2）子程序在没有加入输出调试时，计算结果严重偏离实际情况，怀疑子程序有问题，希望有大神可以帮忙解决一下
现附上子程序文件和inp文件，拜托了，急求，多谢了！ 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=1.d0/3.d0,half=0.5d0,op5=1.5d0,toler=10e-3,
*twothird=2.d0/3.d0,NEWTON=10,fifty=50,hundred=100)
dimension strial(ndir+nshr),ssa(ndir+nshr),etrie(ndir+nshr)
1 ,ssb(ndir+nshr)
integer k,i,iii
double precision e,xnu,A,B,C,en,pcrush,evcrush,plock,evlock,tc,fc,
* k1,k2,k3,SFmax,EFmin,Eodot,dShearModu,d2G,alamda,dBulkModu,
* d3bulk,ddamage,pn,kav,dmu,ptrial,c1,dmup,pt,dT,mua,ev3,D1,D2,
* varj2t,qstar,edoteq,c2,sigmay,pstar,phi,kfactor,dgama
double precision strial,ssa,etrie,ssb,failure
c
c stateNew(k,1)=总应变
c stateNew(k,2)
c stateNew(k,3)
c stateNew(k,4)
c stateNew(k,5)
c stateNew(k,6)
c stateNew(k,7)=弹性应变
c stateNew(k,8)
c stateNew(k,9)
c stateNew(k,10)
c stateNew(k,11)
c stateNew(k,12)
c stateNew(k,13)=ddamage
c stateNew(k,14)=dgama
c
iii=ndir+nshr
e=props(1)
xnu=props(2)
A=props(3)
B=props(4)
C=props(5)
en=props(6)
pcrush=props(7)
evcrush=props(8)
plock=props(9)
evlock=props(10)
tc=props(11)
fc=props(12)
D1=props(13)
D2=props(14)
k1=props(15)
k2=props(16)
k3=props(17)
SFmax=props(18)
EFmin=props(19)
Eodot=props(20)
C
dShearModu=e/(two*(one+xnu))
d2G=two*ShearModu
alamda=xnu*d2G/(one-two*xnu)
dBulkModu=alamda+twothird*dShearModu
d3bulk=three*dBulkModu
c
do k=1,nblock
c
write(7,110)'k', k
110 format(i6,/)
ddamage=stateOld(k,13)
pn=-third*(stressOld(k,1)+stressOld(k,2)+stressOld(k,3))
write(7,120)'pn', pn
120 format(f10.7,/)
c !n平均应力Pn,压为+
kav=(one-ddamage)*dBulkModu+ddamage*k1
dmu=-(strainInc(k,1) + strainInc(k,2) + strainInc(k,3))
write(7,150) 'dmu',dmu
150 format(f10.7,/)
c ! 体积应变增量,压为+
ptrial=pn+kav*dmu
dT=fc*(one-ddamage)
if((ptrial+dT).lt.zero) then
pt=-dT !拉为-,压为+
dmup=zero
else if((ptrial .lt. pn) .or. (ptrial .le. pcrush)) then
c !卸载 or 第一阶段
pt=ptrial
dmup=zero
else if(ptrial .lt. plock) then ! 非卸载 第二阶段
c1=(plock-pcrush)/evlock ! H=(Pl-Pc)/Ul
dmup=(ptrial-pn)/(kav+c1) !塑性体积应变增量
pt=pn+kav*(dmu-dmup)
else ! 第三阶段
mua=dmu-(stateOld(k,1)+stateOld(k,2)+stateOld(k,3))
c !n+1总的体积应变,压为+
c1=(mua-evlock)/(one+evlock)
pt=c1*(k1+c1*(k2+c1*k3))
dmup=zero
endif
c 体积应变计算结束,得到n+1时刻的平均应力pt,塑性体积应变增量dmup
c 计算弹性体积应变/3-->ev3,符号仍是压为+
ev3=pt/d3bulk
c 计算n+1时刻的偏应力.首先计算反对称应变张量产生的偏应力strial(6)
c n时刻的偏应力ssa(6)
ssa(1)=stressOld(k,1)+pn
ssa(2)=stressOld(k,2)+pn
ssa(3)=stressOld(k,3)+pn
ssa(4)=stressOld(k,4)
c if(nshr.gt.one)then
ssa(5)=stressOld(k,5)
ssa(6)=stressOld(k,6)
c endif
c Sx=2(S12*dW12-S13*dW31)Sy=2(-S12*dW12+S23*dW23)Sz=2(-S23*dW23+S31*dW31)
strial(1)=two*(-relSpinInc(k,3)*ssa(4)
*+relSpinInc(k,2)*ssa(6))
strial(2)=two*(relSpinInc(k,3)*ssa(4)
*-relSpinInc(k,1)*ssa(5))
strial(3)=two*(relSpinInc(k,1)*ssa(5)
*-relSpinInc(k,2)*ssa(6))
c Sxy=(S22-S11)dW12+S31*dW23-S23*dW31
strial(4)=(ssa(1)-ssa(2))*relSpinInc(k,3)
1 -ssa(6)*relSpinInc(k,1)+ssa(5)*relSpinInc(k,2)
c Syz=-S31*dW12+（S33-S22)*dW23+S12*dW31
strial(5)=ssa(6)*relSpinInc(k,3)
1 +(ssa(2)-ssa(3))*relSpinInc(k,1)-ssa(4)*relSpinInc(k,2)
cSzx=S23*dW12-S12*dw23+(S11-S33)* dW31
strial(6)=-ssa(5)*relSpinInc(k,3)
1 +ssa(4)*relSpinInc(k,1)-(ssa(3)-ssa(1))*relSpinInc(k,2)
c 记录Trial弹性偏应变,后面仍用到.（n时刻弹性应变+应变增量）
c1=third*(stateOld(k,7)+stateOld(k,8)+stateOld(k,9)-dmu)
etrie(1)=stateOld(k,7)+strainInc(k,1)-c1
etrie(2)=stateOld(k,8)+strainInc(k,2)-c1
etrie(3)=stateOld(k,9)+strainInc(k,3)-c1
etrie(4)=(stateOld(k,10)+strainInc(k,4))/two
etrie(5)=(stateOld(k,11)+strainInc(k,5))/two
etrie(6)=(stateOld(k,12)+strainInc(k,6))/two
c采用弹性Trial偏应变etrie()预估Trial偏应力
do i=1,iii
strial(i)=strial(i)+d2G*etrie(i)
end do
varj2t=etrie(4)**2+etrie(5)**2+etrie(6)**2
1+half*(etrie(1)**2+etrie(2)**2+etrie(3)**2)
qstar=sqrt(three*varj2t)/fc
write(6,180) 'qstar',qstar
180 format(f10.7,/)
c 计算等效偏应变率
c1=dmu/three
ssb(1)=strainInc(k,1)+c1
ssb(2)=strainInc(k,2)+c1
ssb(3)=strainInc(k,3)+c1
ssb(4)=strainInc(k,4)
ssb(5)=strainInc(k,5)
ssb(6)=strainInc(k,6)
if(dt .gt. zero) then
edoteq=sqrt(two*third*(ssb(1)**2+ssb(2)**2+ssb(3)**2
1 +two*(ssb(4)**2+ssb(5)**2+ssb(6) **2)))/(dt*Eodot)
else
edoteq=zero
endif
write(7,*)"strainrate",edoteq
pstar=pt/fc
if(pstar .gt. zero) then
sigmay=(A*(one-ddamage)+B*(pstar**en))*(one+C*log(edoteq))
write(6,190) 'sigmay',sigmay
190 format(f10.7,/)
if((sigmay-SFmax).gt.zero)then
sigmay=SFmax
end if
else
c1=A*(one-ddamage)*(one+C*log(edoteq))
c2=tc/fc*(one-ddamage)
if(abs(c2) .lt. elimit) then
sigmay=zero
else
sigmay=c1/c2*abs(c2-abs(pstar))
end if
end if
phi=qstar-sigmay
if(phi.gt.zero)then
c 塑性屈服
if(qstar .gt. elimit) then
kfactor=sigmay/qstar
else
kfactor=zero
end if
dgama=phi/(three*dShearModu)!等效塑性应变增量
c stateNew(k,14)=stateOld(k,14)+dgama
c计算屈服后的应力
stressNew(k,1)=kfactor*strial(1)-pt
stressNew(k,2)=kfactor*strial(2)-pt
stressNew(k,3)=kfactor*strial(3)-pt
stressNew(k,4)=kfactor*strial(4)
stressNew(k,5)=kfactor*strial(5)
stressNew(k,6)=kfactor*strial(6)
c 计算弹性应变
kfactor=kfactor/d2G
stateNew(k,7)=kfactor*strial(1)-eve3 !弹性体积应变更新
stateNew(k,8)=kfactor*strial(2)-eve3
stateNew(k,9)=kfactor*strial(3)-eve3
stateNew(k,10)=two*kfactor*strial(4)
stateNew(k,11)=two*kfactor*strial(5)
stateNew(k,12)=two*kfactor*strial(6)
else
c 偏应力为弹性,更新应力
stressNew(k,1)=strial(1)-pt
stressNew(k,2)=strial(2)-pt
stressNew(k,3)=strial(3)-pt
stressNew(k,4)=strial(4)
stressNew(k,5)=strial(5)
stressNew(k,6)=strial(6)
c 偏应力为弹性,更新弹性应变:偏应变为弹性,但体积应变不见得是弹性的
stateNew(k,7)=etrie(1)-eve3 !弹性体积应变更新
stateNew(k,8)=etrie(2)-eve3
stateNew(k,9)=etrie(3)-eve3
stateNew(k,10)=two*etrie(4)
stateNew(k,11)=two*etrie(5)
stateNew(k,12)=two*etrie(6)
end if
c结束弹塑性判断,更新总应变
do i=1,iii
stateNew(k,i)=stateOld(k,i)+strainInc(k,i)
end do
c 更新损伤变量
c1=D1*((pt+tc)/fc)**D2
if(c1.lt.EFmin) then
c1=EFmin
end if
ddamage=(dgama+dmup)/c1
if(ddamage .gt. one) then
ddamage=one
end if
write(7,*)"D",ddamage
stateNew(k,13)=ddamage
c
c
if((stateNew(k,13)).lt.(0.d85))then
failure=one
else
failure=zero
end if
stateNew(k,14)=failure
C Update the specific internal energy -
C
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))
enerInternNew(k) = enerInternOld(k) + stressPower / density(k)
C
C Update the dissipated inelastic specific energy -
C
c plasticWorkInc = half * ( yieldOld + yieldNew ) * dgama
c enerInelasNew(k) = enerInelasOld(k)
c * + plasticWorkInc / density(k)
end do
C
return
end
自己顶一个。。。。。。。。。。 顶一个、、、、、、、、、、、、、、、、、、、 inp文件呢？ 你好，请问你有相关参考资料么？ 这是什么材料的本构啊啊
你好，请问一下你调试的方法是怎么样呢？ 可以交流一下吗，我也是做的hjc，计算不收敛，不知道哪不对，可以帮忙看看嘛 有人会吗 可以交流一下吗
页:
[1]