水平同軸轉柱體間隙區域間水動力數值模擬研究

談雯倩，張東來，邵宇陽*，郝林杰，諸裕良

(1.河海大學 港口海岸與近海工程學院，南京 210098；2.中交上海航道勘察設計研究院有限公司，上海 201208；3.中水北方勘測設計研究有限責任公司，天津 300222)

細顆粒泥沙為主的河口海岸的航道往往會有泥沙淤積礙航的困擾，這一問題會嚴重影響相關流域的生產和經濟發展。其中細顆粒泥沙的絮凝過程直接關系到河口地區的泥沙沉降-淤積狀況[1]，是河口地形演變、航槽清淤治理、維護水體質量等問題中的一個重要的研究對象。前人針對絮凝過程開展了大量研究，其中絮凝實驗常用的裝置有攪動裝置[2-3]、泰勒庫特流裝置[4-5]等，泰勒庫特流裝置因其操作便捷，通過改變同軸轉柱間隙內流體的流速梯度來研究對泥沙絮凝的影響，控制精準，而被廣泛使用。

流體在旋轉的兩個同軸圓柱之間的流動，即泰勒庫特流，Couette[6]最早在流體粘性系數研究實驗中對此展開研究，之后由Taylor[7]進一步研究該流動的穩定性理論，泰勒庫特流是流體力學中一個重要的研究課題。國內外學者關于泰勒庫特流的研究有很多，早期的研究主要是實驗，Cloe[8]研究了不同半徑比下的泰勒渦和波動泰勒渦出現所需的臨界轉速，且得出了產生波動泰勒渦的臨界轉速與內外圓柱的間隙大小成反比的結論。Werely[9]等人試驗研究了僅內圓柱轉動的情況下間隙內泰勒庫特流。Akonur[10]等利用PIV粒子圖像測速法研究了僅內圓柱旋轉的條件下泰勒庫特流的子午面速度，實驗結果發現泰勒渦的強度隨著雷諾數的增加而增強。鮑鋒[11]等人通過流動顯示實驗和PIV粒子圖像測速法對間隙流場進行可視化和定量化的研究，分析流體運動的周期性規律以及間隙內的雷諾數分布。近年來數值模擬逐漸成為常用的研究方法，Hwang[12]等人以Werely[9]的實驗結果為研究基礎，建立數學模型比較了環形間隙內有無軸向流動對流態的影響，驗證了軸向流動對流場的穩定作用。馮俊杰[13]將實驗和數值模擬相結合得出內圓柱轉速、旋轉雷諾數與流場渦形態轉變的關系。Berghout[14]將Monin-Obukhov曲率長度的概念應用于泰勒庫特流的研究，得到了湍流狀態下的平面速度剖面。Cheng[15]等人利用大渦模擬捕捉泰勒庫特流流動，并提出粗糙度模型以作為今后泰勒庫特流設計小尺度均勻表面粗糙度實驗的指導。泰勒庫特流的數值模擬是近年來研究的熱點問題，泰勒庫特流模型在軍工、航空、動力工程、生物工程、水處理、防洪、生態環境保護和膜分離等領域有著廣泛的應用，因此對它的研究具有重要的理論和實踐意義。

綜合目前的研究來看，對泰勒庫特流的研究一般為豎直放置，泥沙在豎直放置的泰勒庫特流裝置中受到重力的作用，會在軸向上存在一個泥沙濃度梯度，有研究證明泥沙濃度亦是絮凝研究的主要影響因素之一[16]。將同軸轉柱水平放置，周向的流動可以使縱向上泥沙攪動均勻，即可有效避免此類影響。目前對水平的泰勒庫特流研究較少，因此本文利用Fluent軟件對三維不可壓縮的、僅內圓柱旋轉的泰勒庫特流流場進行模擬研究，對比前人的實驗結果進行模型驗證，論證數值模擬的可行性，并在此基礎上建立模型模擬水平放置的同軸轉柱體間隙區域間流場，研究其水動力特征，為今后的絮凝實驗研究打好基礎。

1 模型建立

Wereley和Lueptow[17]曾經對泰勒庫特流進行了實驗研究，所使用的裝置軸向長度為41 cm，內圓柱直徑為8.68 cm，外圓柱直徑為10.46 cm，間隙寬度為0.89 cm。裝置內的液體為水、甘油和碘化鈉的混合液，液體密度為1 620 kg/m3，運動粘性系數為3.15×10-6m2/s。在實驗中，僅內圓柱旋轉，通過改變內圓柱轉速，研究泰勒庫特流裝置內的流場變化，得到了間隙流態由層流轉變為層流泰勒渦的臨界雷諾數與雷諾數為124時的層流泰勒渦流場。本文根據其實驗結果建立模型，并進行模型驗證。

