VUMAT基本知识汇总.docx

1、VUMAT基本知识汇总NBLOCK:在铜用Vumat吋需要用到的材料点的数量Ndir：对称张量中直接应力的数量(sigmal 1 ,sigma22,sigma33)Nshr：对称张量中间接应力的数x(sigmal2, sigmal3, sigma23)Nstafcv:与材料类型相关朕的用户定义的状态变量的数目Nficldv：用户定义的外场变量的个数Nprops:用户自定义材料属性的个数Laimcal:指示是否在退火过程中被讷用例程的标志。LanncaI=0,指示在常规力 学性能增量，例程被讽用。Laimcall表示，这是退火过程，你应该重新初始化 内部状态变量，stcpTimc:步骤开始后的数

2、值totalTimc:总时间Dt：时间增量值Cmnamc:用户自定义的材料名称，左对齐。它是通过字符串传迷的。一些内部 材料模型是以“AEQ_”字符串开头给定的名称。为了避免冲突，你不应该在“cmnamc”中使用“ABQ_”作为领先字符串。c(x)rdMp(nblock,*):材料点的坐标值。它是壳单元的中层面材料点，梁和管(pipe) 单元的质心。char Length (nblock):特征元素长度，是基于几何平均数的默认值或用户子程序VUCHARLENGTH中 定义的用户特征元长度。props(nprops):用户使用的材料属性 density(nblock):中层结构的物质点的当前密度

3、rclSpinlnc (nblock, nshr):在随转系统中定义的每个物质点处熠加的相对旋转矢量 system. Defined as 一 where W is the antisymmetric part of the velocitygradient, L, and fl = R - RT Stored in 3D as ( 32,13,21) and in 2D as( 21)_tcmp()ld(nblock):物质点开始增加时的温度。dcfgradOld (nb()ck,ndir+2*nshr):在增量开始时，每个物质点出的变形梯度张量， 在 3d 中形为(Fll, F22, F3

4、3, F12, F23, F31, F21, F32, F13)，在 2d 中形为(Fll, F22, F33, F12, F21)strctchOld (nblock, ndir+nshr)Stretch tensor, U, at each material point at the beginning of the increment defined from the polar decomposition of the deformation gradient by F = R U.ficld()ld (nblock, nficldv):在增量开始时，每个物质点处用户定义场变量的值 s

5、trcssOld (nblock, ndir+nshr):在增量开始时，每个物质点处的应力张量： stalcOld (nblock, nstatcv):在增量开始时，每个物质点处的状态变量： tcmpNcw(nbk)ck):在增量结束时，每个物质点处的温度dcfgradNcw (nbk)ck,ndir+2*nshr):在增量结束吋，每个物质点出的变形梯度张量， 在 3d 中形为(Fll, F22, F33, F12, F23, F31, F21, F32, F13)，在 2d 中形为(Fll, F22, F33, F12, F21) ficldNcw (nblock, nficldv):在增量

6、开始时，每个物质点处用户定5L长变量的值As a simple example of the coding of subroutine VUMAT, consider the generalized plane strain case for an elastic/plastic material with kinematic hardening The basic assumptions and definitions of the model are as followsLet Obe the current value of the stress, and define Sto be t

