建筑屋蓋表面及其周圍積雪分布研究

徐 楓， 周高照， 肖儀清， 歐進(jìn)萍,2

(1.哈爾濱工業(yè)大學(xué) 深圳研究生院，深圳 518055；2.哈爾濱工業(yè)大學(xué) 土木工程學(xué)院，哈爾濱 150090)

風(fēng)吹雪是氣流挾帶起分散的雪粒在近地面運(yùn)行的多相流，對自然積雪有重新分配的作用，它所形成的積雪深度可達(dá)到自然積雪厚度3倍～10倍。在一些強(qiáng)降雪地區(qū)，雪荷載是結(jié)構(gòu)設(shè)計(jì)的控制荷載，降雪后由于風(fēng)力作用雪顆粒發(fā)生復(fù)雜的漂移堆積運(yùn)動造成屋面局部積雪過大，從而導(dǎo)致梁、柱扭曲，屋面倒塌的事件層出不窮[1-4]。另一方面，建筑周邊積雪不均勻分布不僅會影響行人和通車，對一些低矮建筑甚至有被積雪覆蓋的危險。因此，準(zhǔn)確預(yù)測建筑屋蓋及其周圍積雪的空間分布是預(yù)防積雪成災(zāi)必不可少的工作。

預(yù)測積雪分布主要有三種方法：現(xiàn)場實(shí)測、風(fēng)洞試驗(yàn)和數(shù)值模擬。現(xiàn)場實(shí)測是進(jìn)行風(fēng)吹雪研究的最直接、可靠的方法，Oikawa[5]和Tsuchiya等[6]分別對邊長為1 m立方體周邊和高低屋屋蓋表面積雪分布進(jìn)行了實(shí)測，給出了積雪的空間分布形式，但由于受自然條件(地理環(huán)境、風(fēng)向等)的約束，很難進(jìn)行大規(guī)模、系統(tǒng)性的研究。風(fēng)洞試驗(yàn)是進(jìn)行風(fēng)吹雪研究的重要研究手段，Iversen[7]對防雪柵周邊風(fēng)吹雪進(jìn)行了試驗(yàn)并與實(shí)測結(jié)果進(jìn)行了比較，Sant’Anna等[8-9]對平屋蓋表面積雪分布進(jìn)行了風(fēng)洞實(shí)驗(yàn)，劉慶寬等[10]分別采用食用鹽、硅砂、明礬、干松木屑代替雪顆粒模擬了高低屋蓋表面積雪分布，其中食用鹽的模擬結(jié)果與Tsuchiya的實(shí)測結(jié)果吻合較好。風(fēng)洞試驗(yàn)雖然彌補(bǔ)了現(xiàn)場實(shí)測受自然條件限制這一難題，但其成本高、周期長，而且很多參數(shù)相似比受到限制。然而，數(shù)值模擬可以方便的調(diào)節(jié)各種研究參數(shù)、成本低、周期短、能實(shí)現(xiàn)全尺度模擬，隨著計(jì)算機(jī)性能的迅速提高，基于CFD技術(shù)的風(fēng)吹雪數(shù)值模擬成為了預(yù)測積雪分布的重要研究方法。Uematsu等[11]首次使用混合長度理論作為湍流模型，通過引入雪質(zhì)量傳輸方程和雪面侵蝕經(jīng)驗(yàn)公式，運(yùn)用CFD軟件模擬了二維雪柵與三維建筑周邊的風(fēng)致雪漂流。隨后Bang等[12]提出了一種新的兩相流飄雪模型—廣義飄雪通量模型，但其未考慮躍移層的影響。Thiis等[13]運(yùn)用該模型對一個簡單建筑周邊的風(fēng)致雪漂移進(jìn)行了模擬，并與全尺寸實(shí)測結(jié)果進(jìn)行了比較。Beyers等[14]運(yùn)用單方程模型和k-ε湍流模型對立方體周邊的積雪分布進(jìn)行了模擬，李雪峰等[15]在Beyers單方程的基礎(chǔ)上考慮了雪顆粒的浮力、慣性力和氣動阻力對空氣相的影響，模擬了立方體周邊積雪分布并與實(shí)測結(jié)果進(jìn)行了比較。孫曉穎等[16]基于 Euler-Euler方法，采用FLUENT中Mixture多相流模型對立方體周邊的風(fēng)雪漂移進(jìn)行了數(shù)值模擬，并與實(shí)測結(jié)果進(jìn)行了對比發(fā)現(xiàn)兩者吻合較好。李躍等[17]采用CFD大型商業(yè)軟件對大跨度球殼屋蓋結(jié)構(gòu)進(jìn)行風(fēng)致積雪分析，并通過比較風(fēng)速、跨度、傾角等因素對積雪侵蝕的影響，擬合出了相應(yīng)的積雪分布系數(shù)。

