楊 晨 ，陳 鑫 ，薛 亮 ，孫 健

（1. 華北水利水電大學(xué) 水利學(xué)院，河南 鄭州 450046；2. 中國(guó)農(nóng)業(yè)大學(xué) 水利與土木工程學(xué)院，北京 100083；3. 中國(guó)石油大學(xué)（北京） 石油工程學(xué)院，北京 102249；4. 清華大學(xué) 水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084）

相對(duì)于歐美國(guó)家來(lái)看，我國(guó)河流含沙量普遍偏高，特別是北方地區(qū)，如黃河流域。另外，我國(guó)有廣大的河口海岸地區(qū)如環(huán)渤海灣、江蘇中部沿長(zhǎng)江口至浙江北部都屬于淤泥質(zhì)海岸，懸移質(zhì)含沙量比較高。由于生產(chǎn)活動(dòng)的需要，我國(guó)許多水利與海洋工程都不可避免地需要直接應(yīng)用高含沙水，例如使用高含沙水吹填島礁或者用于農(nóng)業(yè)灌溉，整治深水航道的疏浚維護(hù)等，即使某些應(yīng)用會(huì)帶來(lái)如增加灌溉成本、損害水力機(jī)械等很多危害，因此研究高含沙條件下的水沙運(yùn)動(dòng)具有重要意義。隨著計(jì)算機(jī)軟硬件技術(shù)的發(fā)展，工程中比較普遍使用ANSYS 軟件模擬高含沙流動(dòng)，以其中的Fluent 通用兩相流模型的應(yīng)用為例，涉及的模型控制方程和紊流模型眾多，邊界條件復(fù)雜，參數(shù)不確定性極大，本文對(duì)這些問(wèn)題進(jìn)行了探討與研究。

1 控制方程

根據(jù)Elghobashi S E[1]的總結(jié)，水沙兩相的耦合模式按照級(jí)別從低到高可以分為單向耦合、雙向耦合和四向耦合三種類(lèi)型。單向耦合模式下，泥沙顆粒對(duì)流體運(yùn)動(dòng)的影響被完全忽略，只需要描述顆粒在紊流中的對(duì)流擴(kuò)散運(yùn)動(dòng)。當(dāng)泥沙體積濃度達(dá)到10-6～10-3時(shí)，需要用雙向耦合模式，此時(shí)考慮了顆粒相對(duì)流體相紊動(dòng)的影響。特別是顆粒自身對(duì)周邊水流繞流渦的影響，以及這種影響對(duì)顆粒運(yùn)動(dòng)的作用。當(dāng)泥沙相濃度更高的時(shí)候，相間的相互作用需要考慮；泥沙顆粒之間的相互碰撞即碰撞效應(yīng)應(yīng)力變得非常重要，并且這種碰撞還強(qiáng)烈影響到自身以及流體的紊動(dòng)特別是繞流渦，以及大量泥沙懸浮后對(duì)更多泥沙懸浮的再抑制的影響，這些因素都考慮的就是四向耦合模式。然而現(xiàn)有的Fluent 軟件多相流模型卻并非針對(duì)高含沙流動(dòng)構(gòu)建，故而需要在運(yùn)用中針對(duì)此問(wèn)題進(jìn)行大量的二次開(kāi)發(fā)，如層移輸沙問(wèn)題、岸灘底層泥沙運(yùn)動(dòng)問(wèn)題和高含沙泥漿抽取及輸送問(wèn)題等。Fluent 軟件包的通用兩相流模型包括VOF 模型、混合模型和歐拉模型。VOF 模型所考慮的問(wèn)題，兩相之間并沒(méi)有充分耦合，而是彼此之間有很明確的連續(xù)整段分界面分層或者自由表面流?；旌夏Ｐ透鶹OF 方法一樣仍然使用單一流體的方法，求解混合物各種物理量的方程，不過(guò)考慮相間的滑移，由于不考慮相間作用力和顆粒粒間應(yīng)力，不適合于高含沙的情況。此處討論的高含沙數(shù)值模擬是體積濃度達(dá)到10-3時(shí)的情形，適合于歐拉模型，是一種四向耦合模式。

紊流時(shí)均的相連續(xù)方程和動(dòng)量方程，軟件中略去了其中紊動(dòng)的時(shí)間導(dǎo)數(shù)項(xiàng)、三階紊動(dòng)項(xiàng)、壓力梯度紊動(dòng)項(xiàng)、應(yīng)力相關(guān)紊動(dòng)項(xiàng)，得到如下控制方程組：

