戴妙林,張 胤,鞠志敏,陳五一,祖 威

(1.河海大學水利水電學院,江蘇南京 210098;2.中國水電顧問集團成都勘測設計研究院,四川成都 610072)

基于FLAC3D和ANSYS的巖質邊坡混凝土置換加固效應的組合網(wǎng)格法

戴妙林1,張 胤1,鞠志敏1,陳五一2,祖 威2

(1.河海大學水利水電學院,江蘇南京 210098;2.中國水電顧問集團成都勘測設計研究院,四川成都 610072)

通過系統(tǒng)分析組合網(wǎng)格模型的計算原理,為達到計算收斂快且精度高的要求,提出結合FLAC3D和ANSYS兩種數(shù)值軟件的組合網(wǎng)格法對巖質邊坡混凝土置換加固效應進行模擬求解,同時采用等效插值方法解決組合網(wǎng)格間邊界不匹配的問題。實例計算結果表明,此方法適用于巖質邊坡混凝土置換加固效應的研究,可充分發(fā)揮FLAC3D與ANSYS軟件各自的功能優(yōu)勢,能減少建模工作量,節(jié)省模型計算時間。

組合網(wǎng)格法;混凝土置換;巖質邊坡;FLAC3D;ANSYS

在巖質邊坡的加固中,混凝土置換是一種常用的加固措施。進行加固效應的數(shù)值計算之前,要建立計算模型的網(wǎng)格,按照通常的置換體模擬方法,會在很大程度上增加建模的工作量,特別對于地形和地質條件復雜的邊坡,尤為如此。Ferket等[1]于1996年開始研究簡單的有限差分法離散化橢圓型方程邊值問題的復合網(wǎng)格。McCormick等[2-3]在有限差分離散的橢圓邊值問題上提出采用離散化的方法統(tǒng)一整體粗網(wǎng)格和局部一個或多個均勻精細網(wǎng)格來達到逼近連續(xù)的解決方案。陳文平等[4-5]提出了整體粗網(wǎng)格求解和局部細網(wǎng)格反復迭代求解的方法,其結果與常規(guī)有限元方法所求得的解相符,為求解大型實際復雜問題提供了一個好的方法和思路。唐菊珍等[6-7]解決了網(wǎng)格局部奇性問題,提出了一種基于有限元自動生成系統(tǒng)(FEPG)的組合網(wǎng)格算法,該算法采用兩套網(wǎng)格求解,在整個求解區(qū)域采用較粗網(wǎng)格,不考慮奇異的影響;而在奇異附近區(qū)域采用較細的網(wǎng)格,考慮奇異的影響;整體粗網(wǎng)格求解和局部細網(wǎng)格求解反復迭代,求得最終結果。2010年孫振波等[8]提出基于重疊型區(qū)域分解原理的快速自適應組合網(wǎng)格(fast adaptive composite grid,FAC)方案,解決了存在地下水時斷裂和防滲帷幕難以模擬的問題,并且能夠在提高求解效率的同時達到預期目標。

針對加固效應數(shù)值計算建模難協(xié)調的問題,有必要探索新途徑和新方法,本文嘗試將組合網(wǎng)格法應用于巖質邊坡置換加固效應的研究中,采用天然狀態(tài)和加固體兩套獨立模型網(wǎng)格,根據(jù)加固力等效的原理,在兩套網(wǎng)格之間進行迭代計算,達到預期的精度時停止計算,巖質邊坡采用FLAC3D計算,加固體采用ANSYS計算,這種方法既能體現(xiàn)置換體的加固效應,又可減少建模的工作量,充分發(fā)揮兩種軟件的功能和綜合優(yōu)勢,是一種值得探討和在工程實際中嘗試應用的方法。

1 算法原理

組合網(wǎng)格法采用兩套網(wǎng)格求解,以巖質邊坡混凝土置換加固效應研究為例,即分別對邊坡巖體和混凝土置換體建立獨立的網(wǎng)格,對邊坡巖體采用較粗網(wǎng)格,并且不考慮置換體的影響,而對混凝土置換體采用較細的網(wǎng)格,在兩套網(wǎng)格之間進行迭代求解。為了闡述組合網(wǎng)格法的理論依據(jù),首先討論以下齊次邊值問題:

