基于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°?；炷林脫Q體采用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