基于FLAC3D和ANSYS的岩质边坡混凝土置换加固效应的组合网格法

戴妙林,张 胤,鞠志敏,陈五一,祖 威

(1.河海大学水利水电学院,江苏南京 210098;2.中国水电顾问集团成都勘测设计研究院,四川成都 610072)

基于FLAC3D和ANSYS的岩质边坡混凝土置换加固效应的组合网格法

戴妙林1,张 胤1,鞠志敏1,陈五一2,祖 威2

(1.河海大学水利水电学院,江苏南京 210098;2.中国水电顾问集团成都勘测设计研究院,四川成都 610072)

通过系统分析组合网格模型的计算原理,为达到计算收敛快且精度高的要求,提出结合FLAC3D和ANSYS两种数值软件的组合网格法对岩质边坡混凝土置换加固效应进行模拟求解,同时采用等效插值方法解决组合网格间边界不匹配的问题。实例计算结果表明,此方法适用于岩质边坡混凝土置换加固效应的研究,可充分发挥FLAC3D与ANSYS软件各自的功能优势,能减少建模工作量,节省模型计算时间。

组合网格法;混凝土置换;岩质边坡;FLAC3D;ANSYS

在岩质边坡的加固中,混凝土置换是一种常用的加固措施。进行加固效应的数值计算之前,要建立计算模型的网格,按照通常的置换体模拟方法,会在很大程度上增加建模的工作量,特别对于地形和地质条件复杂的边坡,尤为如此。Ferket等[1]于1996年开始研究简单的有限差分法离散化椭圆型方程边值问题的复合网格。McCormick等[2-3]在有限差分离散的椭圆边值问题上提出采用离散化的方法统一整体粗网格和局部一个或多个均匀精细网格来达到逼近连续的解决方案。陈文平等[4-5]提出了整体粗网格求解和局部细网格反复迭代求解的方法,其结果与常规有限元方法所求得的解相符,为求解大型实际复杂问题提供了一个好的方法和思路。唐菊珍等[6-7]解决了网格局部奇性问题,提出了一种基于有限元自动生成系统(FEPG)的组合网格算法,该算法采用两套网格求解,在整个求解区域采用较粗网格,不考虑奇异的影响;而在奇异附近区域采用较细的网格,考虑奇异的影响;整体粗网格求解和局部细网格求解反复迭代,求得最终结果。2010年孙振波等[8]提出基于重叠型区域分解原理的快速自适应组合网格(fast adaptive composite grid,FAC)方案,解决了存在地下水时断裂和防渗帷幕难以模拟的问题,并且能够在提高求解效率的同时达到预期目标。

针对加固效应数值计算建模难协调的问题,有必要探索新途径和新方法,本文尝试将组合网格法应用于岩质边坡置换加固效应的研究中,采用天然状态和加固体两套独立模型网格,根据加固力等效的原理,在两套网格之间进行迭代计算,达到预期的精度时停止计算,岩质边坡采用FLAC3D计算,加固体采用ANSYS计算,这种方法既能体现置换体的加固效应,又可减少建模的工作量,充分发挥两种软件的功能和综合优势,是一种值得探讨和在工程实际中尝试应用的方法。

1 算法原理

组合网格法采用两套网格求解,以岩质边坡混凝土置换加固效应研究为例,即分别对边坡岩体和混凝土置换体建立独立的网格,对边坡岩体采用较粗网格,并且不考虑置换体的影响,而对混凝土置换体采用较细的网格,在两套网格之间进行迭代求解。为了阐述组合网格法的理论依据,首先讨论以下齐次边值问题:

式中:L为椭圆微分算子;f为中间变量;Ω为L对应的计算区域;Γ为Ωr和Ωc共用的模型边界(图1);v、u为方程求解变量;下标r、c分别代表粗网格中岩体和细网格中混凝土置换体。对Ω分别进行网格剖分,得到规则或不规则粗细网格Tc和Tr(Tc比Tr网格尺寸小很多),并有相应的有限元空间Sr和Sc。

图1 边坡岩体粗网格和置换体细网格

同时有对应的虚功方程为

2 算法实现

2.1 组合网格法迭代算法

在式(3)的基础上,可实现组合网格法的迭代计算,具体计算步骤如下:

以上为组合网格法迭代计算的计算步骤和数学表达形式,下面结合边坡加固研究,具体说明实施方法。

a.进行边坡岩体整体区域网格模型计算,其中不含加固体网格。

