郝 哲,朱一飛,王鐵男

(1.沈陽大學(xué) 建筑工程學(xué)院,遼寧 沈陽 110044; 2.東北大學(xué) 資源與土木工程學(xué)院,遼寧 沈陽 110004)

基于差分法的土石壩數(shù)值模擬研究

郝 哲1,朱一飛2,王鐵男1

(1.沈陽大學(xué) 建筑工程學(xué)院,遼寧 沈陽 110044; 2.東北大學(xué) 資源與土木工程學(xué)院,遼寧 沈陽 110004)

對有限差分法的基本原理及分析過程進行了闡述,并對FLAC數(shù)值計算軟件的功能、處理過程、應(yīng)用范圍以及與其他軟件的比較等方面進行了系統(tǒng)論述.基于現(xiàn)場調(diào)研收集的相關(guān)壩體的資料,應(yīng)用FLAC程序?qū)Ω沸码姀S四灰場主壩進行位移、變形、應(yīng)力分析,以判斷壩體結(jié)構(gòu)的穩(wěn)定狀態(tài).結(jié)果表明,該壩體在應(yīng)力、位移、變形上是滿足要求的,是穩(wěn)定的,在長時間內(nèi)不會出現(xiàn)問題.

阜新電廠;土石壩;差分法;FLAC;數(shù)值模擬

用當?shù)赝亮虾蜕?、砂礫、卵礫、石渣、石料等筑成壩,依靠壩體自身重量維持壩的穩(wěn)定,是一種古老而至今還廣泛使用的大壩形式,稱為土石壩[1].《阜新發(fā)電廠四灰場主壩壩體穩(wěn)定性監(jiān)測研究》是沈陽大學(xué)承擔的橫向課題.課題組對電廠四灰場主壩進行現(xiàn)場調(diào)研,了解了大壩歷經(jīng)20多年后的運行現(xiàn)狀,對其外型、邊坡養(yǎng)護、水庫水位、周圍環(huán)境、邊坡植被等情況進行了考察,并收集了壩體設(shè)計的相關(guān)數(shù)據(jù)資料、場地地質(zhì)資料及物理力學(xué)參數(shù)指標等,作為本文分析的基礎(chǔ)資料.

阜新電廠四灰場主壩位于阜新市北11 km處,是該電廠粉煤灰及廢水的主要排泄場地[2].場區(qū)地面標高界于219.40～417.60 m之間,為低山丘陵區(qū).整個大壩呈橢圓形,壩的南側(cè)為一“U”形谷,谷底較平坦,寬約100 m.溝谷西側(cè)為黃土陡崖,崖高6～8 m,近于直立.東側(cè)為居民區(qū),壩肩部分為巖石山坡,部分為土山坡,山坡坡度約為15°,山坡上植被較差.西側(cè)壩肩為巖石山坡,山坡坡度約為12°.大壩北側(cè)為一灰渣庫.在主壩東北角處有一副壩,呈橢圓形,壩長約為50 m,壩寬約20 m.主壩壩體為均質(zhì)土石壩,壩高35 m,壩長448.14 m,頂寬為5 m,底寬為45 m.壩基標高為224 m,壩體外側(cè)有一層碎石護坡,厚度約為0.3 m;內(nèi)側(cè)部分為碎石護坡,部分為草皮護坡.壩內(nèi)水位很低,水主要集中在北側(cè),南側(cè)部分可見壩底.

1 計算程序

1.1 拉格朗日差分法

所謂差分法,是把基本方程和邊界條件(一般均為微分方程)近似地改用差分方程(代數(shù)方程)來表示,把求解微分方程的問題改換成求解代數(shù)方程.

