非線性彈性全量模型(江見鯨模型) 脆性斷裂混凝土本構(gòu)程序

1、非線性彈性全量模型(江見鯨模型)+脆性斷裂混凝土本構(gòu)程序TypeDef !定義混凝土材性模塊type : typ_Concretereal*8 fc, ft,E0, ENU, EPS_Crush ; !抗壓強度+，抗拉強度+，初始彈性模量，初始泊松比，壓碎應變-real*8 Es, ENUs !割線模量，割線泊松比real*8 T(6,6) !坐標轉(zhuǎn)換矩陣integer NCrack (3), Pre_NCrack(3), Pre_inc, Pre_incsub; !開裂記錄，前次迭代開裂記錄，前次增量步，前次增量子步real*8 SIG(6), EPS(6),dEPS(6); !開始時應力，

2、應變，應變增量real*8 StressP(6), StrainP(6); !主應力，主應變real*8 Stress(6), Strain(6) !結(jié)束時應力，應變real*8 Beta,Pre_Beta !非線性指標，前次迭代非線性指標real*8 D(6,6), Dela(6,6), Ds(6,6) !剛度矩陣，彈性剛度矩陣，割線剛度矩陣end type typ_Concreteend module TypeDef module My_MOD !開辟公共變量空間use TypeDeftype(typ_Concrete) : My_Con(1000,8) !定義混凝土數(shù)組end modul

3、esubroutine Get_DS(D,G,E,DE,S,TEMP0,1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR,inc,incsub,ncycle)! D(6x6) 迭代本構(gòu)矩陣(out)! G(6) 由于狀態(tài)改變引起的應力變化，不用(out)! E(6) 開始時刻的應變(in)! DE(6) 應變增量(in)! S(6) 開始時刻的應變(in & out)! Temp0 溫度(in)! DTEMP 溫度變化(in)! NGENS 應變維數(shù)(in)! N(2) 單元編號(in)! NN 積分點編號(in)! KC 層號(in)! MATS 材料編號(

4、in)! NDI 正應力維數(shù)(in)! NSHEAR 剪應力維數(shù)(in)! inc 當前增量步(in)! incsub 當前增量子步(in)! ncycle 當前循環(huán)數(shù)(in)use IMSL !引用IMSL函數(shù)庫use typedef use My_Modimplicit noneinteger : ngens,nn,kc,mats,ndi,nshear,inc,incsub,ncyclereal*8 : e(ngens),de(ngens),temp0(1),dtemp(1),g(ngens)1 ,d(ngens,ngens),s(ngens) integer : n(2)type(typ

5、_concrete) : Creal*8 Beta1,strain_mreal*8 s_m,J2,J3,r,sita,TempA, TempB, TempC;integer NSubStep !子步積分步數(shù)integer I,J,K1,K2C=My_Con(n(1),nn) !得到內(nèi)存中保留的數(shù)據(jù)c -c 初次計算，清零并賦值if(inc=0.and.incsub=0.and.ncycle=0) thenopen(77,file='debug.txt')write(77,*)close(77)C%fc=30.; C%ft=3.; C%E0=30e3; C%ENU=0.18;C%

6、EPS_Crush=-0.0033; C%T=0.; C%NCrack=0; C%Pre_NCrack=0; C%Beta=0; C%Pre_Beta=0;C%Pre_inc=0; C%Pre_incsub=0; C%SIG=0.; C%EPS=0.; C%Stress=0.; C%Strain=0.;end ifc -c 如果新的增量步開始，則更新相應變量if(inc>C%Pre_inc .or. incsub>C%Pre_incsub) thenC%Pre_inc=inc; C%Pre_incsub=incsubC%NCrack=C%Pre_NCrack; !修正裂縫狀態(tài)C%B

