材料失效及其模拟(页 1) - A04：ABAQUS - SimWe仿真论坛 |Simwe.com|仿真|设计|有限元|虚拟仪器

chunyu 发表于 2005-2-1 10:25

材料失效及其模拟

在ABAQUS里，模拟材料失效主要有两种方法：一是使用已有的失效材料模型，在定义材料时指明最大剪切应变或最大拉应力失效准则，并给出相关的参数，详细内容可参见ABAQUS Analysis User's Manual /11.2.8 Failure models:dO!_9vY*u d
；另一种方法是使用VUMAT，在该子程序中定义一个状态变量（state variable）来表征材料是否失效。为0，失效，为1，正常。然后根据自己定义的失效准则，来给这个变量赋值。子程序方法比较灵活，可以定义自己的失效准则和材料本构方程。详细内容参见ABAQUS Analysis User's Manual /25.3.4 VUMAT。下面几幅截图是我自己的一个VUMAT产生的结果，用来模拟聚合物材料如PMMA的剪切及断裂破坏，本构方程及断裂准则均为自己定义，考虑了弹性，粘弹性，塑性，粘塑性。网格比较粗，示意而已。试件单向拉伸。@E g*Ry!u:I|T
（1）应力集中

chunyu 发表于 2005-2-1 10:27

Re:材料失效及其模拟

裂纹形成

chunyu 发表于 2005-2-1 10:27

Re:材料失效及其模拟

裂纹扩展

chunyu 发表于 2005-2-1 10:28

Re:材料失效及其模拟

进一步扩展

chunyu 发表于 2005-2-1 10:28

Re:材料失效及其模拟

完全断裂

qiuzhixue 发表于 2005-2-1 11:08

Re:材料失效及其模拟

SALUTE。。。
能否发个INP和。F来看看？

sr 发表于 2005-2-1 17:25

Re:材料失效及其模拟

做的不错！5Hw}9MA
请chunyu 老大贴上你的inp和.f看看。

gsut123 发表于 2005-2-2 15:15

Re:材料失效及其模拟

完全断裂前就应该算不下去了，不知这么老兄怎么设置的？

ganyx 发表于 2005-2-3 22:37

Re:材料失效及其模拟

好东西，kq\ {h2I:kcSX
估计大家更想知道的是更详细实现的过程

chunyu 发表于 2005-2-4 17:04

Re:材料失效及其模拟

谢谢大家的关注和鼓励。inp文件就不多说了，一个正方形的板，中间钻个洞。
至于材料破坏的算法，这里采用了B.P. Gearing, L. Anand / International Journal of Solids and Structures 41 (2004) 3125–3150中推荐的方法。原理很简单:W0Cx L"B,kC6Pk,a
（1）根据应力状态，计算主应力。如果hydrostatic stress为正（就是受拉了），最大主应力为正且大于一个值（材料参数），破坏开始形成（即initiation），并在最大主应力方向上产生一个塑性变形，文章里给出了一个公式计算塑性变形的大小（记为P1），当然自己也可以提出合适的公式。否则剪切变形占主导，采用一般的本构方程计算应力应变(等效塑性应变记为P2)。@*Udad0A M7n}"X
（2）当p1>p1_cr或p2>p2_cr时（p1_cr,p2_cr均为材料参数），积分点失效，在VUMAT中通知ABAQUS删除就行了（这是原方法），当一个单元内所有的积分点失效后，该单元自动失效。
（3）为了计算稳定，并且更符合实际破坏的过程，我做了一些修改，不是在p1>p1_cr或p2>p2_cr时立即将单元删除，而是定义了一个我自己称之为质量因子（q_factor）的自变量，当p1>p1_cr或p2>p2_cr后，q_factor很快衰减，材料的模量等抵抗变形的参数也随之快速衰减，直至没有能力承受载荷，然后删除积分点，F;B^8Y+_L4}
（4）几点说明S8sR&\i:J&D-L
(a)破坏的initiation,growth直至材料breakdown，基本上没有统一的标准，所以可以根据实际材料的性能自己提出标准；
(b)这种failure model只能在abaqus/explicit中使用，所以计算效率可能不令人满意；
(c)过年了，少上网，多休息：），恭祝各位鸡年大吉，万事如意！

fracture 发表于 2005-9-8 21:39

Re:材料失效及其模拟

做的太好了，能否发个INP和。F来看看？非常感谢！

都市猎人 发表于 2005-9-9 17:42

Re:材料失效及其模拟

上面的大虾能给俺发歌INPUT文件吗damage@lut.cn.xiexie !!!!+

lololocad 发表于 2005-9-9 17:48

Re:材料失效及其模拟

欧也要啊，hmpen23@126.com,先谢了！

dudreams 发表于 2005-9-10 18:25

Re:材料失效及其模拟

ding,我看大伙都很关注，我建议楼主把inp传到论坛上吧

zzqqzz19770621 发表于 2005-10-28 14:51

Re:材料失效及其模拟

你好，我也正在做材料的失效及断裂，能否将您的inp和.f给我看一下呢？M"M)@R9d.x1ZH3E3O
zzq19770621@126.com
不胜感激！

andy77 发表于 2005-12-19 18:55

Re:材料失效及其模拟

