典型形式大跨度屋蓋風(fēng)雪漂移的數(shù)值模擬

孫曉穎， 洪財濱， 武 岳

（1．哈爾濱工業(yè)大學(xué) 結(jié)構(gòu)工程災(zāi)變與控制教育部重點實驗室，哈爾濱 150090；2．哈爾濱工業(yè)大學(xué) 土木工程學(xué)院，哈爾濱 150090；3．中國航天建設(shè)集團有限公司，北京 100071）

近年來，因風(fēng)雪運動導(dǎo)致的工程事故接連不斷，造成大量的經(jīng)濟財產(chǎn)損失和人員傷亡［1］。表1列出了近幾年國內(nèi)外一些雪致災(zāi)害的工程事故。

這些工程事故中，大跨度屋蓋結(jié)構(gòu)是主要的結(jié)構(gòu)形式之一。導(dǎo)致這些工程事故的原因，很大程度上是因為積雪的局部堆積造成的，而積雪的局部堆積主要是由于風(fēng)致雪漂移所引起的。從1954年第一本荷載規(guī)范《荷載暫行規(guī)范》到2002年頒布的現(xiàn)行規(guī)范《建筑結(jié)構(gòu)荷載規(guī)范》（GB20009－2001），我國荷載規(guī)范中關(guān)于雪荷載的部分一共進行了四次修訂［2］。但歷次的修訂都把重點放在基本雪壓這方面，而對于雪荷載的分布系數(shù)則只是給出了幾種簡單屋蓋形式的雪荷載分布系數(shù)，并不適合大量涌現(xiàn)的大跨度屋蓋結(jié)構(gòu)。因此，積極開展大跨度屋蓋結(jié)構(gòu)風(fēng)雪運動的相關(guān)研究，對雪災(zāi)預(yù)測、防御具有重要的意義。

表1 國內(nèi)外雪致災(zāi)害Tab．1 Some domestic and foreign snow disaster

目前風(fēng)致雪漂移的研究方法主要有實地觀測、試驗方法及數(shù)值模擬三種方式［1］。實地觀測可獲得最準確的資料，但是受自然條件影響較大。風(fēng)洞試驗受相似準則的約束大，且周期長，成本高。CFD數(shù)值模擬技術(shù)近年來被應(yīng)用于研究建筑物周圍的風(fēng)雪運動。Alhajraf［2］基于兩相流理論，建立了防雪欄周圍風(fēng)雪漂移的數(shù)值模型，并與實測和風(fēng)洞試驗進行比較。Ue-matsu［3］對防雪欄、風(fēng)穴和小山周圍的風(fēng)雪漂移進行數(shù)值模擬研究，考慮了風(fēng)雪漂移的躍移運動和懸移運動，其模擬結(jié)果與實地觀測較為接近。

本文基于CFD數(shù)值模擬軟件Fluent，進行了二次開發(fā)，采用Mixture多相流模型和k－kl－ω湍流模型，對幾種典型形式大跨度屋蓋的風(fēng)雪漂移進行了數(shù)值模擬，并對積雪分布系數(shù)進行了探討。因文章篇幅有限，關(guān)于所用方法及參數(shù)設(shè)置的驗證部分未在文章中給出，詳見文獻［4］。

1 計算方法

本文基于兩相流理論，認為空氣相和雪相均為連續(xù)相［3］，采用Mixture多相流模型模擬風(fēng)雪流，利用有限體積法來求解風(fēng)雪流的控制方程。

1.1 風(fēng)雪流的控制方程

在Mixture模型中，空氣相和雪相共用一套連續(xù)方程和動量方程［5］：

（1）連續(xù)方程：

式中k＝1時表示空氣相，k＝2時表示雪相，以下類似，不再重復(fù)敘述。

（2）動量方程（N－S方程）：

其中，p為流場壓力；μm為空氣相和雪相的混合動力粘度；vdr，k＝vk－vm為相對漂移速度。