本文基于歐拉-歐拉方法，假定空氣相和雪相均為連續(xù)相，基于計(jì)算流體力學(xué)軟件FLUENT 17.0，采用Mixture模型和k-kl-ω湍流模型模擬了立方體周邊和高低屋蓋表面的風(fēng)致雪漂移，并與實(shí)測和風(fēng)洞試驗(yàn)結(jié)果進(jìn)行了對比分析。

1 風(fēng)吹雪數(shù)值模型

積雪表面的雪顆粒受到風(fēng)力作用，當(dāng)雪面的摩擦速度達(dá)到閾值摩擦速度時，雪顆粒將克服雪面約束發(fā)生運(yùn)動，從而導(dǎo)致積雪不均勻分布。雪粒運(yùn)動主要有三種形式：蠕移、躍移與懸移。蠕移對風(fēng)吹雪的質(zhì)量輸送貢獻(xiàn)很小，僅可能在局部行成一些小的波紋，風(fēng)吹雪研究中一般忽略蠕移的影響，躍移層高度即躍移層與懸移層的分界面距雪面的豎向高度，躍移層高度公式

(1)

式中:g為重力加速度，取9.81 m/s2。數(shù)值模擬中僅考慮躍移與懸移的影響。

1.1 控制方程

基于Euler-Euler方法的多相流模型有三種，分別為：VOF模型、Mixture模型和Eulerian模型。相比VOF模型，Mixture模型允許空氣相與雪相之間有相對滑移速度和相互穿插，在各相間耦合關(guān)系不明確時，Mixture模型能與Eulerian模型有同等精度，且計(jì)算量小、穩(wěn)定性好。文獻(xiàn)[16]模擬了防雪柵和立方體周邊的風(fēng)致雪漂移，計(jì)算結(jié)果與實(shí)測結(jié)果吻合較好，驗(yàn)證了Mixture模型模擬風(fēng)吹雪運(yùn)動的可行性。

Mixture模型求解空氣相和雪相組成的混合相的連續(xù)方程、混合的動量方程和第二相體積分?jǐn)?shù)方程，從而得到整個風(fēng)場和體積分?jǐn)?shù)場，其中空氣相為主相，雪相為第二相。

(1) 空氣相控制方程

連續(xù)方程

(2)

動量方程

(3)

(2) 雪相控制方程

(4)

(3)k-kl-ω湍流模型

Walters等[18]在k-ω模型基礎(chǔ)上引入了層流動能kl，提出了k-kl-ω三方程渦黏模型，該模型能夠準(zhǔn)確地預(yù)測層流至湍流的轉(zhuǎn)捩過程，其表達(dá)式如下

(5)

式中:PKT,PKL為湍流、層流動能產(chǎn)生項(xiàng);R為層流動能kl轉(zhuǎn)化湍動能kT的產(chǎn)生率;DT、DL為近壁面耗散項(xiàng)。

1.2 積雪侵蝕與沉積

氣流攜帶著雪粒經(jīng)過地面，當(dāng)摩擦速度μ*小于閾值摩擦速度μ*t時，雪顆粒將在此處發(fā)生沉積；當(dāng)摩擦速度μ*大于閾值摩擦速度μ*t時，雪顆粒將被風(fēng)吹走發(fā)生侵蝕，閾值摩擦速度μ*t取值0.15，摩擦速度μ*表達(dá)式