我也想要一个, 谢谢!N4go:`*w _ y1~
panxd@126.com

skipperyhf 发表于 2006-1-12 20:24

Re:材料失效及其模拟

看到了大家都是作材料失效的阿，我们做失效的可不可以专门开个论坛阿，我是做结构失效分析的，现在在学abaqus,不知道该如何做了，向各位请教，我要摹拟支撑在地震作用小的失效，主要是支撑的疲劳破坏，轴向位移达到一定的程度就认为它破坏了，现在这个疲劳失效准则我不知道该在abaqus中如何考虑，望各位高人指教K-RJ n SZ @!m
还有搂主把你的给我发一份可以么，
skipperyhf@tom.com

步行去远方 发表于 2006-3-13 21:01

Re:材料失效及其模拟

楼主可以给发个INP学习一下吗？8c(s$G,g*lgH f
一直没太学会怎么定义失效，着急呀。
谢谢了！6z_|7a!APK
信箱：lmsalms@sohu.com

zhuhanhan 发表于 2006-7-25 08:36

我也是做失效的8W%O"u7k5fO%z
我也要一份：[email]hanfue@163.com[/email]

wuthing2008 发表于 2006-7-25 12:04

我也要一份，谢谢了！
[email]wuthing@sohu.com[/email]

zhuhanhan 发表于 2006-8-3 17:47

请问chunyu (chunyu):你是用的11.2.8 Failure models中的模型？

我想采用材料Progressive Damage and Failured?S6nRe&W6v
但是我的abaqus6.5材料库里面怎么没有该材料？我该怎么看到explicit库中的材料？

xzp831125 发表于 2006-8-24 16:56

呵呵楼主做的很不错Rd)W7T$~+\!Q2_
给我也来一份吧yDQMve~
[email]xuzhiping1125@gmail.com[/email]

zhuhanhan 发表于 2006-8-28 11:31

咱们这里做失效的同志们，有谁作出来了？可不可以将你的成功经历给我们大家说一下，我觉得失效很难做。

步行去远方 发表于 2006-8-28 15:41

楼主好象没再看过这个贴子,现在这么火了,却不见来了,眼下做失效的朋友好象也没几个出来帮忙喔,期待呀!

wjjerry 发表于 2006-8-30 10:22

哪天把自己做的习题也放上来。

chunyu 发表于 2006-8-30 13:02

这是子程序，因为最近焦头烂额，就不解释了

subroutine vumat(v.GzT2K.?)w~v0h
C Read only -k'?Vnt{
   1  nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
   2  stepTime, totalTime, dt, cmname, coordMp, charLength,
   3  props, density, strainInc, relSpinInc,
   4  tempOld, stretchOld, defgradOld, fieldOld,l f8G.A$Q
   3  stressOld, stateOld, enerInternOld, enerInelasOld,H~iz;}+j-m i
   6  tempNew, stretchNew, defgradNew, fieldNew,
C Write only -
   5  stressNew, stateNew, enerInternNew, enerInelasNew )
include 'vaba_param.inc'

C
C MODEL FOR THE SHEAR DEFORMATION AND CRACK OF POLYMERS
C SHEAR MODEL: SPRING+BURGERS ELEMENT+PLASTIC ELEMENT(D-P & EYRING)
C CRACK MODEL: S11>S11_CR WHEN SKK>0Zx1CVi+hc8oP(\
C {m-L$iV t
C WRITTEN BY  ZHANG CHUNYU([email]CHUNYU@NUS.EDU.SG[/email])
C
C All arrays dimensioned by (*) are not used in this algorithmPf r;O/en

      dimension props(nprops), density(nblock),
   1  coordMp(nblock,*),5|@gT:m D&|
   2  charLength(*), strainInc(nblock,ndir+nshr),
   3  relSpinInc(*), tempOld(*),w7@-fjb|
   4  stretchOld(*), defgradOld(nblock,ndir+nshr),
   5  fieldOld(*), stressOld(nblock,ndir+nshr),
   6  stateOld(nblock,nstatev), enerInternOld(nblock),
   7  enerInelasOld(nblock), tempNew(*),
   8  stretchNew(*), defgradNew(nblock,ndir+nshr), fieldNew(*),
   9  stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
   1  enerInternNew(nblock), enerInelasNew(nblock)
C    strain components stored as state variablesr9huBas!T4l8B&K0@,U
   dimension eelas(ndir+nshr),eplas(ndir+nshr),deplas(ndir+nshr)
dimension evisco(ndir+nshr),devisco(ndir+nshr),veint(ndir+nshr)JJ!|9rc*K(Tq,H
dimension effstrn(ndir+nshr),ps(ndir),an(ndir,ndir),s(ndir+nshr)
   dimension dr(3,3),f_old(6),f_new(6) ~+v8B-O,y

   data newton,toler/10,1.0E-6/
CV/Y#av2c8pY\
   character*80 cmname
C0B:`-Z Nr5a
   parameter( zero = 0., one = 1., two = 2., three = 3.,hNJ~&lOl @
   1  third = one/three, half = .5, twoThirds = two/three,W3Ckp!E)Rg
   2  threeHalfs = 1.5 )6ro*S*T7y
C -----------------------------------------------------------
C    PROPS(1) - EfAR*fzbk&DH
C    PROPS(2) - NU
C    PROPS(3) - E1P(@nqu}.x3i3l2A
C    PROPS(4) - ETA1TDzany5?C9K
C    PROPS(5) - ETA0
C    PROPS(6) - PHAI
C    PROPS(7) - PHAI2 FOR NON-ASSOCIATIVE FLOW
C    PROPS(8) - B    FOR EYRING MODEL
C    PROPS(9) - C1X7d Lb@+h9W$J
C    PROPS(10)- C2
C    PROPS(11)- M\ f{ {&PX}&q/J|
C    PROPS(12)- SCRAZE
C    PROPS(13)- ZETA02hg"Q,b x)b*?
C    PROPS(14)- CRAZE STRAIN LIMIT
C    PROPS(15)- SHEAR STRAIN LIMIT
C    PROPS(16)- Y0 :w t1MmZ.jg9s
C    PROPS(17)... PLASTIC CURVE
C    CALLS AHARD FOR CURVE OF SYIELD VS. PEEQ
C -----------------------------------------------------------{"m]n6t h Z}L+IS
Ced'O?a"W
   if (ndir .NE. 3) then
      write(6,1)RG q0PA G.p u
1    format(//,30X,'***ERROR - THIS VUMAT MAY ONLY BE USED FOR ',
   1       'ELEMENTS WITH THREE DIRECT STRESS COMPONENTS')%xR*M:Ke']
   endif`JV)NP"]UaT5T
C
C    ELASTIC PROPERTIESsX-a-{#N r
   ntens=ndir+nshr

   bniu=props(2)
   if(bniu.GT.0.4999.AND.bniu.LT.0.5001) bniu=0.499

C PARAMETERS FOR THE KELVIN UNIT "t"k/}9c#CG)z*]
7{f5m6n,qko:@-w4G2}
   phai=props(6)$^;O2a.[2IYw P
phai2=props(7)s)xX1elie`8Z
B=props(8)
c1=props(9)S8wjq cerz"N@
c2=props(10)
c3=three*bniu/(one+bniu)3l-_]9?s(fh
bm=props(11)!sbg[5cU,QoXO
scraze=props(12) B"Eh*ai^{Y}
zeta0=props(13) v'Ji Y{,J
craze_cr=props(14)1jl~/G9R O ai6X;v
shear_cr=props(15)
c loop the material block
jz9O3i?SP%V5p o

do 350 i = 1,nblock
#p ]!P\{5G0X5@z
c_factor=stateOld(i,4*ntens+4))H4P}$eV(\JxMm
if(c_factor .EQ. 0)  c_factor=1.0
   q_factor=c_factor
AVz \\ eFK
e0=props(1),H/g#Re9_,OK
e1=props(3)9M"Q @Q8OJ2tg
eta1=props(4)iL@k#DCGW
eta0=props(5)*}y5VrFG:c9ot8l
7`7Rc&miQKt
   e0=e0*q_factor
   bk=e0/(one-two*bniu)/threey%aJw6fgr9J
   g=e0/(one+bniu)/two
e1=e1*q_factorVCjt]6u#KD(T
g1=e1/(2.0*(1.0+bniu))
eta1=eta1*q_factorV)WO Wi~%W`0w
eta0=eta0*q_factor9r2EX$t/lE+_-\Z h6m
coef1=1.0+eta1/eta0+1.0/4.0*g1/eta0*dt?"[UM \htq
   coef2=g1+2.0*eta1/dt
coef3=1.0/4.0*g1/eta0*dtN&@4v6D.]
coef=1.0+g*coef1/coef2}H8f0Ro7@
C RECOVER AND ROTATE DEVIATORIC SHEAR STRAIN & TOTAL STRESS8uT-?xa
c    CALL ROTSIG(STATEV(1),DROT,EELAS,2,NDI,NSHR)
c    CALL ROTSIG(STATEV(NTENS+1),DROT,EPLAS,2,NDI,NSHR)9v5F#~di:R4S
c    CALL ROTSIG(STATEV(2*NTENS+3),DROT,EVISCO,2,NDI,NSHR)
   do 10 j=1,ntens(`&nonaq:SqKm1[
   eelas(j)=stateOld(i,j)9~|$V M8qe(cb:d*d
   eplas(j)=stateOld(i,j+ntens)
   evisco(j)=stateOld(i,j+3+2*ntens)g*S!`9w"_l'tW
   veint(j)=stateOld(i,j+3+3*ntens)
   deplas(j)=0.0'D_;{a*F-B[
   f_old(j)=defgradOld(i,j)XA$Z,UY5vC
   f_new(j)=defgradNew(i,j)g?,NL Z Q{
10 continue
c get rotation tensor and rotate old strain tensor
   call getdr(f_old,f_new,dr)
   call ROTTENSOR(eelas,dr,3,3)*Ib%Np l eJ@^
   call ROTTENSOR(eplas,dr,3,3){;{e A7\d)_z
call ROTTENSOR(evisco,dr,3,3)
   eqplas_shr=stateOld(i,1+2*ntens)
eqplas_crz=stateOld(i,2+2*ntens)-Qj]@@K^9dy
      vplas=stateOld(i,3+2*ntens)

C CALCULATE DEVIATORIC STRESS4f D X]CE'tOo
   oldstresskk=0.0Jv+g Pp YJ0z5S L ] J