目前，工程上常采用對N－S方程平均的方法來計算流體的湍流運動。將N－S方程逐項平均可得基于雷諾應(yīng)力的時均N－S方程，但是方程組會多出一組新的未知變量—雷諾應(yīng)力方程不再封閉，需要引入湍流模型進行求解［6］。

（3）湍流方程

湍流模型選用k－kl－ω模型。k－kl－ω模型是一種三方程的粘渦模型，包含湍動能kT的輸運方程，層流動能kL的輸運方程和逆湍流時間尺度w的輸運方程［5］：

其中，DT、DL為擴散項；PKT、PKL為湍流、層流動能產(chǎn)生項；Cw1、Cw2、Cw3為常數(shù)項；

R為考慮氣流由層流到湍流內(nèi)波動的影響。

k－kl－ω模型能有效描述邊界層中由層流轉(zhuǎn)變到湍流的轉(zhuǎn)捩，而風(fēng)致雪漂移的運動過程中，從躍移運動到懸移運動的轉(zhuǎn)變，類似層流到湍流的轉(zhuǎn)變，因而該模型能更好地模擬近地面的風(fēng)雪流。

（4）雪相濃度方程［5］

其中，mpq為雪相（p相）轉(zhuǎn)化成空氣相（q相）的質(zhì)量；mqp為空氣相（q相）轉(zhuǎn)化成雪相（p相）的質(zhì)量。由于不考慮雪的升華，因此mpq和 mqp的取值都為零。

對控制微分方程進行離散時，采用有限體積法將微分方程分解成一系列關(guān)于多個變量的非線性耦合代數(shù)方程組。在運用有限體積法時，采用二階迎風(fēng)格式離散對流項，對擴散項采用具有二階精度的中心差分格式；對于壓力—速度耦聯(lián)方程，采用SIMPLE算法實現(xiàn)各聯(lián)立方程的解耦及壓力場和速度場的校正。

1.2 積雪的侵蝕與沉積

壁面摩擦速度是風(fēng)致雪漂移運動中非常重要的參數(shù)。壁面摩擦速度與壁面剪切應(yīng)力直接相關(guān)，當(dāng)壁面剪切應(yīng)力克服雪顆粒內(nèi)部作用力和重力時，雪顆粒便開始運動。壁面摩擦速度的表達式為：

亦可用對數(shù)律的風(fēng)速剖面公式近似計算壁面摩擦速度：

式中，ρ為空氣密度；u（z）為z高度處的風(fēng)速；κ為馮卡門常數(shù)，取 0．4；z0為壁面粗糙度高度，取 3．0×10－5［7］。

當(dāng)摩擦速度u*超過閾值速度u*t時，雪顆粒在風(fēng)力的作用下離開積雪面，積雪面發(fā)生侵蝕現(xiàn)象；當(dāng)摩擦速度u*小于閾值速度u*t時，漂移中的雪顆粒便堆積在積雪面，發(fā)生沉積現(xiàn)象。侵蝕量qero和沉積量qdep的表達式分別為［8］：

其中，Aero為常系數(shù)，取7．0×10－4；wf為雪的沉降速度；u*t為閾值速度；Фs為雪的質(zhì)量濃度（Фs＝ρsf，f為雪的體積分數(shù)，ρs為雪顆粒密度）。

受風(fēng)的影響，積雪面高度會時刻變化，這必然會使計算域的邊界產(chǎn)生變化，從而影響計算結(jié)果。而李雪峰［9］在對二維臺階式屋蓋模型采用動態(tài)網(wǎng)格技術(shù)模擬積雪面高度變化的結(jié)論中指出，積雪面變化的影響不到0．1%。另外，三維動態(tài)網(wǎng)格對計算機性能的要求非常高，因此本文未考慮積雪面變化對計算域的影響。

2 典型屋蓋的數(shù)值模擬

2.1 計算參數(shù)設(shè)置

（1）計算域及網(wǎng)格劃分

