FLAC_3D_粘弹性模型的二次开发及其应用_徐平

收稿日期:2003O 11O 08基金项目:国家自然科学基金委员会与清江水电开发有限公司联合资助重点项目(59939190);水利部创新项目(SCX2002-20)作者简介:徐 平(1965O ),男,湖南华容人,长江科学院教授级高级工程师,博士,主要从事岩石力学试验与数值分析研究.文章编号:1001O 5485(2004)02O 0010O 04FLAC 3D粘弹性模型的二次开发及其应用徐 平1,李云鹏2,丁秀丽1,王芝银2(1.长江科学院岩基研究所,湖北武汉 430010; 2.石油大学,北京 102249)摘要:针对FL AC 3D (V ersion 2.0)提供的粘弹性本构模型难以描述岩体变形随时间发展趋于稳定的现象,对FL AC 3D 的粘弹性模拟功能进行了开发,导出广义Kelvin 模型的中心差分格式,并采用FL AC 3D 内嵌功能较强的FISH 语言和FLA C 3D 的相关命令研制了广义Kelvin 模型的接口程序。

通过检验算例说明了开发思路和接口程序的正确性,同时还指出了二次开发时应注意的一些问题,并对基坑开挖进行了实例分析,结果符合基坑开挖变形规律。

关 键 词:粘弹性;二次开发;广义K elvin 模型;FLA C 3D 软件中图分类号:T P311.52;T U 452 文献标识码:A1 概 述FLAC 3D 软件是一个用于分析和解决工程力学问题的三维显式差分法程序,可以有效地模拟岩土体及其它材料结构的工程力学行为,与有限元法等应用软件相比各自具有不同的优缺点。

在FLAC 3D软件中,具有能满足岩土工程结构分析的基本单元模式,包含了进行弹性分析的各向同性、异性和横观各向同性的弹性本构模型,塑性模型有Mohr O Coulomb 弹塑性、Drucker O Prager 弹塑性模型等;描述材料蠕变特性的有Maxw ell 、幂指数等5种蠕变本构模型。

对于静态、动态、蠕变、渗流、温度效应及其相互耦合的工程问题,通过选择适当的模式可以实现对不同类型问题的求解。

目前FLAC 3D (Ver -sion 2.0)软件在我国尤其是在岩土工程界得到了广泛的认可和应用[1~5],并取得了具有工程指导意义的应用成果。

尽管如此,FLAC 3D 软件在应用过程中仍然暴露出了不少的缺点和不足。

如FLAC 3D的前处理部分仍采用不太直观的文本命令格式,通过数据文件与软件连接,虽然FLAC 3D也提供了内嵌式的FISH 语言,但在处理复杂网格时仍给使用者带来了很大的不便。

文[6]就其前处理部分的不足,借助于For -tran 语言对此项进行了功能开发与研制;文[7]为了更细致地了解岩石变形局部化现象,研制了局部化图形再现后的处理程序。

由于岩体力学特性的复杂性,其变形往往呈现出明显的时间性,FLAC 3D程序提供了5种适用于不同特性岩体介质的粘弹塑性分析模型,如Maxw ell 模型描述的是变形随时间变化而不断增大的情况,即属于非稳定流变模型,其对岩体变形随时间的发展趋于稳定状态的情况就很难描述。

FLAC 3D 中其它流变模型均有其特定的适用范围,这就明显影响了FLAC 3D 的广泛应用。

故此,本文以FLAC 3D(Version 2.0)为基础,利用其内嵌的FISH 语言工具,开发研制了广义Kelvin 模型接口程序,从而扩大了FLAC 3D 程序的适用性。

通过算例分析表明,所研制的接口程序开发思路正确,可用于分析工程实际问题。

2 广义Kelvin 模型及其差分格式工程实践表明,许多岩体及支护结构在工程挠动后产生一定变形,但随着时间的推移变形逐渐趋于稳定。

为使FLAC 3D软件能够有效地模拟这一变形特性,本文给出了广义Kelvin 模型的本构方程,并导出图1 广义Kelvin 模型Fig.1 Kelvin O Voigt model适应于FLAC 3D 软件的差分格式。

2.1 广义Kelvin 模型广义Kelvin 模型如图1所示,其微分型的本构关系可表示为[8]R +p 1¤R =q 0E +q 1¤E ,(1)式中:p 1=G 1E +E 1;q 0=EE 1E +E 1;q 1=E G 1E +E 1;G 为材料粘性系数;E,E 1分别为弹性模量和粘弹性模量。

引入算子,式(1)可改写为P c R =Q c E ,(2)第21卷第2期长 江 科 学 院 院 报Vol.21No.22004年4月Journal of Yangtze River Scientific Research Institute A pr.2004式中:P c =(p 1d d t +1);Q c =(q 1dd t+q 0)。

将式(2)扩充到三维状态,得到广义Kelvin 模型的三维本构关系为P d S ij =2Q d e ij ,(3)式中:P d =(p 1c d d t +1);Q d =(q 1c dd t +q 0c );q 1c =G 0G 1G 0+G 1;q 0c =G 0G 1G 0+G 1;p 1c =G 1G 0+G 1;(4)S ij 为偏应力张量;e ij 为偏应变张量;G 0,G 1分别为弹性剪切模量和粘弹性剪切模量。

