基于Mixture兩相流模型和準靜態方法的積雪分布數值模擬

薛銘乾 黃慧萱

(西南交通大學土木工程學院，四川 成都 610031)

?

基于Mixture兩相流模型和準靜態方法的積雪分布數值模擬

薛銘乾 黃慧萱

(西南交通大學土木工程學院，四川 成都 610031)

在風致雪運動的研究中，采用基于Fluent提供的Mixture兩相流模型，并結合準靜態方法，對二維防雪欄開展數值仿真，分析了防雪欄離地面的高度對地面積雪分布的影響，驗證了該方法的可行性與準確性。

風致雪漂移，數值模擬，兩相流，準靜態方法

0 引言

風致雪運動的數值模擬開始于20世紀90年代。2002年Alhajraf[2]使用Flow-3D軟件對復雜建筑的周邊積雪進行了數值模擬，并使用了FAVOR技術對積雪邊界進行調整；2009年Thiis[3]使用了兩相流的方法進行了三維大跨曲面屋蓋的積雪模擬，在模擬中采用了瞬態計算方法；2011年莫華美[4]采用兩相流理論對一些典型屋面進行了二維的數值模擬，在模擬中使用Fluent中的Mixture模型，并使用了Fluent的動網格技術；2016年周晅毅[5]提出了采用準靜態方法來模擬。

本文采用Euler-Euler兩相流模型，使用Fluent中的Mixture模型完成數值模擬。由于瞬態計算十分耗時，尤其是三維模型的模擬，因此本文采用準穩態方法模擬，并考慮了積雪休止角的影響[6]。

1 風致雪漂移數值模擬方法

在風致雪運動的研究中，通常采用Euler-Euler兩相流理論，本文采用Fluent兩相流模型中的Mixture模型，并配合雪邊界調整。

1.1 Mixture模型

在Fluent中，Mixture模型是一種簡化的多相流模型，是通過求解混合相的連續方程、動量方程和能量方程，次相的體積分數運輸方程和各相之間的相對速度來模擬多相流的運動。

體積分數方程：

1.2 湍流模型

在現階段的研究里，湍流的模擬上沒有定論，然而湍流模型事關風場模擬的精度與準確度，因此湍流模型的選擇十分關鍵，也依然是現在的研究熱點。由于k-ε模型具有較好的普適性，因此往往是工程計算的普遍選擇，因此本文采用Realizablek-ε湍流模型。

1.3 雪的沉積與積雪的侵蝕

雪的沉積與侵蝕和積雪邊界上的剪切風速u*與剪切風速閾值u*t有關，同時也和流場內雪相濃度相關，并忽略溫度、壓強等次要因素。當地的剪切風速由流場特性確定，剪切風速閾值的相關因素就相對復雜，自然條件下雪的剪切風速閾值大致在0.15～0.40之間。當剪切風速u*大于剪切風速閾值u*t時，表現為積雪的侵蝕，反之則為雪的沉積，積雪侵蝕通量qero和雪的沉積通量qdep分別為：

因此積雪邊界上的通量改變值qs=qero+qdep，從而可以得到積雪邊界高度改變量Δh：

其中，Δt為時間步長；γ為積雪中雪相的最大體積分數，一般取值為0.62。

2 算例分析

防雪欄可以改變局部空氣流動規律，可以改變積雪在防雪欄附近的堆積形態，常用于鐵路和公路的風雪防護中。由于防雪欄在長度方向遠大于寬度、高度，因此可以被視為二維問題。

2.1 幾何建模及計算參數

該算例來自于Uematsu(1991)的數值模擬，其中的實測數據來自于Takeuchi(1989)。如圖1所示，其計算流域大小為40m×140m，其中防雪欄上游長度為40m，下游長度為100m，防雪欄的尺寸為高3.41m、寬0.1m。設置兩組模型，控制防雪欄到地面的高度，分別為0.1m和0.4m，其他參數一致。

雪的密度ρs=150kg/m3，雪顆粒的下降速度wf=0.20m/s，雪粒子直徑取為0.15mm，剪切速度閾值u*t=0.20m/s。

由于該模型模擬的是在防雪欄作用下的地面積雪，因此在建模及計算中沒有考慮休止角，不預鋪雪，主要考慮雪的沉積效應。

2.2 邊界條件

1)入流邊界。

使用速度入口邊界，大氣邊界層的水平風速呈對數分布，v=w=0：

其中，κ為馮卡門常數；z0為氣動粗糙長度。

雪相體積分數的入口條件中，需要區分躍移運動和懸移運動，躍移層的高度為：

入口處躍移層和懸移層的雪相體積分數分別為：

2)出流邊界。

采用完全發展的出流邊界。

3)計算流域的頂面與側面。

設置對稱邊界條件。

4)地面及防雪欄表面。

采用無滑移邊界。

2.3 計算方法及流程

在自然條件下，風吹雪現象會持續幾個小時到數天不等，若使用常規手段進行非定常計算，計算量極大且非常耗時。因此考慮到計算的精度及效率，而采用準靜態方法。該方法的計算流程如下：

1)根據計算模型參數及其他信息確定數值模擬時間T，并分成n個時間段，第i個時間長度表示為Δti。

2)通過Fluent多相流模型中的Mixture模型計算風場及其雪相濃度分布，從而得到積雪表面的剪切風速u*和雪相體積分數f。

3)由步驟2)計算得到的結果，并根據公式確定雪相通量的改變值qs，并計算積雪邊界的高度改變量Δhi。

4)將步驟3)計算所得的高度改變量Δhi更新到計算網格中，其中包括調整邊界及邊界附近網格，形成下一個時間段的計算網格。

