王海波，沈立群，張晉鋒，武英豪

（1.湖北省水利水電規(guī)劃勘測設(shè)計院，湖北 武漢 430064；2.武漢大學(xué)水利水電學(xué)院，湖北 武漢 430000）

0 引 言

得益于碾壓混凝土技術(shù)（RCC）的應(yīng)用，臺階式溢洪道已經(jīng)成為國內(nèi)外水利工程中通用的泄流消能方式。由于水流沿程逐級摻氣、減速、消能，與光滑溢洪道相比，臺階式溢洪道能夠顯著減輕下游的消能壓力?？諝鈸饺胨餍纬傻乃嗔鲝?fù)雜運動是臺階式溢洪道數(shù)值模擬計算的難點。

隨著CFD 領(lǐng)域和計算機技術(shù)的發(fā)展，數(shù)值誤差評估方法也在不斷更新迭代，目前最可靠的是RE（Richardson Extrapolation）法［1，2］。自從它的創(chuàng)始人于1911年首次提出該方法以來，大量學(xué)者對其進行研究并應(yīng)用到數(shù)值誤差評估中［3，4］。美國機械工程師協(xié)會（American Society of Mechanical Engineers，ASME）的流體工程部一直致力于CFD 數(shù)值誤差的檢測和評估的工作。CFD 可靠性分析領(lǐng)域的第一套質(zhì)量控制措施是由該部門的Roache 等人［5］于1986年發(fā)布，并由Freitas［6］于1993年進行修訂。它就是基于RE 法的計算網(wǎng)格收斂指數(shù)法（Grid Convergence Index，GCI），該方法已應(yīng)用在數(shù)百個CFD 案例中［7］。采用從致密至粗糙3 種網(wǎng)格尺寸，利用Flow 3D?軟件對臺階式溢洪道上的水汽二相流進行數(shù)值模擬，提供評價其計算精度的一般方法。采用從致密至粗糙3 種網(wǎng)格尺寸，利用Flow 3D?軟件對臺階式溢洪道上的水汽二相流進行數(shù)值模擬，提供評價其計算精度的一般方法。

1 試驗方案與方法

本次數(shù)值模擬的計算域包括上游水庫段、寬頂堰段、臺階式溢洪道段和尾水渠段。為保障計算精度的同時加快計算速率，數(shù)值模擬采用二維網(wǎng)格進行計算。以水庫底部起點為坐標(biāo)原點，上游至下游水平方向為X軸正方向，豎直向上為Y軸正方向，計算模型的全局體型見圖1。為保障上游水庫水位穩(wěn)定，水庫長度取為37.5 m；為保障臺階段入流平順，在其前方設(shè)置長度為20.625 m 的寬頂堰；尾水渠段長度取為12.5 m。臺階高度為1 m，步長為1.25 m，坡度為38.66°，共25級臺階。

圖1 計算模型全局示意圖（單位：m）Fig.1 Global schematic diagram of the computational model

Flow 3D?軟件是基于混合框架并兼具模擬以上4 種物理現(xiàn)象的能力。它采用GMRES 方法［8］求解離散方程，采用Favor 技術(shù)［9，10］進行網(wǎng)格處理，可實現(xiàn)用簡單的結(jié)構(gòu)化網(wǎng)格描述復(fù)雜的幾何體型［11］，并開發(fā)了Tru-VOF 方法［12］，大幅改進了其追蹤自由液面的速度和精度。因此，本文選用Flow 3D?軟件對臺階上的水汽二相流進行數(shù)值模擬計算。

2 網(wǎng)格收斂指數(shù)計算步驟

一般而言，網(wǎng)格越密集，GCI 越小，代表計算結(jié)果越精確。GCI 法至少需要從致密到粗糙的3 種網(wǎng)格，本試驗計算網(wǎng)格收斂性分析取網(wǎng)格尺寸為：0.037 5 m×0.037 5 m（網(wǎng)格1：h1×h1）、0.05 m×0.05 m（網(wǎng)格2：h2×h2）、0.066 7 m×0.066 7 m（網(wǎng)格3：h3×h3），即網(wǎng)格加密因子r取為：