(6)

式中:τ為壁面剪切應(yīng)力，ρa(bǔ)為空氣密度。

雪質(zhì)量通量q(kg·m-1·s-1)指單位時間內(nèi)，豎直剖面單位面積內(nèi)通過的雪的質(zhì)量。積雪侵蝕沉積目前沒有統(tǒng)一的公式，本文采用應(yīng)用較多的Mohamed等[19]提出的侵蝕、沉積通量表達(dá)式

(7)

(8)

式中:wf為雪粒沉積速度;A為常數(shù)，取-7.0×10-4;φs=f·ρs為近壁面雪質(zhì)量濃度。

根據(jù)侵蝕沉積通量可以計(jì)算積雪面厚度的改變，假設(shè)dt時間內(nèi)積雪厚度的改變?yōu)閐h，其表達(dá)式如下

(9)

式中:qs為沉積或侵蝕質(zhì)量通量;γ為積雪堆積率，取0.62。

積雪沉積與侵蝕的發(fā)生將導(dǎo)致積雪邊界發(fā)生變化，從而對計(jì)算域風(fēng)場有一定程度的影響，李雪峰[20]運(yùn)用動網(wǎng)格對二維的高低屋蓋積雪分布進(jìn)行了研究，分析了積雪邊界對積雪分布的影響，結(jié)果表明積雪邊界的變化對結(jié)果的影響不到0.1%。本文由于考慮到三維動網(wǎng)格對計(jì)算機(jī)的要求較高，且計(jì)算時間很長，未考慮積雪邊界的動態(tài)變化。

1.3 入流邊界

1.3.1 空氣相入流邊界

空氣相入流面采用速度入(Velocity-Inlet)邊界，指數(shù)型風(fēng)剖面，即：

(10)

式中:v(z)為離地高度z處平均風(fēng)速，z0為參考高度，一般取10 m;v0為參考高度處平均風(fēng)速;α為地面粗糙系數(shù)，對應(yīng)A、B、C、D四類場地，分別取0.12、0.16、0.22、0.3。

湍流入口均參考王福軍[21]書中公式

(11)

(12)

其中Cμ=0.09，I、l為湍流強(qiáng)度和湍流積分尺度，參考日本建筑荷載規(guī)(2004)[22]取值如下

(13)

(14)

式中:ZG為梯度風(fēng)高度，取450 m;Zb為截?cái)喔叨龋? m;α2取值0.15。

1.3.2 雪相入流邊界

雪相邊界入口給定雪相體積分?jǐn)?shù)豎向剖面，對于躍移層和懸移層給定不同的函數(shù)公式，大多數(shù)是實(shí)測到的經(jīng)驗(yàn)公式。

躍移層的體積分?jǐn)?shù)分布采用Pemeroy等[23]提出的傳輸率公式

(15)

懸移層體積分?jǐn)?shù)采用Pomeroy等[24]理論推導(dǎo)的公式

(16)

式中：z為離地高度；ρs為積雪密度。

2 立方體周邊風(fēng)致積雪模擬結(jié)果

建筑對氣流的阻礙作用改變建筑周邊的風(fēng)場，從而引起積雪在建筑周圍的不均勻分布。建筑周圍的不均勻積雪會給人們出行帶來不便，在降雪嚴(yán)重的地區(qū)，甚至可能出現(xiàn)積雪覆蓋低矮建筑的現(xiàn)象。Oikawa對邊長為1 m立方體周邊的積雪分布進(jìn)行了實(shí)測，本文模擬立方體周圍的積雪分布，并與Oikawa實(shí)測結(jié)果、Tominaga等[25]新模型和孫曉穎等[26]Mixture模型等的數(shù)值模擬結(jié)果進(jìn)行比較，驗(yàn)證了Mixture模型對風(fēng)吹雪模擬的可行性。

2.1 幾何建模與計(jì)算參數(shù)

