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

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

0 引 言

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

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

1 試驗(yàn)方案與方法

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

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

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

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

一般而言，網(wǎng)格越密集，GCI 越小，代表計(jì)算結(jié)果越精確。GCI 法至少需要從致密到粗糙的3 種網(wǎng)格，本試驗(yàn)計(jì)算網(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 法的計(jì)算步驟為：

（1）計(jì)算收斂精度p。

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

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

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

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

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

3.1 水面線

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

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

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

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

3.2 流 速

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

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

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

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

3.3 湍動(dòng)能

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

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

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

圖7 網(wǎng)格2計(jì)算湍動(dòng)能對(duì)應(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)

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

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

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

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

4.1 水面線驗(yàn)證

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

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

圖10 數(shù)值模擬計(jì)算與物模試驗(yàn)實(shí)測(cè)的臺(tái)階段上水面線對(duì)比Fig.10 Comparison between numerical simulation calculation and physical model test measurement of water surface line on stepped spillway

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

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

4.2 流速驗(yàn)證

hk/h=2.0 流量工況下，分別取1.25 m 均勻步長(zhǎng)臺(tái)階方案在5、10、15、20 號(hào)臺(tái)階頂豎直斷面的數(shù)值模擬計(jì)算流速與物模試驗(yàn)實(shí)測(cè)流速進(jìn)行對(duì)比（圖11）。

圖11 數(shù)值模擬計(jì)算與物模試驗(yàn)實(shí)測(cè)的臺(tái)階頂豎直斷面流速對(duì)比Fig.11 Comparison between numerical simulation calculation and physical model test measurement of vertical section of convex corner

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

4.3 初始摻氣點(diǎn)驗(yàn)證

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

圖12 物模試驗(yàn)實(shí)測(cè)的初始摻氣點(diǎn)（hk/h=1.8）Fig.12 The initial aeration point measured by the physical model test （hk/h=1.8）

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

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

表1 數(shù)值模擬計(jì)算與物模試驗(yàn)實(shí)測(cè)的初始摻氣點(diǎn)位置對(duì)比Tab.1 Comparison of the initial aeration point between numerical simulation and physical model test

5 結(jié) 論

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

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

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

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