FLOW-3D氣泡流模型關(guān)鍵參數(shù)敏感性探究

董宗師，蔡 芳

（1.長江設(shè)計集團，湖北 武漢 430010；2.四川大學(xué)水利水電學(xué)院，四川 成都 610065）

0 引 言

流動摻氣尤其是水面自摻氣在明渠、溢洪道、水工隧洞等引水、輸水、泄水建筑物中非常常見，對工程的安全運行有著重要的影響［1-3］。具體而言，摻氣會增加流道中的水深［4-6］，影響射流的沖刷效果［7］，還可以影響壩下魚類等水生生物的生存環(huán)境［8-10］；同時，結(jié)構(gòu)面附近的氣泡對減免高速水流的空蝕破壞也有著至關(guān)重要的作用［11，12］。當(dāng)前，對于工程水力學(xué)的研究多依據(jù)重力相似準(zhǔn)則而采用模型試驗的方法，但在摻氣水流中表面張力和黏性力的作用明顯，因此，試驗研究成果與工程實際相比存在較為明顯的比尺效應(yīng)［13］。原型觀測也存在成本高、測量參數(shù)有限、儀器結(jié)構(gòu)強度不足等缺點［14］。而隨著計算機算力的飛速提升和計算物理學(xué)的快速發(fā)展，計算流體力學(xué)（CFD）在工程水力學(xué)研究中得到了越來越廣泛的應(yīng)用。

相比帶有自由面的單相流模擬，摻氣水流的模擬要復(fù)雜的多，需要附加自由水面摻氣、氣泡尺寸計算、水-氣相對速度估計、兩相湍流模擬以及氣泡的變形、破碎及聚合等多種其他模型［15］。目前，Ansys Fluent 和OpenFOAM 等流行的CFD 軟件中雖然含有多個氣∕液、氣∕固、液∕固多相流曳力和湍流模型，但無法考慮界面處的湍流摻氣，也就無法完成對于摻氣水流的模擬［16，17］。近年來，多個學(xué)者使用FLOW-3D 軟件對臺階［18］、寬尾墩［19］、明渠［20］等泄水建筑物中的摻氣水流進行了數(shù)值模擬，開創(chuàng)了在泄洪模擬中考慮摻氣作用的先河。然而，F(xiàn)LOW-3D 軟件官方文檔中并未對以下模型細節(jié)進行說明，限制了用戶對軟件參數(shù)的選擇和準(zhǔn)確性分析。①在計算氣泡直徑時采用的模型方程；②臨界毛細管數(shù)中的參數(shù)含義以及臨界韋伯?dāng)?shù)和臨界毛細管數(shù)對于不同流動的主導(dǎo)作用；③初始氣泡直徑的取值建議以及其通過何種方式影響計算結(jié)果。此外，模型中含有兩個重要的參數(shù)——拖曳系數(shù)以及Richardson-Zaki 調(diào)節(jié)系數(shù)，雖然軟件給出了默認值，但其對計算結(jié)果的影響也不明確。本文對FLOW-3D 摻氣流模型中臨界毛細管數(shù)、臨界韋伯?dāng)?shù)、氣泡初始直徑、拖曳阻力系數(shù)以及Richardson-Zaki 調(diào)節(jié)系數(shù)等模型參數(shù)對計算結(jié)果的影響進行了敏感性分析，以供相關(guān)的數(shù)值模擬參考。

1 研究對象與模型設(shè)置

1.1 研究對象

本文選取了Straub & Anderson 的明槽均勻流試驗［21］中45°坡度下流量為0.269 m3∕s 的工況為驗證算例。試驗中所使用的明槽寬和高分別為0.45 m 和0.3 m，長度為15 m，槽底平均突起高度為0.71 mm。在模擬中，忽略了兩側(cè)邊墻的影響將展向網(wǎng)格設(shè)置為1 進行二維模擬。同時，為了保證水流達到均勻流狀態(tài)，將泄槽長度延長至50 m，并在x=48 m 處進行數(shù)據(jù)采集。計算體型和設(shè)置示意見圖1。