5)重復步驟2)～3)，完成整個數值模擬。

在這個方法中，采用有限個定常計算來取代非定常計算，在計算效率上來說是可取的。風吹雪現象持續的時間很長，邊界上的改變量是一個有限的量，當邊界高度改變量較小時，邊界改變引起的流場變化可以忽略，因此將整個模擬分成n個階段，在計算精度上是可以接受的。

準靜態方法里的核心問題在于階段數目和每個階段持續時間的確定。階段數目取決于模擬的總時長以及流場的穩定程度。由于在每個階段內被當成定常問題來處理，因此倘若流場對邊界變動較為敏感時，單階段的模擬時長不宜過長。

3 計算分析

考慮防雪欄距離地面高度對積雪形態的影響，其中離地高度h分別為0.1m和0.4m，其中第二個模型為驗證算例。

圖2a)為防雪欄距離地面h0=0.4m情況下初始階段的水平風速分布，防雪欄后部存在一個較大的回流旋渦，在防雪欄后部較遠(x>15)的地面上存在更大的回流旋渦。由于防雪欄距離底部還有0.4m，有大量氣流通過，風速較大。圖2b)則為防雪欄距離地面 0.1m的水平風速分布，兩者風速分布規律相似，但存在少量差異：防雪欄后部的旋渦比前者小一些；地面上的回流旋渦位置向前挪動了(大概位置在x>10)；由于防雪欄距離地面的高度變小，底部通過的氣流變少，底部風速比前者低。

將整個數值仿真過程分成三個階段，分別做兩次網格調整。

圖3為驗證算例的數值仿真結果，取防雪欄高度H=3.41m為特征高度，對坐標X和積雪高度S做歸一化(同下)。從三個階段的仿真結果來看，對比分析三個階段中積雪堆積位置以及積雪增量，認為流場中的風場對積雪邊界的變化不敏感。

對比數值仿真和實測結果可以發現，防雪欄上游，數值仿真的積雪堆積與實測結果吻合的較好，而防雪欄后部則吻合的較差。防雪欄上游的空氣流動沒有受到擾動，流場的仿真結果與實際相差不大，但下游則受到防雪欄的擾動，空氣流動變得復雜，脈動特征增強，湍流模型難以捕捉流場特征，因此流場的模擬效果相對較差；同時由于Takeuchi等人所使用的防雪欄是存在14%的孔隙率，然而在數值仿真中并未考慮，因此下游區域的流場仿真與實際有較大出入。在模擬中，防雪欄下游區存在積雪堆積，大致位于2.8≤x/H≤4.5，其規律與實測相似?？傮w上而言，數值仿真結果可以接受，因此認為本文提出的數值仿真方法可行。

圖4為防雪欄距離地面0.1 m的三階段數值仿真結果，其中反映出來的規律與驗證模型一致，存在少許差異(見圖5)：防雪欄上游，兩者差距不大，但是靠近防雪欄位置的積雪坡度隨高度h0的減小而增大；防雪欄下游的積雪堆積位置要比驗證模型靠前，大致位置在1.3≤x/H≤3.0。

4 結語

本文提出了一種基于Fluent的Mixture兩相流模型和準靜態方法的數值模擬方法，能夠提高風致雪漂移數值模擬的計算效率。通過對二維防雪欄地面積雪分布的數值仿真，驗證這種數值模擬方法的可行性，并分析防雪欄距離地面高度不同所帶來的差異。本文得出如下幾個結論：

1)將二維防雪欄的數值仿真結果與文獻實測數據進行對比，驗證了這種數值模擬方法的可行性。

2)隨著防雪欄距離地面的高度的減小，上游的積雪堆積變陡，下游的積雪堆積向前移動。

[1] Sato T,Uematsu T,Nakata T,et al.Three dimensional numerical simulation of snowdrift[J].Journal of Wind Engineering and Industrial Aerodynamics,1993(46):741-746.

[2] Alhajraf S.Numerical simulation of sand and snow drift at porous fences[A].Proceedings of the Fifth International Conference on Aeolian Research and the Global Change and Terrestrial Ecosystem-Soil Erosion Network[C].2002:208-213.

[3] Thiis T K,Potac J,Ramberg J F.3D numerical simulations and full scale measurements of snow depositions on a curved roof[A].5th European & African Conference on Wind Engineering,Florence,Italy[C].2009.

[4] 莫華美.典型屋面積雪分布的數值模擬與實測研究[D].哈爾濱：哈爾濱工業大學,2011.

[5] Zhou X,Kang L,Gu M,et al.Numerical simulation and wind tunnel test for redistribution of snow on a flat roof[J].Journal of Wind Engineering and Industrial Aerodynamics,2016(153):92-105.

[6] 康路陽,周晅毅,顧 明.考慮積雪休止角的屋面積雪漂移數值模擬方法[J].同濟大學學報(自然科學版),2016(1):11-15.

Numerical simulation of snow distribution based on Mixture two-phase flow model and quasi-static method

Xue Mingqian Huang Huixuan

(School of Civil Engineering, Southwest Jiaotong University， Chengdu 610031， China)

In the research of wind-induced snow movement, a numerical simulation method based on Mixture two-phase flow model from Fluent combined with quasi-static method is adopted. Through the numerical simulation of two-dimensional snow-proof bars, this method is tested and the influence of the height of the snow protection bar on the ground snow distribution is analyzed.

snowdrift, numerical simulation, two-phase flow, quasi-static method

1009-6825(2017)16-0041-03

2017-03-22

薛銘乾(1991- )，男，在讀碩士； 黃慧萱(1962- )，男，副教授

TP319

A