模擬泰勒庫特流裝置，選擇初始條件和邊界條件，所建立的泰勒庫特流模型流場由靜止開始啟動，內壁面為旋轉壁面，以設定的轉速均勻轉動，外壁面為靜止壁面。

1.1 網格模型

通過數值模擬計算控制方程確定解時，需要將控制方程離散化。Fluent中采用有限體積法對控制方程組進行空間和時間上的離散，然后使用分離式解法對離散后的方程組進行求解。

本文模擬的泰勒庫特流裝置，計算區域為環形，選擇均勻分布的結構化六面體核心網格進行劃分。模擬同軸圓柱之間的流體流動，需同時考慮徑向、軸向和周向三個方向上的網格數量，本文在三個方向上劃分的網格數量依次為32、128和256個。

在劃分網格時，考慮到邊界層的影響，在內、外壁面和上、下端面共四個壁面內設置了膨脹層，近壁面起始網格高度設置為0.1 mm，膨脹系數為1.1。

1.2 湍流模型

Realizablek-ε兩方程模式適用于旋轉剪切流、含有射流、混合流的自由流動、管道內流動、邊界層流動、帶有分離的流動和有明顯旋流的流動。有對比分析認為Realizablek-ε湍流模型可以很好地模擬泰勒庫特流的湍流狀態[18]。Realizablek-ε兩方程模式的湍流能量輸運方程和能量耗散方程如下

(1)

(2)

式中：右端三項分別代表生成項、耗散項和壁面項。

2 模擬結果驗證

Wereley等人的實驗重點關注了泰特庫特流在雷諾數為103以及124時的流動狀態，本文將參照其實驗結果，進行數值模擬計算。

2.1 臨界雷諾數的確定

同軸圓柱間隙內流體由層流轉變為層流泰勒渦時的雷諾數，稱為臨界雷諾數。Wereley等人的實驗中，當內圓柱的轉速達到0.827 rad/s時，流態發生變化，并得知此時的雷諾數為103。根據此數據，本文模擬雷諾數為103時的流場，發現因端部效應僅兩端出現渦狀，其余中間部分仍處于層流狀態，總體流態基本可以認為是層流。加大轉速，模擬雷諾數為104時的流動狀態，開始出現了層流泰勒渦的特征，有微弱的渦狀，而當雷諾數Re=105時，層流泰勒渦的特征更為明顯，流場出現規則分布的渦狀，間隙內的流動狀態完全由層流轉換為完全的層流泰勒渦。

2.2 臨界雷諾數時流場模擬結果

2.2.1 臨界雷諾數時流態分布

臨界雷諾數子午面速度矢量圖如圖1所示。上實線代表泰勒庫特流裝置內壁，下實線為外壁，內部曲線為速度等值線，等值線之間的差值為0.1倍內壁轉速。

圖1 臨界雷諾數子午面速度矢量圖Fig.1 Velocity vector of meridian plane at critical Reynolds number

圖1-a為Wereley等通過實驗得到的結果，本文模擬所得的臨界雷諾數在104和105附近，則此時的子午面流速矢量圖如圖1-b、1-c所示。

通過圖1-b、圖1-c可以看出，雷諾數為104時，水流泰勒渦狀不明顯，雷諾數為105時，泰勒渦充分發展，流態渦型明顯。

子午面速度矢量圖體現了徑向流動造成的速度變化。圖1-a、圖1-c等值線變形較為一致，內壁處高速流體向外流動，外壁處低速流體向內流動，從而形成渦狀。明顯的流體向外流動的速度大于向內流動，因而向外流動區域的速度等值線的變化大于向內流動區域速度等值線的變化。

2.2.2 臨界雷諾數時徑向流動特性

圖2為子午面間隙中軸線上的徑向速度對比圖。實驗利用PIV粒子圖像測速法測得結果。數值模擬雷諾數為104和105時的徑向速度，其值與實驗結果有一定的差距，進而模擬雷諾數為106時的結果，與實驗結果相似。

圖2 子午面中軸線徑向速度對比圖Fig.2 Comparison of radial velocity on the central axis of meridian plane

