郝 哲,朱一飛,王鐵男

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

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

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

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

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

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

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

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

1 計算程序

1.1 拉格朗日差分法

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

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

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

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

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

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

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

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

2.1 計算模型

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

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

圖1 壩體計算模型

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

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)格結點表示出來.

圖2 網(wǎng)格劃分

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

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

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

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

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

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

故該層物理力學參數(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)很小,甚至可以忽略不計,故其邊界可以看做為被固定的剛性邊界(即沒有變形,也沒有位移),對于壩頂邊界來說是沒有約束的.

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

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

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

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

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

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

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

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

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

2.5 計算結果分析

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

2.5.1 應力分析

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

圖3 主應力矢量(ST)

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

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

2.5.2 位移分析

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

圖6 位移矢量(D)

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

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

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

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

3 結 論

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

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

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

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

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

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

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

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

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

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

[10] 方朝陽,夏富洲.土石壩安全監(jiān)測綜合評判準則研究[J].武漢大學學報,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-),男,遼寧沈陽人,沈陽大學教授,博士后研究人員.