沙粒內稟尺度對道床宏細觀力學特性影響研究

肖 宏，張智海，王宏閣，崔旭浩

(1.北京交通大學 土木建筑工程學院， 北京 100044；2.北京交通大學 軌道工程北京市重點實驗室， 北京 100044)

我國西北地區荒漠化嚴重，鐵路道床積沙現象普遍。散體道床在列車荷載作用下，沙粒會不斷侵入道床填充道砟空隙，改變原有的道床顆粒組分及級配。隨著道床中沙粒含量的增加，軌道結構易出現拱道、道床板結、鋼軌低接頭、三角坑等工務病害。這不僅使軌道幾何形位難以保持，還會加大線路養護維修工作量[1-2]。此外，細沙的灌入輕則降低道床服役性能，重則導致扣件銹蝕斷裂及鋼軌永久變形，直接影響行車安全[3]。因此，開展風沙區鐵路道床宏細觀力學特性研究非常必要。

沙粒侵入引發一系列線路服役狀態問題，受到諸多鐵道工程研究者的關注。由于物理試驗成本較高且試驗結果離散性大、影響因素難以控制[4]，有學者采用離散單元法(DEM)來研究沙粒與道砟之間的宏細觀力學特性。嚴穎等[5]建立了桶狀離散元試樣，對道砟和細沙顆粒進行了數值模擬，分析了不同含沙率下有效變形模量差異性，但在研究過程中將沙粒粒徑放大了20倍。高亮等[6]建立了三維道砟箱剪切模型，研究了黃沙對道床剪切性能的影響，將黃沙粒徑直接視為4～9.5 mm進行模擬，并未深究沙粒合理模擬尺度的選擇問題。以上研究表明，針對沙粒灌入，對道砟材料力學特性影響的數值試驗研究忽略了沙粒粒徑尺度效應的影響，未能準確模擬沙粒與道砟相互作用關系及顆粒之間接觸力傳遞規律，因而不能從本質上揭示道床內部多元混合顆粒之間的相互作用原理。可見，要探明含沙道床力學特性，首先要解決數值模擬中沙粒粒徑的合理選擇問題。

對于顆粒材料的尺度研究，目前主要關注宏觀、細觀、微觀3個層次[7]。不同粒徑尺度是相對的，不是固定不變的，應視研究問題而定。在細、微觀層次，顆粒材料可以被量化為不同粒徑的細小顆粒集合體，其力學性質隨顆粒粒徑量變逐漸引發顆粒集合體特性的質變。在研究沙粒與道砟之間的相互作用關系時，由于沙粒粒徑較小，全尺度道床模型中沙粒數量可能達到數千萬甚至上億，現有的計算條件根本無法精確模擬。為模擬細小顆粒的力學特性，許多學者利用顆粒縮放法來處理細小顆粒尺寸效應問題，并將其廣泛應用于工程研究中。任建莉等[8]利用相似理論和量綱分析法，對鑄鐵煤粉的放大顆粒進行了微觀參數標定，驗證了顆粒縮放理論的正確性。李永祥等[9]利用顆粒縮放原理，將小麥粉粒徑放大了6倍，進行了“Hertz-Mindlin with JKR”接觸模型參數標定。Sakai等[10]在流化床數值試驗中，將細顆粒體用粗顆粒進行替換，結果表明，粗顆粒模型可以較為準確的模擬原始顆粒材料的力學特性。Weinhart等[11]建立了粗顆粒筒倉離散元仿真模型，提出了在高度局部化系統中空間粗顆粒寬度選擇的重要性。但以上研究都是利用顆粒縮放法來表征單一顆粒力學特性，而風沙道床屬于多元混合顆粒體系，內部顆粒之間傳力方式復雜多樣，其粒徑縮放原理是否仍然適用還需進一步驗證。尤其在沙粒粒徑縮放對道床力學特性的影響研究方面，目前少有涉及。