大于推薦的最小網(wǎng)格加密因子1.3［13］。GCI 法的計算步驟為：

（1）計算收斂精度p。

式中：φk為第k種網(wǎng)格下數(shù)值模擬計算的離散解。本試驗中r21=r32=1.333，可得q(p)=0，進而簡化收斂精度p的計算。

（2）計算網(wǎng)格2相對于細(xì)密網(wǎng)格（網(wǎng)格1）的相似誤差。

（3）計算網(wǎng)格2相對于細(xì)密網(wǎng)格（網(wǎng)格1）的網(wǎng)格收斂指數(shù)。

3 計算網(wǎng)格收斂性分析

臺階式溢洪道中，水面線及流速是描述水流的基本參數(shù)，湍動能是決定水流自摻氣是否發(fā)生的關(guān)鍵水力參數(shù)，在數(shù)值模擬計算中必須保證這3個參數(shù)的準(zhǔn)確性。本節(jié)以1.25 m均勻步長臺階方案在hk∕h=2.0 流量工況下的水面線、流速和湍動能作為特征離散解，采用GCI方法進行計算網(wǎng)格收斂性分析。

3.1 水面線

Flow 3D?軟件在3 種網(wǎng)格下計算所得的水庫及寬頂堰水面線［圖2（a）］、臺階段水面線［圖2（b）］均貼合緊密，即3 種網(wǎng)格的整體水面線之間相差甚微，僅在溢洪道后半段觀察到輕微的離散。將5、10、15、20 號臺階上方的水面線進行局部放大［圖2（c）～（f）］，可知3 條水面線之間在15 號臺階前貼合緊密，在20號臺階上開始觀察到輕微的離散：細(xì)密網(wǎng)格（網(wǎng)格1）的水面線最高，粗糙網(wǎng)格（網(wǎng)格3）的水面線最低，網(wǎng)格2 的水面線介于兩者之間。網(wǎng)格越細(xì)致，摻氣模型捕捉到的摻氣就越多，摻氣水深增高，水面線被抬升，而臺階式溢洪道的摻氣充分發(fā)展區(qū)一般位于后半段，這才使得3 種網(wǎng)格的水面線之間在20 號臺階之后出現(xiàn)輕微的離散。

圖2 不同網(wǎng)格尺寸的計算水面線對比Fig.2 Comparison of water surface lines calculated with different grid sizes

網(wǎng)格2 計算水面線與細(xì)密網(wǎng)格（網(wǎng)格1）計算水面線之間的偏差已十分微小。以20號臺階上水面線數(shù)據(jù)為離散特征解，采用GCI方法對其進行數(shù)值誤差評估（圖3），得出網(wǎng)格2的全局網(wǎng)格收斂指數(shù)小于5%。最大GCI21為4.01%，出現(xiàn)在x=82.125 m，對應(yīng)細(xì)密網(wǎng)格計算水深為2.23 m。

圖3 網(wǎng)格2計算水面線對應(yīng)的網(wǎng)格收斂指數(shù)（20號臺階）Fig.3 The GCI corresponding to the water surface line calculated by NO.2 grid

3.2 流 速

Flow 3D?軟件采用3種網(wǎng)格計算所得的5、10、15、20號臺階頂（凸角處）豎直斷面上的流速分布如圖4，均呈現(xiàn)由底部向水面逐漸增大的正確分布規(guī)律。網(wǎng)格2斷面流速曲線與細(xì)密網(wǎng)格（網(wǎng)格1）斷面流速曲線之間貼合緊密，而與粗糙網(wǎng)格（網(wǎng)格3）斷面流速曲線之間存在明顯偏差。細(xì)密網(wǎng)格能夠更加充分考慮水流摻氣對流速發(fā)展的影響，更加精確捕捉空間上流速的變化。一般而言，同一位置處細(xì)密網(wǎng)格計算流速略大于粗糙網(wǎng)格計算流速，且細(xì)密網(wǎng)格的斷面流速梯度大于粗糙網(wǎng)格的斷面流速梯度。

