Vumat代表各向同性完美塑料

你好,

我用fortran 77编写了VUMAT子程序，以获得完美的塑性，与米斯塑性统一。当对一个简单的立方体模型使用单向拉伸模拟进行分析时…在我看来，它只经历了弹性变形，而不是塑性变形……我试图找出是什么错误..但是找不到..有人能找到这个代码中的错误是什么..我已经附上了VUMAT代码供参考…非常感谢……

子例程vumat (

c只读-

1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal，

2 stepTime, totalTime, dt, cmname, coordMp, charLength，

3个道具，密度，strainInc, relSpinInc，

4 tempOld, stretchOld, defgradOld, fieldOld，

5 stress sold, stateOld, enerInternOld, enerInelasOld，

6 tempNew, stretchNew, defgradNew, fieldNew，

c只写-

7 . stressNew, statennew, enerInternNew, enerInelasNew)

c

包括“vaba_param.inc”

C

维道具(nprops),密度(nblock),

1 coordMp(nblock，*)， nElement(nblock)， nMatPoint(nblock)，

2 charLength(*)， strainInc(nblock,ndir+nshr)，

3 relSpinInc(*)， tempOld(*)，nLayer(nblock)，

4 . stretchOld(*)， defgradOld(*)，nSecPoint(nblock)，

5 fieldOld(*)， stressOld(nblock,ndir+nshr)，

6 stateOld(nblock,nstatev)， enerInternOld(nblock)，

7 enerInelasOld(nblock)， tempNew(*)，

8 stretchNew(*)， defgradNew(*)，

9 fieldNew(*)， stressNew(nblock,ndir+nshr)，

1 stateNew (nblock nstatev),

2 enerInternNew(nblock)， enerInelasNew(nblock)

C

字符* 80 cmname

C

参数(0 = 0。D0, 1 = 1。D0, 2 = 2。D0，一半=。5d0，

& 3 = 3。D0, op5 = 1.5d0)

C在分析开始时假设纯弹性材料

C

E =道具(1)

Xnu =道具(2)

Yield1 = props(3)

Pls1 =道具(4)

twomu = E /(1 + xnu)

alambda = xnu * twomu /(两人* xnu)

Thremu = op5 * twomu

c If (steptime .eq. 0) then

做K = 1,nblock

跟踪= strainInc (k, 1) + strainInc (k, 2) + strainInc (k, 3)

stressNew (k, 1) = stressOld (k, 1)

& + twomu*strainInc(k,1)+alambda*trace

stressNew (k, 2) = stressOld (k, 2)

& + twomu*strainInc(k,2)+alambda*trace

stressNew (k, 3) = stressOld (k, 3)

& + twomu*strainInc(k,3)+alambda*trace

stressNew (k, 4) = stressOld (k, 4) + twomu * strainInc (k, 4)

stressNew (k, 5) = stressOld (k, 5) + twomu * strainInc (k, 5)

stressNew (k, 6) = stressOld (k, 6) + twomu * strainInc (k, 6)

最后做

c其他

做K = 1,nblock

跟踪= strainInc (k, 1) + strainInc (k, 2) + strainInc (k, 3)

C试验压力

s11 = stressOld (k, 1) + twomu * strainInc (k, 1)

& + alambda*trace

s22 = stressOld (k, 2) + twomu * strainInc (k, 2)

& + alambda*trace

s33 = stressOld (k, 3) + twomu * strainInc (k, 3)

& + alambda*trace

s12 = stressOld (k, 4) + twomu * strainInc (k, 4)

向= stressOld (k, 5) + twomu * strainInc (k, 5)

s23 = stressOld (k, 6) + twomu * strainInc (k, 6)

C

Smean = third*(s11+s22+s33)

C

S11 = S11 - smean

S22 = S22 - smean

S33 = S33 - smean

C

Vmises = sqrt(op5*((s11*s11)+(s22*s22)+(s33*s33))

& + (2 * s12 * s12) +(2 *向*向)+ (2 * s23 * s23)))

C print*，'value is'，vmises

C

sigdif = vmises-yieldNew

facyld = 0

如果(sigdif.eq.vmises)然后

因子= 1

C更新重音

其他的如果(sigdif.lt.vmises)然后

yieldNew = yield1

facyld =一个

Deqps = facyld * sigdif /thremu

C更新重音

factor = yieldNew / (yieldNew + thremu * deqps)

c

stressNew(k,1) = s11 * factor + smean

stressNew(k,2) = s22 * factor + smean

stressNew(k,3) = s33 * factor + smean

stressNew(k,4) = s12 * factor

stressNew(k,5) = s13 * factor

stressNew(k,6) = s23 * factor

如果

C

C更新状态变量

C

statennew (k,1) = stateOld(k,1) + deqps

C

C更新比内能

C

C

压力功率=一半* (

&(应力sold (k,1) +应力new (k,1)) *应力inc (k,1) +

&(应激sold (k,2) +应激new (k,2)) *应激inc (k,2) +

&(应激sold (k,3) +应激new (k,3)) *应激inc (k,3))+

& (stress sold (k,4) + stressNew(k,4)) * strainInc(k,4) +

& (stress sold (k,5) + stressNew(k,5)) * strainInc(k,5) +

& (stress sold (k,6) + stressNew(k,6)) * strainInc(k,6)

C

enerInternNew(k) = enerInternOld(k) + (stressPower/density(k))

c

C更新耗散的非弹性比能

C

plasticWorkInc = yield1* deqps

enerInelasNew(k) = enerInelasOld(k)

& + plasticWorkInc /密度(k)

最后做

C结束if

返回

结束

C

问候,