考慮顆粒形狀的面板堆石壩振動臺模型試驗DEM模擬

井向陽， 楊利福， 馬 剛， 周 偉

(1. 中國電建集團(tuán)成都勘測設(shè)計研究院有限公司，成都 610072； 2. 欽州市水利局，廣西 欽州 535000；3. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072)

粗粒料因取材方便、抗震性能好、施工快捷等優(yōu)點被大量用于巖土工程中，尤其在超大體積填筑工程中，粗粒料往往被設(shè)計者優(yōu)先考慮為工程主體填筑材料，如碎石樁、鐵(公)路路基、建筑基礎(chǔ)、堆石壩等。近年來，在我國西南強(qiáng)震區(qū)已建、在建或擬建一大批高土石壩(水布埡、天生橋一級、洪家渡、三板溪、糯扎渡、雙江口、兩河口、如美、馬吉等)，其中許多堆石壩壩高達(dá)到了200～300 m級，超過現(xiàn)行規(guī)范規(guī)定的壩高。由于壩高、庫大，一旦遭遇超強(qiáng)震而失事，不僅會造成重大經(jīng)濟(jì)損失，而且所形成的次生災(zāi)害將嚴(yán)重威脅下游人民生命財產(chǎn)安全。因此，有必要對堆石壩動力條件下的力學(xué)響應(yīng)進(jìn)行深入研究，為堆石壩抗震設(shè)計提供參考依據(jù)。

汶川地震以來，堆石壩震害引起了國內(nèi)外學(xué)者廣泛關(guān)注，并開展了大量試驗研究。在堆石壩動力響應(yīng)試驗研究方面，楊正權(quán)等[1-3]對雙江口高心墻堆石壩、兩河口高心墻堆石壩和猴子巖高面板堆石壩進(jìn)行了大型振動臺模型試驗研究；孔憲京等[4]進(jìn)行了地震作用下面板堆石壩面板錯臺振動臺模型試驗研究，并在土石壩壩坡模型振動臺破壞試驗基礎(chǔ)上，通過數(shù)值分析研究了壩坡動力特性[5]；Yuan等[6]制備了心墻堆石壩和面板堆石壩室內(nèi)縮尺模型，并通過振動臺上施加不同動力激勵研究了堆石體的動力響應(yīng)；Torisu等[7]通過大型振動臺模型試驗和室內(nèi)剪切試驗預(yù)測了堆石壩的動力殘余變形；徐澤平等[8]采用離心模型試驗研究了新疆察汗烏蘇砂礫石面板壩在施工期和蓄水期壩體動力響應(yīng)；程嵩等[9-10]通過動力離心模型試驗研究面板堆石壩震動響應(yīng)及變形規(guī)律；Kim等[11]通過動力離心模型試驗研究了心墻堆石壩和面板堆石壩的地震響應(yīng)。堆石壩震害現(xiàn)象觀察和振動臺動力模型試驗表明，堆石壩地震破壞始于下游坡面頂部的堆石體松動、滾落，最終導(dǎo)致壩頂坍塌。

在數(shù)值試驗研究方面，受計算機(jī)處理能力限制，堆石壩動力響應(yīng)較多采用連續(xù)介質(zhì)理論進(jìn)行數(shù)值分析。劉漢龍等[12]研究了云鵬心墻土石壩的地震易損性；張銳等[13]在研究了土石壩地震加速度動態(tài)分布、壩坡抗震能力；Arici等[14-15]在地基振動不同步、不同方向振動荷載對堆石壩動力響應(yīng)的影響做了大量研究。目前，也有一些學(xué)者采用非連續(xù)介質(zhì)力學(xué)方法對堆石壩動力響應(yīng)機(jī)制進(jìn)行研究?？讘椌┑萚16]結(jié)合振動臺物理試驗，采用非連續(xù)變形分析方法(Discontinuous Deformation Analysis,DDA)研究了面板堆石壩地震響應(yīng)特征及破壞過程；劉漢龍等[17]利用離散元軟件PFC2D模擬了土石壩振動臺模型試驗。

