基于DEM 模擬的顆粒材料熱固結(jié)試驗分析*

應(yīng)濤濤，魏良針，楊 磊，夏念興，龔禮岳

（1.溫州市鐵路與軌道交通投資集團有限公司，浙江 溫州 325000；2.溫州大學 建筑工程學院，浙江 溫州 325000）

基于節(jié)能環(huán)保可持續(xù)的發(fā)展原則，如能源樁[1]和地源熱泵[2]等與溫度相關(guān)的能源工程[3-5]開展得越來越多，溫度對能源結(jié)構(gòu)周邊土體的影響已經(jīng)成為工程在使用過程中不可忽略的重要因素。Sang 等[6]分別在室內(nèi)和高溫下進行了飽和粘土的熱排水固結(jié)試驗，發(fā)現(xiàn)溫度可以加快排水固結(jié)的固結(jié)速度。Xiong 等[7]開發(fā)了一種溫度控制的三軸試驗裝置，研究了飽和軟粉質(zhì)粘土在不同溫度下的動態(tài)不排水特征，發(fā)現(xiàn)溫度對軟粘土的累積塑性變形、孔隙水壓力等特性有很大的影響。此外，溫度土體的強度[8-10]、空隙特征[11]、變形[12]等特性也有較大的影響。土體特性的改變必然會造成土體與結(jié)構(gòu)接觸面特性的變化。Maghsoodi 等[13]建立了非等溫土-結(jié)構(gòu)界面模型，討論了溫度對粘土-結(jié)構(gòu)界面空隙率的影響。Donna 等[14]開發(fā)了一種直接剪切裝置，研究了不同溫度下樁-土界面的響應(yīng)，發(fā)現(xiàn)粘土-混凝土界面的強度隨溫度的升高而增加。此外，研究還表明溫度對土-結(jié)構(gòu)界面的接觸力[15]、應(yīng)力-應(yīng)變[16]和抗剪強度[17-19]等特性都有影響。而土-結(jié)構(gòu)接觸面特性的變化勢必會引起結(jié)構(gòu)的不均勻沉降甚至失穩(wěn)，所以需進一步深入土體的溫度效應(yīng)的研究。但是，上述試驗中都只是分析了溫度對土體在宏觀上的如強度、應(yīng)力和應(yīng)變等影響，缺少了如各項異性、內(nèi)部接觸力分布等微觀層面的研究。

離散元法[20]（DEM）作為一種能應(yīng)用于工程領(lǐng)域的數(shù)值計算方法，將巖體視為不連續(xù)的離散塊體能有效地進行顆粒接觸間的微觀模擬。Pan 等[21]通過離散元模擬了粗粒土的熱誘導(dǎo)變形，發(fā)現(xiàn)粗粒土在加熱時發(fā)生膨脹，在冷卻時發(fā)生收縮，且收縮的幅度大于膨脹的幅度；同時還發(fā)現(xiàn)了由于顆粒間的滑動，顆粒重排的增加。Coulibaly 等[22]通過離散元法對粒狀材料進行了熱循環(huán)的模擬，發(fā)現(xiàn)了溫度循環(huán)對顆粒材料的結(jié)構(gòu)和性能都有顯著的影響。Zhao 等[23]則發(fā)現(xiàn)熱循環(huán)在粒狀材料內(nèi)部會引起顯著的應(yīng)力松弛。但是，目前離散元法進行的熱力耦合模擬大多采用的是顆粒體體積膨脹法，較少對真實傳熱過程的模擬。

所以，本文分別采用熱傳遞法和顆粒體積膨脹法，通過離散元軟件PFC2D來模擬單次溫度循環(huán)下的熱固結(jié)實驗。宏觀上分析試樣密實化和應(yīng)力衰減的規(guī)律，在微觀方面研究了顆粒間接觸的各項異性受溫度的影響規(guī)律，最后對比分析了兩種模擬方法在模擬過程中的同異性。

1 熱力耦合模型的機理

1.1 顆粒體積膨脹法

顆粒體積膨脹法是指不考慮顆粒間的熱能傳遞，顆粒溫度的變化直接通過參數(shù)的設(shè)置來實現(xiàn)。在模擬過程中，通過公式（1）實現(xiàn)溫度對顆粒體積的影響。

