超高堆石壩分層填筑的顆粒流模擬

黃小波

（永州市水利水電勘測設計院　永州市　425000）

超高堆石壩分層填筑的顆粒流模擬

黃小波

（永州市水利水電勘測設計院永州市425000）

針對當前堆石壩發(fā)展趨勢及壩體變形預測方法的不足，研究了細觀顆粒流方法用于模擬超高堆石壩應力變形特征的可行性。首先，利用FISH語言二次開發(fā)，實現(xiàn)在復雜邊界區(qū)域快速建立密實的堆石料模型，并充分體現(xiàn)堆石料的級配及破碎效應；然后，以水布埡面板工程為例，基于該建模方法對超高堆石壩的填筑過程進行仿真模擬。結果表明，顆粒流模型的壩體變形與壩體實測變形比較接近。由此可見，顆粒流方法應用于模擬超高堆石壩分層分期填筑等大尺度復雜邊界問題是切實可行的，從而為超高堆石壩變形預測提供新的途徑。

超高堆石壩顆粒流分層填筑

引言

20世紀以來，國際上相繼開工建設了一批200 m級的高面板堆石壩，取得了寶貴的經驗，但從已建成工程的初步運行狀況來看，這些高壩也出現(xiàn)了一些問題，反映了高堆石壩的新特征。目前，我國面板堆石壩的發(fā)展面臨著從200 m級高壩向300 m級超高壩的突破，對于300 m級超高堆石壩，其運行特性不可能簡單地通過對已建工程的觀測數(shù)據(jù)進行適當外推而得到，而必須對現(xiàn)有的設計準則和計算方法進行調整，以適應超高壩的實際情況。

變形控制是面板堆石壩設計階段最重要的考慮因素，對于300 m級的超高堆石壩，壩體堆石的變形控制將是決定其在技術上是否可行的關鍵因素。目前，基于有限元方法的數(shù)值模擬是預測堆石壩變形量的主要手段，然而，有限元方法雖然能夠定性地得到壩體和面板的應力變形規(guī)律，但是計算結果尚不能達到定量分析的精度，這主要是由于堆石體是一種典型的不連續(xù)介質，在荷載作用下，堆石料會出現(xiàn)破碎、滑移及充填，表現(xiàn)出非常復雜的細觀力學特性，而連續(xù)介質力學方法只能從宏觀層面計算壩體的應力變形，難以反映堆石料的上述細觀力學特性，計算結果往往與壩料的真實變形存在較大差距。

近些年來，離散元方法的提出和發(fā)展，為從細觀層面研究堆石壩的應力變形特征提供了新的思路與方法。但是,從目前的發(fā)展現(xiàn)狀來看，在堆石壩的模擬方面，離散元方法還僅僅只是用于研究壩料特性、碾壓及流變的模擬，研究對象主要是針對小尺度的物理試樣或低壩的簡單模擬，而對于高堆石壩的分層填筑等大尺度復雜工程問題，尚未進行深入研究。

對離散元方法而言，計算效率是制約其工程應用的重要因素之一，隨機多邊形模型雖能較為真實地反映堆石料的幾何特性，但是接觸識別十分復雜，嚴重影響計算效率，而顆粒流方法將顆粒簡化為圓形，接觸識別簡單，計算效率相對較高，這是本文采用顆粒流程序的重要原因。另外，計算過程中采用多核并行算法，模型的計算效率將進一步提高。因此,本文擬采用顆粒流程序PFC2D對超高堆石壩分層分區(qū)填筑模擬的可行性進行研究，從細觀層面對超高堆石壩的宏觀變形特性進行初步探索。

1　超高堆石壩分層填筑顆粒流模擬的基本原理

本文利用二維顆粒流方法對超高堆石壩的填筑進行模擬，重點在于實現(xiàn)壩體的分區(qū)分層填筑過程，目前暫不考慮壩體碾壓及流變效應。為了實現(xiàn)這一目的，首先需建立堆石料的顆粒流模擬方法，并充分考慮堆石料的復雜特性，然后，基于該方法建立填筑層模型，并最終實現(xiàn)壩體全斷面的填筑。

1.1堆石體模型的顆粒流建模方法

實驗研究表明，堆石體具有明顯的級配特性，幾何形狀為隨機不規(guī)則的多面體，在一定荷載作用下會出現(xiàn)破裂。本文利用FISH語言二次開發(fā)，提出了一種建立堆石料顆粒流模型的方法，在模型區(qū)域內快速生成顆粒，并同時綜合考慮堆石料的級配、隨機形狀及破碎效應，主要建模步驟如下所示：

