基于FLOW-3D的臺(tái)階式溢流堰摻氣水流數(shù)值模擬

董宗師，蔡 芳

（1.長(zhǎng)江設(shè)計(jì)集團(tuán)有限公司，武漢 430010；2.四川大學(xué)水利水電學(xué)院，成都 610065）

0 引 言

流動(dòng)摻氣現(xiàn)象廣泛存在于自然、工業(yè)和工程之中，一種較具代表性的形式即為水面自摻氣［1，2］。在這種流動(dòng)中，空氣從水面摻入、形成大小不同的氣泡分布于水流之中［3］。20世紀(jì)80年代末，碾壓混凝土筑壩技術(shù)的發(fā)展促進(jìn)了臺(tái)階式溢流堰的工程應(yīng)用及其水力特性研究。同時(shí)，得益于新型探測(cè)技術(shù)的發(fā)明和應(yīng)用，人們對(duì)摻氣水流特性的認(rèn)知在過去幾十年間得到了明顯提升。如美國(guó)的Straub 和Anderson［4］使用電感式探針測(cè)量了明渠自摻氣均勻流中的摻氣濃度分布，并嘗試將均勻流區(qū)的摻氣濃度分布和水深、流量、坡度等參數(shù)關(guān)聯(lián)。Chanson［5］使用電感式雙頭探針在40 kHz 的頻率下探究了渠底法向上的摻氣濃度和流速分布，以及氣泡的弦長(zhǎng)和重現(xiàn)頻率。Boes 和Hager［6，7］使用一種高頻光纖探針對(duì)臺(tái)階式溢流堰上流速、摻氣濃度和氣泡大小等水氣兩相流特性進(jìn)行了探究。Pfister 和Hager［8］使用高速攝像機(jī)拍攝了臺(tái)階式溢流堰初始摻氣點(diǎn)附近的流態(tài)并對(duì)摻氣的發(fā)生和氣泡的輸運(yùn)機(jī)理進(jìn)行了分析。Felder 和Chan‐son［9］測(cè)量了26.6°的臺(tái)階式溢流堰上的湍動(dòng)強(qiáng)度、湍流時(shí)間和空間尺度等湍流參數(shù)。Valero 和Bung［10］在一個(gè)1.74 m 高、坡度為1∶2的溢洪道上測(cè)量了水面上空氣層內(nèi)的流速分布并基于此提出了一個(gè)初始摻氣點(diǎn)的估算公式。湯榮才等［11］提出根據(jù)像素亮度計(jì)算細(xì)觀氣泡運(yùn)動(dòng)的高分辨率光流場(chǎng)的方法，彌補(bǔ)了傳統(tǒng)的基于粒子圖像測(cè)速技術(shù)的不足，加深了人們對(duì)摻氣發(fā)生、氣泡擴(kuò)散等的認(rèn)知。

當(dāng)前，VOF 方法可以較好地解決大尺度自由表面（如未摻氣水面）追蹤的問題；且對(duì)于流道中的高雷諾數(shù)水流，k-ε類雙方程渦黏模型可以較好的反映湍流對(duì)時(shí)均動(dòng)量輸運(yùn)的影響。如關(guān)大瑋和程香菊［12］使用FLOW-3D對(duì)高壩閘孔泄流進(jìn)行了三維數(shù)值模擬，驗(yàn)證了相關(guān)模型的準(zhǔn)確性；董宗師等［13］也使用VOF+RNGk-ε模型對(duì)低Fr數(shù)下非完全寬尾墩的消能機(jī)理進(jìn)行了探究。然而，為了實(shí)現(xiàn)對(duì)水流自摻氣的模擬，仍需附加一些模型來考慮摻氣對(duì)水流的影響。當(dāng)前，考慮兩相貫穿的水氣兩相流的數(shù)值模擬多基于雙歐拉框架，如張?jiān)瑢幍龋?4］采用這類模型對(duì)泄洪霧化的機(jī)理及規(guī)律進(jìn)行了探究，鄧斌等［15］也采用此類模型對(duì)波浪作用下直立結(jié)構(gòu)物附近強(qiáng)湍動(dòng)摻氣流體運(yùn)動(dòng)進(jìn)行了數(shù)值模擬。

