-C UMAT程序实例+david注解: g; F* H' {8 Z% @. j% C$ N4 M/ n# f* w" B" D- a& c( w5 S
, _6 f) Z! w/ i( v: G3 ?SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
UMAT子程序(应力张量矩阵,状态变量矩阵,雅可比矩阵[应变增量引起的应力增量],弹性应变能,塑性耗散,蠕变耗散,5 v% V/ p; _! A$ ?3 a1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
1--没有什么实质性的意义吧,大致是要给自己提个醒,这一行的一些量在定义该本构模型中的温度变化时会用到.
RPL--每增量步结束时,由于材料的机械加工[理解成变形比较好]单位时间所产生的体积热[根据帮助文件翻译,大抵是这个本构模型中考虑到材料变形过程中塑性功产生热量,这个参数是对这个热量的一个衡量吧,我不是很懂]
DDSDDT--温度引起的应力增量的变化量
DRPLDE--应变增量引起的RPL的变化量
DRPLDT--温度引起的RPL的变化量
[RPL,DDSDDT,DRPLDE,DRPLDT,四个参数只有在一个完全耦合的热应力分析中才需要定义,而该本构涉及到此,故而引入定义]
STRAN--应变矩阵
DSTRAN--应变增量矩阵
- s! c2 H/ T3 W, w- d$ l8 E% `/ j7 X! o1 {) |4 i
6 @/ j$ q- A: a2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
2--
TIME--时间
DTIME--时间增量
TEMP--增量步开始时的温度
DTEMP--温度增量
PREDEF--增量步开始时,该点处基于结点处的值内插得到的预定义场变量值的矩阵
DPRED--预定义场变量的增量矩阵
MATERAL--==MATERIAL材料名称
NDI--直接应力分量的个数
NSHR--剪切应力分量的个数
NTENS--总应力分量的个数= NDI--直接应力分量的个数+NSHR--剪切应力分量的个数
+ ^4 ^# L0 Y8 F3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
3--
NSTATV--状态变量矩阵的维数
PROPS--材料常数的矩阵
NPROPS--材料常数的个数
COORDS--当前积分点的坐标值
DROT--转动增量矩阵
PNEWDT--新的时间增量对当前时间增量的建议比例
CELENT--特征单元长度
4 `1 W# M/ B% S2 I. s+ N7 v; i4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
7 o, U$ r" c; i) d0 R" v
4--
DFGRD0--增量步开始时的变形梯度数组
DFGRD1--增量步结束时的变形梯度数组
NOEL--单元编号
NPT--积分结点编号
KSLYA--==LAYER[原帮助文件中]层编号[对复杂的壳体和分层固体]
KSPT--当前层的截面点编号
KSTEP--分析步编号
KINC--增量步编号
9 @1 D$ b8 X3 a6 w' {: E1 Y& e: i/ q7 W$ ^( ^- S, K
---------------------------------------------------以上是变量声明
C$ s8 G2 \) a, Z/ k3 C% ^8 xINCLUDE'ABA_PARAM.INC'----[将此文件包含进来,该文件时定义精度的文件,包含双精度和单精度定义], i% K6 i6 U! |# j6 bC! j8 D* A0 u( j% t* U 2 X# f2 Z9 q6 O ~# q2 E9 e/ x4 O* M& b6 K--------------------------------------------------以上是引入包含文件' t3 s! g' D% s# @' ^8 a, [ 7 r3 E/ @3 G" K& a6 n * A. z+ d( w* c 8 f0 g! ]. [/ W! r! Z5 e2 v2 I $ u/ B* }7 v$ h) J& A7 x! @# j 6 k; Z% V. N; A' U5 |$ C f! c0 `4 s1 w+ ^$ `; F1 s' ]CHARACTER*80 MATERL--[定义80个字符长度的字符型变量,命名为MATERL,跟前面对应]6 s ]+ O4 Q: yDIMENSION STRESS(NTENS),STATEV(NSTATV),--[定义数组:应力矩阵[维数],状态变量矩阵[维数],[Fortran程序中用数组存储矩阵数据]& c* V! _, q2 H" I7 X. D1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTEN)--[雅可比矩阵[行数,列数],参照上面]3 H( K U8 q( D. M2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1)--[数组维数说明符的下标下届为1时,可省略]- ^! c' I' c/ h" F* U+ c: X3 PROPS(NPROPS),COORDS(3),DROT(3,3),. y, Z9 D- q) ~$ r4 DFGRD0(3,3),DFGRD1(3,3) q; |3 M W2 S2 N9 g' G * ?( j% i# c( p* E b9 v2 j/ E ! e* V, n/ K |, ] t2 Y$ x# ^2 T& }* d( [( ?9 D6 |& w' ~' `: y0 {. t1 s8 W+ ~. |. d4 [% r - D, ]4 d e. n6 P/ B7 ^4 m--------------------------------------------以上是变量类型定义( @ f: v. D$ [. W: q. W f, f2 s' v4 b1 |! [ N( L# t; z9 d( A+ }3 D! k# F& ]3 M; T. C7 |" d1 r# v: Y5 l S V! |6 `% C. s+ a" p+ d. e% `' `7 R; r# O: b" pC& b1 V I9 \4 \8 n- hDIMENSION EELAS(6),EPLAS(6),FLOW(6)--[弹性应变,塑性应变,流动?]0 j1 _% ~) g B! h7 ZPARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0,HALF=0.5d0)--[PARAMETER语句,指定某些不变量]' T( v1 ^! w: y S0 G1 }* UDATA NEWTON,TOLER/40,1.D-6/--[DATA语句,指定程序中的某些变量或数组的初值]4 P8 g$ D+ a+ e" w& t9 qC& e; s; T: x3 B8 `. E8 Y9 b4 n3 _3 {. h4 G, D6 h$ m% Y1 m3 Z: x# L! n0 v4 r4 {-------------------------------------------用户自定义变量参数5 J( I+ o& n2 W & [, T" |( l' o% ~$ \5 Z9 J- I4 [, g# ]" v3 p& ]9 S) u$ W# y$ m6 a4 [. H p9 s# E 9 |# U4 B- N3 A" |- F! A) b. h9 C b% ]1 AC-----------------------------------------------------------8 B9 p4 b; r% s) R/ C* n2 [& VC UMAT FOR JOHNSON-COOK MODEL4 f; Q$ Y! ] m1 yC-----------------------------------------------------------( W6 t5 A+ v8 N# gC PROPS(1)-YANG'S MODULUS--[杨氏模量]# Y8 C# `2 `$ w( HC PROPS(2)-POISSON RATIO--[泊松比]1 o- ?- B5 I/ M* I! H; _: }) lC PROPS(3)-INELASTIC HEAT FRACTION--[塑性耗散比,表示塑性功转化成为热量的比例]# _& j/ Y( n7 R0 GC PARAMETERS OF JOHNSON-COOK MODEL:--[JOHNSON-COOK模型中的常量声明]( v7 j+ X2 e" x0 S7 a( zC PROPS(4)-A3 w6 ]8 \) A5 E) _ uC PROPS(5)-B2 w# ?& ~% m$ X, A1 C$ dC PROPS(6)-n2 m- D. ~# L6 x% ]- oC PROPS(7)-C/ K5 m9 ^' j$ @+ \9 F" mC PROPS(8)-m+ A9 P; [. [. G* WC-----------------------------------------------------------6 ~- }" I2 t+ QC# e0 j; k; R6 U* p* Q3 l4 f IF(NDI.NE.3)THEN* f0 Y3 d* T: ]4 W4 q4 Q5 S1 u WRITE(6,1)--[输出设备6,输出格式1]; a- V! C. O6 b w1 FORMAT(//,30X,'***ERROR-THIS UMAT MAY ONLY BE USED FOT'--[输出格式\内容]3 H; g' e# d6 z4 t% X 1 'ELEMENTS WITH THREE DIRECT STRESS @/ O9 V6 o$ ? ENDIF-----------------------[终止条件定义]2 z# j1 P" D" z5 X* _# J: T % B5 o% j- [' u . }4 W6 U# z b* \+ W " N) t. {& S' Q / I/ B: O/ F7 E-----------------------------用户子程序头文件,说明程序一些相关信息0 v- t% L$ Q4 l$ u5 f9 v+ |9 c# A* d4 c& ~) Z& x % o7 k1 x( p* P ' m/ I: \5 T) Z% Q9 W% k 0 L5 K3 M5 j/ f9 J. f+ b 4 i, D( N. }' e7 |4 ]) R1 L" @! R' o7 s' }, c1 }7 A k, M+ Y( i6 S! F; |' T* m& Q$ s: F. o. A" B% F7 y; Y) N, kC( z# N; {" n. V0 VC ELASTIC PROPERTIES--[弹性性质]) n- ~0 R" @7 X, |$ G2 \2 JC$ E( h' k2 S/ d) \3 I( GEMOD=PROPS(1)--[将弹性模量赋予EMOD]# J6 [ b" X( C/ X5 FENU=PROPS(2)---[将泊松比赋予ENU]* A' B0 X% B9 c0 b1 c, l3 [IF(ENU.GT.0.4999.AND.ENU.LT.0.5001) ENU=0.499)----------------[定义泊松比的取值]5 J2 d6 b a: d7 K6 Q& sEBULK3=EMOD/(ONE-TWO*ENU))------------------[定义EBULK3=E/(1-2v),为弹性体积膨胀系数K*3,对应弹性力学公式,[弹性体积膨胀系数]体积模量K=1/3*E/(1-2v)]2 B0 `! \! m' y1 e- m* m; cEG2=EMOD/(ONE+ENU)-----------------[定义EG2=E/(1+v)]/ W5 J6 K( |, ~1 z* q: b1 MEG=EG2/TWO-------------------[剪切弹性模量定义,对应于弹性力学中的剪切弹性模量[拉梅弹性常数之一μ=]G=E/2(1+v)] G2 Q4 O$ G" ]1 PEG3=THREE*EG-------------------[--定义EG3=3*EG=3*E/2(1+v)=3G,不知何意??]8 C2 |( K/ z: M$ K/ U9 \ t7 F. GELAM=(EBULK3-EG2)/THREE--------------------[定义ELAM=1/3*{E/(1-2v)-E/(1+v)}=K-2G/3,对应于弹性力学中的拉梅弹性常数入]+ i) m# j4 O K9 {4 H& o% n3 p8 d; A% S2 h- c/ h. A N; L; s, v0 Z# b* c/ F4 {/ g------------------------------以上是根据已经给出的弹性模量和泊松比求解拉梅弹性常数μ和入2 b; f5 t! g" C 3 s3 T' y1 T# f0 W2 n: |( H: B- I# F/ v) b+ r8 u4 h0 F' z , K6 R9 v# S- b z e5 I0 a V7 p& @& Z1 }' b7 a+ i4 @, c: u, O7 q0 x5 v" y6 t 9 F2 ^8 X4 Z3 g# m7 e/ DC/ i3 Z0 r: w* z& l( O/ N9 u/ ^. SC ELASTIC STIFFNESS-------------------[弹性刚度矩阵], c9 c! i2 i( s/ g( \C7 u6 L- Y$ y6 gDO 20 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;20代表此行循环代号]) I# a$ g0 D3 D/ j7 W% p" hDO 10 K2=1,NTENS-------------------[DO循环,初值1,终值为传递过来的NTENS,步长为1,省略;10代表此行循环代号]/ k5 I: L! ?" y' l) j9 ~DDSDDE(K2,K1)=0.0------------------[循环体,将DDSDDE中的全部值赋值为0.0,K2代表行,K1代表列,按列赋值]; a& h3 h/ a& x6 ~' c/ m10 CONTINUE-------------------[代号为10的循环执行结束--continue为终端语句,作为循环终端的标志]# a; G3 [# @5 q2 B0 d+ W20 CONTINUE-------------------[代号为20的循环执行结束--continue为终端语句,作为循环终端的标志]7 C5 c- q- X2 g. |/ a) B0 x! n0 YC/ F$ u! P y6 S0 w; F. c% sDO 40 K1=1,NDI-------------------[DO循环,初值1,终值为NDI[直接应力分量的个数],步长为默认1,省略;40代表此行循环代号]: o3 _2 p) J! N4 {- D! u4 RDO 30 K2=1,NDI-------------------[DO循环,初值1,终值为NDI[直接应力分量的个数],步长为默认1,省略;30代表此行循环代号]: `& x* A9 M% G b2 O, F1 ^DDSDDE(K2,K1)=ELAM------------------[循环体,将DDSDDE中NDI行NDI列的直接应力分量[对应于正应力分量]全部值赋值为拉梅弹性常数入,K2代表行,K1代表列,按列赋值]6 e( i) h$ _- d$ c* P. f30 CONTINUE-------------------[代号为30的循环执行结束--continue为终端语句,作为循环终端的标志]1 A$ h6 \, l* Z( @# g2 y& X7 ?% R* ADDSDDE(K1,K1)=EG2+ELAM------------------[循环体,将DDSDDE中对角线上的NDI个直接应力分量[对应于正应力分量]全部值赋值为拉梅弹性常数EG2+入=E/(1+v),K2代表行,K1代表列,按列赋值]" h/ m6 |. f# v. Q% i! _40 CONTINUE-------------------[代号为40的循环执行结束--continue为终端语句,作为循环终端的标志]4 Y0 m. I- i8 A c# ~8 d8 i9 j3 _* @DO 50 K1=NDI+1,NTENS-------------------[DO循环,初值NDI+1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;50代表此行循环代号]& V! |. ?" ?) A: u" U X8 i- {7 e& \DDSDDE(K1,K1)=EG------------------[循环体,将DDSDDE中NDI---->NTENS的对角线元素的值赋为EG=G,,按列赋值]- `; z4 z, t" V" O% @' A50 CONTINUE-------------------[代号为50的循环执行结束--continue为终端语句,作为循环终端的标志], I: {4 E5 F4 i5 q& W8 S" U1 M4 S- P 1 d4 M( z# s% S- \6 s6 q) n; T" K) P; e# C 2 |$ `0 B/ l" |7 w8 e( I9 J/ |-------------------以上是为了获得弹性本构方程中的弹性本构\弹性模量矩阵,参见弹性力学相关知识[陈惠发弹塑性力学书P105]2 a# s2 y: T! b: A7 ~5 f' r8 Z0 Q% V2 ]$ o8 l' N" {$ Y! c( y. V% ]) @, f" x . Y ?$ b7 _& [( _" X ( Z, v$ }: {" h# ?& ?/ X/ Q& t % j! K; n9 w7 S1 ~$ aC/ } ?4 q1 ?% h8 [/ B7 ]C CALCULATE STRESS FROM ELASTIC STRAINS-------------------[以弹性应变计算应力]. q1 E% X. Y, P UC6 C' @5 t6 z$ @DO 70 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;70代表此行循环代号]3 Z4 S& D" d5 D$ Z J' m2 f% e. J" pDO 60 K2=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;60代表此行循环代号]" ?/ X0 N/ i& `* }, ]0 }) f4 s- J) E" [( V9 XSTRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)-------------------[循环体,三维矩阵形式的本构关系:{应力}=[弹性模量矩阵]{应变};此处的循环很好的实现了,应力列变量的每一个元素,是弹性模量矩阵中对应的一行元素*应变列阵中的一列元素,不错!]; @% v# l8 D3 `& m6 x" P, y s& n# X+ M+ b- j60 CONTINUE-------------------[代号为60的循环执行结束--continue为终端语句,作为循环终端的标志]5 m% w: G" \4 x70 CONTINUE-------------------[代号为70的循环执行结束--continue为终端语句,作为循环终端的标志]6 y7 N" r. W: T* P+ w6 O l3 H% h4 s8 a2 j1 o6 o2 K " _; O, Z" ]) }8 A' C1 w! f1 V: O, M( H: f' T6 D7 n $ o! R6 Q. W% ?7 A-------------------以上是通过弹性阶段的本构方程,由所给的应变计算应力,参见弹性力学相关知识[陈惠发弹塑性力学书P104]5 W) _0 S4 F: g: k0 z* S q' N9 z, D# g! I4 `% r( e: E' D$ o- x6 h7 c - n5 Q: Z' }* [, x) L) u3 Q) t, J0 p+ U' `0 d: j! s; ^5 Q' U: q4 v" J' R6 C. JC8 t# b5 w; |! Z2 i: VC RECOVER ELASTIC AND PLASTIC STRAINS-------------------[弹性和塑性应变覆盖/更新--david自己猜的]6 }: ~0 M+ o* [% J6 ~5 cC) y7 D9 V6 P' ?1 q8 x6 XDO 80 K1=1,NTENS-------------------[DO循环,初值1,终值为NTENS[总的应力分量的个数],步长为默认1,省略;80代表此行循环代号]$ {: P/ ?9 f1 G6 ?EELAS(K1)=STATEV(K1)+DSTRAN(K1)-------------------[弹性应变=状态变量+应变增量,状态变量矩阵中1-6存储弹性应变]8 i; c$ b4 g; K5 U# ]8 OEPLAS(K1)=STATEV(K1+NTENS)-------------------[状态变量矩阵中的7-12存储塑性应变]% H2 n, s7 k3 a+ y, o/ U80 CONTINUE- M( E+ c0 J6 Y$ M% d" P7 B2 o4 pEQPLAS=STATEV(1+2*NTENS)-------------------[状态变量矩阵中的13存储等效塑性应变]' P" ?( I% @3 O# @C7 Q% s6 U1 L a, X5 y4 F9 }7 L+ z: n" W: r- H+ N, R: y- B# X-------------------状态变量在增量步开始时将数值传递到UMAT中.在增量步结束时必须更新状态变量矩阵中的数据' o1 h& B5 g1 }/ Z8 N& v- e; t2 V- Y# Z- F! E6 D7 Y: _- R2 t# _* w. n1 ?C CALCULATE MISES STRESS-------------------[计算MISE屈服准则]' s z' ]; l! G" `( R8 d# N( }C, U9 ]; A1 U' \IF(NPROPS.GT.5.AND.PROPS(4).GT.0.0)THEN-------------------[如果(材料常数的个数>5并且材料常数(4)>0.0),那么]9 n( ~. U2 Z& T- u" ~( y( YSMISES=(STRESS(1)-STRESS(2))*(STRESS(1)-STRESS(2))+" [3 w- M; \ M& [3 a; a1(STRESS(2)-STRESS(3))*(STRESS(2)-STRESS(3))++ z- `6 U8 o- O, a1(STRESS(3)-STRESS(1))*(STRESS(3)-STRESS(1))) A- B! U( E6 U d G4 @DO 90 K1=NDI+1,NTENS( ?$ m S, X' J1 t9 U* B) K" fSMISES=SMISES+SIX*STRESS(K1)*STRESS(K1)& _7 `3 U% a3 C1 z, H: i3 n. a90 CONTINUE' E- d p! Q( z$ A* h" a& L" XSMISES=SQRT(SMISES/TWO)-------------------[MISE屈服应力计算]' m2 U* W9 z% @7 pC$ d8 b( K" X) h/ u6 O" |/ G2 E2 J" h$ N: q$ w5 X, S0 W3 W1 V H-------------------MISE屈服应力计算/ m& Z: {* j; `/ B7 @5 H: s. m$ ?2 |8 B / h( l3 l+ ]+ Z5 O5 E8 z+ e6 [ Y8 K0 w4 X% FC CALL USERHARD SUBROUTINE,GET HARDENING RATE AND YIELD STRESS-------------------[调用用户硬化子程序,获取硬化率和屈服应力]' L* {: O0 P# Z1 L# N; CC3 X$ o: r, V: Y$ V5 |/ _9 a9 \C0 C. C& @% u# D1 n& uCALL USERHARD(SYIEL0,HARD,EQPLAS,PROPS(4))-------------------[调用用户硬化子程序(虚参列表,SYIEL0??)]9 y2 v( `( n2 v1 G& \& IC DETERMINE IF ACTIVELY YIELDING-------------------[确定是否屈服]- z4 y9 M, L9 l% d6 ^" Z* A pC+ Z+ {( z; {& s9 E. [% ^IF(SMISES.GT.(1.0+TOLER)*SYIEL0)THEN: r2 }6 Z2 j2 x3 y, g- i6 EC O* h( J7 w g. j4 cC MATERIAL RESPONSE IS PLASTIC,DETERMINE FLOW DIRECTION-------------------[条件成立,则材料进入塑性,定义流动法则]% r5 E* ?, c) D& h3 L8 CC6 W( y% A, o2 u, U! X+ s" q3 NSHYDRO=(STRESS(1)+STRESS(2)+STRESS(3))/THREE1 X& O" B0 w, {* T( u! AONESY=ONE/SMISES- \+ _: P: g: U$ }5 \' GDO 110 K1=1,NDI3 Y& v6 K6 A3 @1 o1 y" FFLOW(K1)=ONESY*(STRESS(K1)-SHYDRO)) b# o4 r& |* d9 D110 CONTINUE* R5 ~+ l( {3 l( t) C: [) v, RDO 120 K1=NDI+1,NTENS/ t" Y) i- Q2 v% } SFLOW(K1)=STRESS(K1)*ONESY9 Z/ P& ^' U6 ]+ e# ?, z120 CONTINUE1 R$ K; Q: F$ \; I, `, W* H1 r 6 ^. |# Z9 r7 i( o& H ' C6 O0 I% x7 e' w9 f ( V2 H; I% g2 _0 o-------------------流动法则定义,参见弹塑性力学有关书籍* A- h" `4 {7 @ 4 w6 I1 j# ]# E1 l 0 Y) v! f$ r( M- g* }! j/ N. K) d& v2 X" I4 KC2 c; l! N6 Q) j2 iC READ PARAMETERS OF JOHNSON-COOK MODEL-------------------[读取JOHNSON-COOK模型参数]6 v$ H* j g8 y; Z( P8 D! t4 ~! @- R: QC6 s2 h' g; f% Z6 t* UA=PROPS(4)4 a% c! O9 t! W7 e6 ]! l$ `B=PROPS(5)4 H' n3 |6 @+ U% p. x" CEN=PROPS(6)/ k( r* Q( ^# T! y! OC=PROPS(7): S+ e4 k; |# F9 b+ P }EM=PROPS(8)& x: Z2 D* [ K% y( aC. R* L7 |% @. D% U+ w- [, F1 K1 M6 ]- |8 c-------------------读取模型参数,需要用户输入 L! A1 O/ I" C4 C0 i / R/ I# q* ?, v( i6 `/ t : y; p- f: j) {; N* ? ' G8 O8 R) @, a6 r# N# C$ w: O" t# v/ f# CC NEWTON ITERATION-------------------[NEWTON迭代法--数值迭代]5 [* ?# T% x [5 T. {C# X3 _ G# u" xSYIELD=SYIEL0+ C: x* j8 l; |DEQPL=(SMISES-SYIELD)/EG3" r6 }6 [$ Z; L' ?DSTRES=TOLER*SYIEL0/EG34 l; s, F2 t% A! L. J8 F! wDEQMIN=HALF*DTIME*EXP(1.0D-4/C)5 u; K, o7 I* h9 g, @1 GDO 130 KEWTON=1,NEWTON; n, _& u- ~. l. `/ A+ x, `DEQPL=MAX(DEQPL,DEQMIN), f& j' ] x8 t, F7 l$ wCALL USERHARD(SYIELD,HARD,EQPLAS+DEQPL,PROPS(4)+ Y7 } F/ ^' G, q% QTVP=C*LOG(DEQPL/DTIME)0 x- U" K; r# b+ T h! W& ?, CTVP1=TVP+ONE) _5 b* l/ z1 K) U" VHARD1=HARD*TVP1+SYIELD*C/DEQPL. r7 q/ c' P# ^! sSYIELD=SYIELD*TVP1# p% D8 X8 b* MRHS=SMISES-EG3*DEQPL-SYIELD2 O& B. |0 f% k% |6 P; K9 Q1 |! ?DEQPL=DEQPL+RHS/(EG3+HARD1)+ a M6 }( J" q* V6 lIF(ABS(RHS/EG3).LE.DSTRES)GOTO 140( G3 v. }: ]$ P/ @4 v130 CONTINUE% F t$ X# P; T: a: f, |' BWRITE(6,2)NEWTON& p6 @; j+ ?9 I" ?5 F; t& h! `2 FORMAT(//,30X,'***WARNING-PLASTICITY ALGORITHM DID N( S7 s- ?# X" j2 z M1'CONVERGE AFTER',I3,'ITERATIONS'). ~7 m5 O4 z# v! m+ h4 n140 CONTINUE0 x( C: }" o0 J8 ?2 uEFFHRD=EG3*HARD1/(EG3+HARD1)( l; @8 ~: k; n4 e) n6 @2 DC2 {% q! `) S% u4 ^; M7 M % x9 b# S0 [3 j( K5 R3 `1 @" x, d$ h4 Q( Z" e' S q( c-------------------牛顿迭代法定义,参见弹塑性力学有关书籍9 K* h4 F$ O, t1 v: {/ I* E# W* J) q! Q; ^9 k# f9 e% Y* z6 Q" F" m' c! s - d+ d `! r9 t, }7 j2 l: D( B: Y, _/ `3 f+ BC CALCULATE STRESS AND UPDATE STRAINS-------------------[重新计算应力,并更新应变]1 m+ J" Z! h/ [) ?7 o, ?C& b* B4 S( M* M/ M3 k# M, DDO 150 K1=1,NDI. h( Q, H# E! i8 F2 H! h* jSTRESS(K1)=FLOW(K1)*SYIELD+SHYDRO8 N; e5 g5 {5 u' M# Y5 s+ d6 a1 I( IEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL/TWO$ T1 J0 F ^' J9 z: y/ GEELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL/TWO0 M/ d& j6 f8 G0 E; o150 CONTINUE$ q& ?6 u: X/ M( F5 v: ?DO 160 K1=NDI+1,NTENS" A5 F% r7 W. W! W: |STRESS(K1)=FLOW(K1)*SYIELD$ @- T9 j$ q1 r4 eEPLAS(K1)=EPLAS(K1)+THREE*FLOW(K1)*DEQPL, B; W+ q: h! \ }8 ^% e) O: jEELAS(K1)=EELAS(K1)-THREE*FLOW(K1)*DEQPL h9 E/ [- ?- f, p. H( C( U160 CONTINUE6 g) U! s t8 F- Z& k. l4 a2 jEQPLAS=EQPLAS+DEQPL5 R5 J9 K5 w$ v& e! \$ RSPD=DEQPL*(SYIEL0+SYIELD)/TWO: Y, C" a d4 ?RPL=PROPS(3)*SPD/DTIME G1 f4 v/ r, R3 HC) O! T7 L) j E; q" v) t- z; S/ O. u 4 A4 x2 F! G4 a# p* i-------------------以上皆是塑性力学的知识,参见弹塑性力学有关书籍8 F9 D+ u) n& A6 R8 j& h7 f5 Q/ a* L% D- P. e) ]0 _5 j; ~9 A2 w C0 ^ ]$ B& `6 a7 y# \, vC JACOBIAN-------------------[计算雅可比矩阵]) X, S' B5 v5 w5 d( G+ f# `* k+ UC) W" t3 \0 y5 Y3 l! U- t8 AEFFG=EG*SYIELD/SMISES( T# t6 |1 e0 u# e- HEFFG2=TWO*EFFG) o$ t; W! r; c2 g2 REFFG3=THREE*EFFG2/TWO& k. h2 y6 E/ s3 G5 bEFFLAM=(EBULK3-EFFG2)/THREE D- g% K2 D& f* iDO 220 K1=1,NDI- x% A- R- D+ s9 }! e; ?DO 210 K2=1,NDI4 {# }9 ?0 S$ O* wDDSDDE(K2,K1)=EFFLAM: H" \7 L$ ^3 o( [& I2 Q0 ^210 CONTINUE; i; E. f# j7 i( @ zDDSDDE(K1,K1)=EFFG2+EFFLAM6 b8 {6 c0 d$ j# E c: M2 k3 m) R220 CONTINUE( L4 ^, p( ?8 A- u) ODO 230 K1=NDI+1,NTENS' q; S( w/ P' y& x% ]; lDDSDDE(K1,K1)=EFFG# h9 [; \' X6 B$ J% _# P2 y230 CONTINUE3 m+ c+ m" ^; [# [& u MDO 250 K1=1,NTENS" x3 _1 E" e9 k5 k1 Q8 xDO 240 K2=1,NTENS" k- `' c; r2 m: z! LDDSDDE(K2,K1)=DDSDDE(K2,K1)+FLOW(K2)*FLOW(K1)6 g h4 }1 _- @% N/ {" H1*(EFFHRD-EFFG3)8 N& ~* u9 u: h7 a240 CONTINUE5 P# Y9 O2 c( M; f c* Z250 CONTINUE- i9 o! [0 ?9 g& _6 A' K6 pENDIF' f1 e1 p+ \" P6 u. x' ~) SENDIF/ |& Z0 b: O9 `# t3 o6 @C i3 s6 Y" U- i6 Y5 {+ m v+ {6 y6 y; [- G ?" S+ z8 M G1 c' X) }: W- j7 C; d6 r O1 s) {. m* o! V* a( A* `% y# I! T4 R' tC STORE STRAINS IN STATE VARIABLE ARRAY-------------------[存储应变到状态变量矩阵]+ L& o2 F" y$ r+ T* o3 [. \C5 v# z0 [; u, u' ?0 v5 M8 x5 v- EDO 310 K1=1,NTENS- _* B6 y; E+ y: C' {STATEV(K1)=EELAS(K1)6 C5 w, a8 J+ i; Q5 k$ ?& T* hSTATEV(K1+NTENS)=EPLAS(K1)) Y+ ?( o ^! a0 \4 J310 CONTINUE$ ~, I7 r5 _, xSTATEV(1+2*NTENS)=EQPLAS5 s9 N$ A; ?/ W9 mC7 P6 d( }8 P; W* G0 aRETURN I, `, L9 t) k$ M% D+ E K) X8 f& [END-------------------[UMAT退出,增量步结束]" ~4 E3 i; T$ k+ c 2 B- D1 @! H$ g+ F9 t/ n l/ y, M O-------------------结束/ x# M6 s2 ]* _, C! `' e0 r1 ?9 _7 c1 d; ~C: h9 g4 e' q' B5 W3 |# ]C/ t/ V: Q1 F* P Z+ a! }6 d- BC! _; ]- z- C: w8 J! ASUBROUTINE USERHARD(SYIELD,HARD,EQPLAS,TABLE)-------------------[用户自定义硬化子程序--与UMAT类似,把相关硬化定义编进去即可,跟所要定义的本构有关]8 x8 N6 m A3 m* j3 v$ b, {: PC9 H* Q( @2 c' Q7 _: T4 ZINCLUDE'ABA_PARAM.INC') w2 J+ c; A' f5 D4 p/ f; r, HC. x e& o P9 Q8 P2 }1 b" | Y/ bDIMENSION TABLE(3)5 k0 x; K( ] B8 \7 i1 S9 Y2 h1 iC3 |1 B4 H, Z2 jC GET PARAMETERS,SET HARDENING TO ZERO-------------------[获取常数,设初始硬化值=0]- ?/ t. S* S& |9 U$ ^C3 v4 b2 c: i) j* r# ~& H) T5 i' a' d+ A, VA=TABLE(1), a k% x: d3 sB=TABLE(2)( B' C+ I2 b; j& H. a9 F" yEN=TABLE(3)) q5 w8 S& v( u, |" xHARD=0.0& ^, N0 p' Z" b4 Z. g( f+ SC; H# k* m$ J$ @* oC CALSULATE CURRENT YIELD STRESS AND HARDENING RATE-------------------[计算当前屈服应力和硬化率]3 ]* e% t% {3 a8 ]C0 k/ u7 J D8 q. m/ f0 ] `IF(EQPLAS.EQ.0.0)THEN9 X1 @# f v- l0 V) H$ s) M iSYIELD=A3 A7 ]; N. V8 Q+ `' K) eELSE( N8 m* O, g s/ ~HARD=EN*B*EQPLAS**(EN-1)9 V6 P0 ^) t/ [9 Q6 cSYIELD=A+B*EQPLAS**EN. c! ^1 p+ ]% x3 H9 kEND IF# O- q) Y7 t) }$ k7 J" x, F8 ^RETURN% j# o8 H, D: k4 i- sEND