圖4 不同網(wǎng)格尺寸的臺階頂豎直斷面計算流速對比Fig.4 Comparison of calculated velocity in vertical section of convex corner with different grid sizes

網(wǎng)格2計算流速與細(xì)密網(wǎng)格計算流速之間的偏差已十分微小。以5、10、15、20 號臺階頂豎直斷面上的流速數(shù)據(jù)為離散特征解，采用GCI 方法對其進行數(shù)值誤差評估（圖5），得出網(wǎng)格2的全局網(wǎng)格收斂指數(shù)小于5%。5號臺階最大GCI21為1.56%，出現(xiàn)在y+=0.3 m，對應(yīng)細(xì)密網(wǎng)格計算流速為8.82 m∕s；10 號臺階最大GCI21為2.33%，出現(xiàn)在y+=0.1 m，對應(yīng)細(xì)密網(wǎng)格計算流速為7.46 m∕s；15號臺階最大GCI21為3.53%，出現(xiàn)在y+=0.2 m，對應(yīng)細(xì)密網(wǎng)格計算流速為8.88 m∕s；20號臺階最大GCI21為3.29%，出現(xiàn)在y+=0.9 m，對應(yīng)細(xì)密網(wǎng)格計算流速為14.58 m∕s。

圖5 網(wǎng)格2計算流速對應(yīng)的網(wǎng)格收斂指數(shù)Fig.5 The GCI corresponding to the velocity calculated by NO.2 grid

3.3 湍動能

Flow 3D?軟件采用3種網(wǎng)格計算所得的5、10、15、20號臺階頂豎直斷面上的湍動能分布如圖6，均呈現(xiàn)由底部向水面逐漸減小的正確分布規(guī)律。與流速計算結(jié)果類似，網(wǎng)格2 斷面湍動能曲線與細(xì)密網(wǎng)格（網(wǎng)格1）斷面湍動能曲線之間緊密貼合，而與粗糙網(wǎng)格（網(wǎng)格3）斷面湍動能曲線之間存在明顯差異。細(xì)密網(wǎng)格能夠通過Favor 技術(shù)更加準(zhǔn)確還原臺階的體型輪廓，更加精確進行湍流模型的計算。同一位置處細(xì)密網(wǎng)格計算湍動能大于粗糙網(wǎng)格計算湍動能。

圖6 不同網(wǎng)格尺寸的臺階頂豎直斷面計算湍動能對比Fig.6 Comparison of calculated turbulent kinetic energy in vertical section of convex corner with different grid sizes

網(wǎng)格2計算湍動能與細(xì)密網(wǎng)格計算湍動能之間的偏差已十分微小。以5、10、15、20 號臺階頂豎直斷面上的湍動能數(shù)據(jù)為離散特征解，采用GCI 方法對其進行數(shù)值誤差評估（圖7），得出網(wǎng)格2 的全局網(wǎng)格收斂指數(shù)小于5%。5 號臺階最大GCI21為4.1%，出現(xiàn)在y+=0.8 m，對應(yīng)細(xì)密網(wǎng)格計算湍動能為0.18 m2∕s2；10 號臺階最大GCI21為3.4%，出現(xiàn)在y+=0.6 m，對應(yīng)細(xì)密網(wǎng)格計算湍動能為1.46 m2∕s2；15號臺階最大GCI21為3.84%，出現(xiàn)在y+=0.6 m，對應(yīng)細(xì)密網(wǎng)格計算湍動能為2.48 m2∕s2；20 號臺階最大GCI21為2.56%，出現(xiàn)在y+=1.0 m，對應(yīng)細(xì)密網(wǎng)格計算湍動能為2.02 m2∕s2。