泄洪水流多采用混合框架，在此框架下模擬摻氣還需附加如下模型：①表面摻氣模型：即在何種條件下水面會(huì)發(fā)生摻氣、能摻入多少的氣體；②氣泡輸運(yùn)模型：即氣體摻入后如何在水中與水相對(duì)運(yùn)動(dòng)，這就進(jìn)一步涉及到氣泡大小的估算、氣泡與水間曳力的計(jì)算等。由于水中氣泡和表面渦體的尺度（10-4～10-1m）比工程尺度（101～103m）小的多，所以在現(xiàn)有的計(jì)算機(jī)水平下通過網(wǎng)格直接解析氣泡與水的相互作用仍非常困難?？尚械姆桨甘菍⑸鲜鰜喐褡樱幢染W(wǎng)格尺度小）尺度的物理作用與基于網(wǎng)格的物理量聯(lián)系起來以實(shí)現(xiàn)方程的閉合。例如，Ma等［16］提出了一種基于界面法向速度梯度的摻氣模型，并將其在船舶領(lǐng)域多相流CFD 代碼CFDshipM 中進(jìn)行了測(cè)試。在水利水電工程領(lǐng)域，更為適用的則為Hirt［17］在2003年提出的模型，它集成于商業(yè)軟件FLOW-3D 中，在過去幾十年內(nèi)得到了多次完善［18，19］，目前可以考慮摻氣引起的密度變化、氣泡的破碎和聚合以及氣泡所受的浮力和曳力等。Valero 和Garcia-Bartual［20］曾針對(duì)光滑溢流堰校驗(yàn)過此模型，Khosro 和Afshin［21］使用其對(duì)兩種特殊體型的臺(tái)階式溢流堰上的流動(dòng)進(jìn)行過模擬，并發(fā)現(xiàn)計(jì)算的斷面平均摻氣濃度與實(shí)驗(yàn)值僅差5%，但未給出模型的具體細(xì)節(jié)。然而，針對(duì)此模型的詳細(xì)介紹及評(píng)價(jià)仍未見報(bào)道，尤其未見對(duì)其計(jì)算氣泡半徑準(zhǔn)確性的評(píng)價(jià)以及對(duì)曳力系數(shù)和Richardson-Zaki系數(shù)乘數(shù)兩個(gè)關(guān)鍵參數(shù)的敏感性分析。

本文使用FLOW-3D 中的摻氣模型對(duì)臺(tái)階式溢流堰上的水氣兩相流進(jìn)行了二維數(shù)值模擬。所模擬的算例是Pfister 和Hager 開展的臺(tái)階式溢流堰試驗(yàn)［8］，其中包括了流態(tài)、初始摻氣點(diǎn)、渠底法向流速及摻氣濃度分布等數(shù)據(jù)。本文通過試驗(yàn)數(shù)據(jù)與模擬數(shù)據(jù)的對(duì)比，對(duì)模型進(jìn)行了系統(tǒng)的評(píng)價(jià)，并對(duì)若干模型參數(shù)的敏感性進(jìn)行了分析，同時(shí)對(duì)模型的后期修正提供了建議。

1 數(shù)學(xué)模型

FLOW-3D 軟件采用的計(jì)算框架可視為混和框架的一種特殊形式，其采用Tru-VOF 技術(shù)［22］在動(dòng)態(tài)追蹤自由水面的同時(shí)，通過在水面施加合適的邊界條件來避免對(duì)氣相的計(jì)算，從而減少計(jì)算時(shí)間。本計(jì)算中采用了RNGk-ε湍流模型進(jìn)行湍流的模擬。

1.1 帶自由面的單相流模型

當(dāng)前，對(duì)水利水電工程中水流進(jìn)行的數(shù)值模擬大都采用帶自由面的單相流模型，其矢量形式的控制方程為：

方程（1）和（2）即為黏性不可壓縮單流體的N-S 方程組。其中U為速度，ρ為流體密度（常量），p為壓強(qiáng)，τ為黏性力，下標(biāo)i和j代表笛卡爾坐標(biāo)系下x、y、z。

除流速和壓強(qiáng)這兩個(gè)最為重要的流動(dòng)參數(shù)之外，還需要?jiǎng)討B(tài)的追蹤自由面的位置。這一任務(wù)通常由VOF 方程完成，其基本形式為：