式中：下標(biāo) k 指代某一項(xiàng)，k=f，s；而 f 代表流體相；s 代表泥沙相；ρ 為相密度；t 代表時(shí)間；x 代表空間坐標(biāo)軸；下標(biāo)i,j=1，2 分別代表水平和垂直方向上的坐標(biāo)，遵循愛(ài)因斯坦求和約定；u 為流速；p 為壓力；τ 為應(yīng)力張量；g 為單位質(zhì)量力；F 為相間作用力，它由相界面上應(yīng)力的積分所得。

運(yùn)用Fluent 連續(xù)介質(zhì)模型，需要強(qiáng)調(diào)的是，模型的建立基于雙流體模型連續(xù)介質(zhì)理論，只適用于高含沙濃度的情況，Hsu T 等[2]認(rèn)為是泥沙體積濃度大于10-4。低濃度條件下的連續(xù)性假設(shè)不是客觀成立的，在實(shí)際問(wèn)題中強(qiáng)行將模型運(yùn)用于低濃度問(wèn)題不一定能獲得收斂解，計(jì)算精度極有可能低于混合模型甚至是普通的擴(kuò)散模型。此外，這個(gè)模型得到的只是兩相的濃度和運(yùn)動(dòng)場(chǎng)，并沒(méi)有方程表達(dá)單個(gè)粒子的運(yùn)動(dòng)，所以不能直接從此模型得到泥沙顆粒運(yùn)動(dòng)軌跡。

在一般濃度兩相流問(wèn)題中，相界面不連續(xù)且相間作用力太過(guò)復(fù)雜，很難對(duì)這些基本方程進(jìn)行封閉求解。對(duì)模型各項(xiàng)的物理意義理解不一樣而導(dǎo)致模型封閉不同，從而使得模型的最終表達(dá)式存在差異。Fluent 的紊動(dòng)關(guān)聯(lián)項(xiàng)中，兩相體積分?jǐn)?shù)與兩相速度相關(guān)項(xiàng)可以根據(jù)梯度輸移假定由式（3）給出，雷諾應(yīng)力項(xiàng)可以根據(jù)Boussinesq 假設(shè)由式（4）給出。

式中：κ=νst/σ 為泥沙紊動(dòng)擴(kuò)散系數(shù)，σ 為泥沙Schmidt 數(shù)；νft和 νst分別為流體相和泥沙相紊動(dòng)粘性系數(shù)。

2 基本應(yīng)力項(xiàng)

Fluent 通用多項(xiàng)模型的基礎(chǔ)是分子運(yùn)動(dòng)論，大量使用動(dòng)理學(xué)理論描述泥沙相運(yùn)動(dòng)。Fluent 軟件完整的固體動(dòng)量方程式（2）包含了固體應(yīng)力項(xiàng)τ，只考慮平移和碰撞從顆粒的動(dòng)量交換中產(chǎn)生的剪切和體積應(yīng)力，源于 Chapman S 等[3]、Syamlal M 等[4]、Gidaspow D 等[5]學(xué)者的理論。基礎(chǔ)是分子運(yùn)動(dòng)論，多用動(dòng)理學(xué)模型建模，主要運(yùn)用于氣固兩相流和顆粒流。與顆粒溫度有關(guān)，隨顆粒溫度增加而增大。顆粒溫度物理意義與紊動(dòng)能相同，與顆粒隨機(jī)運(yùn)動(dòng)的動(dòng)能成比例。應(yīng)力粘性包括碰撞粘性（Collisional viscosity）、動(dòng)力粘度（Kinetic viscosity）、體積粘度（Bulk viscosity）和摩擦粘度（Friction viscosity）。

這些源于分子運(yùn)動(dòng)論和動(dòng)理學(xué)理論的應(yīng)力項(xiàng)運(yùn)用于氣固兩相流和顆粒流獲得了較好的效果，但是在高含沙水流中使用仍然有待完善，因?yàn)榧词刮蓜?dòng)微弱甚至無(wú)紊動(dòng)的高含沙流動(dòng)中顆粒的碰撞和摩擦也是存在的。在河流海洋極限含沙區(qū)域，如底床區(qū)域，流動(dòng)微弱，其適用性值得進(jìn)一步商榷。此外體積粘度是顆粒壓縮和擴(kuò)張的抵抗力；摩擦粘度與可壓縮機(jī)制下顆粒流動(dòng)的固體壓力正相關(guān)，而固體壓力隨顆粒溫度或者是紊動(dòng)能增加而增大。水沙兩相流動(dòng)中的假設(shè)是泥沙顆粒不可壓縮，因此在高含沙的計(jì)算中不適合使用這兩項(xiàng)。

Fluent 軟件并不包含顆粒粒間應(yīng)力 (Intergranular stress)，也叫離散應(yīng)力(Dispersive stress)，這個(gè)是高含沙模擬的關(guān)鍵。為了獲得更貼合實(shí)際的高精度模擬結(jié)果，可考慮在二次開(kāi)發(fā)中將此項(xiàng)編譯于源項(xiàng)中。Bagnold R A[6-7]最早對(duì)它進(jìn)行了相關(guān)研究，其后還有 Ahilan R V 等[8]和 Fredsoe 等[9]，其等價(jià)粘性的一個(gè)基本形式如式（5）和式（6）：

