复合材料层合板三维Hashin失效准则VUMAT子程序调试

wzx901221 · 发表于 2017-3-20 20:12:27

本人根据三维Hashin失效准则和Chang-Chang退化准则编写了一段VUMAT子程序，但一直未能调试成功，求论坛里各位大神指点帮助啊！做过类似工作的我们也可以互相交流！

glove831 · 发表于 2017-4-6 20:10:20

楼主可以去请教一下yiqinuli楼主，他比较擅长VUMAT
http://home.simwe.com/space-uid-619943.html

wzx901221 · 发表于 2017-3-20 21:00:19

想交流朋友的可以私信我，我会把我害存在错误的程序分享，共同学习！

wzx901221 · 发表于 2017-3-21 23:31:13

本帖最后由 wzx901221 于 2017-3-22 22:22 编辑

复制代码

kqvtkdss · 发表于 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

复制代码

阿幺 · 发表于 2017-6-9 17:37:50

你计算的有什么问题？ABAQUS自带的有三维hashin的uvmat

hitery · 发表于 2017-8-25 10:04:38

阿幺发表于 2017-6-9 17:37
你计算的有什么问题？ABAQUS自带的有三维hashin的uvmat

自带的三维hashin准则umat在哪里啊，怎么调用？

upczhou · 发表于 2017-8-30 10:38:56

hitery 发表于 2017-8-25 10:04
自带的三维hashin准则umat在哪里啊，怎么调用？

链接：http://pan.baidu.com/s/1ge7IKgz 密码：sljs

optimus103 · 发表于 2017-10-15 21:07:26

kqvtkdss 发表于 2017-3-22 12:02

这是基于参数退化的三维Hashin VUMAT渐进失效吧？

simon1106 · 发表于 2018-5-5 09:59:54

请问用Vumat的话必须删除单元吗？这样的话就看不到纤维失效的云图了（直接删除了），只能看到基体和层间失效云图

Linda010 · 发表于 2018-5-21 00:31:48

楼主程序调通了吗？？

godlovedy · 发表于 2018-10-23 18:19:16

楼主程序调通了吗？？

侠在江湖 · 发表于 2019-1-24 14:13:48

楼主程序怎么样了啊

wfylzy · 发表于 2019-2-7 01:08:45

阿幺发表于 2017-6-9 17:37
你计算的有什么问题？ABAQUS自带的有三维hashin的uvmat

你好请问你说的这个是unifiber吗

somehow · 发表于 2019-3-2 12:04:58

最近在学编子程序楼主可以给个联系方式吗？有问题请教

莫与陌与沫 · 发表于 2019-4-16 10:55:33

hitery 发表于 2017-8-25 10:04
自带的三维hashin准则umat在哪里啊，怎么调用？

层主，自带的三维Hashin准则umat子程序你找到了吗？能不能分享一下呀

IsaacHuo · 发表于 2019-8-30 11:34:51

这代码的应变更新有问题吧，还有剪应变没有改为工程应变。

nuaa2013 · 发表于 2019-9-9 20:15:37

upczhou 发表于 2017-8-30 10:38
链接：http://pan.baidu.com/s/1ge7IKgz 密码：sljs

请问有CAE的模型文件吗？我使用这个子程序调用出了问题，想试一下，谢谢啦。(有奖回复哦)
qq有效联系：269102349@qq.com

bluebluesky · 发表于 2021-3-29 21:41:39

不错的分享

蝶恋花 · 发表于 2021-4-6 12:10:24

VUMAT会删除单元

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

[复合材料] 复合材料层合板三维Hashin失效准则VUMAT子程序调试