圖1 計算體型與設(shè)置示意圖Fig.1 Schematic diagram of simulated geometry and simulation setup

1.2 模擬設(shè)置

在計算中，為了保證FAVOR 界面捕捉技術(shù)對槽底的解析精度，明槽與x軸平行（即水平）放置，并根據(jù)明槽45°坡度來調(diào)整重力分量，流向和垂直槽底向下的重力加速度均為6.937 m∕s2。計算使用了均勻的結(jié)構(gòu)化網(wǎng)格，單元尺寸為Δx=0.017 m，Δz=0.025 m，并將槽底第一層網(wǎng)格高度設(shè)置為0.75 mm 以保證壁面函數(shù)計算湍流參數(shù)的準(zhǔn)確性（實際計算得到的槽底y+值為97）。計算域邊界條件設(shè)置如下：底部設(shè)置為無滑移固壁邊界，頂部給定零壓邊界，進口根據(jù)流量和相同水力條件下的均勻流水深h0=0.049 m 給定均勻流速邊界，并假定湍流強度為5%，出口則設(shè)置為自由出流。所有計算結(jié)果均在計算域內(nèi)水體的總體積變化率小于0.5%時采集。

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

FLOW-3D 軟件因其特有的Tru-VOF 自由面追蹤方法而在水利水電、港口工程等領(lǐng)域得到了廣泛的應(yīng)用。Tru-VOF 方法在每個時間步完成自由水面追蹤之后，即將空氣區(qū)域移出計算區(qū)域而只對壁面和水面圍成的水體進行計算。這種技術(shù)依賴于對壁面和水面處施加恰當(dāng)?shù)倪吔鐥l件［22］。而對于湍流的模擬，考慮到所模擬流動的高雷諾數(shù)而采用了RNGk-ε湍流模型。對氣泡流的模擬則使用了表面摻氣模型、變密度模型、漂移流模型、氣泡尺度模型等系列模型來實現(xiàn)。上述模型相應(yīng)的方程如下。

2.1 納維-斯托克斯方程組

式中：u為流速；f為體積力；ρ和υ分別為水的密度和運動黏度系數(shù)。

2.2 湍流模型

式中：k為湍動能，ε為湍動耗散率；Pk為由于平均運動梯度引起的湍動能產(chǎn)生項；Deffk和Deffε分別為k和ε的擴散系數(shù)，其計算方法及其他常數(shù)的取值詳見文獻［15］。

2.3 表面追蹤模型

式中：f為水的體積分數(shù)；Sa為摻氣源項［詳見公式（9）］；Vc為網(wǎng)格單元體積。

此方程為一純對流方程，為了避免產(chǎn)生數(shù)值擴散，F(xiàn)LOW-3D采用施主-受主技術(shù)來幾何求解和重構(gòu)界面。

2.4 表面摻氣和逸出模型

FLOW-3D 中的摻氣模型［23］假設(shè)摻氣的發(fā)生是由于湍動能在垂直水面的分量大于表面張力和重力兩個抗力，其控制方程為：

式中：LT為湍流特征長度；gn是重力在水面法向的分量；σ為表面張力系數(shù)，計算時取0.073 N∕m；Sa為摻氣速率；ka為一用戶定義的系數(shù)，默認值為0.5；As為水面處網(wǎng)格內(nèi)水面的面積，湍流模型常數(shù)Cμ=0.085。

FLOW-3D 還允許氣泡從水面逸出，其實現(xiàn)方法為用戶可以定義一個最大摻氣濃度cmax，當(dāng)水面處的摻氣濃度大于cmax時會被重新設(shè)定為cmax。

2.5 變密度模型

變密度模型用以考慮摻氣引起的流體密度變化。摻入的氣泡與水體的相對運動可被抽象為摻氣濃度的對流-擴散方程：

式中：c為摻氣濃度；Ua為氣相的速度；Dc為摻氣擴散系數(shù)；Sa為式（5）和式（9）中的摻氣源項；Vc為各網(wǎng)格單元的總體積。