圖7 網(wǎng)格2計算湍動能對應(yīng)網(wǎng)格收斂指數(shù)Fig.7 The GCI corresponding to the turbulent kinetic energy calculated by NO.2 grid

3.4 網(wǎng)格尺寸確認(rèn)

通過分別以水面線、流速和湍動能作為離散特征解進行的網(wǎng)格收斂性分析可知，網(wǎng)格2 計算結(jié)果與細(xì)密網(wǎng)格（網(wǎng)格1）計算結(jié)果之間的偏差較小，3 種水力參數(shù)對應(yīng)的全局GCI21均小于5%；粗糙網(wǎng)格（網(wǎng)格3）計算水面線與細(xì)密網(wǎng)格（網(wǎng)格1）計算水面線之間的偏差也較小，但兩者的計算流速和湍動能之間存在較大偏差。為保障計算精度的同時縮減計算時長，本文所涉及的數(shù)值模擬均采用網(wǎng)格2（0.05 m×0.05 m）尺寸進行計算，即每個臺階在高度方向上被均分為20個網(wǎng)格，如圖8。

圖8 推薦計算網(wǎng)格尺寸示意圖Fig.7 Schematic diagram of recommended calculation grid size

4 數(shù)值模擬精確性分析

本文以張文傳［14］進行的1.25 m 均勻步長臺階方案的物理模型試驗測量數(shù)據(jù)作為參照，選用水面線、流速和初始摻氣點作為比較參數(shù)，評價數(shù)值模擬計算結(jié)果的精確度。物理模型長度比尺為1∶12.5，模型臺階高度和步長分別為8 cm和10 cm。

4.1 水面線驗證

hk/h=2.0流量工況下，1.25 m均勻步長臺階方案通過數(shù)值模擬計算和物理模型試驗測量所得的寬頂堰上水面線、臺階段上水面線分別如圖9、圖10，其中B表示寬頂堰長度。

圖9 數(shù)值模擬計算與物模試驗實測的寬頂堰上水面線對比Fig.9 Comparison between numerical simulation calculation and physical model test measurement of water surface line on wide-top weir

圖10 數(shù)值模擬計算與物模試驗實測的臺階段上水面線對比Fig.10 Comparison between numerical simulation calculation and physical model test measurement of water surface line on stepped spillway

寬頂堰上數(shù)值模擬計算水面線與物模試驗實測水面線之間整體貼合緊密，僅在前半部分出現(xiàn)少許離散，后半部分及溢洪道入流部分幾乎無偏差。兩者的最大偏差出現(xiàn)在x∕B=0.1處，為2.8%，對應(yīng)物模試驗實測水深為2.613 m。

臺階段上數(shù)值模擬計算水面線與物模試驗實測水面線之間整體貼合緊密，兩者在臺階段前半段幾乎無偏差，卻在后半段出現(xiàn)少許離散。前半段水流未摻氣，為純水流流動，或摻氣濃度較低，數(shù)值模擬計算精度較高，同一位置處計算水面線略低于物模試驗實測水面線，兩者偏差均不超過5%；后半段水流摻氣逐漸發(fā)展，水汽二相流的運動特性更加復(fù)雜，數(shù)值模擬計算精度有所下降，同一位置處計算水面線略高于物模試驗實測水面線，兩者最大偏差出現(xiàn)在x=89.375 m 處，為12%，對應(yīng)物模實測水深為1.5 m。

4.2 流速驗證

hk/h=2.0 流量工況下，分別取1.25 m 均勻步長臺階方案在5、10、15、20 號臺階頂豎直斷面的數(shù)值模擬計算流速與物模試驗實測流速進行對比（圖11）。

圖11 數(shù)值模擬計算與物模試驗實測的臺階頂豎直斷面流速對比Fig.11 Comparison between numerical simulation calculation and physical model test measurement of vertical section of convex corner