式中:L為橢圓微分算子;f為中間變量;Ω為L對應的計算區(qū)域;Γ為Ωr和Ωc共用的模型邊界(圖1);v、u為方程求解變量;下標r、c分別代表粗網(wǎng)格中巖體和細網(wǎng)格中混凝土置換體。對Ω分別進行網(wǎng)格剖分,得到規(guī)則或不規(guī)則粗細網(wǎng)格Tc和Tr(Tc比Tr網(wǎng)格尺寸小很多),并有相應的有限元空間Sr和Sc。

圖1 邊坡巖體粗網(wǎng)格和置換體細網(wǎng)格

同時有對應的虛功方程為

2 算法實現(xiàn)

2.1 組合網(wǎng)格法迭代算法

在式(3)的基礎上,可實現(xiàn)組合網(wǎng)格法的迭代計算,具體計算步驟如下:

以上為組合網(wǎng)格法迭代計算的計算步驟和數(shù)學表達形式,下面結合邊坡加固研究,具體說明實施方法。

a.進行邊坡巖體整體區(qū)域網(wǎng)格模型計算,其中不含加固體網(wǎng)格。

b.采用a的計算結果,插值計算得到加固體區(qū)域的邊界條件。

c.根據(jù)b的加固體邊界條件,采用加固體區(qū)域網(wǎng)格,以及加固體的材料參數(shù),計算影響反力。

其中燃料煤平均低位發(fā)熱量QL根據(jù)環(huán)境統(tǒng)計數(shù)據(jù)取值22051 kJ/kg；過量空氣系數(shù)α根據(jù)《火電廠大氣污染物排放標準》燃煤鍋爐取值1.4；發(fā)電煤耗g采用《浙江省電力行業(yè)節(jié)能環(huán)保白皮書》中發(fā)電標準煤耗284 g/kWh。根據(jù)公式(4)、式(5)和式(6)計算得出300、750和1000 MW機組的廢氣排放量M分別為4.87×109、1.22×1010和1.62×1010 m3。

d.根據(jù)b的加固體邊界條件,采用加固體區(qū)域網(wǎng)格,以及加固體所在部位的天然巖體的材料參數(shù),計算影響反力。

e.將c得到的邊界反力減去d的邊界反力,得到加固體區(qū)域邊界的反力差值,并通過插值作用到邊坡巖體整體區(qū)域上,轉至a,進入下一流程計算,若滿足收斂精度,即停止計算。采用現(xiàn)流程反力差值與上一個流程反力差值的差值范數(shù)來表達收斂精度。

2.2 邊界節(jié)點位移和力的計算

在粗細兩套網(wǎng)格間進行數(shù)值傳遞時,為保證邊界節(jié)點的位移及節(jié)點力的連續(xù)性,本文通過坐標變換,在每個單元上建立局部坐標系,實現(xiàn)粗網(wǎng)格邊界單元節(jié)點與細網(wǎng)格對應節(jié)點的轉換,將粗網(wǎng)格某一單元的邊界節(jié)點信息設置為細網(wǎng)格對應單元系列的范圍信息,在其范圍內節(jié)點均根據(jù)各自所處位置進行加權平均插值處理;反之,在細網(wǎng)格邊界單元上節(jié)點也能在粗網(wǎng)格單元內找到唯一的對應點。為建立這種坐標變換,最方便最直觀的方法是將坐標變換式表示成為關于整體坐標的插值函數(shù),以8節(jié)點六面體單元為例有

式中:uk(k=x,y,z)為細網(wǎng)格邊界節(jié)點3個方向的位移;Ni(ξ,η)為單元形函數(shù),其中ξ、η為局部坐標軸對應數(shù)值;vki為細網(wǎng)格邊界節(jié)點落在粗網(wǎng)格單元內第i節(jié)點的位移。

