邓肯张模型FORTRAN子程序源代码
- 格式:doc
- 大小:41.00 KB
- 文档页数:4
(1) 根据邓肯等人总结的经验公式计算参数a 、b :b =1(σ1−σ3)ult =(ε1σ1−σ3)95%−(ε1σ1−σ3)70%(ε1)95%−(ε1)70%()()111195%70%13131395%70%112a 1i a a a ultE p p p εεεεσσσσσσ==⎛⎫⎛⎫⎛⎫⎡⎤+-+ ⎪ ⎪ ⎪⎣⎦---⎝⎭⎝⎭⎝⎭()131313ult()()-ff fR b σσσσσσ-==-计算得到表一如下。
f 80117.03321=++=Rf Rf Rf Rf又因为a 为起始变形模量E i 的倒数,即E i =1a可得表二,并绘制lg (Ei/Pa) 与lg(σ3/Pa)的试验关系图如图一所示。
表二图一:承德中密砂lg (Ei/Pa) 与lg(σ3/Pa)的试验关系图对图一中的试验点进行拟合,得到lg (Ei/Pa) 与lg(σ3/Pa)的直线关系:y=0.8033x+2.2914.根据公式:E i=Kp a(σ3 p a )n可知K、n分别代表lg (Ei/Pa) 与lg(σ3/Pa)直线的截距和斜率,故可得K=2.2914;n=0.8033。
E-ν法在常规三轴试验中,轴向应变ε1与侧向应变—ε3之间也存在双曲线关系,经变换之后可得如下公式:−ε3ε1=f−Dε3由上式知—ε3/ε1与—ε3为直线关系,但实际上,二者并不是严格的直线关系,需先对试验结果进行取舍,然后选取某一区间进行拟合。
本文中选取试验曲线的后半部分进行拟合,得到不同围压下相应的拟合曲线,如下图所示。
图二:—ε3/ε1与—ε3关系曲线对应不同围压下的拟合曲线分别为:σ3=100kpa时,y=2.8211x+0.4719;σ3=300kpa时,y=2.8809x+0.4381;σ3=500kpa时,y=3.258x+0.4177.f和D分别为—ε3/ε1与—ε3直线的截距和斜率,结果如下表所示。
又因为νi=f=G -Flg (σ3/Pa )故可做νi—lg (σ3/Pa )关系曲线如下所示。
文章编号:1007 2284(2010)03 0076 04邓肯 张非线性模型研究及其在ANSYS 中的实现宿 辉1,2,党承华2,崔佳佳2(1.北京科技大学土木与环境学院,北京100083;2.河北工程大学水电学院,河北邯郸056021)摘 要:对工程领域使用广泛的邓肯 张非线性本构模型进行了研究,总结了国内外的研究现状及理论成果,针对其无法判定因结压力降低时的加载情况,提出了相应的变形模量的计算方法,同时考虑中主应力、土体抗拉强度的影响对模型进行了修正。
利用AP DL 编写程序实现了A N SYS 的材料本构模型的二次开发,运用重启动方法实现单元应力修正后数据重写入数据库,通过试验模拟对比分析验证了模型的适用性。
关键词:邓肯 张模型;抗拉强度;中主应力;应力分析 中图分类号:T U 470+.3 文献标识码:ADuncan Chang Nonlinear Elastic Model and Realization in ANSYSSU Hui 1,2,DANG Cheng hua 2,CUI Jia jia2(1.Schoo l of Civ il and Env iro nmen Engineering ,Beijing U niver sity o f Science and T echnolog y,Beijing 100083,China;2.Scho ol o f Water Resources and H ydropow er,Hebei U niv ersity of Engineering ,H andan 056021,H ebei Pro vince ,China)Abstract:Based on the pr esent research situation and theor etical r esults,research is do ne o n Duncan Chang no nlinear elastic mo del,w hich is applied w idely in engineer ing,.A cco rding t o the fact that Duncan Chang model can't judg e the lo ading situatio n w hile co n solidat ion pressure decr eases .T he model is mo dified cosidering the effect of int ermediate principal stress and tensile strength of so il.T he seco ndar y dev elo pment of mater ial co nstitut ive model in A N SYS is accomplished by utilizing the A PDL lang uag e.T hen t he r e start method is used to modify element str ess,the co rrected data is database rew ritten.A ser ies of compliance tests verify the accur a cy and applicability of the modified Duncan Chang nonlinear elastic mo del embedded AN SY S.Key words:Duncan Chang no nlinear elastic mo del;tensile str eng th;inter mediat e pr incipal stress;str ess analy sis 收稿日期:2009 05 18作者简介:宿 辉(1972 ),男,副教授,博士研究生,主要从事岩土工程及水工结构教学及研究工作。
有限元编程算例(Fortran)本程序通过Fortran语言编写,程序在Intel Parallel Studio XE 2013 with VS2013中成功运行,程序为《计算力学》(龙述尧等编)一书中的源程序,仅作研究学习使用,省去了敲写的麻烦.源程序为:!Page149COMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)OPEN(5,FILE=’DATAIN’)!OPEN(6,FILE=’DATAOUT’,STATUS=’NEW')CALL DATAIF(IND.EQ.0)GOTO 10EO=EO/(1.0—UN*UN)UN=UN/(1。
0—UN)10 CALL TOTSTICALL LOADCALL SUPPORCALL SOLVEQCALL STRESSPAUSE!STOPENDSUBROUTINE DATACOMMON/X1/NJ,NE,NZ,NDD,NPJ,IND,NJ2,EO,UN,GAMA,TE,AECOMMON/X2/JM(100,3),NZC(50),CJZ(100,2),PJ(100,2),B(3,6),D(3,3),S(3,6),TKZ(200,20),EKE(6,6),P(200)READ(5,*)NJ,NE,NZ,NDD,NPJ,INDNJ2=NJ*2NPJ1=NPJ+1READ(5,*)EO,UN,GAMA,TEREAD(5,*)((JM(I,J),J=1,3),I=1,NE)READ(5,*)((CJZ(I,J),J = 1,2),I=1,NJ)!Page150READ(5,*)(NZC(I),I=1,NZ)READ(5,*)((PJ(I,J),J=1,2),I=1,NPJ1)WRITE(6,10)(I,(CJZ(I,J),J=1,2),I=1,NJ)10 FORMA T(4X,2HNO,6X,1HX,6X,1HY/(I6,2X,F7。
ADINA软件中邓肯—张模型的二次开发与应用
陈志坚;江涛;陈松
【期刊名称】《水电能源科学》
【年(卷),期】2007(25)2
【摘要】根据AD INA提供的用户子程序,用FORTRAN编写了邓肯—张E-v本构模型,给出了模型开发的方法和步骤,并用常规三轴压缩试验对模型进行了检验。
检验结果表明,开发的模型正确合理、计算精度较理想,可应用于实际工程计算以解决复杂的非线性土工问题。
【总页数】3页(P65-67)
【关键词】ADINA;邓肯-张模型;二次开发
【作者】陈志坚;江涛;陈松
【作者单位】河海大学土木工程学院
【正文语种】中文
【中图分类】TU411;TP311.5
【相关文献】
1.邓肯-张E-B模型的ANSYS二次开发及应用 [J], 孙明权;陈姣姣;刘运红
2.Duncan-Chang E-υ模型在ADINA软件中的开发与应用 [J], 陈志波;朱俊高;陈浩锋
3.邓肯-张模型开发及其在面板坝计算中的应用 [J], 宋志宇;李斌;董莉莉;梁春雨
4.邓肯-张EB模型参数求解的二次优化法 [J], 陈立宏
5.一种改进邓肯张模型及其在土石坝数值模拟中的应用 [J], 邵东琛
因版权原因,仅展示原文概要,查看原文内容请购买。
应用MATLAB确定邓肯-张双曲线模型中的K,n参数简介:接合承德中密砂常规三轴试验数据,介绍应用Matlab语言编写计算及绘图程序来处理试验数据的方法,可显著提高试验研究的数据处理效率和结果的可视化程度。
关键字:Matlab 三轴试验邓肯-张模型1 前言基于广义胡克定律的线弹性理论形式简单,参数少,物理意义明确,而且在工程界有广泛深厚的基础,得以应用于许多工程领域中。
早期土力学中的变形计算主要是基于线弹性理论的,只有在计算机得到迅速发展之后,非线性理论模型才得到较广泛的应用。
邓肯-张模型是建立在增量广义胡克定律基础之上的变模量的弹性模型,可以反映土变形的非线性,并在一定程度上反映土变形的弹塑性,很容易为工程界所接受,加之所用参数和材料参数不多,物理意义明确,只需用常规三轴压缩试验即可确定这些参数及材料常数适应的土类比较广,所以该模型为岩土工程界所熟知,并得到了广泛的应用,成为土的最为普及的本构模型之一。
本文主要是应用MATLAB编写计算及绘图程序来处理承德中密砂常规三轴试验数据。
2 基于MATLAB的计算过程实现现场的观测数据经过采集和整理后,按照一定的格式把数据存储在数据文件中,然后可以使用MATLAB丰富的数值运算功能可以非常容易地编制出数据处理程序,先用函数fope n()打开数据文件,fid=fopen(‘filename’,’r’)再用fscanf 函数依次从文件中读取格式化数据来完成对各变量地赋值,其使用语法为:matrix=fscanf(fid,format)。
本文由于数据不是太多,所以在计算过程中没有采取调用存储文件地形式。
直接在计算过程中输入试验数据计算。
2.1 数据的处理对第一组数据,通过编写Matlab语言,由轴向应变和应力差的试验数据可以作出~()和~双曲线关系图形,主要用到的MATLAB命令为:plot(x1,y1);axis([0 0.04 0 3]) ;hold on%(1)plot(x1, x1./y1);a=polyfit(x1, x1./y1,1);t1=0:0.001:0.07;plot(x1, x1./y1,'.',t1,a(1)*t1 +a(2))%(2)其中x1代表第一组轴向应变,x2代表第一组应力差。
Fortran 讲解FORTRAN是科学计算语言。
很多年不碰它了(87年用)。
你怎么不去FOR TRAN论坛?试着说说:SUBROUTINE KE(IO,NE,NWE,T,A1,A2,V,EK,BCA) *子例行程序KE(可以CALL调用)DIMENSION B(7),BCA(7,NE),EK(6,6) *定义1个一维数组和2个两维数组DO 10 I=1,7 *循环,10是下面的标号B(I)=BCA(I,IO) *给一维数组B(I)赋值。
但BCA函数我没见过,外部函数?10 CONTINUE *未完成7次循环继续A=A1/B(7)*T *表达式计算结果赋值给ADO 20 I=1,3 *下面是双重循环——外循环DO 20 J=I,3 *内循环I1=2*IJ1=2*JEK(I1-1,J1-1)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3)) *2维数组通过表达式计算后赋值EK(I1-1,J1)=A*(V*B(I)*B(J+3)+A2*B(I+3)*B(J))EK(I1,J1-1)=A*(V*B(I+3)*B(J)+A2*B(I)*B(J+3))EK(I1,J1)=A*(B(I+3)*B(J+3)+A2*B(I)*B(J))20 CONTINUE *循环未完成继续DO 30 I=3,6 *有是一个双重循环DO 30 J=1,IEK(I,J)=EK(J,I)30 CONTINUE *循环未完成继续IF(NWE.EQ.0) GOTO 60 *如果NWE=0 转标号60处WRITE(6,40) IO *输出IO,6=显示器或打印机,40是表控格式,就是由40标号语句控制输出格式40 FORMAT(/1X,'EK NE='I5) *格式说明,似乎1X前多个/WRITE(6,50) EK *输出EK50 FORMAT(1X,6E11.4) *同上60 RETURN *返回操作系统END *程序结束EK(I1-1,J1-1)=A*(B(I)*B(J)+A2*B(I+3)*B(J+3))假如I=2,J=3I1=2*I=4,J1=2*J=6,这样上面的语句就相当于:EK(3,5)=A*(B(2)*B(3)+A2*B(5)*B(6))呵呵,就这意思。
邓肯-张EB模型参数求解的二次优化法陈立宏【摘要】邓肯-张非线性弹性模型是土石坝工程中最常用的本构模型.水利行业《土工试验规程》中根据应力水平75%和90%两点法进行计算时,得到的结果往往并不合理,有时n值还可能出现负数.一般的适线法仅仅对单个试样结果进行优化,而并不是针对整组试验结果,因此无法得到最优结果.提出了一种二步优化的参数计算方法,首先对每级围压下单个试样的试验成果采用适线法优化,得到每级围压下的参数a、b.在此基础上,计算得到参数K、n、Rf的初值.然后以邓肯-张理论为基础,根据获得的参数初值针对整组试验成果进行二次优化,以理论计算与试验的应力应变曲线差的平方和最小为目标函数,从而得到EB模型的主要参数.该方法简单实用,能够快速和准确地获得邓肯-张模型参数,并结合糯扎渡大坝堆石料三轴试验数据,对方法进行了验证.%Duncan-Chang nonlinear elastic constitutive model is the most used one in embankment dam engineering.The Specification of Soil Test in hydraulic industry proposes a computational method based on the values of two points from the stress-axial strain curve of the triaxial testing results.The stress levels of these two points are 75% and 90%respectively.However the proposed method cannot obtain reasonable results all the times,and sometimes even the parameter n maybe negative.Curve fitting methods make some progress,but still could not gain the optimal value for the parameters because these methods only based on single sample result.A two step optimization method for acquiring the optimal values of Duncan-Chang model is presented herein.First,the traditional curve fitting method is adopted to obtain thevalues of parameters a and b under each confining pressure.Then the parameters K,n and Rf are ing these parameters as initial values,a second optimization procedure is carried out to fit all the resultsof triaxial test to gain the parameters of Duncan-Chang model,in which,the minimum square sum of the differences of stress and strain curves of theoretical calculation and experiment is taken as the objectivefunction.The method is simple and practical,and can quickly and accurately obtain the parameters of DuncanZhang model.The method is validated based on the triaxial test data of Nuozhadu Dam.【期刊名称】《水力发电》【年(卷),期】2017(043)008【总页数】5页(P52-55,75)【关键词】堆石料;邓肯-张模型;优化方法;土石坝【作者】陈立宏【作者单位】北京交通大学土建学院,北京100044【正文语种】中文【中图分类】TU413堆石料作为高土石坝工程的主体填料,其工程特性和本构模型参数一直为大家所关注。
DIMENSION X(31),Y(31),DYCS(30,4),KB(93,94),DRTA(93)REAL KBINTEGER JDCS(31,4),TXFLAG(31),CSFLAGOPEN(1,FILE='E:\\INPUT.TXT',STATUS='OLD')OPEN(2,FILE='E:\\OUTPUT.TXT',STATUS='UNKNOWN')CALL INPUT(NJ,X,Y,DYCS,JDCS,RZ,E,TKH,TKV,QH,QV,HA,TXFLAG)NTRY=011 CSFLAG=0NTRY=NTRY+1CALL SUBKB(NJ,DYCS,JDCS,TXFLAG,X,Y,QH,QV,TKH,TKV,E,RZ,HA,KB)N=3*NJM=N+1CALL FCQJ(KB,N,M,DRTA)WRITE(2,100) NTRYWRITE(*,100) NTRY100 FORMAT(1X,I3,'STEP')WRITE(2,200)WRITE(*,200)200 FORMAT(1X,'NO.',5X,'X-DISP',5X,'Y-DISP',1X,'X-SPRING',1X,'R-DISP') DO 21 I=1,NJDISP=SQRT(DRTA(3*I-2)**2+DRTA(3*I-1)**2)WRITE(2,300) I,DRTA(3*I-2),DRTA(3*I-1),TXFLAG(I),DISPWRITE(*,300) I,DRTA(3*I-2),DRTA(3*I-1),TXFLAG(I),DISP21 CONTINUE300 FORMAT(1X,I2,2E11.3,I5,E11.3)DO 31 I=1,NJFLAG=DRTA(3*I-2)IF(FLAG.LE.0.0.AND.TXFLAG(I).EQ.0) THENTXFLAG(I)=1CSFLAG=1END IFIF(FLAG.GT.0.0.AND.TXFLAG(I).EQ.1) THENTXFLAG(I)=0CSFLAG=1END IF31 CONTINUEIF(CSFLAG.NE.0) GO TO 11CALL SUBFE(NJ,DYCS,E,DRTA)ENDSUBROUTINE INPUT(NJ,X,Y,DYCS,JDCS,RZ,E,TKH,TKV,QH,QV,HA,TXFLAG) DIMENSION X(31),Y(31),DYCS(30,4),JDHD(31)INTEGER JDCS(31,4),TXFLAG(31)REAL JDHDREAD(1,*) NJ,RZ,E,QH,QV,TKH,TKV,HADO 10 I=1,NJREAD(1,*) X(I),Y(I),(JDCS(I,J),J=1,4),JDHD(I)TXFLAG(I)=JDCS(I,1)10 CONTINUEDO 20 I=1,NJ-1DX=X(I+1)-X(I)DY=Y(I+1)-Y(I)DYCS(I,1)=SQRT(DX*DX+DY*DY)DYCS(I,2)=(JDHD(I+1)+JDHD(I))/2.0DYCS(I,3)=DY/DYCS(I,1)DYCS(I,4)=DX/DYCS(I,1)20 CONTINUEENDSUBROUTINE SUBKT(NJ,Y,TKH,KT)REAL Y(31),KT(31)DO 10 I=1,NJKT(I)=010 CONTINUEDO 20 I=2,NJ-1KT(I)=TKH*ABS(Y(I+1)-Y(I-1))/220 CONTINUEKT(1)=TKH*ABS(Y(2)-Y(1))/2KT(NJ)=TKH*ABS(Y(NJ)-Y(NJ-1))/2ENDSUBROUTINE SUBKC(HA,TKV,KC)REAL KC(3,3)DO 10 I=1,3DO 10 J=1,3KC(I,J)=0.010 CONTINUEKC(1,1)=0KC(2,2)=TKV*HAKC(3,3)=(TKV*HA**3)/12.0ENDSUBROUTINE JDHZ(NJ,X,Y,RZ,QH,QV,DYCS,HZ) DIMENSION HZ(31,3),X(31),Y(31),DYCS(30,4)DO 10 I=1,NJDO 10 J=1,3HZ(I,J)=010 CONTINUEDO 20 I=1,NJ-1PH=(Y(I+1)-Y(I))*QH/2HZ(I,1)=HZ(I,1)+PHHZ(I+1,1)=HZ(I+1,1)+PH20 CONTINUEDO 30 I=1,NJ-1IF(X(I).LT.X(I+1)) THENPV=-(X(I+1)-X(I))*QV/2HZ(I,2)=HZ(I,2)+PVHZ(I+1,2)=HZ(I+1,2)+PVEND IF30 CONTINUEDO 40 I=1,NJ-1G=-DYCS(I,1)*DYCS(I,2)*RZ/2HZ(I,2)=HZ(I,2)+GHZ(I+1,2)=HZ(I+1,2)+G40 CONTINUEENDSUBROUTINE TT(SI,CO,T)REAL T(6,6)DO 10 I=1,6DO 10 J=1,6T(I,J)=010 CONTINUET(1,1)=COT(1,2)=SIT(2,1)=-T(1,2)T(2,2)=T(1,1)T(3,3)=1.0DO 20 I=1,3DO 20 J=1,320 T(I+3,J+3)=T(I,J)RETURNENDSUBROUTINE KSE(E,BL,D,KE)REAL KE(6,6)DO 10 I=1,6DO 10 J=1,6KE(I,J)=010 CONTINUEKE(1,1)=E*D/BLKE(1,4)=-KE(1,1)KE(2,2)=E*D**3/BL**3KE(2,5)=-KE(2,2)KE(2,3)=E*D**3/2/BL**2KE(2,6)=KE(2,3)KE(3,3)=E*D**3/3/BLKE(3,5)=-KE(2,6)KE(3,6)=KE(3,3)/2KE(4,4)=KE(1,1)KE(5,5)=KE(2,2)KE(5,6)=KE(3,5)KE(6,6)=KE(3,3)DO 20 I=2,6DO 20 J=1,I-1KE(I,J)=KE(J,I)20 CONTINUEENDSUBROUTINE ESM(E,BL,D,SI,CO,KE)REAL KE(6,6),T(6,6),A(6,6),C(6,6)CALL TT(SI,CO,T)CALL KSE(E,BL,D,KE)DO 10 I=1,6DO 10 J=1,6A(I,J)=T(J,I)10 CONTINUEDO 20 I=1,6DO 20 J=1,6C(I,J)=0DO 20 K=1,6C(I,J)=C(I,J)+A(I,K)*KE(K,J)20 CONTINUEDO 30 I=1,6DO 30 J=1,6KE(I,J)=0DO 30 K=1,6KE(I,J)=KE(I,J)+C(I,K)*T(K,J)30 CONTINUEENDSUBROUTINE SUBKB(NJ,DYCS,JDCS,TXFLAG,X,Y,QH,QV,TKH,TKV,E,RZ,HA,KB) REAL KE(6,6),KC(3,3),KB(93,94),KT(31),HZ(31,3),DYCS(30,4),X(31),Y(34)INTEGER JDCS(31,4),TXFLAG(31)DO 10 I=1,3*NJDO 10 J=1,3*NJ+1KB(I,J)=010 CONTINUEDO 20 K=1,NJ-1BL=DYCS(K,1)HD=DYCS(K,2)SI=DYCS(K,3)CO=DYCS(K,4)CALL ESM(E,BL,HD,SI,CO,KE)DO 20 I=1,6DO 20 J=1,6KB(3*(K-1)+I,3*(K-1)+J)=KB(3*(K-1)+I,3*(K-1)+J)+KE(I,J) 20 CONTINUECALL SUBKT(NJ,Y,TKH,KT)DO 30 I=1,NJIF(TXFLAG(I).NE.0) THENKB(3*I-2,3*I-2)=KB(3*I-2,3*I-2)+KT(I)END IF30 CONTINUECALL SUBKC(HA,TKV,KC)DO 40 I=1,3KB(I,I)=KB(I,I)+KC(I,I)40 CONTINUECALL JDHZ(NJ,X,Y,RZ,QH,QV,DYCS,HZ)DO 50 I=1,NJDO 50 J=1,3KB(3*(I-1)+J,3*NJ+1)=KB(3*(I-1)+J,3*NJ+1)+HZ(I,J)50 CONTINUEDO 60 I=1,NJDO 60 K=2,4IF(JDCS(I,K).EQ.0) THENDO 70 J=1,3*NJ+170 KB(3*(I-1)+K-1,J)=0DO 80 L=1,3*NJ80 KB(L,3*(I-1)+K-1)=0KB(3*(I-1)+K-1,3*(I-1)+K-1)=1.0END IF60 CONTINUEENDSUBROUTINE FCQJ(A,N,N1,X)DIMENSION A(93,94),X(93)DO 50 K=1,ND=0.0DO 20 I=K,NIF (D-ABS(A(I,K))) 10,20,2010 D=ABS(A(I,K))L=I20 CONTINUEIF (L.EQ.K) GO TO 30DO 25 J=K,N1T=A(L,J)A(L,J)=A(K,J)25 A(K,J)=T30 T=1.0/A(K,K)K1=K+1DO 40 J=K1,N1A(K,J)=A(K,J)*TDO 40 I=K1,N40 A(I,J)=A(I,J)-A(I,K)*A(K,J)50 CONTINUEDO 60 IK=2,NI=N1-IKI1=I+1DO 60 J=I1,N60 A(I,N1)=A(I,N1)-A(I,J)*A(J,N1)DO 70 I=1,N70 X(I)=A(I,N1)ENDSUBROUTINE SUBFE(NJ,DYCS,E,D)DIMENSION DYCS(30,4),D(93),T(6,6),DD(6),KE(6,6),FF(6)REAL KEWRITE(2,100)100 FORMAT(1X,'NO',8X,'NI',8X,'QI',8X,'MI',8X,'NJ',8X,'QJ',8X,'MJ') DO 10 I=1,NJ-1BL=DYCS(I,1)HD=DYCS(I,2)SI=DYCS(I,3)CO=DYCS(I,4)CALL TT(SI,CO,T)DO 20 J=1,6DD(J)=0DO 20 K=1,6DD(J)=DD(J)+T(J,K)*D((I-1)*3+K)20 CONTINUECALL KSE(E,BL,HD,KE)DO 30 J=1,6FF(J)=0DO 30 K=1,6FF(J)=FF(J)+KE(J,K)*DD(K) 30 CONTINUEWRITE(2,200) I,(FF(J),J=1,6) 200 FORMAT(1X,I2,6F10.3)10 CONTINUEEND。
邓肯张模型FORTRAN子程序源代码
SUBROUTINE UMA T(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,
1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN,
2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MATERL,NDI,NSHR,NTENS,
3 NSTATV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT,
4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC)
C
INCLUDE 'ABA_PARAM.INC'
C
CHARACTER*80 MATERL
DIMENSION STRESS(NTENS),STATEV(NSTA TV),
1 DDSDDE(NTENS,NTENS),DDSDDT(NTENS),DRPLDE(NTENS),
2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),PREDEF(1),DPRED(1),
3 PROPS(NPROPS),COORDS(3),DROT(3,3),
4 DFGRD0(3,3),DFGRD1(3,3)
C
DIMENSION PS(3),DSTRESS(NTENS),TDSTRESS(NTENS),TSTRESS(NTENS) PARAMETER (ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,SIX=6.0D0)
K=PROPS(1)
N=PROPS(2)
RF=PROPS(3)
C=PROPS(4)
FAI=PROPS(5)/180.0*3.1415926
G=PROPS(6)
D=PROPS(7)
F=PROPS(8)
KUR=PROPS(9)
PA=PROPS(10)
DFAI=PROPS(11)/180.0*3.1415926
S1S3O=STATEV(1)
S3O=STATEV(2)
SSS=STATEV(3)
CALL GETPS(STRESS,PS,NTENS)
FAI=FAI-DFAI*LOG10(S3O/PA)
CALL GETEMOD(PS,K,N,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O,G,D,F
1 ,SSS,S1S3O)
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
TDSTRESS=0.0
CALL GETSTRESS(DDSDDE,TDSTRESS,DSTRAN,NTENS)
DO 701 I1=1,NTENS
TSTRESS(I1)=STRESS(I1)+TDSTRESS(I1)*0.5
701 CONTINUE
CALL GETPS(TSTRESS,PS,NTENS)
CALL GETEMOD(PS,K,N,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O,G,D,F,
1 SSS,S1S3O)
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
DSTRESS=0.0
CALL GETSTRESS(DDSDDE,DSTRESS,DSTRAN,NTENS)
DO 702 I1=1,NTENS
STRESS(I1)=STRESS(I1)+DSTRESS(I1)
702 CONTINUE
CALL GETPS(STRESS,PS,NTENS)
CALL GETEMOD(PS,K,N,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O,G,D,F,
1 SSS,S1S3O)
EBULK3=EMOD/(ONE-TWO*ENU)
EG2=EMOD/(ONE+ENU)
EG=EG2/TWO
EG3=THREE*EG
ELAM=(EBULK3-EG2)/THREE
CALL GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
IF(PS(3).GT.S3O)S3O=PS(3)
IF((PS(1)-PS(3)).GT.S1S3O)S1S3O=PS(1)-PS(3)
IF(S.GT.SSS)SSS=S
STA TEV(1)=S1S3O
STA TEV(2)=S3O
STA TEV(3)=SSS
END
SUBROUTINE GETPS(STRESS,PS,NTENS)
INCLUDE 'ABA_PARAM.INC'
DIMENSION PS(3),STRESS(NTENS)
CALL SPRINC(STRESS,PS,1,3,3)
DO 310 I=1,2
DO 320 J=I+1,3
IF(PS(I).GT.PS(J))THEN
PPS=PS(I)
PS(I)=PS(J)
PS(J)=PPS
END IF
320 CONTINUE
310 CONTINUE
DO 330 K1=1,3
PS(K1)=-PS(K1)
330 CONTINUE
RETURN
END
SUBROUTINE GETEMOD(PS,K,EN,RF,C,FAI,ENU,PA,KUR,EMOD,S,S3O
1 ,G,D,F,SSS,S1S3O)
INCLUDE 'ABA_PARAM.INC'
DIMENSION PS(3)
S=(1-SIN(FAI))*(PS(1)-PS(3))
IF(PS(3).LT.0) THEN
PSFEI=0.1
ELSE
PSFEI=PS(3)
END IF
S=S/(2*C*COS(FAI)+2*PSFEI*SIN(FAI))
IF(S.GE.0.99) THEN
S=0.99
END IF
AA=D*(PS(1)-PS(3))
AA=AA/(K*PA*((S3O/PA)**N))
AA=AA/(1-RF*S)
ENU=G-F*LOG10(S3O/PA)
ENU=ENU/(1-AA)/(1-AA)
IF(ENU.GT.0.49)ENU=0.49
EMOD=K*PA*((S3O/PA)**N)*((1-RF*S)**2)
IF(S.LT.SSS.AND.(PS(1)-PS(3)).LT.S1S3O)THEN
EMOD=KUR*PA*((S3O/PA)**N)
END IF
END
SUBROUTINE GETFEI(STRESS,SMISES,NTENS,NDI)
INCLUDE 'ABA_PARAM.INC'
DIMENSION STRESS(NTENS)
SMISES=STRESS(1)**2+STRESS(2)**2+STRESS(3)**2+
1STRESS(4)**2+STRESS(5)**2+STRESS(6)**2
SMISES=SQRT(SMISES)
RETURN
END
SUBROUTINE GETDDSDDE(DDSDDE,NTENS,NDI,ELAM,EG2,EG)
INCLUDE 'ABA_PARAM.INC'
DIMENSION DDSDDE(NTENS,NTENS)
DO 20 K1=1,NTENS
DO 10 K2=1,NTENS
DDSDDE(K2,K1)=0.0
10 CONTINUE
20 CONTINUE
DO 40 K1=1,NDI
DO 30 K2=1,NDI
DDSDDE(K2,K1)=ELAM
30 CONTINUE
DDSDDE(K1,K1)=EG2+ELAM
40 CONTINUE
DO 50 K1=NDI+1,NTENS
DDSDDE(K1,K1)=EG
50 CONTINUE
RETURN
END
SUBROUTINE GETSTRESS(DDSDDE,STRESS,DSTRAN,NTENS)
INCLUDE 'ABA_PARAM.INC'
DIMENSION DDSDDE(NTENS,NTENS),STRESS(NTENS),DSTRAN(NTENS) 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。