式中：λ 是泥沙相的線性濃度；αsm是泥沙相最大體積濃度，一般取值是 0.6～0.65。Liu H 等[10]、Li M等[11]、Bakhtyar R 等[12]在高濃度水沙兩相模型中充分使用并發(fā)展了粒間應(yīng)力，F(xiàn)luent 模擬中可考慮編譯此項(xiàng)。

3 紊流模型

單相流紊流模型僅有速度紊動(dòng)相關(guān)項(xiàng)，形式相對(duì)簡(jiǎn)單；而兩相流動(dòng)量方程由于包含體積分?jǐn)?shù)、壓力、速度的紊動(dòng)相關(guān)項(xiàng)，使用的紊流模型遠(yuǎn)遠(yuǎn)要復(fù)雜?；贔luent 可以使用或者二次開(kāi)發(fā)的兩相紊流模型跟單相紊流一樣有很多種。在模擬高含沙水流時(shí)所有的常規(guī)紊流模型以及較為先進(jìn)的大渦模擬和直接數(shù)值模擬都可以選擇使用。由于紊流理論及先進(jìn)的紊流模型本身也一直在發(fā)展中，綜合考慮計(jì)算成本、穩(wěn)定性、成熟性以及精度，其中最常用的是基于雙方程紊流模型的?；?。無(wú)論使用何種紊流模型，在高含沙狀態(tài)下均需要合理考慮水沙的特性，特別是大量泥沙運(yùn)動(dòng)對(duì)水流的影響、顆粒碰撞的影響和顆粒繞流渦的影響。以經(jīng)典常用的k-ε 模型為基礎(chǔ)，F(xiàn)luent 提供了3 種方法模擬多相流中的紊流。不過(guò)其既有兩相紊流模型主要源于氣固兩相流動(dòng)，有關(guān)固體相的紊動(dòng)理論主要依靠分子運(yùn)動(dòng)論，在水沙兩相流動(dòng)中的適用性值得進(jìn)一步完善。

最完整和容易理解的是每相紊流模型（Turbulence Model for Each Phase）。它為水相和泥沙相都求解一套k 和ε 輸運(yùn)方程。同時(shí)考慮體積分?jǐn)?shù)、壓力和流速的紊動(dòng)后，兩相流動(dòng)方程按雷諾展開(kāi)平均后推導(dǎo)所得到的紊動(dòng)能方程和紊動(dòng)能耗散率項(xiàng)數(shù)非常大[13]，分別具有三十多項(xiàng)和六十多項(xiàng)。當(dāng)紊流傳遞在相間起重要作用時(shí)，這個(gè)紊流模型是合適的選擇。然而此模型需要模化參數(shù)太多，計(jì)算成本巨大并且不易收斂，還有待進(jìn)一步發(fā)展以提高結(jié)果的可靠性，周立行[14]進(jìn)行了大量研究?？紤]到兩相最大速度差一般不會(huì)超過(guò)最大流速2%[11]，在高含沙水流動(dòng)模擬中的應(yīng)用需要慎重考慮。

Fluent 默認(rèn)的是混合紊流模型（Mixture Turbulence Model）。它源于混合動(dòng)量方程的時(shí)均化，僅求解混合物共同的k-ε 方程，一般認(rèn)為它適合于相分離、分層或接近分層的多相流，或者兩相之間的密度比接近1 的情況。不去考慮泥沙相跟水相紊動(dòng)區(qū)別的情形下，使用混合屬性和混合速度捕獲紊流的重要特征是足夠的。比較理想的運(yùn)用情況是泥沙相的濃度不會(huì)太高，在一般的水力機(jī)械問(wèn)題里面應(yīng)當(dāng)是合適的；而在河流海洋具有低床的情形，由于床面附近含沙濃度達(dá)到了極限狀態(tài)，模擬效果會(huì)差很多。

最普通的是分散紊流模型（Dispersed Turbulence Model）。對(duì)水相使用k-ε 方程求解，在此基礎(chǔ)上采用顆粒追隨理論獲得泥沙相紊動(dòng)特征。這種情形下，對(duì)泥沙相隨機(jī)運(yùn)動(dòng)的起支配作用的是水相紊流的影響，顆粒間的碰撞可忽略。此時(shí)泥沙相的紊動(dòng)量根據(jù)水相的平均特征、顆粒弛豫時(shí)間和粒子相互作用時(shí)間的旋渦給出。當(dāng)明顯地水相是主要連續(xù)相而泥沙相是分散稀釋的第二相時(shí)，這個(gè)模型是比較適用的。總體上看水力機(jī)械中很多高含沙水流都是屬于這種情形，再考慮到此模型相對(duì)成熟和完善，計(jì)算成本不大，因此采用此模型用于水力機(jī)械高含沙水流計(jì)算是比較合適的。