FLAC3D是由Itasca公司研發(fā)推出的連續(xù)介質力學分析軟件,能夠模擬并計算三維巖土體及其他介質中工程結構的受力與變形形態(tài)。內置網(wǎng)格單元生成器提供了多種基本形狀的網(wǎng)格單元,網(wǎng)格單元生成器可以將不同形狀的網(wǎng)格單元連接起來,以得到期望的幾何模型,并且還能適應不同的邊界形狀,適合較大較復雜的單元網(wǎng)格計算。由于FLAC3D采用了一個“顯示解”方案,所以非線性的應力應變關系的求解所花費的時間幾乎與線性本構關系相同,比隱式求解方案花費時間少;此外,由于采用迭代法求解不用存儲剛度矩陣,因此采用中等容量的內存就可以較快求解大型結構,模擬大變形問題所用時間與模擬小變形問題差不多,在非線性計算方面,具有收斂快等優(yōu)勢。

采用有限元軟件對混凝土力學特性模擬是非常普遍的,作為世界有限元軟件的杰出代表,ANSYS軟件中的混凝土本構模型的功能很強大,有許多種材料模型可供選擇,如多線性隨動強化模型,在合理選取參數(shù)后可以比較全面描述及模擬混凝土破壞曲線的下降段,能夠反映混凝土的軟化特性;其中用于混凝土結構的非線性Solid65單元可以很好地考慮和體現(xiàn)混凝土的壓潰和開裂,在模擬混凝土力學特性方面具有很大優(yōu)勢。

2.4 基于FLAC3D和ANSYS的組合網(wǎng)格法計算流程

考慮到FLAC3D與ANSYS兩種軟件的特點,前者在單元類型和非線性計算方面具有較強的優(yōu)勢,后者在模擬混凝土力學特性方面的功能較強,綜合應用兩種軟件進行組合網(wǎng)格法計算,粗網(wǎng)格的邊坡巖體模型采用FLAC3D計算得到一組細網(wǎng)格對應節(jié)點位移,再通過ANSYS的Solid65單元計算細網(wǎng)格的混凝土置換體得到加固與不加固的邊界節(jié)點反力差,中間采用Fortran語言編程,循環(huán)調用FLAC3D 和ANSYS進行粗細網(wǎng)格計算,相互迭代至某一精度后停止循環(huán)計算。具體流程見圖2。

圖2 組合網(wǎng)格模型計算流程

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的置換洞?；瑒芋w、滑動面薄層、下盤巖體和置換體的體積模量分別為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 模型輪廓尺寸與分區(qū)示意圖

3.2 網(wǎng)格剖分

常見的數(shù)值計算模型在采用通常的方法進行加固效應的數(shù)值計算時,邊坡和置換體處于同一個體系,采用一套網(wǎng)格,混凝土置換體邊界節(jié)點與相連巖體的周邊節(jié)點位置是一一對應的。計算模型中的置換體剖分成圖4(a)所示形式的網(wǎng)格。

在采用組合網(wǎng)格法進行加固計算時,需要建立兩套獨立的網(wǎng)格,第一套為由FLAC3D軟件進行計算的邊坡巖體粗網(wǎng)格,第二套為由ANSYS軟件進行計算的置換體細網(wǎng)格。本實例計算采用兩種計算方案:第一種計算方案,第一套網(wǎng)格與上述通常的數(shù)值計算模型網(wǎng)格相同,第二套網(wǎng)格為圖4(a)的形式,兩套網(wǎng)格在邊界上的節(jié)點位置是一一對應的;第二種計算方案,第一套網(wǎng)格與第一種計算方案相同,但第二套網(wǎng)格則采用圖4(b)的形式,其中置換體的網(wǎng)格全部為規(guī)則的長方體,與第一套網(wǎng)格中的被置換區(qū)域部分的網(wǎng)格形態(tài)(圖4(a))是不相同的,單元數(shù)量增加近9倍,而且兩套網(wǎng)格的邊界節(jié)點位置并不完全相同,其中第二套網(wǎng)格的邊界節(jié)點位移采用公式(5)計算。

圖4 混凝土置換體網(wǎng)格剖分

3.3 加固效應分析

3.3.1 天然狀態(tài)計算分析

