邓肯张模型FORTRAN子程序源代码

合集下载

邓肯-张模型参数求取

（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 ）关系曲线如下所示。

Fortran子程序.ppt

PROGRAM FLOWER DO I=100, 999 IF (NUM(I,100)**3+NUM(I,10)**3+NUM(I,1)**3==I) THEN
WRITE(*,*) I ENDIF ENDDO END
FUNCTION NUM(N, I)
SELECT CASE(I)
CASE(100)
第八章子程序
语句函数函数子程序子例行程序程序单元之间的数据传递递归调用数据共用存储单元与数据块子程序内部子程序模块
8.1 语句函数语句函数——
PROGRAM EXAM1A X=1 FX=5*X**3-2*X**2+7*X+6
用一个语句定义的函数。WRITE(*,*) 'f(', X, ')=', FX
PARAMETER (N=20)
DIMENSION A(N)
INTEGER A
A(1)=17
！该语句及随后的三个语句生成数组A
DO I=2, 20
A(I)=MOD(19*A(I-1), 1024)
ENDDO
DO I=1, N-1 ！二重循环完成对A的排序
DO J=I+1, N
CALL SWAP(A(I), A(J))
当未进行说明时，函数值的类型遵守I-N规则。虚参的类型在函数体中说明，否则遵守I-N规则。
二、函数子程序的调用
P130 例8.4 用函数子程序的方法设计一个程序，求50-100内的所有素数及其和。
分析：设计一个函数子程序 PRIME(N)，函数 PRIME的值定义如下：
1 当n是素数时 prim e(n) 0 当n非素数时
注意：
函数用一条语句可以定义时，才用语句函数的形式定义函数。

邓肯_张非线性模型研究及其在ANSYS中的实现

文章编号: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中的子程序

一、子例行程序的定义
子例行程序是以Subroutine 语句开头，并以End语句结束的一个程序段，其定义的一般格式为：
Subroutine 子例行程序名（虚参表）
子例行程序体 end
2015/12/20 12
注意：（1）子例行程序的命名方法与变量相同。虚参由变量、数组名（不能是数组元素，常数、表达式）充当，当有多个虚参时，之间用“ ,” 分隔，没有虚参时子例行程序名后面的一对括号可以省略；（2）子例行程序的设计方法与函数子程序相同，但不能有对子例行程序名的赋值语句（因为其名字没有值）。
第六讲 Fortran中的子程序
实际中的程序由若干个程序单元组成，但是有且只有一个主程序（ Program main ），其它的都是子程序。子程序是构造大型实用程序的有效工具，设计程序要善于利用子程序，因此，本讲学习Fortran中的子程序：函数子程序和子例行程序。此外，在Fortran中还有一种类似与函数子程序的语句函数。
2015/12/20 15
分析：牛顿迭代公式为：
program newton f(x)=7*x**4+6*x**3-5*x**2+4*x+3 f ( xn ) df(x)=28*x**3+18*x**2-10*x+4 xn 1 xn ' f ( xn ) read*,x0,e !x0=-1.0,e=0.0001 x1=x0 首先要给个迭代初值 x0 ， x2=x1-f(x1)/df(x1) 代入公式的右边计算 x1 ， do while(abs(x1-x2)>e) 再把 x1 代入计算出 x2 ， … ， x1=x2 x2=x1-f(x1)/df(x1) 当第n次和第(n+1)次计算 end do 结果差不多时（在规定的 print*,x1,f(x1) 精度内），则计算结果就 end

本构模型之邓肯张模型