式中：α為水相的體積分?jǐn)?shù)，則氣相的為1-α。不難發(fā)現(xiàn)，此方程為一純對(duì)流方程，即可視為速度場(chǎng)對(duì)水面的對(duì)流輸運(yùn)。

1.2 湍流模型

本文使用RNGk-ε湍流模型來模擬湍流對(duì)時(shí)均動(dòng)量輸運(yùn)的影響，其湍動(dòng)能k方程和耗散率ε方程分別為：

式中：Pk為湍動(dòng)能的產(chǎn)生項(xiàng)；為相比標(biāo)準(zhǔn)模型改進(jìn)后的系數(shù)。

Pk和的計(jì)算公式為：

其中Deffk和Deffε一般均采用υ+νt/σ計(jì)算。湍流黏度的計(jì)算公式為模型參數(shù)為Cμ=0.084 5，σk=0.719 42，σε=0.719 42，C1ε=1.42，C2ε=1.68，η0=4.38，β=0.012。值得強(qiáng)調(diào)的是，當(dāng)不開啟摻氣模型時(shí)，流動(dòng)為純水流動(dòng)，此時(shí)方程（4）和（5）中密度為常數(shù)，可以省略。

1.3 表面摻氣模型

FLOW-3D 中的摻氣模型［17］由湍流邊界層發(fā)展至水面導(dǎo)致?lián)綒獾睦碚摪l(fā)展而來［23，24］。其假設(shè)自由表面摻氣由失穩(wěn)力和穩(wěn)定力來控制，其中失穩(wěn)力與湍動(dòng)能線性相關(guān)，而穩(wěn)定力來源于表面張力和重力兩部分，且此兩力均和湍流的長(zhǎng)度尺度密切相關(guān)。當(dāng)失穩(wěn)力大于穩(wěn)定力時(shí)，摻氣即可發(fā)生，其控制方程為［17］：

式中：LT代表湍流長(zhǎng)度尺度；Cμ=0.085 為湍流模型常數(shù)；k和ε分別代表湍動(dòng)能及其耗散率；ρ為流體密度；gn是水面法向的重力分量；σ為表面張力系數(shù)；Sa為網(wǎng)格內(nèi)單位時(shí)間內(nèi)摻入的氣體體積；ka為一可調(diào)節(jié)的比例系數(shù)；As為網(wǎng)格內(nèi)自由水面的面積。

腫瘤組織miR-320a與CYLD mRNA表達(dá)水平呈正相關(guān)(r=0.607，P<0.001)(圖6)。

此模型用來模擬自由水面處氣體摻入的過程，在氣體摻入以后，它在水體內(nèi)的擴(kuò)散和運(yùn)動(dòng)由下文介紹的氣泡輸運(yùn)模型和漂移流模型進(jìn)行控制。

1.4 氣泡輸運(yùn)模型

摻氣發(fā)生后，氣泡被水體進(jìn)行平流輸運(yùn)和湍流擴(kuò)散［25］。此運(yùn)動(dòng)可被描述為氣體體積分?jǐn)?shù)的對(duì)流-擴(kuò)散方程：

式中：c為摻氣濃度；Ua為氣相的運(yùn)動(dòng)速度；Dc為氣濃擴(kuò)散系數(shù)，一般通過一常數(shù)分子擴(kuò)散系數(shù)和湍流黏度計(jì)算；Sa為方程（11）中的摻氣源項(xiàng)；Vc為網(wǎng)格體積。

兩相流的平均密度為：

式中：ρw和ρa(bǔ)分別為水和空氣的密度；ρ為兩相的平均密度。

1.5 漂移流模型

浮力、相間阻力和氣泡之間的相互作用是影響兩相流運(yùn)動(dòng)的最重要因素。這些影響在FLOW-3D 中通過漂移流模型來反映。在此模型中，由于水氣密度差較大，所以兩相間的滑移速度被假定為常量，進(jìn)而水中氣體的運(yùn)動(dòng)方程就被簡(jiǎn)化為：

式中：K為相間阻力系數(shù)；Ur為滑移速度；K可以通過單氣泡的阻力系數(shù)Kp計(jì)算得到：

式中：Ap為氣泡的斷面面積；Ur為滑移速度Ur的模；Cd為用戶定義的阻力系數(shù)；ρc和μc分別為連續(xù)相的密度和動(dòng)力黏度；Vp是單個(gè)氣泡的體積；Rp為氣泡半徑。

