基于有限元法的心墻礫石土壓實特性分析

崔 博，王毓杰，王佳俊，余 佳

(天津大學水利工程仿真與安全國家重點實驗室，天津300350)

0 引 言

礫石土料經過壓實后有較高的密度和強度以及較好的防滲性能，目前作為填筑材料廣泛應用于高心墻土石壩中[1]。嚴格控制現場施工工藝以保證礫石土壓實質量，其中，碾壓作業是心墻礫石土施工的關鍵環節，開展壓實特性(物理特性和變形特性)研究對于控制碾壓施工質量有著非常關鍵的作用[2-3]。在碾壓過程中，高振狀態使礫石土產生較大的塑性變形，碾輪與礫石土的接觸面不斷變化，導致礫石土的壓實特性也發生了顯著變化，準確把握礫石土壓實特性的變化對控制心墻礫石土的壓實質量有著重要的研究意義。

碾壓過程中，由于振動荷載下土體塑性變形的滯后性和碾輪與土體接觸面的不確定性，碾輪與土體的動力相互作用表現出材料非線性、幾何非線性和接觸非線性[4]。采用傳統的靜力學或者動力學理論分析時會很復雜繁瑣，且難以得到精確的解析結果，因此需采用數值模擬方法對這類問題進行數值分析。數值模擬是揭示土體壓實機理的重要方法，目前已成功應用于瀝青路面[3，5-6]、碾壓混凝土[7]、黃土[8]、堆石體[9]、宕渣[10]等領域[11-14]的壓實研究。從以上研究可知，利用數值模擬方法通過分析土體的壓實特性、力學響應以及各種因素對壓實效果的影響等方面，能夠有效地研究土體在碾壓過程中的壓實行為，不過尚未發現數值模擬方法在礫石土壓實特性的相關研究。

礫石土是由一定級配的巖石顆粒和土料組成的無凝聚性土石混合體[15]。喬蘭、Xiao等[16-17]對礫石土進行室內擊實試驗，分析了在擊實作用下礫石土的顆粒破碎規律，并分析了礫石土壓實特性的變化規律；李朝政、馮瑞玲等[18-19]通過室內擊實試驗和現場碾壓試驗結合確定礫石土的碾壓參數，以此控制工程施工質量。同時，近年來一些學者也利用數值模擬方法來研究土石混合體的物理力學特性。王鵬飛等[20]對土體混合體進行數值模擬，從含石量、塊石空間分布方面來研究其力學特性的變化；崔博等[21]基于離散元法進行土石混合體的三軸數值模擬試驗，對其力學性質進行研究。從上述分析可知，礫石土的壓實特性大多數是通過室內擊實試驗開展相關研究，并且通過現場碾壓試驗確定碾壓參數；而數值模擬方面的研究是通過結合三軸壓縮與離散元法研究其力學性質。

國內外研究中，尚未發現將振動凸塊碾與心墻礫石土結合起來作為耦合體系研究碾壓過程中礫石土體壓實特性的變化；同時，考慮到現場試驗有耗資大、試驗工作量大、周期長及數據離散型性高等缺陷，因此有必要在傳統試驗工作的基礎上，引入數值模擬的手段來彌補室內和現場試驗方面的不足。由于凸塊碾與礫石土相互作用時存在復雜的非線性關系以及礫石土有非均質性和各向異性等復雜工程性質，目前沒有精確的本構關系來描述其力學特性。ABAQUS作為有限元計算分析軟件的代表，具有良好的非線性解析功能和收斂性，常用于求解各種復雜材料的非線性問題，而且具有能真實反映土體復雜性狀的本構模型，易于處理礫石土這類土石混合料的非均質性和各向異性問題。

本文基于有限元法，構建“凸塊碾-礫石土”三維有限元模型，以研究心墻礫石土的壓實特性。首先，根據現場碾壓的振動凸塊碾尺寸，將振動凸塊碾簡化為碾輪并建立其有限元模型；其次，確定礫石土的模型尺寸、本構模型以及基本參數以構建礫石土的三維有限元模型；再者，通過“凸塊碾與礫石土”模型的接觸分析，模型作用荷載的施加進行模型的仿真計算，并與實測值對比證明該模型的可行性；最后，從沉降量、豎向應力幅值、密度等方面研究碾壓過程中礫石土的壓實特性，分析礫石土的應力分布規律，為心墻礫石土的壓實質量研究提供科學的理論依據。

