固化子程序
subroutineusdfld(field,statev,pnewdt,direct,t,celent,
1 time,dtime,cmname,orname,nfield,nstatv,noel,npt,layer,
2 kspt,kstep,kinc,ndi,nshr,coord,jmac,jmtyp,matlayo,laccfla)
c redefine fiel
d variables at a material point.
include 'aba_param.inc'
c
character*80 cmname,orname
character*3 flgray(15)
dimension field(nfield),statev(nstatv),direct(3,3),
1 t(3,3),time(2)
dimension array(15),jarray(15),jmac(*),coord(*) !jmatyp(*),
c this subroutine must call utility routine GETVRM to access material point data.
c Get temperatures from previous increment
field(1)=statev(1)
return
end
C
SUBROUTINE UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,
2 STRAN,DSTRAN,TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,CMNAME,
3 NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,
4 CELENT,DFGRD0,DFGRD1,NOEL,NPT,LAYER,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*8 CMNAME
DIMENSION STRESS(NTENS),STATEV(NSTATV),DDSDDE(NTENS,NTENS),
1 DDSDDT(NTENS),DRPLDE(NTENS),STRAN(NTENS),DSTRAN(NTENS),
2 TIME(2),PREDEF(1),DPRED(1),PROPS(NPROPS),COORDS(3),DROT(3,3),
3 DFGRD0(3,3),DFGRD1(3,3)
DIMENSION PS(NDI),temp(2)
C-----------------------------------------------------------------------------------
C AK即k,AN即n,RF破坏比,c即凝聚力,FA即内摩擦角,PA大气压力,VKB即Kb,VNB即m,AUR 为卸荷Kur/k
C-----------------------------------------------------------------------------------
C statev(1)=temp(1)/100
c 此处的statev(1)=temp(1)/100只是为了调试，真实的固化度不是等于这个
c write(6,*)nodel,temp(1),statev(1)
common uvr,uer,spred1 !定义的全局变量，各子程序之间可传递，每个子程序开头都要定义saa=1.53e5
see=0.665e5
c
srr=8.31
c 普适气体常数
smm=0.813
snn=2.74
scc=43.1
saco=-1.684
sact=5.475e-3
sdensity=1300
svf=0.4
c 纤维体积含量
shr=77500
c 单位质量树脂固化反应放热
IF(Time(2).LT.1.)THEN
DPRED(1) = 0.0
PREDEF(1)=0.055
ELSE
sabc1=((1.-PREDEF(1))**snn)*(PREDEF(1)**smm)
sabc2=scc*(PREDEF(1)-(saco+sact*(TEMP(1)+273)))
sabc3=1+EXP(sabc2)
DPRED(1) = saa*EXP(-see/(srr*(TEMP(1)+273)))*sabc1/sabc3 PREDEF(1) = PREDEF(1)+ DPRED(1)* DTIME
statev(1)=PREDEF(1)
spred1=DPRED(1)* DTIME
C OPEN(UNIT=1,FILE="E:\abaqusTemp\quanguhua\2.TXT")
C write(1,*)dpred(1)
END IF
c 树脂弹性系数参数
UER0=4.67e6
UER1=4.67e9
UTC1A=-45.7
UTC1B=0
UTC2=-12
UTG0=268
UaTG=220
UT0=20
UaER=0
c 树脂弹性模量
UTC1=UTC1A+UTC1B*(TEMP(1)+273)
if(temp(1).lt.50.and.time(2).lt.5000)then
UT=-50
else if(temp(1).ge.50.or.time(2).gt.17782)then
UT=(UTG0+UaTG*statev(1))-(TEMP(1)+273)
end if
if(UT.lt.UTC1)then
UER=UER0
else if(UT.ge.UTC1.and.UT.lt.UTC2)then
UER=UER0+(UT-UTC1)*(UER1-UER0)/(UTC2-UTC1)
else if(UT.ge.UTC2)then
UER=UER1
end if
UVR1=0.37 !树脂泊松比系数
UVR=0.5*(1-UER*(1-2*UVR1)/UER1) !树脂泊松比
UGR=UER/(2*(1+UVR)) !树脂剪切模量
c 纤维弹性系数参数
UEF1=210E9
UEF2=17.24E9
UEF3=17.24E9
UVF13=0.2
UVF12=UVF13
UVF23=0.25
UGF13=27.6E9
UGF23=UEF3/(2*(1+UVF23))
UKR=UER/(2*(1-UVR-2*UVR**2)) !树脂体积模量
UKF=UEF1/(2*(1-UVF12-2*UVF12**2)) !此处纤维体积模量不知道是否用这个参数
UKK11=4*(UVR-UVF13**2)*UKF*UKR*UGR*(1-svf)*svf
UKK12=(UKF+UGR)*UKR+(UKF-UKR)*UGR*svf
UKK1=UKK11/UKK12
E11=svf*UEF1+(1-svf)*UER+UKK1 !层合板1方向弹性模量
UKK21=(UGF13+UGR)+(UGF13-UGR)*SVF
UKK22=(UGF13+UGR)-(UGF13-UGR)*SVF
UKK2=UKK21/UKK22
G12=UGR*UKK2 !层合板12剪切模量
G13=G12
UKK31=UGR*(UKR*(UGR+UGF23)+2*UGF23*UGR+UKR*(UGF23-UGR)*SVF)
UKK32=UKR*(UGR+UGF23)+2*UGF23*UGR-(UKR+2*UGR)*(UGF23-UGR)*SVF
G23=UKK31/UKK32 !层合板23剪切模量
UKK41=(UVR-UVF13)*(UKR-UKF)*UGR*(1-SVF)*SVF
UKK42=(UKF+UGR)*UKR+(UKF-UKR)*UGR*SVF
UKK4=UKK41/UKK42
V12=UVF13*SVF+UVR*(1-SVF)+UKK4 !V12
V21=V12*E23/E11
V13=V12
V31=V13*E23/E11
G31=E11/(2*(1+V31))
UKT1=(UKF+UGR)*UKR+(UKF-UKR)*UGR*SVF
UKT2=(UKF+UGR)-(UKF-UKR)*SVF
UKT=UKT1/UKT2
E23=1/((0.25/UKT)+(0.25/G23)+(V12**2/E11)) !层合板2方向弹性模量
V23=(2*E11*UKT-E11*E23-4*E23*UKT*V13**2)/(2*E11*UKT)
V32=V23
D=(1-V12*V21-V23*V32-V13*V31-2*V21*V32*V13)/(E11*E23*E23)
C if(PREDEF(1).gt.0.829)then
C OPEN(UNIT=1700,FILE="E:\abaqusTemp\quanguhua\2.TXT")
C write(1700,*)e11,e23,g12,g23,v12,v23
C end if
DO 10 K1=1,NTENS
DO 20 K2=1,NTENS
DDSDDE(K1,K2)=0
20 END DO
10 END DO
DDSDDE(1,1)=(1-V23*V32)/(E23*E23*D)
DDSDDE(1,2)=(V12+V32*V13)/(E11*E23*D)
DDSDDE(1,3)=(V13+V23*V12)/(E11*E23*D)
DDSDDE(2,2)=(1-V13*V31)/(E11*E23*D)
DDSDDE(2,1)=DDSDDE(1,2)
DDSDDE(2,3)=(V23+V21*V13)/(E11*E23*D)
DDSDDE(3,1)=DDSDDE(1,3)
DDSDDE(3,2)=DDSDDE(2,3)
DDSDDE(3,3)=(1-V12*V21)/(E11*E23*D)
DDSDDE(4,4)=G23
DDSDDE(5,5)=G31
DDSDDE(6,6)=G12
DO 70 K1=1,NTENS
DO 60 K2=1,NTENS
STRESS(K2)=STRESS(K2)+DDSDDE(K2,K1)*DSTRAN(K1)
60 CONTINUE
70 CONTINUE
RETURN
END
C
subroutineuexpan(expan,dexpandt,temp,time,dtime,predef,dpred,
$ statev,cmname,nstatv,noel)
c
include 'aba_param.inc'
c
character*80 cmname
c
dimension expan(*),dexpandt(*),temp(2),time(2),predef(*),
$ dpred(*),statev(nstatv)
c
common uvr,uer,spred1 !定义的全局变量，各子程序之间可传递，每个子程序开头都要定义ra=0.173 !常数A
rvp=0.099 !固化完成时的单位体积树脂固化收缩量
rvf=0.4 !纤维体积含量
rvr=0.37 !树脂泊松比
ragel=0.469 !凝胶点固化度
rac1=0.055 !树脂固化收缩开始点
rac2=0.670 !树脂固化收缩结束点
rac=(predef(1)-rac1)/(rac2-rac1)
rv=1+ra*rac+(rvp-ra)*rac**2
rv2=rv**2
rv3=ra/(rac2-rac1)+2*(rvp-ra)*rac/(rac2-rac1)
rv4=3*rv2**(1./3.)
chexpanr1=0. ! 树脂1方向上固化收缩应变系数
chexpanr2=(rv3/rv4)*(1-rvf)*(1+0.37) ! 树脂1方向上固化收缩应变系数，次数statev（3）是umat里面计算出来的树脂泊松比
rexp11=0.6e-6 !单层板1方向上在温度t1时的热胀系数
rexp12=0.6e-6 !单层板1方向上在温度t2时的热胀系数
rexp21=28.6e-6 !单层板2方向上在温度t2时的热胀系数
rexp22=28.6e-6 !单层板1方向上在温度t2时的热胀系数
rexp111=5.04e-6 !纤维在1方向上的热胀系数
rexp112=5.04e-6 !纤维在2方向上的热胀系数
rexp211=7.2e-5 !树脂在1方向上的热胀系数
rexp212=7.2e-5 !树脂在2方向上的热胀系数
rt1=0 !温度t1
rt2=100 !温度t2
c 复合材料三个方向上热胀系数
thexpan1=(rexp12-rexp11)*(temp(1)-rt1)/(rt2-rt1)+rexp11 thexpan2=(rexp22-rexp21)*(temp(1)-rt1)/(rt2-rt1)+rexp21 C thexpan1=rvf*rexp111+(1-rvf)*rexp211
C thexpan2=rvf*rexp112+(1-rvf)*rexp212
UVF13=0.2 !纤维13方向泊松比
UEF1=210E9 !纤维1方向弹性模量
rexpan21=(UVF13*rvf+uvr*(1-rvf))*chexpanr2*UER rexpan22=rexpan21*(1-rvf)/(UEF1*rvf+UER*(1-rvf)) rexpan2=chexpanr2*(1+uvr)*(1-rvf)-rexpan22
chexpan2=rexpan2
if(predef(1).gt.ragel.and.predef(1).lt.rac2)then
expan(1) = thexpan1*temp(2)
expan(2) = -chexpan2*spred1+thexpan2*temp(2)
expan(3) = expan(2)
else if(predef(1).ge.rac2)then
expan(1) = thexpan1*temp(2)
expan(2) = thexpan2*temp(2)
expan(3) = expan(2)
end if
predef(2)=chexpanr2
predef(3)=spred1
predef(4)=chexpan2*spred1
C predef(5)=thexpan2*temp(2)
C predef(6)=expan(2)
C if(time(2).ge.16000)then
C expan(1) = chexpan1*dpred(1)
C expan(2) = chexpan2*dpred(1)
C expan(3) = expan(2)
C end if
return
end