為描述風沙道床數值模擬過程中不同沙粒尺寸與道砟粒徑的關系，將數值模型中沙粒最小顆粒粒徑與道砟最大顆粒粒徑的比值，定義為內稟尺度。該指標重在反映可變沙粒粒徑與不變道砟顆粒粒徑的關系，是風沙道床數值模型內部顆粒跨尺度作用的一種體現。

基于上述分析，本文將從不同沙粒粒徑的填充方法入手，基于顆粒縮放原理提出一種新的研究多元混合顆粒粒徑尺度效應的方法——“胞分法”。該方法是基于細胞分裂分化的思想，通過閾值人為控制胞體的分裂分化數量與次數，嚴格保證局部區域內新產生的胞體與原胞體所占空間大小相同。本文借助胞分原理來保證了不同沙粒尺寸下道砟空隙中沙粒質量守恒，進而研究沙粒粒徑尺寸對道床宏細觀力學特性的影響，揭示多元混合顆粒之間力鏈傳遞及演變規律，給出數值試驗沙粒的合理尺寸，以期為風沙區鐵路道床力學特性研究提供先決條件。

1 胞分法分析模型

1.1 胞分法理論模型

風沙道床是多元混合材料相互作用構成的復雜顆粒體系，其力學特性取決于顆粒體系的細、微觀作用力。不同沙粒尺寸對道床結構產生不同的力學效應，而沙粒比表面積和質量隨粒徑的改變而逐漸變化，會影響道砟與沙粒之間接觸的產生[12]。顆粒縮放法[8-10]雖能模擬細小顆粒，但僅用于單一顆粒力學特性的研究，而針對多元混合顆粒研究未必適用，且風沙道床顆粒體系復雜多變，需保證道砟空隙含沙量不變，減少計算誤差，確保計算結果真實可靠。

風沙道床顆粒數量較多，利用PFC3D模擬，存在小顆粒懸浮、模型平衡困難、計算效率低等問題，且在分析多根軌枕長度的線路時，存在局限性。已有文獻[13-15]表明，二維模型可以反映道砟顆粒之間的接觸特性，因此本文擬采用PFC2D模型進行分析。

為考慮沙粒粒徑對道床受力狀態的影響，保證道砟空隙含沙量不變，道床體系沙粒質量守恒，本文將單個沙粒看成一個胞體，新產生的沙粒視為胞元，引入閾值含沙量γ來建立原胞體與新胞元之間的聯系，即

(1)

(2)

(3)

(4)

(5)

式中：ms為道床中沙粒總質量；mb為道砟顆粒總質量；Ns為沙粒數量；Nb為道砟顆粒數量；ρs為沙粒密度；Vis為第i個沙粒的體積；ρb為道砟密度；Vjb為第j個道砟的體積；Sis為第i個沙粒的面積；Sjb為第j個道砟的面積；Ts為沙粒模擬單元的厚度；Tb為道砟模擬單元的厚度。

若沙粒采用單位圓盤模擬，垂直于平面方向道砟顆粒厚度Tb與沙粒厚度Ts均為單位厚度。由式(4)和式(5)可得，單個顆粒體積與圓盤面積數值相等，含沙量γ與沙粒數量、顆粒所占面積正相關，且第i個沙粒面積Sis與沙粒粒徑d2成正比。因此要保證道砟空隙中含沙量不變，需確定沙粒數量與沙粒粒徑之間的關系。

本文擬通過控制分裂次數與新生胞元的半徑來逐漸逼近沙粒的真實細、微觀作用效果。在模擬時，首先采用先張后縮法將沙粒粒徑放大10倍并視為原胞體，也可稱為一分胞體，然后利用fish語言來保證每次分裂產生的新胞元粒徑相同，且胞元之間為最密排列方式，減小由于顆粒錯動、重新排列造成的計算誤差，縮減離散元模型的平衡時間。