3
G , F 为试验常数，其确定见图(c)。
将式（16）微分，得：
d v 1 D ) f D f 3 ( i 1 1 v t 2 2 d ( 1 D ) ( 1 D ) 1 1 1
(19)
将式（4）、式（8）、式（10）和式（18）代入式（19），则得到任一应力（σ1，σ3）时的泊松比的邓肯－张计算公式：
R f 值一般在0.75~1.0之间
R 1 f b ( ) ( )f 1 3 u l t 1 2
(8)
将式（8）、式（4）代入式（3）中，得
1 1 Et Rf Ei 1 1 E ( ) 1 3 f i
t
(4)
这表明代表的是在这个试验中的起始变形模量(初始切线模量)的倒数。
在式（1）中，如果 1 ，则：
(1 3 )ult 1 b
(5) (6)
或者 :
1 b (1 3 )ult
由此可看出b代表的是双曲线的渐进线所对应的极限偏差应力 (1 3)ult 的倒数。
在试验中 1 不可能无限大，求取 (1 3)ult ；对于有峰值 ) ( ) 3 ) ( 3 ) 点的情况，取 ( 1 3 f 1 3 1 f 1 u l t 峰，这样 ( 定义破坏比 R f 为： (1 3 ) f Rf (7) (1 3 )ult
本构模型之邓肯张模型
主要是根据试验成果拟合推导得出
邓肯-张双曲线模型
• 该模型是一种建立在增量广义虎克定律基础上的非线性弹性模型，可经反映应力～应变关系的非线性，模型参数只有８个，且物理意义明确，易于掌握，并可通过静三轴试验全部确定，便于在数值计算中运用，因而，得到了广泛地应用。

有限元编程算例(fortran)【范本模板】

有限元编程算例（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。

ABAQUS二次开发Duncan-Chang模型在粗粒料三轴试验中的应用

2004 级江涛硕士的毕业论文《大直径超长钻孔灌注桩群桩效应的三维非线性有限元分析》中使用的 2 次开发模型，模型可靠性已得到验证[2]。
3.粗粒料三轴实验
粗粒料通常是指块石、碎石(或砾卵石)、石屑、石粉等粗颗粒组成的无粘性
2

混合料，或是粘性土中含有大量粗颗粒的混合土。《土工试验规程》（SL237－
含量(％)
46.2
38.4 15.4
试样直径 101mm，高度 200mm，采用击实制样，控制干密度取 ρ = 1.9g / cm3 。
进行 3 组常规三轴试验，围压分别为 200kPa、500kPa、800kPa，采用固结排水
剪，剪切速率取 0.001mm/s，固结和剪切过程均按《土工试验方法标准》（GB/T
Duncan-Chang 模型是目前岩土工程方面运用十分广泛的岩土体本构模型，其能够很好地反映土体非线性变形特征。其参数具有比较明确的物理意义，可教容易由常规三轴试验得到。本文将用 2 次开发邓肯—张模型来模拟高围压三轴排水实验。
2. Duncan-Chang 模型
1970 年Duncan、Chang提出了邓肯E-V双曲线非线性弹性模型，该模型属于变弹性模型，其理论基础是基于广义虎克定律基础上的增量线弹性模型，这种模型假定：土体变形虽然是弹性的，但是在微小变形时，变形可以看作是线性的，并且服从增量线性的各向同性广义虎克定律。其基本假设依然保留了弹性理论中
的斜率或截距得到。在拟合过程中，得到的参数是综合 3 组试验数据的结果，所
4

以高围压 0.8MPa 下模拟的曲线大部分在试验曲线以下，而围压 0.2MPa、0.5MPa 的情况下模拟曲线大部分在试验曲线以上。如果试验的过程严格安照规范的规定，操作过程在没有失误的情况下，排除了人为因素，一定能得到很好的模拟效果。

fortran子程序【精选】

parameter (n1=3,n2=3,n3=3)
dimension a(n1,n2)
dimension b(n2,n3),c(n1,n3)
open(5,file='input.dat')
call getmat(a,n1,n2)
call getmat(b,n2,n3)
call matpro(a,n1,n2,b,n3,c)
20 continue end
real w(5,5) call readin(w)
输入矩阵 1,2,3,4,5 2,3,4,5,1
call opp(w,x1,x2) write(*,100) x1,x2 100 format(1x,'The two sum of',
3,4,5,1,2
$ ' oppusite angles elements:',
subroutine sub(ch) charact...er*(*) ch end
11
（4）如果实参是变量或数组元素，在调用子程序时，对应的虚实参数实际上将共用同一存储单元。
program main
subroutine sub(x,a)
integer a,c(3)
integer x,a
data c/3*0/
4,5,1,2,3 5,1,2,3,4
$ /1x,'x1=',f8.2,' x2=',f8.2) end
9
主要区别：
1. 形式不一样（function、subroutine、call），无虚参时括号使用不一样。
2. 子程序名的意义不一样。函数子程序名代表函数值，通过其传递数据，需要类型说明；而子例行程序名仅为了调用使用，通过虚实参传递数据。

ADINA软件中邓肯—张模型的二次开发与应用

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确定邓肯－张双曲线模型中的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代表第一组应力差。

ADAMS FORTRAN程序

Fortran subroutine to run in ADAMS/Solver FORTRAN
1
Fortran subroutine to run in ADAMS/Solver FORTRAN
! PARAMETER (AROD=3.8D−4) PARAMETER (ATUBE=8.04D−4) PARAMETER (MRAROD=78D−6) PARAMETER (TEM=301) PARAMETER (VSTAT=89D−6) PARAMETER (VISCO=24) ! ! Identified parameters: ! ! Rebound parameters ! PARAMETER (FOR=80) PARAMETER (FPR=510) PARAMETER (KLR=8200) PARAMETER (KHR=7) PARAMETER (KPR=120) PARAMETER (KSR=840) PARAMETER (KTR=3.5) ! ! Compression parameters ! PARAMETER (FOC=95) PARAMETER (FPC=265) PARAMETER (KLC=1400) PARAMETER (KHC=1) PARAMETER (KPC=23) PARAMETER (KSC=640) PARAMETER (KTC=2.2) ! ! Defining lower damper mount ! MARKER/01, QP=267,600,180 ! ! Defining upper damper mount ! MARKER/02, QP=307,500,680 ! ! Determine displacement ! DIS, FUNCTION=DM(01,02) ! ! Determine velocity ! VEL, FUNCTION=VM(01,02) ! ! Determine acceleration ! ACC, FUNCTION=ACCM(01,02) ! !$$$$$$$$$$$$$$$$$$$$$$$$ EXECUTABLE CODE $$$$$$$$$$$$$$$$$$$$$$$$$$$

Fortran 讲解-平面有限元源程序

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))呵呵，就这意思。