式中：kRZ是Richardson-Zaki 系數(shù)乘數(shù)，其本質(zhì)上為一個(gè)對(duì)φ0的線性修正系數(shù)，原始公式和程序推薦的值均為1；φ0為Richard‐son-Zaki系數(shù)，其值取決于氣泡的雷諾數(shù)Rb=dbur/υw（其中db為氣泡半徑，υw為水的運(yùn)動(dòng)黏度）；當(dāng)1＜Rb＜500 時(shí)，φ0=4.45/R0.1b，當(dāng)Rb＞500時(shí)，φ0取2.39。

2 試驗(yàn)簡(jiǎn)況及模擬設(shè)置

Pfister 和Hager 的試驗(yàn)在瑞士蘇黎世聯(lián)邦理工學(xué)院水力學(xué)水文學(xué)及冰川學(xué)研究所進(jìn)行。試驗(yàn)所用的臺(tái)階式泄槽詳見圖1：其寬度為0.5 m，長(zhǎng)度為3.4 m，共有25 級(jí)臺(tái)階，每級(jí)高度均為0.093 m，寬度均為0.078 m，相應(yīng)的堰面坡度為50°；該泄槽首部為標(biāo)準(zhǔn)WES 堰，設(shè)計(jì)單寬流量為0.864 m2/s，設(shè)計(jì)水頭為0.533 m。在數(shù)值模擬中，泄槽下游被水平延長(zhǎng)至末級(jí)臺(tái)階后4.8 m，以使水流在水平段充分發(fā)展并在出口處施加自由出游的邊界條件。本模擬計(jì)算了單寬流量為0.215 m2/s 的工況，其對(duì)應(yīng)的堰頂水頭hc為0.167 m。此工況下，試驗(yàn)中觀察到的初始摻氣點(diǎn)位置在x=1.53 m，也就是大約在第6級(jí)臺(tái)階處。

圖1 試驗(yàn)臺(tái)階泄槽剖面圖Fig.1 Longitudinal view of stepped chute in the experiment

在數(shù)值模擬中，計(jì)算域長(zhǎng)度被定為9.5 m（堰頂上游2 m，下游7.5 m）。上邊界被設(shè)置為相對(duì)壓強(qiáng)邊界，相對(duì)壓強(qiáng)為0。上游水頭為0.227 m，下游設(shè)置為自由出流邊界。臺(tái)階結(jié)構(gòu)由FLOW-3D 的FAVOR 技術(shù)［26］捕捉，且其表面被設(shè)置為無滑移壁面。計(jì)算中共使用三個(gè)網(wǎng)格塊，其中第二塊網(wǎng)格塊包絡(luò)了全部泄槽結(jié)構(gòu)。網(wǎng)格塊1 和3 分別在上、下游，三塊網(wǎng)格塊的位置如圖2所示。在模擬坐標(biāo)系（x＇-z＇）中網(wǎng)格塊1 中的網(wǎng)格尺寸為0.025 m×0.016 m，網(wǎng)格塊2 和3 的尺寸則分別為0.002 m×0.002 m和0.025 m×0.013 m。

圖2 網(wǎng)格劃分、算例設(shè)置和臺(tái)階解析度的示意圖Fig.2 Schematic diagram of mesh，simulation setup and step resolution

模擬中的數(shù)學(xué)模型設(shè)置如下：允許氣泡溢出，即當(dāng)水面處摻氣濃度大于允許的最大摻氣濃度Cm時(shí)直接將其設(shè)置為Cm。氣泡初始大小被設(shè)定為1 mm，臨界韋伯?dāng)?shù)和毛細(xì)管數(shù)分別被設(shè)定為1.6 和1。拖曳阻力系數(shù)和Richardson-Zaki 系數(shù)乘數(shù)被分別設(shè)定為0.5和1。水相的最小和最大體積分?jǐn)?shù)被設(shè)定為0.1和1，即允許的最大摻氣濃度為0.9。空氣的密度和黏度被分別設(shè)定為1.225 kg/m3和1.7 e-5kg/m/s，以環(huán)境溫度為15 ℃考慮。

3 結(jié)果分析

3.1 網(wǎng)格無關(guān)性驗(yàn)證