在利用胞分法建立不同沙粒粒徑離散元模型時，將沙粒粒徑縮放區間離散化，文獻[5]擬定了沙粒粒徑縮放上限值為10。考慮計算效率和模擬精度等問題，將分裂次數n的上限值取為8，使第8次分裂后的沙粒粒徑處于實際粗沙粒徑1 mm附近。分裂次數和沙粒胞元半徑的關系為

(6)

式中：R1為一分胞體的沙粒半徑；Rn為第n次分裂后沙粒半徑。

沙粒胞體分裂過程見圖1，圖1中黑色虛線表示原胞體輪廓，綠色虛線是各個胞元質心的連線。通過局部坐標反演的方法，利用質心線構成正三角形、正方形等幾何體，推算出以原胞體質心為中心的周圍各個胞元質心的位置坐標。借助fish語言將每次分裂的沙粒質心坐標導入PFC中進行顆粒質心更新計算，實現沙粒的分裂分化過程，使沙粒粒徑逐漸減小，向真實沙粒粒徑逐漸逼近，表1是道床內部沙粒粒徑與分裂次數的對應關系，其中，d為中細沙粒的真實顆粒粒徑，一般在0.25～0.5 mm之間取值。

表1 沙粒胞元顆粒信息統計

圖1 沙粒胞體分裂過程

1.2 顆粒流模型的構建

1.2.1 軌道結構離散元模型

在建立風沙區軌道結構離散元模型之前，需建立精細化的道砟顆粒模型。已有研究表明[16]，道砟顆粒外形、級配、棱角的各向異性等是影響道砟精確模擬的關鍵。為考慮道砟顆粒棱角及外在形態特征[17-18]，本文利用三維逆向工程軟件提取了道砟三維廓形[19]，沿著最長軸線進行投影得到道砟顆粒二維廓形文件，通過自編fish語言導入PFC中生成不規則道砟顆粒，并按照既有線一級碎石道砟顆粒級配生成道砟顆粒，粒徑范圍為25～63 mm。

圖2 風沙區軌道結構離散元模型

結合現場實際軌道鋪設情況，道床厚度取0.35 m，道砟用clump單元模擬，沙粒采用ball單元模擬。一般而言，車輪荷載影響區域范圍[20]為5根軌枕，為消除邊界效應的影響，特利用塊體疊加拼裝法生成7根軌枕長度的風沙道床仿真模型。軌枕采用ⅩⅡ型混凝土枕，以顆粒簇塊體方式導入道床模型中；扣件采用兩個半徑為35 mm的ball單元模擬；綜合考慮60 kg/m鋼軌的實際界面尺寸，離散元模型中采用直徑為0.15 m的28個小球規則排列黏結而成，利用等質量法[13]計算單個顆粒密度約為510 kg/m3。考慮鋼軌與列車之間的相互作用，根據一個轉向架的距離來設置兩個加載顆粒的間距，利用fish語言將荷載時程曲線分別施加在3號和4號軌枕中間及7號軌枕左側的加載顆粒上，來模擬車輪與鋼軌之間的荷載傳遞。采用“胞分法”建立不同沙粒粒徑的多尺度軌道結構離散仿真模型。

1.2.2 本構模型

表2 離散元模型細觀參數

2 現場試驗及模型驗證

2.1 現場試驗

為研究沙粒作用對道床力學特性的影響，在風沙嚴重的甘萬鐵路昌吉高勒站和巴音杭蓋站區段開展了現場輪軌力與道床加速度動測試驗，分別布置鋼軌垂向力測點和道床加速度測點。測得鋼軌垂向力最大為117.69 kN，道床最大加速度為1.41g。

2.2 列車荷載施加

為研究實際列車荷載作用下風沙道床沙粒多尺度效應，以現場速度58 km/h試驗列車C70E的實測輪軌垂向力為依據，利用程序語言對實測曲線相鄰數據點進行了加密處理，生成了荷載數據文件，采用table命令將荷載數據文件導入到PFC中，并將其施加到圖2(b)中的車輪上。現場實測輪軌垂向力時程曲線見圖3。為準確捕捉圖3中的荷載曲線峰值數據，在數值加載過程中將計算時步調至10-7～10-8之間。考慮列車循環荷載特性及離散元模型的計算效率問題，本文選取了兩節車廂8軸荷載模式的時程曲線來模擬實際列車荷載，見圖4。