離散單元法(Discrete Element Method, DEM)是研究顆粒系統(tǒng)宏細(xì)觀力學(xué)特性的數(shù)值模擬方法之一，其能夠從細(xì)觀角度對顆粒體的各種力學(xué)特性進(jìn)行分析，為研究地震條件下堆石壩動力特性提供了一條有效途徑。本文借助有限元網(wǎng)格劃分技術(shù)，提出了復(fù)雜形狀顆粒DEM模擬方法，實現(xiàn)了包括了凹多邊形在內(nèi)的不規(guī)則隨機(jī)多邊形顆粒模擬；考慮堆石體顆粒復(fù)雜形狀，模擬了堆石壩模型振動臺試驗，研究了地震峰值加速度對面板堆石壩動力響應(yīng)的影響，并揭示了堆石壩細(xì)觀動力特性。

1 復(fù)雜形狀顆粒DEM模擬

在有限元-離散元耦合分析方法(Combined Finite-Discrete Element Method, FDEM)中，可以采用隨機(jī)多邊形(或多面體)來模擬顆粒形狀[18-19]，然而，當(dāng)顆粒形狀拓?fù)浣Y(jié)構(gòu)比較復(fù)雜(如凹多邊形)時，F(xiàn)DEM接觸檢索方法將無法準(zhǔn)確識別接觸顆粒。在DEM模擬中，顆粒的形狀多為二維圓盤或三維圓球，無法考慮顆粒之間的嵌入咬合作用，數(shù)值模擬結(jié)果與實際顆粒的力學(xué)性質(zhì)存在一定差異。為此，許多學(xué)者將形狀簡單的圓盤或圓球通過一定的方式組合在一起形成“顆粒簇”(cluster)或“超級顆?！?clump)來反映顆粒的不規(guī)則形狀[20-21]。此外，在DEM方法中也有以橢圓或者橢球為基本顆粒的，通過改變橢圓的長短軸之比來模擬不同顆粒形狀[22]。上述DEM研究只能模擬一些形狀簡單且規(guī)則的顆粒，對于巖土工程中的復(fù)雜形狀顆粒往往需要借助圖像獲取技術(shù)才能實現(xiàn)真實顆粒形狀的再現(xiàn)。Wang等[23-24]采用CT掃描技術(shù)獲得顆粒表面信息后，在離散元模擬中用“超級顆?！敝貥?gòu)了真實顆粒形狀，但CT掃描方法能提供的顆粒形狀很有限，無法反映真實巖土顆粒形狀的隨機(jī)性和多樣性。

常曉林等[25]在橢球面上隨機(jī)布點后生成三維凸多面體顆粒，通過改變橢球長短軸之比和顆粒點數(shù)來控制顆粒形狀。本文提出的隨機(jī)多邊形模擬方式可以生成包括凹多邊形在內(nèi)的復(fù)雜顆粒形狀。復(fù)雜多邊形顆粒由邊數(shù)、極角、圓半徑表征確定，如圖1所示, 具體生成方法包括六個步驟。

圖1 復(fù)雜多邊形顆粒的表征Fig.1 Characterization of polygon particles

步驟1在指定區(qū)域內(nèi)按照給定級配曲線生成隨機(jī)圓形顆粒，此時，顆粒圓心位置(x0，y0)和半徑r已知。

步驟2指定最大邊數(shù)和最小邊數(shù)，生成多邊形隨機(jī)邊數(shù)。

n=nmin+int[(nmax-nmin)rand]

(1)

式中:n為生成隨機(jī)多邊形的邊數(shù)；nmin為最小邊數(shù)； int為取整函數(shù)；nmax為最大邊數(shù)； rand為0～1的隨機(jī)數(shù)。

θk=2π[1+(2rand-1)δ]/n

(2)

(3)

步驟4獲取多邊形頂點坐標(biāo)。