從式（8）～（11）可知，湍動(dòng)能k和湍動(dòng)耗散率ε控制著表面摻氣的發(fā)生。但在臺(tái)階式斜槽中，由于壁面的位置由FAVOR技術(shù)處理決定，因此臺(tái)階面的實(shí)際位置與網(wǎng)格精度（Cs）相關(guān)，k和ε的值預(yù)計(jì)也將受網(wǎng)格精度影響。因此，在本模擬中，對(duì)網(wǎng)格塊2 使用了3 種尺度正方形網(wǎng)格進(jìn)行了解析，分別為5、3 和2 mm。這3種網(wǎng)格精度下流速模值、湍動(dòng)能和湍流特征長(zhǎng)度如圖3所示。在網(wǎng)格無關(guān)性驗(yàn)證的計(jì)算中，與摻氣相關(guān)的模型被關(guān)閉，以免除它們對(duì)湍流計(jì)算的影響，其中在堰面上的位置被無量綱化為實(shí)驗(yàn)坐標(biāo)系中的x/xi（見圖1，xi為初始摻氣點(diǎn)的x坐標(biāo)，模擬工況中為1.54 m）。

從圖3可見，隨著網(wǎng)格的加密，速度模的收斂性要明顯好于湍動(dòng)能和湍流特征長(zhǎng)度。而且計(jì)算結(jié)果的差別主要體現(xiàn)在斷面底部、臺(tái)階附近。在自由表面處，2 mm網(wǎng)格和3 mm網(wǎng)格所給出的3 個(gè)變量的計(jì)算結(jié)果非常接近，因此這兩種網(wǎng)格下最終得到的初始摻氣點(diǎn)位置也基本相同。整體而言，臺(tái)階泄槽上的流動(dòng)參數(shù)對(duì)網(wǎng)格較為敏感，尤其是湍流參數(shù)。最明顯的差別出現(xiàn)在x/xi=0.57斷面的湍流長(zhǎng)度尺度分布上，此斷面上該值不僅在水流中有明顯差別，在水面的值的差別也非常顯著。分析可知，這種敏感性與FLOW-3D 的FAVOR 解析技術(shù)緊密相關(guān)，因?yàn)镕AVOR 解析出來的固體壁面位置會(huì)直接影響壁面函數(shù)，表現(xiàn)為對(duì)無量綱壁面距離y+的影響。在本次模擬中，2 mm 和3 mm 網(wǎng)格所對(duì)應(yīng)的y+分別為60 和70 左右，處于對(duì)數(shù)律層。3 mm 網(wǎng)格得到xi=1.54 m，2 mm 網(wǎng)格得到xi=1.48 m，且與試驗(yàn)值xi=1.54 m 也很接近。因此，本節(jié)后面的結(jié)果均使用2 mm 網(wǎng)格的算例得出。此外，圖中最靠近壁面處的湍動(dòng)能不為0，這是因?yàn)榱鲃?dòng)雷諾數(shù)很高，壁面邊界層非常薄，在圖中并不可見。

圖3 3種網(wǎng)格下垂直偽底面方向的3個(gè)斷面中流速模、湍動(dòng)能以及湍流特征長(zhǎng)度Fig.3 Velocity magnitude，turbulent kinetic energy，and characteristic length of turbulence at three cross-sections perpendicular to the pseudo-bottom under three grid schemes

3.2 初始摻氣點(diǎn)附近的湍動(dòng)能和摻氣濃度

初始摻氣點(diǎn)位置的準(zhǔn)確預(yù)估對(duì)于準(zhǔn)確模擬臺(tái)階溢流堰上的摻氣水流非常重要。目前，學(xué)者們最廣為接受的摻氣初生理論仍是摻氣點(diǎn)的位置是湍流邊界層與水面的交點(diǎn)［23，27］。在本模型中，水面的失穩(wěn)力與湍動(dòng)能k線性相關(guān)［見式（9）］，因此可以用來分析湍流邊界層的發(fā)展過程。圖4展示了試驗(yàn)中的流態(tài)側(cè)視圖以及初始摻氣點(diǎn)附近湍動(dòng)能和摻氣濃度的分布。在圖4（b）中，初始摻氣點(diǎn)正好是湍動(dòng)能k=0.05 m2/s2的等值線與水面的交點(diǎn)，這與模型基本理論和試驗(yàn)觀測(cè)是基本一致的。