圖3 現場實測輪軌垂向力時程曲線

圖4 試驗現場C70E敞車示意(單位：mm)

2.3 模型可靠性驗證

為驗證模型的正確性，將圖4中現場實測的輪軌垂向力荷載譜施加到本文所建的離散元模型中，計算得到風沙區道床振動加速度最大峰值為1.44g，現場實測道床振動加速度最大峰值為1.41g，試驗值與離散元模擬值僅相差2.13%，數值吻合度高。現場實測道床加速度時程曲線見圖5，離散元仿真道床加速度時程曲線見圖6。對比圖5和圖6可知，在列車荷載作用下二者的曲線走勢及增長幅值變化規律幾乎相同，實測加速度值因受現場環境、輪軌噪聲信號、測試信號線固定不到位、線路不平順等因素的干擾，基線較粗，與數值分析結果細部形態略有差異，但不影響道床加速度幅值隨時間的變化趨勢。綜上所述，本文從數值大小和圖形規律都驗證了模型的可靠性，也說明了離散元模型中微觀接觸參數可以用于沙粒尺寸的研究。

圖5 現場實測道床加速度時程曲線

圖6 離散元仿真道床加速度時程曲線

3 道床受力特性分析

為揭示細沙顆粒大小對道床力學行為的影響，以及純粹的顆粒粒徑尺度縮放是否會引起道床受力體系的改變，論文利用胞分法建立了8種不同沙粒粒徑的離散元仿真模型，從宏細觀角度分析了不同沙粒粒徑作用下道床的力學特性，數值試驗方案見表3。

表3 不同沙粒粒徑計算方案

3.1 道床細觀接觸力演變規律分析

3.1.1 平均接觸力分析

道砟顆粒平均接觸力是衡量道砟接觸狀態的重要指標，可以反映道床內部道砟顆粒的平均受力情況。為研究分析方便，將8軸荷載的各軸按照從左到右依次編號，列車軸次與道砟平均接觸力的關系曲線見圖7。

圖7 道砟平均接觸力變化

由圖7可知，當沙粒處于一分胞體狀態時，道砟平均接觸力最大為139.557 N；當沙粒處于八分胞體狀態時，道砟平均接觸力最大為83.306 N，相對減小了40.31%。隨著沙粒分裂次數的增加，沙粒內稟尺度逐漸減小，平均接觸力相對減小百分比逐漸趨于0，道砟平均接觸力逐漸收斂趨于定值。六分胞體時，道砟平均接觸力變化幅度較小，因而可以反映實際沙粒的作用效果。這是由于在含沙量不變的情況下，隨著沙粒粒徑的減少，道床中沙粒數量成倍增加，使道床內部顆粒接觸增多，力鏈的傳遞和延伸方向逐漸增多，道砟顆粒之間的接觸力逐漸減小，平均接觸力也逐漸減小。

3.1.2 接觸力概率密度分析

為從細觀角度揭示不同沙粒粒徑作用下道床內部顆粒接觸力分布規律，采用張科芬等[4]、Liu等[27]提出的度量顆粒間力大小的概率密度曲線，衡量在一定范圍內顆粒接觸力出現的概率，借此直接觀察顆粒體系中強接觸力與弱接觸力分布規律。對于風沙道床散粒體系，內部顆粒接觸力大小是離散且有限的，因此可利用fish語言遍歷顆粒之間的接觸力，接觸力的概率密度為

(7)

式中：x為接觸力的離散值，介于[0，b]之間，b值取決于接觸力區間上限值；v(x)為顆粒體系接觸力小于x的個數；N為顆粒體系總的接觸對；F(x)為接觸力分布函數；p(x)為接觸力概率密度函數。接觸力概率密度統計結果見圖8。