本文采用的計算方法為拉格朗日差分法,這是一種利用拖帶坐標系分析大變形問題的數(shù)值方法.模型經(jīng)過網(wǎng)絡(luò)劃分,物理網(wǎng)格映射成數(shù)學(xué)網(wǎng)格,數(shù)學(xué)網(wǎng)格上的某個結(jié)點就與物理網(wǎng)格上的相應(yīng)結(jié)點坐標相對應(yīng).對某一個結(jié)點而言,在每一時刻它受到來自其周圍區(qū)域的合力影響.如果合力不等于零,結(jié)點就具有了失穩(wěn)力,就要產(chǎn)生運動,根據(jù)牛頓定律,結(jié)點就要產(chǎn)生加速度,進而可以在一個時步中求得速度和位移的增量.對每一個區(qū)域而言,可以根據(jù)其周圍結(jié)點的運動速度求得它的應(yīng)變率,然后根據(jù)材料的本構(gòu)關(guān)系求得應(yīng)力的增量.由應(yīng)力增量求出 t和 t+Δt時刻各個結(jié)點的不平衡力和各個結(jié)點在 t+Δt時的變加速度.對加速度進行積分,即可得結(jié)點的新坐標值.由于物體的變形,單元發(fā)生局部的平均旋轉(zhuǎn)或整旋,只要計算相應(yīng)的應(yīng)力改正值,通過應(yīng)力疊加就可以得到新的應(yīng)力值,到此計算為一個循環(huán).然后按時步進行下一輪的計算,如此一直進行到問題收斂[3,4].

1.2 數(shù)值模擬軟件選取

關(guān)于數(shù)值計算方法,目前較為常用的有:有限元法、邊界元法、有限差分法、加權(quán)余量法、離散元法、剛體元法、不連續(xù)變形分析法、流形元法等.前四種方法是基于連續(xù)介質(zhì)力學(xué)的方法,隨后三種方法是基于非連續(xù)介質(zhì)力學(xué)的方法,最后一種方法具有這兩大類方法的共性.有限單元法是最常用的,但一般用有限元計算的位移通常較實測的結(jié)果要小,有時甚至差一兩個數(shù)量級,究其原因,主要是在計算中忽略了非線性的大變形和沿弱面的不連續(xù)變形.有限差分法是將問題的基本方程和邊界條件以簡單、直觀的差分形式來表達,專門求解巖土力學(xué)非線性大變形問題,能得到較滿意的位移解,易于在實際工程中應(yīng)用.

目前,面向各學(xué)科專業(yè),針對不同的數(shù)值計算方法,已開發(fā)出眾多的數(shù)值計算軟件,其中有通用程序,也有專用程序.表1列出了目前國內(nèi)外較知名的大、中型軟件.

表1 各種數(shù)值計算軟件功能對比

其中,FLAC程序(Fast Language Analysis of Continue)是基于拉格朗日差分法的國際知名大型數(shù)值計算軟件,由美國 ITASCA Consulting Group.Inc.開發(fā)研制.主要應(yīng)用于巖土力學(xué)分析,能對其應(yīng)力、位移進行精確計算,它采用顯式拉格朗日算法及混合離散單元劃分技術(shù),代替了原先廣泛使用的隱式有限元法[5].能夠精確地模擬材料的塑性流動和破壞,對靜態(tài)系統(tǒng)模型也可采用動態(tài)方程來進行求解.因為不需要形成剛度矩陣,故占用微機內(nèi)存小,便于求解大型工程問題.與有限元計算相比,FLAC解線性問題較慢,而解非線性問題較快,特別是在大變形和幾何非線性的情況下,具有明顯優(yōu)勢.FLAC程序具有計算靈活的特點,是其他有限元程序無法比擬的,特別是它通過類似批處理文件的前處理可以任意改變原始參數(shù),同時還能保持前后計算的連續(xù)性.也就是說后面參數(shù)改變,可以繼續(xù)在原基礎(chǔ)上計算,而且這種變化十分方便.這樣就為變彈模法模擬時間因素的影響創(chuàng)造了條件,而其他大多數(shù)有限元程序很難做到這一點.FLAC程序可以用于分析多種巖石材料、斷層和節(jié)理,可以模擬不同加載條件下的地應(yīng)力場生成、邊坡或地下硐室開挖、大壩和地基的誘發(fā)破壞、混凝土襯砌、錨桿或錨索設(shè)置、地震工程和巖石爆破、地下滲流等.本文的數(shù)值計算采用FLAC 3.0程序.

2 阜新電廠四灰場主體壩體FLAC分析

2.1 計算模型