數(shù)值模擬計算與物模試驗實測的臺階頂豎直斷面流速分布規(guī)律相同，即自底部至水面的流速逐漸增大、流速梯度逐漸減小。在水流底部，物模試驗實測流速大于數(shù)值模擬計算流速，在水流表面，卻正好相反，即在同一豎直斷面上，物模試驗實測的流速梯度小于數(shù)值模擬計算的流速梯度，與前人研究結(jié)果一致［11］。物模試驗實測流速與數(shù)值模擬計算流速之間在5號臺階頂豎直斷面上的偏差很小，不超過5%，后在10、15、20號臺階頂豎直斷面上的偏差逐漸增大。數(shù)值模擬計算未充分考慮水流摻氣導(dǎo)致水流底部與臺階面之間的摩擦減小這一效果，隨著水流摻氣的沿程發(fā)展，這一缺陷會導(dǎo)致數(shù)值模擬計算流速與物模試驗實測流速之間的偏差被逐漸拉大。兩者之間的最大流速偏差出現(xiàn)在20 號臺階底部，為16%，對應(yīng)物模試驗實測流速為9.0 m∕s。

4.3 初始摻氣點驗證

初始摻氣點作為水流表面摻氣的起點，是描述水流自摻氣的關(guān)鍵參數(shù)，需確保CFD 軟件對其位置模擬的可靠性。hk/h=1.2、hk/h=1.4、hk/h=1.6、hk/h=1.8、hk/h=2.0 五種流量工況下，分別取1.25 m 均勻步長臺階方案數(shù)值模擬計算的初始摻氣點位置與物模試驗實測的初始摻氣點位置進行對比，以評價Flow 3D?軟件內(nèi)嵌摻氣模型的計算精度。初始摻氣點后水流迅速摻氣，湍動加劇，空氣和水混合，呈現(xiàn)為乳白色的水汽二相流［14］（圖12）。物模試驗中可將水流表面開始變白的位置作為初始摻氣點。數(shù)值模擬可繪制水流的摻氣濃度云圖（圖13），可以清晰觀察到初始摻氣點的位置［15，16］。

圖12 物模試驗實測的初始摻氣點（hk/h=1.8）Fig.12 The initial aeration point measured by the physical model test （hk/h=1.8）

圖13 數(shù)值模擬計算的初始摻氣點（均勻1.25 m步長）Fig.13 Initial aeration point calculated by numerical simulation

統(tǒng)計各工況下數(shù)值模擬計算的初始摻氣點和物模試驗實測的初始摻氣點對應(yīng)臺階級數(shù)如表1，隨著單寬流量的增大，初始摻氣點的位置逐漸下移。數(shù)值模擬計算初始摻氣點與物模試驗實測初始摻氣點的位置偏差在1 個臺階長度范圍內(nèi)，F(xiàn)low 3D?軟件內(nèi)嵌摻氣模型的計算精度較高。

表1 數(shù)值模擬計算與物模試驗實測的初始摻氣點位置對比Tab.1 Comparison of the initial aeration point between numerical simulation and physical model test

5 結(jié) 論

本文提供了評價臺階式溢洪道上水汽二相流數(shù)值模擬計算精度的一般方法，包括計算網(wǎng)格的收斂性分析和計算結(jié)果的精確性分析。主要結(jié)論有：

（1）采用GCI 方法進行網(wǎng)格收斂性分析至少需要從致密到粗糙的3種網(wǎng)格，推薦使用相同的網(wǎng)格加密因子，以簡化網(wǎng)格收斂指數(shù)的計算。

（2）本文選用了計算水面線、流速及湍動能作為特征離散解，得出網(wǎng)格2 的全局GCI 均小于5%，推薦選用此網(wǎng)格尺寸進行數(shù)值模擬計算，在保障計算精度的同時縮減計算時長。

（3）通過與物理模型試驗的對照，F(xiàn)low 3D?軟件對臺階式溢洪道上關(guān)鍵水力學(xué)參數(shù)（水面線、流速及初始摻氣點位置）的計算精度較高，且對純水流的計算精度優(yōu)于對水汽二相流的計算精度，兩者均位于較高水準(zhǔn)。