平屋面摩擦速度的數值模擬研究

辛林桂，周晅毅，顧 明

(同濟大學土木工程防災國家重點實驗室，上海 200092)

平屋面是一種典型的屋蓋形式，其跨度對建筑周圍的流場有著重要的影響[1-5]。風經過平屋面時會出現復雜的繞流、再附現象，在氣流作用下雪顆粒將發生復雜的飄移運動，從而造成屋面和建筑周圍積雪的不均勻分布[6-9]。大跨度的平屋面結構對于不均勻分布的雪荷載很敏感，雪荷載的不均勻分布形式往往是結構的最不利荷載分布之一。因此，準確預測不同跨度屋面上的積雪分布和傳輸特征有利于保障結構的安全。

目前對平屋面積雪分布和傳輸的研究方法主要有現場觀測[10-11]、風洞或水槽試驗[12-13]和數值模擬[14-15]三種?，F場觀測記錄了平屋面在自然環境中的積雪分布和傳輸規律，反映了最為真實積雪飄移的運動。THIIS 和GJESSING[10]在挪威對平屋面進行了屋面積雪的實測研究，結果表明屋面尺寸對屋面積雪的飄移和分布有顯著的影響；ZHANG 等[11]采用現場觀測研究了不同比例的平屋面和其他屋面模型的縮尺比對屋面積雪分布的影響，觀測結果表明幾何縮尺比對屋面積雪分布影響甚微。與現場測量相比，風洞或水槽試驗可以提供一個穩定可控的試驗平臺來研究屋面的積雪飄移。O’ROURKE 等[12]在水槽中利用核桃殼對平屋面積雪傳輸率進行了研究，結果表明，水槽模型與全比例模型的輸運率吻合較好；QIANG 等[13]在低溫風洞中使用人工雪顆粒對積雪在平屋頂上的飄移進行了模擬，研究了降雪期間與無降雪期間積雪飄移的異同。不同于現場實測和風洞或水槽試驗，對屋面積雪的數值模擬更偏向于機理研究。LIU 等[14]利用拉格朗日方法建立了一個數值模型研究了平屋面的積雪飄移，探討了擊濺夾帶對平屋頂積雪輸運的作用及其貢獻。QIANG 等[15]介紹了一種積雪傳輸模型，用于計算平屋面的風致積雪傳輸量，在該模型中考慮了降雪量和閾值摩擦速度等變量的差異，模擬結果與前人的實測結果較為一致。上述的研究采用不同方法分析了平屋面積雪的分布特征，但忽略了屋面跨度可能對積雪傳輸造成的影響；同時，O’ROURKE 等[16]通過研究發現飽和積雪傳輸率應按傳輸距離進行折減。由于目前對平屋面的積雪研究中，鮮見有學者考慮屋面跨度對屋面雪荷載的影響，因此結合屋面積雪傳輸率與傳輸距離分析跨度對平屋面雪荷載的影響是非常必要的。

積雪野外實測的研究結果表明，傳輸率與摩擦速度之間存在指數關系[17-18]。并且一些研究結果表明，建筑屋面的摩擦速度是決定其風致遷移雪荷載的重要因素[19]，如：康路陽等[20]采用分階段的歐拉-歐拉方法模擬了有休止角的積雪表面的摩擦速度，并結合摩擦速度預測了平屋面的積雪分布；CAO 等[21]通過風洞試驗和數值模擬研究了低坡度屋面的積雪質量傳輸率與摩擦速度之間的關系。因此，準確獲得平屋面的摩擦是研究屋面積雪傳輸的前提條件。

本文首先采用CFD 方法模擬了平屋面的摩擦速度，并且為了研究不同湍流模型的模擬精度，利用Irwin 探頭風洞試驗與數值模擬結果進行了對比。接著，根據CFD 模擬的結果分析了不同跨度平屋面流場和摩擦速度的分布特征。最后，基于數值模擬得到的摩擦速度，計算了平屋面積雪傳輸率，并討論了不同來流風速下屋面跨度對積雪傳輸率和侵蝕通量的影響。

1 數值模型

1.1 研究對象

為了獲得不同跨度平屋面的摩擦速度，本文采用RANS 方法對二維平屋面的周圍的流場進行模擬。研究對象的跨度L分別為30 m、60 m、90 m和120 m，高度H為15 m，屋面跨度與高度之比L/H分別為2、4、6 和8，如圖1 所示。

圖1 研究對象Fig.1 Research object

1.2 計算域及網格方案

CFD 數值模擬采用ANSYS FLUENT 軟件進行，數值模型的幾何縮尺比為1∶100，圖2 為計算域和網格劃分示意圖(未按比例畫出)，計算域的長度為16L，高度為20H。建筑迎風面到入口(velocity inlet)邊界的距離為5L，出口(outflow)邊界距離建筑背風面為10L， 如圖2(a)所示。以模型L/H=4 為例，高分辨率的結構化網格如圖2(b)所示，最小網格尺寸約為0.03H，網格的橫向和縱向的增長率均為1.05，總網格數量約為5.5 萬。