關(guān)于紊流模型，還需要注意到高含沙濃度對(duì)水流紊動(dòng)的抑制作用。Elghobashi S E 等[13]最初紊流模型包含體積紊動(dòng)相關(guān)量，但并沒(méi)有具體應(yīng)用于高含沙問(wèn)題，其參數(shù)也沒(méi)有含沙影響。大量?jī)上嗄Ｐ筒恢苯佑?jì)算高濃度處的紊動(dòng)能量和耗散率，而像單相流動(dòng)利用壁函數(shù)在底部一定高度給定高含沙條件下的紊動(dòng)能和紊動(dòng)能耗散率[12，15-16]。

4 邊界條件

水相的邊界條件中，最需要注意的是壁面邊界條件。傳統(tǒng)的壁函數(shù)表達(dá)式僅僅適合于清水狀況，與泥沙（粒徑、濃度）沒(méi)有關(guān)系。當(dāng)大量泥沙存在時(shí)，粗糙高度、粘性底層和過(guò)渡層的厚度都會(huì)發(fā)生較大的變化；流速分布跟清水條件下的分布也有很大不同，如圖 1 所示[17]。其中：uf與式（1）一樣代表水相流速；u*代表摩阻流速；y' 代表相對(duì)水深，y'=0 為水底，y'=1 為水面。而壁面邊界條件的選擇又極大影響紊流模型及動(dòng)量方程的求解精度，然后影響泥沙擴(kuò)散系數(shù)和濃度分布計(jì)算的精度，最終影響的是泥沙輸送率的估計(jì)。關(guān)于這方面的論著已有不少[18-20]，可以考慮將最先進(jìn)的成果用于Fluent 高含沙模擬中，改善數(shù)值結(jié)果。

泥沙相的邊界條件，最大的問(wèn)題是現(xiàn)有的泥沙運(yùn)動(dòng)理論體系都源于恒定流動(dòng)，在考慮時(shí)間變化時(shí)并非總是適用的，極有可能給計(jì)算結(jié)果帶來(lái)巨大誤差。以周期性的流動(dòng)為例，見(jiàn)圖2 二階Stokes 波流動(dòng)過(guò)程中層移輸沙問(wèn)題泥沙參考濃度的變化過(guò)程：點(diǎn)代表實(shí)驗(yàn)值，點(diǎn)畫(huà)線和虛線代表兩個(gè)經(jīng)典恒定流動(dòng)公式，實(shí)線代表非恒定流動(dòng)理論推導(dǎo)的理論解[21]，t 即式（1）的時(shí)間；T 代表二階 Stokes 波流動(dòng)的周期。傳統(tǒng)恒定流動(dòng)的邊界條件不受相位差影響，在流動(dòng)轉(zhuǎn)向時(shí)候t/T～0.42 瞬時(shí)流速很小，根據(jù)經(jīng)典的Engelund F 等[22]和 Zyserman J A 等[23]公式（即圖 2中EF76 和ZF94）獲得的參考濃度很小，而實(shí)際情況下已懸浮泥沙含量來(lái)不及沉降會(huì)保持一定的數(shù)值；不考慮大量已懸浮泥沙對(duì)再懸浮的抑制，兩個(gè)經(jīng)典公式的解隨著流動(dòng)速度增加一直增長(zhǎng)到很大的水平t/T～0.21，與實(shí)驗(yàn)背離。揚(yáng)沙率或者挾沙能力等泥沙參數(shù)的變化也跟參考濃度類(lèi)似。由于現(xiàn)實(shí)中大部分流動(dòng)都是時(shí)間變化的非恒定流動(dòng)，因此可以考慮在運(yùn)用Fluent 軟件時(shí)重新編譯泥沙邊界條件。

5 結(jié)論

圖1 清水和飽和含沙水流速分布

圖2 非恒定流動(dòng)中泥沙濃度邊界條件

高含沙兩相流動(dòng)的精確數(shù)值模擬在我國(guó)水利工程中具有重要的意義。本文針對(duì)較為普遍使用的Fluent 商業(yè)軟件通用兩相模型，在控制方程、基本應(yīng)力項(xiàng)、紊流模型和邊界條件這幾個(gè)要素的比選上展開(kāi)了討論。在此基礎(chǔ)上本文提出了在軟件模擬仿真中可能出現(xiàn)問(wèn)題的原因及改進(jìn)的方向，并期望Fluent 軟件在高含沙流動(dòng)模擬這方面能得到更好的發(fā)展和應(yīng)用。