本文對邊長為1 m的立方體模型進(jìn)行分析，計(jì)算域長、寬、高分別為20 m、10 m、8 m(如圖1)，阻塞率為1.25%，采用結(jié)構(gòu)化網(wǎng)格進(jìn)行離散，邊界層最小網(wǎng)格尺寸為0.005 m，網(wǎng)格總數(shù)為160萬。雪密度取150 kg/m3，雪粒直徑取0.000 15 m，閾值摩擦速度取0.15 m/s。采用Mixture多相流模型，湍流模型為k-kl-ω模型，雙精度求解，控制方程殘差降到10-6。

圖1 計(jì)算域網(wǎng)格

2.2 邊界條件

入流面：采用速度入口邊界條件(Velocity-Inlet)，由于忽略雪相對空氣相的影響，雪相和空氣相采用相同的入口速度，指數(shù)風(fēng)剖面見式(10)，其中參考高度取立方體高度，參考風(fēng)速取5 m/s，地面粗糙系數(shù)α取0.25，湍流動能k和湍流耗散率ω分別由式(11)和(12)給出。雪相體積分?jǐn)?shù)入口分為躍移層和懸移層分別由式(15)和(16)給出，兩者分界面高度即躍移層高度。

出流面：采用充分發(fā)展邊界條件。

積雪面和立方體壁面：采用無滑移wall邊界條件。

計(jì)算域其他邊界：采用對稱邊界。

2.3 網(wǎng)格無關(guān)性驗(yàn)證

網(wǎng)格無關(guān)性的檢驗(yàn)是數(shù)值模擬中非常重要的部分，在網(wǎng)格數(shù)量的增加對數(shù)值模擬結(jié)果的影響可以忽略不計(jì)時，網(wǎng)格無關(guān)性才滿足要求。本文采用4種不同數(shù)量的結(jié)構(gòu)化網(wǎng)格來驗(yàn)證網(wǎng)格無關(guān)性(如表1)，由于研究建筑周邊的積雪分布，采用立方體(邊長H=1 m)迎風(fēng)向中心線地面的摩擦速度進(jìn)行比較，其結(jié)果如圖2所示。

表1 計(jì)算方案

四組Case計(jì)算所得摩擦速度μ*的總體趨勢一致，隨著第一層網(wǎng)格高度逐漸減小、網(wǎng)格數(shù)量逐漸增加，摩擦速度總體呈減小趨勢，當(dāng)?shù)谝粚泳W(wǎng)格高度h減小到0.005和0.002 5(即Case3和Case4)時計(jì)算結(jié)果基本一致，本文認(rèn)為Case3的網(wǎng)格達(dá)到網(wǎng)格無關(guān)性要求，在后續(xù)的積雪模擬中采用此套網(wǎng)格。

2.4 時間步長獨(dú)立性驗(yàn)證

本文采用Case3網(wǎng)格對時間步長獨(dú)立性進(jìn)行驗(yàn)證，根據(jù)特征長度與特征速度的比值，選取時間步長dt分別為0.02、0.01、0.005、0.002 5 s進(jìn)行時間步長獨(dú)立性驗(yàn)證。通過比較單位時間(h)內(nèi)積雪侵蝕沉積通量來分析時間步長對積雪模擬的影響，運(yùn)用EXECUTE_AT_END宏實(shí)現(xiàn)積雪通量的累加，并與穩(wěn)態(tài)計(jì)算的結(jié)果進(jìn)行了對比，其結(jié)果如圖3所示。時間步長的改變對積雪通量基本沒有影響，且和穩(wěn)態(tài)計(jì)算結(jié)果一致，而現(xiàn)實(shí)中的風(fēng)吹雪運(yùn)動持續(xù)時間可達(dá)數(shù)天，采用穩(wěn)態(tài)計(jì)算可以大大減少計(jì)算量，本文后續(xù)計(jì)算參考文獻(xiàn)[17,20]采用穩(wěn)態(tài)計(jì)算。

圖2 網(wǎng)格無關(guān)性檢驗(yàn)

圖3 時間步長獨(dú)立性檢驗(yàn)

2.5 平均風(fēng)壓系數(shù)比較

