多因素下開洞口煤倉風致積雪CFD模擬及擬合研究*

殷志祥，姜昕彤，崔翰博

(遼寧工程技術大學 力學與工程學院，遼寧 阜新 123000)

0 引言

隨著大跨度建筑的跨度越來越大，由于屋面積雪分布不均勻且局部雪壓過大，風雪運動下結構屋面發(fā)生損毀事故增多[1-4]。我國《建筑結構荷載規(guī)范》(GB 50009-2012)分別對拱形屋面和大跨度屋面(跨度>100 m)的積雪分布系數做出規(guī)定，但規(guī)定較為簡單，并不適用于開口煤倉的球殼屋蓋，且在風雪共同作用下此類開口結構表面的積雪運動變化將更為復雜。因此，尋找適合風雪作用下開洞口煤倉球殼屋蓋的積雪分布系數十分必要。

Murakami等[5]分別采用k-ε湍流模型和大渦模型(LES)對立方塊周圍的湍流分布情況進行數值模擬，并與風洞試驗數據對比，發(fā)現 LES 模擬結果較為準確； Uematsu等[6]對防雪欄、風穴和小山3種模型模擬時考慮了風致雪漂的躍移和懸移運動，得到結果更接近實地觀測； Beyers等[7]在對立方體模擬時細化影響因素，精確得到風雪作用下雪顆粒運動軌跡；周暄毅等[8]模擬計算了北京南站屋蓋表面的積雪不均勻分布情況；孫芳錦等[9]模擬研究風向角對膜結構表面積雪漂移運動的影響，得到膜結構表面不均勻雪壓分布情況；孫曉穎等[10]對比不同湍流模型下防雪欄周圍風致雪漂移運動，得出k-kl-ω模型對風雪運動模擬較為準確；李躍等[11]研究矢跨比和風速對屋面雪壓的影響，并求得近似公式。以上研究都以具體模型展開分析，沒有適合此類開口結構的風雪研究。本文以開洞口圓形煤倉為模型，研究風雪共同作用下屋蓋不均勻雪壓分布，根據所得模擬結果推導出此類開口結構屋面積雪值公式，可直接應用于此后煤倉工程設計中，對今后準確預測此類開洞口結構表面雪壓分布有一定的參考價值。

1 風致積雪CFD理論

1.1 空氣相及雪相的控制方程

采用Reynolds平均法，風可近似為不可壓縮性粘性流體，質量守恒方程及動量守恒方程[12]分別為：

(1)

(2)

(3)

式中：ρs為積雪密度，kg/m3；f為單位體積里雪相所占組分；uj為風速矢量；μt為空氣相湍流粘性系數，此系數代表空氣相對雪相的影響程度；uR，j為雪相對空氣的運動速度，m/s，取0.3[13]。

1.2 雪的侵蝕和沉積

摩擦速度u*(風速)是體現屋蓋積雪是否發(fā)生侵蝕或沉積的決定性因素。u*的表達式為：

u*=u(z)K/ln(z/zs)

(4)

式中：u(z)為近壁面的風速，m/s；K為Karman常數，取0.4；z為積雪面到近壁的第一格柵中心高度，m；zs為積雪面粗糙度，m，取0.000 5[9]。屋蓋表面積雪發(fā)生侵蝕時摩擦速度u*大于閾值速度u*t(0.20 m/s[8])，沉積時摩擦速度小于閾值速度。侵蝕量和沉積量[7]分別如式(5)～(6)所示：

(5)

(6)

式中：Aero為比例常數，取7.0×10-4；C為單位體積內雪的質量，kg，C=fρs其中，f為雪的體積分數，ρs為積雪密度，kg/m3；wf為降雪速度，m/s。

2 開洞口煤倉風致積雪CFD模擬

2.1 計算模型與網格劃分

本文以某電廠圓形煤倉為工程背景，模型如圖1所示，煤倉筒倉高度H為18 m[14]，煤倉帽頂直徑d為20 m，高度h為2 m，煤倉跨度D分別設置為110 m和50 m，球殼屋面上正方形洞口模擬棧橋洞口。110 m跨度煤倉洞口尺寸為4 m×4 m[14]，50 m跨度煤倉洞口尺寸為2 m×2 m，門洞尺寸為4.5 m×4.5 m[15]。