計算域的設(shè)置如圖1所示，其中計算域各參數(shù)設(shè)置如表2所示。各屋蓋模型的網(wǎng)格劃分如圖所示。

圖1 計算域示意圖Fig．1 Size of numerical simulation field

表2 計算域各參數(shù)設(shè)置Tab．2 Com putational field setting

表3 網(wǎng)格劃分Tab．3 Grid meshing

（2）邊界條件

空氣相邊界條件：采用速度入流的邊界條件，入口風(fēng)速以 v（z）＝v0（z／z0）α的形式給出，其中 α取值為0．15，參考高度 z0取10 m，參考高度處風(fēng)速 v0取為5 m／s。

來流湍流特性通過直接給定湍流動能k和湍流耗散率ε值的方式來定義：

式中：l是湍流特征尺度；I為湍流強度，根據(jù)《建筑結(jié)構(gòu)荷載規(guī)范》（GB50009－2012）進行取值［13］：

式中：I10＝0．14，α＝0．15。

雪相邊界條件：受空氣相的影響，雪相入口的速度與空氣相一致。而雪相的入口濃度在躍移層和懸移層內(nèi)則不同［11－12］：

式中：g為重力加速度；up為雪相運動速度；hs為躍移層和懸移層的臨界高度。

出口邊界條件：采用壓力出口的邊界條件。

流體域頂部和兩側(cè)：自由滑移的壁面條件。

建筑物表面和地面：無滑移的壁面條件。

（3）計算結(jié)果的表示形式

根據(jù)《建筑結(jié)構(gòu)荷載規(guī)范》（GB50009－2001）［14］中的規(guī)定，屋面水平投影面上雪荷載標準值的計算公式為：

式中：Sk為雪荷載標準值；μr為屋蓋積雪分布系數(shù)；S0為基本雪壓，可由規(guī)范查得。

CFD計算結(jié)果將以積雪分布系數(shù)μr表示，μr＞1．0表示積雪沉積；μr＜1．0表示積雪被侵蝕。

2.2 數(shù)值模擬結(jié)果

2．2．1 平屋蓋

平屋蓋的計算模型如圖2所示，屋蓋跨度L取80 m，屋蓋離地面高度H取20 m，高跨比H／L為1／4。

圖2 計算模型Fig．2 Sketch ofmodel

圖3 不同風(fēng)向角下平屋蓋的積雪分布系數(shù)Fig．3 Snow distribution factor on flat roofs in differentwind directions

圖3給出了不同風(fēng)向角下平屋蓋積雪分布系數(shù)云圖，可以看出：

（1）屋蓋表面的積雪分布以侵蝕為主，靠近屋蓋前沿的部位侵蝕較大，隨著距離的增大，侵蝕逐漸減小，局部出現(xiàn)沉積，最大沉積發(fā)生在距迎風(fēng)面前沿約L／3處。在0°風(fēng)下，風(fēng)雪流垂直于平屋蓋外邊緣，來流在屋面前沿分離形成柱狀渦，在漩渦區(qū)風(fēng)速大，故積雪侵蝕大；當(dāng)風(fēng)雪流以某一斜向角度（45°）經(jīng)過建筑物時，屋蓋表面的積雪受錐形渦影響的區(qū)域侵蝕量大，隨著距離的增大侵蝕也逐漸減小。總體看來，平屋蓋迎風(fēng)側(cè)屋面的積雪受柱狀渦或錐形渦的作用發(fā)生較大的侵蝕，局部有少量沉積，在其他區(qū)域由于尾流作用積雪高度變化不大。

（2）風(fēng)向角對積雪分布的影響較大。不同風(fēng)向下，由于來流和漩渦作用均有較大的不同，使得屋蓋表面的風(fēng)速場產(chǎn)生較大的差異，因而積雪的侵蝕和沉積最大值出現(xiàn)的位置也不同。