圖2 計算域和網格劃分示意圖Fig.2 Schematics of computational domain and computational grid

1.3 邊界條件和參數設置

計算域入口邊界采用速度入口，參考《建筑結構荷載規范》對B 類地貌的規定，并按照幾何縮尺比(1∶100)縮放后在數值模擬中進行應用。入口平均風剖面U(y)如下所示[22]：

式中：風剖面指數α 取0.15；為了探究風速對屋面摩擦速度和積雪傳輸率的影響，本文取U0.1(距地面0.1 m 高度處風速)的風速大小分別為6 m/s、8 m/s 和10 m/s。

湍流強度剖面I(y)為[22]：

根據風速剖面U(y)和湍流強度剖面I(y)，湍動能剖面k(y)和湍流耗散率ε(y)剖面為[23]：

式中，Cμ為0.09。

下邊界和建筑表面的壁面邊界條件采用無滑移標準壁面函數。上邊界設置為對稱邊界，出口邊界設置為0 梯度的出流邊界。采用穩態RANS方法，對Standardk-ε、 RNGk-ε 和Realizablek-ε三種湍流模型進行了比較分析。所有的傳輸方程均采用二階迎風離散格式，壓力速度耦合方法為SIMPLE 算法，收斂準則為10-6。

1.4 網格無關性分析

采用三種尺寸的網格進行網格敏感性分析，分別為細網格、中網格和粗網格。細網格的最小網格尺寸為0.01H，屋面y+的范圍為25～50；中網格的最小網格尺寸為0.03H，屋面y+的范圍為31～70；粗網格的最小網格尺寸為0.05H，屋面y+的范圍為38～95。圖3 為模型尺寸L/H=4 時，采用Realizablek-ε 模型和三種網格劃分方案模擬的屋面摩擦速度對比。從圖3 可以看出，粗網格與細網格的差異較大，而中網格與細網格的差異較小，因此本文在下文中選擇中網格的模擬結果進行分析。

圖3 網格無關性分析Fig.3 Grid-sensitivity analysis

2 湍流模型對模擬結果的影響

為了驗證Standardk-ε、 RNGk-ε 和Realizablek-ε 三種湍流模型的模擬精度，并對比不同湍流模型對模擬結果的影響：首先，對Irwin 探頭進行了標定；接著，在不同來流風速下采用標定后Irwin探頭測量了平屋面的摩擦速度；最后，將不同湍流模型的數值模擬和風洞試驗進行了對比。風洞試驗的風壓測量采用美國Scanvival 公司生產的ZOC33 電子壓力掃描模塊，壓力掃描閥的采樣頻率為312.5 Hz，采樣時間為30 s。

2.1 Irwin 探頭的標定

本文采用的Irwin 探頭為無方向探頭[24]，尺寸如圖4 所示，測壓孔K1到K2的垂直距離h為5 mm。

圖4 Irwin 探頭尺寸 /mmFig.4 Dimensions of the Irwin sensor

FERREIRA 等[25]和FARIA 等[26]認為壁面切應力與Irwin 探頭測得壓差存在指數關系，即：

式中：a、b為Irwin 探頭測量壁面切應力時的標定系數；ΔP=P0-Ph，為Irwin 探頭測得的壓差。由摩擦速度u*=， ρa為空氣密度，可得u*與ΔP的關系為：

式中，m、n為Irwin 探頭測量摩擦速度時的標定系數。標定試驗在同濟大學TJ-2 大氣邊界層風洞的空風場中，根據WANG 等[27]的方法對Irwin 探頭測量摩擦速度進行了標定，如圖5(a)所示。假設標定試驗的近壁面風速符合對數分布：

圖5 Irwin 探頭測量摩擦速度的標定示意圖Fig.5 Calibration diagram of friction velocity measured by the Irwin sensor

式中：U(z)為z/m 高度處的風速，標定試驗中不同風速下的風剖面U(z)分布如圖5(b)所示；Von Karman 常數κ為0.42；z0為標定平臺的氣動粗糙度高度。通過式(7)的關系對風剖面進行擬合后，可得到z0為1.173×10-5m，五種不同風速下的u*分別為0.17 m/s、0.23 m/s、0.32 m/s、0.41 m/s 和0.50 m/s。于是再根據式(6)對高度h=5 mm 的Irwin 探頭測得的 ΔP和u*進行標定，通過重復性和方向性檢驗后，可得到h=5 mm 的Irwin 探頭在測量摩擦速度時的標定系數m=0.183、n=0.338。

