基于DPM模型的高速列車轉向架區域風雪流數值模擬

謝菲，高廣軍，何侃，陳功

(1.湖南鐵道職業技術學院，湖南 株洲 412001；2.中南大學 交通運輸工程學院，湖南 長沙 410075；3.中車株洲機車電力機車有限公司，湖南 株洲 412001)

高速列車在寒冷的冬季運行時，軌道上的積雪會受到列車風影響，在列車的底部和尾部飛揚[1]。此外，列車轉向架上的電機、制動鉗夾等設備在運行時會產生熱能。當空氣中運動的雪粒子接觸這些發熱設備時，將會融化為水，進一步在高寒低溫條件下凝結成冰。因此，列車在長期運行過程中，積雪會在轉向架區域持續堆積，最終形成厚重的冰雪混合物，如圖1所示。這一現象將增加轉向架配重，危害轉向架各關鍵部件的工作性能，進而對列車行車安全產生極大隱患。針對列車積雪結冰問題，JEMT[2]采用對環境無害的加熱丙二醇作為材料加入除冰設備中，來改善瑞典等北歐國家的軌道積雪問題。PAULUKUHN等[3]針對俄羅斯冬季低溫、多雪環境，為保證列車能夠經得起無法避免的冰雪侵害，在除冰后的列車底部噴射水和丙烯類混合液，從而減少結冰。BETTEZ[4]針對日本鐵路線的列車運行對高效除冰的要求，采用噴灑器作為除冰裝置，有效解決了因積雪結冰而導致列車延誤的問題。綜上所述，來自于冬季處于高寒環境國家的研究人員提出了一系列的列車防積雪應對策略。然而，針對轉向架區域的積雪結冰研究鮮有文獻提及。高廣軍等[5-8]采用化學發熱材料或物理結構改進等方式來改善轉向架區域的積雪情況。丁叁叁等[9-11]均對列車轉向架區域積雪結冰問題開展了流場數值模擬研究。然而，這些研究僅從流場進風量的角度間接預測轉向架區域的積雪情況，并未考慮雪粒子在空氣中的分布。雪粒子雖對流場運動具有一定跟隨性，但氣流出現擾動或者流速降低時，雪粒子的運動慣性將使得其對流場的跟隨性減弱。在轉向架區域，由于復雜幾何模型影響，氣流較為紊亂，存在大量湍流結構，且流速相對車身周圍流場較低。與真實情況相比，僅從流場的角度考慮轉向架區域積雪分布有失偏頗。因此，本文試圖從風雪兩相流的角度探索轉向架區域積雪分布情況。關于風雪兩相流，FUJII 等[12]針對日本鐵路積雪情況進行研究，采用實驗的方式分析了雪粒子運動的臨界風速與雪密度及含水量之間的關系，提出雪粒子初始密度范圍為36～132 kg/m3。TOMINAGA 等[13]采用數值仿真對雪粒子的堆積、侵蝕等情況進行模擬，并通過實驗驗證得出雪粒子的密度為100 kg/m3時數值仿真與實驗結果具有良好的一致性。桑嘉賓[14]建立風雪流的數值模型，研究了混合粒徑，得出雪粒子直徑的影響因素，指出雪粒子直徑一般分布在0～0.6 mm 之間，且主要集中在0.2 mm 附近。莫華美[15]針對屋面積雪分布建立了風雪流數值模型，利用該數值模型對屋面上的積雪分布規律進行數值模擬，指出雪粒子直徑主要分布在0～0.4 mm 之間。TOMINAGA 等[16]采用CFD 數值模擬仿真的方法模擬雪粒子的堆積深度，并且采用風洞試驗驗證了模擬方法的可行性。BEYERS等[17]采用數值仿真模擬了雪粒子在建筑物周圍逐步形成堆積的過程，有助于在建筑設計時對潛在的積雪情況進行處理。SERINE 等[18]將兩維空間上的雪粒子數值模擬方法和風洞試驗對比，并將其推廣到三維空間上，從而確定雪粒子數值仿真的可行性。張強等[19]基于歐拉兩相流理論和擴展傳熱模型，提出一種模擬飛機機翼積雪結冰的方法，并在實驗中得到驗證。GORDON 等[20-21]針對嚴寒天氣的積雪現象，對雪粒子的形狀、大小、速度及密度等進行研究，得出雪粒子的相關物理性質。SATO 等[22]針對在強風情況下的雪粒子破碎問題，采用風洞試驗對積雪的密度進行研究，并在實地觀測結果中驗證其合理性。ZHOU 等[23]采用數值模擬技術對機場航站樓屋頂的雪粒子現象進行模擬，并用風洞試驗對比驗證數值模擬技術的可靠性，從而得出湍流模型、風速及風向等對積雪的影響。BEYERS等[24]采用數值方法模擬立方體周圍的積雪情況，并通過實驗進行對比驗證。韓俊臣等[25]采用三維非定常雷諾時均模擬高速動車組轉向架區域流場和雪粒子分布情況。婁振[26]采用準瞬態的方法模擬高速列車拖車轉向架區域的積冰過程，研究轉向架區域的積冰分布特征。上述文獻表明，采用風雪兩相流數值模擬仿方法可以再現積雪分布實際情況，具有較高的可行性。因此，本文采用數值模擬方法，復現列車轉向架區域的流場流動情況，結合DPM 模型預測了雪粒子在空氣中的運動軌跡，揭示了雪粒子在轉向架上的分布情況，并從雪粒子隨流場運動的角度描述了轉向架區域積雪形成的原因，為高速列車轉向架區域積雪結冰研究提供參考。