由圖可以看出，所有曲線走向相同，近似正弦曲線，實驗結果所得徑向速度大于數值模擬結果。這是由于PIV法測得的是空間上的平均值，而實驗中不可避免的存在振動和二次流的影響，因此在尋求與實驗結果相似的數模結果時，所使用的雷諾數要大于臨界雷諾數。且圖中也可以看出，向外壁的流速要大于向內壁的流速，空間上向內流動的區域寬度要大于向外流動的區域。

2.3 雷諾數Re=124時流場模擬結果

2.3.1Re=124時流態分布

圖3為雷諾數等于124時的子午面速度矢量圖，上實線代表泰勒庫特流裝置內壁，下實線為外壁，內部曲線為速度等值線，等值線之間的差值為0.1倍內壁轉速。

圖3-a為實驗結果、圖3-b為數值模擬結果。由圖可以看出液體總的流動走向與臨界雷諾數時結果相似，近內壁處的高流速流體向外壁流動，近外壁處的較低流速的流體向內壁方向流動。且雷諾數為124時，圖中泰勒渦明顯增強，渦狀流動更為顯著。徑向流速越大，等值線變化幅度越大，故兩個相鄰的泰勒渦之間等值線變化幅度最大。

2.3.2Re=124時徑向流動特性

圖4為雷諾數Re=124時，子午面中軸線上徑向速度對比圖。圖中模擬結果和Wereley實驗的結果擬合較好。

由圖4所示，不論在走向上還是數值大小上，數值模擬結果與實驗結果都較為相似。雷諾數Re=124時，其模擬結果仍滿足由內壁向外壁流動徑向流速最大，但該區域較小。由外壁向內壁流動區域流速較小，但該區域明顯大于向外流動區域。這一結果符合質量守恒定律，向外壁的流速較大則可能是受到了離心力的影響。當雷諾數Re=124時，最大徑向流速約是雷諾數Re=105時的三倍，即隨著雷諾數的增大，子午面徑向流速隨之增大。

3 水平同軸轉柱體間隙區域間流場數值模擬

3.1 水平同軸轉柱體模型建立

為便于今后的制作考慮，本文數值模擬的水平同軸轉柱體外圓柱直徑設為184.1 mm，內圓柱直徑設為165.1 mm，間隙寬度設為9.5 mm，軸向長度設為400 mm。對水平放置的同軸轉柱體進行數值模擬，建立新的數學模型，在徑向、軸向和周向上網格數量確定為32個、128個和256個，內壁面、外壁面和兩個端面共四個面設置膨脹層，近壁面起始網格高度設置為0.1 mm，膨脹系數為1.1。

線型穩定性分析和能量理論常被用來計算臨界泰勒數[19]，由此計算泰勒庫特流裝置的臨界泰勒數，進而轉化為臨界雷諾數。經計算該裝置由層流轉變為層流泰勒渦時的臨界雷諾數為124，流體完全發展為湍流的臨界雷諾數為3 844，并以此來控制流態變化的模擬。

3.2 庫特泊肅葉流模擬結果

當雷諾數低于臨界雷諾數時，裝置中的流場表現為層流狀態，此時稱為庫特泊肅葉流。本文數值模擬雷諾數Re=65時同軸圓柱間的流場，模擬中將內圓柱的轉速設置為0.082 7 rad/s，進行穩態層流的數值模擬，模擬結果子午面速度云圖如圖5。

由圖可見，當內圓柱的轉速較小，即所設雷諾數低于臨界雷諾數時，整個流場簡單而穩定，除了周向，其他方向上的流速基本為零，流體沒有出現擾動，所產生的影響可忽略不計。此時同軸圓柱間隙間的流場呈現出典型的層流特征，滿足庫特泊肅葉流的要求。

3.3 層流泰勒渦模擬結果

將內圓柱轉速提高到0.191 7 rad/s，此時的雷諾數為150，大于所計算出的臨界雷諾數。間隙內的流動狀態隨之發生改變，流體受到擾動影響，由庫特泊肅葉流發展為完全的層流泰勒渦。

3.3.1 層流泰勒渦流動特性

當雷諾數增加至150時，內外圓柱間隙內流場子午面速度云圖結果如圖6所示。

此時的流態發生了變化，整體出現了規則的擾動，速度等值線發生了規則的彎曲變形。由圖可以看出，流場內流速仍有明顯的分層，且較為均衡穩定。間隙內流場表現出層流泰勒渦的所有特征，圖6直觀地展示了層流泰勒渦在間隙內的流動形式。層流泰勒渦在運動時，流場內出現了規則分布的漩渦，但在此雷諾數下，渦狀流速較小，周向流速仍為主要流動。