由圖8可知，隨著顆粒間接觸力的逐漸增大，接觸法向力與接觸切向力概率密度逐漸減小，且小的接觸力占比更高。道床結構不受外力時，隨沙粒胞分次數的增加，沙粒粒徑逐漸減小，無論是法向接觸力還是切向接觸力，概率密度分布曲線都越來越陡，最大接觸力b值越來越小，可以明顯看出三分胞體后，接觸力概率密度曲線變化趨勢逐漸趨于穩定，出現大接觸力概率愈小，小接觸力概率愈大；峰值荷載作用下最大接觸力概率密度隨沙粒粒徑減小而逐漸減小，道床接觸法向力與切向力大小趨于均勻化。對比可知，法向接觸力在第五次分裂后概率密度曲線趨于穩定，而切向接觸力在第六次分裂后概率密度才趨于穩定，說明沙粒粒徑對切向接觸力分布影響更大。峰值荷載作用下道床內部法向和切向接觸力離散值是自重作用下的5倍，表明列車荷載作用是道床內部產生強接觸力的主要原因，弱接觸力是道床顆粒體系接觸力的主要表現形式。

綜上可知，沙粒內稟尺度會影響道床顆粒間的接觸力分布規律，且對峰值荷載作用下接觸力分布影響更為顯著。為減小沙粒粒徑對數值試驗的影響，可選擇六分胞體沙粒粒徑進行模擬計算，以保證計算結果的可靠性，揭示沙粒對道床受力影響的內在機理。

3.1.3 接觸力方向演變規律分析

為更加直觀地分析沙粒內稟尺度對道床內部顆粒體系接觸力各向異性演變規律的影響，分組統計了峰值荷載作用下各方向接觸力的數目及各接觸力絕對值累計值(接觸力合力)，利用自編CAD小程序繪制了道砟顆粒法向接觸合力和切向接觸合力各向異性玫瑰圖，見圖9和圖10。圖9和10中同心圓半徑大小代表不同等級接觸合力量值，圓盤外圍角度表示各組接觸合力方向。由于繪圖前，設置了繪圖基準力，圖9和10中所顯示的是各組接觸合力量值，該量值是各組接觸力絕對值累計值除以基準力得到，因此無量綱。圖9的繪圖基準力為100 kN，圖10的繪圖基準力為50 kN。

各組法向接觸合力分布；平均法向接觸合力。

各組切向接觸合力分布；平均切向接觸合力。

3.2 道床應力場細觀分析

為了量化不同沙粒粒徑作用下，道床應力場的變化規律，選取了列車荷載作用影響最大區間，通過自編fish語言提取單個顆粒的平均接觸應力，按照一定比例繪制了峰值荷載作用下不同沙粒粒徑的道床應力場分布云圖，見圖11。

圖11 不同沙粒尺寸下道床應力場分布

由圖11可知，列車荷載作用下風沙道床內部應力分布較為均勻，低應力(綠色區域)分布面積較大，高應力(紅色區域)僅分布在道床局部區域，且軌枕底部應力較大，而枕側區域平均應力較小，這與軌枕傳力方式、枕底顆粒接觸數目密切相關。通過觀察發現，隨著沙粒內稟尺度的逐漸減小，沙粒數量倍增，道床中高應力區域逐漸減小，低應力區域占比越來越高，道床應力均勻化程度顯著增強。一分胞體道床內部高應力區域占比較高，但隨著沙粒的分裂次數的增加，高應力區域明顯減少，直到六分胞體時，顆粒之間的強接觸力減少，道床應力衰減幅度明顯降低，道床內部高應力區域占比趨于穩定，處于低應力狀態。這是由于列車荷載與顆粒接觸數目是影響道床細觀應力狀態的主要因素，在特定的列車荷載作用下，顆粒接觸數目受顆粒數量的影響較大。一般而言，顆粒數目越多，細觀接觸面越大，顆粒之間接觸配對相對容易，顆粒之間的接觸對愈多，作用力相對分散，強接觸力不易產生，道床內部趨于低應力狀態。

