子程序subroutinevumat.doc

《子程序subroutinevumat.doc》由会员分享，可在线阅读，更多相关《子程序subroutinevumat.doc（17页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、!-subroutine vumat(1 k4 A$ K- v. c3 F! YC Read only -9 n) Z5 Z3 D3 v* C |0 F F 1nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,. n, O. _1 b* c# K& 2stepTime, totalTime, dt, cmname, coordMp, charLength, p3 k+ : C( b 3props, density, strainInc, relSpinInc,8 w8 |/ W& ?0 4 K4 B/ x- T# 4tempOld, s

2、tretchOld, defgradOld, fieldOld,1 m f7 r- z f; O f7 p& K9 i 3stressOld, stateOld, enerInternOld, enerInelasOld,! a. s+ r7 % Z& s; F. r e 6tempNew, stretchNew, defgradNew, fieldNew,; J# H( & - s! yC Write only -0 p; 0 N1 & Z4 j( V 5stressNew, stateNew, enerInternNew, enerInelasNew )+ _6 e1 m% - K/ B$

3、 f include vaba_param.inc% R# b# z% + e; _8 m8 C S$ z! E3 Q6 T+ _ Z- PC! * Q: i- t; ! V6 b LC MODEL FOR THE SHEAR DEFORMATION AND CRACK OF POLYMERS! R$ M/ o! C- Z+ A5 v C SHEAR MODEL: SPRING+BURGERS ELEMENT+PLASTIC ELEMENT(D-P & EYRING)$ % l# O g q8 m7 i3 _C CRACK MODEL: S11S11_CR WHEN SKK0 A9 p! B5

4、 / U+ N! IC : W0 T/ J3 F. G3 R9 6 hC WRITTEN BYZHANG CHUNYU(CHUNYUNUS.EDU.SG) b6 w3 e. ( yC , N3 I) 2 H8 Z) y4 l* u2 ZC All arrays dimensioned by (*) are not used in this algorithm+ S s3 Rp, y6 G+ T - t% X( r$ m; e- T, P! h dimension props(nprops), density(nblock),$ p5 ! t8 i2 m9 e C 1coordMp(nblock

5、,*),) F/ P/ ! Y) R9 O. c 2charLength(*), strainInc(nblock,ndir+nshr),0 Y d- l9 n: D. P% S$ r 3relSpinInc(*), tempOld(*),8 9 h# H/ w6 Q. G: A 4stretchOld(*), defgradOld(nblock,ndir+nshr), T( W _$ O5 P) g2 U 7 X 5fieldOld(*), stressOld(nblock,ndir+nshr),$ F, e. n/ z6 6stateOld(nblock,nstatev), enerInt

6、ernOld(nblock),* d6 H! A, Q# P 7enerInelasOld(nblock), tempNew(*), m6 P8 K3 K/ f 8stretchNew(*), defgradNew(nblock,ndir+nshr), fieldNew(*),$ K, H! N3 ; 2 L# 9stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),7 N2 bQ f- H% * c4 H8 a 1enerInternNew(nblock), enerInelasNew(nblock) l8 $ V! ?0 TC stra

7、in components stored as state variables! P& I) N2 b4 C, h! c& H# F dimension eelas(ndir+nshr),eplas(ndir+nshr),deplas(ndir+nshr) c; z9 y3 e N/ K* 1 j5 h dimension evisco(ndir+nshr),devisco(ndir+nshr),veint(ndir+nshr)$ y, c J/ Q2 z dimension effstrn(ndir+nshr),ps(ndir),an(ndir,ndir),s(ndir+nshr)8 r2

8、x! A5 F2 dimension dr(3,3),f_old(6),f_new(6)8 r& Q+ z# 4 D5 S9 m# N0 1 e( |$ Y6 , A$ b data newton,toler/10,1.0E-6/2 dU/ & e2 J: G2 b- A* C$ Z9 x Y# b& x* R6 Z! 1 P character*80 cmname1 y; a5 L. w0 a& Y0 a v( DC, w2 E% Y8 W5 c9 W8 U parameter( zero = 0., one = 1., two = 2., three = 3.,7 V! 0 d2 h; a