在對加固效應計算前,首先分析天然狀態(tài)下的安全度。如圖5所示,使用FLAC3D和ANSYS兩種軟件分別計算,得到安全系數(shù)均為1.20,與剛體極限平衡法計算結果相同。折減系數(shù)為1.20時坡肩x向位移組合網(wǎng)格法計算為0.950 mm,小于ANSYS計算的3.626 mm。此現(xiàn)象可能是由于這兩種軟件的收斂方式與計算原理不同造成的,但從整體角度研究模型安全度的大小時,所得安全系數(shù)是相同的。

圖5 天然狀態(tài)下模型x向位移與折減系數(shù)關系曲線

3.3.2 加固分析

從組合網(wǎng)格模型的兩種計算方案結果來看(圖6),兩者的結果非常接近,坡肩的x向位移突變點均在折減系數(shù)為1.45處,即安全系數(shù)均為1.45。第二種計算方案與第一種計算方案的加固效應計算結果基本相同,表明加固體的網(wǎng)格剖分可以獨立于邊坡巖體的網(wǎng)格剖分,對置換體采用較細的網(wǎng)格,有助于對其應力、變形等進行重點研究,體現(xiàn)了組合網(wǎng)格法的優(yōu)越性。為了與通常的計算方法進行比較,采用第一種計算方案的網(wǎng)格形態(tài),在同一套網(wǎng)格模型中,考慮置換加固效應,物理力學參數(shù)設置同前,采用ANSYS計算,坡肩的x向位移隨折減系數(shù)的關系曲線也示于圖6中,安全系數(shù)亦為1.45,與組合網(wǎng)格法分析結果相同。折減系數(shù)為1.45時,組合網(wǎng)格法第一種計算方案計算的坡肩x向位移為0.209mm,第二種計算方案為0.217 mm,小于ANSYS計算的2.875mm,這是由于在組合網(wǎng)格法計算中,邊坡巖體的粗網(wǎng)格模型采用FLAC3D計算,所以,與天然狀態(tài)計算位移的規(guī)律相同。

圖6 加固狀態(tài)下模型x向位移與折減系數(shù)關系曲線

4 結 語

利用Fortran語言編程實現(xiàn)了FLAC3D與ANSYS兩種軟件在巖質邊坡混凝土加固中的組合使用,既能充分發(fā)揮軟件各自的優(yōu)勢,提高模型計算效率,又能回避巖體邊坡網(wǎng)格與置換體的協(xié)調問題,使得基于FLAC3D與ANSYS組合網(wǎng)格法在計算收斂性與混凝土力學特性方面均有較好的體現(xiàn)。同時,該組合網(wǎng)格法的加固體與巖體可以采用兩套不相同的獨立網(wǎng)格,既能突出重點研究的置換加固體,又能減少建模的工作量,故可以進一步嘗試將組合網(wǎng)格法應用于復雜巖質邊坡多處置換加固效應的研究。

[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]陳文平,馬昌鳳,蔣利華.組合網(wǎng)格法在攪拌摩擦焊接數(shù)值模擬中的應用[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]郭勝山,李德玉,唐菊珍.水工結構動接觸問題的組合網(wǎng)格算法[J].水力發(fā)電,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的組合網(wǎng)格法與接觸問題的Lagrange乘子法[D].桂林:桂林電子科技大學,2006.

[7]唐菊珍,馬昌鳳.基于有限元自動生成系統(tǒng)(FEPG)的組合網(wǎng)格算法[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]孫振波,朱國榮,江思珉,等.地下水模型的快速自適應組合網(wǎng)格求解方法[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]戴妙林,黃天成,張胤,等.單滑動面巖質邊坡數(shù)值計算影響因素分析[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]周桂云.基于強度折減的邊坡穩(wěn)定安全系數(shù)有限元迭代解法[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);水利部公益性行業(yè)科研專項(201001035);江蘇省普通高校研究生科研創(chuàng)新計劃(CXLX12_0252)

戴妙林(1964—),男,浙江溫嶺人,副教授,碩士,主要從事水工結構和巖土工程數(shù)值計算及試驗研究。E-mail:mldai@hhu.edu.cn