abaqus vumat 金属塑性断裂代码

macgician · 发表于 2019-3-23 00:11:43

论文里面看到的代码，想要在abaqus中实现，调试运行了很久，始终没有应力返回，请大神看看是否代码本身有问题啊？整个模型只有kinetic 能量有变化，其他能量都是零。
（后处理中显示selected primary variable is not available...）

SUBROUTINE VUMAT(nblock,ndir,nshr,nstatev,nfieldv,
1 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,
7 stressNew,stateNew,enerInternNew,enerInelasNew)
C
INCLUDE 'VABA_PARAM.INC'
C VARIABLE DEFINITION
C --------------------------------------------------------------------
dimension props(nprops), density(nblock),
1 coordMp(*),
2 charLength(*), strainInc(nblock,ndir+nshr),
3 relSpinInc(*), tempOld(*),
4 stretchOld(*), defgradOld(nblock,ndir+2*nshr),
5 fieldOld(*), stressOld(nblock,ndir+nshr),
6 stateOld(nblock,nstatev), enerInternOld(nblock),
7 enerInelasOld(nblock), tempNew(*),
8 stretchNew(nblock,ndir+2*nshr),
9 defgradNew(nblock,ndir+2*nshr),
9 fieldNew(*),
9 stressNew(nblock,ndir+nshr),
9 stateNew(nblock,nstatev),
9 enerInternNew(nblock), enerInelasNew(nblock)
CHARACTER*80 CMNAME
Double Precision nue, E, Cxy, Cnn, Css, sig1,sig2,sig3,
1 sig4,sig5,sig6,f, dfs1, dfs2, dfs3, dfs4,misesqr,
2 dfs5,dfs6,FF,GG,HH,NN,ep1,ep2,ep3,gamp1,gamp2,gamp3,
3 gampt1,gampt2,gampt3,ept1,ept2,ept3,e0,
4 dflam,lam,eqmises,dD,ds1lam,ds2lam,ds3lam,ds4lam,
5 peeq, peeqt, dpeeq, ds5lam,ds6lam,temp4,mexp,
6 s1, s, st, H1, table(2,25),A,N,Fra,onethird,
7 DamageC,xflag,ef,xi,thetab,ThirdI,temp1,temp2,temp3,
8 C1,C2,C3,f1,f2,f3,smean,eta,thetap
integer km,i
C
C material parameters
C--------------------------------------------------------------------
C Piecewise input of hardening curve [ table(2,*) for plastic strain points,table(1,*) for corresponding yielding stresses]
table(2,1)=0
table(2,2)=0.010
table(2,3)=0.015
table(2,4)=0.012
table(2,5)=0.023
table(2,6)=0.029
table(2,7)=0.034
table(2,8)=0.038
table(2,9)=0.043
table(2,10)=0.048
table(2,11)=0.053
table(2,12)=0.058
table(2,13)=0.063
table(2,14)=0.068
table(2,15)=0.0728
table(2,16)=0.078
table(2,17)=0.083
table(2,18)=0.089
table(2,19)=0.095
table(2,20)=0.10
table(2,21)=0.11
table(2,22)=0.112456
table(2,23)=0.118
table(2,24)=0.3
table(2,25)=0.6
C
table(1,1)=537.33
table(1,2)=552.42
table(1,3)=565.015
table(1,4)=580.95
table(1,5)=593.97
table(1,6)=603.93
table(1,7)=613.67
table(1,8)=620.78
table(1,9)=626.39
table(1,10)=637.22
table(1,11)=644.566
table(1,12)=647.818
table(1,13)=658.494
table(1,14)=664.096
table(1,15)=665.88
table(1,17)=672.064
table(1,19)=682.71
table(1,20)=682.22
table(1,21)=690.669
table(1,22)=694.155
table(1,23)=694.347
table(1,24)=750.279
table(1,25)=803.213
C
C elasticity
E = props(1)
nue = props(2)
C plasticity
A = props(3)
eO = props(4)
N = props(5)
C MMC fracture
C1 = props(6)
C2 = props(7)
C3 = props(8)
C post failure
xflag = props(9)
mexp= props(10)
C plastic anisotropy
FF = props(11)
GG = props(12)
HH = props(13)
NN = props(14)
LL = props(15)
MM = props(16)
C 3D C matrix
Cxy=E*nue/((1-2*nue)*(1+nue))
Cnn=E*(1-nue)/((1-2*nue)*(1+nue))
Css=E/(2*(1+nue))
C Fracture definition
if ( xflag .eq. 1.) then
DamageC = 1.001
else if (xflag .gt. 1.) then
DamageC = xflag
else
DamageC = 2.
end if
C
C BEGIN(integration point specific)
do 1000 km = 1,nblock
C
C initialize state variables
if(totalTime .le. 0.) then
do j=4,10
stateOld(km,j)=0.
end do
stateOld(km,11)=table(1,1)
stateOld(km,3)=1.
stateOld(km,1)=0.
stateOld(km,2)=0.
stateOld(km,12)=0.
stateOld(km,13)=0.
stateOld(km,14)=0.
end if
if(stateold(km,3).eq. 0.) then
statenew(km,3)=0
statenew(km,2)=0
statenew(km,4) =0
statenew (km, 12) =0
statenew (km, 13) =0
statenew (km, 14) =0
goto 1000
end if
peeq = stateOld(km,1)
epl = stateOld(km,5)
ep2 = stateOld(km,6)
ep3 = stateOld(km,7)
gamp1 = stateOld(km,8)
gamp2 = stateOld(km,9)
gamp3 = stateOld(km,10)
s = stateOld(km,11)
C
C--------------------------------------------------------------------
C calculate trial stress state
sigl=stressOld(km,1) + Cnn*strainInc(km,1)
1 + Cxy*(strainInc(km,2) + strainInc(km,3))
sig2=stressOld(km,2) + Cnn*strainInc(km,2)
1 + Cxy*(strainInc(km,1) + strainInc(km,3))
sig3=stressOld(km,3) + Cnn*strainInc(km,3)
1 + Cxy*(strainInc(km,1)+ strainInc(km,2))
sig4=stressOld(km,4) + Css*2*strainInc(km,4)
sig5=stressOld(km,5) + Css*2*strainInc(km,5)
sig6=stressOld(km,6) + Css*2*strainInc(km,6)
C
f=sqrt(FF*(sig2-sig3)**2 + GG*(sig3-sig1)**2
1 + HH*(sig1-sig2)**2 + 2*NN*sig4**2 +3*sig5**2
2 + 3*sig6**2) - s
if (f.lt.1.e-2) then
C elastic step
ept1 = ep1
ept2 = ep2
ept3 = ep3
gampt1 = gamp1
gampt2 = gamp2
gampt3 = gamp3
st = s
peeqt = peeq
dpeeq = 0.
else
C plastic step
C ----------------------------------------------------------------------
C begin return mapping
C initialize iteration variables
lam =0.
dfs1 =0.
dfs2 =0.
dfs3 =0.
dfs4 =0.
dfs5 =0.
dfs6 =0.
C
do i=1,100
C calculate plastic strains using associated flow rule and plastic
C multiplier 'lam'
ept1 = ep1 + lam*dfs1
ept2 = ep2 + lam*dfs2
ept3 = ep3 + lam*dfs3
gampt1 = gamp1 + lam*dfs4
gampt2 = gamp2 + lam*dfs5
gampt3 = gamp3 + lam*dfs6
C calculate stresses
sig1=stressOld(km,1) + Cnn*(strainInc(km,1) - lam*dfs1)
1 + Cxy*(strainInc(km,2)-lam*dfs2+ strainInc(km,3)-lam*dfs3)
sig2=stressOld(km,2) + Cnn*(strainInc(km,2)-lam*dfs2)
1 + Cxy*(strainInc(km,1)-lam*dfs1+strainInc(km,3) -lam*dfs3)
sig3=stressOld(km,3) + Cnn*(strainInc(km,3)-lam*dfs3)
1 + Cxy*(strainInc(km,1)-lam*dfs1+strainInc(km,2) -lam*dfs2)
sig4=stressOld(km,4) + Css*(2*strainInc(km,4)-lam*dfs4)
sig5=stressOld(km,5) + Css*(2*strainInc(km,5)-lam*dfs5)
sig6=stressOld(km,6) + Css*(2*strainInc(km,6)-lam*dfs6)
C hardening modulus
dpeeq = lam
peeqt = peeq + dpeeq
call UHARD(peeqt, table, 25, s1, H1)
if((stateOld(km,2).gt.1.).and.(xflag .gt. 1.)) then
st=s1*((DamageC-stateOld(km,2))/(DamageC-1.))**mexp
C Hl=Hl*((DamageC-stateOld(km,2))/(DamageC-1.))**0.5
C -s1*0.5*((DamageC-stateOld(km,2))/(DamageC-1.))**(-0.5)
C /(ef*(DamageC-1.))
else
st=s1
end if
C evaluate yield condition
f=sqrt(FF*(sig2-sig3)**2 + GG*(sig3-sig1)**2
1 + HH*(sig1-sig2)**2 + 2*NN*sig4**2
1 + 3*sig5**2+3*sig6**2) - st
if (abs(f).lt.1.e-2) goto 12
C calculate partial derivatives of f : dfs1=df/dsigl, dfs2=df/dsig2...
dfs1 = (-GG*(sig3-sig1)+HH*(sig1-sig2))/(f+st)
dfs2 = (FF*(sig2-sig3)-HH*(sig1-sig2))/(f+st)
dfs3 = (-FF*(sig2-sig3)+GG*(sig3-sig1))/(f+st)
dfs4 = 2*NN*sig4/(f+st)
dfs5 = 3*sig5/(f+st)
dfs6 = 3*sig6/(f+st)
C calculate partial derivatives of s : ds1lam=ds1/dlam, ds2lam=ds2/dlam.
ds1lam=-Cnn*dfs1-Cxy*(dfs2+dfs3)
ds2lam=-Cnn*dfs2-Cxy*(dfs1+dfs3)
ds3lam=-Cnn*dfs3-Cxy*(dfs1+dfs2)
ds4lam=-Css*dfs4
ds5lam=-Css*dfs5
ds6lam=-Css*dfs6
C calculate dflam = df/dlam=df/ds * ds/dlam
dflam = - (H1) + (dfs1*ds1lam+dfs2*ds2lam+
1 dfs3*ds3lam+dfs4*ds4lam +dfs5*ds5lam+dfs6*ds6lam)
C calculate updated plastic multiplier
lam = lam -f/dflam
C end return mapping
C----------------------------------------------------------------------
end do
12 continue
end if
C--------------------------------------------------------------------
C UPDATE
C--------------------------------------------------------------------
C stresses
stressNew(km,1) = sig1
stressNew(km,2) = sig2
stressNew(km,3) = sig3
stressNew(km,4) = sig4
stressNew(km,5) = sig5
stressNew(km,6) = sig6
C other state variables
stateNew(km,1) = peeqt
stateNew(km,5) = ept1
stateNew(km,6) = ept2
stateNew(km,7) = ept3
stateNew(km,8) = gampt1
stateNew(km,9) = gampt2
stateNew(km,10)= gampt3
stateNew(km,11)= st
C--------------------------------------------------------------------
C UPDATE for fracture criteria
C--------------------------------------------------------------------
misesqr=0.5*(sig2-sig3)**2 + 0.5*(sig3-sig1)**2
1 + 0.5*(sig1-sig2)**2 + 3*sig4**2
2 +3*sig5**2+3*sig6**2
if(misesqr.le.0.) then
eqmises=0.001
else
eqmises=sqrt(misesqr)
end if
smean=(sig1+sig2+sig3)/3
temp1=(sig1-smean)*((sig2-smean)*(sig3-smean)-sig5**2)
temp2=sig4*(sig4*(sig3-smean)-sig5*sig6)
temp3=sig6*(sig4*sig5-sig6*(sig2-smean))
temp4=(27./2.)*(temp1-temp2+temp3)
if(temp4.lt.0) then
ThirdI=-(-temp4)**(1.0/3.0)
else
ThirdI=(temp4)**(1.0/3.0)
end if
if(eqmises.eq.0.) then
eta=1./3.
xi=0
else
eta=smean/eqmises
xi=(ThirdI/eqmises)**3.
end if
if(xi.lt.(-1))then
xi=-1
else if(xi.gt.1.)then
xi=1
end if
thetab=1-2*acos(xi)/3.1415927
stateNew(km,12)=smean
stateNew(km,13)=eqmises
stateNew(km,14)=eta
stateNew(km,15)=thetab
stateNew(km,16)=xi
C--------------------------------------------------------------------
C--Apply Fracture Criteria
C--------------------------------------------------------------------
if((xflag .ge. 1.).and.(peeqt.gt.0.)) then
if ((stateOld (km,2) .lt. DamageC).and.
1 (stateOld(km,3).eq.1.)) then
C Mohr-Coulomb fracture criteria
f1 = cos(thetab*3.1415927/6)
f2 = sin(thetab*3.1415927/6)
f3 = C3+(sqrt(3.)/(2.-sqrt(3.)))*(1-C3)*(1./f1-1.)
ef =((A/C2)*f3*(f1*sqrt((1+C1**2)/3.)
1 +C1*(eta+f2/3.))) **(-1./N)
C Cut-off value of stress triaxiality
if (eta .lt. -1) ef=10000.0
C Fracture locus term
Fra = 1./ef
dD = Fra * dpeeq
stateNew(km,2)= stateOld(km,2) + dD
stateNew(km,4)= ef
end if
C ---------------------------------------------
C Check element deletion and delete damaged element
if ((stateNew(km,2) .ge. DamageC) .and.
1 (stateOld(km,3) .eq. 1.)) then
stateNew(km,3) = 0.
stressNew(km,1)= 0.
stressNew(km,2)= 0.
stressNew(km,3)= 0.
stressNew(km,4)= 0.
stressNew(km,5)= 0.
stressNew(km,6)= 0.
end if
C end of fracture coding
end if
C--------------------------------------------------------------------
C USER SUBROUTINE CODE ENDS HERE
C--------------------------------------------------------------------
1000 continue
RETURN
END
C ----------------------------------------------------------------------
subroutine UHARD(eqplasr, table, nvalue, syieldw, hardw)
include 'vaba_param.inc'
C
C variable definition
dimension table(2,nvalue)
Double Precision eqplas0, eqplas1, syield0, syield1, eqplasr,
1 syieldw,hardw,table
C
C eqplasr ............. equivalent plastic strain (read only)
C table ............... array of stress-strain pairs (read only)
C nvalue .............. number of stress-strain pairs (read only)
C syieldw ............. von Mises yield stress (write only)
C hardw ............... slope of yield surface at eqplas (write only)
C
C table(1,1), table(2,1), table(1,2), table(2,2)
C stress_1, eps_1, stress_2, eps_2C
do k1=1, nvalue-1
eqplas1 = table(2,k1+1)
if(eqplasr.lt.eqplas1) then
eqplas0 = table(2,k1)
syield0 = table(1,k1)
syield1 = table(1,k1+1)
hardw = (syield1 - syield0) / (eqplas1 - eqplas0)
syieldw = syield0 + hardw * (eqplasr - eqplas0)
goto 11
end if
end do
hardw = 0.d0
syieldw = table(1,nvalue)
11 continue
return
end