1 “凸塊碾-礫石土”三維有限元模型

1.1 建立振動凸塊碾模型

振動凸塊碾是一種作用力組合式的壓實機械，是將振動平碾、凸塊的壓實原理進行復合，在其碾壓滾筒上有多排矩形塔狀凸塊，按滾筒軸線交錯布置。振動凸塊碾在碾壓過程中不僅借助碾輪重力與激振力的作用，同時也借助凸塊插入土體時使土體受到擠壓和揉搓的聯合作用，從而使土體達到良好的壓實效果。本文研究的凸塊碾機械參數、礫石土模型參數、碾壓試驗數據資料均依托于西南某高心墻堆石壩工程，其中凸塊碾模型為心墻碾壓中的三一重工26 t自行式振動凸塊碾，其機械參數如表1所示。

表1 凸塊碾機械參數

由于振動凸塊碾碾壓心墻礫石土時，主要是其碾輪部分對礫石土有壓實作用，綜合前人振動壓實數值模擬的相關研究[9-10]，加之考慮到其結構的復雜性，也為了減少計算量，將振動凸塊碾簡化為碾輪，利用ABAQUS建立光滑碾輪模型和凸塊模型。將建立好的碾輪模型和凸塊模型導入到SolidWorks軟件進行排列裝配，再將其導入到ABAQUS得到凸塊碾模型如圖1所示。由于碾輪相對土體而言剛度很大，而且也并非研究對象，將其視為剛體，不參與變形計算。

圖1 振動凸塊碾模型示意

1.2 建立礫石土模型

礫石土體屬于半空間結構，前后左右及下方是無限延伸的，在進行有限元分析時可取長方體礫石土塊來代替理論上的半空間結構[22]。在接觸表面除外的5個面(四個側面和底面)上施加固定邊界約束，并將礫石土有限元模型分為2層，上層為碾壓層，厚30 cm；下層為已壓實的下墊層，厚25 cm，與現場施工情況吻合。

根據現場碾壓作業的施工情況，凸塊碾在碾壓心墻礫石土時行駛速度要嚴格控制在2～3 km/h，碾輪的直徑是1.7 m，當凸塊碾以最低速度2 km/h滾動行駛時，滾動一周所需的時間最長，經計算最大周期9.048 s。

為了模擬礫石土在整個周期內的碾壓過程，必須使碾壓一遍的仿真時間大于最大周期9.048 s，因此本文模型在模擬心墻礫石土的碾壓施工時，碾壓一遍仿真時間選取為10 s。由于凸塊碾的碾輪寬度為2.17 m，可將礫石土模型寬度設為2.5 m，本文仿真中凸塊碾的最高行駛速度為3 km/h，仿真時間為10 s，則模型所需最短長度8.3 m，將模型長度設為10 m。本文的礫石土模型分為兩層，上層為碾壓層，模型長10 m、寬2.5 m、高0.3 m，下層為壓實后的下墊層，模型長10 m、寬2.5 m、高0.25 m，得到礫石土模型如圖2所示。

圖2 礫石土模型示意

針對礫石土本構模型的選取，由徐鼎平等[23]的研究可知Mohr-Coulomb模型能夠有效表達礫石土在受壓作用下的應力-應變關系，因此本文選取ABAQUS里面的Mohr-Coulomb模型進行模擬，模型初始的各項計算參數見表2。

表2 礫石土模型參數

1.3 模型接觸作用分析

礫石土在凸塊碾的滾動過程中受到碾輪凸塊的擠壓和揉搓作用，加上凸塊碾的振動作用使其力學性能不斷發生變化，加上凸塊碾與礫石土在碾壓過程中接觸面積的變化導致二者之間存在復雜的非線性問題。ABAQUS中解決非線性線性問題的求解器包括ABAQUS/Standard和ABAQUS/Explicit，對于高速復雜的非線性行為應用ABAQUS/Explicit求解器較高效準確。凸塊碾與礫石土存在復雜的接觸行為且凸塊碾的振動頻率較高，故本文選擇ABAQUS/Explicit求解器進行有限元計算分析。對于存在復雜拓撲關系的模型之間的接觸分析，考慮到碾壓過程中凸塊碾對礫石土的擠壓和揉搓作用導致接觸面不斷發生變化，從而產生非線性接觸，因此需要使用通用接觸來施加接觸。通用接觸可以通過軟件自動分析模型來定義模型之間的接觸關系，從而節約時間成本。