圖1 開洞口煤倉模型Fig.1 Opening round coal bunker model diagram

110 m跨度煤倉計算流域為1 000 m×500 m×300 m，50 m跨度煤倉計算流域為600 m×350 m×270 m，滿足阻塞率要求?；旌暇W格劃分模型，非結構化網格劃分建筑物，結構化網格劃分外部區(qū)域。最小網格尺寸0.5 m，網格總數約90～150萬個，近壁面網格劃分如圖2所示。

圖2 近壁面網格劃分Fig.2 Near wall meshing

2.2 邊界條件與求解方法

入流面空氣相采用速度入口邊界條件，風速u(z)為：

u(z)=ub(z/zb)α

(7)

式中：zb為標準參考高度，m，取10；ub為zb對應的平均風速，m/s；z為任一高度，m；α為地面粗糙度指數，取0.15[16]。湍流動能k及湍流耗散率ε為：

k=1.5(u(z)·I)2，ε=0.090.75·k1.5/l

(8)

式中：I為湍流強度；l為湍流積分尺度，m，I和l分別為：

(9)

按日本規(guī)范取II類地貌，zb取5m；zG取350 m。

(10)

入口處雪相邊界條件為：

(11)

式中：μp為雪相運動速度，m/s；g為重力加速度，m/s2；hs為躍移層和懸移層臨界高度，m，為簡化計算取為0.5。

出流面為完全發(fā)展出流邊界條件，模型表面與地面為無滑移的壁面條件，流域上部和兩側為自由滑移避免條件。將空氣和雪顆粒之間關系假設為單項耦合，采用Euler-Euler方法，認為空氣相和雪相都是連續(xù)相，選擇計算量小且計算穩(wěn)定性較好的Mixture多相流模型，采用可實現k-ε兩方程模型。

3 開洞口煤倉風致積雪參數影響分析

模擬在不同風向角(θ=0°，20°，45°)、不同積雪厚度(h=100，300，500 mm)、下跨度為50 m和110 m的開洞口圓形煤倉表面積雪分布情況，煤倉矢跨比分別設置為α=0.14，0.32，0.45的共54組工況(模型1～54)，模型編號如表1所示。

3.1 開洞口煤倉表面測點設置

在煤倉表面設4個測點，測點A和B分別位于煤倉中軸線背風面和迎風面最外側，測點C位于垂直風向方向中軸線最外側，即表面積雪侵蝕量最大區(qū)域，因棧橋洞口迎風側積雪侵蝕量大于背風側，故測點D設置在迎風側棧橋洞口邊緣，如圖3所示。

圖3 測點設置Fig.3 Measuring point setting diagram

3.2 風向角對開口煤倉積雪分布的影響及測點分析

取模型9、模型18、模型27的3種工況進行分析。風向角角度定義如圖4(a)所示，因篇幅有限，僅以模型9為例，得出的開口煤倉表面單位時間單位屋蓋面積上的積雪侵蝕或沉積量如圖4(b)所示，來流風向在圖中標出，其他工況與之相同。由圖4得出，煤倉表面積雪均發(fā)生侵蝕，迎風面到背風面積雪侵蝕量逐漸增大后變小，迎風面積雪侵蝕量略大于背風面，中心部位積雪侵蝕量達到最大，這與文獻[11]得出結論相吻合，迎風面邊緣及由于帽頂遮擋作用背風面中心部位積雪侵蝕量極小。棧橋洞口背風側比迎風側積雪侵蝕量小，是因為棧橋洞口可以吸收部分來流風，減小了風對表面積雪的侵蝕作用。

表1 模型編號與參數Table 1 Model number and parameter

圖4 不同風向角下單位時間單位面積積雪侵蝕或沉積量Fig.4 Snow erosion or sedimentation volume per unit area per unit time at different wind direction angles

圖5 煤倉表面積雪量分布(模型9，18，27)Fig.5 Coal bunker surface snow volume distribution curve(model 9，18，27)