9、% l7 6 ) b 1third = one/three, half = .5, twoThirds = two/three,( |. U* |) d: i$ j a 2threeHalfs = 1.5 ). w* o6 M! F, L; rC -5 % F8 vU5 OC PROPS(1) - E% E, & G _M; o0 h) m) GC PROPS(2) - NU( L( i2 e1 G6 _ d# CC PROPS(3) - E15 a7 Y* V& a3 E0 iC PROPS(4) - ETA1- k1 k Z8 e; Q 7 X1 D. i+ lC PROPS(5) - E

10、TA02 L/ 6 E8 P/ 5 p1 g: lC PROPS(6) - PHAI2 ?9 m! K9 5 hC PROPS(7) - PHAI2 FOR NON-ASSOCIATIVE FLOW. R0 j) F d/ h. e0 X* D2 S) d: JC PROPS(8) - B FOR EYRING MODEL7 K4 / g: F7 r# p% SC PROPS(9) - C1/ n2 S U( s% : Y! lC PROPS(10)- C2/ f, n# v! i6 D2 p; cC PROPS(11)- M% X/ z! D/ E$ j T& h; S5 aC PROPS(

11、12)- SCRAZE8 x5 P$ N+ 1 r) X; m) 5 2 yC PROPS(13)- ZETA0% w( v$ d- A3 H) z) L2 SC PROPS(14)- CRAZE STRAIN LIMIT8 P3 0 i& m: ; m/ N* mC PROPS(15)- SHEAR STRAIN LIMIT, J4 L/ C4 _2 d: R, W7 JC PROPS(16)- Y0 & g6 e7 I j4 G6 X. iC PROPS(17). PLASTIC CURVE; w1 % |% x) : E- ( rC CALLS AHARD FOR CURVE OF SY

12、IELD VS. PEEQ: Z2 ?9 L; a; C -: y; X4 a/ W9 t$ E: lC& h! N# o- Z& L4 i if (ndir .NE. 3) then$ g8 z3 Z! t- m write(6,1)$ R; 0 f* n5 v6 o1 u 1 format(/,30X,*ERROR - THIS VUMAT MAY ONLY BE USED FOR ,) O, D7 V/ K5 E; M T3 C 1 ELEMENTS WITH THREE DIRECT STRESS COMPONENTS)+ _6 w+ y& A# v# k9 j1 D+ i endif

13、+ V* J+ D- m: w ) S- C- e4 l) f: Z+ s- aC ELASTIC PROPERTIES7 b, p. i. P: % T% l7 4 W ntens=ndir+nshr $ N0 K5 5 A( Z 6 q) T Q! c/ 0 v bniu=props(2), e d* ! 8 O; s V c& A% if(bniu.GT.0.4999.AND.bniu.LT.0.5001) bniu=0.4990 P q, E5 2 3 s, T8 va - sd/ 2 h- u& T$ |! 1 z* ) d% pC PARAMETERS FOR THE KELVIN

14、 UNIT 3 ! d, O: D) w; 1 2 Z9 V8 W# l( H1 L phai=props(6)6 , E: $ k+ f8 s, E: D phai2=props(7)z; J6 W# / y T2 f N B=props(8)+ ! M, ! $ a& u6 U& t* y c1=props(9)3 p. & B/ s, b c2=props(10)% M( a3 & P: _, H: g c3=three*bniu/(one+bniu)6 o4 , B1 * U! q( 4 h& Y bm=props(11)8 H- s7 F+ E5 k2 P& r; B 9 z& v,

15、 J* B scraze=props(12)& V s% k% R; y zeta0=props(13) 5 a4 t) 6 R# k6 # H% I9 J% v2 h& Q craze_cr=props(14)- z, - + |8 G2 I/ a2 shear_cr=props(15)% j, Q; w7 g3 g4 N7 o * i ?c loop the material block; , m! m# Q1 * E7 d$ ) T b( . g% |- c. _( s2 b 4 2 w# y4 v( C$ do 350 i = 1,nblock% L- N- A9 A* N1 W5 T