式中，ΔT 是顆粒溫度的變化量；α 是顆粒的熱膨脹系數(shù)；R0是顆粒變溫前的半徑；ΔR 是顆粒在變溫時產(chǎn)生的半徑的變化量。

顆粒體積的改變使顆粒間接觸的重疊量δn隨之發(fā)生變化，然后通過接觸力的計算公式（2）實現(xiàn)熱力耦合的模擬。

在采用顆粒體積膨脹法的溫度固結(jié)試驗?zāi)M過程中，全部顆粒的溫度都是相同的，且都是統(tǒng)一發(fā)生變化的。由于不需要考慮溫度的傳遞，采用顆粒體積膨脹法的熱力耦合在PFC2D軟件上的模擬更加快速。因此在進行如溫度循環(huán)試驗這類溫度需要長期持續(xù)變化的模擬中顆粒體積膨脹法的運用是必不可少的。

1.2 熱傳遞法

在采用熱傳遞法的熱力耦合模型中，每一個顆粒都代表一個蓄熱池。在兩個蓄熱池顆?；ハ嘟佑|后，中心會形成一根連接這兩個蓄熱池的傳熱管道簡稱熱管。通過熱管就可以實現(xiàn)熱能在兩個蓄熱池之間的轉(zhuǎn)移。

在模擬過程中，將邊界墻體作為恒定的熱源對內(nèi)部的顆粒進行熱能的供應(yīng)。且邊界墻體不隨溫度的變化而產(chǎn)生變形。在傳熱過程中，邊界墻體的熱能會從較近的顆粒向較遠的顆粒傳輸。而顆粒的溫度在邊界熱能的作用下會不斷升高直到與邊界墻體的溫度相同。與此同時，在顆粒溫度升高的過程中，顆粒體積則會根據(jù)公式（1）產(chǎn)生膨脹。并通過公式（2）實現(xiàn)溫度對試樣力學性質(zhì)的影響。

通過熱傳遞法進行的熱力耦合模擬更加接近現(xiàn)實中能源結(jié)構(gòu)對周邊土體的影響。但是由于其在模擬過程中需要兼顧熱傳遞和力學計算，會導(dǎo)致模擬速度的降低，不利于開展如多周期溫度循環(huán)固結(jié)試驗等有長時間傳熱試驗的模擬。

2 模型的建立及模擬方案

2.1 數(shù)值模型的建立

圖1 為試樣模型的顆粒粒徑分布曲線，曲線中的相關(guān)數(shù)據(jù)見表1。按照這一級配曲線，生成相對密實度Dr為0.2 的松散試樣其對應(yīng)的孔隙比為0.2376。

圖1 顆粒級配曲線

模型中力學參數(shù)和熱力學參數(shù)的選取主要參考石英砂的數(shù)值如表2 所示。

在顆粒生成后，循環(huán)計算至試樣內(nèi)部平均不平衡系數(shù)αf（見公式（3））小于1e-5。將平衡后的試樣進行側(cè)限壓縮固結(jié)。在固結(jié)過程中的邊界條件為在上下墻體施加100kPa 的恒定荷載，左右墻體保持靜止不動。在后續(xù)的溫度固結(jié)試驗過程中會保持這個邊界條件不變。側(cè)向固結(jié)后的試樣如圖2 所示，此時試樣的孔隙比e0為0.2312，將此狀態(tài)下的試樣視為溫度固結(jié)試驗中的初始狀態(tài)。

圖2 數(shù)值試樣模型

式中N 是顆粒的總數(shù)；C 是總的接觸數(shù)量；Ci是第i 個顆粒的接觸數(shù)量；F(ij)是第i 個顆粒的第j 個接觸的接觸力；Fb(i)是第i 個顆粒的體積力；Fc(i)是第i 個接觸的接觸力。

2.2 模擬方案

分別采用熱傳遞法和顆粒體積膨脹法進行溫度固結(jié)試驗的模擬。在模擬過程中，設(shè)置顆粒的初始溫度T0為20℃，溫度變化過程為20℃～60℃～20℃。