b.采用a的计算结果,插值计算得到加固体区域的边界条件。

c.根据b的加固体边界条件,采用加固体区域网格,以及加固体的材料参数,计算影响反力。

其中燃料煤平均低位发热量QL根据环境统计数据取值22051 kJ/kg；过量空气系数α根据《火电厂大气污染物排放标准》燃煤锅炉取值1.4；发电煤耗g采用《浙江省电力行业节能环保白皮书》中发电标准煤耗284 g/kWh。根据公式(4)、式(5)和式(6)计算得出300、750和1000 MW机组的废气排放量M分别为4.87×109、1.22×1010和1.62×1010 m3。

d.根据b的加固体边界条件,采用加固体区域网格,以及加固体所在部位的天然岩体的材料参数,计算影响反力。

e.将c得到的边界反力减去d的边界反力,得到加固体区域边界的反力差值,并通过插值作用到边坡岩体整体区域上,转至a,进入下一流程计算,若满足收敛精度,即停止计算。采用现流程反力差值与上一个流程反力差值的差值范数来表达收敛精度。

2.2 边界节点位移和力的计算

在粗细两套网格间进行数值传递时,为保证边界节点的位移及节点力的连续性,本文通过坐标变换,在每个单元上建立局部坐标系,实现粗网格边界单元节点与细网格对应节点的转换,将粗网格某一单元的边界节点信息设置为细网格对应单元系列的范围信息,在其范围内节点均根据各自所处位置进行加权平均插值处理;反之,在细网格边界单元上节点也能在粗网格单元内找到唯一的对应点。为建立这种坐标变换,最方便最直观的方法是将坐标变换式表示成为关于整体坐标的插值函数,以8节点六面体单元为例有

式中:uk(k=x,y,z)为细网格边界节点3个方向的位移;Ni(ξ,η)为单元形函数,其中ξ、η为局部坐标轴对应数值;vki为细网格边界节点落在粗网格单元内第i节点的位移。

FLAC3D是由Itasca公司研发推出的连续介质力学分析软件,能够模拟并计算三维岩土体及其他介质中工程结构的受力与变形形态。内置网格单元生成器提供了多种基本形状的网格单元,网格单元生成器可以将不同形状的网格单元连接起来,以得到期望的几何模型,并且还能适应不同的边界形状,适合较大较复杂的单元网格计算。由于FLAC3D采用了一个“显示解”方案,所以非线性的应力应变关系的求解所花费的时间几乎与线性本构关系相同,比隐式求解方案花费时间少;此外,由于采用迭代法求解不用存储刚度矩阵,因此采用中等容量的内存就可以较快求解大型结构,模拟大变形问题所用时间与模拟小变形问题差不多,在非线性计算方面,具有收敛快等优势。

采用有限元软件对混凝土力学特性模拟是非常普遍的,作为世界有限元软件的杰出代表,ANSYS软件中的混凝土本构模型的功能很强大,有许多种材料模型可供选择,如多线性随动强化模型,在合理选取参数后可以比较全面描述及模拟混凝土破坏曲线的下降段,能够反映混凝土的软化特性;其中用于混凝土结构的非线性Solid65单元可以很好地考虑和体现混凝土的压溃和开裂,在模拟混凝土力学特性方面具有很大优势。

2.4 基于FLAC3D和ANSYS的组合网格法计算流程

考虑到FLAC3D与ANSYS两种软件的特点,前者在单元类型和非线性计算方面具有较强的优势,后者在模拟混凝土力学特性方面的功能较强,综合应用两种软件进行组合网格法计算,粗网格的边坡岩体模型采用FLAC3D计算得到一组细网格对应节点位移,再通过ANSYS的Solid65单元计算细网格的混凝土置换体得到加固与不加固的边界节点反力差,中间采用Fortran语言编程,循环调用FLAC3D 和ANSYS进行粗细网格计算,相互迭代至某一精度后停止循环计算。具体流程见图2。

图2 组合网格模型计算流程

3 实例计算分析

3.1 边坡计算模型