Silsoe Research Institute在空曠地區(qū)建造了邊長6 m的立方體，立方體表面的中線上布置測點(diǎn)，測量其壓力值和周圍的風(fēng)速值,用來研究鈍體周圍的壓力場和速度場。建筑物表面的風(fēng)壓分布與周圍的風(fēng)速場存在一定的關(guān)系，而建筑物近表面的風(fēng)速場主要決定了積雪遷移的結(jié)果。

本文采用k-kl-ω模型對立方體中心線上的風(fēng)壓進(jìn)行了模擬，并將模擬結(jié)果與Silsoe Research Institute實(shí)測結(jié)果和文獻(xiàn)[15] k-ε模型模擬結(jié)果進(jìn)行了比較，結(jié)果如圖4所示。從圖中可以看出：數(shù)值模擬平均風(fēng)壓系數(shù)與實(shí)測值趨勢一致，本文k-kl-ω模型計(jì)算結(jié)果相比文獻(xiàn)[15]結(jié)果在立方體迎風(fēng)面和頂面相對較好，在立方體背風(fēng)面兩數(shù)值模擬結(jié)果較為一致且數(shù)值上都略大于實(shí)測結(jié)果，由此可知k-kl-ω模型對風(fēng)場和平均風(fēng)壓系數(shù)的模擬要優(yōu)于k-ε模型，本文采用k-kl-ω模型模擬風(fēng)致雪漂移。

圖4 平均風(fēng)壓系數(shù)

2.6 立方體周邊積雪分布

2.6.1 立方體周圍風(fēng)場與積雪分布

立方體周圍風(fēng)場由簡化圖5表示，在立方體迎風(fēng)面由于立方體的阻擋，在A區(qū)形成駐渦，駐渦區(qū)風(fēng)速大，導(dǎo)致靠近立方體附近地面的積雪被風(fēng)吹走發(fā)生侵蝕；立方體背風(fēng)面B區(qū)，氣流繞過立方體在背風(fēng)面形成回流，但由于該旋渦風(fēng)速較小，因此積雪侵蝕很小；在立方體橫風(fēng)向的C區(qū)和D區(qū)，由于來流沿立方體兩側(cè)迎風(fēng)角分離形成馬蹄渦，速度很大發(fā)生較大侵蝕。

(a) 立方體迎風(fēng)截面流線圖

(b) 流場簡化圖

Oikawa實(shí)測、Tominaga和本文給出的立方體周圍無量綱水平積雪分布如圖6所示。Tominaga基于其以前的研究提出了一種新的雪漂模型，并給出了積雪水平分布圖6(b)(顏色越深，積雪厚度越大)。本文與Oikawa 和Tominaga都在兩側(cè)迎風(fēng)角附近發(fā)生較大侵蝕，但Tominaga結(jié)果侵蝕范圍相對較大，而本文模擬結(jié)果與Oikawa實(shí)測值僅在迎風(fēng)角附近小范圍內(nèi)雪深為0，Oikawa實(shí)測則出現(xiàn)了兩迎風(fēng)角侵蝕范圍不對稱現(xiàn)象，主要是實(shí)測過程中風(fēng)速大小和風(fēng)向發(fā)生變化導(dǎo)致的。在立方體迎風(fēng)區(qū)域由于Tominaga的模型未考慮空氣相與雪相的相互作用，未發(fā)生積雪侵蝕與實(shí)測結(jié)果不符合；橫風(fēng)向靠近立方體附近Tominaga的模擬過高的估計(jì)了積雪沉積，積雪深度達(dá)到2左右，大于實(shí)測結(jié)果，這主要是由于其模型過高的估計(jì)了非平衡躍移層的雪粒傳輸率，導(dǎo)致地面附近雪濃度過大所致，相對而言本文結(jié)果與實(shí)測值較為一致。

(a) Oikawa(2003)實(shí)測結(jié)果

(b) Tominaga(2011)模擬結(jié)果

(c) 本文模擬結(jié)果

2.6.2 積雪分布曲線比較

積雪分布曲線如圖7所示，圖中橫坐標(biāo)為采用立方體邊長H量綱一化后距立方體中心的距離，0表示立方體中心，-0.5和0.5表示立方體表面距中心的距離；縱坐標(biāo)為采用初始積雪深度hs(不受建筑影響區(qū)域的平均積雪深度)量綱一化后的積雪深度，大于1代表積雪沉積，小于1代表積雪侵蝕。