圖1 轉向架區域冰雪混合物堆積Fig.1 Accumulation of the ice and snow mixture in the bogie area

1 研究方法

1.1 物理模型、計算域及邊界條件設置

本文基于轉向架區域流場，建立了如圖2所示的帶有單個轉向架區域的高速列車車體模型。考慮到實際運行中，頭車的結構和形狀與其他車不一致，不具有普遍性。因此，選擇中間車型作為研究對象，在實踐中能更好地估計一般轉向架上的積雪結冰情況。該簡化模型忽略了列車頭、尾車等流線型結構對真實列車宏觀流場產生影響的因素，進而著眼于轉向架區域。車體的高度為3.6 m，寬度為3.2 m，長度為17.6 m。轉向架區域距離車體迎風側端面距離為6.9 m，使得車體底部流場在進入轉向架區域之前得到充分發展，符合實際運行的列車轉向架區域流場情況。同時，轉向架區域距離車體背風側端面同樣為6.9 m，使得轉向架區域的尾流得到充分發展，避免尾流對轉向架區域流場產生干擾。此外，為真實模擬高速列車運行時車體底部的流場結構，增加了離地高度為0.2 m的鋼軌模型。

圖2 車體模型Fig.2 Geometry of the train

為模擬車體周圍流場，建立了長60 m，寬20 m，高20 m 的計算區域，如圖3 所示。車體上游處(平面ABCD)設置為速度入口，給定速度為55.56 m/s，使得雷諾數為1.325×107。車體下游處(平面EFGH)設置為壓力出口，給定壓力為0 Pa。車體兩側(平面ADHE和平面BCGF)及頂面(平面ABFE)均設置為固定壁面，模擬風洞環境。地面(平面DCGH)設置為滑移地面，給定與車速相同的滑移速度。車體及轉向架壁面設置為固定壁面。車體距離氣流入口15 m，車體兩側與計算域側邊界距離等寬。其中，模型的阻塞比約為2.8%，少于5%，滿足一般鈍體車輛計算需求[27]。計算時間步長設置為0.001 s，使得庫朗數平均值小于1。

圖3 計算區域Fig.3 Domain for simulation

1.2 網格設置