（1）根據(jù)顆粒級配在模型區(qū)域內生成初始模型，并設置一個合適的粒徑d，將級配中粒徑超過d的粒料分為粗粒料，小于d的粒料分為細粒料，記錄初始模型中所有粗粒料的位置及半徑，并刪除全部顆粒。

（2）按照級配中細粒料及其相對含量在模型區(qū)域內再次生成顆粒模型。

（3）根據(jù)（1）中記錄的粗粒料形心坐標和粒徑，在（2）中生成的計算模型里確定對應的區(qū)域，將落入該區(qū)域的顆粒粘結在一起，構成顆粒簇，并賦予模型細觀力學參數(shù)。

圖1是按照根據(jù)上述方法建立的堆石料的顆粒流模型，其中黃色顆粒為細粒料，不可破碎，藍色顆粒為粗粒料，在一定應力條件下可發(fā)生破碎。

圖1　堆石料顆粒流模型

模型中粗粒料由若干彼此相鄰的細顆粒粘結在一起構成，圖2為模型中隨機選取的幾種不同形狀的粗粒料。

圖2中，紅圈表示初始模型中記錄的粗粒料位置，黃線表示顆粒簇內部的平行粘結，相鄰顆粒簇間的接觸為線性無粘結接觸。

圖2　粗粒料的隨機構成

上述顆粒流模型中，粗粒料的等效粒徑按照如下式進行計算：

式中n——構成某粗粒料的細顆粒數(shù)量；

ri——細顆粒的粒徑。

可以看出，該式按照面積等效的原則可以較為方便地計算得到顆粒簇的等效粒徑，但是沒有考慮顆粒簇內部的空隙面積，計算得到的顆粒簇的等效粒徑比實際粒徑偏小，由此得到的級配曲線中細粒料的含量會偏大。若忽略這一誤差，在截斷物理試樣級配細小粒徑的條件下，按照上述建模方法，數(shù)值模型的級配和物理試樣的級配是十分接近的，如圖3所示。

圖3　堆石體顆粒流模型的粒料級配

可以看出，通過上述建模方法，可基本保證數(shù)值模型的級配與壩料真實級配一致，并且粗粒料可發(fā)生破碎，充分體現(xiàn)了壩料的破碎效應。雖然模型中粗粒料的形狀與壩料真實形狀存在一定差異，但是該建模方法不需要定義粗粒料的形狀，也避免了粗粒料初始位置的搜索與判別，極大地簡化了建模過程，這對于實現(xiàn)超高堆石壩分層分區(qū)填筑過程的模擬，具有十分重要的意義。

1.2二維顆粒流模型的細觀特性

與三維模型相比，二維模型相對更加密實，且在垂直模型平面方向不受力或位移約束，這就使得二維模型的宏觀力學特性與三維模型相比存在較大差異。為了使二維模型能夠表現(xiàn)出與三維模型基本一致的力學特性，需要對二維顆粒流模型的細觀特性進行修改。

（1）孔隙率。由上可知，在相同條件下，二維模型模型相對更加密實，在荷載作用下易表現(xiàn)出過量的剪脹現(xiàn)象，難以較為真實的反映三維模型的體積變形特性，因此在建立二維顆粒流模型時，需要提高模型的孔隙率，以增強模型的可壓縮性。但是，若在初始建模時將模型的孔隙率設置得過大，會導致顆粒間接觸不良，或無法填滿整個模型區(qū)域，因此，本文所采取的方法是：首先，在生成顆粒時孔隙率取合適值，使顆粒充滿模型區(qū)域，又不至于使得顆粒間重疊量過大；然后，賦予顆粒接觸細觀力學參數(shù)，使得模型構成一個穩(wěn)定的受力體系；最后，隨機刪除模型中的部分細小顆粒，提高模型的孔隙率。

（2）顆粒轉動特性。顆粒流方法將顆粒簡化為理想的圓形，降低了顆粒間的咬合作用，對于二維顆粒流模型而言，由于缺少中間主應力方向的約束，顆粒的咬合作用進一步降低，因此在加載條件下，顆粒出現(xiàn)過量的轉動，使得模型宏觀上表現(xiàn)出明顯的應變軟化現(xiàn)象，這與三維堆石料模型的應變硬化特性不符。目前，在基于二維顆粒流模型模擬堆石料力學特性的研究成果中，為了解決這一問題，一種方法是將增強基本顆粒幾何形狀的復雜性，另一種方法是在接觸關系中賦予顆?？罐D特性[16]。顯然，后者建模過程更加方便，更是適合于模擬超高堆石壩分層填筑等復雜邊界問題。

1.3填筑層的模擬方法