從數(shù)值模擬結(jié)果可以看出，積雪在風(fēng)力作用下呈現(xiàn)明顯的不均勻分布，并且積雪的分布以侵蝕為主，大部分區(qū)域的積雪分布系數(shù)小于1．0，但在局部地方有積雪堆積，最大積雪分布系數(shù)為1．27。而我國規(guī)范對于平屋蓋的積雪分布系數(shù)給定的值為1．0，即積雪是均勻分布的，并沒有平屋蓋的局部堆雪效應(yīng)。

2．2．2 球面屋蓋

球殼是的計算模型如圖4所示，四周封閉，跨度L為60 m，高度H為12m，矢跨比 f／L分別為1／5、1／8和1／10。

圖4 計算模型Fig．4 Sketch ofmodel

圖5給出了球面屋蓋積雪分布系數(shù)云圖，可以看出：

（1）迎風(fēng)側(cè)屋蓋上的積雪出現(xiàn)大面積的侵蝕，且侵蝕最大的地方主要分布在屋蓋的拐角區(qū)域，且拐角區(qū)域積雪高度的變化梯度較大，而在背風(fēng)區(qū)域出現(xiàn)沉積。從圖6可以看出，迎風(fēng)側(cè)屋蓋上風(fēng)雪流速較大，積雪在風(fēng)力作用下發(fā)生大面積的侵蝕；風(fēng)雪流在屋蓋拐角處發(fā)生分離，使得背風(fēng)區(qū)域風(fēng)速較小，故背風(fēng)側(cè)積雪出現(xiàn)沉積現(xiàn)象。

（2）矢跨比是影響積雪分布的重要因素之一。隨著矢跨比的增加，迎風(fēng)側(cè)積雪的侵蝕量和侵蝕范圍逐漸減小，背風(fēng)側(cè)的積雪沉積量和沉積范圍則逐漸增加。

圖5 不同矢跨比球面屋蓋在不同風(fēng)向角下的積雪分布系數(shù)Fig．5 Snow distribution factor on spherical dome with different rise-span ratio

圖7 計算模型Fig．7 Sketch ofmodel

2．2．3 鞍形屋蓋

鞍形屋蓋的計算模型如圖7所示，建筑結(jié)構(gòu)四周封閉，跨度L為80 m，高度H為10 m，矢跨比f／L分別為 1／8、1／12和 1／16。

圖8給出了鞍形屋蓋積雪分布系數(shù)云圖，可以看出：

（1）在屋蓋的邊沿和拐角區(qū)域，屋蓋上的積雪發(fā)生較大的侵蝕；局部有積雪沉積。0°風(fēng)向角和90°風(fēng)向角下，受風(fēng)雪流錐形渦影響區(qū)域氣流流速大，積雪的侵蝕亦較大；45°風(fēng)向角下，風(fēng)雪流垂直于平屋蓋外邊緣，在橫風(fēng)向屋面邊沿（LH′邊）出現(xiàn)了最大的侵蝕，而在風(fēng)雪流分離區(qū)和背風(fēng)側(cè)則出現(xiàn)局部堆積。

（2）風(fēng)向角對積雪分布的影響較大。不同風(fēng)向角下，屋蓋表面的漩渦作用非常不同，屋蓋表面的風(fēng)速場也有較大的差異，因而積雪的分布規(guī)律也有不同的表現(xiàn)。

（3）鞍形屋蓋的積雪分布形式既體現(xiàn)了與平屋蓋類似的分布形式，又隨著屋蓋曲率的變化體現(xiàn)出差異性，這表明積雪分布受矢跨比的影響較大。在0°風(fēng)向角下，屋蓋順風(fēng)向呈現(xiàn)向下凹的形式，隨著矢跨比的增加，在屋蓋頂點附近的積雪侵蝕量先增加后減少，其他區(qū)域的侵蝕量則逐漸減少；45°風(fēng)向角下，隨著矢跨比的增加，積雪沉積量先減少后增加，而侵蝕范圍則不斷擴大；90°風(fēng)向角下，屋蓋拱起，隨著矢跨比的增加，侵蝕量減少，且侵蝕范圍也逐漸減小。