計算壩體結(jié)構(gòu)時,為了使地基對結(jié)構(gòu)應(yīng)力的影響反映出來,必須把和結(jié)構(gòu)相連的一部分地基取為彈性體,與結(jié)構(gòu)一起作為計算對象.按照彈性力學(xué)中關(guān)于接觸應(yīng)力的理論,所取地基范圍大小應(yīng)視結(jié)構(gòu)底部的寬度而定(與結(jié)構(gòu)高度無關(guān)).早期文獻[6,7]中,一般都建議在結(jié)構(gòu)的兩邊和下方把地基范圍取為大致等于結(jié)構(gòu)底部寬度B.但后來一些文獻[8]中,大都把所取的范圍擴大為L= 2B,個別文獻[9]還把它擴大為 L=4B.此外,還有一些文獻[10]認為,應(yīng)當把地基范圍取為矩形區(qū)域,以便將鉸支座改為連桿支座,減少對地基的人為約束.筆者認為在地基比較均勻而且結(jié)構(gòu)與地基彈性相差不大的情況下,沒有必要使 L超過2B,用鉸支座比連桿支座更接近實際情況,地基范圍的形狀影響也不大.

對于阜新電廠主壩,由于為土壩,故其結(jié)構(gòu)與地基結(jié)構(gòu)相當,取 L=2B即可.為便于計算,取地基為矩形區(qū)域比較符合程序要求.得計算模型如圖1所示.

圖1 壩體計算模型

FLAC軟件在給定的區(qū)域或范圍內(nèi)可賦予的本構(gòu)模型有:橫觀各向同性彈性模型(Anisotropic)、彈性本構(gòu)模型(Elastic)、摩爾-庫侖彈塑性模型(Mohr-Coulomb)、零模型(Null)、應(yīng)變軟化模型(SS)、徹體節(jié)理模型(Ubiquitous).由于土體有一定的流塑性,巖石也有一定的彈塑性,故選擇摩爾-庫侖彈塑性本構(gòu)模型.

2.2 網(wǎng)格劃分

根據(jù)程序要求,用相隔等間距 h且平行于坐標軸的兩組平行線組成網(wǎng)格.劃分橫向網(wǎng)格數(shù)為225 m/5 m=45,縱向網(wǎng)格數(shù)為(90 m+35 m)/5 m =25,網(wǎng)格總個數(shù)為45×25=1 125.網(wǎng)格劃分如圖2所示,且壩體橫截面的邊緣輪廓均能由網(wǎng)格結(jié)點表示出來.

圖2 網(wǎng)格劃分

2.3 物理力學(xué)參數(shù)分析

根據(jù)阜新電廠四灰場地質(zhì)勘測資料,按三層考慮:壩基土層厚度約為5.4 m;砂石層厚度為10.1 m;底層為花崗巖.各層物理力學(xué)參數(shù)如下:

(1)對于壩體本身.由勘測資料可知,筑壩材料為黃土.該土為黃～黃褐色,含少量碎石、角礫,稍濕,硬塑;最大干容重為1.80～1.81 g/cm3,最佳含水質(zhì)量分數(shù)為15.4%～15.8%,可溶鹽質(zhì)量分數(shù)為0.03%,密度為1.4～1.61 g/cm3,內(nèi)摩擦角為22.3°,彈性模量為104kPa.據(jù)此計算剪切模量S=3.79 MPa,體積模量 K=9.26 MPa,粘聚力COH=0.049 MPa.

故該層物理力學(xué)參數(shù)為:S=3.79 MPa,K= 9.26 MPa,DENS=1.47×10-3kg/cm3,FRIC= 22.3°,COH=4.9×10-2MPa.

(2)花崗巖層.根據(jù)勘測資料可知,主壩地段花崗巖實驗結(jié)果如下:吸水率 0.73%,容重2.56 g/cm3,比重 2.62 g/cm3,風干抗壓強度69.20 MPa,軟化系數(shù)0.91,內(nèi)摩擦角37°,泊松比0.384.故彈性模量 E=1.774×107kPa.剪切模量 S=6.4×103MPa,體積模量 K=2.54×104MPa,粘聚力COH=21 MPa.

故該層物理力學(xué)參數(shù)為:S=6.4×103MPa, K =2.54 ×104MPa,DENS = 2.56 × 10-3kg/cm3,FRIC=37°,COH=21 MPa.

(3)地基土石層.由于該層是由土層和砂石層組成,且每層厚度較小,在程序中不能一一表示出來,所以由土層和砂石層的各項力學(xué)參數(shù)取其平均值可得出結(jié)果.土層參數(shù)可與筑壩黃土相仿,砂石層力學(xué)參數(shù)可根據(jù)相關(guān)資料查得.密度取2.0 g/cm3,摩擦角取30°.剪切模量 K=5.7 ×102MPa,體積模量 K=5.6×102MPa,粘聚力取COH=10 MPa.