dstrankk=0.0
strankk=0.09oY7V|+r,q:x
do 20 j=1,ndirF#Ra } c,}6J
   oldstresskk=oldstresskk+stressOld(i,j)
   dstrankk=dstrankk+strainInc(i,j)
      strankk=strankk+eelas(j)+eplas(j)+evisco(j)+strainInc(i,j)"}n:v^m.Q
20 continue
C CALCULATE EFFECTIVE STRAIN
   do 30 j=1,ndir
   effstrn(j)=eelas(j)+(strainInc(i,j)-dstrankk/3.0)-
1    (0.5*(coef1+2.0*coef3)*(stressOld(i,j)-oldstresskk/3.0)+
1    g1/eta0*veint(j)-2.0*g1*evisco(j))/coef2
30 continue#q!ej^!}C/Y
   do 31 j=ndir+1,ntens0v@Eu'P)s&z
   effstrn(j)=eelas(j)+strainInc(i,j)/2.0-sx|(Aj3A8wB@4l
1             (0.5*(coef1+2.0*coef3)*stressOld(i,j)+
1             g1/eta0*veint(j)-2.0*g1*evisco(j))/coef2"UAX/R+D\,EG c(V
31 continue7TX2g8_:u%`R;N2\
C UPDATED DEVIATORIC STRESS
   do 40 j=1,ntens N#J8Qiq~
   stressNew(i,j)=2.0*g/coef*effstrn(j)
40 continue d}j@3Kef
C CALCULATE BULK STRESS & PRESSURE'y.x,r3r'Z9N
   stresskk=oldstresskk+3.0*bk*dstrankk2CZ{l5J:\PuQ
p=-stresskk/3.0 +Q'd6z{AF
C UPDATE DIRECT STRESSo+J*\%^%v bf&XqW
   do 45 j=1,ndira7Hi-n'F1Z#~bIx
stressNew(i,j)=stressNew(i,j)-p)@~:i4hYR o
45 continue
C DETERMINE SHEAR YILED OR CRAZE YIELDf0A3j'`s6EY`Y
   do 50 j=1,ntens
   s(j)=stressNew(i,j)yFEY~!?+v
50 continue
C PRINCIPAL STRESS AND ,wMFcMx?
   call PRINCIPAL(s,ps,an,ndir,nshr,1.0E-6)
ps1=ps(1)
nmax=1
if(ps(2) .GT. ps1) then
   ps1=ps(2))e^p4@'PR.E
   nmax=2Q;j)Y&GZ
endif
if(ps(3) .GT. ps1) then
   ps1=ps(3)
   nmax=3K*EUk6F:J
endif
bms=-pi2L.Lz2CQ!g
cs=c1+c2/bms+c3*bms RBv`fg
if (bms .GT. 0 .AND. ps1 .GT. 0 .AND. ps1 .GT. cs) then
      flow=oneUEs6F:EX[axT?
else @;MQNu\4S
      flow=zeroKu6ga7H`#L-r|5Xb
endif
8V~+c[Qd1y.XeVG
C CRAZE FLOW,calculate deplas
   if (flow .EQ. one .AND. c_factor .GT. 1.0E-4) then*W _C#T9B9Nx
   do 60 j=1,ndir
         deplas(j)=zeta0*(ps1/scraze/q_factor)**(1.0/bm)*dt*EG$T"fC#i l{
1             an(j,nmax)**2
60    continue&F;A3{(P3T^K)B
      do 61 j=ndir+1,ntens
         n1=j-ndir:m Jf:D+nnj?
         n2=j-ndir+1
         if(n2 .GT. ndir) thenzb/\8y#^!l2}"e@xv
   n2=n2-ndir
      endif .j\N/f$amI g^2_
         deplas(j)=zeta0*(ps1/scraze/q_factor)**(1.0/bm)*dt*
1             an(n1,nmax)*an(n2,nmax)*A1E8gaPb4n
61    continue
   eqplas_crz=eqplas_crz+zeta0*(ps1/scraze/q_factor)**(1.0/bm)*dtz+N#v.B!cL8R`
if(eqplas_crz .GT. craze_cr .AND. c_factor .GT. zero) then1@f3U i`}IX
c_factor=c_factor*0.1*craze_cr/eqplas_crz8P"\v,K7l/Wq
endif
dU*MI$P)oJtZ*vsv s
c    update stress

   dvplaskk=deplas(1)+deplas(2)+deplas(3)(b?p'oB%?+zrh
   vplas=vplas+dvplaskk
   stresskk=oldstresskk+3.0*bk/q_factor*c_factor*_+H"fR~._
1          (dstrankk-dvplaskk)fI'w6B/yZVw
      do 70 j=1,ndir
      deplas(j)=deplas(j)-dvplaskk/3.0.ClL5{j
         stressNew(i,j)=2.0*g/q_factor*c_factor*+n I)DXp n5x
1                   (effstrn(j)-deplas(j))/coef
         stressNew(i,j)=stressNew(i,j)+stresskk/3.0V/SvWI1u0gj
70    continue
      do 71 j=ndir+1,ntens
         stressNew(i,j)=2.0*g/q_factor*c_factor*v*|2d |mD c
1                   (effstrn(j)-deplas(j))/coefH#Le9DyP
71    continue    pExV7{t
endif
V"S-g-S!G6s$GfB^
C SHEAR YIELD FLOW,CALCULATE DEPLAS
c test crack-->shear yielding
c if crack and shear yielding can not coexist,delete the following the command
c    flow=zero
   if(nprops.GT.15 .AND. props(16).GT. 0.0 .AND. flow .EQ. zero) then,]W[:{vME'ra
C
C    MISES STRESSQ_[C9p
Cs A8F"K!M]I/[d+TW&\h
      smises=(s(1)-s(2))*(s(1)-s(2)) +
   1       (s(2)-s(3))*(s(2)-s(3)) +
   1       (s(3)-s(1))*(s(3)-s(1))
      do 90 j=ndir+1,ntensP4C$d\M l
            smises=smises+6.0*s(j)*s(j)]8`U_@2rH.k YM
90    CONTINUE
      smises=sqrt(smises/2.0)ym"D3w#RL;g{
C    SHEAR STRAIN RATE.W(b|q!H uiN
      effrate=0.0
   do 95 j=1,ndir
      effrate=effrate+(strainInc(i,j)-dstrankk/3.0)*rm6y"kj5SB
   1                   (strainInc(i,j)-dstrankk/3.0)
95    continue
      do 100 j=ndir+1,ntens
      effrate=effrate+2.0*strainInc(i,j)*strainInc(i,j)5B3SFcz#}0a^6R8CN
100    CONTINUE
   effrate=sqrt(effrate/2.0)/dt
   if (effrate .LT. 1.0E-6) then/WkN/Go7X!v
         effrate=1.0E-6,^1g.O:_%m
   endif
C    HARDENING CURVE, GET YIELD STRESS
C    CALCULATE HOW MANY POINTS FOR THE PLASTIC CURVE
C       COEF5  DEPQ=COEF5*LAMTA.vM3\e@
   coef5=sqrt(1.0+phai2*phai2/2.0)I5|5rI U2Z-R
      NVALUE=nprops/2-7
      CALL AHARD(syiel0,HARD,eqplas_shr*coef5,props(16),NVALUE)
C    DETERMINE IF ACTIVELY YIELDING
C
   effyield=syiel0+B*log10(effrate)+phai*p
   if (smises .GT. (1+toler)*effyield ) then?P:G"B.t3@P
C    CALCULATE EQUIVALENT STRAIN 8C8m"wK3u-i#oX
      equie=0.0  a s @#w;Eb R2W
      do 110 j=1,ndirH |^ N](bJg
         equie=equie+effstrn(j)*effstrn(j)
110    continueEu$QA&ZcWFs
      do 111 j=ndir+1,ntens2NX} q4x }O9D
         equie=equie+2*effstrn(j)*effstrn(j)
111    continueD~d.dK#Q\
      equie=sqrt(2.0/3.0*equie)^[5}l&h7S3_f
C
C    SOLVE FOR EQUIV STRESS, NEWTON ITERATIONz'Q"n}Je
C    Maybe unnecessary for explicit solver
      syield=syiel0z@;pp:T9@!uVh
      deqpl=0.0I1~KA._F%Q
      do 130 j=1,newtonmV.B3U2D.}"P4jS
         rhs=3.0*g*(equie-deqpl)-coef*(syield+B*log10(effrate)
1       -phai*bk*(strankk-vplas-phai2*deqpl))
         deqpl=deqpl+rhs/(3.0*g+coef*(HARD*coef5+bk*phai*phai2))
   call AHARD(syield,HARD,(eqplas_shr+deqpl)*coef5,props(16),NVALUE)
      if(abs(rhs).LT.toler*effyield) goto 140