7、eta=C%Pre_Beta; !修正非線性指標狀態(tài)! 判斷是否壓壞strain_m=(C%EPS(1)+C%EPS(2)+C%EPS(3)/3.if(Strain_m>0.) Strain_m=0.if(minval(C%EPS(1:3)-Strain_m)<C%EPS_Crush) thenC%NCrack=100 !徹底破壞end ifend ifc 數(shù)據(jù)賦值 open(77,file='debug.txt',position='append')C%SIG=s; C%EPS=e; C%dEPS=de;C%Pre_NCrack=C%NCrackC

8、%Pre_Beta=C%BetaNSubStep=4c -c 計算彈性矩陣C%Dela=0.do K1=1, 3do K2=1, 3C%Dela(K2,K1)=C%ENUend doC%Dela(K1,K1)=1.-C%ENUEND dodo K1=4,6C%Dela(K1,K1)=(1.-2.*C%ENU)*0.5end doC%Dela=C%Dela*C%E0/(1.+C%ENU)/(1.-2.*C%ENU)c -c 如果已經(jīng)壓碎，應力清零，剛度為很小值，結(jié)束計算if(maxval(C%NCrack)=100) thenC%D=0.0001*C%DelaC%SIG=0.;C%Stress=

9、0.;s=0.returnend ifC 計算主應力和割線剛度C%Stress=C%SIGdo I=1, NSubSteps_m=(C%Stress(1)+C%Stress(2)+C%Stress(3)/3. !計算平均應力s(1:3)=C%Stress(1:3)-s_ms(4:6)=C%Stress(4:6) !計算應力偏量J2=-s(1)*s(2)-s(2)*s(3)-s(3)*s(1)+s(4)*2+s(5)*2+s(6)*2 !計算J2J3=s(1)*s(2)*s(3)+2.*s(4)*s(5)*s(6) ! 計算J31 -s(1)*s(5)*2-s(2)*s(6)*2-s(3)*s(

10、4)*2r=sqrt(4.*J2/3.)if(r.ne.0.) thensita=acos(4.*J3/r*3)/3.elsesita=0.end ifif(maxval(abs(C%Pre_NCrack)=0) then !沒有裂縫call Get_T_Matrix(C%SIG,C%T) !計算坐標轉(zhuǎn)換矩陣end ifC%StressP=matmul(transpose(C%T),C%Stress); !計算主應力C%StrainP=matmul(transpose(C%T),C%Strain); !計算主應變if(C%StressP(1)<0.05*C%fc.and.1 maxval

11、(abs(C%Pre_NCrack)=0) then !沒有裂縫TempA=1.2856/C%fc*2;TempB=(1.4268+10.2551*cos(sita)/C%fc;TempC=3.2128*s_m*3./C%fc-1.;Beta1=-TempB+sqrt(TempB*2-4.*TempA*TempC)Beta1=Beta1/2./TempABeta1=sqrt(J2)/Beta1 !根據(jù)江見鯨模型求解Betaif(Beta1>C%Beta) then !Beta應該始終增大(對于全量模型)C%Pre_Beta=Beta1elseBeta1=C%Betaend if! 計算割

12、線剛度和泊松比if(Beta1>1.) Beta1=1.C%Es=C%E0/2.*(1.+sqrt(1.-Beta1)if(Beta1<0.8) C%ENUs=C%ENUif(Beta1.ge.0.8) C%ENUs=0.42-(0.42-C%ENU)*1 sqrt(1.-(Beta1-.8)/.2)*2)!計算割線剛度矩陣C%Ds=0.do K1=1, 3do K2=1, 3C%Ds(K2,K1)=C%ENUsend doC%Ds(K1,K1)=1.-C%ENUsend dodo K1=4,6C%Ds(K1,K1)=(1.-2.*C%ENUs)*0.5end doC%Ds=C%D

13、s*C%Es/(1.+C%ENUs)/(1.-2.*C%ENUs)else !如果處于開裂控制區(qū)C%Ds=C%Delado K1=1,3if(C%StressP(K1)>C%ft .OR. C%Pre_NCrack(K1)>0) then !按開裂處理C%Pre_NCrack(K1)=1;Call Crack_Open(C,K1) !計算開裂矩陣end ifend doC%Ds=matmul(C%T,matmul(C%Ds,transpose(C%T); !計算割線剛度矩陣end ifif(Beta1<0.99999d0) then !如果沒有達到極限應力C%Strain=C

14、%EPS+C%dEPS*real(I/NSubStep);C%Stress=matmul(C%Ds,C%Strain);else !達到極限應力后應力不變C%Stress=C%SIGC%Ds=1.d-6*C%Dsend ifend docccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 設置迭代剛度矩陣 D=c%Ds+1.d-6*C%Delas=C%StressMy_Con(N(1),nn)=C;c c write(77,*) 'Step: ',inc, incsub, ncy

15、clec write(77,*) Beta1c write(77,*) C%Pre_NCrackc write(77,*) C%Stress(1:3)c write(77,*)close(77)returnend subroutineC 根據(jù)開裂修正剛度矩陣subroutine Crack_Open(C,I0)use typeDefimplicit nonetype (typ_Concrete) : Cinteger,intent(in) : I0real*8 FACTOR,ECRA1,ECRA2,ECRA12C%Ds(I0,:)=0.C%Ds(:,I0)=0.ECRA1=C%StrainP(

16、I0);ECRA2=0.;ECRA12=maxval(abs(C%StrainP(4:6)c 計算裂面剪力傳遞系數(shù)call USHRET0 (FACTOR,ECRA1,ECRA2,ECRA12)C%Ds(4,4)=FACTOR*C%Ds(4,4)C%Ds(5,5)=FACTOR*C%Ds(5,5)C%Ds(6,6)=FACTOR*C%Ds(6,6)returnend subroutine c 計算裂面剪力傳遞系數(shù)SUBROUTINE USHRET0 (FACTOR,ECRA1,ECRA2,ECRA12)IMPLICIT REAL *8 (A-H, O-Z)factor=0.4;RETURNEN

17、Dc 計算主應力及坐標轉(zhuǎn)換矩陣subroutine Get_T_Matrix(olds,T)use IMSLimplicit nonereal*8 olds(6), T(6,6)real*8 SIG(3,3),EVAL(3), EVEC(3,3)real*8 l(3),m(3),n(3)real*8 SIGP(3)integer ISIG(1,1)=olds(1)SIG(2,2)=olds(2)SIG(3,3)=olds(3)SIG(1,2)=olds(4); SIG(2,1)=olds(4);SIG(2,3)=olds(5); SIG(3,2)=olds(5);SIG(1,3)=olds(6

18、); SIG(3,1)=olds(6);call DEVCSF(3, SIG, 3, EVAL, EVEC, 3)SIGP(1)=maxval(EVAL)do I=1,3if(EVAL(I)=SIGP(1) thenl=EVec(:,I)EVAL(I)=minval(EVAL)-10;exitend ifend doSIGP(2)=maxval(EVAL)do I=1,3if(EVAL(I)=SIGP(2) thenm=EVec(:,I)EVAL(I)=minval(EVAL)-10;exitend ifend doSIGP(3)=maxval(EVAL)do I=1,3if(EVAL(I)=

19、SIGP(3) thenn=EVec(:,I)EVAL(I)=minval(EVAL)-10;exitend ifend dodo I=1,3T(I,:)=(/l(I)*2,m(I)*2,n(I)*2,l(I)*m(I),1 m(I)*n(I),n(I)*l(I)/)end doT(4,:)=(/2.d0*l(1)*l(2),2.d0*m(1)*m(2),2.d0*n(1)*n(2),1 l(1)*m(2)+l(2)*m(1),m(1)*n(2)+m(2)*n(1),n(1)*l(2)+n(2)*l(1)/)T(5,:)=(/2.d0*l(2)*l(3),2.d0*m(2)*m(3),2.d0

20、*n(2)*n(3),1 l(2)*m(3)+l(3)*m(2),m(2)*n(3)+m(3)*n(2),n(2)*l(3)+n(3)*l(2)/)T(6,:)=(/2.d0*l(3)*l(1),2.d0*m(3)*m(1),2.d0*n(3)*n(1),1 l(3)*m(1)+l(1)*m(3),m(3)*n(1)+m(1)*n(3),n(3)*l(1)+n(1)*l(3)/)returnend subroutinec marc 接口程序SUBROUTINE HYPELA(D,G,E,DE,S,TEMP0,1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR)impli

21、cit real*8 (a-h,o-z)INCLUDE './common/concom' ! 通過concom模塊得到當前的計算步數(shù)integer : ngens,nn,kc,mats,ndi,nshearreal*8 : e(1),de(1),temp0(*),dtemp(*),g(1),d(ngens,ngens),s(1) integer : n(2) if(mats=1) then !如果材料編號是1call Get_DS(D,G,E,DE,S,TEMP0,1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR,inc,incsub,ncycle)end ifreturnend subroutinec 后處理子程序 subroutine plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,* nshear,jpltcd)c* * * * * *cc select a variable contour plotting (user subroutine).cc v variablec s (idss) stress arrayc sp stresses in preferred directionc etot t