本文數值計算中采用六面體主導的網格。計算域采用開源計算流體力學分析軟件OpenFoam 進行SnappyHex 網格離散，網格單元總數2 000 萬。網格采用切割單元(cut-cell)方式，以級數遞增方式從計算域邊界向車體壁面逐層加密。車體壁面網格尺寸為0.02 m，轉向架壁面網格尺寸為0.01 m。為準確模擬邊界層流場，在車體和轉向架壁面上布置了8 層附面層網格，近壁面第1 層網格高度為0.6 mm，使得壁面無量綱化尺寸y+處于30～300 之間。圖4展示了轉向架表面網格以及車身周圍計算域剖面網格。

圖4 數值計算網格Fig.4 Mesh for simulation

為了明確本文數值模擬不受網格密度的影響，基于相同的壁面法向y+，在流向和展向使用了不同的網格尺度，并采用了中等密度網格和精細密度網格進行對比。其中，中等網格和精細網格數量分別是2 000萬和4 000萬。如圖5所示，為從車體底面中心線的2 個網格獲得的時均壓力系數Cp。在圖中可以發現中等網格和精細網格的壓力系數基本保持一致，說明本文采用的中等網格已經達到計算精度的收斂區間，具有較好的準確性。

圖5 車體底面中心線壓力系數對比Fig.5 Time-averaged Cp along the centerline of the bottom surface of the train body

1.3 流場數值模擬方法

針對高速列車周圍流場，HE 等[28]指出采用大渦模擬(LES)能較好地捕捉湍流特性，但采用改進的延遲分離渦(IDDES)或者非穩態雷諾平均(URANS)方法同樣能較為準確地模擬流場速度及壓力分布?？紤]到LES 或IDDES 對于網格分辨率要求較高，不適應于較高雷諾數條件下的數值仿真。此外，本文研究內容除了對流場運動方程求解之外，還需要額外求解雪粒子運動方程，計算量較大。因此，本文采用對網格尺度要求更低、在工程尺度得到廣泛認可及應用的URANS 方法對流場部分進行求解。Realizablek-?湍流模型能較好地捕捉流動分離，保持雷諾應力與真實湍流一致，使計算結果更符合真實性。因此，為適配壁面網格尺寸，選用了Realizablek-?湍流模型，其控制方程見參考文獻[29]。

1.4 雪粒子數值模擬方法

風吹雪是一種氣固兩相流，根據現代多相流計算理論，計算模型主要有歐拉兩流體模型、混合模型以及拉格朗日粒子軌跡模型等[30]。拉格朗日粒子模型是另一種氣固兩相流計算方法，目前風吹雪數值模擬中應用較少。拉格朗日模型的優點是物理概念明確、結果直觀，可計入粒子的軌道經歷效應和歷史效應，能得到粒子的詳細運動信息[31]。因此，本文采用拉格朗日隨機模型方法，對轉向架區域的積雪分布進行數值模擬。

在本文中，通過對粒子的力平衡方程進行積分，來預測離散相粒子的運動軌跡。其中，作用在粒子上的外力主要包括氣動力、重力、曳力等[32]。除此之外，其余力對粒子運動影響相對較小，在分析中可忽略不計[30]，其運動方程見式(1)：

對方程(1)積分就得到了顆粒軌道上每一個位置上的顆粒速度，見式(2)：

其中：x為粒子位移矢量，t為時間。粒子的軌跡是在離散的時間步長上逐步進行積分運算得到的，對式(2)進行積分，就可以得到粒子的軌跡。

對于車體、輪對等不易積雪部件，設定雪粒子邊界為反射條件；對于構架、電機、端板、制動鉗夾等易積雪部件，設定雪粒子邊界為捕捉條件，認為轉向架區域雪粒子全部被捕捉，以此來考慮積雪最嚴重情況。其中，雪粒子密度為100 kg/m3，雪粒子直徑采用均勻分布為0.2 mm，雪粒子材料為water-liquid，列車運行速度為200 km/h。