在熱傳遞法中，先將邊界墻體的溫度設(shè)置為60℃并保持溫度固定，在試樣的平均溫度升高到60℃后再將邊界墻體的溫度改為20℃，使試樣的平均溫度重新降低到20℃，以此實現(xiàn)溫度循環(huán)固結(jié)試驗的模擬。

在顆粒體積膨脹法中，以固定的溫度增量ΔT=1℃進行顆粒溫度的設(shè)置。試樣平均溫度具體的變化過程為20℃～21℃～22℃～......～60℃～59℃～58℃～......～20℃，如圖3所示。

圖3 溫度變化示意圖

最后，將兩種熱力耦合方法的模擬數(shù)據(jù)在宏觀角度（孔隙比、側(cè)向應(yīng)力）和微觀角度（各項異性）上進行對比分析。

3 試驗結(jié)果

3.1 孔隙比及側(cè)向應(yīng)力

為方便在數(shù)據(jù)圖中進行標注，將熱傳遞法簡稱為方法Ⅰ，將顆粒體積膨脹法簡稱為方法Ⅱ。圖4 為兩種方法下的孔隙比在溫度循環(huán)固結(jié)模擬下的對比結(jié)果。當溫度從20℃升高到60℃時，兩種方法分別從初始孔隙比0.23122降低到了0.23082 和0.23081，當溫度從60℃降低到20℃時，孔隙比又分別升高到了0.23114 和0.23117?？梢钥闯?，在經(jīng)過一次溫度循環(huán)后，試樣的孔隙比相較于初始狀態(tài)時是降低的。這是由于在恒定的邊界條件作用下，溫度的升高使顆粒的體積不斷膨脹降低了試樣的孔隙比，而溫度的降低使顆粒體積逐漸收縮，導(dǎo)致試樣孔隙比的增加，顆粒體積在不斷變化過程中產(chǎn)生了塑性變形，使試樣體積無法還原到初始狀態(tài)。這種孔隙比的變化規(guī)律與Liu 等[30]關(guān)于不同溫度下的孔隙率變化規(guī)律相同。同時也可以看出，兩種熱力耦合方法在孔隙比上的變化是高度吻合的。

圖4 孔隙比在溫度循環(huán)過程中的變化

由于在邊界條件中上下墻體施加了恒定的荷載，試樣在豎直方向上的應(yīng)力是恒定的，因此對于應(yīng)力的研究主要針對側(cè)向應(yīng)力。圖5 為試樣的側(cè)向應(yīng)力在熱固結(jié)試驗?zāi)M中的變化情況。隨著溫度的升高，側(cè)向應(yīng)力從初始應(yīng)力41.54kPa 分別增加到了48.34kPa 和48.11kPa。隨著溫度的降低，側(cè)向應(yīng)力也隨之降低到了40.44kPa 和40.59kPa。這也是由于側(cè)向墻體對試樣的約束力會隨著顆粒體積的膨脹而增加，隨著顆粒體積的收縮而降低的原因造成的。此外，可以看出模型Ⅰ與模型Ⅱ的側(cè)向應(yīng)力曲線十分相近的。

圖5 溫度循環(huán)過程中側(cè)圍壓力的變化

3.2 各項異性

參考Rothenburg 和Bathurst[31]關(guān)于各向異性的擬合公式（4）分別從法向接觸力和切向接觸力兩個方面對比方法Ⅰ和方法Ⅱ?qū)υ嚇拥挠绊憽?/p>

式中：fn和ft分別為法向接觸力和切向接觸力在θ（θ∈[0，2π]）角度上的大??；f0和fs為法向接觸力和切向接觸力的平均值；an和as分別表示了法向接觸力和切向接觸力的方向變化量，其數(shù)值越大，接觸力的各向異性就越明顯，由公式（5）計算所得；θn為最大平均法向接觸力的傾角；θt為平均切應(yīng)力為零時的傾角。

式中：E（θ）為接觸密度函數(shù)；ai表示分布的各向異性大??；θi表示i 對應(yīng)的各向異性方向，其中i 分別為n 和t。