2.2 Irwin 探頭風洞試驗

采用Irwin 探頭在同濟大學TJ-2 大氣邊界層風洞對平屋面的摩擦速度進行測量，風洞的試驗段長15 m、寬3.0 m、高2.5 m，如圖6 所示。屋面模型的幾何縮尺比為1∶100，跨度為0.6 m，高度為0.15 m，阻塞比小于5%。在進行風洞試驗時僅模擬風場的二維效應，模型的兩側增加了側擋板，并且為了減小擋板對風場的影響，擋板端部進行了削尖處理。Irwin 探頭在風洞試驗模型上的布置如圖7 所示?？紤]到Irwin 探頭相互遮擋可能會影響測量精度，故在屋面上對Irwin 探頭進行了交錯布置，兩排Irwin 探頭的間距為50 mm，Irwin探頭距離迎風側和背風側屋檐的距離均為10 mm，其他的Irwin 探頭在屋面上均勻排布。

圖6 風洞試驗試驗段示意圖 /mFig.6 Schematic diagram of wind tunnel test section

圖7 風洞試驗模型測點布置 /mmFig.7 Layout of measuring points of wind tunnel test mode

在測量平屋面的摩擦速度時，定義屋面高度H為參考高度，平屋面Irwin 探頭風洞試驗采用了三種風速：參考風速U(H)分別為6 m/s、8 m/s 和10 m/s，參考高度處的湍流強度I(H)約為16%，風速和湍流強度剖面如圖8 所示。參考了康路陽等[20]對風洞試驗風速的擬合方法，風速剖面可由式(7)進行擬合，氣動粗糙度高度z0為5.362×10-5m，湍動能由式(2)進行計算，湍流耗散率的定義為：

圖8 風洞試驗風速和湍流強度剖面Fig.8 Wind speed and turbulence intensity profile of the wind tunnel test

式中，Cμ為經驗常數，取0.09。

2.3 不同湍流模型的對比

為了對比Standardk-ε、 RNGk-ε、Realizablek-ε 三種湍流模型對數值模擬結果的影響，將CFD 模擬獲得的屋面摩擦速度u*與風洞試驗結果進行了對比。以參考風速U(H)為8 m/s 為例，數值模擬與風洞試驗的對比如圖9 所示。

圖9 數值模擬與風洞試驗的摩擦速度對比(U(H)=8 m/s)Fig.9 Comparison of friction velocity between numerical simulation and wind tunnel test (U(H)=8 m/s)

在三種風速下，采用FAC2(the factor of 2)、線性相關系數R和平均百分比誤差APD 對數值模擬結果進行評價。FAC2、R和APD 的定義如下所示[28 - 30]：

式中：Pi和Mi(i=1, 2,…,N;N=48)分別為CFD 模擬的模擬值和風洞試驗測量值； σP和 σM分別為CFD 模擬和風洞試驗結果的標準差；上劃線代表平均值，R理想值為1.0，APD 理想值為0%。

根據上述三種評價指標，不同湍流模型的數值模擬結果評價如表1 所示。在FAC2 評價中，Realizablek-ε 模型的摩擦速度的評價結果與理想值最為接近，這表明Realizablek-ε 模型的整體表現最好，很少發生明顯過高或過低的預測。對比線性相關系數R可得，Realizablek-ε 模型模擬的摩擦速度與風洞試驗測量結果相比相關性最好。在平均百分比誤差APD 評價中，Realizablek-ε 模型的摩擦速度評價結果與理想值0%最接近。綜上所述，Realizablek-ε 模型對摩擦速度的預測表現最好。因此，在后文的研究中選擇Realizablek-ε模型的數值模擬結果進行分析。

表1 三種風速下不同湍流模型模擬結果對比Table 1 Comparisons of simulation results of different turbulence models under three wind speeds

3 結果討論

3.1 不同跨度平屋面的流場和摩擦速度分布

以U10=8 m/s 為例，不同跨度平屋面的風速分布如圖10 所示。由圖10 可得，來流在屋蓋迎風側的前緣發生分離，迎風前緣最大風速約為來流風速U10(8 m/s)的1.5 倍。來流在迎風前緣發生強烈的分離，形成旋渦后發生再附，并在屋面后方形成一個回流區。屋面跨度L=2H時(圖10(a))，分離產生的旋渦幾乎覆蓋了整個屋面；隨著屋面跨度的增大(圖10(b)～圖10(d))，旋渦的長度也在不斷增大，但旋渦的長度與屋面跨度的比值在不斷的縮小。同時，旋渦后方近屋面的風速也隨著跨度的增大不斷增大。

圖10 不同跨度平屋面周圍的風速分布Fig.10 Wind speed distribution around flat roofs with different spans