圖7 立方體周圍積雪分布

迎風(fēng)向截面，本文模擬結(jié)果與實(shí)測結(jié)果趨勢基本一致，迎風(fēng)面由于駐渦的影響靠近立方體附近發(fā)生侵蝕，最大侵蝕位置發(fā)生在0.75H左右，相比武岳，本文結(jié)果在-1.0H至1.0H之間與實(shí)驗(yàn)值更加吻合；實(shí)測值在-1.2H處出現(xiàn)沉積，而本文在-1.5H處出現(xiàn)沉積，這主要由于k-kl-ω模型過高的估計(jì)了湍流對雪顆粒的侵蝕作用；而Tominaga的模擬在迎風(fēng)側(cè)未出現(xiàn)侵蝕，且沉積量大于實(shí)測值。在背風(fēng)面受到背風(fēng)旋渦的影響發(fā)生侵蝕，除緊靠立方體地面，本文結(jié)果與武岳結(jié)果基本一致，最大侵蝕發(fā)生在1.2H處，Tominaga結(jié)果侵蝕較小，實(shí)測最大侵蝕位置在1.5H處，這主要由于數(shù)值模擬未考慮空氣中雪顆粒撞擊雪面引起的侵蝕和實(shí)測過程中風(fēng)向和風(fēng)速的改變導(dǎo)致的。

橫風(fēng)向截面，由于立方體的阻擋兩側(cè)形成馬蹄渦風(fēng)速很大，從而橫風(fēng)向截面主要發(fā)生侵蝕，靠近立方體處有少許沉積，本文結(jié)果與武岳較為一致，最大侵蝕位置發(fā)生在-0.9H和0.9H處；Tominaga的模擬結(jié)果侵蝕較大，在-0.75H至-1.35H之間積雪深度為0；實(shí)測結(jié)果C區(qū)發(fā)生-1.1H處，D區(qū)發(fā)生在1.0H處，由于實(shí)測過程中風(fēng)向的改變，兩側(cè)出現(xiàn)了不對稱的形狀。

整體而言，本文模擬結(jié)果趨勢基本與實(shí)測結(jié)果一致，在靠近立方體附近與實(shí)測結(jié)果吻合的很好，但由于數(shù)值模擬簡化了顆粒之間相互撞擊的影響和未考慮顆粒對雪面撞擊引起的侵蝕，在遠(yuǎn)離立方體區(qū)域模擬結(jié)果有些不足，而研究主要關(guān)注的也是靠近建筑附近的積雪情況，因此采用Mixture模型進(jìn)行風(fēng)吹雪的模擬在一定程度上是可行的。

3 高低屋蓋積雪模擬結(jié)果

高低屋蓋是我國屋蓋中比較常見的一種形式，一般由于高屋蓋對低屋蓋的阻擋作用，低屋蓋容易形成不均勻積雪，局部荷載過大導(dǎo)致屋蓋破壞。

圖8 高低屋蓋模型

Tsuchiya等在平坦的地貌對高低屋蓋模型積雪分布進(jìn)行了實(shí)測，并對其縮尺模型進(jìn)行了風(fēng)洞實(shí)驗(yàn)，其原型迎風(fēng)向截面尺寸如圖8所示，其寬度為6H，其中H等于0.9 m。

3.1 幾何建模與網(wǎng)格劃分

本文采用Tsuchiya實(shí)測模型模擬低屋蓋迎風(fēng)時積雪分布情況，H取0.9 m，計(jì)算域長、寬、高分別為100 m、43.2 m、15 m。計(jì)算域采用結(jié)構(gòu)化網(wǎng)格進(jìn)行離散，邊界層最小網(wǎng)格尺寸為0.005 m，網(wǎng)格總數(shù)為260萬，網(wǎng)格劃分如圖9所示。入流邊界參考高度取低屋蓋高度H=0.9 m，風(fēng)速取5 m/s，其余邊界與2.2節(jié)描述一致。