3.3.2 層流泰勒渦渦場

圖7為雷諾數為150時子午面速度矢量圖。上、下兩條直線分別代表外壁和內壁，軸向上的曲線為速度等值線，每兩條等值線的差值均為0.1倍的內圓柱轉速。在軸向方向上可以看出有三個完整的渦狀，流速等值線隨著流動而發生變形。

如圖所示，相鄰的兩個泰勒渦間速度等值線變形最大，即相鄰的兩個泰勒渦間徑向流速最大，軸向流速最小。向外壁的最大徑向速度要大于向內壁的徑向速度，根據質量守恒定律，向內壁流動的區域大于向外壁流動的區域，即向內壁徑向流動的兩個泰勒渦之間的間距要大于向外壁徑向流動的渦間的間距。

3.4 湍流泰勒渦模擬結果

進一步增大內圓柱的轉速，當內圓柱轉速增至6.390 7 rad/s時，雷諾數增加至5 000，此時間隙內的流動形態為湍流泰勒渦。

3.4.1 湍流泰勒渦流動特性

圖8為流態變為湍流泰勒渦時子午面速度云圖。由圖可見，同軸圓柱間隙內流體出現了規則的擾動，速度更加均勻，不再具有明顯的分層。

子午面速度云圖直觀地展示出間隙內湍流泰勒渦的流動形式。當同軸圓柱間隙內的主流變化為湍流后，軸向上仍然出現了規則分布的泰勒渦，與較低雷諾數下的子午面速度云圖相比，渦的長度明顯變長，間隙內泰勒渦的總數變少。

3.4.2 湍流泰勒渦渦場

圖9為雷諾數為5 000時子午面速度矢量圖，上、下兩條直線分別代表外壁和內壁，軸向上的曲線為速度等值線，每兩條等值線的差值均為0.1倍的內圓柱轉速。圖中包含三個完整渦的軸向部分，流速等值線隨著流動變形明顯。

圖9 雷諾數Re=5 000時子午面速度矢量圖Fig.9 Velocity vector of meridian plane at Re=5 000

如圖所示，徑向速度的最大值出現在相鄰的兩個泰勒渦間的交界區域，在一個泰勒渦的兩側的交界區域上徑向速度方向相反。徑向速度的最小值出現在通過渦心的徑向線上，其值基本為零，而軸向速度則在此處達到了極值。與較低雷諾數下的子午面速度矢量圖相比，二次流發展更加充分，周向速度更加均勻，在間隙內中心線附近區域周向速度大小幾乎相同，在不同軸向位置徑向線上，軸向速度的零點出現在不同的半徑位置上，這說明了在湍流泰勒渦下，紊動更加劇烈。

3.4.3 湍動特性

湍動特性是湍流泰勒渦最基本的特性，本文采用單位質量內的湍動動能kt來描述湍動的大小。

(3)

圖10為不同半徑位置上的湍動動能沿軸向變化圖，軸向上選取170～215 mm段，該段包含了一對泰勒渦，徑向上選取的半徑為84 mm、86 mm、87 mm、88 mm和90 mm位置處進行對比。

圖10 不同半徑位置上湍動動能沿軸向變化圖Fig.10 Axial variation of turbulent kinetic energy at different radii

軸向上190～200 mm處，徑向流動向外，內壁區域的湍動動能達到最大，軸向上170～175 mm和210～215 mm處，徑向流動向內，外壁區域的湍動動能達到最大，且內壁附近的最大湍動動能要高于外壁附近的最大湍動動能。在徑向流動向外的區域，半徑越大的位置，湍動動能越小，而徑向流動向內的區域與之相反，半徑越大的位置處湍動動能越大。在泰勒渦中心部分，除了外壁附近區域，沿徑向位置的湍動動能趨于一致，都達到最小值。

本節模擬水平放置的同軸轉柱體間隙區域間流體流動，模擬較低雷諾數下的庫特泊肅葉流與層流泰勒渦兩種典型流場，以及高雷諾數下的湍流泰勒渦流場，庫特泊肅葉流為周向層流，流場簡單且較為穩定；隨著雷諾數的增大，層流泰勒渦階段，周向速度在徑向上仍具有明顯分層，同軸圓柱間隙內在軸向上出現規則分布的泰勒渦；在湍流泰勒渦階段紊動劇烈，流速更為均勻，不再分層明顯，泰勒渦中心部分，湍動動能基本一致。從總體上看，豎直與水平放置下的流場中，流態變化過程基本一致，模擬結果與葉立[20]等對泰勒庫特流流態過渡變化的數值模擬結果相符。