130    continue
      write(6,2) newton
2       format(//,30X,'***WARNING - PLASTICITY ALGORITHM DID NOT ',
   1       'CONVERGE AFTER ',I3,' ITERATIONS')'T'u|+F'r
140    continueu!PMtZwH

      stresskk=3.0*bk*(strankk-vplas-phai2*deqpl)
p=-stresskk/3.0:Hhx!gZ
      effyield=syield+phai*p+B*log10(effrate)$}U$Pe@I$[
coef4=3.0*g*deqpl/effyield+coef?W~0e({rL^Nl
      do 150 j=1,ndir
         stressNew(i,j)=2.0*g*effstrn(j)/coef4
         deplas(j)=3.0*stressNew(i,j)*deqpl/(2.0*effyield)
         stressNew(i,j)=stressNew(i,j)+stresskk/3.05y,or~O b/f
150    continuev5d4B8c Q!~
      do 160 j=ndir+1,ntens
         stressNew(i,j)=2.0*g*effstrn(j)/coef4b.rI5zl^
         deplas(j)=3.0*stressNew(i,j)*deqpl/(2.0*effyield)z&wF_B3VhW8g
160    continue)\Xn!m/j {t
      eqplas_shr=eqplas_shr+deqpl%f:Y\ ]9\7\
      vplas=vplas+phai2*deqpl qG"N[bI4}
if(eqplas_shr .GT. shear_cr .AND. c_factor .GT. zero) theneX@/u!x\Z3}
      c_factor=c_factor*0.1*shear_cr/eqplas_shr
      endif
      endif
endif-j|2W n\Tr h1m+N(ru@
C UPDATE STATE VARIABLES
C CALCULATE DEVISCO
   do 260 j=1,ndirR@0iw8s/[5}Ff
   devisco(j)=(0.5*coef1*(stressOld(i,j)-oldstresskk/3.0+
1 stressNew(i,j)-stresskk/3.0)+coef3*(stressOld(i,j)-r]Mn^||*O
   1 oldstresskk/3.0)+g1/eta0*veint(j)-2*g1*evisco(j))/coef2
260 continuecKu4AN8a3Z}-_
   do 265 j=ndir+1,ntens3_7oM.J9X.mtI
   devisco(j)=(0.5*coef1*(stressOld(i,j)+stressNew(i,j))+
1             coef3*stressOld(i,j)+g1/eta0*veint(j)-
1    2*g1*evisco(j))/coef2
265 continue
C STORE STATE VARIABLE ARRAY
do 300 j=1,ndir
   stateNew(i,j)=eelas(j)+strainInc(i,j)-dstrankk/3.0-devisco(j)-"h{;_*Hw$l,ox~
   1             deplas(j)
   stateNew(i,j+ntens)=eplas(j)+deplas(j)
   stateNew(i,j+3+2*ntens)=evisco(j)+devisco(j)m-|7l})o;h#o]
      stateNew(i,j+3+3*ntens)=veint(j)+0.5*dt*(stressOld(i,j)-
1          oldstresskk/3.0+stressNew(i,j)-stresskk/3.0)FR0m M*o
300 continue
do 310 j=ndir+1,ntens
      stateNew(i,j)=eelas(j)+strainInc(i,j)/2.0-devisco(j)-deplas(j)(Ao;j{Ik#[
   stateNew(i,j+ntens)=eplas(j)+deplas(j)
   stateNew(i,j+3+2*ntens)=evisco(j)+devisco(j).M%GRg!`AUW
      stateNew(i,j+3+3*ntens)=veint(j)+0.5*dt*(stressOld(i,j)+`6X`0vT0y&yS
1                                     +stressNew(i,j))
310  continueHq$@Fh"`(HE
   stateNew(i,1+2*ntens)=eqplas_shr
   stateNew(i,2+2*ntens)=eqplas_crz
   stateNew(i,3+2*ntens)=vplas
stateNew(i,4*ntens+4)=c_factor
if(c_factor .LT. 1.0E-4) then
   stateNew(i,4*ntens+5)=zero'Kbfz-^ R+M7}$?
else3D2w2FB#F:ri[
      stateNew(i,4*ntens+5)=oneEq-F;{O%@(mn
endif
350  continue_!Y;w @BH
C5I1O.xlD0XC%OF
   return
   end7M$mn/}%J;O"u'Ub)P

C
c    CALCULATE EQUIVALENT STRESS4_DwhfD
   SUBROUTINE AHARD(SYIELD,HARD,EQPLAS,TABLE,NVALUE)
C
   INCLUDE 'ABA_PARAM.INC'
   DIMENSION TABLE(2,NVALUE)
C {"M0_zI"Bj\
C SET YIELD STRESS TO LAST VALUE OF TABLE, HARDENING TO ZERO
   SYIELD=TABLE(1,NVALUE)R]hg%DL^
   HARD=0.0
CYb,[(S a#?(\
C IF MORE THAN ONE ENTRY, SEARCH TABLE
CG$mBHyY+r6N4e(W a
   IF(NVALUE.GT.1) THEN
      DO 10 K1=1,NVALUE-1
         EQPL1=TABLE(2,K1+1)3| I@YB};EQ
         IF(EQPLAS.LT.EQPL1) THEND1J3{f)OCE1e
         EQPL0=TABLE(2,K1)W0z\s8`8g1o
         IF(EQPL1.LE.EQPL0) THEN
            WRITE(6,1)
1             FORMAT(//,30X,'***ERROR - PLASTIC STRAIN MUST BE ',wW7TgP6r&r
   1                'ENTERED IN ASCENDING ORDER')
C             CALL XIT0|8o*SOb5@y
            ENDIF
C
C          CURRENT YIELD STRESS AND HARDENINGQh@a%~D
C
         DEQPL=EQPL1-EQPL0
         SYIEL0=TABLE(1,K1)
         SYIEL1=TABLE(1,K1+1)
         DSYIEL=SYIEL1-SYIEL0
         HARD=DSYIEL/DEQPL
         SYIELD=SYIEL0+(EQPLAS-EQPL0)*HARD
         GOTO 20Mhi9}2caS.FEG
         ENDIF
10       CONTINUE)i#b5}jkq['O
20       CONTINUE
   ENDIF+K~r1Y%qz2D0c
   RETURN
   END!IOb(T2?6l
C ROUNTINE TO CALCULATE EIGENVALUES AND EIGENVECTORS
   SUBROUTINE PRINCIPAL(S,PS,V,NDIR,NSHR,TL)b0z?Qx
C    IMPLICIT REAL*8 (A-H,O-Z)    uQ{'YH"p&Z
   Y*X#TW6D5_?2O
   DIMENSION A(NDIR,NDIR),V(NDIR,NDIR),S(NDIR+NSHR),PS(NDIR)CL`#G DS|.N

   NEQ=NDIR
C    EIGENVALUE SOLUTION BY JACOBI METHOD -
C    A - MATRIX (ANY RANK) TO BE SOLVED ---:q/P$M\&f1|N}i
C       EIGENVALUES ON DIAGONAL'A)Swf3~~2m S
C    V - MATRIX OF EIGENVECTORS PRODUCEDb^(l FcSXyStm
C    TL- NUMBER OF SIGNIFICANT FIGURES W7K-m:w{(\s
C---- INITIALIZATION -----------------------
   ZERO = 0.0D0T`4V+E0mHg
   SUM = ZERO
   TOL = ABS(TL).M^!U{}n
A(1,1)=S(1)
A(2,2)=S(2)
A(3,3)=S(3):L7`6Xn&I-~h
   A(1,2)=S(NDIR+1)vz~F6iop,uc
A(2,1)=S(NDIR+1)
IF(NSHR .GT. 1) THEN
A(1,3)=S(NDIR+3)
A(3,1)=S(NDIR+3)
A(2,3)=S(NDIR+2)kT@,V0fih'|+_ v
A(3,2)=S(NDIR+2)
ENDIFCy v\I.Y~
C---- SET INITIAL EIGENVECTORS -------------.RkI DIjl?%w
   DO 200 I=1,NEQ$k6Z$PEO#W@!fX
   DO 190 J=1,NEQ
      IF (TL.GT.ZERO) V(I,J) = ZERO
190    SUM = SUM + ABS(A(I,J))+y(_'^#Ml
   IF (TL.GT.ZERO) V(I,I) = 1.0 }(H^f.@ qa
200  CONTINUE:D,I X"pn"D%O3h
C---- CHECK FOR TRIVIAL PROBLEM -----------5t7x7U6}y0n;~;o
   IF (NEQ.EQ.1) RETURN0ERIN vLZ'J
      IF (SUM.LE.ZERO) RETURNp#A7v } p7@_
      SUM = SUM/FLOAT(NEQ*NEQ),V8I8Q2nc%qg$t$R8Z!c-w
C-------------------------------------------\CC3Lpit
C---- REDUCE MATRIX TO DIAGONAL ------------2p ZoUh I+j
C-------------------------------------------
400  SSUM = ZERO-[i5H:tsj J [G-X
   AMAX = ZERO-tZ2s7qcj*y
   DO 700 J=2,NEQ
   IH = J-1
   DO 700 I=1,IH:oc1Q,U"BPv)\"oO
C---- CHECK IF A(I,J) IS TO BE REDUCED ----- p3aT"TS:d~e+R
   AA = ABS(A(I,J))b6s FH0IVs
   IF (AA.GT.AMAX) AMAX = AAxg/J#LaR)C
   SSUM = SSUM + AAmu D+c~m%B6b
   IF (AA.LT.0.1*AMAX) GO TO 700
C---- CALCULATE ROTATION ANGLE ----------
   AA=ATAN2(2.0*A(I,J),A(I,I)-A(J,J))/2.04WDK$LeVz
   SI = SIN(AA)
   CO = COS(AA)
C---- MODIFY "I" AND "J" COLUMNS --------hhP/k o#a$l3[o
   DO 500 K=1,NEQIb.F{#Z I(\s&jb
   TT = A(K,I)e8zO T9U*Ko
   A(K,I) = CO*TT + SI*A(K,J)$jAyFLz
   A(K,J) = -SI*TT + CO*A(K,J)
   TT = V(K,I)?B1s!]4T
   V(K,I) = CO*TT + SI*V(K,J)2XU;xtz"i3io9{
500 V(K,J) = -SI*TT + CO*V(K,J)YY&t+o]0Be"T$o
C---- MODIFY DIAGONAL TERMS -------------
   A(I,I) = CO*A(I,I) + SI*A(J,I)
   A(J,J) =-SI*A(I,J) + CO*A(J,J)
   A(I,J) = ZERO
C---- MAKE "A" MATRIX SYMMETRICAL -------/m6lq%l;KK/i
   DO 600 K=1,NEQ]a;Um1]#v(_'\ ~
   A(I,K) = A(K,I)%Ooj#K n7XBd
   A(J,K) = A(K,J)Om5]~/R}%K
600 CONTINUE:`K#N&^#N2NjS]Y
C---- A(I,J) MADE ZERO BY ROTATION ------0YeY8wJ \-l9X#X L
700 CONTINUE5]iD8t?tzX
C---- CHECK FOR CONVERGENCE -------------o,`0c+u0p)pDn
   IF(ABS(SSUM)/SUM .GT.TOL)GO TO 400v_ X$@\ @ YP_
DO 900 I=1,NDIR&yK&Eh W{0}?i7D9V
   PS(I)=A(I,I)j.|}Z([DL;HUY
900 CONTINUE
   RETURN!I't~3q~ BI-W
   END2^*gA(MV2x%{]
C****************************************************************************
   SUBROUTINE ROTTENSOR(TENSOR,DR,NDIR,NSHR)
C    ROTATE A TENSOR:  TENSOR=DR.TENSOR.trans(DR)x l#GV}(|DO
C    TENSOR IS STORED IN A VETOR(NDIR+NSHR)
C****************************************************************************,E(R%wf%g0H:q
c    IMPLICIT REAL*8(A-H,O-Z)
   DIMENSION TENSOR(NDIR+NSHR),DR(NDIR,NDIR)(S:UI\6E`
DIMENSION STENSOR(3,3),TEMP(3,3),TDR(3,3)#m2_{1u;r.}8yc8N
   CALL TRANS(DR,TDR)@2Va5x~
   CALL GETU(TENSOR,STENSOR,NDIR,NSHR)
CALL MUL(DR,STENSOR,TEMP,NDIR)^(mkHI_\]5Hh-|
CALL MUL(TEMP,TDR,STENSOR,NDIR)
C    REARRANGE AND STORE IN A VECTOR
   TENSOR(1)=STENSOR(1,1)
TENSOR(2)=STENSOR(2,2)
TENSOR(3)=STENSOR(3,3)!l@c5Bh
TENSOR(4)=STENSOR(1,2)PHU/G&^p8^
TENSOR(5)=STENSOR(2,3)9V2?kvQ` Ma/q;|
TENSOR(6)=STENSOR(1,3)5J1yJsg
'q-} ZE5ro3K
   RETURN
END

C****************************************************************************
   SUBROUTINE GETDR(FO,FN,DR)8vUnx,{\:U2|
C    CALCULATE DELTA_R  DF=DR.DU
C    DU=TRANSPOSE(DU) DR.TRANSPOSE(DR)=I4WKr~ s)\Dio
C****************************************************************************
c    IMPLICIT REAL*8(A-H,O-Z)#kIa&tL|/c-\
DIMENSION FOLD(6),FNEW(6),DR(3,3)$WV'K.tb?\{;O
   DIMENSION FO(3,3),FN(3,3),FOINV(3,3),DF(3,3),DC(3,3),TRANSDF(3,3)
DIMENSION V(3,3),DU(3,3),DUINV(3,3)2ER*cH0r&Re}g
CALL GETF(FOLD,FO,3,3)c8o$Q7mlu0j
CALL GETF(FNEW,FN,3,3)5RTL$My#c%P qI
C    FOLD^-1+Gx@BJ@JY(y:]2p#[
   CALL KMINV(FO,FOINV,DET)
C    DF=FNEW.FOLD^-1#dL!sJK5t
   CALL MUL(FN,FOINV,DF,3)T,X:m[?tG"h
C    TRANSPOSE(DF)QH;vb6QU*]
   CALL TRANS(DF,TRANSDF)
C    CALCULATE TRANSPOSE(DF).DFdT-@{&dF
   CALL MUL(TRANSDF,DF,DC,3)^BY1_v*qI(~j
C    CALCULATE THE EIGENVALUES AND EIGENVECTORS OF DC
   CALL EIGEN(DC,V,3,1.0e-6)7qLV{"u'^ h@ b
C    CALCULATE DU Z@3lV,E'{#M
   DO 20 I=1,3
   DO 10 J=1,3 Zp,WI"v7lZo.B
      DU(I,J)=SQRT(DC(1,1))*V(I,1)*V(J,1)+6Ensf&CCp.h#a i
1          SQRT(DC(2,2))*V(I,2)*V(J,2)+-QP1FWVT,Q
1          SQRT(DC(3,3))*V(I,3)*V(J,3)T$z2EX{
10    CONTINUE
20 CONTINUE
C    CLACULATE DU^-1
   CALL KMINV(DU,DUINV,DET)
C    DR=DF.DU^-1f4Na4ams4s QZi
CALL MUL(DF,DUINV,DR,3)
RETURN 0J5W P2b:j,}]2i
END.i;f K#Ez{f7P
C****************************************************************************
   SUBROUTINE EIGEN(A,V,NEQ,TL){%[|*hG8b
c    IMPLICIT REAL*8 (A-H,O-Z)'\ r^XM;`
   DIMENSION A(NEQ,NEQ),V(NEQ,NEQ)W&@#kd+K:mK
C    EIGENVALUE SOLUTION BY JACOBI METHOD -
C    A --MATRIX (ANY RANK) TO BE SOLVED ---BZ2F;lw&tR5~
C       EIGENVALUES ON DIAGONAL_b:k A5N2^O9@)Q
C    V - MATRIX OF EIGENVECTORS PRODUCEDZ%HJ4x)b:H
C    TL- NUMBER OF SIGNIFICANT FIGURES
C---- INITIALIZATION -----------------------2L`l2S!{mK}
   ZERO = 0.0D0
   SUM = ZERO_*}5|Z&k;F'p7o2[xS
   TOL = ABS(TL)5Q(X7I!fyaa
C---- SET INITIAL EIGENVECTORS --------------`:g4u0t?9F#?
   DO 200 I=1,NEQi d-T%` zT
   DO 190 J=1,NEQ
   IF (TL.GT.ZERO) V(I,J) = ZERO
190 SUM = SUM + ABS(A(I,J))
   IF (TL.GT.ZERO) V(I,I) = 1.0
200 CONTINUE
C---- CHECK FOR TRIVIAL PROBLEM -----------'cc@6d5Q]
   IF (NEQ.EQ.1) RETURN
   IF (SUM.LE.ZERO) RETURN
   SUM = SUM/FLOAT(NEQ*NEQ) Lz7lx`GJn
C-------------------------------------------
C---- REDUCE MATRIX TO DIAGONAL ------------
C-------------------------------------------Rn)U4? _%Ue-Y
400 SSUM = ZERO
   AMAX = ZERO
   DO 700 J=2,NEQT*N8W'E!Oyo
   IH = J -1
   DO 700 I=1,IH)BaOhg4?V
C---- CHECK IF A(I,J) IS TO BE REDUCED -----
   AA = ABS(A(I,J))"tj*Yo G
   IF (AA.GT.AMAX) AMAX = AAP;Fk'y3pU!}M7G&x
   SSUM = SSUM + AA`%hQ3H3G9S
   IF (AA.LT.0.1*AMAX) GO TO 700)aR JG/p7X
C---- CALCULATE ROTATION ANGLE ----------$a0@mI;[
   AA=ATAN2(2.0*A(I,J),A(I,I)-A(J,J))/2.0.QqNcD5p|0F
   SI = SIN(AA)(p/~/Pg-r i0t d
   CO = COS(AA)
C---- MODIFY "I" AND "J" COLUMNS --------VE9h*V0G;} @
   DO 500 K=1,NEQ
      TT = A(K,I)
      A(K,I) = CO*TT + SI*A(K,J)
      A(K,J) = -SI*TT + CO*A(K,J)MiC#Pd'v
      TT = V(K,I) etbP pKf-}
      V(K,I) = CO*TT + SI*V(K,J) kv F G~/s5@
500    V(K,J) = -SI*TT + CO*V(K,J)F0?&{5h K
C---- MODIFY DIAGONAL TERMS -------------
      A(I,I) = CO*A(I,I) + SI*A(J,I)
      A(J,J) =-SI*A(I,J) + CO*A(J,J)3@NKjR VW
      A(I,J) = ZERO
C---- MAKE "A" MATRIX SYMMETRICAL -------
   DO 600 K=1,NEQ
      A(I,K) = A(K,I)
      A(J,K) = A(K,J)
600 CONTINUE%x1A ^ q3D(H8_
C---- A(I,J) MADE ZERO BY ROTATION ------(e'l:^&ts%g4[B
700 CONTINUEU$k1M!aU}+r
C---- CHECK FOR CONVERGENCE -------------,m4^R7IFh)m3Xt
   IF(ABS(SSUM)/SUM .GT.TOL)GO TO 400"]:veg6i
,n od7` Zd%Q
   RETURN
   END
-iZ[Y:J c
C****************************************************************************d(d/G1^Sl
xK^w(r \/BY
C****************************************************************************||g4i*G-R e
SUBROUTINE KMINV(A,AINV,DET_AINV)ESK+w4K:s
Li#@P4e
C This subroutine calculates the inverse of a {3 x 3} matrix and the
C determinant of the inverse+kJ[}8k$P
C**************************************************************************** +h'mSx9Uda
c    IMPLICIT REAL*8(A-H,O-Z)
A#i+f_v!aU#{R
DIMENSION A(3,3), AINV(3,3)

PARAMETER(ZERO=0.D0, ONE=1.D0)9T,Xy aixG @(E&U
E^$v)k#vj0z
DET_A = A(1,1)*(A(2,2)*A(3,3) - A(3,2)*A(2,3)) -I1MnWIXb2g
   +          A(2,1)*(A(1,2)*A(3,3) - A(3,2)*A(1,3)) +C!m;IcdU\n(p
   +          A(3,1)*(A(1,2)*A(2,3) - A(2,2)*A(1,3))
$i*G2SSHJ/^7SI%\6`
IF (DET_A .LE. ZERO) THEN
   WRITE(80,*) 'WARNING: DET OF MAT IS ZERO/NEGATIVE !!';Y:vSAh3^2a3W7Iy
ENDIF
#iV,_1w)mg
DET_AINV = ONE/DET_Aa%[F3}}4MK
#f*_`q8vf$V6X,dU
AINV(1,1) = DET_AINV*(A(2,2)*A(3,3) - A(3,2)*A(2,3))5ju#lD,fjw
AINV(1,2) = DET_AINV*(A(3,2)*A(1,3) - A(1,2)*A(3,3))
AINV(1,3) = DET_AINV*(A(1,2)*A(2,3) - A(2,2)*A(1,3))
AINV(2,1) = DET_AINV*(A(3,1)*A(2,3) - A(2,1)*A(3,3))
AINV(2,2) = DET_AINV*(A(1,1)*A(3,3) - A(3,1)*A(1,3))
AINV(2,3) = DET_AINV*(A(2,1)*A(1,3) - A(1,1)*A(2,3))7dg:MyE$].c f
AINV(3,1) = DET_AINV*(A(2,1)*A(3,2) - A(3,1)*A(2,2))
AINV(3,2) = DET_AINV*(A(3,1)*A(1,2) - A(1,1)*A(3,2))
AINV(3,3) = DET_AINV*(A(1,1)*A(2,2) - A(2,1)*A(1,2))

RETURNR3Nb|(gI e%f _
END_-NX[v6X^
C*******************************************************************
SUBROUTINE MUL(A,B,C,N)%r%U.~#[$A@
C MULTIPLICATION OF SQUARE MATRIX A,B. THE RSULTE IS STORED IN CtEZ/j;]
C N IS THE DIMENSION OF THE MATRIX;xRKv2|-\&d)qu
C*******************************************************************
c    IMPLICIT REAL*8(A-H,O-Z)o@l GG.Yac,m)^ d
   DIMENSION A(N,N),B(N,N),C(N,N)
DO 30 I=1,N#g)y#]uSB-`3S
   DO 20 J=1,N
      C(I,J)=0 T hq0IKZ
      DO 10 K=1,N\"P0A7koG?C
         C(I,J)=C(I,J)+A(I,K)*B(K,J))S.Z&Au/o
10       CONTINUE@0KR6oj,~
20 CONTINUEJ:x+vzZ#Ov"Z
30  CONTINUE u qUYL
$TBRF%w1H)_
   RETURN
END
C**********************************************************************!`lH3}.Q I)X8\$w"`
SUBROUTINE TRANS(A,ATRANS)6r@2Sy*p$T`

C THIS SUBROUTINE CALCULATES THE TRANSPOSE OF AN 3 BY 3 s1u9[(]#C,A {q|
C MATRIX [A], AND PLACES THE RESULT IN ATRANS.
C**********************************************************************

c REAL*8 A(3,3),ATRANS(3,3)
   dimension A(3,3),ATRANS(3,3)
DO 1 I=1,3
   DO 2 J=1,3h%oMtm@(s8T
      ATRANS(J,I) = A(I,J)F@;b']+A|F
2    CONTINUE
1  CONTINUE
RETURN
END
$_\Z*t/P0f@F2Q^
C**********************************************************************M&}1I&cKhM
X*q!d&M*Z,Y5~
C****************************************************************/P*E8X;v(Z
   SUBROUTINE GETF(F,SF,NDIR,NSHR)Id*e4jh&Q!L R
C    REARRANGE DEFORMATION GRADIENT STORED IN A ONE-DIMENSION VETOR #q1}!Jd6qty$iFr
C    TO A SQUARE MATRIXG2ND,^0?n
C    F(NDIR+2*NSHR)    SF(NIDR,NDIR)
C****************************************************************
c    IMPLICIT REAL*8(A-H,O-Z)
   DIMENSION F(NDIR+2*NSHR),SF(NDIR,NDIR)
   IF (NSHR .LT. 3) THEN
   WRITE(80,*) 'WARNING: CAN ONLY BE USED IN 3D MODEL '
ENDIF
SF(1,1)=F(1)
SF(1,2)=F(4)[ `3pNdbtY$I7|*b O
SF(1,3)=F(9)&[}8_"CW EC
SF(2,1)=F(7)
SF(2,2)=F(2)}RAD1L
SF(2,3)=F(5)
SF(3,1)=F(6)
SF(3,2)=F(8)
SF(3,3)=F(3)erkLRj zs
.^9qSfwm^
RETURN
END1M$J6yS+Y0^t2A
C****************************************************************
   SUBROUTINE GETU(U,SU,NDIR,NSHR)L!C;ckF/WG0A*`
C    REARRANGE STRETCH TENSOR STORED IN A ONE-DIMENSION VETOR
C    TO A SQUARE MATRIXIew$u^
C    U(NDIR+NSHR)    SU(NIDR,NDIR)x+Q,olNrj
C****************************************************************
c    IMPLICIT REAL*8(A-H,O-Z)Sd;j(g1D|.\
   DIMENSION U(NDIR+NSHR),SU(NDIR,NDIR)
   IF (NSHR .LT. 3) THEN*J1]*FmqB|0l+{
   WRITE(80,*) 'WARNING: CAN ONLY BE USED IN 3D MODEL '
ENDIF
SU(1,1)=U(1)
SU(1,2)=U(4)
SU(1,3)=U(6)
SU(2,1)=U(4)
SU(2,2)=U(2).sp?"uy$Z6M
SU(2,3)=U(5)
SU(3,1)=U(6)| d)vU8N~Ab
SU(3,2)=U(5)N2D Y-{!US&SJD
SU(3,3)=U(3)]jwM(lww(V"V
G}K&\)FN
RETURN
END
仅供参考，请不要另作他用。如果需要，可与本人联系。

hanqinhu 发表于 2006-9-1 13:15

真仗义，呵呵，论坛里的哥们都这样就好了，向你学习，当然了，俺是有心无力，低手一个，呵呵

luanluan 发表于 2006-9-17 10:48

谢谢了

blarry 发表于 2006-9-17 11:38

呵呵，什么时候我能看懂就好了

hby087 发表于 2006-9-18 22:07

麻烦给我发一份，谢谢了
[email]hby087@gmail.com[/email]

wxwswsw 发表于 2006-9-21 16:22

太谢谢楼主，之前看过帖子，觉得没有源程序，实在是太遗憾了，现在看到chunyu兄把vumat都贴出来了，真是感动啊！这帖子太好了，大家顶啊

lbxsnh 发表于 2006-10-25 15:15

Re:材料失效及其模拟

你好，我也正在做材料的失效及断裂，能否将您的inp和.f给我看一下呢？
[email]lbxsnh@163.com[/email]-`7y/O6j EDLcmWa
不胜感激！

wwwaba 发表于 2006-10-26 07:38

您好，我正在做材料的失效及断裂，能否将您的inp和.f给我看一下
用[email]ytdoc@163.com[/email]

yaooay 发表于 2006-11-5 07:11

你好，我也正在开始做材料的失效及断裂，能否将您的inp和.f给我看一下呢？
[email]yaooay@gmail.com[/email] ga5y~mrw:J
非常感谢！

herojoe 发表于 2006-11-5 11:07

高手何不上传一个inp文件或cae 文件呢？想参考一下，谢谢！！

herojoe 发表于 2006-11-10 10:57

哪位大哥要到了楼主的cae文件，想参考一下，谢谢！！

loveluckycs 发表于 2006-12-11 21:52

非常谢谢！！！！！！！！！！！！！！！！！！

lingke125 发表于 2007-3-21 10:59

非常感谢 !&aE m+~}4fC
可是我还是不很明白！1u{3]3] ~d
楼主能给我发一份inp文件么？LZ X"HD-X
[email]lingke125@163.com[/email]
谢谢！

drkxzw 发表于 2007-4-29 13:18

偶也要一份:blue-[email]green-eye@163.com[/email],多谢!

cip05yh 发表于 2007-4-30 02:18

我也要一份，谢谢了！'B*TR~~-A
ZB)R$g8Nvw

[email]huyings@yahoo.com.cn[/email]

页: [1] 2 3

SimWe仿真论坛's Archiver