16、: N8 A! n! a$ f H- E1 S, l: G c_factor=stateOld(i,4*ntens+4)7 _/ t$ a0 , e) z( M( z if(c_factor .EQ. 0)c_factor=1.0+ Y6 ?$ . Y3 q. Q) I( m5 l U* E q_factor=c_factor M- O8 e7 A9 & i1 L5 4 U7 l2 Q1 x7 z- h- a e0=props(1)ks: v* c* Q: _- Y e1=props(3)& O8 J; . P: h: b- b3 m, F |6 eta1=props(4)a F e; Z7

17、L0 ei0 A* 8 eta0=props(5)6 g& x+ E- $ W% 2 W/ Wr7 s i6 x3 N3 O e0=e0*q_factor. VL; z# C& a7 E9 t. Y& to bk=e0/(one-two*bniu)/three% b+ P8 w/ P ?0 l$ Yo g=e0/(one+bniu)/two0 y$ w: h1 y- _/ z e1=e1*q_factor+ ) o0 xd- A( j3 p g1=e1/(2.0*(1.0+bniu): n i f5 k! q; A) p/ h eta1=eta1*q_factor; ) t) I( |( B2

18、 j$ % D eta0=eta0*q_factor7 H# w& B/ W) ? coef1=1.0+eta1/eta0+1.0/4.0*g1/eta0*dt; v i( M8 W% P0 coef2=g1+2.0*eta1/dtC6 ! c/ J* a coef3=1.0/4.0*g1/eta0*dt0 % l4 4 p) e! H8 A$ b% ) e coef=1.0+g*coef1/coef21 n( P( J8 C& f& F* N5 s |; JC RECOVER AND ROTATE DEVIATORIC SHEAR STRAIN & TOTAL STRESS g% p9 A:

19、 ?6 $ nc CALL ROTSIG(STATEV(1),DROT,EELAS,2,NDI,NSHR)( H8 X) 7 F$ y5 d, z9 nc CALL ROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR) S# b+ c4 ?2 M# ?c CALL ROTSIG(STATEV(2*NTENS+3),DROT,EVISCO,2,NDI,NSHR) 9 y* X! _- & a+ m& Q) K6 U do 10 j=1,ntens/ . h8 M- J; i. W/ B/ N eelas(j)=stateOld(i,j) q/ X& e6 %

20、 F! T5 I( r5 E; k eplas(j)=stateOld(i,j+ntens) B$ D( C! l! F: % |2 p0 Y: p: t evisco(j)=stateOld(i,j+3+2*ntens)* |5 t7 T5 & D% o7 B veint(j)=stateOld(i,j+3+3*ntens)& ! Y; O3 : I R V deplas(j)=0.0# e( X4 p M; T# _. m1 F, v8 E f_old(j)=defgradOld(i,j)- D% f2 W. c& l: T+ o8 G q f_new(j)=defgradNew(i,j)

21、9 u A( w& q; f 10 continue ) E/ z2 / f( U5 sc get rotation tensor and rotate old strain tensor H$ _% r) n+ P2 X5 call getdr(f_old,f_new,dr), / p* y) U6 v9 D call ROTTENSOR(eelas,dr,3,3)1 q, V9 F0 V! |! / m- a8 ? call ROTTENSOR(eplas,dr,3,3) Q/ X1 r: m3 s* j call ROTTENSOR(evisco,dr,3,3)0 r& R3 G4 !

22、, K8 q6 J% C eqplas_shr=stateOld(i,1+2*ntens)% G8 f. P6 e( Q eqplas_crz=stateOld(i,2+2*ntens)2 w9 b/ _* - 8 d7 q8 w% vplas=stateOld(i,3+2*ntens)I9 o0 z |* t4 v0 r4 T$ n ( |S3 q9 I6 W+ U7 zC CALCULATE DEVIATORIC STRESS$ U7 Pz0 t% o oldstresskk=0.01 r- n& a: P3 j dstrankk=0.00 g8 V1 b- t2 a+ M strankk

23、=0.0 L5 n! Y3 : X do 20 j=1,ndir5 B4 B. C7 x) Q. C; R oldstresskk=oldstresskk+stressOld(i,j)A! j! m4 s& n9 X4 e dstrankk=dstrankk+strainInc(i,j). S1 r9 1 a- ?6 S- a3 U/ F# T strankk=strankk+eelas(j)+eplas(j)+evisco(j)+strainInc(i,j)4 + d& Q9 oC0 * 20 continue8 x; c! P) 9 J, n) TC CALCULATE EFFECTIVE

24、 STRAIN * R: y* q6 h8 W( + w) x do 30 j=1,ndir$ n. Z3 n2 n3 a5 0 Y v4 f/ c effstrn(j)=eelas(j)+(strainInc(i,j)-dstrankk/3.0)-x# uD. t* ?7 z4 W, u$ z 1 (0.5*(coef1+2.0*coef3)*(stressOld(i,j)-oldstresskk/3.0)+ L9 y% Z! l. 1 g1/eta0*veint(j)-2.0*g1*evisco(j)/coef29 7 b7 P. P5 v; F5 O% a30 continue# f0

25、H) : g$ - v2 h# o do 31 j=ndir+1,ntensX) wG+ Ej. F. |/ effstrn(j)=eelas(j)+strainInc(i,j)/2.0-3 X2 Z/ c1 p: i$ v6 : U 1 (0.5*(coef1+2.0*coef3)*stressOld(i,j)+; o F; e2 m. Y+ c* h 1 g1/eta0*veint(j)-2.0*g1*evisco(j)/coef27 w, U, k1 P; I& ?3 n( ? 31 continue; O3 r2 k7 Z) |) R! S! S& NC UPDATED DEVIATO

26、RIC STRESS 0 , c& w+ Q2 yn do 40 j=1,ntens, u0 t# m( m/ * R, b: d stressNew(i,j)=2.0*g/coef*effstrn(j)2 t t, a- d8 k) Y1 40 continue 9 V- j$ P- M2 a: _rC CALCULATE BULK STRESS & PRESSURE8 X8 0 R7 D M8 U2 ( b d stresskk=oldstresskk+3.0*bk*dstrankk+ N EZ# j* e8 H7 C p=-stresskk/3.0 # P, B/ C; S$ J* T&

27、 C h3 o$ w1 NC UPDATE DIRECT STRESS- 6 q: q$ p8 T6 n6 d) P& o do 45 j=1,ndir D; i: X c) X stressNew(i,j)=stressNew(i,j)-p1 K G/ h9 1 $ e45 continue5 n. L- Y0 j# A; F0 e0 v1 pC DETERMINE SHEAR YILED OR CRAZE YIELD& p6 i3 t. m6 s1 o0 j do 50 j=1,ntens3 D m8 |8 H3 F( B s(j)=stressNew(i,j), L% d. b# % s