為更直觀地分析風向角對開洞口煤倉的影響，匯總數據見圖5，以與順風向平行的中軸線AB為軸，每隔13.75 m設置1個測點，對測點1～ 9進行分析，測點1與B點重合，測點9與A點重合。由圖5可知，隨風向角增大，迎風面(測點1～ 5)積雪侵蝕量明顯增加，背風面(測點5～ 9)增長幅度小。測點3和C點在45°風向角下侵蝕量比0°風向角分別增大80%和69.23%，說明風向角對開口煤倉迎風面風致積雪漂移影響尤為明顯；45°風向角下D點比0°風向角增大66.67%，說明風向角對棧橋洞口附近積雪變化影響也較大，在此后煤倉設計中應重點考慮風向角對煤倉表面積雪的影響。

3.3 積雪厚度對開口煤倉積雪分布影響及測點分析

取模型30、模型33、模型36為例進行分析，煤倉表面積雪量分布曲線，如圖6所示。

圖6 煤倉表面積雪量分布(模型30，33，36)Fig.6 Coal bunker surface snow volume distribution curve(model 30，33，36)

由圖6可知，3種工況煤倉表面積雪侵蝕變化規(guī)律大致相同，積雪侵蝕量隨積雪厚度增大而增大。位于迎風面的測點2、測點3在500 mm積雪厚度的侵蝕量比100 mm積雪厚度分別增大28.57%和25%，背風面測點7、測點8則分別增大33.33%和50%，位于侵蝕量最大區(qū)域的C點增大29.41%，測點D則增大28.57%，由此得出積雪厚度對開洞口煤倉表面的風致積雪漂移影響也很大。

3.4 矢跨比對開口煤倉積雪分布的影響及測點分析

對110 m跨度下不同矢跨比、風向角、積雪厚度的煤倉模型進行數值模擬后發(fā)現，矢跨比為0.14時，煤倉表面積雪變化規(guī)律不同于其他情況，取模型25說明，如圖7所示。由圖7可知，除棧橋洞口背風側邊緣區(qū)域積雪發(fā)生侵蝕，其余區(qū)域積雪全部沉積，因0.14矢跨比的煤倉網殼坡度極小，相當于平屋蓋，故積雪全部沉積在表面，棧橋洞口背風側附近風速較大，一部分積雪發(fā)生侵蝕。

圖7 單位時間單位面積積雪侵蝕或沉積量Fig.7 Snow erosion or sedimentation volume per unit area

以50 m跨度的模型28、模型29、模型30進行分析，得出煤倉表面積雪量分布曲線如圖8所示。由圖8可知，50 m跨度的煤倉積雪侵蝕量和棧橋洞口邊緣測點D的積雪侵蝕量隨矢跨比增大而增加，且風向角的改變不影響積雪變化規(guī)律?？缍葹?10 m，風向角0°下3種工況(模型4、模型5、模型6)，風向角20°下3種工況(模型13、模型14、模型15)模擬得出的結果如圖9所示。由圖9可知，110 m跨度、矢跨比為0.14的煤倉表面積雪大部分發(fā)生沉積，且積雪隨風向角變化及其微弱。由圖9(a)可知，0°風向角下、矢跨比為0.32的煤倉表面和測點C，D的積雪侵蝕量小于矢跨比為0.45的積雪侵蝕量；由圖9(c)看出，風向角20°下此種情況正好相反。通過圖8和圖9得出，煤倉設計時不應僅考慮矢跨比這1種參數的影響，應把風向角對不同矢跨比的影響考慮其中。

圖8 煤倉表面積雪量分布(模型28,29,30)Fig.8 Coal bunker surface snow volume distribution curve(model 28,29,30)

圖9 不同工況下煤倉表面積雪量分布Fig.9 Coal bunker surface snow volume distribution curve in different conditon

4 開洞口煤倉積雪侵蝕量與沉積量公式擬合

4.1 開洞口煤倉表面積雪侵蝕量與沉積量公式擬合

由積雪量分布曲線圖可以看出，煤倉中軸線背風面測點6～9的積雪值基本一致，且與測點1大致相同，故僅對迎風面測點1～6和測點C，D進行公式擬合。煤倉球殼屋面積雪呈等值線分布，故進行等值線擬合，球殼屋面上的積雪量q可以看作是中軸線測點坐標x的函數，發(fā)現測點1～6積雪值分布符合二次函數分布規(guī)律，故進行最小二乘法擬合。