故該層物理力學(xué)參數(shù)為:S=5.7×102MPa, K=5.6×102MPa,DENS=2.0×10-3kg/cm3, FRIC=30°,COH=10 MPa.

2.4 邊界條件及外部載荷的加入

由計算模型可知,由于所選計算區(qū)域的邊界所受擾動已經(jīng)很小,甚至可以忽略不計,故其邊界可以看做為被固定的剛性邊界(即沒有變形,也沒有位移),對于壩頂邊界來說是沒有約束的.

由《施工圖階段工程、水文地質(zhì)勘察任務(wù)書》可知,發(fā)展遠景壩頂標高為276 m,發(fā)展遠景最高水位標高為271 m,故最高水位離壩頂5 m.水深h為30 m,則水對壩體底部側(cè)壓應(yīng)力為

式中,σ水為水的密度,kg/m3;g為重力加速度, m/s-2.

由勘測資料可知,地下水位為2 m,即在壩體有水側(cè)面的水深為30 m,壩基15 m為土砂層和水層,15 m以下為花崗巖層.對于壩體無水側(cè)面,壩基15 m為土砂層,地下水高為13 m,15 m以下為花崗巖層,該壩體左側(cè)為有水側(cè)面.則對于網(wǎng)格(0,0)點處的應(yīng)力為

同理,網(wǎng)格(0,15)點處的應(yīng)力為

網(wǎng)格(0,18)點處的應(yīng)力為:

σx=σy=1.0×103×10×30=0.3 MPa.式中,h1為土壩中水頭高度(m),h2為地基土石層中水頭高度(m),h3為花崗巖層中水頭高度(m),μ為泊松比,γ為重度(N/m3).

對于壩體無水側(cè)面,由于地下水位為2 m,在程序中較小,無法表示出來,可認為從地基土層到頂部應(yīng)力為均勻變化的.則對于網(wǎng)格處(45,0)點處的應(yīng)力為

對于網(wǎng)格(45,15)點處的應(yīng)力為

網(wǎng)格(45,18)點處的應(yīng)力為

2.5 計算結(jié)果分析

將壩體及壩基各層材料的物理力學(xué)參數(shù)、邊界條件和外部荷載按照FLAC程序所需格式組成數(shù)據(jù)文件,輸入程序運算.

2.5.1 應(yīng)力分析

對于所取計算模型的主應(yīng)力矢量、水平方向應(yīng)力的等值線、豎直方向應(yīng)力的等值線,由程序計算可得結(jié)果如圖3～圖5所示.從圖中可知,主應(yīng)力變化是比較均勻的.對于整個模型來說,最大主應(yīng)力在模型底層,主要是由各層材料的自重力產(chǎn)生.而對于壩體應(yīng)力從上往下變化幾乎是線性的,其最大主應(yīng)力也在壩體底部.從圖4、圖5可知,壩體對其周圍的應(yīng)力影響在豎直方向較大,水平方向較小,水平方向的影響主要集中在土砂層.計算結(jié)果的最大主應(yīng)力為2.814 MPa,在允許范圍內(nèi).

圖3 主應(yīng)力矢量(ST)

圖4 水平方向應(yīng)力等值線(SXX)

圖5 豎直方向應(yīng)力等值線(SYY)

2.5.2 位移分析

對于所取計算模型的位移矢量、水平方向位移等值線、豎直方向位移等值線,由程序計算可得結(jié)果如圖6～圖8.

圖6 位移矢量(D)

圖7 水平方向位移等值線(XD)

圖8 豎直方向位移等值線(YD)

從圖中可知,對于整個模型總位移主要集中在壩體,壩基的位移由上往下越來越小.對于水平方向的位移也是主要集中在壩體,而壩基幾乎沒有.豎直方向上的位移在壩體上尤為突出:一是由于水壓力的豎向分力,二是由于壩體的自重.壩基有一部分沉降,也是由于自重和水壓力.