圖8 不同矢跨比鞍形屋面在不同風(fēng)向角下的積雪分布系數(shù)Fig．8 Snow distribution factor on shaddle－shaped roofs with different rise-span ratio under differentwind directions

2．2．4 拱形屋蓋

如圖9所示，拱形屋面計算模型為一個菱形平面拱形屋蓋結(jié)構(gòu)，四周封閉，跨度L為30 m，高度H為12 m，矢跨比 f／L分別為1／3、1／5和 1／10。

圖9 計算模型Fig．9 Sketch ofmodel

圖10給出了拱形屋蓋的積雪分布系數(shù)云圖，可以看出：

（1）迎風(fēng)側(cè)屋蓋上的積雪出現(xiàn)大面積的侵蝕，且侵蝕最大的地方主要分布在屋蓋的拐角區(qū)域。屋蓋拐角區(qū)域積雪高度的變化梯度較大，并在背風(fēng)區(qū)域出現(xiàn)沉積。0°風(fēng)向角下，積雪分布與球殼屋蓋類似，迎風(fēng)側(cè)屋蓋上風(fēng)雪流速較大，積雪在風(fēng)力作用下大面積侵蝕，風(fēng)雪流在屋蓋拐角處發(fā)生分離，使得背風(fēng)區(qū)域風(fēng)速較小，故背風(fēng)側(cè)積雪出現(xiàn)沉積現(xiàn)象；90°風(fēng)向角下，積雪分布則與平屋蓋類似，局部出現(xiàn)堆積。

（2）風(fēng)向角對積雪分布的影響較大。不同風(fēng)向角下，屋蓋表面的漩渦作用非常不同，屋蓋表面的風(fēng)速場也有較大的差異，因而積雪的分布規(guī)律也有不同的表現(xiàn)。

（3）矢跨比是影響積雪分布的重要因素之一。隨著矢跨比的增加，積雪的侵蝕量和侵蝕范圍逐漸減小，沉積量和沉積范圍則逐漸增加。

表4 規(guī)范中拱形屋蓋的積雪分布系數(shù)劃分Tab．4 The snow distribution factor of cylinder roof in Chinese Code、Eurocode and NBCC

圖10 不同矢跨比拱形屋蓋在不同風(fēng)向角下的積雪分布系數(shù)Fig．10 Snow distribution factor on cylinder roofs with different rise-span ratio under differentwind directions

圖11 數(shù)值模擬結(jié)果與規(guī)范對比（矢跨比1／5）Fig．11 Comparison of the simulated results and codes

對于矢跨比為1／5的拱形屋蓋，按中國規(guī)范，積雪呈均勻分布，積雪分布系數(shù)為5／8，而從數(shù)值模擬結(jié)果可以看出，積雪在不同風(fēng)向下呈現(xiàn)明顯的不均勻分布，在多個區(qū)域里出現(xiàn)的分布系數(shù)都大于5／8，且5／8小于歐洲規(guī)范［15］的 0．8和 NBCC［16－17］的 1．0，偏于不安全。圖11給出了數(shù)值模擬結(jié)果與歐洲規(guī)范與NBCC中不均勻分布情況的對比圖。可以看出：各國規(guī)范積雪不均勻分布情況不盡相同，數(shù)值模擬結(jié)果積雪分布趨勢與加拿大規(guī)范相似。在數(shù)值上，數(shù)值模擬結(jié)果比各國規(guī)范小，但在屋檐處則較各國規(guī)范大，這主要是由于初始假設(shè)屋檐處能堆積的積雪厚度與其他區(qū)域相同。盡管數(shù)值模擬結(jié)果與各國規(guī)范略有出入，但不可否認積雪確實存在不均勻分布的情況，而我國規(guī)范并未考慮這點。