摻氣水流的宏觀密度則直接計算為體積平均：

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

2.6 漂移流模型

浮力、水-氣間相對運動的阻力和其他相互作用通過漂移流模型來模擬［24，26］。考慮到水體和空氣的密度相差近1 000倍，F(xiàn)LOW-3D 中假定水-氣滑移在極短時間內(nèi)達到平衡，因此相間的滑移速度可通過下式進行計算：

式中：K基于網(wǎng)格單元的阻力系數(shù)；ur為水-氣間的滑移速度；K的值通過單氣泡的阻力系數(shù)Kp來計算：

式中：Ap為單個氣泡的截面積；Ur為ur的模；ρc和μc為連續(xù)相（在本文中為水）的密度和動力黏度；Vp為單個氣泡的體積；Rp為氣泡半徑；Cd為用戶定義的拖曳系數(shù)，默認值為0.5。

當(dāng)摻氣濃度較高時，網(wǎng)格內(nèi)氣泡量較多，此時還需要考慮氣泡間的相互作用對滑移速度的影響，F(xiàn)LOW-3D 采用文獻［25］中的修正公式來考慮這一影響：

式中：CRZ為Richardson-Zaki調(diào)節(jié)系數(shù)，軟件推薦其值為1（即不進行調(diào)節(jié)修正）；ζ為Richardson-Zaki 系數(shù)，其值與水流的流動狀態(tài)有關(guān)，可通過氣泡雷諾數(shù)Reb計算。當(dāng)1 ＜Reb≤ 500 時ζ=4.45∕Re0.1b，當(dāng)Reb＞500時，ζ=2.39。

在計算中，可將氣泡直徑設(shè)為定值或由程序根據(jù)流動狀態(tài)計算。FLOW-3D 中使用臨界韋伯?dāng)?shù)和臨界毛細管數(shù)來估計氣泡的尺寸［26］，計算時采用的方程如下：

式中：μm為混和相的分子動力黏度；μt，m為混和相的湍流動力黏度；由量綱分析可知，e的量綱為［T-1］。但官方文檔和相關(guān)文獻并未對其進行說明?？紤]到此值與湍流相關(guān)，因此在本文中取為湍流時間尺度的倒數(shù)。

如上文所述，模型中有若干用戶自定義參數(shù)。氣泡半徑采用動態(tài)計算，這些參數(shù)的默認值如下：臨界韋伯?dāng)?shù)Wec=1.6，臨界毛細管數(shù)Ca，c=1，拖曳系數(shù)Cd=0.5，Richardson-Zaki 調(diào)節(jié)系數(shù)Crz=1。在計算中，水體的體積分數(shù)范圍設(shè)定為0.1～1，即當(dāng)摻氣濃度高于90%時會被自動設(shè)置為90%?？諝獾拿芏群宛ざ确謩e取為1.225 kg∕m3和1.7×10-5kg∕（m·s），且水相被一直認為連續(xù)相。

3 計算成果與分析

3.1 氣泡尺寸的主要控制參數(shù)

如上文所述，F(xiàn)LOW-3D 中的氣泡尺寸模型主要通過臨界韋伯?dāng)?shù)Wec、臨界毛細管數(shù)Ca，c和氣泡聚合模型來實現(xiàn)，但未給出上述參數(shù)的閾值和模型方程。這些細節(jié)直接影響了氣泡的大小，對于算例的設(shè)置非常重要。雖然在軟件中并未直接提供韋伯?dāng)?shù)和毛細管數(shù)的值，但可以通過一定的方法進行二次計算。在均勻流區(qū)域（即本文的數(shù)據(jù)采集斷面），流向的壓強梯度與垂向相比可以忽略，因此可據(jù)式（12）～（17）計算韋伯?dāng)?shù)和毛線管數(shù)的值：①根據(jù)z 向壓強梯度、摻氣濃度c和氣泡直徑db通過式（12）～（14）求得滑移速度Ur；②根據(jù)式（15）對Ur進行修正得到Ueff r；③根據(jù)Ueff r用式（16）和（17）計算得到We和Ca。圖2展示了用上述方法計算得到的毛細管數(shù)和韋伯?dāng)?shù)，其中藍色實線為軟件默認值。