為進一步探究沙粒內稟尺度與道床內部應力的作用關系，統計了峰值荷載作用下，不同沙粒粒徑的道床最大應力和平均應力，見圖12。由圖12可知，不同沙粒粒徑的道床內部應力差異性較大，且隨著沙粒粒徑的減小，即沙粒胞元個數的增加，道床最大應力和平均應力也逐漸減小；自五分胞體后，道床最大應力在375.124～380.225 kPa之間，平均應力在32.206～32.501 kPa之間，應力幅值變化不大，沙粒粒徑對道床應力的作用效果減弱，這與圖11道床應力分布云圖得出的結論一致。

圖12 不同沙粒粒徑作用下道床應力值

道砟顆粒之間的接觸力大小和接觸數目是影響道床內部平均應力的關鍵，而沙粒的侵入改變了原有道床的受力體系，顆粒之間的無規則排列，導致接觸力分布方向的各向異性，且隨著沙粒粒徑的減小，沙粒數量愈多，道床內部接觸數目愈多，作用力的傳遞及力鏈延伸方向多樣，道床內部的強接觸力愈少，顆粒之間的最大接觸力也逐漸變小。因此，道床的最大應力和平均應力隨沙粒粒徑的減小而逐漸減小，但當沙粒粒徑減小到一定程度時，沙粒內稟尺度對道床作用減弱，出現了道床應力變化穩定點，這為沙粒侵入道床數值試驗最小顆粒尺度的選擇，提供了理論依據。

4 結論

為研究沙粒內稟尺度對道床宏細觀力學特性的影響，確定風沙區道床數值試驗沙粒粒徑的合理尺度，考慮道砟空隙中沙粒質量守恒和無相互重疊原則，首次采用“胞分法”建立了不同沙粒粒徑尺度的精細化軌道結構離散元模型，分析發現沙粒粒徑對道床宏細觀力學特性有明確影響。

(1)沙粒粒徑會影響道砟的平均接觸力大小。沙粒處于一分胞體狀態時，道砟顆粒最大平均接觸力為139.557 N；沙粒處于八分胞體時，道砟最大平均接觸力為83.306 N，相對減小了40.31%；隨著沙粒胞分次數的增加，沙粒粒徑逐漸減小，道砟平均接觸力也逐漸減小。

(2)沙粒內稟尺度會影響道床顆粒間的接觸力分布規律，對峰值荷載作用下接觸力分布影響更為顯著。隨著顆粒間接觸力的逐漸增大，接觸法向力與接觸切向力概率密度逐漸減小，小接觸力占比更高。為減小沙粒粒徑對數值試驗的影響，建議選擇六分胞體沙粒粒徑進行道床受力分析，以保證研究結果的可靠性。

(3)沙粒內稟尺度主要影響切向接觸合力量值分布方向各向異性，而對呈“蝴蝶狀”分布的法向接觸合力量值分布方向影響較小。隨著沙粒內稟尺度愈減，即逐漸逼近真實沙粒粒徑時，法向接觸合力量值各向異性增強，而切向接觸合力量值在不同沙粒粒徑作用下接觸方向分布形態各異，接觸力量值主分布方向不斷變化，各向異性更為顯著。另外，接觸合力方向玫瑰圖也再現了隨著沙粒粒徑增大，道砟顆粒接觸力由弱到強的演變過程。

(4)沙粒粒徑與道床內部應力分布存在關聯。隨著沙粒粒徑的減小，道床最大應力和平均應力都逐漸減小。當沙粒粒徑減小到一定程度時，沙粒內稟尺度對道床作用減弱，出現了道床應力變化穩定變化點。

(5)通過對不同沙粒尺寸下道床的宏細觀力學特性進行分析，采用球形顆粒模擬沙粒時，沙粒半徑建議取為0.51～1.02 mm，即沙粒內稟尺度比為0.016 2，可以較為真實的模擬沙粒與道砟的實際作用效果。