本帖最后由 david840530 于 2010-1-13 23:15 编辑
J-C UMAT程序实例+david注解
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
UMAT子程序(应力张量矩阵,状态变量矩阵,雅可比矩阵[应变增量引起的应力增量],弹性应变能,塑性耗散,蠕变耗散,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 1--没有什么实质性的意义吧,大致是要给自己提个醒,这一行的一些量在定义该本构模型中的温度变化时会用到. RPL--每增量步结束时,由于材料的机械加工[理解成变形比较好]单位时间所产生的体积热[根据帮助文件翻译,大抵是这个本构模型中考虑到材料变形过程中塑性功产生热量,这个参数是对这个热量的一个衡量吧,我不是很懂] DDSDDT--温度引起的应力增量的变化量 DRPLDE--应变增量引起的RPL的变化量 DRPLDT--温度引起的RPL的变化量 [RPL,DDSDDT,DRPLDE,DRPLDT,四个参数只有在一个完全耦合的热应力分析中才需要定义,而该本构涉及到此,故而引入定义] STRAN--应变矩阵 DSTRAN--应变增量矩阵
2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
2-- TIME--时间 DTIME--时间增量 TEMP--增量步开始时的温度 DTEMP--温度增量 PREDEF--增量步开始时,该点处基于结点处的值内插得到的预定义场变量值的矩阵 DPRED--预定义场变量的增量矩阵 MATERAL--==MATERIAL材料名称 NDI--直接应力分量的个数 NSHR--剪切应力分量的个数 NTENS--总应力分量的个数= NDI--直接应力分量的个数+NSHR--剪切应力分量的个数
3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 3-- NSTATV--状态变量矩阵的维数 PROPS--材料常数的矩阵 NPROPS--材料常数的个数 COORDS--当前积分点的坐标值 DROT--转动增量矩阵 PNEWDT--新的时间增量对当前时间增量的建议比例 CELENT--特征单元长度
4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
4-- DFGRD0--增量步开始时的变形梯度数组 DFGRD1--增量步结束时的变形梯度数组 NOEL--单元编号 NPT--积分结点编号 KSLYA--==LAYER[原帮助文件中]层编号[对复杂的壳体和分层固体] KSPT--当前层的截面点编号 KSTEP--分析步编号 KINC--增量步编号
---------------------------------------------------以上是变量声明 C
INCLUDE'ABA_PARAM.INC'----[将此文件包含进来,该文件时定义精度的文件,包含双精度和单精度定义]
C
--------------------------------------------------以上是引入包含文件
CHARACTER*80 MATERL--[定义80个字符长度的字符型变量,命名为MATERL,跟前面对应]
DIMENSION STRESS(NTENS),STATEV(NSTATV),--[定义数组:应力矩阵[维数],状态变量矩阵[维数],[Fortran程序中用数组存储矩阵数据]
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTEN)--[雅可比矩阵[行数,列数],参照上面]
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1)--[数组维数说明符的下标下届为1时,可省略]
3 PROPS(NPROPS),COORDS(3),DROT(3,3),
4 DFGRD0(3,3),DFGRD1(3,3)
--------------------------------------------以上是变量类型定义
C
DIMENSION EELAS(6),EPLAS(6),FLOW(6)--[弹性应变,塑性应变,流动?]
PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0,HALF=0.5d0)--[PARAMETER语句,指定某些不变量]
DATA NEWTON,TOLER/40,1.D-6/--[DATA语句,指定程序中的某些变量或数组的初值]
C
-------------------------------------------用户自定义变量参数
C-----------------------------------------------------------
C UMAT FOR JOHNSON-COOK MODEL
C-----------------------------------------------------------
C PROPS(1)-YANG'S MODULUS--[杨氏模量]
C PROPS(2)-POISSON RATIO--[泊松比]
C PROPS(3)-INELASTIC HEAT FRACTION--[塑性耗散比,表示塑性功转化成为热量的比例]
C PARAMETERS OF JOHNSON-COOK MODEL:--[JOHNSON-COOK模型中的常量声明]
C PROPS(4)-A
C PROPS(5)-B
C PROPS(6)-n
C PROPS(7)-C
C PROPS(8)-m
C-----------------------------------------------------------
C
IF(NDI.NE.3)THEN
WRITE(6,1)--[输出设备6,输出格式1]
1 FORMAT(//,30X,'***ERROR-THIS UMAT MAY ONLY BE USED FOT'--[输出格式\内容]
1 'ELEMENTS WITH THREE DIRECT STRESS
ENDIF-----------------------[终止条件定义]
-----------------------------用户子程序头文件,说明程序一些相关信息
C
C ELASTIC PROPERTIES--[弹性性质]
C
EMOD=PROPS(1)--[将弹性模量赋予EMOD]
ENU=PROPS(2)---[将泊松比赋予ENU]
IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499)----------------[定义泊松比的取值]
EBULK3=EMOD/(ONE-TWO*ENU))------------------[定义EBULK3=E/(1-2v),为弹性体积膨胀系数K*3,对应弹性力学公式,[弹性体积膨胀系数]体积模量K=1/3*E/(1-2v)]
EG2=EMOD/(ONE+ENU)-----------------[定义EG2=E/(1+v)]
EG=EG2/TWO-------------------[剪切弹性模量定义,对应于弹性力学中的剪切弹性模量[拉梅弹性常数之一μ=]G=E/2(1+v)]
EG3=THREE*EG-------------------[--定义EG3=3*EG=3*E/2(1+v)=3G,不知何意??]
ELAM=(EBULK3-EG2)/THREE--------------------[定义ELAM=1/3*{E/(1-2v)-E/(1+v)}=K-2G/3,对应于弹性力学中的拉梅弹性常数入]
------------------------------以上是根据已经给出的弹性模量和泊松比求解拉梅弹性常数μ和入
C
C ELASTIC STIFFNESS-------------------[弹性刚度矩阵]
C
DO 20 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;20代表此行循环代号]
DO 10 K2=1,NTENS-------------------[DO循环,初值1,终值为传递过来的NTENS,步长为1,省略;10代表此行循环代号]
DDSDDE(K2,K1)=0.0------------------[循环体,将DDSDDE中的全部值赋值为0.0,K2代表行,K1代表列,按列赋值]
10 CONTINUE-------------------[代号为10的循环执行结束--continue为终端语句,作为循环终端的标志]
20 CONTINUE-------------------[代号为20的循环执行结束--continue为终端语句,作为循环终端的标志]
C
DO 40 K1=1,NDI-------------------[DO循环,初值1,终值为NDI[直接应力分量的个数],步长为默认1,省略;40代表此行循环代号]
DO 30 K2=1,NDI-------------------[DO循环,初值1,终值为NDI[直接应力分量的个数],步长为默认1,省略;30代表此行循环代号]
DDSDDE(K2,K1)=ELAM------------------[循环体,将DDSDDE中NDI行NDI列的直接应力分量[对应于正应力分量]全部值赋值为拉梅弹性常数入,K2代表行,K1代表列,按列赋值]
30 CONTINUE-------------------[代号为30的循环执行结束--continue为终端语句,作为循环终端的标志]
DDSDDE(K1,K1)=EG2+ELAM------------------[循环体,将DDSDDE中对角线上的NDI个直接应力分量[对应于正应力分量]全部值赋值为拉梅弹性常数EG2+入=E/(1+v),K2代表行,K1代表列,按列赋值]
40 CONTINUE-------------------[代号为40的循环执行结束--continue为终端语句,作为循环终端的标志]
DO 50 K1=NDI+1,NTENS-------------------[DO循环,初值NDI+1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;50代表此行循环代号]
DDSDDE(K1,K1)=EG------------------[循环体,将DDSDDE中NDI---->NTENS的对角线元素的值赋为EG=G,,按列赋值]
50 CONTINUE-------------------[代号为50的循环执行结束--continue为终端语句,作为循环终端的标志]
-------------------以上是为了获得弹性本构方程中的弹性本构\弹性模量矩阵,参见弹性力学相关知识[陈惠发弹塑性力学书P105]
C
C CALCULATE STRESS FROM ELASTIC STRAINS-------------------[以弹性应变计算应力]
C
DO 70 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;70代表此行循环代号]
DO 60 K2=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;60代表此行循环代号]
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)-------------------[循环体,三维矩阵形式的本构关系:{应力}=[弹性模量矩阵]{应变};此处的循环很好的实现了,应力列变量的每一个元素,是弹性模量矩阵中对应的一行元素*应变列阵中的一列元素,不错!]
60 CONTINUE-------------------[代号为60的循环执行结束--continue为终端语句,作为循环终端的标志]
70 CONTINUE-------------------[代号为70的循环执行结束--continue为终端语句,作为循环终端的标志]
-------------------以上是通过弹性阶段的本构方程,由所给的应变计算应力,参见弹性力学相关知识[陈惠发弹塑性力学书P104]
C
C RECOVER ELASTIC AND PLASTIC STRAINS-------------------[弹性和塑性应变覆盖/更新--david自己猜的]
C
DO 80 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;80代表此行循环代号]
EELAS(K1)=STATEV(K1)+DSTRAN(K1)-------------------[弹性应变=状态变量+应变增量,状态变量矩阵中1-6存储弹性应变]
EPLAS(K1)=STATEV(K1+NTENS)-------------------[状态变量矩阵中的7-12存储塑性应变]
80 CONTINUE
EQPLAS=STATEV(1+2*NTENS)-------------------[状态变量矩阵中的13存储等效塑性应变]
C
-------------------状态变量在增量步开始时将数值传递到UMAT中.在增量步结束时必须更新状态变量矩阵中的数据
C CALCULATE MISES STRESS-------------------[计算MISE屈服准则]
C
IF(NPROPS.GT.5.AND.PROPS(4).GT.0.0)THEN-------------------[如果(材料常数的个数>5并且材料常数(4)>0.0),那么]
SMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2))+
1(STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3))+
1(STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1))
DO 90 K1=NDI+1,NTENS
SMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)
90 CONTINUE
SMISES=SQRT(SMISES/TWO)-------------------[MISE屈服应力计算]
C
-------------------MISE屈服应力计算
C CALL USERHARD SUBROUTINE,GET HARDENING RATE AND YIELD STRESS-------------------[调用用户硬化子程序,获取硬化率和屈服应力]
C
C
CALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4))-------------------[调用用户硬化子程序(虚参列表,SYIEL0??)]
C DETERMINE IF ACTIVELY YIELDING-------------------[确定是否屈服]
C
IF(SMISES.GT.(1.0+TOLER)*SYIEL0)THEN
C
C MATERIAL RESPONSE IS PLASTIC,DETERMINE FLOW DIRECTION-------------------[条件成立,则材料进入塑性,定义流动法则]
C
SHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE
ONESY=ONE/SMISES
DO 110 K1=1,NDI
FLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)
110 CONTINUE
DO 120 K1=NDI+1,NTENS
FLOW(K1)=STRESS(K1)*ONESY
120 CONTINUE
-------------------流动法则定义,参见弹塑性力学有关书籍
C
C READ PARAMETERS OF JOHNSON-COOK MODEL-------------------[读取JOHNSON-COOK模型参数]
C
A=PROPS(4)
B=PROPS(5)
EN=PROPS(6)
C=PROPS(7)
EM=PROPS(8)
C
-------------------读取模型参数,需要用户输入
C NEWTON ITERATION-------------------[NEWTON迭代法--数值迭代]
C
SYIELD=SYIEL0
DEQPL=(SMISES-SYIELD)/EG3
DSTRES=TOLER*SYIEL0/EG3
DEQMIN=HALF*DTIME*EXP(1.0D-4/C)
DO 130 KEWTON=1,NEWTON
DEQPL=MAX(DEQPL,DEQMIN)
CALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4)
TVP=C*LOG(DEQPL/DTIME)
TVP1=TVP+ONE
HARD1=HARD*TVP1+SYIELD*C/DEQPL
SYIELD=SYIELD*TVP1
RHS=SMISES-EG3*DEQPL-SYIELD
DEQPL=DEQPL+RHS/(EG3+HARD1)
IF(ABS(RHS/EG3).LE.DSTRES)GOTO 140
130 CONTINUE
WRITE(6,2)NEWTON
2 FORMAT(//,30X,'***WARNING-PLASTICITY ALGORITHM DID N
1'CONVERGE AFTER',I3,'ITERATIONS')
140 CONTINUE
EFFHRD=EG3*HARD1/(EG3+HARD1)
C
-------------------牛顿迭代法定义,参见弹塑性力学有关书籍
C CALCULATE STRESS AND UPDATE STRAINS-------------------[重新计算应力,并更新应变]
C
DO 150 K1=1,NDI
STRESS(K1)=FLOW(K1)*SYIELD+SHYDRO
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO
150 CONTINUE
DO 160 K1=NDI+1,NTENS
STRESS(K1)=FLOW(K1)*SYIELD
EPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL
EELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL
160 CONTINUE
EQPLAS=EQPLAS+DEQPL
SPD=DEQPL*(SYIEL0+SYIELD)/TWO
RPL=PROPS(3)*SPD/DTIME
C
-------------------以上皆是塑性力学的知识,参见弹塑性力学有关书籍
C JACOBIAN-------------------[计算雅可比矩阵]
C
EFFG=EG*SYIELD/SMISES
EFFG2=TWO*EFFG
EFFG3=THREE*EFFG2/TWO
EFFLAM=(EBULK3-EFFG2)/THREE
DO 220 K1=1,NDI
DO 210 K2=1,NDI
DDSDDE(K2,K1)=EFFLAM
210 CONTINUE
DDSDDE(K1,K1)=EFFG2+EFFLAM
220 CONTINUE
DO 230 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EFFG
230 CONTINUE
DO 250 K1=1,NTENS
DO 240 K2=1,NTENS
DDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)
1*(EFFHRD-EFFG3)
240 CONTINUE
250 CONTINUE
ENDIF
ENDIF
C
C STORE STRAINS IN STATE VARIABLE ARRAY-------------------[存储应变到状态变量矩阵]
C
DO 310 K1=1,NTENS
STATEV(K1)=EELAS(K1)
STATEV(K1+NTENS)=EPLAS(K1)
310 CONTINUE
STATEV(1+2*NTENS)=EQPLAS
C
RETURN
END-------------------[UMAT退出,增量步结束]
-------------------结束
C
C
C
SUBROUTINE USERHARD(SYIELD,HARD,EQPLAS,TABLE)-------------------[用户自定义硬化子程序--与UMAT类似,把相关硬化定义编进去即可,跟所要定义的本构有关]
C
INCLUDE'ABA_PARAM.INC'
C
DIMENSION TABLE(3)
C
C GET PARAMETERS,SET HARDENING TO ZERO-------------------[获取常数,设初始硬化值=0]
C
A=TABLE(1)
B=TABLE(2)
EN=TABLE(3)
HARD=0.0
C
C CALSULATE CURRENT YIELD STRESS AND HARDENING RATE-------------------[计算当前屈服应力和硬化率]
C
IF(EQPLAS.EQ.0.0)THEN
SYIELD=A
ELSE
HARD=EN*B*EQPLAS**(EN-1)
SYIELD=A+B*EQPLAS**EN
END IF
RETURN
END |
IP卡
狗仔卡
显身卡
发表于 2010-1-13 22:01:25
］
