IJMS期刊喷丸诱导晶粒细化的程序分享
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 常量声明 const parameter
real*8,parameter:: x0=1.0d-30,PRECISION_GAMMA=5.0d-3,
1 PRECISION_RHO=5.0d-3,DPEEQ_SHRINK=1.0d-5
parameter(ZERO=0.0d0,ONE=1.0d0,TWO=2.0d0,THREE=3.0d0,HALF=0.5d0,
1 SIX=6.0d0,MAXSTEP=50)
C
real*8 E,mu,Rho_c0,Rho_w0,diameter,dpeeq,peeq,Rho,rho_w,rho_c,
1 rhoc_old,rhow_old,pre_dpeeq,dpeeq_trial,volFrac,K_frac
c 状态变量内容
c state(*,1) 屈服应力 Yeild Stress
c state(*,2) 等效塑性应变 PEEQ
c state(*,3) 位错壁密度 Dislocation Wall Density
c state(*,4) 位错胞密度 Dislocation Cell Density
c state(*,5) 位错胞直径 Dislocation Cell Diameter
c state(*,6) 位错壁体积分数 Volume Fraction
c state(*,7) 总位错密度 Total Dislocation Density
c 用户输入参数赋予 Parameters about Material
E = props(1) !弹性模量
Mu = props(2) !泊松比
Sig1 = props(3) !初始应力S1
Alpha_star = props(4) !α*
Beta_star = props(5) !β*
K_c = props(6) !kc
K_w = props(7) !kw
N_c = props(8) !nc
N_w = props(9) !nw
M_star = props(10) !m*
Burger = props(11) !b 伯氏矢量
Taylor = props(12) !M 泰勒因子
Eta = props(13) !η 常数
volFrac_zero = props(14) !f_0 体积分数参数
volFrac_infin = props(15) !f_infin 体积分数参数
Gamma0 = props(16) !参考剪应变率
K_zero = props(17) !K参数
K_infin = props(18) !K参数
Peeq_wave = props(19) !常数
Beta = props(20) !常数
Rho_c0 = props(21) !初始位错胞密度
Rho_w0 = props(22) !初始位错壁密度
c 拉梅常数 Lame constant
La = (Mu*E)/((ONE+Mu)*(ONE-TWO*Mu)) !lambda
Lb = E/(TWO*(ONE+Mu)) !G
if(stepTime .eq. ZERO)then !初始时间步 initial step
do k=1,nblock
trace = strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
stressNew(k,1) = TWO*Lb*strainInc(k,1) + trace*La
stressNew(k,2) = TWO*Lb*strainInc(k,2) + trace*La
stressNew(k,3) = TWO*Lb*strainInc(k,3) + trace*La
stressNew(k,4) = TWO*Lb*strainInc(k,4)
stressNew(k,5) = TWO*Lb*strainInc(k,5)
stressNew(k,6) = TWO*Lb*strainInc(k,6)
end do
else !非初始时间步 following step
do k=1,nblock !计算试应力 calculating trial stress
trace = strainInc(k,1)+strainInc(k,2)+strainInc(k,3)
s1 = stressOld(k,1) + TWO*Lb*strainInc(k,1) + trace*La
s2 = stressOld(k,2) + TWO*Lb*strainInc(k,2) + trace*La
s3 = stressOld(k,3) + TWO*Lb*strainInc(k,3) + trace*La
s4 = stressOld(k,4) + TWO*Lb*strainInc(k,4)
s5 = stressOld(k,5) + TWO*Lb*strainInc(k,5)
s6 = stressOld(k,6) + TWO*Lb*strainInc(k,6)
!计算mise应力，判断是否屈服 calculating mise stress then check Is yeild?
s m = (s1+s2+s3)/THREE
s1 = s1 - s m
s2 = s2 - s m
s3 = s3 - s m
mise = sqrt(1.5d0 *
1 ((s1*s1 + s2*s2 + s3*s3) + TWO*(s4*s4 + s5*s5 +s6*s6)))
!传递量 variable could be changed in process
sig_y = stateOld(k,1)
peeq = stateOld(k,2)
rho_w = stateOld(k,3)
rho_c = stateOld(k,4)
diameter = stateOld(k,5)
volFrac = stateOld(k,6)
!不传递量 variable havn't changed in process
Rho = StateOld(k,7)
dpeeq = ZERO
if(sig_y .EQ. ZERO)then !在未屈服时为相关变量赋值 Assign values to related variables when unyielding
peeq = ZERO
rho_w = rho_w0
rho_c = rho_c0
!相关变量初值计算 Calculation of initial values of relevant variables
K_frac = K_zero
volFrac = volFrac_zero
Rho = volFrac*rho_w + (ONE-volFrac)*rho_c
diameter = K_frac/sqrt(Rho)
sig_y = sig1 + Taylor*Eta*Lb*Burger*
1 (volFrac*sqrt(rho_w)+(ONE-volFrac)*sqrt(rho_c))
end if
if(mise .LE. sig_y)then !未屈服 试探应力就是实际应力 The test stress at unyielding is the actual stress
!应力计算 stress calculation
stressNew(k,1) = s1 + s m
stressNew(k,2) = s2 + s m
stressNew(k,3) = s3 + s m
stressNew(k,4) = s4
stressNew(k,5) = s5
stressNew(k,6) = s6
!能量计算 energy calculation
strainEnergy = HALF * (
1 (stressNew(k,1)+stressOld(k,1))*strainInc(k,1) +
2 (stressNew(k,2)+stressOld(k,2))*strainInc(k,2) +
3 (stressNew(k,3)+stressOld(k,3))*strainInc(k,3) ) +
4 (stressNew(k,4)+stressOld(k,4))*strainInc(k,4) +
5 (stressNew(k,5)+stressOld(k,5))*strainInc(k,5) +
6 (stressNew(k,6)+stressOld(k,6))*strainInc(k,6)
enerInternNew(k) = enerInternOld(k) +
1 strainEnergy/density(k)
enerInelasNew(k) = enerInelasOld(k)
else !屈服 按位错密度本构模型计算 According to dislocation density constitutive model when yeild
!计算试探等效塑性应变初值 Calculating the initial value of equivalent plastic strain
rhoc_old = rho_c
rhow_old = rho_w
d_old = diameter
!迭代计算增量步开始时的等效应变增量(按牛顿迭代法求解)
!The equivalent strain increment at the beginning of the increment step is computed iteratively (solved by Newton iteration method).
!该增量步作为为初始试探值
!This incremental step is taken as the initial trial stress value
!迭代常数声明 Iterative constant declaration
jstep = ZERO
dpeeq = (mise-Sig1)/(THREE*Lb) !由mise控制的假设迭代初值 mise controlled initial value of hypothetical iteration
!计算迭代初值 Calculating the initial iteration value
f = ONE
do while(f .GT. ZERO)
dpeeq = dpeeq*DPEEQ_SHRINK
f = Sig1-mise+THREE*Lb*dpeeq+Taylor*Eta*Lb*Burger*
1 (volFrac*dsqrt(rho_w)+(ONE-volFrac)*dsqrt(rho_c))*
2 (((Taylor*dpeeq)/(Gamma0*dt))**(ONE/M_star))
end do
pre_dpeeq = -dpeeq
!开始用牛顿迭代计算等效应变增量 The equivalent strain increment is calculated by Newton iteration
do while(abs((pre_dpeeq-dpeeq)/pre_dpeeq) .GT. PRECISION_GAMMA)
f = Sig1-mise+THREE*Lb*dpeeq+Taylor*Eta*Lb*Burger*
1 (volFrac*dsqrt(rho_w)+(ONE-volFrac)*dsqrt(rho_c))*
2 (((Taylor*dpeeq)/(Gamma0*dt))**(ONE/M_star))
df = THREE*Lb + Taylor*Eta*Lb*Burger*
1 (volFrac*dsqrt(rho_w)+(ONE-volFrac)*dsqrt(rho_c))*
2 (((Taylor*dpeeq)/(Gamma0*dt))**((ONE/M_star)-ONE))*
3 (Taylor/(Gamma0*dt))*(ONE/M_star)
pre_dpeeq = dpeeq
dpeeq = dpeeq - f/df
!计数器
jstep = jstep + ONE
if(jstep .GT. MAXSTEP)then
EXIT
end if
end do
dpeeq_trial = -dpeeq
!迭代计算试探位错密度与屈服应力
!The dislocation density and yield stress were tested by iterative calculation
!设置最多迭代步数，以防迭代不收敛时卡死程序
!Set the maximum number of iteration steps to prevent the program from getting stuck when iteration does not converge
kstep = ZERO
do while(kstep .LT. MAXSTEP)!迭代循环 iterative loop
dpeeq_trail = dpeeq
!体积分数 Volume Fraction
volFrac = volFrac_infin + (volFrac_zero-volFrac_infin)*
1 exp((-Taylor*(peeq+dpeeq))/Peeq_wave)
!位错胞密度 Dislocation Wall Density
rho_c = rhoc_old + (
1 dsqrt(rhow_old/THREE)*((Taylor*Alpha_star)/Burger) -
2 (SIX*Beta_star*Taylor)/
3 (Burger*d_old*((ONE-volFrac)**(ONE/THREE))) -
4 K_c*Taylor*rhoc_old*(((Taylor*dpeeq)/(dt*Gamma0))**(-ONE/N_c))
5 ) * dpeeq
!位错壁密度 Dislocation Cell Density
rho_w = rhow_old + (
1 (dsqrt(THREE*rhow_old)*Beta_star*Taylor*(ONE-volFrac))/
2 (volFrac*Burger) +
3 (SIX*Beta_star*Taylor*((ONE-volFrac)**(TWO/THREE)))/
4 (Burger*d_old*volFrac) -
5 K_w*Taylor*rhow_old*(((Taylor*dpeeq)/(dt*Gamma0))**(-ONE/N_w))
6 ) * dpeeq
!总位错密度 Total Dislocation Density
Rho = volFrac*rho_w + (ONE-volFrac)*rho_c
!晶胞直径表达式中的分子 Molecules in the cell diameter expression
K_frac = K_infin + (K_zero-K_infin)*
1 exp(-Beta*Taylor*(peeq+dpeeq))
!晶胞直径 Dislocation Cell Diameter
diameter = K_frac/dsqrt(Rho)
!计算位错更新后的等效塑性应变增量
!Calculate the equivalent plastic strain increment after dislocation updating
!迭代常数声明
!Iterative constant declaration
jstep = ZERO
dpeeq = (mise-Sig1)/(THREE*Lb) !由mise控制的假设迭代初值 mise controlled initial value of hypothetical iteration
!计算迭代初值 Calculate the initial iteration value
f = ONE
do while(f .GT. ZERO)
dpeeq = dpeeq*DPEEQ_SHRINK
f = Sig1-mise+THREE*Lb*dpeeq+Taylor*Eta*Lb*Burger*
1 (volFrac*dsqrt(rho_w)+(ONE-volFrac)*dsqrt(rho_c))*
2 (((Taylor*dpeeq)/(Gamma0*dt))**(ONE/M_star))
end do
pre_dpeeq = -dpeeq
!开始用牛顿迭代计算等效应变增量
!The equivalent strain increment is calculated by Newton iteration
do while(abs((pre_dpeeq-dpeeq)/pre_dpeeq) .GT. PRECISION_GAMMA)
f = Sig1-mise+THREE*Lb*dpeeq+Taylor*Eta*Lb*Burger*
1 (volFrac*dsqrt(rho_w)+(ONE-volFrac)*dsqrt(rho_c))*
2 (((Taylor*dpeeq)/(Gamma0*dt))**(ONE/M_star))
df = THREE*Lb + Taylor*Eta*Lb*Burger*
1 (volFrac*dsqrt(rho_w)+(ONE-volFrac)*dsqrt(rho_c))*
2 (((Taylor*dpeeq)/(Gamma0*dt))**((ONE/M_star)-ONE))*
3 (Taylor/(Gamma0*dt))*(ONE/M_star)
pre_dpeeq = dpeeq
dpeeq = dpeeq - f/df
!计数器 counter
jstep = jstep + ONE
if(jstep .GT. MAXSTEP)then
EXIT
end if
end do
if(abs((dpeeq_trial-dpeeq)/dpeeq_trial)
1 .GT. PRECISION_GAMMA)then
dpeeq_trial = dpeeq
else
EXIT
end if
!计数器 counter
kstep = kstep + ONE
end do!迭代循环结束 End of iteration loop
!屈服应力 Yeild Stress
if(dpeeq .EQ. ZERO)then
sig_y = mise
else
sig_y = Sig1 + (Taylor*Eta*Lb*Burger)*
1 (volFrac*dsqrt(rho_w)+(ONE-volFrac)*dsqrt(rho_c))*
2 (((Taylor*dpeeq)/(dt*Gamma0))**(ONE/M_star))
end if
!根据迭代得出的等效塑性应变计算折减系数
!The reduction coefficient is calculated according to the equivalent plastic strain obtained by iteration
factor = sig_y/mise
!计算应力 Calculating Stress
stressNew(k,1) = s1*factor + s m
stressNew(k,2) = s2*factor + s m
stressNew(k,3) = s3*factor + s m
stressNew(k,4) = s4*factor
stressNew(k,5) = s5*factor
stressNew(k,6) = s6*factor
!计算能量 Calculating Energy
strainEnergy = HALF * (
1 (stressNew(k,1)+stressOld(k,1))*strainInc(k,1) +
2 (stressNew(k,2)+stressOld(k,2))*strainInc(k,2) +
3 (stressNew(k,3)+stressOld(k,3))*strainInc(k,3) ) +
4 (stressNew(k,4)+stressOld(k,4))*strainInc(k,4) +
5 (stressNew(k,5)+stressOld(k,5))*strainInc(k,5) +
6 (stressNew(k,6)+stressOld(k,6))*strainInc(k,6)
plasEnergy = dpeeq*sig_y
enerInternNew(k) = enerInternOld(k) +
1 strainEnergy/density(k)
enerInelasNew(k) = enerInelasOld(k) +
2 plasEnergy/density(k)
end if
stateNew(k,1) = sig_y
stateNew(k,2) = peeq + dpeeq
stateNew(k,3) = rho_w
stateNew(k,4) = rho_c
stateNew(k,5) = diameter
stateNew(k,6) = volFrac
c
stateNew(k,7) = Rho
end do
end if
return
end
