基于CFD-DEM的有砟道床多孔介質(zhì)模型研究

殷 浩，高 亮

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

有砟道床由散粒體的道砟顆粒堆積密實而成，隨著行車速度的提高，會出現(xiàn)道砟顆粒飛離道床，擊打車輛、軌道及線路周邊設施的情況，即道砟飛濺現(xiàn)象。近年來，高速鐵路道砟飛濺現(xiàn)象持續(xù)增多[1-4]，其對軌道結構和列車具有巨大危害，極大增加了線路和車輛的養(yǎng)護維修成本，威脅行車安全。

道砟飛濺受軌道結構流場特征影響顯著，為分析其發(fā)生和發(fā)展機理，國內(nèi)外學者開展了大量研究工作。基于數(shù)值方法的有砟軌道流場特征分析，在飛砟機理研究中起到了不可或缺的作用，通過計算還原軌道流場，能夠準確分析道砟飛濺易發(fā)生區(qū)域及其發(fā)展規(guī)律。Sima等[5]較早采用計算流體動力學(CFD)方法建立了車廂底部與軌道結構的數(shù)值模型，對基于軌道流場特征的道砟飛濺機理進行仿真研究。Rocchi等[6]采用高速鐵路線路測試數(shù)據(jù)對所建立的有砟軌道CFD數(shù)值模型進行驗證，并分析了軌道流場與道砟所受荷載之間的關系。殷浩等[7]建立了軌道與車廂結構的數(shù)值模型，并精確考慮二者外形因素影響，分析了行車速度和道床尺寸等參數(shù)對道砟飛濺的影響。

有砟道床的散粒體特征對軌道流場影響顯著，但既有仿真研究多采用連續(xù)壁面單元(wall單元)對散體道床進行模擬，無法反映道砟級配和形狀等顆粒特征對軌道流場的影響，不能準確分析道砟飛濺的發(fā)生和發(fā)展機理。García等[8-9]通過設置道床壁面單元的等效粗糙度參數(shù)，模擬道床表面的不連續(xù)特征。然而，該方法仍是對道床表面顆粒特征的過度簡化，且無法反映散體道床對流經(jīng)其內(nèi)部流體的影響。

郄錄朝等[10]在開展道砟飛濺CFD仿真分析時，將有砟道床考慮為多孔介質(zhì)模型，通過現(xiàn)場試驗數(shù)據(jù)對模型參數(shù)進行標定。多孔介質(zhì)模型能夠反映結構內(nèi)部的孔隙特征，真實還原散粒體道砟對流經(jīng)道床流體的阻力影響，該方法有效提高了流場的還原精度。然而，道床多孔介質(zhì)模型參數(shù)在不同線路條件下的差異顯著，采用現(xiàn)場試驗的參數(shù)獲取方法成本較高，且無法得到既有線路之外的道床模型參數(shù)。

針對上述問題，本文提出有砟道床的計算流體動力學-離散元(CFD-DEM)耦合分析方法，采用基于CFD-DEM耦合的精細計算方法，對局部道床多孔介質(zhì)模型中所涉及的關鍵參數(shù)進行獲取和驗證。該參數(shù)獲取方法既能反映道砟級配和形狀等顆粒特征，還能體現(xiàn)道床與流體的相互作用，實現(xiàn)了軌道流場的真實還原。為精確分析基于大尺度的有砟軌道流場特征的飛砟發(fā)生機理，以及基于道砟飛濺運動特征的飛砟發(fā)展機理提供了必要的技術手段。

1 有砟道床CFD-DEM耦合分析方法

采用DEM分析模塊對有砟道床進行模擬，能夠反映道砟形狀、級配、位置等顆粒特征；采用CFD分析模塊計算流經(jīng)道床的軌道流場，能夠反映散體道床與流體的相互作用關系。通過編制CFD-DEM耦合算法自定義函數(shù)，描述有砟道床與流體之間的動量和質(zhì)量數(shù)據(jù)傳遞，真實還原有砟軌道流場特征，對道床多孔介質(zhì)模型參數(shù)進行獲取和驗證。