由圖4（c）可見，連續(xù)摻氣的開始位置大約在初始摻氣點(diǎn)一個(gè)臺(tái)階以后，這與流動(dòng)的波動(dòng)有關(guān)。此后，摻氣的空氣被輸運(yùn)到下游并最終在2 個(gè)臺(tái)階后在全斷面中得以分布。然而，這種輸運(yùn)速度似乎比試驗(yàn)測(cè)量值要慢。例如，文獻(xiàn)［7］建議可以使用底摻氣濃度1%作為初始摻氣點(diǎn)的位置。但數(shù)值模擬結(jié)果中初始摻氣點(diǎn)后第三級(jí)臺(tái)階表面的摻氣濃度仍小于0.1%。

圖4 試驗(yàn)流態(tài)及模擬中初始摻氣點(diǎn)附近的湍動(dòng)能和摻氣濃度分布Fig.4 Flow pattern in the experiment and turbulent kinetic energy and air concentration distribution near the inception point

3.3 摻氣濃度分布

在文獻(xiàn)中，斷面摻氣濃度的分布常作為從偽底面起算的無量綱化的摻氣水深z/z90的函數(shù)加以描述，其中z90為從偽底面起算并與之垂直方向上摻氣濃度為90%處的水深。為了將數(shù)值模擬結(jié)果與試驗(yàn)值進(jìn)行比較，數(shù)值模擬中的自由表面高程被視為z90。在本節(jié)中，使用zsur來同時(shí)表征試驗(yàn)和數(shù)值模擬中的摻氣水深。

圖5展示的是不同斷面不同無量綱摻氣水深處的摻氣濃度分布，圖5（b）中水面處還同時(shí)標(biāo)注了水面上計(jì)算得到的摻氣濃度值。從圖5（a）可以發(fā)現(xiàn)，在z/zsur＜0.4 范圍內(nèi)，摻氣濃度隨z/zsur相對(duì)緩慢，而在上半部分，它則快速增長(zhǎng)到0.9。從x/xi=1.37 開始，摻氣濃度的分布曲線幾乎不再變化，預(yù)示著水流達(dá)到均勻流狀態(tài)。均勻流區(qū)域的底摻氣濃度約為0.2 左右。在數(shù)值模擬中大約從x/xi=1.21 開始?xì)馀荼惠斶\(yùn)到整個(gè)斷面。這顯然比x/xi=1.05的實(shí)驗(yàn)值要更靠下游。同時(shí)，表面處的摻氣濃度也并沒有達(dá)到0.9，而是從x/xi=1.21 的0.25 逐漸增加至x/xi=1.85 的0.8。由于實(shí)驗(yàn)中缺乏摻氣水深z90數(shù)據(jù)，當(dāng)前模型計(jì)算摻氣水深的準(zhǔn)確性仍然需要進(jìn)一步研究。從數(shù)值模擬結(jié)果中可以發(fā)現(xiàn)摻氣濃度約從z/z90=0.4 位置開始迅速增加，這與實(shí)驗(yàn)得到的結(jié)果是一致的。然而，數(shù)值模擬中的底摻氣濃度比試驗(yàn)中增加更為迅速，底摻氣濃度的最大值出現(xiàn)在x/xi=2.01 處，其值約為0.4?？傮w來說，F(xiàn)LOW-3D 中的摻氣模型目前對(duì)于臺(tái)階泄槽上初始摻氣點(diǎn)的模擬較為準(zhǔn)確，但對(duì)于摻氣濃度計(jì)算的準(zhǔn)確性仍然欠佳。

圖5 不同斷面的摻氣濃度分布Fig.5 Air concentration distribution at different cross-sections

3.4 流速分布