摩擦速度可以確定雪顆粒運動狀態，在數值模擬中屋面摩擦速度u*可由屋面剪切應力τw計算獲得，即u*=。受近屋面風速的影響，不同跨度屋面的摩擦速度分布如圖11 所示。由于來流在屋蓋前緣分離的影響，摩擦速度在屋蓋前緣處較大，并且不同跨度屋面的摩擦速度在屋蓋前緣處大小分布基本相等；而后經過分離產生的旋渦覆蓋區域，摩擦速度先增大后減??；最后旋渦覆蓋區域后方風速平穩發展，摩擦速度也逐漸增大。參考孫曉穎等[31]對雪顆粒閾值摩擦速度的取值，本文取雪顆粒的閾值摩擦速度u*t為0.15 m/s。當屋面跨度與高度之比L/H=2 時，屋面摩擦在背風側屋檐小于閾值，屋面其他位置均大于閾值，屋面整體呈現侵蝕狀態。隨著屋面跨度的增大，屋面摩擦速度的最小值所在的位置逐漸向迎風側屋檐靠近，即屋面上產生積雪沉積的位置也逐漸靠近迎風側屋檐。

圖11 不同跨度屋面的摩擦速度分布Fig.11 Friction velocity distribution of roofs with different spans

圖12 不同跨度的平屋面表面的平均摩擦速度Fig.12 Mean friction velocity of roofs with different spans

3.2 屋面跨度對積雪傳輸和侵蝕的影響

本節采用CFD 模擬獲得的平屋面摩擦速度的分布討論屋面跨度對屋面積雪傳輸的影響。屋面積雪傳輸率Qu*與數值模擬獲得的平均摩擦速度的關系可由OWEN[17]和KIND[18]的研究結果獲得：

根據屋面積雪傳輸率Qr，屋面積雪沿跨度的平均侵蝕通量qr可以表示為Qr與L比值：

圖13 為不同風速下跨度對平屋面的積雪傳輸率和侵蝕通量的影響。由于來流風速為屋面積雪的傳輸和侵蝕提供了動力，對比圖13(a)、圖13(b)可得，屋面跨度相同時，風速越大屋面積雪的侵蝕就越劇烈，屋面積雪傳輸率和侵蝕通量也就越大。圖13(a)顯示屋面的積雪傳輸率隨著屋面跨度的增大而增大。同時，從圖13(b)可見，屋面積雪的沿跨度平均侵蝕通量隨著屋面跨度的增大而減小。

圖13 不同跨度平屋面的積雪傳輸率和侵蝕通量Fig.13 Snow transport rate and erosion flux of flat roofs with different spans

3.3 屋面跨度對雪荷載暴露系數的影響

為了進一步分析屋面跨度對屋面雪荷載的影響，本文參考ASCE 荷載規范[34]對平屋面雪荷載Pf的定義：

式中：Ce為暴露系數；Ct為熱力系數；Is為重要系數；Pg為地面雪荷載，kg/m2。將Ct與Is取1.0 后，Ce可以表示為：

式中：Gg/(kg/(m2·s))為地面積雪的沉積通量，即單位時間內在單位面積上的降雪質量，新降積雪密度為100 kg/m3時，這里將Gg取為1.5×10-4kg/(m2·s)，對應的氣象條件為暴雪天氣[35]；T/s 為屋面風吹雪的持續時間；qr/(kg/(m2·s))為屋面積雪的侵蝕通量。圖14 為不同風速下屋面跨度對雪荷載暴露系數的影響。由圖可得，屋面跨度相同時，雪荷載暴露系數隨著來流風速的增大而減??；并且在同樣的氣候條件(風速)下，跨度越大的屋蓋結構，來流風速對其屋面雪荷載的折減就越小。

4 結論

本文采用CFD 方法模擬了平屋面的摩擦速度，并結合其分布特征分析了屋面積雪傳輸率、侵蝕通量和暴露系數隨跨度的變化趨勢，主要結論如下：

(1) 本文根據風洞試驗對Irwin 探頭測量摩擦速度進行了標定，結合Irwin 探頭風洞試驗對比不同湍流模型的CFD 模擬結果發現，Realizablek-ε模型的模擬結果與風洞試驗最為接近。

(2) 對比不同跨度平屋面的CFD 模擬結果可得，隨著屋面跨度的增大，屋面摩擦速度的最小值所在的位置和積雪沉積的區域均逐漸向迎風側屋檐靠近。不同跨度平屋面的平均摩擦速度隨著風速的增大而增大，并且屋面跨度與高度之比L/H=4 時，屋面的平均摩擦速度在本文選取的4 種跨度的屋面中最小。

(3) 屋面積雪的平均侵蝕通量隨著跨度的增大而減小。并且在同樣的氣候條件(風速)下，跨度越大的屋蓋結構，來流風速對其屋面雪荷載的折減就越小。