1.1 有砟道床DEM分析模塊

采用DEM模擬有砟道床結構，將道砟顆?？紤]為相互獨立的單元個體。采用Hertz-Mindlin接觸模型計算顆粒之間接觸力，根據(jù)牛頓第二定律建立運動方程，計算顆粒在迭代時步內(nèi)的位移，進而對道砟顆粒的位置信息進行更新。道砟顆粒單元的運動方程為

(1)

式中：m為顆粒質(zhì)量；u為顆粒位移；Fi為道砟顆粒與流體之間的相互作用力；F為顆粒間作用力。

不同形態(tài)特征的道砟顆粒在流場中的受力和運動特征差異顯著，本文采用三維激光掃描法，對道砟顆粒的真實外形進行獲取[11-12]，并將其作為模板導入DEM計算模塊。該道砟輪廓模板僅以虛擬面的形式存在，無法考慮道砟顆粒之間的相互作用，故采用單元球模型對其進行填充。根據(jù)道砟輪廓模板的形狀和尺寸，確定填充球單元的位置和粒徑。不同形狀特征的道砟輪廓模板，以及填充后的道砟單元模型見圖1。

圖1 道砟顆粒DEM模型示意圖

在計算空間內(nèi)隨機生成不同形態(tài)狀特征的道砟顆粒仿真單元，構成有砟道床離散元數(shù)值模型。根據(jù)TB 10621—2014《高速鐵路設計規(guī)范》[13]，設置道床頂面寬度3.6 m，道床厚度0.35 m，邊坡坡度1∶1.75。模型中的結構單元體均滿足剛性假設，參考文獻[14-15]開展的有砟道床離散元仿真分析，并結合課題組既有研究，確定道砟模型細觀參數(shù)見表1。

表1 道砟DEM模型參數(shù)

在道床模型建立過程中，需分層生成道砟顆粒，并采用墻體單元以較慢的速度對其逐層壓實，確保每層道砟均達密實穩(wěn)定狀態(tài)。建立的有砟道床DEM分析模塊見圖2。

圖2 有砟道床DEM分析模塊

1.2 軌道流場CFD分析模塊

流體的流動需滿足質(zhì)量守恒定律，即單位時間流體微元體增大的質(zhì)量，與流入微元體的凈質(zhì)量相等。在此基礎上，通過在守恒方程加入額外的顆粒單元體積分數(shù)項ε，考慮道砟顆粒對流體的影響，流體質(zhì)量守恒方程為

(2)

式中：u、v和w是速度矢量分別在x、y和z方向上的分量；ρ為流體密度；t為時間。

流體流動還需滿足動量守恒定律，即微元體中流體動量對時間的變化率與微元體所受合外力相等，動量守恒方程也需考慮顆粒單元體積分數(shù)項ε的影響，表達式為

(3)

(4)

(5)

式中：μ為流體動力黏度；P為流體微元體壓力；Su、Sv和Sw為動量守恒方程廣義源項。

采用多面體單元對軌道流場CFD分析模塊的計算域進行網(wǎng)格劃分，該模塊中的軌道結構尺寸與DEM模塊一致，所建立的網(wǎng)格模型見圖3。在CFD-DEM耦合算法中，CFD模塊與DEM模塊中數(shù)據(jù)的交互傳遞是以全局坐標為基礎開展的，故需確保二者具有相同的全局坐標，進而實現(xiàn)道砟顆粒與其所處區(qū)域軌道流場的相互關聯(lián)。

圖3 軌道流場CFD分析模塊

1.3 CFD-DEM耦合算法

(6)

式中：ρf為流體密度；Ap為顆粒投影面積；Vf和Vp分別為流體和顆粒的體積；Cd、κ分別為阻力系數(shù)及冪指系數(shù)，計算式分別為

(7)

(8)

其中，Re為雷諾數(shù)。

(9)

(10)

(11)

式中：Cb為附加體力系數(shù)，取0.5。