28、9 Q9 D% p% I( $ F50 continue* D% B n* i3 DC PRINCIPAL STRESS AND * L# o# g c2 Y call PRINCIPAL(s,ps,an,ndir,nshr,1.0E-6)0 c% s0 U7 5 8 _ ps1=ps(1)9 Z6 c% f; 4 ( n, Ft nmax=1# s/ U! 4 E# i/ N! X if(ps(2) .GT. ps1) then* y s9 w. G# ?2 F. ? ps1=ps(2). o. o2 ) N c$ X nmax=2( o6 |. h: 2 f endif) / u/ D+

29、g* j$ f! R Y# x. K5 y if(ps(3) .GT. ps1) then# Y; K7 F7 P: : u- p% Z2 e3 G6 q2 F ps1=ps(3)1 p& P# S A9 o9 Z# T/ l: m* I I5 A nmax=35 W3 X* z; s! M# 7 endif1 h8 l) M9 y2 O! : R0 H% o% _ bms=-p9 j) x* F! o8 E! X cs=c1+c2/bms+c3*bms) d+ K( Q! D6 s+ q if (bms .GT. 0 .AND. ps1 .GT. 0 .AND. ps1 .GT. cs) t

30、hen: N6 P. Y4 J- H2 / b$ c flow=one% R9 |* M/ N! A$ X1 K( U else ?2 E4 7 X+ J E flow=zero# F/ Q% F1 J$ J0 z, v: P; i& o + L endif# T; z) 5 H O: s$ O3 C9 P! u/ T1 O4 N: V8 aC CRAZE FLOW,calculate deplas% G. wt6 F- i5 5 t) _( o* j if (flow .EQ. one .AND. c_factor .GT. 1.0E-4) then3 B! ) S: e% I& # pd7

31、 f do 60 j=1,ndir- h( U9 # |. f$ j0 ! 3 b) deplas(j)=zeta0*(ps1/scraze/q_factor)*(1.0/bm)*dt*% a# _9 # |! d; , r 1 an(j,nmax)*2# q# G+ Q+ o5 % G4 1 X9 r3 r5 C 60 continue w# ) M$ K. p/ C0 + q; F l2 t do 61 j=ndir+1,ntens+ t% S( & l0 d# H# n1=j-ndir# b1 ?! k+ G: l/ G* y n2=j-ndir+1* u S5 & |0 o+ t6 J

32、0 Y, l& F if(n2 .GT. ndir) then% - ) s) . T, W- Nd6 I0 R n2=n2-ndirZ3 X, _ r& V2 |4 w endif - ; r2 ! P* T( deplas(j)=zeta0*(ps1/scraze/q_factor)*(1.0/bm)*dt*% 7 R u q6 K) D9 Pg 1 an(n1,nmax)*an(n2,nmax), 3 _ t5 K( L: R7 8 i* D61 continue3 t; q* Q: Z4 c- a9 h* t: pn eqplas_crz=eqplas_crz+zeta0*(ps1/s

33、craze/q_factor)*(1.0/bm)*dt8 l! b3 R5 Q2 t/ x- M# l if(eqplas_crz .GT. craze_cr .AND. c_factor .GT. zero) then5 i5 J+ Q2 Z0 B& + f, ?5 Q c_factor=c_factor*0.1*craze_cr/eqplas_crz$ w8 W+ r y7 K. F endif6 x% p& q& K1 z2 g7 A) c + l d; p9 M6 K8 j( mc update stress* A% 4 b1 FP7 g9 c* x. u* z4 r M) v3 M!

34、 M$ e_ dvplaskk=deplas(1)+deplas(2)+deplas(3)5 Y1 - Q w0 x2 F, w c vplas=vplas+dvplaskk9 j% D* O% b stresskk=oldstresskk+3.0*bk/q_factor*c_factor*9 * V/ a4 t, q4 e2 U5 f: A& P 1 (dstrankk-dvplaskk); g a# s* p9 g! M o6 f do 70 j=1,ndir/ F0 O5 z. : Kh deplas(j)=deplas(j)-dvplaskk/3.0- P$ P( r 4 F stre

35、ssNew(i,j)=2.0*g/q_factor*c_factor*4 w- Y. F7 n& Z2 O! W2 I# C; D 1 (effstrn(j)-deplas(j)/coef, u4 u% i$ S: H X8 n_7 ; I stressNew(i,j)=stressNew(i,j)+stresskk/3.0 T! & Q0 c* G7 eR0 f+ Y 70 continue9 E2 A4 X# B, E5 g; do 71 j=ndir+1,ntens$ % t9 R/ h! z. uF, 1 3 K4 n stressNew(i,j)=2.0*g/q_factor*c_f

36、actor*( b0 X. X% r4 T# l0 k8 W3 X0 u 1 (effstrn(j)-deplas(j)/coef# 0 Y) u4 B6 w0 w X9 71 continue 7 |8 S5 V$ Q2 V, q! U endif+ Y- C$ H( d+ Q% G7 Q E . w+ 6 S1 i9 y# R& P3 WC SHEAR YIELD FLOW,CALCULATE DEPLASg8 T1 z q9 j. Kc test crack-shear yielding+ G3 . Y8 E0 E+ Z& g+ c if crack and shear yielding

37、 can not coexist,delete the following the command6 w0 T$ x- e) R# s: K5 lc flow=zero / F. n; x* S# r8 | if(nprops.GT.15 .AND. props(16).GT. 0.0 .AND. flow .EQ. zero) then7 K9 8 x- G0 Z. r* + 7 vC3 _9 _8 B4 a! $ H* q$ N. M6 vC MISES STRESS g( o+ O+ B& x C1 P# i# x1 |4 F3 D& ?9 % T5 M smises=(s(1)-s(2

38、)*(s(1)-s(2) +& L* k9 w8 W0 9 d2 g 1 (s(2)-s(3)*(s(2)-s(3) +3 vo) 6 k* Q2 A 1 (s(3)-s(1)*(s(3)-s(1)9 e9 SB, J1 A5 K% _, M( b do 90 j=ndir+1,ntens( d# G. i1 |: P: z* y3 d smises=smises+6.0*s(j)*s(j), c3 P+ cD) d; W 90 CONTINUE+ c/ a r# i! k4 M5 U, E smises=sqrt(smises/2.0)$ d7 l RE4 p7 e% zC SHEAR ST

39、RAIN RATE9 L$ k& p& S, ! r9 K( d x effrate=0.0 Z7 H; j! c, 7 H, R& E do 95 j=1,ndir8 g7 B9 f6 b9 0 F effrate=effrate+(strainInc(i,j)-dstrankk/3.0)*5 q3 N; Z 9 RJ8 S9 n 1 (strainInc(i,j)-dstrankk/3.0)6 % R/ k , m3 + a! k4 % 95 continue M1 j& m3 n( D4 e O: g do 100 j=ndir+1,ntens; F/ m, ) L; y/ O& p: C$ B effrate=effrate+2.0*strainInc(i,j)*strainInc(i,j)# l$ Q7 E9 k6 p4 e3 7 100 CONTINUE ) c- o/ J$ p! m: b D0 | effrate=sqrt(effrate/2.0)/dt* _: p! p/ m; l! f( M