此外，時間步長Δt由庫朗數Co最大值控制，定義為Co=vΔt/Δx。其中，Δx是網格尺寸，v是網格內的流體速度。此外，隱式算法中的庫朗數范圍為0～100。如圖6 所示，在本文數值模擬計算中，庫朗數分布為0～11.85，且絕大部分集中在0和1 之間，基本滿足計算要求。當t=2.00 s 時，氣流多次通過該區域(風速55.56 m/s)，轉向架中的風雪流場得到相對穩定和充分發展。并且，基于課題組多年計算轉向架積雪的經驗，2 s 的時間足以使得風夾帶的雪粒子數量在轉向架區域達到相對穩定。因此，本文的計算時間設置為2 s，時間步長為0.001 s，共計2 000步。

圖6 庫朗數等值線Fig.6 Contours of the courant number

2 數值模擬方法驗證

關于本文中流場數值模擬方法的驗證請參考文獻[33]。本文僅針對數值模擬中雪粒子在轉向架區域的分布情況進行驗證分析，開展離散相風洞試驗。該試驗在中南大學高速列車研究中心的風洞實驗室高速段進行。實驗段的尺寸為3.4 m×1.0 m×0.8 m(長×寬×高)。由于真實雪粒子存在難制備、難保存等問題，試驗中采用了與雪粒子密度相似的木屑粒子(密度約為100 kg/m3)作為替代物。木屑均經過篩網嚴格過濾，挑選出與雪粒子直徑(0.2 mm)大小類似的粒子進行試驗。試驗段主要安裝的裝置包括高速列車原始轉向架模型、高速列車原始車體模型、固定裝置，其中固定裝置又包括支撐桿、支撐座、分流板以及安裝架等構成，如圖7所示。

圖7 模型的安裝與固定Fig.7 Installation and fixation of the model

本次兩相流風洞實驗需要模擬雪粒子在轉向架區域各關鍵部件的黏附情況，為此，在轉向架和車體模型表面涂抹蜂蜜以黏附吹過的木屑。如圖8所示，為涂刷蜂蜜后的轉向架以及轉向架各個細節部位的示意圖。

圖8 涂刷蜂蜜后的轉向架Fig.8 Bogie painted with honey

為驗證數值模擬方法的可靠性，且保障計算與試驗的一致性，依據風洞試驗室和試驗模型的大小，建立了和風洞試驗模型完全相同的數值模擬模型，且計算入口風速與風洞試驗自由流風速一致。網格離散原則及數值模擬方法與本文第1節描述一致，由于簡化了車體模型，網格總數為1 400萬。

圖9展示了風洞試驗和仿真計算的轉向架粒子堆積情況。從圖9可以看出，風洞試驗中木屑分布特點與數值仿真計算的積雪分布規律一致，粒子的堆積均產生在轉向架各部件的迎風側。因此，本文采用的數值仿真計算方法能較為準確地模擬雪粒子在轉向架區域的分布情況。

圖9 轉向架底部表面積雪結果對比Fig.9 Comparison of snow results on the bottom surface of bogie

3 數值仿真結果分析

本節采用第1節所述流場及雪粒子數值模擬方法，針對轉向架區域進行數值模擬仿真，得出轉向架區域壓力分布、流線分布、雪粒子分布及轉向架表面積雪分布情況。

3.1 轉向架區域壓力分布

針對轉向架區域主要部件的積雪結冰現象進行分析，在轉向架區域作縱截面空間切片，如圖10 所示，切片1 和4 在轉向架輪對、制動鉗夾位置，切片2和3在轉向架電機位置。

圖10 轉向架區域切片位置Fig.10 Position of slices in the bogie area

如圖11 所示，展示了轉向架區域時均壓力云圖，來流方向為從左至右。由于切片1 和4 較為對稱，為簡略表達僅展示切片1，2 和3 上的時均壓力分布。從圖11 可以看出，輪對、制動鉗夾、電機等突出部件在迎風側由于受氣流沖刷出現明顯正壓，而在這些部件的背風側則為負壓。轉向架區域的壓力分布體現了復雜幾何結構對流場分布具有較大影響。因此，在進一步探討雪粒子軌跡之前，需對流線分布進行分析。

圖11 轉向架區域壓力分布Fig.11 Pressure distribution in the bogie area

3.2 轉向架區域流線分布