DEM模塊與CFD模塊通過網(wǎng)格單元進行數(shù)據(jù)交換，顆粒單元體積分數(shù)項ε的獲取也在網(wǎng)格單位開展。計算耦合模塊的顆粒體積分數(shù)項時，首先對顆粒單元的質(zhì)心坐標進行檢索，若顆粒質(zhì)心位于目標網(wǎng)格單元區(qū)域內(nèi)，則認為該顆粒屬于目標網(wǎng)格單元，顆粒體積也被計入目標網(wǎng)格單元的體積分數(shù)項的計算中[18]。網(wǎng)格單元顆粒體積分數(shù)項ε可計算為

(12)

式中：Vmesh為網(wǎng)格單元體積；Vi為質(zhì)心位于所計算流體網(wǎng)格單元內(nèi)部的顆粒單元體積；n為該網(wǎng)格單元中顆粒數(shù)量。

在每個計算時步內(nèi)，CFD模塊首先檢索計算域內(nèi)流場的初始參數(shù)，以及顆粒單元的位置信息。在此基礎上，對流場數(shù)據(jù)進行求解計算，收斂后進一步計算流場對處于其內(nèi)部的每枚顆粒單元受到的作用力Fi。DEM模塊通過顆粒運動狀態(tài)以及顆粒與其他顆粒單元或邊界單元之間的接觸狀態(tài)，計算DEM模塊中顆粒所受到的接觸力。并進一步考慮CFD模塊對顆粒產(chǎn)生的作用力，由牛頓運動定律，計算顆粒的位置和速度等運動信息。在下一個時步內(nèi)，將DEM模塊中更新顆粒運動信息后的顆粒模型數(shù)據(jù)導入CFD計算模塊中，再次重新檢索計算域內(nèi)流場參數(shù)和顆粒位置信息，對流場數(shù)據(jù)和顆粒受力進行求解，依次開展迭代計算，直至滿足計算條件要求。在計算域邊界位置，采用面網(wǎng)格單元對DEM中的計算邊界進行描述，進而實現(xiàn)與CFD網(wǎng)格邊界單元的耦合對應。

基于CFD-DEM耦合分析方法的有砟軌道模型計算精度，由DEM模塊計算精度、CFD模塊計算精度，以及二者之間的耦合算法共同控制。其中，DEM模塊計算精度會受到分析時間步長、接觸模型、道砟顆粒球簇單元填充數(shù)量和質(zhì)量等因素的影響；CFD模塊計算精度會受到分析時間步長、網(wǎng)格類型和尺寸、控制模型等因素的影響；CFD-DEM耦合算法主要通過控制兩個模塊之間的數(shù)據(jù)交換頻率影響分析計算精度。CFD-DEM耦合算法流程見圖4。

圖4 CFD-DEM耦合算法流程

2 有砟道床多孔介質(zhì)模型

在采用CFD-DEM耦合算法開展軌道流場分析時，在每個迭代時步內(nèi)均需對大量道砟顆粒單元的位置和受力信息進行遍歷檢索。若采用該方法開展局部范圍的軌道流場分析，其計算精度和效率均較高；而在涉及大尺度的高速列車-有砟軌道的流場分析時，計算效率會受到DEM模塊算法的限制而顯著不足。故本文將有砟道床考慮為能反映道床對流體阻力影響的多孔介質(zhì)模型，采用CFD-DEM耦合算法對道床多孔介質(zhì)模型參數(shù)進行獲取和驗證。為研究基于大尺度高速列車-有砟軌道道砟飛濺流場特征和運動特征的飛砟機理提供技術手段。

2.1 多孔介質(zhì)數(shù)學模型

有砟道床是由道砟顆粒和空氣組成的多相物質(zhì)；道砟顆粒間的孔隙相對較小，道砟占據(jù)道床結構的絕大部分區(qū)域；顆粒間孔隙能夠相互連通，使空氣能夠在其內(nèi)部和表面流動。有砟道床的上述結構特征符合Bear[19]、Warrick等[20]對多孔介質(zhì)結構的定義，故在分析道床流場特征時，可將其考慮為多孔介質(zhì)結構。