3.2 計(jì)算結(jié)果分析與比較

本文模擬了低屋蓋迎風(fēng)時積雪分布，由圖10可見,氣流由于低屋蓋的阻擋在地面形成回流，同時在低屋面前緣出現(xiàn)強(qiáng)分離流動, 較大的回流旋渦在低屋面上方形成，在建筑物后方形成較大回流旋渦，建筑物后面的地面拐角處存在較小的回流旋渦。

圖11給出了低屋蓋積雪分布三維視圖，在迎風(fēng)向0.5H處由于較強(qiáng)回流漩渦，積雪被風(fēng)吹走發(fā)生較大侵蝕；由于高屋蓋的阻擋，氣流沿橫風(fēng)向兩側(cè)繞流，在兩側(cè)角部位置積雪全部吹走，積雪深度為0；橫風(fēng)向積雪中間位置積雪較多，向兩側(cè)減少。

圖9 高低屋蓋網(wǎng)格劃分

圖10 高低屋蓋迎風(fēng)截面流線圖

圖11 低屋蓋積雪分布圖

本文將低屋面中線上的積雪深度與劉慶寬風(fēng)洞實(shí)驗(yàn)(采用食用鹽的風(fēng)洞實(shí)驗(yàn)結(jié)果)、Tsuchiya實(shí)測和王衛(wèi)華數(shù)值模擬等結(jié)果進(jìn)行了比較(圖12)，圖中橫坐標(biāo)是采用低屋蓋高度H量綱一體化，0表示高低屋蓋交界處，3代表低屋蓋邊緣處，縱坐標(biāo)為采用初始積雪深度ds(不受干擾的地面平均積雪深度)量綱一化后的積雪深度。在靠近高屋蓋附近0H至0.5H之間，本文結(jié)果與實(shí)測值吻合較好，兩者最大侵蝕都發(fā)生在0.35H左右；0.5H至1.2H之間數(shù)值模擬結(jié)果有些不足，主要原因是數(shù)值模擬無法考慮積雪休止角的影響導(dǎo)致；1.2H至3H本文結(jié)果介于實(shí)測和風(fēng)洞試驗(yàn)之間。

4 結(jié) 論

圖12 低屋蓋中線積雪分布

本文基于Euler-Euler方法，對計(jì)算流體力學(xué)軟件FLUENT17.0進(jìn)行了二次開發(fā)，采用Mixture多相流模型和湍流模型k-kl-ω三方程模型對立方體周邊和高低屋蓋表面的風(fēng)致積雪漂移進(jìn)行了模擬，積雪分布形式與實(shí)測和風(fēng)洞結(jié)果吻合較好。k-kl-ω三方程模型相比k-ε模型，對立方體中線的平均風(fēng)壓系數(shù)模擬更為準(zhǔn)確，其模型中引入了層流動能來考慮轉(zhuǎn)捩過程的影響。開展了時間步長獨(dú)立性驗(yàn)證并與穩(wěn)態(tài)計(jì)算結(jié)果進(jìn)行了對比，結(jié)果顯示瞬態(tài)與穩(wěn)態(tài)計(jì)算的單位時間積雪侵、蝕沉積通量基本一致，風(fēng)吹雪運(yùn)動往往持續(xù)時間達(dá)數(shù)天，采用穩(wěn)態(tài)計(jì)算可以大大減小計(jì)算量并取得合理結(jié)果。

Mixture多相流模型求解共同的連續(xù)方程和動量方程，在一定程度上考慮了空氣與雪相的相互作用，通過與實(shí)測和風(fēng)洞結(jié)果的對比，數(shù)值模擬結(jié)果積雪侵蝕、沉積趨勢一致。目前Euler-Euler方法無法考慮顆粒之間相互撞擊對積雪的影響和顆粒對雪面撞擊引起的積雪侵蝕，并且忽略了積雪休止角、雪粒升華等因素的影響，導(dǎo)致數(shù)值模擬結(jié)果在遠(yuǎn)離立方體處誤差較大，在今后的研究中應(yīng)進(jìn)一步完善。