接觸屬性包括3方面：法向作用、切向作用和阻尼作用。法向作用是在“凸塊碾-礫石土”模型中，當凸塊碾與礫石土不發生接觸行為時，二者之間不存在拉應力，故選取硬接觸；切向作用為利用靜-動摩擦指數衰減形式來定義接觸面的切向摩擦，該摩擦模型還需定義彈性滑移剛度，假設碾壓過程中不存在滑移現象，故彈性滑移剛度設為無限大；阻尼作用為接觸阻尼力在法向方向上與滑移率成比例，在切向方向與彈性滑移率成比例，切向阻尼的定義是通過定義其和法向阻尼系數的比值來定義的，ABAQUS/Explicit中默認該比值大小為1。

1.4 模型作用力

振動碾壓過程中，“凸塊碾-礫石土”模型施加的荷載為

P=G+F

(1)

F=F0sin(ωt)

(2)

式中，P為礫石土的作用力；G為機械自重；F0為激振力；F為F0的分力；ω為振動頻率；t為時間。

“凸塊碾-礫石土”模型施加作用力后如圖3所示。

圖3 模型作用力示意

1.5 “凸塊碾-礫石土”有限元模型驗證

在振動壓實過程中，隨著土體不斷地被壓實，土體的相關特性參數也在不斷發生變化。根據前人關于振動壓實數值模擬的研究，有學者[22]通過改變不同遍數下的模型參數來進行模擬碾壓過程，也有學者[3，13]考慮到壓實前與壓實后的模型參數變化不大，假設模型參數不變的情況下采用初始的模型參數來進行模擬碾壓過程，均取得了不錯的效果。本文采取這2種方式對“凸塊碾-礫石土”模型進行模擬，一方面能驗證該模型的可行性，另一方面為后續的研究奠定了基礎，結合現場試驗數據資料得到相應的模型參數見表3。將實測的沉降量與模擬計算得到的沉降量進行對比，對比結果見圖4，其中模擬曲線1為模型參數改變時沉降量的模擬結果曲線，模擬曲線2為模型參數不變時沉降量的模擬結果曲線。

表3 模型參數

圖4 模擬結果與實測沉降曲線對比

由圖4可知，3條曲線的沉降量隨碾壓遍數的變化基本上保持一致，即初始狀態下沉降量增加較快，隨著碾壓遍數的增加，沉降量逐漸趨于穩定。可以明顯看出，模型參數不變時，與實測的沉降量曲線的變化規律較為吻合，而模型參數改變時的沉降量與實際沉降量數值上相差較大。模擬曲線1與實測值的總沉降誤差為1 cm，此時的總沉降的相對誤差為14.7%，最大相對誤差為16%，最小相對誤差為3.67%；模擬曲線2與實測值的總沉降誤差為0.18 cm，此時總沉降的相對誤差為1.47%，最大相對誤差為16%，均滿足數值模擬的要求[7]，說明基于有限元法構建的“凸塊碾-礫石土”有限元模型研究心墻礫石土的壓實過程是可行的。模擬曲線1與實際沉降曲線相差較大是由于礫石土壓實后與壓實前的模型參數變化并不是很大，改變模型參數在高強度的激振力作用下模型進行有限元數值計算時，礫石土仍會產生少量豎向位移，不同模型參數下產生的豎向位移不斷疊加，這才導致了模擬曲線1相較于模擬曲線2與實測沉降的總沉降誤差較大；而在模型參數不變時，凸塊碾的振動作用使礫石土在前期產生較大變形，隨著碾壓遍數的增加，變形會逐漸減小并趨于穩定，與實測沉降的變化規律也比較相符。綜上，本文假定模型的初始參數不變進行“凸塊碾-礫石土”有限元模型的數值模擬。

2 基于“凸塊碾-礫石土”模型的礫石土壓實特性分析

2.1 沉降量分析