3 結(jié) 論

（1）大跨度屋蓋表面的積雪分布以侵蝕為主，局部出現(xiàn)積雪沉積。風(fēng)雪流在屋蓋前沿及拐角處產(chǎn)生的漩渦脫落（分離和再附）是影響屋面積雪分布規(guī)律的主要因素。

（2）屋蓋表面的積雪分布受風(fēng)向角的影響較大。風(fēng)向角不同，來流和漩渦作用均有較大的不同，使得屋蓋表面的風(fēng)速場產(chǎn)生較大的差異。

（3）矢跨比是影響積雪分布的重要因素之一。矢跨比直接影響著屋蓋的形狀，進而影響屋蓋表面的風(fēng)速場，使不同矢跨比的屋蓋上的積雪形成不同的分布形式。。

（4）通過各國規(guī)范、數(shù)值模擬結(jié)果與我國荷載規(guī)范對比發(fā)現(xiàn)：① 風(fēng)致雪漂移將使平屋蓋產(chǎn)生局部堆雪效應(yīng)，而我國荷載規(guī)范并未考慮這點；② 我國規(guī)范中拱形屋面均勻分布的積雪分布系數(shù)取值較歐洲規(guī)范和加拿大規(guī)范都偏小，且未考慮積雪的不均勻分布形式。

［1］范峰．雪致建筑結(jié)構(gòu)災(zāi)害調(diào)查及國內(nèi)外雪荷載規(guī)范對比［R］．2010，杭州．

［2］莫華美．典型屋面積雪分布的數(shù)值模擬與實測研究［D］．哈爾濱：哈爾濱工業(yè)大學(xué)，2011年．

［3］Bintanja R．Snowdrift suspension and atmospheric turbulence．Part I：Theoretical background and model description［J］．Springer Netherlands，2000，95（3）：343－368．

［4］洪財濱．大跨度屋蓋風(fēng)致雪漂移的數(shù)值模擬研究［D］．哈爾濱：哈爾濱工業(yè)大學(xué)，2012．

［5］Ansys Inc．ANSYSFLUENT 14．0：Theory Guide［Z］．2011．

［6］王福軍．計算流體動力學(xué)分析－CFD軟件原理與應(yīng)用［M］．北京：清華大學(xué)出版社，2004．

［7］張相庭．工程抗風(fēng)設(shè)計計算手冊［M］．北京：中國建筑工業(yè)出版社，1998．

［8］Beyers J H M，Sundsbo P A，Harms T M．Numerical simulation of three-dimensional， transient snow drifting around a cube［J］．Journal of Wind Engineering and Industrial Aerodynamics，2004，92（9）：725－747．

［9］李雪峰．風(fēng)致建筑屋蓋表面及其周邊積雪分布研究［D］．上海：同濟大學(xué)，2011．

［10］AIJ．Recommendations for loads on Buildings（English Version）［S］．2004．

［11］ Pomeroy JW，Gray D M．Saltation of snow［J］．Water Resources Research，1990，26（7）：1583－1594．

［12］Pomeroy JW，Male D H．Steady-state suspension of snow［J］．Journal of Hydrology，1992，136（1－4）：275－301．

［13］中華人民共和國建設(shè)部．GB 50009－2012建筑結(jié)構(gòu)荷載規(guī)范［S］．北京：中國建筑工業(yè)出版社，2012．

［14］中華人民共和國建設(shè)部．GB 50009－2001建筑結(jié)構(gòu)荷載規(guī)范［S］．北京：中國建筑工業(yè)出版社，2001．

［15］BSI．BS EN1991－1－3：2003 Eurocode 1－Actions on structures－Part 1－3：General actions－Snow loads［S］．London：2003．

［16］NBCC．NBCC 2005National Building Code of Canada 2005［S］．Ottawa：2005．

［17］ NBCC．User’s Guide－NBC 2005 Structural Commentaries（Part4 of Division B）［M］．Ottawa：2005．