umat程序代码及注释文档格式.docx
《umat程序代码及注释文档格式.docx》由会员分享，可在线阅读，更多相关《umat程序代码及注释文档格式.docx（16页珍藏版）》请在冰豆网上搜索。
C输入材料的强度参数，用于失效判断。
xt=props（10）
xc=props（11）
yt=props（12）
yc=props（13）
s12=props（14）
s13=props（15）
s23=props（16）
zt=props（17）
zc=props（18）
eta=props（19）粘性系数
gc1=props（20）复合材料的1方向的临界能量释放率
gc2=props（21）
gc3=props（22）
!
CFULL无损伤下的刚度阵
u21=（u12/e1）*e2
u31=（u13/e1）*e3
u32=（u23/e2）*e3
fact=1-u12*u21-u23*u32-u13*u31-
&
2*u21*u32*u13
factor=fact/（e1*e2*e3）
CFULL（1,1）=（1-u23*u32）/（e2*e3*factor）
CFULL（2,2）=（1-u13*u31）/（e1*e3*factor）
CFULL（3,3）=（1-u12*u21）/（e1*e2*factor）
CFULL（1,2）=（u12+u32*u13）/（e1*e3*factor）
CFULL（2,3）=（u32+u12*u31）/（e1*e3*factor）
CFULL（1,3）=（u31+u21*u32）/（e3*e2*factor）
CFULL（3,1）=CFULL（1,3）
CFULL（3,2）=CFULL（2,3）
CFULL（2,1）=CFULL（1,2）
CFULL（4,4）=g12
CFULL（5,5）=g13
CFULL（6,6）=g23
C
CRECOVERELASTICSTRAINS更新应变
DO10K1=1,NTENS
STRANT（K1）=STATEV（K1）+DSTRAN（K1）
10CONTINUE
CSTORESTRAINSINSTATEVARIABLEARRAY存储应变值，用于下一步的应变更新
DO20K1=1,NTENS
STATEV（K1）=STRANT（K1）
20CONTINUE
damagestrain计算各种失效情况下发生失效时的应变极限值
epsilonxt=xt/CFULL（1,1）
epsilonxc=xc/CFULL（1,1）
epsilonyt=yt/CFULL（2,2）
epsilonyc=yc/CFULL（2,2）
epsilonzt=zt/CFULL（3,3）
epsilonzc=zc/CFULL（3,3）
epsilons12=s12/CFULL（4,4）
epsilons13=s13/CFULL（5,5）
epsilons23=s23/CFULL（6,6）
fibertensedamage失效准则
if（STRANT
（1）.GE.0.0）theneft=（STRANT
（1）/epsilonxt）**2+（STRANT（4）/epsilons12）**2+
（STRANT（5）/epsilons13）**2
efc=0.0
else
eft=0.0
fibercompressbuckling
efc=（STRANT
（1）/xc）**2
endif
matrixcrack
epm=STRANT
（2）+STRANT（3）
if（epm.GE.0.0）then
emt=（STRANT
（2）+STRANT（3））**2/（epsilonyt*epsilonzt）-
STRANT
（2）*STRANT（3）/epsilons23**2+
（STRANT（4）/epsilons12）**2+
（STRANT（5）/epsilons13）**2+
（STRANT（6）/epsilons23）**2
emc=0.0
emt=0.0
matrixcruck
emc=（STRANT
（2）+STRANT（3））**2/（epsilonyc*epsilonzc）+
（（STRANT
（2）+STRANT（3））/epsilonyc）*
（epsilonyc/2.0/epsilons12-1.0）-
endif
!
delemination
if（STRANT（3）.GE.0.0）then
edt=（STRANT（3）/epsilonzt）**2+（STRANT（6）/epsilons23）**2+
edc=0.0
edt=0.0
edc=（STRANT（3）/epsilonzc）**2+（STRANT（6）/epsilons23）**2+
Open（16,File='
D:
/1.txt'
status='
unknown'
）
write（16,*）ef,em,ed
ef=eft+efc
em=emt+emc
ed=edt+edc
damageornot判断是否发生失效
fiberdamage
if（ef.GT.1.0）then纤维发生损伤
ff=SQRT（ef）计算ff
if（STRANT
（1）.GE.0.0）then若为纤维拉断
df=1.0-exp（（-1.0）*xt*epsilonxt*CELENT*（ff-1.0）/gc1）/ff计算纤维的损伤变量df
dfdff=xt*epsilonxt*CELENT/gc1*计算
，更新雅克比矩阵用
exp（（-1.0）*xt*epsilonxt*CELENT*（ff-1.0）/gc1）/ff+
exp（（-1.0）*xt*epsilonxt*CELENT*（ff-1.0）/gc1）/ff**2
DOi=1,6
DFFDE（i）=0.0计算
ENDDO
DFFDE
（1）=1.0/ff*STRANT
（1）/epsilonxt**2
DFFDE（4）=1.0/ff*STRANT（4）/epsilons12**2
DFFDE（5）=1.0/ff*STRANT（5）/epsilons13**2
if（df.GT.STATEV（7））then因为损伤不可逆，保存df的最大值，用到粘性规律中
STATEV（7）=df
else
df=STATEV（7）
endif
else若为纤维压缩屈曲
df=1.0-exp（（-1.0）*xc*epsilonxc*CELENT*（ff-1.0）/gc1）/ff计算纤维的损伤变量df
dfdff=xc*epsilonxc*CELENT/gc1*计算
exp（（-1.0）*xc*epsilonxc*CELENT*（ff-1.0）/gc1）/ff+
exp（（-1.0）*xc*epsilonxc*CELENT*（ff-1.0）/gc1）/ff**2
DOi=1,6计算
DFFDE（i）=0.0
DFFDE
（1）=1.0/ff*STRANT
（1）/epsilonxc**2
if（df.GT.STATEV（7））then因为损伤不可逆，保存df的最大值，用到粘性规律中
else纤维没有发生损伤，纤维的损伤变量等都置零
df=0.0
dfdff=0.0
matrixdamage基体损伤，思路与纤维损伤相同
if（em.GT.1.0）then
fm=SQRT（em）
if（epm.GE.0.0）then
dm=1.0-exp（（-1.0）*yt*epsilonyt*CELENT*（fm-1.0）/gc2）/fm
dmdfm=yt*epsilonyt*CELENT/gc2*
exp（（-1.0）*yt*epsilonyt*CELENT*（fm-1.0）/gc2）/fm+
exp（（-1.0）*yt*epsilonyt*CELENT*（fm-1.0）/gc2）/fm**2
DFMDE（i）=0.0
DFMDE
（2）=1.0/2.0/fm*
（2.0*（STRANT
（2）+STRANT（3））/epsilonyt/epsilonzt-
STRANT（3）/epsilons23**2）
DFMDE（3）=1.0/2.0/fm*
（2.0*（STRANT
（2）+STRANT（3））/epsilonyt/epsilonzt-
STRANT
（2）/epsilons23**2）
DFMDE（4）=1.0/fm*STRANT（4）/epsilons12**2
DFMDE（5）=1.0/fm*STRANT（5）/epsilons13**2
DFMDE（6）=1.0/fm*STRANT（6）/epsilons23**2
if（dm.GT.STATEV（8））then
STATEV（8）=dm
dm=STATEV（8）
dm=1.0-exp（-1.0*yc*epsilonyc*CELENT*（fm-1.0）/gc2）/fm
dmdfm=yc*epsilonyc*CELENT/gc2*
exp（（-1.0）*yc*epsilonyc*CELENT*（fm-1.0）/gc2）/fm+
exp（（-1.0）*yc*epsilonyc*CELENT*（fm-1.0）/gc2）/fm**2
（2.0*（STRANT
（2）+STRANT（3））/epsilonyc/epsilonzc+
1.0/epsilonyc*（1.0/2.0*epsilonyc/epsilons12-1）-
STRANT（3）/epsilons23**2）
DFMDE（4）=1.0/fm*
STRANT（4）/epsilons12**2
DFMDE（5）=1.0/fm*
STRANT（5）/epsilons13**2
DFMDE（6）=1.0/fm*
STRANT（6）/epsilons23**2
endif
dm=0.0
dmdfm=0.0
delamination分层损伤，思路同上
if（ed.GT.1.0）then
fd=SQRT（ed）
dd=1-exp（（-1.0）*zt*epsilonzt*CELENT*（fd-1.0）/gc3）/fd
dddfd=zt*epsilonzt*CELENT/gc3*
exp（（-1.0）*zt*epsilonzt*CELENT*（fd-1.0）/gc3）/fd+
exp（（-1.0）*zt*epsilonzt*CELENT*（fd-1.0）/gc3）/fd**2
DFDDE（i）=0.0
DFDDE（3）=1.0/fd
*STRANT（3）/epsilonzt**2
DFDDE（5）=1.0/fd
*STRANT（5）/epsilons13**2
DFDDE（6）=1.0/fd
*STRANT（6）/epsilons23**2
if（dd.GT.STATEV（9））then
STATEV（9）=dd
dd=STATEV（9）
dd=1.0-exp（（-1.0）*zc*epsilonzc*CELENT*（fd-1.0）/gc3）/fd
dddfd=zc*epsilonzc*CELENT/gc3*
exp（（-1.0）*zc*epsilonzc*CELENT*（fd-1.0）/gc3）/fd+
exp（（-1.0）*zc*epsilonzc*CELENT*（fd-1.0）/gc3）/fd**2
*STRANT（3）/epsilonzc**2
dd=0.0
dddfd=0.0
if（df.GT.1.0）then
df=1.0
if（dm.GT.1.0）then
dm=1.0
if（dd.GT.1.0）then
dd=1.0
dfv=df*DTIME/（eta+DTIME）+STATEV（10）*eta/（eta+DTIME）施加粘性规律，得dfv等
dmv=dm*DTIME/（eta+DTIME）+STATEV（11）*eta/（eta+DTIME）
ddv=dd*DTIME/（eta+DTIME）+STATEV（12）*eta/（eta+DTIME）
STATEV（10）=dfv保存dfv
STATEV（11）=dmv
STATEV（12）=ddv
write（16,*）dfv
CALCULATECDFULL计算折减后的刚度阵
CDFULL（1,1）=（1.0-dfv）*（1.0-dfv）*CFULL（1,1）
CDFULL（2,2）=（1.0-dmv）*（1.0-dmv）*CFULL（2,2）
CDFULL（3,3）=（1.0-ddv）*（1.0-ddv）*CFULL（3,3）
CDFULL（1,2）=（1.0-dfv）*（1.0-dmv）*CFULL（1,2）
CDFULL（1,3）=（1.0-dfv）*（1.0-ddv）*CFULL（1,3）
CDFULL（2,3）=（1.0-dmv）*（1.0-ddv）*CFULL（2,3）
CDFULL（2,1）=CDFULL（1,2）
CDFULL（3,1）=CDFULL（1,3）
CDFULL（3,2）=CDFULL（2,3）
CDFULL（4,4）=（1.0-dfv）*（1.0-dmv）*CFULL（4,4）
CDFULL（5,5）=（1.0-dfv）*（1.0-ddv）*CFULL（5,5）
CDFULL（6,6）=（1.0-dmv）*（1.0-ddv）*CFULL（6,6）
Open（16,File='
write（16,*）CDFULL（1,1）,CDFULL（2,2）
CALCULATEDCDDFV计算
，用于更新雅克比矩阵
DOj=1,6
DCDDFV（i,j）=0.0
ENDDO
DCDDFV（1,1）=（2.0*dfv-2.0）*CFULL（1,1）
DCDDFV（1,2）=（dmv-1.0）*CFULL（1,2）
DCDDFV（1,3）=（ddv-1.0）*CFULL（1,3）
DCDDFV（2,1）=DCDDFV（1,2）
DCDDFV（3,1）=DCDDFV（1,3）
DCDDFV（4,4）=（dmv-1.0）*CFULL（4,4）
DCDDFV（5,5）=（ddv-1.0）*CFULL（5,5）
CALCULATEDCDDMV计算
DCDDMV（i,j）=0.0
DCDDMV（2,2）=（2.0*dmv-2.0）*CFULL（2,2）
DCDDMV（1,2）=（dfv-1.0）*CFULL（1,2）
DCDDMV（2,3）=（ddv-1.0）*CFULL（2,3）
DCDDMV（2,1）=DCDDMV（1,2）
DCDDMV（3,2）=DCDDMV（2,3）
DCDDMV（4,4）=（dfv-1.0）*CFULL（4,4）
DCDDMV（6,6）=（ddv-1.0）*CFULL（6,6）
CALCULATEDCDDDV计算
DCDDDV（i,j）=0.0
DCDDDV（3,3）=（2.0*ddv-2.0）*CFULL（3,3）
DCDDDV（1,3）=（dfv-1.0）*CFULL（1,3）
DCDDDV（2,3）=（dmv-1.0）*CFULL（2,3）
DCDDDV（3,1）=DCDDDV（1,3）
DCDDDV（3,2）=DCDDDV（2,3）
DCDDDV（5,5）=（dfv-1.0）*CFULL（5,5）
DCDDDV（6,6）=（dmv-1.0）*CFULL（6,6）
CCALCULATESTRESSFROMSTRAINS通过折减后的刚度阵及应变更新应力
UpdatetheStrsee
stress
（1）=CDFULL（1,1）*STRANT
（1）
+CDFULL（1,2）*STRANT
（2）
+CDFULL（1,3）*STRANT（3）
stress
（2）=CDFULL（2,1）*STRANT
（1）
+CDFULL（2,2）*STRANT
（2）
+CDFULL（2,3）*STRANT（3）
stress（3）=CDFULL（1,3）*STRANT
（1）
+CDFULL（2,3）*STRANT
（2）
+CDFULL（3,3）*STRANT（3）
stress（4）=CDFULL（4,4）*STRANT（4）
stress（5）=CDFULL（5,5）*STRANT（5）
stress（6）=CDFULL（6,6）*STRANT（6）
DOi=1,6计算
TEMP1（i）=TEMP1（i）+DCDDFV（i,j）*STRANT（j）
ENDDO
ENDDO
DOi=1,6计算
TEMP2（i）=TEMP2（i）+DCDDMV（i,j）*STRANT（j）
DOi=1,6计算
TEMP3（i）=TEMP3（i）+DCDDDV（i,j）*STRANT（j）
TEMP4（i）=ddfdff*DFFDE（i）
TEMP5（i）=ddmdfm*DFMDE（i）计算
TEMP6（i）=ddddfd*DFDDE（i）
计算（
+
升级会员 