4 紊動作用下的絮凝研究討論

水流動力環境是細顆粒泥沙絮凝的關鍵影響因素之一，水流的紊動作用能夠帶動和加強泥沙顆粒間的碰撞從而促進絮凝的發生；另一方面，水流的強紊動作用會破壞絮團結構從而抑制絮凝的發生[21]。水流紊動可被看成各種不同尺度的渦旋運動與平均流速的疊加，渦旋運動產生的剪切力和慣性離心力是泥沙顆粒絮凝發生和破壞的主要作用力[22]。泰勒庫特流裝置具有剪切力小、對絮體破壞小、流動與停留時間易調等優點，因而常用于絮凝實驗研究。

柴朝暉[4]利用同軸旋轉雙筒以層流轉變為層流泰勒渦時的臨界速度ω1、流體完全發展為湍流時的臨界速度ω2為分界點，劃分流態及水流強度。實驗結果證明，當水流強度較小、流態以層流為主時，泥沙顆粒運動以絮凝為主；當水流強度較大、水流為湍流時，紊動作用強，阻礙泥沙的絮凝。韓曉婷[23]數值模擬泰勒庫特流裝置并做絮凝試驗，結果表明泰勒渦渦旋結構清晰時絮凝效果最佳。內筒轉速較小時，湍動動能較小，不足以使絮體顆粒發生有效碰撞，隨著轉速的增加，過大的湍動動能則會破壞絮體結構。

本文的數值模擬結果如表1所示，隨著內圓柱轉速的增加，同軸轉柱體內流速、湍動動能逐級增加，流態隨之變化，經歷了由層流到湍流泰勒渦的過渡過程，本文模擬結果整體趨勢與韓曉婷[23]數值模擬結果相似。豎直放置的同軸轉柱體不可避免的存在泥沙濃度梯度，濃度梯度的存在必然會對絮凝產生一定的影響。本文提出將同軸轉柱水平放置，利用周向流消除濃度梯度影響，更直觀研究流態、水流強度的變化對絮凝的影響，因而更適用于絮凝實驗，當然具體的結果還需要后續根據本文數值模擬結果制作物理模型并開展進一步論證研究。

表1 不同轉速下數值模擬結果Tab.1 Numerical simulation results at different speeds

5 結論

(1)數值模擬結果與實驗結果對比，數值模擬得到的結果與實驗所得基本相同，兩者差約為0.9%。考慮到誤差存在的必然性，可以認為本文所建立的數學模型能較好地模擬出泰勒庫特流裝置內的流場特性，其結果較為精確。

(2)泰勒庫特渦流場中，內圓柱附近流體流速較高，外圓柱附近流速較低，高速流體向外圓柱流動，低速流體向內圓柱流動，從而形成渦狀。子午面徑向流速最大值出現在由內圓柱向外圓柱流動的區域。液體由內向外流動的區域較小，由外向內流動的區域較大。

(3)對水平的同軸轉柱間隙區域間流場進行模擬結果表明，隨著內圓柱轉速的提高，雷諾數逐漸提高，當雷諾數大于臨界雷諾數時，同軸圓柱間的流動不再是穩定的層流，在軸向上會出現明顯的均勻分布的泰勒渦，這一現象與豎直放置的泰勒庫特流裝置結果是相似的。

(4)隨著雷諾數的增加，泰勒渦發展逐漸充分，渦狀特征逐漸明顯，水流流態也由層流變為渦流。雷諾數越大，水流渦動越強，徑向流速越大，對流體的影響更為明顯。但不同雷諾數下，徑向流速與轉速相比只占2%～7%，流場中轉向流速仍起主導作用。

(5)徑向流動向外的渦間區域，內圓柱附近湍動動能最大，沿徑向向外半徑越大的位置，湍動動能越小；與之相反，徑向流動向內的渦間區域，外圓柱附近區域湍動動能最大，湍動動能由外向內逐漸減小；渦中心處，湍動動能基本一致。

(6)同軸轉柱研究水流紊動對細顆粒泥沙絮凝的影響，流態為層流時，水流強度較小，泥沙絮凝；流態為湍流時，水流強度較大，泥沙絮凝受阻。水平放置的同軸轉體，能有效避免豎直方向上的泥沙濃度梯度對絮凝的影響，可進一步開展論證研究，為今后水平同軸轉柱體利用于絮凝實驗提供支撐基礎。