圖2 x=48 m斷面上的毛細管數(shù)和韋伯?dāng)?shù)垂向分布Fig.2 Vertical distribution of Capillary number and Weber number at x=48 m

由圖2可見計算得到的毛細管數(shù)和韋伯?dāng)?shù)與設(shè)定值均有一定差別。毛細管數(shù)的計算值在水體中部為一恒定值，且僅比設(shè)定值大15%左右，但算得的韋伯?dāng)?shù)與設(shè)定值則有著顯著差別，雖然在水面處較為接近，但在底部相差9 個量級。由此可見，F(xiàn)LOW-3D 的氣泡流模型中對于明渠這種強紊動摻氣水流，氣泡尺寸主要是受臨界毛細管數(shù)控制。從物理作用機理的角度來看，因為湍流對氣泡破碎的作用強于水-氣間的慣性力，因此這一算法是合理的。至于毛細管數(shù)計算值與設(shè)定值的差別，則于滑移速度修正、湍流黏度計算等因素均有一定關(guān)系，由于這些均為軟件的黑匣子細節(jié)，因此難以進一步定論。

3.2 臨界毛細管數(shù)的影響

為了進一步考察FLOW-3D 摻氣水流特性計算結(jié)果對臨界毛細管數(shù)Ca，c的敏感性，分別將Ca，c取為1、2、4 三個不同的值探進行模擬探究。圖3展示了3個Ca，c取值下的湍動能、湍流時間尺度以及湍流黏度的計算結(jié)果。由圖3可見在不同的Ca，c下湍動能、湍流時間尺度以及湍流黏度的垂向分布曲線幾乎完全重合，表明將臨界毛細管數(shù)提高至默認值的4 倍對水流的湍流特性沒有影響。

圖3 3種臨界毛細管數(shù)下的湍動能、湍流時間尺度和湍流黏度Fig.3 Turbulent kinetic energy，turbulent time scale and turbulent viscosity under three critical capillary numbers

圖4展示了Ca，c分別設(shè)定為1、2、4 時計算得到的摻氣濃度、氣泡直徑、滑移速度及毛細管數(shù)。從式（17）來看，對于本文模擬的強烈紊動的明渠水流，相比湍流黏度而言分子黏性可以忽略不計。因此，在同一湍流時間尺度下，氣泡直徑與臨界毛細管數(shù)線性相關(guān)。然而，本文測試的3 個最大相差4 倍的臨界毛細管數(shù)所得到的氣泡的直徑相差不大，同時增加臨界毛細管數(shù)使滑移速度和毛細管數(shù)也小幅增加。這主要是由于氣泡聚合作用強于臨界毛細管數(shù)的影響，使臨界毛細管數(shù)對計算結(jié)果的影響并不明顯。

圖4 3種臨界毛細管數(shù)下的摻氣濃度、氣泡直徑、滑移速度和毛細管數(shù)Fig.4 Air concentration， bubble diameter，slip velocity and capillary number under three critical capillary numbers

3.3 臨界韋伯?dāng)?shù)的影響

本節(jié)對計算結(jié)果對臨界韋伯?dāng)?shù)Wec值的敏感性進行探究。在FLOW-3D 中，Wec的默認值為1.6，本文將其值分別取為3.2和6.4來考察Wec對計算結(jié)果的影響。類似地，從圖5可見，將臨界韋伯?dāng)?shù)Wec取為默認值的2倍和4倍并不會對湍動能、湍流時間尺度以及湍流黏度等湍流特性產(chǎn)生顯著影響。

圖5 3種臨界韋伯?dāng)?shù)下的湍動能、湍流時間尺度和湍流黏度Fig.5 Turbulent kinetic energy，turbulent time scale and turbulent viscosity under three critical Weber numbers

