基于FISH语言的邓肯-张模型在FLAC3D中的开发与实现
基于FISH语言的邓肯-张模型在FLAC3D中的开发与
实现
尚波，刘忆，周建波
河海大学水利水电学院，南京（210098）
E-mail：sdg414@https://www.wendangku.net/doc/4f16473172.html,
摘要：根据FLAC3D软件中提供的二次开发平台，利用FISH语言在FLAC3D软件中实现了邓肯-张模型的开发，并提供一算例对所开发的模型的正确性进行了验证。
关键词：邓肯-张模型，FLAC3D，FISH语言，二次开发
1. 引言
1986年，美国ITASCA咨询集团公司研制开发了FLAC(Flac Lagrangian Analysis of Continua) 程序，它基于联系介质显示拉格朗日差分法，适合求解非线性的大变形问题，FLAC3D是三维有限差分程序。在国外，FLAC/FLAC3D软件在众多行业中得到了较大的应用[1]，如工程地质学、构造地质学、大陆动力学、成矿学等。上世纪90年代初该软件被引入中国，主要应用于岩土力学分析，例如矿体滑坡、煤矿开采沉陷预测、水利枢纽掩体稳定性分析、采矿巷道稳定性研究等。目前FLAC/FLAC3D软件已广泛应用于岩土工程、采矿工程、水利工程、地质工程等行业和部门。
FLAC3D具有强大的适合模拟岩土材料的本构模型及结构模型，FLAC提供了十种基本的本构模型，分别归类到零模型、弹性模型和塑性模型之中。其中零模型通常用来表示被移除或开挖掉的材料。弹性模型包括各向同性弹性模型和横观各向同性弹性模型。除此，还包括：摩尔-库仑模型、德鲁克-普拉格模型、节理化模型、应变硬化-软化模型、双线性应变硬化-软化的节理化模型、双屈服模型、修正的剑桥粘土模型以及霍克-布朗模型。但国内岩土工程领域中广泛采用的邓肯-张本构模型没有在FLAC3D中实现，这使得FLAC3D在国内岩土数值分析中的应用受到了限制。为了弥补不足，本文利用FLAC3D提供的二次开发平台，用FISH语言成功的实现了邓肯-张模型在FLAC3D软件中的开发，并通过一算例对模型的正确性进行了验证。
2. FLAC3D计算步骤及FISH语言
2.1 计算步骤
图1表明了FLAC所包含的一般计算过程。这个过程首先调用运动方程从应力和外力导出了新的速度和位移。然后，根据速度导出应变速率，以及有应变速率得出新的应力。对于循环的每一周期，我们采用一个时步。值得注意的是图1中的每个方框根据已知值更新另外自身的网格变量，而这些已知值在方框内操作时是保持恒定的。例如，下部的方框取一组已算出的速度值，对每个单元计算新的应力。方框内应力-应变关系运算是可以将速度值假设为定值，也就是说，最新计算的应力并不影响速度。由于我们所选的时步很小，则信息在如此小的时间间隔里就不会从一个单元传到另一个单元，信息在所以物质中传播速度都有个极限。由于每个循环圈占用一个小的时步，相邻单元在计算过程中的确不能互相影响，那么关于速度恒定的假定也就是合理的了[2]。
图1 基本显式计算循环
2.2 FISH 语言
FISH 是FALC 中自带的编程语言，被用户用于定义变量和新建函数，从而拓展了FLAC 的用途，增加了用户自定义功能。FISH 可以用于编写新的本构模型。自定义模型中可拥有类似FLAC 内置参数一样的用户自己命名的特性参数，可以用命令PROPERITY 设置，PRINT 和PLOT 命令打印输出参数。用户自定义的模型的运行速度比内置模型要慢。模型优化后，正常的运行速度为内置模型的1/4到1/3。自定义模型虽增加了运行时间，却节省了人力，克服了内置模型的不适用问题[2]。
3. 邓肯-张模型的数学表达式
对于邓肯-张E-B 模型：切线弹性模量
2
3123(1sin )()1()2cos 2sin n t f a
a
E R kp c p σ?σσ?σ???
??=???+?? （1）
切线体积模量
3
(
)m
t b a a
K k p p σ= （2）
由于v 只能在0 ~ 0.49之间变化，K 要有所限制，须限制在（0.37 ~ 17）t E 之间。
对于卸载情况，回弹模量的表达式为
3
(
)n
ur ur a a
E k p p σ= （3）
同时，为了考虑颗粒体材料的抗剪强度随围压增大而降低，内摩擦角改用为
1
0log(
)
a p σ???=?V （4）
式（1），（2），（3）和（4）中，c 、?、?V
、f
R 、k 、n 、
b k 、m 、ur k 均由
常规三轴实验得出[3]。
4. 程序流程
用FISH 编制的邓肯-张模型的程序流程[4]见图2。
图2 程序流程图
5. 算例验证
为了对所开发的邓肯-张模型的正确性进行验证，本文提供了一个简单的算例。为简便起见，该模型为一均质土石坝，覆盖层厚度为20米，长410米，宽120米。坝体高度为50米，坝顶宽10米，上下游坝坡为1：2，如图3。本文只考虑施工期，施工期荷载分期加载，坝体共分十级加载，考察其竖向沉降和水平向位移。其基本材料参数见表1。计算结果：竖向位移见图4，水平位移见图5。
表1 材料参数
图3 坝体网格
图4 竖向位移
从图4坝体竖向沉降的规律看，竖向最大位移在坝体的高度（1/3 ~ 1/2）附近，图5水平位移左右对称，与一般工程实践相符合，因此，证明了模型开发的合理性。
6. 结论
根据FLAC3D软件提供的二次开发平台，成功的实现了邓肯-张本构模型的开发。本文的主要目的是为了验证该模型的正确性，所以设计了一个简单的算例进行数值模拟计算。结
果表明，模型的计算结果比较合理。
参考文献
[1] 龚纪文，席先武，王岳军等．应力与变形的数值模型方法——数值模拟软件 FLAC 介绍[J]，华东地质学院学报，2002，25（3）： 220-227．
[2] 刘波，韩彦辉（美国）．《FLAC原理、实例与应用指南》[M]，北京：人们交通出版社，2005.9．
[3] 钱家欢，殷宗泽．《土工原理与计算》[M]，北京：中国水利水电出版社，1996 ：54-92．
[4] 花加凤．土石坝膜防渗结构问题探讨[D]，河海大学工学硕士学位论文，2006．
Development and implementation of
Duncan-Changconstitutive model based with FISH language
in FLAC3D
Shang Bo，Liu Yi，Zhou Jianbo
Institute of Water Conservancy and Hydropower Engineering，Hohai University，Nanjing
（210098）
Abstract
According to the further developing platform, the Duncan-Chang constitutive model is developed in FLAC3D with FISH language. The provided example is performed to verify the correctness of the defined model.
Keywords：Duncan-Chang constitutive model，FLAC3D，FISH language，further development
邓肯-张模型研究认识
塑性力学读书报告 邓肯-张模型研究认识 学院：建设工程 姓名：王吉亮 学号：2006631011 专业：地质工程
教师：金英玉
邓肯-张模型研究认识 王吉亮（83分） 摘 要:从邓肯-张模型的本源开始，分析研究了邓肯-张模型与E-B 模型的建立过程和模型中参数如何确定的问题，结合对该模型的认识，提出该模型具有的缺点与不足。 关键词:邓肯-张模型；E-B 模型；参数确定 CONGNITION ON THE STUDY OF DUNCAN-CHANG MODEL Wang Jiliang Abstract: rom the parent of Duncan-Chang model, studing the establish procedure of Duncan-Chang model and E-B model, introducing the problem of how to define the indexes in the model. Associate the congnition on this model, present the shortcomings. Keywords: Duncan-Chang model; E-B model; indexes define 1 引言 邓肯－张模型是一个非线性本构模型，既然是一个本构模型，可想而之他反应的是应力与应变之间的关系。说它是非线性的，那么反映应力应变关系的模量就不是一个常数Ｅ那么简单。在介绍该模型之前，先要介绍一个概念，就是反映非线性关系的增量广义胡克定律： 1123()t t t v d d d d E E σεσσ= -+ （1） 1963年，康纳（Kondner ）根据大量土的三 轴试验的应力应变关系曲线，提出可以用双曲线拟合出一般土的三轴试验13()~a σσε-曲线，即： 13a a a b εσσε-= + （2） 其中，a 、b 为试验常数。对于常规三轴压缩试验，1a εε=。邓肯等人根据这一双曲线应力应变关系提出了一种目前被广泛的增量弹性模型， 一般被称为邓肯－张（Duncan-Chang ）模型。 在常规三轴压缩试验中，13a a a b εσσε-=+可以写成： 1113 a b εεσσ=+- （3） 将常规三轴压缩试验的结果按 11 13 ~εεσσ-的关系进行整理，则二者近似成线性关系(见图1)。其中，a 为直线的截距；b 为直线的斜率。 在常规三轴压缩试验中，由于 230d d σσ==，所以切线模量为 ε1/(σ1 -σ3 ) -σ3 )ult 图1 1113 ~εεσσ-线性关系图 132 11()() t d a E d a b σσεε-= =+ （4） 在试验的起始点，10ε=，t i E E =，则： 1 i E a = ，这表明a 表示的是在这个试验中的起始变形模量E i 的倒数。如果1ε→∞，则： 131 ()ult b σσ-= （5）
邓肯张模型参数(精)
5.370569 应力差(б1-б3/100kPa轴向变形ε1 体应变εv ε3ε1/(б1-б3 0.5080.002250.00074-0.0007554.42913E-051.0020.004490.0013- 0.0015954.48104E-051.4630.006740.00176-0.002494.60697E-051.8490.008980.00223-0.0033754.85668E-052.1490.011230.00223-0.00455.22569E-052.3310.013480.00214-0.005675.78293E-052.4770.015720.00176-0.006986.34639E-052.5870.017970.00158-0.0081956.94627E-052.6650.020210.00111-0.009557.58349E-052.730.022460.00056-0.010958.22711E-052.7770.024710.00009-0.012318.89809E-052.8120.026950.00003-0.013469.58393E-052.8450.03032-0.00012-0.015220.0001065732.8780.03369-0.00195-0.017820.000117062.8890.03706-0.00288-0.019970.000128282.8940.04043-0.00381-0.022120.0001397032.8890.04492-0.00474-0.024830.0001554862.8790.04941-0.00567-0.027540.0001716222.8690.0539-0.0065-0.03020.000187872.8530.05839-0.00752- 0.0329550.0002046622.8360.06289-0.00827-0.035580.0002217562.809 0.06738 -0.00891 -0.038145 0.000239872 应力差(б1-б3/100kPa轴向变形ε1 体应变εv ε3ε1/(б1-б3 0.9090.001250.00075-0.000250.1375137512.2860.00350.00151- 0.0009950.1531058623.5330.005750.00236-0.0016950.1627512034.3410.0080.00302-
邓肯-张模型公式推导 高土
邓肯－张模型是一个非线性本构模型，既然是一个本构模型，可想而之他反应的是应力与应变之间的关系。说它是非线性的，那么反映应力应变关系的模量就不是一个常数Ｅ那么简单。在介绍该模型之前，先要介绍一个概念，就是反映非线性关系的增量广义胡克定律： 1123()t t t v d d d d E E σεσσ= -+ （1） 1963年，康纳（Kondner ）根据大量土的三 轴试验的应力应变关系曲线，提出可以用双曲线拟合出一般土的三轴试验13()~a σσε-曲线，即： 13a a a b εσσε-= + （2） 其中，a 、b 为试验常数。对于常规三轴压缩试验，1a εε=。邓肯等人根据这一双曲线应力应变关系提出了一种目前被广泛的增量弹性模型，一般被称为邓肯－张（Duncan-Chang ）模型。 在常规三轴压缩试验中，13a a a b εσσε-=+可以写成： 1113 a b εεσσ=+- （3） 将常规三轴压缩试验的结果按 11 13 ~εεσσ-的关系进行整理，则二者近似成线性关系(见图 1)。其中，a 为直线的截距；b 为直线的斜率。 在常规三轴压缩试验中，由于 230d d σσ==，所以切线模量为 ε1 /(σ1 -σ3 ) 1 b=1/(σ1 -σ3 )ult a =1/E i 图1 1113 ~εεσσ-线性关系图 132 11()() t d a E d a b σσεε-= =+ （4） 在试验的起始点，10ε=，t i E E =，则： 1i E a = ，这表明a 表示的是在这个试验中的起始变形模量E i 的倒数。如果1ε→∞，则： 131 ()ult b σσ-= （5） 由此可以看出b 代表的是双曲线的渐近线所对应的极限偏差应力13()ult σσ-的倒数。 在土的试样中，如果应力应变曲线近似于双曲线关系，则往往是根据一定的应变值（如 115%ε=）来确定土的强度13()f σσ-，而不可 能在试验中使1ε无限大，求取13()ult σσ-；对于有峰值点的情况，取1313()()f σσσσ-=-峰， 这样1313()()f σσσσ--ult <。定义破坏比R f 为： 1313()()f f R σσσσ-=-ult （6） 而 13131 ()()f f R b σσσσ== --ult （7） 将上式与1 i E a = 代入 132 11()() t d a E d a b σσεε-==+ （8） 得到：
邓肯张模型FORTRAN子程序源代码
邓肯张模型FORTRAN子程序源代码 SUBROUTINE UMA T(STRESS,STA TEV,DDSDDE,SSE,SPD,SCD, 1 RPL,DDSDDT,DRPLDE,DRPLDT,STRAN,DSTRAN, 2 TIME,DTIME,TEMP,DTEMP,PREDEF,DPRED,MA TERL,NDI,NSHR,NTENS, 3 NSTA TV,PROPS,NPROPS,COORDS,DROT,PNEWDT,CELENT, 4 DFGRD0,DFGRD1,NOEL,NPT,KSLAY,KSPT,KSTEP,KINC) C INCLUDE 'ABA_PARAM.INC' C CHARACTER*80 MA TERL DIMENSION STRESS(NTENS),STA TEV(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=STA TEV(1) S3O=STA TEV(2) SSS=STA TEV(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)
邓肯张模型模拟
研究生课程作业邓肯张模型参数计算 学生姓名李俊 学科专业岩土工程 学号201420105614 任课教师周小文教授 作业提交日期2014年12月
1．计算轴向应变 c h h ?∑= 1ε 式中 1ε－轴向应变； h ?∑－固结下沉量，由轴向位移计测得 0h －土样初始高度 c h —按实测固结下沉的试样高度 c h ?—试样固结下沉量 2．计算按实测固结下沉的试样高度，面积： 式中 Ac －按实测固结下沉的试样面积 0V －土样初始体积 3．计算剪切过程中试样的平均面积： 式中 a A －剪切过程中平均断面积 c V －按实测固结下沉的试样的体积 i V ?－排水剪中剪切时的试样体积变化 按体变管或排水管读数求得 1h ?－固结下沉量，由轴向位移计测得 3. 计算主应力差 c i c h V V A ?-= 01 h h V V A c i c a ?-?-= C c c A h V ?=
103 1?=-a A CR σσ 式中 31σσ- － 主应力差 1σ―大主应力 3σ－小主应力 C －测力计率定系数 R －测力计读数 2 数据处理 2.1 3σ=100kPa 数据初步计算 当3σ=100kPa 时，各数据初步计算如表1所示。 围压100kPa 数据初步计算表 表1
2.1.1 由切线模量计算数据 对公式 ) (311 σσε-=a +b 1ε进行直线拟合，如图1所示。 图1 1131 /()~εσσε-拟合曲线 a =0.0002，1 i E a = =5000kPa b ==0.0028，()131 ult b σσ-= =263.16kPa ()13f σσ-=204.26kPa ，()()1313f f ult R σσσσ-= -=0.7762 2.1.2 由泊松比计算数据 对公式()313/f D εεε-=+-进行直线拟合，如图2所示。
在ansys中导入自定义本构模型---邓肯-张模型
邓肯-张模型的关键点是材料的弹性模量随大小主应力差及小主应力(围压)的变化而变化，用APDL实现之的基本思路是：给每个单元定义一个材料号，分级施加荷载，在每个荷载步结束时提取出各单元的大小主应力，据此计算出下个荷载步的弹性模量Et,修改各单元之MP，用于下一步计算。 以下是一个简单算例，copy出去可直接运行。 !!!常规三轴试验模拟 !!!by taomingxing,NWPU !!!2003.7.16 FINISH /CLEAR /TITLE,Numerical Simulation of three axes testing of soils /PREP7 *dim,SUy,array,50 !Settlement records *dim,MaxPs,array,120 !Max history p1-p3 *dim,MaxDs,array,120 !Max history Ds !*dim,EEt,array,50 !Et of elememt !!!Duncan-Chang Model !!!Symbols：c-粘滞力，Fai-内摩擦角，Sf-破坏强度(p1-p3)f，Ds-应力水平，Pa-大气压,P3-围压 *CREATE,Duncan-Chang !Creat Macro file *afun,deg !Unit of angle *set,Pa,1e5 *set,P1,-ArrS3(i) !注意：岩土工程中应力为拉负压正
*set,P3,-ArrS1(i) *if,P3,LT,0.1*Pa,then P3=0.1*Pa !围压最小取值 *endif Sf=2*(c*cos(Fai)+P3*sin(Fai))/(1-sin(Fai)) !Mohr-Coulomb破坏强度(p1-p3)f Ds=(P1-P3)/Sf !应力水平， *if,Ds,GT,0.95,then Ds=0.95 !应力水平最大取值 *endif !判断加卸荷，如果(P1-P3)小于历史最大值视为卸荷－再加荷过程 *if,MaxPs(i),LT,P1-P3,then Ei=k*Pa*(P3/Pa)**n Et=Ei*(1-Rf*Ds)**2 !加荷情况的切线模量 MaxPs(i)=P1-P3 !保存历史最大应力 *elseif,MaxPs(i),GE,P1-P3 Et=Kur*Pa*(P3/Pa)**n !卸荷模量 *endif mp,ex,i,Et !修改单元i的Et mp,nuxy,i,Mu *END !!!单元类型 et,1,42 !平面四节点单元
ANSYS邓肯张材料模型
ANSYS邓肯-张材料模型 楼主给的在ANSYS上实现邓肯－张模型的方法很有用， 但其中还有几点需要修正的，这也是楼上的兄弟们有疑问的原因。我把楼主的代码运行了一下，然后对照作了修改，现在上传一下，有问题的兄弟可以仔细对照一下，在这里我对其中几个比较明显的问题说明一下： 1.MP命令不能直接给单元加材料，这是对的。在这里，楼主遗漏了一下命令：MPCHG，具体见下面的修改过的代码。 2.关于密度的问题。这些要在宏中定义，每修改一种材料(即调用一次邓肯－张子程序)就要修改一次材料的密度，其他有关材料的问题可以类推。 3.关于施加重力的问题。要在调用宏后，在同一个循环中重新定义一下重力。 以下是我修改过的楼主的代码，希望对兄弟们有所帮助。 !用APDL得到初步成果，贴于此供感兴趣的朋友参考，不当之处敬请指正，!欢迎加以完善。 !基本思路： !邓肯-张模型的关键点是材料的弹性模量随大小主应力差 !及小主应力(围压)的变化而变化，用APDL实现之的基本思路是： !给每个单元定义一个材料号，分级施加荷载，在每个荷载步结束时提取出各 !单元的大小主应力，据此计算出下个荷载步的弹性模量Et,修改各单元之MP，!用于下一步计算。 !以下是一个简单算例，copy出去可直接运行。 !!!常规三轴试验模拟
!********************************************************** FINISH /CLEAR /TITLE,Numerical Simulation of three axes testing of soils /PREP7 *dim,SUy,array,50!Settlement records *dim,MaxPs,array,120!Max history p1-p3 *dim,MaxDs,array,120!Max history Ds !*dim,EEt,array,50!Et of elememt !!!Duncan-Chang Model !!!Symbols：c-粘滞力，Fai-内摩擦角，Sf-破坏强度(p1-p3)f， !Ds-应力水平，Pa-大气压,P3-围压 !********************************************************************** *CREATE,Duncan-Chang!Creat Macro file *afun,deg!Unit of angle *set,Pa,1e5 *set,P1,-ArrS3(i)!注意：岩土工程中应力为拉负压正 *set,P3,-ArrS1(i) *if,P3,LT,0.1*Pa,then P3=0.1*Pa!围压最小取值 *endif Sf0=2*(c0*cos(Fai)+P3*sin(Fai))/(1-sin(Fai))!Mohr-Coulomb破坏强度
邓肯张模型
以土的常三轴实验学习Duncan-Chang本构关系模型 一、实验过程 1、试样制备 试验土样取自于南水北调焦作段一处工程，取回后，人工制成含水量15%的土体。在实验制样过程中，由于含水量较高，所以在通过制样器后，土柱未能成型，于是在原来土样基础上，添加了较干的土，再在制样器侧壁涂抹凡士林。最后制成高度7厘米，直径3.5厘米的土柱实验样品 2、不固结不排水（UU）剪切试验 试验是在土木工程学院深部矿井重点实验室进行的，试验装置如图1所示。 图1 常三轴实验仪
主要试验步骤为 (1)记录体变管的初始读数； (2)对试样加周围压力，并在周围压力下固结。当孔隙水压力的读数接近零时，说明固结完成，记下排水管的读数； (3)开动马达，合上离合器，按0.0065%/min的剪切应变速率对试样加载。按百分表读数为0，30，60，90，120，150，180，210，240，300，360，420，480，540，600，660，?，的间隙记读排水管读数和量力环量表读数，直到试样破坏为止。 二、邓肯张双曲线模型 到目前为止，国内外学者提出的土体本构模型不计其数，但是真正广泛用于工程实际的模型却为数不多，邓肯-张模型为其中之一。该模型是一种建立在增量广义虎克定律基础上的非线性弹性模型，可经反映应力～应变关系的非线性，模型参数只有８个，且物理意义明确，易于掌握，并可通过静三轴试验全部确定，便于在数值计算中运用，因而，得到了广泛地应用。 1、邓肯-张双曲线模型的本质 邓肯-张双曲线模型的本质在于假定土的应力应变之间的关系具有双曲线性质，见图2（a)。
图2(a ） 12()~a σσε- 双曲线 图2（b) 1131/()~εσσε-关系 图2 三轴试验的应力应变典型关系理论图