选用文献[9]中边坡模型作为计算分析实例,该模型xOz平面示意图见图3。假设模型的求解为平面应变问题,即y向的厚度取2m,设置1层单元。边坡高度H=120 m,坡角60°,坡脚至右侧边长距离R=1.5H,坡顶宽度L=2.5H,坡顶至底部高度B= 2H。潜在滑动体的高度h=80 m,滑动面倾角为30°,滑动面设置厚度为0.2 m的薄层单元且下部约1/3处位置布置尺寸为4m×4m的置换洞。滑动体、滑动面薄层、下盘岩体和置换体的体积模量分别为14.67GPa、1.80GPa、24.40GPa和16.70 GPa,剪切模量分别为8.80 GPa、0.60GPa、16.80 GPa和12.50 GPa,密度分别为2 700 kg/m3、2 700 kg/m3、2 850 kg/m3和2400 kg/m3。下盘岩体与滑动体均设为弹性体,滑动面薄层采用Mohr-Coulomb屈服准则,其黏聚力为104 kPa,摩擦角为26.6°。混凝土置换体采用ANSYS中的Solid65单元模拟,黏聚力为2 000 kPa,摩擦角为46.6°,抗拉强度ft=1.43 MPa,抗压强度fc=14.30 MPa,等压双轴抗压强度fcb=1.2fc,静水压力下的双轴抗压强度f1=1.45fc,单轴抗拉强度f2=1.725fc。

图3 模型轮廓尺寸与分区示意图

3.2 网格剖分

常见的数值计算模型在采用通常的方法进行加固效应的数值计算时,边坡和置换体处于同一个体系,采用一套网格,混凝土置换体边界节点与相连岩体的周边节点位置是一一对应的。计算模型中的置换体剖分成图4(a)所示形式的网格。

在采用组合网格法进行加固计算时,需要建立两套独立的网格,第一套为由FLAC3D软件进行计算的边坡岩体粗网格,第二套为由ANSYS软件进行计算的置换体细网格。本实例计算采用两种计算方案:第一种计算方案,第一套网格与上述通常的数值计算模型网格相同,第二套网格为图4(a)的形式,两套网格在边界上的节点位置是一一对应的;第二种计算方案,第一套网格与第一种计算方案相同,但第二套网格则采用图4(b)的形式,其中置换体的网格全部为规则的长方体,与第一套网格中的被置换区域部分的网格形态(图4(a))是不相同的,单元数量增加近9倍,而且两套网格的边界节点位置并不完全相同,其中第二套网格的边界节点位移采用公式(5)计算。

图4 混凝土置换体网格剖分

3.3 加固效应分析

3.3.1 天然状态计算分析

在对加固效应计算前,首先分析天然状态下的安全度。如图5所示,使用FLAC3D和ANSYS两种软件分别计算,得到安全系数均为1.20,与刚体极限平衡法计算结果相同。折减系数为1.20时坡肩x向位移组合网格法计算为0.950 mm,小于ANSYS计算的3.626 mm。此现象可能是由于这两种软件的收敛方式与计算原理不同造成的,但从整体角度研究模型安全度的大小时,所得安全系数是相同的。

图5 天然状态下模型x向位移与折减系数关系曲线

3.3.2 加固分析

从组合网格模型的两种计算方案结果来看(图6),两者的结果非常接近,坡肩的x向位移突变点均在折减系数为1.45处,即安全系数均为1.45。第二种计算方案与第一种计算方案的加固效应计算结果基本相同,表明加固体的网格剖分可以独立于边坡岩体的网格剖分,对置换体采用较细的网格,有助于对其应力、变形等进行重点研究,体现了组合网格法的优越性。为了与通常的计算方法进行比较,采用第一种计算方案的网格形态,在同一套网格模型中,考虑置换加固效应,物理力学参数设置同前,采用ANSYS计算,坡肩的x向位移随折减系数的关系曲线也示于图6中,安全系数亦为1.45,与组合网格法分析结果相同。折减系数为1.45时,组合网格法第一种计算方案计算的坡肩x向位移为0.209mm,第二种计算方案为0.217 mm,小于ANSYS计算的2.875mm,这是由于在组合网格法计算中,边坡岩体的粗网格模型采用FLAC3D计算,所以,与天然状态计算位移的规律相同。

图6 加固状态下模型x向位移与折减系数关系曲线

4 结 语

利用Fortran语言编程实现了FLAC3D与ANSYS两种软件在岩质边坡混凝土加固中的组合使用,既能充分发挥软件各自的优势,提高模型计算效率,又能回避岩体边坡网格与置换体的协调问题,使得基于FLAC3D与ANSYS组合网格法在计算收敛性与混凝土力学特性方面均有较好的体现。同时,该组合网格法的加固体与岩体可以采用两套不相同的独立网格,既能突出重点研究的置换加固体,又能减少建模的工作量,故可以进一步尝试将组合网格法应用于复杂岩质边坡多处置换加固效应的研究。

