基于Mixture兩相流模型和準靜態(tài)方法的積雪分布數(shù)值模擬

薛銘乾 黃慧萱

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

?

基于Mixture兩相流模型和準靜態(tài)方法的積雪分布數(shù)值模擬

薛銘乾 黃慧萱

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

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

風致雪漂移，數(shù)值模擬，兩相流，準靜態(tài)方法

0 引言

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

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

1 風致雪漂移數(shù)值模擬方法

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

1.1 Mixture模型

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

體積分數(shù)方程：

1.2 湍流模型

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

1.3 雪的沉積與積雪的侵蝕

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

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

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

2 算例分析

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

2.1 幾何建模及計算參數(shù)

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

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

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

2.2 邊界條件

1)入流邊界。

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

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

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

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

2)出流邊界。

采用完全發(fā)展的出流邊界。

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

設置對稱邊界條件。

4)地面及防雪欄表面。

采用無滑移邊界。

2.3 計算方法及流程

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

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

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

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

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

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

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

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

3 計算分析

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

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

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

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

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

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

4 結語

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

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

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] 莫華美.典型屋面積雪分布的數(shù)值模擬與實測研究[D].哈爾濱：哈爾濱工業(yè)大學,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] 康路陽,周晅毅,顧 明.考慮積雪休止角的屋面積雪漂移數(shù)值模擬方法[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