分層填筑中一個難點是分層填筑過程的處理問題?；陬w粒流方法模擬壩體分層填筑時，首先需以填筑層輪廓為邊界生成墻體（walls），構建封閉的區(qū)域，然后在該封閉區(qū)域內生成顆粒，這就意味著，在填筑層中生成顆粒時，填筑層的邊界必須是已知的。但是實際上，在填筑層區(qū)域生成顆粒后，經過碾壓和重力作用，該填筑層會發(fā)生一定程度的沉降和側向變形，變形后填筑層的輪廓與其初始輪廓示意圖如圖4所示。由圖4可知，填筑層1碾壓沉降至穩(wěn)定狀態(tài)后，其頂部實際高程與設計高程間存在高度差。在模擬后續(xù)填筑層2時，在碾壓及自重作用下，填筑層2中的顆粒會首先將該高度差的區(qū)域填滿，同時加上自身的壓縮變形，填筑層2碾壓沉降至穩(wěn)定狀態(tài)后，其頂部實際高程與設計高程間會形成更加明顯的高度差，進而對后續(xù)的填筑層產生更加顯著的影響。由此，在利用顆粒流程序模擬壩體分層填筑時，必須使得各填筑層碾壓沉降至穩(wěn)定狀態(tài)后，其輪廓恰好與設計輪廓重合。

圖4　填筑層變形示意圖

為了滿足這一要求，本文提出了一種解決方法，即在模擬填筑層1時，在設計厚度的基礎之上增加一定厚度，以彌補該填筑層由于碾壓及自重而引起的壓縮變形，同時對上下游坡面施加法向約束；待填筑層沉降至穩(wěn)定狀態(tài)后，將該填筑層頂部設計高程以上顆粒全部刪除，并解除上下坡面的法向約束。

由于加高層的厚度遠小于填筑層厚度，當刪除設計輪廓線以上顆粒及解除坡面法向約束時，填筑層內會發(fā)生一定程度的卸荷，但卸荷量很小，可以忽略填筑層邊界由于卸荷而發(fā)生的改變。1.4蓄水作用的模擬方式

圖5　堆石壩填筑層的顆粒流模擬方法示意圖

超高面板堆石壩在蓄水時，在面板運行狀態(tài)良好的情況下，上游水頭基本以面力形式作用在面板上，并傳遞給壩體。因此本文在考慮蓄水作用時，將水壓力作為面力進行考慮。

假定上游水面高程為H，上游坡腳為α，對于面板上游面任意顆粒，設其半徑為R，形心坐標（x,y），厚度為t，則作用在該顆粒上的等效水壓力為：

2　超高堆石壩分層分區(qū)填筑的顆粒流模擬

基于上述基本方法與原理，本文以水布埡面板堆石壩為例，對顆粒流方法在超高堆石壩分層分期填筑模擬上的可行性進行研究。首先，通過對主、次堆石料的壓縮試驗進行仿真模擬，通過反演獲取壩料的細觀力學參數(shù)；然后，結合壩體施工的實際過程，利用顆粒流程序對壩體填筑過程進行仿真模擬，并將壩體的變形結果與實測資料及有限元方法進行對比，以驗證顆粒流方法在模擬超高堆石壩分層分區(qū)填筑問題上的有效性。

2.1壩料細觀力學參數(shù)的反演

顆粒流模型的宏細觀參數(shù)間不存在明確的對應關系，通常只能通過反演來確定與壩料宏觀特性相適應的細觀力學參數(shù)。本文中，通過對水布埡主、次堆石料的壓縮試驗進行仿真模擬，并通過反演的方式確定了壩料的細觀力學參數(shù)，如表1所示。

表1　水布埡堆石料細觀力學參數(shù)

表中各變量的物理含義：kn，線性接觸法向剛度，（N/m）；ks，線性接觸切向剛度，（N/m）；fric，細觀摩擦系數(shù)；pb_kn，平行連接法向剛度，（Pa/m）；pb_ks，平行連接切向剛度,（Pa/m）；pb_ten，細觀抗拉強度，（Pa）；pb_coh，細觀粘聚力，（Pa）；pb_fa，細觀摩擦角，（°）；dp_nratio，法向臨界阻尼系數(shù)；dp_sratio，切向臨界阻尼系數(shù)。

按照上述細觀力學參數(shù)，主、次堆石料的宏觀應力變形關系如圖6所示。

圖6　堆石料顆粒流模型宏觀應力應變關系

同時根據(jù)設計資料，水布埡堆石料的干密度為2 160 kg/m3，由于顆粒流模型的初始孔隙率為0.20，因此顆粒的密度取2 700 kg/m3。