計(jì)算得到的流速模在若干斷面上的分布如圖6所示，其中，試驗(yàn)的流速模分布曲線以粗黑實(shí)線畫出用以和模擬值對(duì)比。之所可以用一條線描述試驗(yàn)中流速分布是因?yàn)樵囼?yàn)中x/xi=0.65以后所有斷面的流速相對(duì)無量綱水深的分布是幾乎相同的。然而，模擬得到的渠底流速值僅為表面流速的25%到30%，而在試驗(yàn)中渠底流速則約為表面流速的60%。這可能與摻氣對(duì)壁面阻力的減緩作用有關(guān)［28］，因?yàn)檫@在數(shù)值模型中是沒有被考慮在內(nèi)的。顯而易見，這種流速分布和梯度的差異會(huì)顯著影響湍流的統(tǒng)計(jì)參數(shù)，進(jìn)而改變摻氣濃度的值。更具體的說，模擬中得到的速度梯度較大，會(huì)導(dǎo)致湍動(dòng)能和湍流黏度更大，并進(jìn)而加劇摻氣在垂直渠底方面上的擴(kuò)散。因此，額外的修正模型對(duì)于臺(tái)階泄槽上水流的流速和摻氣濃度的準(zhǔn)確模擬就非常重要。

圖6 不同斷面處計(jì)算得到的流速分布曲線與試驗(yàn)對(duì)比Fig.6 Comparison of calculated and experimental velocity distribution at different cross-sections

3.5 氣泡半徑

本模擬中計(jì)算所得的氣泡大小如圖7所示?？梢姡?jì)算所得的氣泡大小約在0～12 mm 之間。渠底的氣泡較小，越靠近水面其尺度越大，且臺(tái)階腔內(nèi)部的氣泡大小還略大于偽底面處的氣泡大小。氣泡大小的變化可以通過壓強(qiáng)來解釋，因?yàn)檩^高的環(huán)境壓強(qiáng)將迫使氣泡破碎分解來產(chǎn)生更高的表面張力來平衡氣泡圍壓。這與試驗(yàn)中的觀察結(jié)果是一致的。文獻(xiàn)［5］發(fā)現(xiàn)氣泡的尺寸范圍很廣，從接近于0 到高達(dá)140 mm，而高概率的氣泡大小在0～3 mm 左右。文獻(xiàn)［29］則發(fā)現(xiàn)預(yù)摻氣的水流中氣泡大小一般在1～40 mm 左右。因此，模擬所得到的氣泡大小是基本合理的，但水面處氣泡尺寸的突降則比較奇怪，這可能與邊界條件有關(guān)。

圖7 氣泡直徑的計(jì)算結(jié)果Fig.7 Calculated air bubble diameter

3.6 阻力系數(shù)和Richardson-Zaki 系數(shù)乘數(shù)的影響

除氣泡尺寸參數(shù)之外，在FLOW-3D 的漂移流模型中還有兩個(gè)非常重要的參數(shù)，即拖曳阻力系數(shù)Cd和Richardson-Zaki系數(shù)乘數(shù)kRZ。本節(jié)對(duì)計(jì)算結(jié)果對(duì)這兩個(gè)參數(shù)的敏感性進(jìn)行了初步分析。除軟件推薦的默認(rèn)值Cd=0.5和kRZ=1之外，Cd也被設(shè)為0.1 和0.9，kRZ也被設(shè)置為了另外兩個(gè)值0.5 和1.5。相應(yīng)的計(jì)算結(jié)果如圖8所示。在z/zsur＜0.5的下部區(qū)域，Cd和kRZ的值對(duì)摻氣濃度的影響很小；但在上部，區(qū)別則較為明顯。但這些差別與模擬值與試驗(yàn)值之間的差別相比是很小的。由于目前的模擬已經(jīng)非常耗時(shí)（在16 核心的Intel Xeon E5-2667 v3，3.2 Ghz處理器上，2 mm 網(wǎng)格下單個(gè)算例的計(jì)算時(shí)間為226 h），筆者并沒有對(duì)這兩個(gè)參數(shù)進(jìn)行進(jìn)一步的研究。同時(shí)，在氣泡尺度沒有被準(zhǔn)確估計(jì)之前，對(duì)這兩個(gè)參數(shù)的校準(zhǔn)意義有限。在后續(xù)研究中，筆者建議先找到一個(gè)可靠性較好的氣泡尺度模型，或者使用氣泡尺度基本恒定的試驗(yàn)來對(duì)Cd和kRZ進(jìn)行校準(zhǔn)。

圖8 不同Cd和Crz下x/xi=1.53斷面的摻氣濃度分布Fig.8 Air concentration distribution at x/xi=1.53 with different Cd and Crz values

4 模型討論與分析