7、he deviatoric part of the stress The center of the yield surface in deviatoric stress space is given by the tensor 0(, which has initial values of zero. The stress difference, , is the stress measured from the center of the yield surface and is given by = S a .The von Mises yield surface is defined

8、aswhere Ois the uniaxial equivalent yield stress The von Mises yield surface is a cylinder in deviatoric stress space with a radius ofFor the kinematic hardening model, /?is a constant The normal to the Mises yield surface can be written asWe decompose the strain rate into an elastic and plastic par

9、t using an additive decomposition：The plastic part of the strain rate is given by a normality conditionwhere the scalar multiplier must be determined A scalar measure of equivalent plastic strain rate is defined byThe stress rate is assumed to be purely due to the elastic part of the strain rate and

10、 is expressed in terms of Hooke, s law bya = Atrace(c f)I + where Aand /are the Lams constants for the material.Theevolution law for CXis given asa =:纭hq、where H is the slope of the uniaxial yield stress versus plastic strain curveDuring active plastic loading the stress must remain on the yield sur

11、face, so that Q ： Q = 1“The equivalent plastic strain rate is related to byThe kinematic hardening constitutive model is integrated in a rate form as follows A trial elastic stress is computed asG盅* = oid + Atrace(Ae)I + 2“ 氐where the subscripts oldand newrefer to the beginning and end of the increm

12、ent, respective!y. If the trial stress does not exceed the yieldstress, the new stress is set equal to the trial stress If the yield stress is exceeded, plasticity occurs in the increment We then write the incremental analogs of the rate equations asW =甬鬻-=毗-2 口 4 Q、2anew = aold 可 HwhereA7 = 7Ai.Fro

13、m the definition of the normal to the yield surface at the end of the increment, Q,This can be expanded using the incremental equations asaoid + + |(70Q = s伫倍一 Aj2Q.Taking the tensor product of this equation with Q, using the yield condition at the end of the increment, and solving for “7:丁 = 2p(l+H

14、/3“)The value for is used in the incremental equations to determinencwy and 仁etrsubroutine vumat(C Read only 一1nblock, ndir, nshr, nstatev, nfieldv, nprops, 1anneal2stepTime, totalTime, dt, cmname, coordMp, charLength,3props, density, straininc, relSpinlnc,Cinclude vaba_param incCC J2 Mises Plastici

15、ty with kinematic hardening for planeC strain caseC Elastic predictor, radial corrector algorithm.C The state variables are stored as:stressstressstressC All arrays dimensioned by (*) are not used in this algorithm dimension props(nprops), density(nblock),1coordMp(nblock, *),2charLength(*), strainin

16、c(nblock, ndir+nshr),3relSpinlnc (*), %eiiy)01d(*),4stretchOld (*), defgrad01d(*),5fieldOla 1 1. 6stateOld(nblock, nstatev), enerlnternOld(nblock),7enerlnelasOld(nblock), tempNew(*),8stretchNew(*), defgradNew(*), fieldNew(*),9stressNew(nblock, ndir+nshr), stateNew(nblock, nstatev),1 enerlnternNew(nb

17、lock), enerlnelasNew(nblock)Ccharacter*80 cmnameCparameter ( zero = 0., one = 1,two = 2.、three = 3.,1third = one/three, half =5， twoThirds = two/three,2threeHalfs 二 1 5 )props (1)xnuprops (2)yield = props (3)hard = props (4)ctwomu = e / ( one + xnu ) thremu = threeHalfs * twomu sixmu = three * twomu

18、* ( e - twomu ) / ( sixmu - two * e ) (twomu * ( one + hard/thremu )do 100 i = 1,nblockC the yield stress and facyld will be zero.dsmag = dsmag + ( one - facyld )cC Calculated increment in gamma (this explicitly includes theC time step)diff = dsmag 一 radiusdgamma = facyld * term * diffcC Update equi

19、valent plastic straindeqps = coni * dgammastateNew(i, 5) = state01d(i, 5) + deqpscC Divide dgamma by dsmag so that the deviatoric stresses areC explicitly converted to tensors of unit magnitude in theC following calculationsdgamma = dgamma / dsmagC*dslds2ds3s4C Update the specific internal energy 一

20、stressPower = half * (1 ( stress01d(i, l)+stressNew(i, 1) )*strainlnc(i, 1)1 + ( stress01d(i, 2)+stressNew(i, 2) )*strainlnc(i, 2)1 + ( stress01d(i, 3)+stressNew(i, 3) )*strainlnc(i, 3)1 + two*( stress01d(i, 4)+stressNew(i, 4) )*strainlnc(i, 4)cenerlnternNew(i) = enerlnternOld(i)1 + stressPower / de

21、nsity(i)CC Update the dissipated inelastic specific energy 一 plasticWorklnc = dgamma * half * (1 ( stress01d(i, 1)+stressNew(i, 1) )*dsl1 + ( stress01d(i, 2)+stressNew(i, 2) )*ds21 + ( stress01d(i, 3)+stressNew(i, 3) )*ds31 + two*( stress01d(i, 4)+stressNew(i, 4) )*s4 )enerlnelasNew(i) = enerlnelasOld(i)1 + plasticWorklnc / density(i)100 continuecreturnend