[1]FERKET P J J,REUSKEN A.Further analysis of the local defect correction method[J].Computing,1996,56(2): 117-139.

[2]McCORMICKSF,THOMASJ.Thefastadaptivecomposite grid method(FAC)for elliptic boundary value problems[J].Math Comp,1986,46:439-456.

[3]WHITCOMB J D.Iterative global/local finite element analysis[J].Computers and Structures,1991,40(4): 1027-1031.

[4]陈文平,马昌凤,蒋利华.组合网格法在搅拌摩擦焊接数值模拟中的应用[J].安徽理工大学学报:自然科学版,2009,29(2):57-61.(CHENWenping,MA Changfeng,JIANG Lihua.Application of composite grid method in numerical simulation of friction stir welding [J].Journal of Huainan Institute of Technology:Natural Science,2009,29(2):57-61.(in Chinese))

[5]郭胜山,李德玉,唐菊珍.水工结构动接触问题的组合网格算法[J].水力发电,2009,35(5):24-59.(GUO Shengshan,LI Deyu,TANG Juzhen.A composite grid methodofdynamicalcontactprobleminhydraulic structure[J].Water Power,2009,35(5):24-59.(in Chinese))

[6]唐菊珍.基于FEPG的组合网格法与接触问题的Lagrange乘子法[D].桂林:桂林电子科技大学,2006.

[7]唐菊珍,马昌凤.基于有限元自动生成系统(FEPG)的组合网格算法[J].广西科学院学报,2007,23(1):7-10.(TANG Juzhen,MA Changfeng.A composite grid method based on the system of finite element program generator(FEPG)[J].Journal of Guangxi Academy of Sciences,2007,23(1):7-10.(in Chinese))

[8]孙振波,朱国荣,江思珉,等.地下水模型的快速自适应组合网格求解方法[J].水文地质工程地质,2010,37 (1):18-21.(SUN Zhenbo,ZHU Guorong,JIANG Simin, et al.Solving methods of groundwater model with fast adaptive composite grids[J].Hydrogeology＆Engineering Geology,2010,37(1):18-21.(in Chinese))

[9]戴妙林,黄天成,张胤,等.单滑动面岩质边坡数值计算影响因素分析[J].水电能源科学,2011,29(2):79-81. (DAI Miaolin,HUANG Tiancheng,ZHANG Yin,et al. Influencing factors analysis of numerical calculation for rocky slopewithsingleslidingsurface[J].Water Resources and Power,2011,29(2):79-81.(in Chinese))

[10]周桂云.基于强度折减的边坡稳定安全系数有限元迭代解法[J].水利水电科技进展,2010,30(3):58-61. (ZHOU Guiyun.Finite element iteration method for safety factor of slope stability based on strength reduction method [J].Advances in Science and Technology of Water Resources,2010,30(3):58-61.(in Chinese))

Composite grid method of FLAC3D and ANSYS for reinforcement effect of concrete replacement for rock slopes//

DAI Miaolin1,ZHANG Yin1,JU Zhiming1,CHEN Wuyi2,ZU Wei2(1.College of Water Conservancy and Hydropower Engineering,Hohai University,Nanjing 210098,China;2.HydroChina Chengdu Engineering Corporation,Chengdu 610072,China)

The principles of the composite grid method were systematically analyzed.In order to achieve fast convergence and high precision of calculations,a composite grid method of FLAC3D and ANSYS was proposed to simulate the reinforcement effect of concrete replacement for rock slopes.Simultaneously,the equivalent interpolation method was employed to solve the non-uniform problem of the boundary between thick and fine grids.The calculated results of a case show that the proposed method is applicable to research on the reinforcement effect of concrete replacement for rock slopes. It can give full play of the advantages of both FLAC3D and ANSYS,reduce the workload of modeling,and save computing time of the model.

composite grid method;concrete replacement;rock slope;FLAC3D;ANSYS

10.3880/j.issn.10067647.2013.03.016

TU753.8

A

10067647(2013)03007205

2012-06-23 编辑:熊水斌)

国家自然科学基金(51079044);水利部公益性行业科研专项(201001035);江苏省普通高校研究生科研创新计划(CXLX12_0252)

戴妙林(1964—),男,浙江温岭人,副教授,硕士,主要从事水工结构和岩土工程数值计算及试验研究。E-mail:mldai@hhu.edu.cn