考慮到本文數(shù)據(jù)結(jié)果表征的不足，課題下一步的重點是修正生態(tài)足跡模型，篩選恰當(dāng)變量指標(biāo)，使數(shù)據(jù)結(jié)果更為精確，希期為寧德市經(jīng)濟(jì)、社會、生態(tài)可持續(xù)發(fā)展提供科學(xué)的理論參考和數(shù)據(jù)借鑒。

(4)

圖2、圖3分別為λ=0和λ=0.2時生成的隨機(jī)多邊形。 當(dāng)λ=0時，所生成的隨機(jī)多邊形為圓內(nèi)接凸多邊形，隨著邊數(shù)的增加，隨機(jī)多邊形的形狀越接近圓。 當(dāng)λ≠0時，本文提出的隨機(jī)多邊形生成方法可以生成凸多邊形和凹多邊形，當(dāng)邊數(shù)較少時，所生成的多邊形為凸多邊形但與圓無內(nèi)接關(guān)系，隨著邊數(shù)的增加，顆粒形狀就越復(fù)雜，與巖土工程中的真實顆粒形狀越接近。

圖2 生成的凸隨機(jī)多邊形(λ=0)Fig.2 Polygon particles considering λ=0

圖3 生成的復(fù)雜隨機(jī)多邊形(λ=0.2)Fig.3 Polygon particles considering λ=0.2

步驟5根據(jù)以上復(fù)雜多邊形的邊數(shù)、極角、圓半徑等幾何信息，在有限元軟件ANSYS中生成多邊形面，并對各個面進(jìn)行網(wǎng)格劃分。

步驟6通過在有限元軟件ANSYS中獲取的網(wǎng)格形心位置及網(wǎng)格面積等信息，在離散元軟件PFC2D中用等面積同形心圓盤替換各個網(wǎng)格，并將同一多邊形內(nèi)的圓盤采用clump命令組構(gòu)成為“超級顆?！?，以此模擬復(fù)雜隨機(jī)多邊形，如圖4所示。在PFC2D中clump作為一個獨立的剛體參與計算，clump內(nèi)部的顆粒重疊將會被忽略。

圖4 有限元網(wǎng)格與DEM計算顆粒Fig.4 Finite element meshes and DEM clumps

2 堆石壩振動臺試驗DEM模擬

2.1 計算模型與參數(shù)

模型壩體斷面為梯形，如圖5所示，模型壩壩高取為0.6 m，上游水位0.55 m，上下游邊坡均為1∶1.4，壩頂寬0.08 m，面板厚度6 mm，壩體采用粒徑0.008～0.024 m的均勻級配，初始孔隙率設(shè)定為0.15。

圖5 振動臺試驗DEM模型Fig.5 Discrete element model of shaking table model test

在振動過程中顆粒間的相互作用會消耗部分能量，為了模擬系統(tǒng)中能量的耗散，本文引入黏滯阻尼與局部阻尼作為耗能機(jī)制，局部阻尼系數(shù)取0.1，法向和切向黏滯阻尼系數(shù)均取為0.2[26]。數(shù)值模擬中涉及到的參數(shù)有: 顆粒密度取為2 650 kg/m3，顆粒間摩擦因數(shù)取0.5，多邊形顆粒間的接觸采用修正的線性接觸模型，即kn=k0×r，kn為法向接觸剛度，r為多邊形顆粒等效半徑，k0取1×105kN/m2， 切向剛度ks=kn。 采用汶川波作為輸入地震荷載，見圖6。為了減小顆粒振動分離對數(shù)值模擬的影響，所施加的地震荷載為水平方向，振動持續(xù)15 s，并設(shè)置了峰值加速度不同的三組輸入地震荷載作為對比試驗，分別為0.12g，0.26g和0.445g。通過時間積分，得到不同峰值加速度下的速度曲線，將經(jīng)過基線校正的速度曲線作為壩體底部墻體的水平運(yùn)動速度。