圖6 為溫度固結(jié)試驗?zāi)M中平均法向接觸力隨溫度的變化曲線。當溫度從20℃升高到60℃時，方法Ⅰ和方法Ⅱ的f0從65.26kN 分別增加到了69.86kN 和70.06kN。當試樣的溫度重新降低到20℃時，f0又分別降低到了65.02kN 和65.05kN。相較于初始狀態(tài)分別降低了0.24kN和0.21kN。

為了進一步觀察法向接觸力在試樣內(nèi)部的分布情況，分別取平均溫度為60℃和完成溫度循環(huán)后（平均溫度為20℃）的試樣進行分析，如圖7 所示。當平均溫度分別為60℃和20℃時，法向接觸力在水平方向上的分布分別約為64.72kN 和52.51kN，在豎直方向上的分布則因為穩(wěn)定的邊界條件保持在76.49kN 左右。當試樣的平均溫度重新降低到20℃時，由于法向接觸力在水平方向上的分布降低，使其整體處于降低的趨勢，與圖6 中的變化趨勢相同。此外，可以看出，熱傳遞法和顆粒體積膨脹法在相同溫度下的法向接觸力的分布高度吻合。

圖6 平均法向接觸力隨溫度的變化

圖7 法向接觸力分布示意圖

圖8 為平均切向接觸力fs在熱固結(jié)試驗?zāi)M中的變化曲線。從圖中可以看出，隨著溫度的升高，兩種方法的平均切向接觸力fs從初始的1.735kN 分別降低到了0.637kN 和0.631kN。隨著溫度的降低，fs分別增加到了1.883kN 和1.855kN。相較于初始狀態(tài)，經(jīng)過溫度循環(huán)后的試樣的切向接觸力是有所增加的。

圖8 平均切向接觸力隨溫度的變化

如圖9 所示，分析了切向接觸力在平均溫度分別為20℃和60℃時的分布。當試樣的平均溫度從60℃降低到20℃時，切向接觸力在方向上分布的最大值從0.4635kN增加到了2.837kN，變化情況與圖8 中的曲線相吻合。同時，當平均溫度為60℃時，方法Ⅰ和方法Ⅱ的切向接觸力方向上的最大值分別為0.4635kN 和0.4237kN，其差異可以忽略不計。因此，兩種方法在切向接觸力的分布上也有高度的相似性。

圖9 平均切向接觸力分布示意圖

4 結(jié)論

為了進一步研究在能源樁的使用過程中，樁身溫度的變化對其周邊土體的影響。本文通過PFC2D軟件，利用離散元法（DEM），分別采用熱傳遞法和顆粒體積膨脹法對石英砂材料進行了熱固結(jié)試驗的模擬。分別從宏觀角度和微觀角度分析了結(jié)構(gòu)接觸面周圍土壤的熱耦合力學性能。并對比分析了兩種熱力耦合方法的模擬情況。得到以下結(jié)論：

（1）通過PFC2D軟件，利用離散元法，無論是基于熱傳遞法還是顆粒體積膨脹法，都能很好地實現(xiàn)溫度固結(jié)試驗的模擬，在模擬過程中，能夠明顯地體現(xiàn)出溫度對顆粒體積的影響。

（2）在宏觀角度上，在上下墻體施加穩(wěn)定100kPa 的荷載，左右墻體保持固定的邊界下，在溫度循環(huán)固結(jié)試驗?zāi)M中，顆粒體積產(chǎn)生塑性變形導(dǎo)致孔隙比的降低，試樣變化更加密實；因此也導(dǎo)致了試樣側(cè)向應(yīng)力的降低。

（3）在微觀角度上，由于在水平方向的分布隨側(cè)向應(yīng)力的降低而降低，法向接觸力在經(jīng)過溫度循環(huán)后是降低的。與之完全相反的是平均切向接觸力是增加的。

（4）熱傳遞法和顆粒體積膨脹法在循環(huán)的熱固結(jié)試驗中的模擬結(jié)果無論是在宏觀方面還是在微觀方面，都具有高度的一致性。

以上研究結(jié)果可為能源樁周土的研究提供參考依據(jù)。