圖12 展示了轉向架區域流線分布，來流方向為從左至右。從圖12 可以看出，車體底部氣流在轉向架區域發生了偏轉，部分氣流在流經轉向架前端時向上揚起進入轉向架區域，并在轉向架區域各個部件處形成多處漩渦，說明轉向架區域流場具有較高湍流度。此外，氣流進入轉向架區域之后受復雜幾何結構影響，速度明顯降低，尤其在轉向架輪對、制動鉗夾以及電機等部位的流速接近為0。當流速降低時，雪粒子跟隨性降低，當流速低于雪粒子起動風速時，重力占雪粒子所受合力比重增大。因此，當雪粒子隨氣流上揚進入轉向架區域時，更易在低速區域附著至轉向架各部件上形成堆積。值得注意的是，在轉向架后端，部分氣流在設備艙端板受阻，再次回流入轉向架區域。

圖12 轉向架區域流線分布Fig.12 Streamline in the bogie area

3.3 轉向架區域雪粒子分布

通過對轉向架區域凈風場流場情況的分析，已經初步判斷出轉向架區易形成雪粒子堆積區域。為進一步直觀體現雪粒子在轉向架區域的分布情況，圖13 展示了雪粒子運動的快照，體現了雪粒子的運動歷程。從圖13 可以看出，在t=0.12 s 時刻，雪粒子進入到轉向架區域的前部；在t=0.16 s時刻，雪粒子進入到轉向架區域的中部，且集中在高度較低區域，對制動鉗夾、電機等發熱元件的沖刷不明顯；在t=0.20 s 時刻，雪粒子經過整個轉向架區域，更多的雪粒子進入到轉向架的內部；當時間到達t=2.00 s 時刻，經過更長物理時間后，大量雪粒子進入轉向架區域。此后，雪粒子分布趨于穩定。

圖13 雪粒子運動時間序列快照Fig.13 Snap shots of the movement of snow particles

整體而言，在列車運行過程中，大量雪粒子跟隨氣流從列車車體底部流入，盡管大多數雪粒子到達并且通過了轉向架區域，仍然有許多雪粒子隨著上揚的空氣進入了轉向架區域。在轉向架區域前方，雪粒子跟隨氣流運動，具有一定的跟隨性；在轉向架區域后方，雪粒子在運動慣性作用下和氣流發生偏離。進入轉向架區域的雪粒子流速降低、跟隨性減弱，因此被轉向架區域各部件捕捉。隨著時間推移，雪粒子將會持續進入轉向架區域，從而在各部件上產生堆積。

3.4 轉向架表面積雪分布

如圖14 所示，為展示轉向架區域主要積雪部位，本文顯示了轉向架區域積雪濃度大于1×10-12kg/m3以上時的積雪覆蓋情況，并采用速度渲染顯示沖擊到轉向架區域時雪粒子的速度情況。圖14 中，來流方向為從左至右。可以看出，在轉向架底部區域產生大量積雪，且積雪更易堆積在轉向架各部件迎風側，而在各部件背風區域積雪相對較少。此外，積雪主要集中在電機、制動鉗夾、輪對以及端板等部件底部迎風區域，而在轉向架部件頂部和背風區域堆積相對較少。

圖14 轉向架表面積雪分布Fig.14 Snow accumulation on the surface of the bogie

4 結論

1) 數值計算結果與風洞試驗觀測結果具有較高一致性，說明本文數值模擬方法準確可靠，適用于轉向架區域積雪結冰研究。

2) 空氣流經轉向架迎風側時會向上揚起，進入轉向架的輪對、電機、制動鉗夾等關鍵部位，空氣流經轉向架背風側時會沖刷轉向架區域的設備艙后端板。

3) 雪粒子運動軌跡和積雪堆積主要集中在輪對、電機、制動鉗夾以及端板等部件底部迎風區域，而在轉向架頂部和背風區域堆積較少。

本文的數值計算方法可適用于模擬高速列車在風雪天氣下運行時轉向架區域積雪結冰情況，對于高寒動車組防積雪結冰研究具有一定參考價值。