為了更加細致地了解礫石土沉降量在碾壓過程中的變化情況，本文將“凸塊碾-礫石土”模型網格劃分出來的單元結點進行編號，通過分析壓實軌跡中心處某觀測點的沉降量在各結點對應的不同深度隨時間的變化曲線，以了解不同深度處的沉降量變化情況。圖5是對土體中心橫向剖切得到的結點圖，在圖中選取該觀測點沿碾壓層土體深度方向(即Y軸方向)一條直線上的4個結點，分別是13 881(Y=0 cm)、35 088(Y=10 cm)、34 889(Y=20 cm)、6 120(Y=30 cm)。圖6為各結點豎向位移的時程曲線。圖7為各結點沉降量隨碾壓遍數的變化曲線。

圖5 橫向剖面結點

圖6 各結點沉降量時程曲線

圖7 不同深度的沉降量隨碾壓遍數變化曲線

由圖6可知，該處土體在接近5 s時受到振動凸塊碾的沖擊作用。由于碾壓一遍的周期為10 s，隨后時間每增加10 s左右都要受到振動凸塊碾的沖擊作用，該時刻結點13 881處礫石土受到的沖擊效果最強烈，結點13 881處土體的沉降量最大，在隨后0.5 s左右的時間內出現少許回彈，此后該結點土體在周期內的沉降量基本上逐步趨于穩定。

綜合圖6和圖7可以看出，結點13 881處被碾壓第1遍后(5 s左右)、第2遍后(15 s左右)、第3遍后(25 s左右)和第4遍后(35 s左右)，該結點土體沉降量上升幅度比較大，隨后隨著碾壓遍數的增加，當碾壓到第7遍時，累積的沉降量逐步穩定在6.9 cm左右，且第8遍與第7遍的沉降差小于0.1 cm，說明此時土體已經達到較高的密實狀態；結點35 088沉降量的變化規律和表層結點13 881一樣，不過小于土體表面發生的沉降量，最終的沉降量在1.6 cm左右；而結點36 459和結點6 695隨著碾壓遍數的增加，沉降量也只是毫米級(即在10 mm以內)，尤其是結點6 695由于離振動源最遠，沉降量已經趨近于0。

此外，由礫石土的表面沉降量變化規律可以看出，在心墻礫石土的碾壓過程中，前期礫石土沉降量的增加速率較快，隨著碾壓遍數的增加，沉降量的增加速率慢慢減小，土體并逐漸趨于穩定狀態；從圖中也可以明顯的看出礫石土的主要沉降量發生在土體表層部分(0～10cm)，土體產生的沉降量沿著深度方向會越小，受到沖擊作用后土體回彈的位移也在沿著深度方向不斷變小，這說明碾壓過程中凸塊碾產生的壓實能量傳遞到土體內部時，會沿著深度方向不斷減少，從而導致礫石土的沉降量也會沿著深度方向逐漸遞減。

2.2 豎向應力幅值分析

振動凸塊碾和心墻礫石土在碾壓過程中產生了復雜的動態響應問題，從而使礫石土產生動壓應力信號，動壓應力信號包含相應的壓實信息。采用礫石土的豎向應力幅值(礫石土受到振動沖擊力時產生的豎向動壓應力最大值)來揭示碾壓過程中礫石土壓實效果的變化規律。圖8、圖9分別為豎向應力幅值隨碾壓遍數、土體深度的變化曲線。

圖8 豎向應力幅值-遍數曲線

圖9 豎向應力幅值隨土體深度的變化曲線

由圖8可知，豎向應力幅值與沉降量隨碾壓遍數的變化趨勢是一致的，即激振力開始作用于礫石土時，礫石土表面產生的豎向應力幅值較小，隨著碾壓遍數的增加，表面豎向應力幅值也在不斷增大，碾壓7遍后，表面豎向應力幅值逐漸穩定在1.3 MPa左右。圖9可以很直觀地反映出豎向應力幅值隨深度變化關系，其中0～30 cm為碾壓層，30～55 cm為下墊層。當礫石土受到振動沖擊作用時，土體表面的豎向應力幅值最大，土體深度從0～55 cm的豎向應力幅值都有不同程度的衰減，碾壓層土體(0～30 cm)的豎向應力幅值沿著深度方向衰減幅度較大，而下墊層土體(30～55 cm)的豎向應力幅值沿著深度方向衰減幅度較小。