圖6展示了3 種Wec值下的摻氣濃度、氣泡直徑及韋伯?dāng)?shù)的計算結(jié)果?？梢园l(fā)現(xiàn)增大臨界韋伯?dāng)?shù)會略微降低底部摻氣濃度，同時使水面附近摻氣濃度小幅增加，即使摻氣濃度的垂向差別更大。進一步觀察氣泡尺寸的分布可知這是因為較大的臨界韋伯?dāng)?shù)的值使氣泡尺寸分布更廣，將Wec變?yōu)槟J值的2倍和4倍后水面氣泡尺寸從4 cm 變?yōu)? cm 和9 cm。綜上所述，臨界韋伯?dāng)?shù)對計算的影響也較小，但略大于比臨界毛細管數(shù)的影響。

圖6 3種臨界韋伯?dāng)?shù)下的摻氣濃度、氣泡直徑及韋伯?dāng)?shù)Fig.6 Air concentration，bubble diameter and Weber number under three critical Weber numbers

3.4 氣泡初始直徑的影響

在FLOW-3D 軟件中，對摻氣水流進行計算時需在漂移流子模型中指定初始氣泡直徑，但相關(guān)的模型方程［即方程（5）～（17）］中并未出現(xiàn)氣泡初始直徑這一參數(shù)，考慮到氣泡直徑由軟件中相關(guān)模型自動計算，初步判斷此值對計算結(jié)果沒有影響，本節(jié)主要對此進行驗證。類似地，圖7和圖8分別展示了1、5 和10 mm 三種初始氣泡直徑d0下計算的湍流特性和氣泡特性，顯然3種不同的d0下各參數(shù)的計算結(jié)果幾乎是相同的，表明初始氣泡直徑確實對計算結(jié)果沒有任何影響。

圖7 不同初始氣泡直徑d0下的湍動能k、湍流時間尺度Tt以及湍流黏度μt分布Fig.7 Distribution of turbulent kinetic energy k，turbulent time scale Tt and turbulent viscosity μt under different initial bubble diameter d0

圖8 不同初始氣泡直徑d0下的摻氣濃度Cair和氣泡直徑dbFig.8 Calculated air concentration Cair and bubble diameter db under different initial bubble diameter d0

3.5 拖曳系數(shù)的影響

從水-氣細觀物理作用機理上來看，水-氣間的拖曳作用與水流流態(tài)、氣泡形狀及大小等均密切有關(guān)。在FLOW-3D 的摻氣水流模型中，拖曳系數(shù)也對水-氣間的滑移速度和水流中的摻氣濃度有著重要影響，其默認值為0.5。由式（13）可以發(fā)現(xiàn)，F(xiàn)LOW-3D 在計算兩相拖曳力時，除了拖曳系數(shù)線性項CdUr外，還包含了氣泡雷諾數(shù)的相關(guān)項同時，考慮到一般認為在湍流流態(tài)中單獨使用CdUr計算拖曳力時Cd取值在2.67左右，本文將Cd分別取值為0.1、0.5、1、1.5 和2.67 五個值來考察拖曳系數(shù)對計算結(jié)果的影響，結(jié)果如圖9和10所示。

由圖9可以發(fā)現(xiàn)，拖曳系數(shù)Cd對湍流特性的影響較為顯著。隨著Cd增大，摻氣水流的湍動能也有所增加，湍流黏度則表現(xiàn)為在下半部分減小，上半部分增加，但Cd對湍流時間尺度的影響不明顯。而且，將Cd設(shè)為比默認值0.5 更大的值比使用較小的值對結(jié)果的影響要小。從圖10可見，將拖曳系數(shù)由默認值0.5 減小至0.1 時，摻氣水深由0.2 m 降低了30%，底部摻氣濃度由70%下降到了40%；而增大Cd直至2.67對摻氣濃度和氣泡直徑的垂向分布沒有太大影響，只使摻氣濃度和摻氣水深小幅增加，對氣泡直徑的影響也僅體現(xiàn)在上部水體內(nèi)。因此，拖曳系數(shù)對于計算結(jié)果的影響較為顯著，將其默認值定為0.5 是合理的。