圖6 汶川地震加速度記錄曲線Fig.6 Acceleration of Wenchuan earthquake time history

2.2 結(jié)果分析

圖7和圖8 所示為不同峰值加速度下面板堆石壩震后變形模擬結(jié)果。由圖可知，當(dāng)峰值加速度為0.12g時，下游壩坡表層顆粒近乎平行原壩坡面滑動；隨著峰值加速度的增大，下游壩坡顆粒滑移面角度變大，壩坡坍塌明顯，這與文獻(xiàn)[27]一致。受面板約束和上游水壓力作用影響，上游壩坡位移較小，但面板與底板出現(xiàn)了滑動，最大位移出現(xiàn)在坡腳附近，表明面板止水部位在地震作用下容易受損。在豎直方向上，不同加速度下壩體內(nèi)部顆粒均出現(xiàn)了向下位移，且隨著加速度的增大，顆粒沉降越明顯。

圖7 震后堆石壩水平向變形Fig.7 Horizontal displacement of CFRD after the earthquake

圖8 震后堆石壩豎直向變形Fig.8 Vertical displacement of CFRD after the earthquake

由地震前后壩體輪廓變化可見，地震導(dǎo)致了壩頂?shù)奶蜕舷掠螇纹孪蛲夤某觯瑸榱硕勘容^峰值加速度對壩體變形的響，統(tǒng)計了不同峰值加速度下面板堆石壩模型震后輪廓特性，如圖9。震前，上下游壩坡斜率均為0.714 3，震后上下游坡腳向外擴(kuò)張，壩坡斜率減小，且下游壩坡比上游壩坡更緩，隨著峰值加速度的增加，下游壩坡斜率從0.623 9變?yōu)?.556 9，壩頂高程從0.559 2 m減小為0.537 8 m。在0.12g峰值加速度汶川波作用壩體斷面面積從0.552 m2最大增大到了

圖9 震后堆石壩輪廓線Fig.9 Skeleton for CFRD after the earthquake

0.567 2 m2，峰值加速度越大，壩體剪脹越明顯。模擬結(jié)果與紫平鋪面板壩汶川地震后實際變形不一致[28]，但與室內(nèi)振動臺模型試驗一致。這主要是堆石體在低圍壓下剪脹變形而高圍壓下顆粒破碎抑制剪脹，振動臺試驗未考慮堆石體顆粒破碎效應(yīng)，模擬結(jié)果(包括物理試驗和數(shù)值試驗)與高堆石壩震后變形不一致。

圖10為不同峰值面板堆石壩震后力鏈分布圖。圖中接觸力大于平均力且被接觸力較小顆粒圍繞的柱狀結(jié)構(gòu)顆粒體系的接觸被定義為強(qiáng)力鏈[29-30]。由圖可知，壩體蓄水后強(qiáng)力鏈基本呈對稱分布。地震過程中，靜力平衡下的強(qiáng)力鏈開始崩塌，并伴有新的強(qiáng)力鏈形成。震后壩體強(qiáng)力鏈總體分布向下游傾斜。

顆粒材料堆積特性與其細(xì)觀組構(gòu)密切相關(guān)。Rothenburg等[31]采用傅里葉函數(shù)來擬合強(qiáng)力鏈法向接觸力的角域分布，其表達(dá)式為

fn(θ)=f0[1+ancos 2(θ-θa)]

(5)

式中:f0為顆粒法向接觸力平均值；θ為顆粒接觸法向與試樣水平面的夾角；θa分別為法向接觸力各向異性

圖10 震后堆石壩力鏈分布圖Fig.10 Force chain distribution of CFRD after the earthquake

的主方向；an為反映法向接觸力各向異性程度的傅里葉系數(shù)。