(12)

式中：q(x)為煤倉屋面各點積雪量，kg·m-2·s-1；x為屋面測點坐標，m；a，b，c為待定常數。

二次函數開口大小和豎向位置分別由橫向因子和豎向因子決定。積雪量q(x)是由風向角、積雪厚度和矢跨比共同決定的，為了準確地反應其比重大小引入此概念，擬合出全部工況下的二次函數，得出不同工況下的橫向因子和豎向因子，橫向和豎向因子會隨著參數不同而發(fā)生變化，故用橫向波動因子來表示橫向因子的變化情況，橫向波動因子如式(13)所示。

(13)

式中：Kwave-hz為橫向波動因子；max(hz)為最大橫向因子；min(hz)為最小橫向因子；n為工況個數；hz(i)為不同工況下的橫向因子。豎向波動因子同理可求。由此得出3種參數的波動因子如表2所示。修正表2各數據，如表3所示。

表2 各參數影響下煤倉屋面積雪值波動因子Table 2 The snow value fluctuation factor of coal bunker under the influence of various parameters

表3 修正后各參數所占比例Table 3 Proportion of parameters after modification

由表3可知，矢跨比和風向角的橫向和豎向波動因子都不可忽略，經過對比分析，將測點1～6公式擬合成式(14)，各測點系數如表4所示。

(14)

式中：q1～6為各測點積雪值，kg·m-2·s-1；α為矢跨比；θ為風向角，°；x為測點坐標，m；h為積雪厚度，m；a～e為待定系數。

測點C積雪值qC擬合公式為：

qC=-9.449×10-5α2·θ-2.400×10-6θ2+0.001α-5.657×10-5h-7.590×10-4

(15)

測點D積雪值qD擬合公式為：

qD=1.569×10-4α2·θ-1.204×10-4θ- 0.003α-4.277×10-5h+0.002 (16)

4.2 公式擬合結果驗證

為驗證擬合公式的準確性，與文獻[11]模擬結果進行對比，文獻為跨度60 m、矢跨比0.22的球殼屋蓋，對比結果如圖10所示，因文獻的球殼屋蓋頂部無帽頂，且根據文獻模擬結果其屋面積雪基本呈對稱分布，故僅擬合半面屋蓋積雪值，因其表面無洞口，故僅對式(14)和式(15)的準確性進行驗證。

由圖10得出，文獻[11]網殼邊緣到中心積雪絕對值逐漸增大，中心部位積雪絕對值最大，與本文模擬結論一致，從公式擬合結果來看，擬合規(guī)律一致，由于文獻中沒有考慮積雪厚度對球殼屋面的影響，故公式擬合結果與文獻模擬數值有一定差別。

圖10 文獻模擬與擬合公式結果對比Fig.10 Comparison of simulation and fitting formula results

5 結論

1)開口煤倉屋面積雪量由迎風面到背風面逐漸增大然后變小，積雪侵蝕量最大位置在煤倉中心部位，隨風向角和積雪厚度增大而變大，風向角對其屋蓋積雪量影響不容忽視，積雪厚度對其也有一定影響；矢跨比變化下的屋蓋積雪值應同時考慮風向角對其影響；煤倉棧橋洞口背風側積雪侵蝕量小于迎風側，在此后煤倉設計中也應注意洞口周圍積雪變化情況。

2)推導出開口圓形煤倉屋面積雪量與風向角、積雪厚度、矢跨比相關公式，通過不同擬合公式可直接求出煤倉表面最不利積雪值、棧橋洞口附近積雪值及屋面其他各點積雪值。

3)擬合公式求得結果與前人模擬結果比較，發(fā)現兩者屋面積雪量變化規(guī)律完全一致，最大積雪侵蝕量都在網殼中心部位，且數值相差不大，擬合效果較為良好。此公式彌補了規(guī)范中沒有涉及到的開洞口圓形煤倉屋面積雪分布系數，可直接應用于工程設計中，對今后準確預測此類開洞口結構表面雪壓分布具有重要的參考價值。