圖9 不同拖曳系數(shù)Cd下的湍動能k、湍流時間尺度Tt以及湍流黏度μt分布Fig.9 Distribution of turbulent kinetic energy k，turbulent time scale Tt and turbulent viscosity μt under different drag coefficient Cd

圖10 不同拖曳系數(shù)Cd下的摻氣濃度Cair和氣泡直徑dbFig.10 Calculated air concentration Cair and bubble diameter db under different drag coefficient Cd

3.6 Richardson-Zaki調(diào)節(jié)系數(shù)的影響

FLOW-3D 的氣泡流模型通過一個乘數(shù)Crz對Richardson-Zaki 系數(shù)進行線性調(diào)節(jié)來模擬氣泡聚合對滑移速度的影響。由式（15）可知，Ur的修正項為一底小于1的指數(shù)函數(shù)，因此當(dāng)Crz＞1 時將減小滑移速度，反之亦然。圖11和圖12分別展示了3個不同的Crz值下計算得到的湍流特性和氣泡特性。由圖11和圖12可以發(fā)現(xiàn)計算結(jié)果對Crz較為敏感：小于1的Crz對會降低湍流程度并減弱摻氣，反之亦然；將Richardson-Zaki 調(diào)節(jié)系數(shù)Crz由默認值1降低至0.5對計算結(jié)果有著明顯影響：槽底摻氣濃度從70%下降43%至40%左右，但Crz對氣泡直徑的影響主要體現(xiàn)在上部區(qū)域，這與拖曳系數(shù)Cd對氣泡大小的影響是相似的。相比之下，將Crz提高至1.5 使滑移速度減小，湍流強度升高，進一步增加了摻氣?？傮w而言計算結(jié)果對Richardson-Zaki 調(diào)節(jié)系數(shù)Crz是較為敏感的。

圖11 3種Richardson-Zaki調(diào)節(jié)系數(shù)下的湍動能、湍流時間尺度和湍流黏度Fig.11 Turbulent kinetic energy，turbulent time scale and turbulent viscosity under three Richardson-Zaki coefficient multiplier values

圖12 不同Richardson-Zaki調(diào)節(jié)系數(shù)Crz下的摻氣濃度Cair和氣泡直徑dbFig.12 Calculated air concentration Cair and bubble diameter db under different Richardson-Zaki coefficient multiplier Crz

4 結(jié) 論

當(dāng)前，F(xiàn)LOW-3D 軟件是唯一一個可以同時考慮大尺度界面摻氣和氣泡空間尺度分布的軟件，值得在摻氣水流的模擬中進行推廣。然而，軟件官方文檔未對氣泡尺度模型的細節(jié)進行詳細說明，這不僅限制了對模型誤差的分析，也限制了軟件的普及應(yīng)用。本文通過數(shù)值試驗對相關(guān)模型的黑匣子細節(jié)進行了逆向探究，同時對關(guān)鍵模型參數(shù)進行了敏感性分析。研究表明：對于明渠自摻氣水流，氣泡大小主要受臨界毛細管數(shù)控制，而臨界韋伯?dāng)?shù)僅在摻氣水面附近對模擬有微弱影響。在本文測試的參數(shù)中，拖曳系數(shù)和Richardson-Zaki 調(diào)節(jié)系數(shù)對計算結(jié)果的影響最為顯著，主要體現(xiàn)在摻氣濃度分布和水體上部的氣泡半徑上?？傮w來看，大于默認值的拖曳系數(shù)和Richardson-Zaki 調(diào)節(jié)系數(shù)會提高湍流強度并增加摻氣，且小于默認值的參數(shù)對計算結(jié)果的影響更為顯著。臨界毛細管數(shù)、臨界韋伯?dāng)?shù)及氣泡初始直徑等參數(shù)對計算結(jié)果的影響很小。在不同算例的模型調(diào)校中，應(yīng)主要對拖曳系數(shù)和Richardson-Zaki 調(diào)節(jié)系數(shù)進行調(diào)校，以提高模型計算的準(zhǔn)確性。