有砟道床多孔介質(zhì)模型能夠反映道床對空氣流場的阻力作用，在數(shù)學模型上表現(xiàn)為流體動量方程的附加動量源項Si，該附加動量源項由黏性損失項和慣性損失項兩部分構成，各向同性多孔介質(zhì)的動量源項Si可表示為

(13)

附加動量源項Si與流經(jīng)多孔介質(zhì)模型的壓力梯度直接相關，故結合式(13)可知，道床上下游的流體壓降與流經(jīng)道床的流體速度之間的關系為

(14)

式中：Δn為沿流體流動方向的道床縱向長度。

2.2 基于CFD-DEM的多孔介質(zhì)參數(shù)獲取

基于有砟道床 CFD-DEM 耦合算法，計算提取流經(jīng)道床的流體壓降，并將其與流場風速擬合，得到壓降與風速擬合關系式，將該關系式與多孔介質(zhì)模型的動量方程附加源項對比分析，得到道床模型的黏性阻力系數(shù)和慣性阻力系數(shù)，實現(xiàn)道床多孔介質(zhì)模型參數(shù)的獲取。

將沿道床模型縱向的邊界條件分別設置為風速入口邊界和壓力出口邊界，其他表面均為無滑移壁面條件，用于限制道床內(nèi)道砟顆粒單元的流出，見圖5。當流場達穩(wěn)定后，分別對道床入口和出口邊界面處的風壓數(shù)據(jù)進行積分，提取流場上游和下游平均風壓，計算得到流經(jīng)道床的空氣流場壓降。當入口邊界的輸入風速為10 m/s時，計算過程中的流場壓降曲線見圖6。由圖6可知流場壓降先逐漸增大，后趨于穩(wěn)定，提取該穩(wěn)定值作為耦合模型的壓降數(shù)據(jù)。

圖5 耦合模型邊界條件

圖6 流體壓降曲線(流速10 m/s)

設置入口風速分別為 5、10、15、20、25、30 m/s，計算模型上下游流場壓降，提取不同輸入風速條件下的壓降數(shù)據(jù)。將輸入風速與壓降數(shù)據(jù)進行二次擬合，擬合曲線見圖7，壓降ΔP與風速v之間的擬合關系式為

圖7 流體壓降與輸入風速關系

ΔP=68.307v2+24.735v

(15)

擬合曲線相關系數(shù)R2=0.999 7，可認為二者完全相關。

將式(15)與式(14)比對分析，可得到道床多孔介質(zhì)模型的黏性阻力系數(shù)與慣性阻力系數(shù)。Δn取沿流體流動方向的道床縱向長度1.8 m，考慮外界條件為標準大氣壓強101.325 kPa，溫度20 ℃，此時空氣密度ρ=1.205 kg/m3，空氣黏度μ=1.79×10-5Pa·s。通過計算，得到道床多孔介質(zhì)模型的黏性阻力系數(shù)1/α=767 678，慣性阻力系數(shù)C′=62.984 9。

2.3 多孔介質(zhì)模型參數(shù)驗證

將由CFD-DEM耦合算法計算得到的黏性阻力系數(shù)1/α與慣性阻力系數(shù)C′施加于有砟道床多孔介質(zhì)CFD模型中，將其設置為與耦合算法相同的邊界條件，開展流場仿真計算。將采用道床多孔介質(zhì)CFD模型與CFD-DEM耦合模型計算得到的有砟軌道風壓分布和壓降數(shù)據(jù)對比分析，對該參數(shù)獲取方法進行驗證。

將CFD-DEM耦合模型與道床多孔介質(zhì)模型的計算域入口風速均設置為10 m/s。耦合模型中流體流經(jīng)道床時的道砟顆粒壓強云圖見圖8。由圖8可知，顆粒模型統(tǒng)一顯示為圓球單元，未體現(xiàn)其外形和粒徑特征。