礫石土在受到激振力的振動作用時產生的壓應力包括彈性反力和慣性力，土體的密實度越大，彈性反力也越大，即豎向應力幅值越大，故豎向應力幅值會隨著碾壓遍數的增加而增加，當達到一定的碾壓遍數時，礫石土表面的豎向應力幅值逐漸趨于穩定。由于凸塊碾產生的壓實能量是沿著深度方向逐漸遞減的，導致豎向應力幅值沿著深度方向逐漸遞減，豎向應力幅值在下墊層的大小是比較明顯的，這說明凸塊碾的激振力對礫石土有較大的影響范圍。

2.3 密度值分析

由于ABAQUS軟件沒有密度的量值輸出，所以必須通過換算進行密度的求解與分析。在碾壓過程中，碾輪與礫石土的作用時間很短，礫石土呈現彈性為主塑形為輔的特征[24]，由此可以得到沖壓前后土體應變ε與密度ρ的關系為

(3)

式中，ρ0為土體的初始密度；ρ為沖壓后的土體密度。

圖10 不同深度的密度-遍數曲線

由圖10可知，礫石土在一定深度(0～20 cm)下的密度值均隨著碾壓遍數的增加不斷變大，當碾壓7遍后，密度值的變化量均很小并逐漸趨于一個穩定值，此時礫石土表面(Y=0 cm左右)的密度值最大，約為2.50 g/cm3；密度值沿著深度方向逐漸遞減，Y=10 cm和Y=20 cm處最終的密度穩定值分別為2.30、2.04 g/cm3，碾壓層底部的密度值幾乎為0。與此同時，深度0～20 cm范圍內的土體均在碾壓第7遍時密度值趨于一個穩定值，這說明此時凸塊碾振動產生的壓實能量能夠有效傳遞到礫石土的一定深度，使礫石土在不同深度均能達到較好的壓實效果。

2.4 豎向應力幅值分布特性分析

從碾輪的行駛方向以及寬度方向，分析豎向應力幅值的分布特性。圖11是激振力作用在礫石土時，在壓實軌跡中心上沿行駛方向均勻選取10個觀測點的豎向應力幅值變化曲線。圖12是激振力作用于某一區域時，沿輪寬方向選取11個觀測點(觀測點6在碾輪寬度方向的中心線上，觀測點1～5和7～11對稱分布于碾輪中心線兩側)的豎向應力幅值變化曲線。

圖11 碾輪行駛方向的豎向應力幅值變化曲線

圖12 碾輪輪寬方向的豎向應力幅值變化曲線

由圖11可知，在激振力的作用下礫石土沿著凸塊碾行駛方向的豎向應力幅值相差較小，振動壓實時在長度方向保證壓實效果的一致性和均勻性。由圖12可知，凸塊碾輪寬中心線上礫石土產生的豎向應力幅值最大，沿輪寬方向兩側的豎向應力幅值較小且呈近似為拋物線分布，這說明土體沿輪寬方向的壓實均勻性不理想，即中間壓實效果較好，而兩邊壓實效果較差，從理論方向上解釋了碾壓作業中兩兩相鄰的條帶之間必須有一定的搭接寬度以保證每一個條帶的壓實質量。

3 結 語

本文利用有限元法，將凸塊碾與礫石土作為耦合體系來研究心墻礫石土在碾壓過程中的壓實特性，利用有限元軟件ABAQUS構建了“凸塊碾-礫石土”三維有限元模型，模擬了心墻礫石土碾壓的整個過程，并將模擬計算的仿真值與實測值作了對比，結果表明該模型具有較好的可靠性。

基于數值模擬結果，分析了心墻礫石土在碾壓過程中的壓實特性，得到以下結論：

(1)礫石土在碾壓過程中受到沖擊作用時一開始會產生較大的塑性變形即土體有較大的沉降量，隨著碾壓遍數的增加土體沉降量逐漸變小并趨于穩定；礫石土表面受到的沖擊效果最好，沉降量也最大，隨著土體深度的增加，沉降量也會逐漸減小，當到達土體底部時，其沉降量幾乎為0。

(2)豎向應力幅值和密度與沉降量變化規律具有一致性，均是達到一定碾壓遍數后逐步趨于穩定；各項數值均是隨著深度的增加而減小，且由豎向應力幅值隨深度的變化可知碾輪的沖擊力對土體有較大的影響范圍。

(3)碾壓過程中礫石土產生的豎向應力幅值相差，即礫石土的壓實效果沿長度方向具有一致性和均勻性；凸塊碾作用下的礫石土中間壓實效果較好，兩邊壓實效果較差，豎向應力幅值近似為拋物線分布。