2.2計算剖面及加載過程

水布埡面板堆石壩壩體最大剖面如圖7所示，根據(jù)壩體的施工組織設計，填筑過程的載荷分級如表2所示。

3　顆粒流模擬結果及分析

按照圖7和表2所示的加載過程，運用顆粒流方法對水布埡面板堆石壩的分層填筑過程進行模擬，可得到不同工況下壩體的宏觀應力變形結果及細觀變化過程。壩體填筑完成后，模型中顆粒數(shù)量為378 066。

圖7　水布埡面板堆石壩壩體分區(qū)及加載過程

3.1壩體變形

根據(jù)模擬結果，在不同的填筑工期，壩體的沉降位移分布如圖8、圖9所示。

表2　顆粒流模型仿真計算荷載及加載過程

從圖可以看出，壩體竣工后壩體沉降主要發(fā)生在次堆石區(qū)，主堆石區(qū)變形量相對較小?？⒐て趬误w最大沉降量約2.11 m，蓄水期壩體最大沉降量約2.24 m，均發(fā)生在壩軸線偏下游約305 m高程處。

3.2壩體受力

堆石體表現(xiàn)出復雜的力學特性，其主要原因就是內部力鏈網(wǎng)絡的復雜演化特性，圖10、圖11反映了不同工況下壩體內部堆石料間接觸力分布。

圖8　一期壩體填筑及蓄水工況下壩體沉降

圖9　三期壩體填筑及蓄水工況下壩體沉降

圖10　四期壩體填筑及蓄水工況下壩體沉降

由圖可以看出，壩體內部力鏈網(wǎng)絡中，強力鏈主要分布在壩軸線上游主堆石區(qū)內，而次堆石區(qū)內部的接觸力則相對較弱，同時，壩體蓄水后，上游壩坡附近堆石料間的接觸力強度明顯上升。

為便于與有限元結果進行比較，可將壩體內部顆粒間的接觸力轉換為顆粒集合體的等效平均應力，該轉換過程可通過測量圓實現(xiàn)。通過在模型中設置大量測量圓，測量圓的相對尺寸為20～33.2倍顆粒平均粒徑，于是可得到測量圓內顆粒集合體的平均應力分量，并合成主應力。通過計算，壩體內部的主應力分布如圖12、圖13所示。

圖中縱坐標為高程（m），應力單位MPa。從圖中可以看出，竣工期壩體最大主應力3.8 MPa，最小主應力1.0 MPa；蓄水期壩體最大主應力4.0 MPa，最小主應力1.1 MPa。蓄水后，壩體上游坡面主應力明顯升高。3.3模擬結果比較及分析

“九五”公關期間，多家單位運用不同模型對水布埡面板堆石壩的宏觀應力變形進行了分析計算，得出的結果如表3所示。

圖11　壩體內部堆石料間接觸力分布(N)

圖12　竣工期壩體內部主應力分布

圖13　蓄水期壩體內部主應力分布

表3　不同模型三維有限元計算壩體變形應力最大值表[19]

同時，汪明元等采用“南水”雙屈服面彈塑性模型對壩體填筑蓄水過程進行了模擬，得到蓄水期最大沉降約2.0 m，最大水平位移約0.50 m；清華大學錯誤!未找到引用源?；谛拚齂-G模型反饋參數(shù)計算得到蓄水期壩體最大沉降約2.39 m，最大水平位移0.55 m。

根據(jù)水布埡面板堆石壩觀測結果[19]，壩體內部最大沉降約為2.45 m，發(fā)生在最大斷面壩軸線偏下游側高程300.0 m處，最大水平位移約0.47 m，發(fā)生在壩軸線之間高程265.0 m處。

根據(jù)本文顆粒流方法的模擬結果，蓄水期壩體最大沉降2.24m，發(fā)生在壩軸線偏下游高程305 m高程處；上游最大水平位移0.27m，發(fā)生在高程287m處，下游最大水平位移0.45 m，發(fā)生在高程305 m處。

通過對比可以發(fā)現(xiàn)，運用顆粒流方法模擬水布埡面板堆石壩分層填筑過程，模型的宏觀變形特征與有限元結果比較接近，這說明利用顆粒流方法模擬超高堆石壩的分層分區(qū)填筑過程是切實可行的。通過與現(xiàn)場觀測對比可以看出，顆粒流模型的變形小于現(xiàn)場觀測值，這主要是由于本文尚未考慮流變的影響，對壩料碾壓考慮得也不充分，可以預見，若在顆粒流模型中考慮堆石料的流變，并對壩體碾壓進行更加細致的模擬，模擬結果將會更加真實地反映壩體的應力變形特征。