VOF+湍流模型只適用于水面未摻氣情況下的流動(dòng)。摻氣發(fā)生后，VOF 捕捉到的自由面下方的水體內(nèi)含有氣泡，其密度不再為常數(shù)，而是一個(gè)場(chǎng)；同時(shí)摻氣水流的黏度和湍流特性也與純水流有著顯著區(qū)別。從數(shù)學(xué)模型的角度來看，進(jìn)行摻氣模擬的第一步就是附加相關(guān)模型以使氣泡摻入過程被考慮在模擬過程之中。本文研究表明，F(xiàn)LOW-3D 中的摻氣模型可以較好的捕捉臺(tái)階式溢流堰上初始摻氣點(diǎn)的位置，此模型的機(jī)理符合“湍流邊界層發(fā)展至水面引起摻氣發(fā)生”的普遍認(rèn)知，較好的解決了摻氣水流模擬的首要難題。但由于摻入氣體的量在試驗(yàn)中難以測(cè)量，也就難以對(duì)模型在這方面的準(zhǔn)確性直接驗(yàn)證，而只能通過摻氣水深、斷面摻氣濃度分布等其它指標(biāo)來衡量。但這些指標(biāo)的準(zhǔn)確計(jì)算建立在對(duì)氣泡大小和摻氣水流的湍流特性的可靠估計(jì)的基礎(chǔ)之上。因此，如何將氣泡大小、湍流模型、拖曳力模型等對(duì)計(jì)算結(jié)果的影響?yīng)毩㈤_來，并一一測(cè)試、驗(yàn)證和改進(jìn)就是一個(gè)非常重要且不可回避的問題。整體來看，F(xiàn)LOW-3D 中的模型從功能性上來說具備了預(yù)測(cè)氣泡大小、估計(jì)拖曳阻力等的能力，其計(jì)算得到的氣泡尺度從量級(jí)上來說較為合理，但對(duì)摻氣水流相關(guān)模型進(jìn)一步的探究和修正仍依賴于對(duì)模型細(xì)節(jié)的更深入的了解。一個(gè)較為理想的研究思路是找到氣泡大小相對(duì)變化不大的案例，或者通過人為預(yù)設(shè)氣泡大小分布的方式排除此關(guān)鍵參數(shù)對(duì)計(jì)算結(jié)果的影響，同時(shí)通過模型修正摻氣區(qū)域的壁面阻力，來以實(shí)現(xiàn)對(duì)拖曳力模型和湍流模型的可靠性評(píng)估。

5 結(jié) 語

本文對(duì)FLOW-3D 中的亞格子摻氣模型對(duì)于臺(tái)階泄槽上的摻氣水流的適用性進(jìn)行了探究。文中所選取的試驗(yàn)數(shù)據(jù)包括初始摻氣點(diǎn)的位置、摻氣濃度以及速度分布曲線。網(wǎng)格收斂性驗(yàn)證表明湍流參數(shù)對(duì)于網(wǎng)格比較敏感，在最精細(xì)的網(wǎng)格（單元約為1/50 臺(tái)階大?。┲?，湍流參數(shù)的收斂性仍然不理想。最精細(xì)網(wǎng)格的計(jì)算中準(zhǔn)確模擬了非摻氣區(qū)域湍流邊界層的發(fā)展和摻氣初生的位置。然而，計(jì)算所得到的摻氣區(qū)的流速分布和摻氣濃度分布與試驗(yàn)值仍然有著較大的差別。同時(shí)，計(jì)算得到的氣泡大小較為合理。本文還發(fā)現(xiàn)拖曳阻力系數(shù)和Richardson-Zaki 系數(shù)乘數(shù)兩個(gè)用戶自定義參數(shù)在本算例中對(duì)摻氣濃度分布的影響不大。因此，包括調(diào)校準(zhǔn)確氣泡尺度模型、修正摻氣區(qū)域的壁面阻力等進(jìn)一步的改進(jìn)對(duì)于準(zhǔn)確模擬水氣兩相流的摻氣濃度、流速分布等參數(shù)而言非常重要。并且，以上模型的修正和調(diào)校最好均在開源代碼中完成，這也是當(dāng)前基于FLOW-3D 這類閉源商業(yè)軟件無法詳細(xì)獲取相關(guān)模型的全部細(xì)節(jié)并對(duì)其進(jìn)行改進(jìn)的主要原因。