基于FLOW-3D的臺階式溢流堰摻氣水流數值模擬

董宗師，蔡 芳

（1.長江設計集團有限公司，武漢 430010；2.四川大學水利水電學院，成都 610065）

0 引 言

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

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

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

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

1 數學模型

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

1.1 帶自由面的單相流模型

當前，對水利水電工程中水流進行的數值模擬大都采用帶自由面的單相流模型，其矢量形式的控制方程為：

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

除流速和壓強這兩個最為重要的流動參數之外，還需要動態的追蹤自由面的位置。這一任務通常由VOF 方程完成，其基本形式為：

式中：α為水相的體積分數，則氣相的為1-α。不難發現，此方程為一純對流方程，即可視為速度場對水面的對流輸運。

1.2 湍流模型

本文使用RNGk-ε湍流模型來模擬湍流對時均動量輸運的影響，其湍動能k方程和耗散率ε方程分別為：

式中：Pk為湍動能的產生項；為相比標準模型改進后的系數。

Pk和的計算公式為：

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

1.3 表面摻氣模型

FLOW-3D 中的摻氣模型［17］由湍流邊界層發展至水面導致摻氣的理論發展而來［23，24］。其假設自由表面摻氣由失穩力和穩定力來控制，其中失穩力與湍動能線性相關，而穩定力來源于表面張力和重力兩部分，且此兩力均和湍流的長度尺度密切相關。當失穩力大于穩定力時，摻氣即可發生，其控制方程為［17］：

式中：LT代表湍流長度尺度；Cμ=0.085 為湍流模型常數；k和ε分別代表湍動能及其耗散率；ρ為流體密度；gn是水面法向的重力分量；σ為表面張力系數；Sa為網格內單位時間內摻入的氣體體積；ka為一可調節的比例系數；As為網格內自由水面的面積。

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

此模型用來模擬自由水面處氣體摻入的過程，在氣體摻入以后，它在水體內的擴散和運動由下文介紹的氣泡輸運模型和漂移流模型進行控制。

1.4 氣泡輸運模型

摻氣發生后，氣泡被水體進行平流輸運和湍流擴散［25］。此運動可被描述為氣體體積分數的對流-擴散方程：

式中：c為摻氣濃度；Ua為氣相的運動速度；Dc為氣濃擴散系數，一般通過一常數分子擴散系數和湍流黏度計算；Sa為方程（11）中的摻氣源項；Vc為網格體積。

兩相流的平均密度為：

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

1.5 漂移流模型

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

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

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

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

2 試驗簡況及模擬設置

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

圖1 試驗臺階泄槽剖面圖Fig.1 Longitudinal view of stepped chute in the experiment

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

圖2 網格劃分、算例設置和臺階解析度的示意圖Fig.2 Schematic diagram of mesh，simulation setup and step resolution

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

3 結果分析

3.1 網格無關性驗證

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

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

圖3 3種網格下垂直偽底面方向的3個斷面中流速模、湍動能以及湍流特征長度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 初始摻氣點附近的湍動能和摻氣濃度

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

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

圖4 試驗流態及模擬中初始摻氣點附近的湍動能和摻氣濃度分布Fig.4 Flow pattern in the experiment and turbulent kinetic energy and air concentration distribution near the inception point

3.3 摻氣濃度分布

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

圖5展示的是不同斷面不同無量綱摻氣水深處的摻氣濃度分布，圖5（b）中水面處還同時標注了水面上計算得到的摻氣濃度值。從圖5（a）可以發現，在z/zsur＜0.4 范圍內，摻氣濃度隨z/zsur相對緩慢，而在上半部分，它則快速增長到0.9。從x/xi=1.37 開始，摻氣濃度的分布曲線幾乎不再變化，預示著水流達到均勻流狀態。均勻流區域的底摻氣濃度約為0.2 左右。在數值模擬中大約從x/xi=1.21 開始氣泡被輸運到整個斷面。這顯然比x/xi=1.05的實驗值要更靠下游。同時，表面處的摻氣濃度也并沒有達到0.9，而是從x/xi=1.21 的0.25 逐漸增加至x/xi=1.85 的0.8。由于實驗中缺乏摻氣水深z90數據，當前模型計算摻氣水深的準確性仍然需要進一步研究。從數值模擬結果中可以發現摻氣濃度約從z/z90=0.4 位置開始迅速增加，這與實驗得到的結果是一致的。然而，數值模擬中的底摻氣濃度比試驗中增加更為迅速，底摻氣濃度的最大值出現在x/xi=2.01 處，其值約為0.4。總體來說，FLOW-3D 中的摻氣模型目前對于臺階泄槽上初始摻氣點的模擬較為準確，但對于摻氣濃度計算的準確性仍然欠佳。

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

3.4 流速分布

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

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

3.5 氣泡半徑

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

圖7 氣泡直徑的計算結果Fig.7 Calculated air bubble diameter

3.6 阻力系數和Richardson-Zaki 系數乘數的影響

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

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

5 結 語

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