4　總結

本文針對顆粒流方法模擬超高堆石壩分層分區(qū)填筑過程的可行性進行了研究，結果表明，顆粒流方法模擬超高堆石壩的填筑過程是切實可行的，從與有限元結果對比可以看出，顆粒流方法的模擬結果是可靠的。在模擬過程中，采用并行算法使得單機所能容納的顆粒數(shù)量明顯增加，計算效率極大提高；同時，本文建立了一套快速有效的建模方法，可在復雜邊界區(qū)域快速生成密實的滿足級配要求的顆粒流模型。這兩個方面的改進與突破，使得顆粒流方法應用于模擬超高堆石壩分層分區(qū)填筑成為現(xiàn)實。

綜上所述，本文的研究成果對于顆粒流方法應用于超高堆石壩分層填筑模擬，具有重要的工程應用前景和理論價值。當然，本文的研究中尚未考慮堆石料的流變，對壩料的碾壓也考慮的不是很充分，對于這些的不足之處，將在后續(xù)的工作中作進一步深入研究，從而使得該模擬方法更加完善。

[1]徐澤平，鄧剛.高面板堆石壩的技術進展及超高面板堆石壩關鍵技術問題探討[J].水利學報，2008，39（10）:1226-1234.

[2]MaG，ZhouW，ChangXL.Modelingtheparticle breakage of rockfill materials with the cohesive crack model [J].Computers&Geotechnics，2014，61（3）:132-143.

[3]周偉，常曉林，周創(chuàng)兵，等.堆石體應力變形細觀模擬的隨機散粒體不連續(xù)變形模型及其應用 [J].巖石力學與工程學報，2009，28（03）:491-499.

[4]Deluzarche R，Cambou B.Discrete numerical modelling of rockfill dams[J].International journal for numerical and analytical methods in geomechanics，2006，30（11）:1075-1096.

[5]Alaei E，Mahboubi A.A discrete model for simulating shear strength and deformation behaviour of rockfill material, onsidering the particle breakage phenomenon[J].Granular Matter，2012，14（6）:707-717.

[6]李曉柱，劉洋，吳順川.堆石壩現(xiàn)場碾壓試驗與離散元數(shù)值分析[J].巖石力學與工程學報，2013，（S2）:3123-3133.

[7]劉洋，李曉柱，吳順川.多塊體形狀堆石體碾壓顆粒破碎數(shù)值模擬[J].巖土力學，2014（11）:3269-3280.

[8]Tran T H，Vénier R,Cambou B.Discrete modelling of rock-ageing in rockfill dams[J].Computers&Geotechnics，2009，36（1）:264-275.

[9]Silvani C，Désoyer T，Bonelli S.Discrete modelling of time-dependent rockfill behaviour[J].International Journal for Numerical&Analytical Methods in Geomechanics，2009，33（5）:665-685.

[10]邵磊，遲世春.堆石料單軸流變試驗的顆粒流模擬[J].巖土力學，2013，（5）:1487-1494.

[11]Cundall P A，Strack O D L.A discrete numerical model for granular assemblies[J].Geotechnique，1979，29（1）:47-65.

[12]柏樹田，崔亦昊.堆石的力學性質[J].水力發(fā)電學報，1997（3）: 21-30.

[13]孫其誠，王光謙.顆粒物質力學導論[M].北京:科學出版社，2009.

[14]Thomas P A,Bray J D.Capturing nonspherical shape ofgranularmediawithdiskclusters[J].Journalof GeotechnicalandGeoenvironmentalEngineering，1999，125（3）:169-178.

[15]Jiang M J,Yu H S,Harris D.A novel discrete model forgranularmaterialincorporatingrollingresistance[J]. Computers and Geotechnics，2005，32（5）:340-357.

[16]Ng T T,Dobry R.A non-linear numerical model for soil mechanics[J].International Journal for Numerical and Analytical Methods in Geomechanics，1992，16（4）:247-263.

[17]Itasca Consulting Group,Inc.（2014a）.PFC（particle flow code in 2 and 3 dimensions），version 5.0[User’s manual].Minneapolis，MN:ICG.

[18]楊啟貴，劉寧，孫役，等.水布埡面板堆石壩筑壩技術[M].北京:中國水利水電出版社，2010.

[19]汪明元，程展林，林紹忠，等.水布埡面板堆石壩的三維彈塑性數(shù)值分析研究[J].巖土力學，2004，（S2）:507-512.

黃小波（1981-），男，大學本科，工程師，從事水利工程設計工作，手機：15907492008。

（2016-02-24）