圖8 CFD-DEM耦合模型道砟顆粒壓強云圖

由圖8可知，不同區(qū)域道砟受到的壓強差異顯著。流場入口附近顆粒多受正壓影響，出口附近顆粒多受負壓影響，且正壓幅值顯著高于負壓幅值。道床上游與下游之間道砟所受壓強均勻變化，該壓強變化即為流體流經(jīng)道床產(chǎn)生的壓降，是由散體道床對空氣流體的黏性阻力和慣性阻力作用所致。同時，道床內(nèi)部也存在所受壓強異常變化的道砟，主要分布在砟肩和邊坡坡底等區(qū)域，這是因為該區(qū)域道床密實程度較道床內(nèi)部其他區(qū)域更低，流場作用下顆粒會產(chǎn)生不同幅度的錯動；同時，道砟顆粒的形狀、尺寸和所處狀態(tài)各不相同，其受到流場作用的影響也存在差異，導致局部區(qū)域道砟所受風壓的異常變化。

道床多孔介質(zhì)CFD模型計算得到的線路中心橫、縱斷面風壓云圖見圖9。由圖9可知，由于將散體道床考慮為多孔介質(zhì)結構，導致道床橫向流場的分布較為均勻，無大幅度的波動情況。在沿線路縱向上，流場的分布和強度特征與CFD-DEM耦合模型的計算結果較為一致，均存在均勻過度的風壓梯度，且在計算域入口附近風壓為正，出口附近為負，在該風壓梯度作用下產(chǎn)生流場壓降。當輸入風速為10 m/s時，道床多孔介質(zhì)模型的流場壓降曲線見圖10，由圖10可知，該風速條件下，流體流經(jīng)道床多孔介質(zhì)模型前后的壓降為7 218 Pa，與采用CFD-DEM耦合模型得到的流場壓降僅相差2.82%。

圖9 道床多孔介質(zhì)模型截面風壓云圖

圖10 多孔介質(zhì)道床模型壓降曲線

計算輸入風速分別為5、10、15、20、25、30 m/s條件下的道床流場壓降。在不同輸入風速條件下，采用有砟道床CFD-DEM耦合模型，以及基于CFD仿真方法的道床多孔介質(zhì)模型達計算收斂后，流經(jīng)道床的流體壓降對比情況見表2。

表2 道床耦合模型與多孔介質(zhì)模型壓降對比

由表2可知，基于道床多孔介質(zhì)CFD模型計算得到的流場壓降，較CFD-DEM耦合模型的計算數(shù)據(jù)非常相近，二者最大誤差不超過5%，即可認為采用CFD-DEM耦合模型計算得到的道床參數(shù)，可應用在基于道床多孔介質(zhì)模型的CFD仿真分析中。同時，采用道床多孔介質(zhì)模型大幅提高了仿真計算效率，為基于大尺度高速列車-有砟軌道流場特征和運動特征的道砟飛濺機理研究，提供了重要的技術手段。

4 結論

本文將有砟道床考慮為能夠反映其對流場黏性阻力和慣性阻力影響的多孔介質(zhì)模型，基于提出的有砟道床CFD-DEM耦合分析方法，對道床多孔介質(zhì)模型關鍵參數(shù)進行獲取和驗證。

(1)有砟軌道CFD-DEM耦合分析方法，既能反映道砟的級配、形狀、位置等顆粒特征，又能體現(xiàn)道床與流體的相互作用關系，實現(xiàn)了有砟軌道結構的準確模擬。

(2)應用CFD-DEM參數(shù)獲取方法建立的有砟道床多孔介質(zhì)模型，能夠準確模擬有砟道床對流體的黏性阻力和慣性阻力影響，模擬最大誤差僅為4.68%，實現(xiàn)了有砟軌道流場特征的真實還原。

(3)基于CFD-DEM耦合的道床多孔介質(zhì)模型參數(shù)獲取及驗證方法，為研究大尺度高速列車-有砟軌道的道砟飛濺流場特征和運動特征提供了技術手段。