當程序計算到800步時位移值已經(jīng)收斂,該點的位移值為3.093 cm,在允許范圍內(nèi)(設(shè)計水平位移在總寬度的0.5%以內(nèi),45 m×0.5%= 22.5 cm).壩體計算的最大水平位移為3.25 cm,最大豎直位移為2.25 cm,在允許范圍內(nèi)(設(shè)計最大沉降為總高度的1%以內(nèi),35 m×1%=35 cm).而對于整個壩體的最大位移為3.359 cm,也在允許范圍內(nèi).

3 結(jié) 論

本文基于差分法對阜新電廠四灰場主壩進行了數(shù)值模擬研究.通過對其進行實地整理分析,提煉出所需物理力學(xué)參數(shù)及荷載作用參數(shù).根據(jù)壩體實際情況,在FLAC程序計算中選擇摩爾—庫侖彈塑性本構(gòu)模型,并且選擇擴展基礎(chǔ)為矩形,有利于網(wǎng)格的劃分.通過應(yīng)用FLAC程序所作的數(shù)值分析可以得出,該壩體在應(yīng)力、位移方面是滿足要求的,是穩(wěn)定的.由于邊坡有一層護石,還有草皮護坡都可以增加壩體的穩(wěn)定性,在計算中比較保守,而且現(xiàn)在壩內(nèi)水位很低,個別地方壩體甚至接觸不到水,按照最高設(shè)計水位計算,壩體是穩(wěn)定的,長時間內(nèi)不會出現(xiàn)問題.對于壩體的滲流方面,由于需要長期監(jiān)測(3～5年,甚至10多年),本文沒有涉及到,但它是壩體監(jiān)測的一個重要組成部分.該壩由于水位低、水少,所以滲流方面也是符合規(guī)定的.

[1] 白永年,吳士寧,王洪恩.土石壩加固[M].北京:水利電力出版社,1992:2-3.

[2] 北京有色冶金設(shè)計研究院.阜新電廠灰渣庫設(shè)施施工圖、階段工程水文地質(zhì)勘測任務(wù)書[R].1983:20-26.

[3] 張東日,陶連金,李風儀,等.拉格朗日元法及其應(yīng)用軟件FLAC[J].礦山壓力與頂板管理,1997,14(3):224-226.

[4] 楊立強,張中杰,林舸,等.FLAC基本原理及其在地學(xué)中的應(yīng)用[J].地學(xué)前緣,2003,10(1):24.

[5] 龔紀文,崔建軍,席先武,等.FLAC數(shù)值模擬軟件及其在地學(xué)中的應(yīng)用[J].大地構(gòu)造與成礦學(xué),2002,26(3):321 -325.

[6] 閻中華.土石壩穩(wěn)定計算方法評述[J].人民黃河,1995, 17(7):33-37.

[7] 吳中如,顧沖時,沈振中,等.大壩安全綜合分析和評價的理論、方法及其應(yīng)用[J].水利水電科技進展,1998,18 (3):2-6,65.

[8] 傅瓊?cè)A,黃真.土石壩安全監(jiān)測及其資料整理分析方法綜述[J].江西水利科技,1997,23(2):124-128.

[9] 方衛(wèi)華.綜論土石壩的安全監(jiān)測[J].紅水河,2002,21 (4):64-67.

[10] 方朝陽,夏富洲.土石壩安全監(jiān)測綜合評判準則研究[J].武漢大學(xué)學(xué)報,2001,34(5):27-31.

Numerical Simulation Study on Earth-Rock Dam Based on Calculus of Difference

HAO Zhe1,ZHU Yifei2,WAN G Tienan1
(1.School of Architectural and Civil Engineering,Shenyang University,Shenyang 110044,China;2.College of Resources and Civil Engineering,Northeastern University,Shenyang 110004,China)

The basic theory and analytic procedure about the calculus of difference are elaborated.The function,disposed process,used realm ofFLAC and comparison to other software are introduced systematically.Based on information collected from Fuxin Electric Power Plant,FLAC was used to calculate the displacement,deformation and stress of the dam to determine the stability of the dam.The results show that the dam can meet the requirements at the aspects of stress,displacement,and it is stable in a long time without problems.

fuxin electric power plant;earth-rock dam;calculus of differences;FLAC;numerical simulation

TV 698.1

A

【責任編輯 楊 敏】

1008-9225(2010)02-0023-05

2009-03-22

郝 哲(1972-),男,遼寧沈陽人,沈陽大學(xué)教授,博士后研究人員.