在ABAQUS中开发实现Duncan_Chang本构模型

降低的影响。
1034
岩土力学
2004 年
对于卸载情况的处理方法，不同文献所考虑的准则有所不同，但卸载弹性模量计算公式应统一表达为
Eur
= kur Pa
⎜⎛ ⎝
−σ 1 Pa
⎟⎞ nur ⎠
(8)
至此，由式 (2) ～ (8) 可以全部确定 DuncanChang E-B 模型刚度系数矩阵的元素，然后，可按照三维问题或平面应变问题的数学表达式形成相应的刚度系数矩阵 [D] ，如对于 Duncan-Chang E-B 三维模型问题有
修改稿收到日期：2003-09-04
作者简介：徐远杰，男，1956 年生，教授，结构破坏与数值仿真研究所所长，从事岩土工程数值仿真研究。
第7期
徐远杰等：在 ABAQUS 中开发实现 Duncan-Chang 本构模型
1033
一不足，本文作者利用 ABAQUS 提供的二次开发用户子程序，完成了 Duncan-Chang 本构模型[1]的接口开发工作，并且对若干模型算例进行了测试，计算结果表明，模型算例的求解效率和精度令人满意，有望成为强大的解决大型复杂土工问题的非线性数值分析工具。同时，开发与维护的工作量将大大降低，是解决三维复杂土工数值分析问题的快捷方案和有效途径。
XU Yuan-jie, WANG Guan-qi, LI Jian, TANG Bi-hua
（School of Civil and Architectural Engineering, Wuhan University, Wuhan 430072, China）
Abstract: The software of ABAQUS, which is developed and supported by HKS Inc., is currently one of the most powerful nonlinear finite element analysis implement in the world. It includes many different kinds of material models. However, the release of ABAQUS currently absents Duncan-Chang constitutive model that has been widely used to simulate some deformation problems in the field of geotechnical engineering. The shortage limits the software’s application scopes in some extents. In this paper, the procedure and outline to develop Duncan-Chang constitutive model in ABAQUS by using User Subroutine were proposed and implemented. Two numerical examples were also solved and compared with conventional triaxial compression test data. All numerical results exhibit that it will be possible to accomplish some complex stress-strain analysis for soil engineering problems by fully making use of the powerful solving platform after ABAQUS is added Duncan-Chang constitutive model. Meanwhile, there are some advantages of faster convergence, higher efficiency, higher precision and more conveniently preprocessing-postprocessing. The programming difficulties to develop soil engineering and workloads of maintenance are also greatly decreased. Key words: finite element; Duncan-Chang constitutive model; ABAQUS; redevelopment

邓肯-张EB模型参数求解的二次优化法

邓肯-张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堆石料作为高土石坝工程的主体填料,其工程特性和本构模型参数一直为大家所关注。

fortran程序坐标转换源代码

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。