- 积分
- 1
- 注册时间
- 2015-1-23
- 仿真币
-
- 最后登录
- 1970-1-1
|
发表于 2017-3-22 12:02:50
|
显示全部楼层
来自 江苏南京
- subroutine vumat(
- C Read only (unmodifiable)variables -
- 1 nblock, ndir, nshr, nstatev, nfieldv, nprops, lanneal,
- 2 stepTime, totalTime, dt, cmname, coordMp, charLength,
- 3 props, density, strainInc, relSpinInc,
- 4 tempOld, strainOld, defgradOld, fieldOld,
- 5 stressOld, stateOld, enerInternOld, enerInelasOld,
- 6 tempNew, stretchNew, defgradNew, fieldNew,
- C Write only (modifiable) variables -
- 7 stressNew, stateNew, enerInternNew, enerInelasNew )
- C
- include 'vaba_param.inc'
- C
- dimension props(nprops), density(nblock), coordMp(nblock,*),
- 1 charLength(*), strainInc(nblock,ndir+nshr),
- 2 relSpinInc(nblock,nshr), tempOld(nblock),
- 3 strainOld(nblock,ndir+nshr),
- 4 defgradOld(nblock,ndir+nshr+nshr),
- 5 fieldOld(nblock,nfieldv), stressOld(nblock,ndir+nshr),
- 6 stateOld(nblock,nstatev), enerInternOld(nblock),
- 7 enerInelasOld(nblock), tempNew(*),
- 8 strainNew(nblock,ndir+nshr),
- 8 defgradNew(nblock,ndir+nshr+nshr),
- 9 fieldNew(nblock,nfieldv),
- 1 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev),
- 2 enerInternNew(nblock), enerInelasNew(nblock)
- C
- character*80 cmname
- parameter (ZERO=0., ONE=1., TWO=2., half=0.5)
- DIMENSION STRAN(nblock,ndir+nshr)
- DIMENSION C(6,6), CD(6,6)
- E11=props(1)
- E22=props(2)
- E33=props(3)
- V12=props(4)
- V13=props(5)
- V23=props(6)
- G12=props(7)
- G13=props(8)
- G23=props(9)
- XT=props(10)
- XC=props(11)
- YT=props(12)
- YC=props(13)
- ZT=props(14)
- ZC=props(15)
- S12=props(16)
- S13=props(17)
- S23=props(18)
- V21=E22*V12/E11
- V31=E33*V13/E11
- V32=E33*V23/E22
- ATEMPS1=(ONE-V12*V21-V23*V32-V13*V31-TWO*V21*V32*V31)
- 1 /(E11*E22*E33)
- DO m=1,6
- DO n=1,6
- C(m,n)=ZERO
- END DO
- END DO
- ** 刚度矩阵
- C(1,1)=(ONE-V23*V32)/(E22*E33*ATEMPS1)
- C(2,2)=(ONE-V13*V31)/(E11*E33*ATEMPS1)
- C(3,3)=(ONE-V12*V21)/(E11*E22*ATEMPS1)
- C(1,2)=(V12+V32*V13)/(E11*E33*ATEMPS1)
- C(1,3)=(V13+V12*V23)/(E11*E22*ATEMPS1)
- C(2,3)=(V23+V21*V13)/(E11*E22*ATEMPS1)
- C(4,4)=G23
- C(5,5)=G13
- C(6,6)=G12
- DO m=2,6
- DO n=1,m-1
- C(m,n)=C(n,m)
- END DO
- END DO
- C-------STRAIN UPDATE
- DO 100 i=1,nblock
- DO j=1,6
- STRAN(i,j)=ZERO
- END DO
- DO j=1,6
- strainNew(i,j)=strainOld(i,j)+strainInc(i,j)
- END DO
-
- C-------STRESS UPDATE
- stressNew(i,1)=C(1,1)*strainNew(i,1)+C(1,2)*strainNew(i,2)
- 1 +C(1,3)*strainNew(i,3)
- stressNew(i,2)=C(1,2)*strainNew(i,1)+C(2,2)*strainNew(i,2)
- 2 +C(2,3)*strainNew(i,3)
- stressNew(i,3)=C(1,3)*strainNew(i,1)+C(2,3)*strainNew(i,2)
- 3 +C(3,3)*strainNew(i,3)
- stressNew(i,4)=TWO*C(4,4)*strainNew(i,4)
- stressNew(i,5)=TWO*C(5,5)*strainNew(i,5)
- stressNew(i,6)=TWO*C(6,6)*strainNew(i,6)
- DO j=1,13
- stateNew(i,j)=ZERO
- END DO
-
- stateNew(i,1)=StressNew(i,1)
- stateNew(i,2)=StressNew(i,2)
- stateNew(i,3)=StressNew(i,3)
- stateNew(i,4)=StressNew(i,4)
- stateNew(i,5)=StressNew(i,5)
- stateNew(i,6)=StressNew(i,6)
- stateNew(i,7)=(stateNew(i,1)/XT)**2+(stateNew(i,4)/S12)**2
- 1 +(stateNew(i,6)/S13)**2
- stateNew(i,8)=(stateNew(i,1)/XC)**2
- stateNew(i,9)=(stateNew(i,2)/YT)**2+(stateNew(i,4)/S12)**2
- 1 +(stateNew(i,5)/S23)**2
- stateNew(i,10)=(stateNew(i,2)/YC)**2+(stateNew(i,4)/S12)**2
- 1 +(stateNew(i,5)/S23)**2
- stateNew(i,11)=(stateNew(i,1)/XT)**2+(stateNew(i,4)/S12)**2
- 1 +(stateNew(i,6)/S13)**2
- stateNew(i,12)=(stateNew(i,3)/ZT)**2+(stateNew(i,6)/S13)**2
- 1 +(stateNew(i,5)/S23)**2
- stateNew(i,13)=(stateNew(i,3)/ZC)**2+(stateNew(i,6)/S13)**2
- 1 +(stateNew(i,5)/S23)**2
-
-
- IF(stateNew(i,1).GT.0) THEN
- IF(stateNew(i,7).GE.1) THEN
- FT=1
- ELSE
- FT=0
- ENDIF
- ENDIF
-
- IF(stateNew(i,1).LT.0) THEN
- IF(stateNew(i,8).GE.1) THEN
- FC=1
- ELSE
- FC=0
- ENDIF
- ENDIF
- IF(stateNew(i,2).GT.0) THEN
- IF(stateNew(i,9).GE.1) THEN
- MT=1
- ELSE
- MT=0
- ENDIF
- ENDIF
-
- IF(stateNew(i,2).LT.0) THEN
- IF(stateNew(i,10).GE.1) THEN
- MC=1
- ELSE
- MC=0
- ENDIF
- ENDIF
-
- IF(stateNew(i,1).LT.0) THEN
- IF(stateNew(i,11).GE.1) THEN
- MTS=1
- ELSE
- MTS=0
- ENDIF
- ENDIF
-
- IF(stateNew(i,3).GT.0) THEN
- IF(stateNew(i,12).GE.1) THEN
- TS=1
- ELSE
- TS=0
- ENDIF
- ENDIF
-
- IF(stateNew(i,3).LT.0) THEN
- IF(stateNew(i,13).GE.1) THEN
- CS=1
- ELSE
- CS=0
- ENDIF
- ENDIF
-
-
- DO m=1,6
- DO n=1,6
- CD(m,n)=C(m,n)
- END DO
- END DO
- TF=FT+FC+MT+MC+MTS+TS+CS
- C---CD为损伤矩阵
- IF(FT.GE.1) THEN
- E11NEW=0.07*E11
- G12NEW=0.07*G12
- G13NEW=0.07*G13
- V21NEW=E22*V12/E11NEW
- V31NEW=E33*V13/E11NEW
- ATEMPS1NEW=(ONE-V12*V21NEW-V23*V32-V13*V31NEW-TWO*V21NEW*V32*V31NEW)/(
- 1 E11NEW*E22*E33)
- CD(1,1)=(ONE-V23*V32)/(E22*E33*ATEMPS1NEW)
- CD(2,2)=(ONE-V13*V31NEW)/(E11NEW*E33*ATEMPS1NEW)
- CD(3,3)=(ONE-V12*V21NEW)/(E11NEW*E22*ATEMPS1NEW)
- CD(1,2)=(V12+V32*V13)/(E11NEW*E33*ATEMPS1NEW)
- CD(1,3)=(V13+V12*V23)/(E11NEW*E22*ATEMPS1NEW)
- CD(2,3)=(V23+V21NEW*V13)/(E11NEW*E22*ATEMPS1NEW)
- CD(4,4)=G23
- CD(5,5)=G13NEW
- CD(6,6)=G12NEW
- DO m=2,6
- DO n=1,m-1
- CD(m,n)=CD(n,m)
- END DO
- END DO
- END IF
-
- IF(FC.GE.1) THEN
- E11NEW=0.14*E11
- G12NEW=0.07*G12
- G13NEW=0.07*G13
- V21NEW=E22*V12/E11NEW
- V31NEW=E33*V13/E11NEW
- ATEMPS1NEW=(ONE-V12*V21NEW-V23*V32-V13*V31NEW-TWO*V21NEW*V32*V31NEW)/(
- 1 E11NEW*E22*E33)
- CD(1,1)=(ONE-V23*V32)/(E22*E33*ATEMPS1NEW)
- CD(2,2)=(ONE-V13*V31NEW)/(E11NEW*E33*ATEMPS1NEW)
- CD(3,3)=(ONE-V12*V21NEW)/(E11NEW*E22*ATEMPS1NEW)
- CD(1,2)=(V12+V32*V13)/(E11NEW*E33*ATEMPS1NEW)
- CD(1,3)=(V13+V12*V23)/(E11NEW*E22*ATEMPS1NEW)
- CD(2,3)=(V23+V21NEW*V13)/(E11NEW*E22*ATEMPS1NEW)
- CD(4,4)=G23
- CD(5,5)=G13NEW
- CD(6,6)=G12NEW
- DO m=2,6
- DO n=1,m-1
- CD(m,n)=CD(n,m)
- END DO
- END DO
- END IF
-
- IF(MT.GE.1) THEN
- E22NEW=0.3*E22
- G12NEW=0.3*G12
- G23NEW=0.3*G23
- V21NEW=E22NEW*V12/E11
- V32NEW=E33*V23/E22NEW
- ATEMPS1NEW=(ONE-V12*V21NEW-V23*V32NEW-V13*V31-TWO*V21NEW*V32NEW*V31)/(
- 1 E11*E22NEW*E33)
- CD(1,1)=(ONE-V23*V32NEW)/(E22NEW*E33*ATEMPS1NEW)
- CD(2,2)=(ONE-V13*V31)/(E11*E33*ATEMPS1NEW)
- CD(3,3)=(ONE-V12*V21NEW)/(E11*E22NEW*ATEMPS1NEW)
- CD(1,2)=(V12+V32NEW*V13)/(E11*E33*ATEMPS1NEW)
- CD(1,3)=(V13+V12*V23)/(E11*E22NEW*ATEMPS1NEW)
- CD(2,3)=(V23+V21NEW*V13)/(E11*E22NEW*ATEMPS1NEW)
- CD(4,4)=G23NEW
- CD(5,5)=G13
- CD(6,6)=G12NEW
- DO m=2,6
- DO n=1,m-1
- CD(m,n)=CD(n,m)
- END DO
- END DO
- END IF
-
- IF(MC.GE.1) THEN
- E22NEW=0.4*E22
- G12NEW=0.4*G12
- G23NEW=0.4*G23
- V21NEW=E22NEW*V12/E11
- V32NEW=E33*V23/E22NEW
- ATEMPS1NEW=(ONE-V12*V21NEW-V23*V32NEW-V13*V31-TWO*V21NEW*V32NEW*V31)/(
- 1 E11*E22NEW*E33)
- CD(1,1)=(ONE-V23*V32NEW)/(E22NEW*E33*ATEMPS1NEW)
- CD(2,2)=(ONE-V13*V31)/(E11*E33*ATEMPS1NEW)
- CD(3,3)=(ONE-V12*V21NEW)/(E11*E22NEW*ATEMPS1NEW)
- CD(1,2)=(V12+V32NEW*V13)/(E11*E33*ATEMPS1NEW)
- CD(1,3)=(V13+V12*V23)/(E11*E22NEW*ATEMPS1NEW)
- CD(2,3)=(V23+V21NEW*V13)/(E11*E22NEW*ATEMPS1NEW)
- CD(4,4)=G23NEW
- CD(5,5)=G13
- CD(6,6)=G12NEW
- DO m=2,6
- DO n=1,m-1
- CD(m,n)=CD(n,m)
- END DO
- END DO
- END IF
-
- IF(TS.GE.1) THEN
- E22NEW=0.3*E22
- E33NEW=0.3*E33
- G13NEW=0.3*G13
- G23NEW=0.3*G23
- V21NEW=E22NEW*V12/E11
- V31NEW=E33NEW*V13/E11
- V32NEW=E33NEW*V23/E22NEW
- ATEMPS1NEW=(ONE-V12*V21NEW-V23*V32NEW-V13*V31NEW-TWO*V21NEW
- 1 *V32NEW*V31NEW)/(E11*E22NEW*E33NEW)
- CD(1,1)=(ONE-V23*V32NEW)/(E22NEW*E33NEW*ATEMPS1NEW)
- CD(2,2)=(ONE-V13*V31NEW)/(E11*E33NEW*ATEMPS1NEW)
- CD(3,3)=(ONE-V12*V21NEW)/(E11*E22NEW*ATEMPS1NEW)
- CD(1,2)=(V12+V32NEW*V13)/(E11*E33NEW*ATEMPS1NEW)
- CD(1,3)=(V13+V12*V23)/(E11*E22NEW*ATEMPS1NEW)
- CD(2,3)=(V23+V21NEW*V13)/(E11*E22NEW*ATEMPS1NEW)
- CD(4,4)=G23NEW
- CD(5,5)=G13NEW
- CD(6,6)=G12
- DO m=2,6
- DO n=1,m-1
- CD(m,n)=CD(n,m)
- END DO
- END DO
- END IF
-
- IF(CS.GE.1) THEN
- E22NEW=0.4*E22
- E33NEW=0.4*E33
- G13NEW=0.4*G13
- G23NEW=0.4*G23
- V21NEW=E22NEW*V12/E11
- V31NEW=E33NEW*V13/E11
- V32NEW=E33NEW*V23/E22NEW
- ATEMPS1NEW=(ONE-V12*V21NEW-V23*V32NEW-V13*V31NEW-TWO*V21NEW
- 1 *V32NEW*V31NEW)/(E11*E22NEW*E33NEW)
- CD(1,1)=(ONE-V23*V32NEW)/(E22NEW*E33NEW*ATEMPS1NEW)
- CD(2,2)=(ONE-V13*V31NEW)/(E11*E33NEW*ATEMPS1NEW)
- CD(3,3)=(ONE-V12*V21NEW)/(E11*E22NEW*ATEMPS1NEW)
- CD(1,2)=(V12+V32NEW*V13)/(E11*E33NEW*ATEMPS1NEW)
- CD(1,3)=(V13+V12*V23)/(E11*E22NEW*ATEMPS1NEW)
- CD(2,3)=(V23+V21NEW*V13)/(E11*E22NEW*ATEMPS1NEW)
- CD(4,4)=G23NEW
- CD(5,5)=G13NEW
- CD(6,6)=G12
- DO m=2,6
- DO n=1,m-1
- CD(m,n)=CD(n,m)
- END DO
- END DO
- END IF
-
- IF(MTS.GE.1) THEN
- E11NEW=0.07*E11
- E22NEW=0.3*E22
- E33NEW=0.3*E33
- G12NEW=0.07*G12
- G13NEW=0.07*G13
- G23NEW=0.3*G23
- V21NEW=E22NEW*V12/E11NEW
- V31NEW=E33NEW*V13/E11NEW
- V32NEW=E33NEW*V23/E22NEW
- ATEMPS1NEW=(ONE-V12*V21NEW-V23*V32NEW-V13*V31NEW-TWO
- 1 *V21NEW*V32NEW*V31NEW)/(E11NEW*E22NEW*E33NEW)
- CD(1,1)=(ONE-V23*V32NEW)/(E22NEW*E33NEW*ATEMPS1NEW)
- CD(2,2)=(ONE-V13*V31NEW)/(E11NEW*E33NEW*ATEMPS1NEW)
- CD(3,3)=(ONE-V12*V21NEW)/(E11NEW*E22NEW*ATEMPS1NEW)
- CD(1,2)=(V12+V32NEW*V13)/(E11NEW*E33NEW*ATEMPS1NEW)
- CD(1,3)=(V13+V12*V23)/(E11NEW*E22NEW*ATEMPS1NEW)
- CD(2,3)=(V23+V21NEW*V13)/(E11NEW*E22NEW*ATEMPS1NEW)
- CD(4,4)=G23NEW
- CD(5,5)=G13NEW
- CD(6,6)=G12NEW
-
- DO m=2,6
- DO n=1,m-1
- CD(m,n)=CD(n,m)
- END DO
- END DO
- END IF
- stressNew(i,1)=CD(1,1)*strainNew(i,1)+CD(1,2)*strainNew(i,2)
- 1 +CD(1,3)*strainNew(i,3)
- stressNew(i,2)=CD(1,2)*strainNew(i,1)+CD(2,2)*strainNew(i,2)
- 2 +CD(2,3)*strainNew(i,3)
- stressNew(i,3)=CD(1,3)*strainNew(i,1)+CD(2,3)*strainNew(i,2)
- 3 +CD(3,3)*strainNew(i,3)
- stressNew(i,4)=TWO*CD(4,4)*strainNew(i,4)
- stressNew(i,5)=TWO*CD(5,5)*strainNew(i,5)
- stressNew(i,6)=TWO*CD(6,6)*strainNew(i,6)
-
- C Update the specific internal energy -
- stressPower = half * (
- 1 ( stressOld(i,1)+stressNew(i,1) )*strainInc(i,1)
- 2 + ( stressOld(i,2)+stressNew(i,2) )*strainInc(i,2)
- 3 + ( stressOld(i,3)+stressNew(i,3) )*strainInc(i,3)
- 4 + two*( stressOld(i,4)+stressNew(i,4) )*strainInc(i,4)
- 5 + two*( stressOld(i,5)+stressNew(i,5) )*strainInc(i,5)
- 6 + two*( stressOld(i,6)+stressNew(i,6) )*strainInc(i,6))
- C
- enerInternNew(i) = enerInternOld(i)
- 1 + stressPower / density(i)
- 100 continue
-
- return
- end
复制代码
|
|