2.2 差分格式将式(4)中P d ,Q d 代入式(3),改为增量形式并化简得p 1c$S ij $t +S ij =2(q c $e ij$t+q 0c e ij ),(5)式中$符号表示相应物理量的增量。

为了有效地与FLAC 3D 软件连接,采用中心差分格式将上式改写为S c ij -S 0ij $t +S c ij +S 0ij2p c 1=2q 1c p 1c $e c ij $t +2q 0cp 1c(e c ij +e 0i j ),(6)式中:S 0ij为t 时刻差分点的偏应力张量;S c ij 为t +$t 时刻差分点的偏应力张量;e 0ij 为t 时刻差分点的偏应变张量;e c ij 为t +$t 时刻差分点的偏应变张量;$e c ij 为在$t 内差分点的偏应变张量增量。

将式(6)整理并解得t +$t 时单元的偏应力张量为S c ij =(C 1S 0ij +C 2$e c ij +C 3e 0ij )#C 4,(7)式中:C 1=1-G 0+G 12G 1$t ;C 2=2G 0+G 0G 1G 1$t ;C 3=2G 0G 1G 1$t ;C 4=(1+G 0+G 12G 1$t )-1。

则t +$t 时单元的应力张量为R c ij =S c ij +R c m D ij ,(8)式中:R c ij 为单元应力张量;R c m 为单元平均应力,即R c m =13R 0kk +K $E c kk ,其中R 0kk =R 0x +R 0y +R 0z ,$E c kk =$E c x +$E c y +$E c z ,K 为弹性体积模量;D ij 为Kronecker 符号。

3 二次开发与关键函数3.1 二次开发流程FLAC 3D的主要计算思路是首先由初应力或外荷载计算获得节点不平衡力,并以等效在节点上的力和质量建立以节点位移速率为未知量的节点运动方程,解得节点位移速率后并进一步求得某一时步的单元应变增量、单元应力增量,以及经累加得到当前时步的单元总应力。

然后再以单元应力等求得节点不平衡力。

如此循环上述过程,最终得到工程结构的应力与变形结果。

根据上述分析思路,并利用FLAC 3D 内嵌的FISH 语言开发的接口程序主要流程如图2所示。

图2 接口程序主流程图F ig.2 M ain flow chart of interface code3.2 关键函数在FLAC 3D 软件系统中内嵌有功能丰富的程序语言FISH,用户可根据不同的要求设定新的变量或函数,实现对单元类型、特殊边界、力学参数、节点位移、速率、单元应力,及其特性的随时修改和伺服控制等。

在广义Kelvin 本构模型的开发中用到了FISH 的节点变量函数、单元变量函数及FLAC 3D 的控制语句和命令语句等。

关键函数、命令语句详见文[9]。

4 算例验证及应用通过上述思路及FLAC 3D 命令与FISH 语言,研制并完成了广义Kelv in 接口程序。

以下给出检验算例和一工程实例,以说明所开发程序的有效性。

4.1 验证算例一粘弹性矩形厚板,尺寸为20m @5m @2m,弹性剪切模量G 0=5@108Pa,体积模量K =1.5@109Pa,粘弹性剪切模量G 1=3.2@109Pa,粘性参数G 1=1.574@1010Pa #d 。

由对称性取1/4模型进行研究,在x =0.0面上x ,y 方向约束,在z =0.0,z =2.0面为z 单向约束。

在x =10.0面上作用有沿x 方向的均匀边界载荷(R xx =20MPa),如图3所示。

用开发的广义Kelv in 模型进行计算,其结果与理论解的比较列于图4和图5中。

模拟流变时间为8.64d,图4中显示的应力曲11第2期徐 平等 FL AC 3D 粘弹性模型的二次开发及其应用图3 计算模型及载荷示意图Fig.3 Illustration o f numerical model andevenly load along x O axis direction图4 单元应力随时间变化规律及等色区图Fig.4 Cur ves of element str ess vs time andstress iso chromatic zonepattern图5 x 向位移随时间变化规律及等色区图Fig.5 Curv es of displacement in the x O direction vstime and disp.isochromatic zone pattern线为中心坐标点(7.5,2.5,0.5)处单元z 向正应力值,其中,依次由上向下第一条曲线为解析曲线,第二条曲线为开发的FLAC 解,其相对误差为8.5%。

图5显示的是坐标点为(7.0,3.0,1.0)处的x 方向的位移值,上方曲线为解析解,下方曲线为FLAC解,相对误差为9.5%。

由应力等值线和位移等值线可见其规律性基本上与理论解分布规律接近。

由图5可知,广义Kelvin 模型所获结果达到稳定状态的时间约为2d 左右。

图6 第1,3,5步开挖后点1和点2的位移随时间的变化曲线Fig.6 Curves of displacement in the No.1and No.2p oints vs time after excavation of step 1,3and 54.2 工程应用某深基坑,由对称性取1/4区域(10m @10m @10m)进行分析研究,每步开挖深度为0.8m,共分5步,开挖区域为1m @4m @2m 。