复制代码

zhaojv920402 · 发表于 2019-3-26 16:38:09

不懂帮顶

juwendy · 发表于 2019-5-1 09:41:21

给你一个赞，very good

慌张的小铁匠 · 发表于 2020-2-24 11:06:32

不懂帮顶

战仙 · 发表于 2020-3-24 12:52:37

你好楼主，你这个程序可以计算材料得疲劳吗，咱可以交流一下

llzyq · 发表于 2020-7-4 16:25:14

求论文链接

SecretZe · 发表于 2020-10-20 18:20:40

楼主，你好，能给论文连接吗

ZXY140910abc · 发表于 2020-12-7 21:07:55

小白问一句，楼主调试出来了吗，问一下编好了Vumat子程序之后，不需要给参数赋值吗

努力学仿真 · 发表于 2020-12-10 22:30:33

求论文链接

pcstudying · 发表于 2021-3-24 10:13:37

求论文链接

Hasol · 发表于 2021-3-25 15:10:30

谢谢分享

jcz1996 · 发表于 2021-4-2 20:27:00

求论文连接

蝶恋花 · 发表于 2021-4-6 12:33:53

学习一下

g423423113 · 发表于 2023-4-26 16:38:29

来学习一下

ductility · 发表于 2023-5-17 11:40:47

哇，第一次见这么好的人，也不要仿真币

ductility · 发表于 2023-7-27 09:48:32

求论文链接

ductility · 发表于 2023-7-27 11:12:51

这些代码都是咋哪找的啊

HEAs · 发表于 2024-3-22 21:17:26

求论文名称

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

[二次开发] abaqus vumat 金属塑性断裂代码