圖11給出了不同峰值加速度下堆石顆粒間接觸法向力各向異性分布圖和相應(yīng)的傅里葉函數(shù)擬合結(jié)果。圖中每6°一個區(qū)間，統(tǒng)計接觸法向落入角度區(qū)間內(nèi)的強(qiáng)力鏈平均法向接觸力。震前，壩體強(qiáng)力鏈玫瑰圖呈“花生狀”分布，法向接觸力主軸方向接近于豎直方向。震后，壩體強(qiáng)力鏈玫瑰圖呈“橢圓狀”分布，法向接觸力主軸方向減小，即向下游偏轉(zhuǎn)，與圖10結(jié)果相似。與震前相比，震后法向接觸力各向異性變小，但隨著峰值加速度的增大，顆粒間接觸法向接觸力各向異性程度逐漸增強(qiáng)。

圖12為地震過程中面板堆石壩內(nèi)部強(qiáng)力鏈數(shù)演化圖。由圖12可知，地震過程中，壩體內(nèi)部強(qiáng)力鏈的失效多于新的強(qiáng)力鏈形成，且隨著峰值加速度的增大，強(qiáng)力鏈?zhǔn)?qiáng)度越大。為了研究顆粒接觸處的摩擦特性，定義摩擦激勵指標(biāo)Im=|ft|/(fntanφu)，IM為顆粒集合體中所有接觸的平均摩擦激勵指標(biāo)，圖13為不同峰值加速度下的壩體平均激勵指標(biāo)IM在地震程中的演化過程。地震初期，平均摩擦激勵有所減小。隨著峰值加速度增加，平均摩擦激勵增大，但峰值加速度為0.26g和0.445g時，壩體平均摩擦激勵相差很小，在0.62處波動，即在不同地震荷載下，壩體破壞時所激發(fā)的平均摩擦激勵是一個定值。說明面板堆石壩在相同初始組構(gòu)條件下抵抗地震破壞的能力是一定的，與加載條件無關(guān)。

圖11 震后堆石壩法向接觸力各向異性分布圖Fig.11 Normal contact forces anisotropy distribution for CFRD after the earthquake

圖12 地震過程中堆石壩內(nèi)部強(qiáng)力鏈數(shù)Fig.12 Number of strong force chain for CFRD during the earthquake

圖13 地震過程中堆石壩摩擦激勵演化規(guī)律Fig.13 Evolutions of average friction mobilization for CFRD during the earthquake

3 結(jié) 論

本文提出了復(fù)雜形狀顆粒DEM模擬方法，模擬了面板堆石壩振動臺模型試驗，從細(xì)觀角度揭示了面板堆石壩動力特性，主要結(jié)論如下:

(1) 提出了復(fù)雜形狀顆粒DEM模擬方法，在PFC2D中實現(xiàn)了包括凹多邊形在內(nèi)的復(fù)雜形狀顆粒模擬，通過參數(shù)控制使得生成的顆粒形狀更加真實。

(2) 峰值加速度較小時，下游壩坡表層顆粒沿原壩坡面方向滑動；隨著峰值加速度的增大，下游壩坡顆粒滑移面角度變大，壩坡坍塌明顯。面板與底板出現(xiàn)了滑動。不考慮顆粒破碎時，地震導(dǎo)致上下游坡腳向外擴(kuò)張，壩頂沉降，但壩體剪脹。

(3) 震前壩體強(qiáng)力鏈玫瑰圖呈“花生狀”分布，法向接觸力主軸方向接近于豎直方向。震后，壩體強(qiáng)力鏈玫瑰圖呈“橢圓狀”分布，法向接觸力主軸方向向下游偏轉(zhuǎn)。隨著峰值加速度的增大，顆粒間接觸法向接觸力各向異性程度逐漸增強(qiáng)。

(4) 地震過程中平均摩擦激勵先減小后增大到一個穩(wěn)定水平。隨著峰值加速度增加，平均摩擦激勵增大，但在較高加速度下，壩體最大平均摩擦激變化不大，說明堆石壩在相同初始組構(gòu)條件下抵抗地震破壞的能力是一定的，與加載條件無關(guān)。