楊航，楊曉勇，白志山，汪營磊，高福磊

（1 華東理工大學機械與動力工程學院，上海 200237；2 西安近代化學研究所，陜西 西安 710000）

泰勒反應器是兩同軸圓筒相對旋轉(zhuǎn)而引起攪拌作用的一類反應器。Taylor[1]于1923 年最早描述了兩個同軸相對旋轉(zhuǎn)圓柱間的泰勒-庫埃特流動。泰勒反應器最早出現(xiàn)于20 世紀70 年代[2]，其具備反應推動力高、比表面積大、傳質(zhì)系數(shù)高和剪切力低等優(yōu)良性能[3]。近年來，泰勒反應器在液-液混合[4]、顆粒制備[5-6]、光催化降解[7-8]以及生物細胞培養(yǎng)[9-10]等許多工程領(lǐng)域得到研究和應用，具有良好的發(fā)展趨勢。

在泰勒反應器中，隨著內(nèi)圓筒轉(zhuǎn)速的增加，可以識別到一系列流型，包括庫埃特流、泰勒渦流、波狀渦流、隨機波狀渦流和湍流泰勒渦流等[11]。在這些流型中，波狀渦流因存在周向波和渦旋間流體輸運而受到了人們的廣泛關(guān)注。Coles[12]首先利用流動可視化技術(shù)研究了波狀渦及其狀態(tài)的多重性，把波狀渦流描述為“a pattern of travelling waves”。隨后，Davey 等[13]、Marcus[14]和Martinand 等[15]分 別 解釋了泰勒渦向波狀渦轉(zhuǎn)變的原因。Bust等[16]研究了轉(zhuǎn)動雷諾數(shù)對波狀渦流的波高和波長的影響，發(fā)現(xiàn)波高和波長都隨著轉(zhuǎn)動雷諾數(shù)的增加而增加。Lueptow 等[17-18]分別使用粒子圖像測速技術(shù)（PIV）來測量環(huán)隙子午面和周向橫截面上波狀渦流的速度場。Kádár和Balan[19]研究了周向波振幅的瞬態(tài)行為以及波狀渦流的波頻率。

上述研究大多數(shù)是在沒有施加軸向流的情況下進行的，然而泰勒-庫埃特反應器總是需要連續(xù)運行，因此研究具有軸流的泰勒庫埃特流對于泰勒反應器的工業(yè)應用具有重要意義。具有軸向流的泰勒-庫埃特流也被稱為泰勒-庫埃特-泊肅葉流[20]，這種流動已經(jīng)引起了越來越多研究者的關(guān)注[21-25]。根據(jù)線性理論[26]，隨著軸向雷諾數(shù)的增加，從庫埃特流轉(zhuǎn)變到泰勒渦流的臨界泰勒數(shù)會變大，因此，軸向流對穩(wěn)定泰勒庫埃特流具有積極作用[27]。Wereley 和Lueptow[28]采用PIV 研究了含軸向流動的波狀渦流在環(huán)隙子午面上的速度場。Hwang 和Yang[29]采用數(shù)值計算的方法研究了子午面上含軸向流動的波狀渦流場，這兩項研究得到了不同泰勒數(shù)和雷諾數(shù)時的速度場，但并未詳細分析周向波動和軸向流對泰勒渦的具體影響。

本文采用計算流體力學（CFD）計算連續(xù)泰勒反應器中的三維瞬態(tài)波狀渦流，并與Wereley 和Lueptow[17,28]的PIV 實驗結(jié)果進行比較，驗證了數(shù)值模擬結(jié)果。研究了周向波動對波狀渦流的特定影響，包括速度分量和渦量的周期性變化，分析了環(huán)隙中不同點的瞬態(tài)行為。探究了軸向流對泰勒渦特性的影響，包括速度分量、渦量和瞬態(tài)行為等。所得的結(jié)果可以為工業(yè)應用中能夠更好地控制泰勒反應器中的相關(guān)流場提供基礎(chǔ)。

1 數(shù)值模擬

1.1 數(shù)學模型

本文主要研究環(huán)隙中的波狀渦流，流體視作不可壓縮的牛頓流體。模型的連續(xù)性方程和動量守恒方程見式(1)、式(2)。

式中，ρ、v、μ、p分別為流體的密度、速度、動力黏度和壓力；g和t分別為重力加速度和時間。

1.2 物理模型與網(wǎng)格劃分

泰勒反應器幾何結(jié)構(gòu)如圖1 所示，其類似于Wereley和Lueptow[28]研究中的結(jié)構(gòu)。流體域內(nèi)徑Ri=43.4mm，外徑Ro=52.3mm，反應器高度L=410mm，間隙寬度d=Ro-Ri=8.9mm，半徑比η=Ri/Ro=0.83，縱橫比Γ=L/d=46.07。標記了幾何結(jié)構(gòu)內(nèi)環(huán)的中心柱面、子午面和緯向環(huán)面。重力加速度為9.81m/s2。從上至下施加低Re的軸向流到環(huán)隙區(qū)域中，其值低于波狀渦流到螺旋波狀渦流的轉(zhuǎn)變雷諾數(shù)。

圖1 泰勒反應器結(jié)構(gòu)

采用軟件ICEM-CFD進行非均勻結(jié)構(gòu)化六面體網(wǎng)格劃分，圖2顯示了模型的網(wǎng)格細節(jié)，由于近壁區(qū)域的速度梯度很高，故對圓筒內(nèi)外壁附近的網(wǎng)格進行了優(yōu)化。

圖2 模型網(wǎng)格劃分和局部放大

對泰勒反應器劃分4套網(wǎng)格，網(wǎng)格方案見表1。采用沿穿過渦旋中心的徑向線的切向速度分布作為評價網(wǎng)格無關(guān)性的標準。切向速度分布、軸向速度分布、徑向速度分布均由內(nèi)圓柱的表面速度RiΩi歸一化。其中，切向速度分布和軸向速度分布均沿著穿過渦旋中心的徑向線，徑向速度分布沿著通過環(huán)隙中心的軸線。如圖3 所示，隨著網(wǎng)格數(shù)量的增加，切向速度分布達到漸近值。網(wǎng)格3與網(wǎng)格4的計算結(jié)果相近，考慮到計算機配置及計算效率的情況下，作為網(wǎng)格獨立性驗證的結(jié)果，選取網(wǎng)格3作為本文的網(wǎng)格方案。

表1 網(wǎng)格方案

圖3 切向速度分布

1.3 條件設(shè)定

泰勒反應器的流動不穩(wěn)定性由量綱為1數(shù)泰勒數(shù)Ta和雷諾數(shù)Re決定。泰勒數(shù)Ta的物理意義為離心力與黏性力之比，如式(3)。

式中，Ri為內(nèi)圓筒半徑；Ω為內(nèi)圓筒轉(zhuǎn)速；d為間隙寬度；υ為流體的運動黏度。

雷諾數(shù)的定義為慣性力與黏性力之比，如式(4)。

式中，w為流體通過環(huán)形通道的平均軸向速度。

本文采用商業(yè)軟件ANSYS-Fluent 14.5對泰勒-庫埃特反應器進行CFD模擬。從頂部施加軸向流進入環(huán)隙區(qū)域，因此設(shè)置模型的頂部平面為速度入口邊界條件，設(shè)置模型的底部出口為壓力出口邊界條件。內(nèi)圓柱壁具有固定的轉(zhuǎn)速，不同的轉(zhuǎn)速對應于特定的Ta。設(shè)置流體域的所有壁面為無滑移邊界條件，設(shè)置靜態(tài)流場為所有情況的初始條件。采用有限體積法離散控制方程，采用三階空間離散格式處理對流項，采用二階中心差分格式處理其他項，采用二階Crank-Nicolson 隱式格式處理時間積分，采用PISO算法處理壓力-速度耦合。時間步長設(shè)置為3×10-2，對應于內(nèi)圓柱旋轉(zhuǎn)2°的時間，方程的收斂準則設(shè)置為10-6。為了獲得時間平均后的結(jié)果，求解進行了20000個時間步長，對應于內(nèi)筒旋轉(zhuǎn)110圈。

2 結(jié)果與討論

2.1 無軸流時的波狀渦流

眾所周知，泰勒渦流的速度場是軸對稱的，而且在周向上沒有速度梯度。但是對于波狀渦流，其速度場比泰勒渦流要復雜得多，并且隨周向位置變化。圖4 顯示了在Ta=139 時不同周向位置的子午面上的速度場和渦量場，渦量ω表征渦流的強度和方向，其表達式如式(5)。

圖4 不同周向位置子午面上的速度場和渦量場

式中，Vr和Va分別表示徑向和軸向速度，渦量的正負值表示漩渦的方向。波狀渦流在周向上有兩個周向波，因此只需要分析從0 到π 上的5 個子午面上的流場。曲線上標出了一對相鄰且旋轉(zhuǎn)方向相反的渦流V1 和V2。可以看出，V1 和V2 的形狀和位置隨周向位置而變化。渦對的形狀左右變形和扭曲，并且位置上下波動。連接渦旋中心的曲線類似于波浪線。在θ=π處，渦對的形狀和位置時又回到了其原始形狀和位置（θ=0 處）。從渦量的角度來看，渦的強度由于渦的變形而改變，波峰（θ=π/4）和波谷（3π/4）時的渦強度比其他位置的大。

圖5(a)為不同子午面上的速度分量分布。由于θ=π處的速度場與θ=0處的速度場相同，因此僅需分析θ=0、π/4、π/2 和3π/4 處的速度分布。圖5(a)展示了沿著線AA′（AA′位置見圖4）的徑向速度分布。對于徑向速度曲線，徑向速度Vr＞0的部分表示流出，Vr＜0的部分表示流入，徑向速度最大值Vr,max的軸向位置表示兩個渦旋的中心邊界?？梢钥闯?，渦對的最大徑向速度和中心邊界隨周向位置而變化。在θ=0時，Vr,max等于0.069RiΩi，相應的軸向位置為0.94d＞0.5λ（λ=1.75d），表明渦流V1被拉伸，渦流V2 被壓縮。在θ=π/4 時，Vr,max增加到0.089RiΩi，中心邊界在0.85d，接近0.5λ。在θ＝π/2處，Vr,max減小至0.071RiΩi，中心邊界在0.77d＜0.85d，表明渦流V1被壓縮而渦流V2被拉伸。θ=3π/4處的Vr,max和中心邊界與θ=π/4 處的相似，并且θ=π 處的Vr,max和中心邊界與θ=0處的重合?？梢缘贸鼋Y(jié)論，周向波的存在導致Vr的周期性變化和中心邊界的振蕩。

圖5 速度分量分布

圖5(b)顯示了沿線BB′（渦V2）和線CC′（渦V1）的軸向速度分布（BB′和CC′的位置見圖4）。對于軸向速度曲線，軸向速度Va＞0 的部分為上升流，Va＜0部分為下降流，Va=0的位置表示渦流的中心。根據(jù)Deng 等[30]的研究，泰勒渦流的渦流對的兩條軸向速度曲線關(guān)于Va=0 幾乎是對稱的。但是對于波狀渦流，由于渦流對的變形，這種對稱性消失了。對于不同的周向位置，軸向速度也顯示出振蕩特性。例如，渦流V1 軸向速度最大值Va,max=0.056RiΩi處的上升流強于軸向速度為Va,min=-0.042RiΩi處的下降流，并且渦流中心位于間隙寬度的中心位置（0.5d）。在θ=π/4時，上升流（Va,max=0.056RiΩi）和下降流（Va,min=-0.042RiΩi）得到了增強。同時，渦流中心（0.53d處）向外圓筒移動。Va,max和Va,min分別增加到0.074RiΩi和-0.069RiΩi，并且渦流中心在θ=π/2時回到0.5d。之后，在θ＝3π/4處，Va,max=0.068RiΩi并且Va,min=-0.065RiΩi，上升流和下降流開始減弱。相應的渦流中心（0.46d處）向內(nèi)圓柱移動。最終，在θ=π處的Va,max和Va,min與在θ=0處的值重復，渦旋中心再次回到0.5d。從模擬結(jié)果可以得出結(jié)論，Va周期性變化，渦流中心在0.5d附近振蕩。

圖6顯示了不同周向位置渦流的最大渦量。由于速度分量的振蕩，渦量在周向上也顯示出波動。如圖4 所示，渦流V2 在θ=0 或π 處渦流變形最大，但渦量最小。在θ=π/4 或3π/4 處，渦流變形最小，但渦量最大，比θ=0處的大24%。渦流V1在θ=π/2處渦流變形最大，但渦量最小。在θ=π/4或3π/4處，渦流變形最小，但渦量最大，比θ=π/2處的大16%，說明渦流V1渦量的變化范圍比渦流V2的小。

圖6 不同周向位置渦流的最大渦量

圖7 為環(huán)隙的一部分中心柱面上的切向速度。為了方便觀察周向波動，將笛卡爾坐標系轉(zhuǎn)換為圓柱坐標。當Ta=131 時，切向速度隨周向位置上下波動，稱為周向波。當Ta=131 和139 時都有兩個周向波，但在Ta=139時的波動比Ta=131時的波更加強烈。當Ta=202時，波數(shù)增加到3個。

如圖8 所示，為了將模擬與實驗結(jié)果進行比較，將衰減泰勒數(shù)ε［定義為(Ta/Tac)-1］作為另一個橫坐標，其中，Tac是從庫埃特流轉(zhuǎn)變?yōu)樘├諟u流的臨界泰勒數(shù)。對于半徑比η=0.83，臨界泰勒數(shù)Tac等于102[31]。波高度H（定義為波的波谷與波峰之間的垂直距離，如圖7所示）是波能量的典型代表。Bust 等[16]發(fā)現(xiàn)，對于η=0.88的情況，波高隨著泰勒數(shù)的增加而增加，直到1.8Tac。本文的模擬結(jié)果也發(fā)生了類似的趨勢。當泰勒數(shù)＜180（≈1.75Tac）時，具有兩個波，并且波高隨著泰勒數(shù)的增加而增加。在Ta=180～202，波高突然下降，同時波數(shù)增加。然后在Ta=253 時，波高降低同時伴隨著波數(shù)增加。同樣，在相似的幾何條件下，模擬的結(jié)果近似于Wereley 和Lueptow[17]的PIV 測量結(jié)果，這也體現(xiàn)了數(shù)值模擬的準確性。

圖7 不同Ta時部分中心柱面上的切向速度

圖8 不同Ta時中心柱面上的波高和波數(shù)

參照圖4，為了量化渦變形的程度，本文定義了一個參數(shù)ξ，如式(6)。

式中，點O為AA'與渦對的中心邊界的交點，ξ＞0 表示渦V1 被拉伸渦V2 被壓縮，ξ＜0 表示渦V1被壓縮渦V2被拉伸。ξ的值越大，表示渦V1被拉伸的程度越大，渦V2被壓縮的程度也就越大。圖9為波數(shù)m=2 時波高對渦變形程度的影響，可知渦V1隨著周向波波高的增加被拉伸的程度越大，渦V2隨著周向波波高的增加被壓縮的程度越大。而且隨著波高的增加渦變形的程度增加得越快。

圖9 波高對渦變形程度的影響

圖10顯示了在Ta=139時不同軸向位置的緯向環(huán)面上的切向速度，3個緯向環(huán)面的位置如圖7(b)所示。對于恒定的θ，切向速度在內(nèi)筒外壁處最大，隨著徑向位置的增加而減小，在外筒內(nèi)壁處減小到零。眾所周知，泰勒渦流的切向速度不會隨周向位置而變化。對于波狀渦流，可以在圖中觀察到切線速度分布隨周向位置的波動。在軸向位置ζ=22.45d處，對應于圖7(b)的波峰，沿周向有兩個波。切向速度在R=0.5d處比在其他位置處的切向速度波動更大，并且在內(nèi)筒外壁或外筒內(nèi)壁附近，波動幾乎消失。在ζ=21.85d處，可以看到在0～2π 之間存在4個波。在對應于波谷的ζ=21.24d處，沿周向也有兩個波，但是波相位與環(huán)面1相差π。

圖10 不同軸向位置緯向環(huán)面上的切向速度

根據(jù)以上內(nèi)容可知，波狀渦流的流場隨周向位置而變化且不穩(wěn)定。因此，泰勒反應器中一點的速度是瞬態(tài)的，研究波狀渦流的瞬態(tài)行為有助于為泰勒反應器的發(fā)展和工業(yè)應用提供有益的見解。由于最大的波動出現(xiàn)在環(huán)隙中心，因此在中心柱面上設(shè)置了3 個監(jiān)視點以檢測切向速度的瞬態(tài)行為。監(jiān)測點的軸向位置分別位于波峰、中點和波谷處。圖11(a)給出了監(jiān)測點處切向速度的時域圖?？梢钥闯?，切向速度隨時間周期性地波動。根據(jù)上述結(jié)論，在Ta＝139處存在兩個從0到2π的周向波，因此，中心表面的旋轉(zhuǎn)周期T為25.64s。在監(jiān)測點1 處，切向速度的周期為T1=12.82s，是監(jiān)測點2 處的T2=6.41s 的兩倍。在監(jiān)測點3 處，周期T3也是12.82s。圖11(b)給出了通過快速傅里葉變換得出的監(jiān)測點的頻域圖，在監(jiān)測點1和3處，主頻率和振幅相同，分別為0.0158Hz和0.0105。主頻為對應于兩倍的周期T1時的頻率0.078Hz。在監(jiān)測點2 上，由于在周期T內(nèi)波數(shù)加倍，主頻為0.0316Hz，其幅度為0.00443，比監(jiān)測點1或3的幅度小57.9%。

圖11 切向速度波動波譜

2.2 有軸流時的波狀渦流

圖12 為Ta=139 和Re=5 時子午面上的速度矢量，其中圖12(a)是數(shù)值模擬的結(jié)果，圖12(b)是Wereley 和Lueptow[28]的PIV 測量結(jié)果，可以看出，計算和實驗結(jié)果吻合良好，這表明計算方法同樣適用于具有軸向流的波狀渦流。

圖12 速度矢量圖

為了進一步驗證數(shù)值模擬結(jié)果的準確性，如表2所示，定量地比較了幾種流動條件下的渦流對波長λ的數(shù)值結(jié)果與Wereley 和Lueptow[28]的PIV 測量結(jié)果，發(fā)現(xiàn)數(shù)值模擬與實驗值之間波長值的最大差小于2%，證明數(shù)值模擬結(jié)果與實驗結(jié)果吻合很好。

表2 不同流動條件下渦對的波長λ

圖13為Ta=139時不同周向位置的子午面上的速度場和渦量場。與圖4 中Re=0 時的情況相比，可以發(fā)現(xiàn)當θ=0時渦的下降流明顯增強，軸向流繞渦運動。同時，由于軸向流的存在，渦心從環(huán)隙的中心向壁面移動，其中渦V1 向內(nèi)圓筒壁面移動，V2 向外圓筒壁面移動。渦的形狀和位置隨著周向位置發(fā)生變化，但是變形的程度和渦心振蕩的幅值均小于無軸向流時。此外，最大渦量的位置也發(fā)生了改變，由渦心附近變?yōu)榱谁h(huán)隙壁面附近。將Re=2.5與Re=5時的速度場和渦量場對比，可以發(fā)現(xiàn)隨著軸向流的增強，渦心的位置向邊壁移動得越厲害，最大渦量的位置也進一步向邊壁移動。Re=5時渦的形狀隨著周向位置的變化比Re=2.5 時的小。這是由于軸向流對于流場具有穩(wěn)定作用。

圖13 不同周向位置子午面上的速度場和渦量場（Ta=139）

圖14 為Re=0 和5 兩種情況下速度分量分布的對比，本文僅分析θ=π/4子午面上的速度分布便可獲取軸流的影響。圖14(a)為徑向速度分布的定量對比，兩條徑向速度曲線非常接近，說明軸流對徑向速度的影響很小。圖14(b)為軸向速度分布的定量對比，可以明顯看出，渦V1和V2的上升流都受到軸流的抑制，而下降流得到了顯著增強。渦V1的中心遷移到了0.35d，渦V2的中心遷移到了0.64d。在圖中，黑色的虛線代表了渦V1 和渦V2 在Re=0和5 時的軸向速度差值Va1-Va2，而紅色的虛線表示層流泊肅葉流動的速度分布，它們基本重合，說明有軸流時泰勒-庫埃特流動的軸向速度為無軸流時泰勒-庫埃特流動的軸向速度與層流泊肅葉流動的軸向速度的疊加。

圖14 速度分量分布

圖15 為Re=5 與Re=0 時 順時針渦V2 在不同周向位置的最大渦量的比較。當Re=5時，渦的最大渦量隨周向位置變化，但其變化范圍僅為10%，小于Re=0時的情況，這也證明了Re=5時渦的畸變較小。

圖15 順時針渦（V2）的最大渦量

軸向流動的存在可以穩(wěn)定泰勒-庫埃特流場，當存在軸流時，轉(zhuǎn)動雷諾數(shù)的臨界值比無軸流時高，因此，軸流是周向波的一個重要影響因素。圖16 為環(huán)隙中心柱面上的切向速度分布?？梢园l(fā)現(xiàn)，圖中周向波的一個明顯特征是當Re從2.5增加到5時，波高明顯下降，但波數(shù)不變。軸向流動對周向波動也有衰減作用，如圖17 所示為不同Re的波高。需要說明的是，本文的Re均不超過6，因為在較大的雷諾數(shù)時，流型可以由波狀渦流轉(zhuǎn)變?yōu)槁菪筒顪u流。如圖可見，隨著Re的增加，波高呈下降趨勢，這說明周向波會隨著軸向雷諾數(shù)的增加而減弱，也證實了軸向流動對流場的穩(wěn)定起著積極的作用。

圖16 不同雷諾數(shù)時環(huán)隙中心柱面上的切向速度

圖17 不同Re周向波波高

在Ta=139和Re=5時，設(shè)置了一個監(jiān)測點來監(jiān)測切向速度的瞬態(tài)行為。圖18(a)為切向速度的時域圖，如圖可見，速度的波動與無軸流時不同，此處存在兩種周期規(guī)律：第一種周期規(guī)律如圖中曲線所示，主要為渦的通過周期，即Tv=7.39s；第二種周期規(guī)律如圖中菱形圖標所示，這是由周向波動和軸向流共同作用導致的，周期為51.73s。由于渦沿著周向旋轉(zhuǎn)，同時渦被軸流推動通過監(jiān)測點。圖18(b)是通過快速傅里葉變換得到的切向速度波譜，非定常的切向速度隨著渦通過頻率0.135Hz及其高次諧波而振蕩。另一方面，一個較低主頻0.019Hz被捕獲，正好對應第2種周期規(guī)律。

圖18 切向速度波動的波譜

3 結(jié)論

本文對連續(xù)泰勒反應器中的波狀渦流進行了數(shù)值模擬研究。為了證明數(shù)值計算方法的可靠性，將數(shù)值模擬結(jié)果與文獻中的PIV 實驗結(jié)果進行了比較，結(jié)果具有良好的一致性。周向波破壞了軸對稱的泰勒渦結(jié)構(gòu)，導致非定常流動，包括渦及渦量隨周向位置發(fā)生周期性變化。切向速度沿著周向波動，最大的波動發(fā)生在環(huán)隙的中心。波高隨著泰勒數(shù)的增加而增大，直到Ta=180，隨著Ta繼續(xù)增加，波數(shù)增加。切向速度隨著時間波動，且在波動中不同軸向位置的波動特性存在明顯差異，軸向流的引入降低了渦及渦量的周期變化，減弱了波高，使流場更加穩(wěn)定。同時還發(fā)現(xiàn)軸向流動也影響著切向速度隨時間的變化，切向速度隨漩渦通過頻率及其高次諧波而振蕩。

符號說明

d——環(huán)隙寬度，m

f——頻率，Hz

H——波高，m

L——反應器長度，m

m——波數(shù)

R——徑向位置，m

Re——軸向雷諾數(shù)

Ri——內(nèi)圓柱半徑，m

Ro——外圓柱半徑，m

T——中心柱面旋轉(zhuǎn)周期

Ta——泰勒數(shù)

Tac——臨界泰勒數(shù)

Ti——監(jiān)測點的旋轉(zhuǎn)周期

Tv——渦的通過周期

t——時間，s

Va——軸向速度

Vr——徑向速度

Vt——切向速度

w——流體通過環(huán)形通道的平均軸向速度

Z——軸向位置，m

Г——縱橫比

ζ——軸向位置，Z/d

η——半徑比

θ——周向位置，rad

λ——軸向平均波長，m

υ——運動黏度，m2/s

ξ——渦變形程度

ρ——密度，kg/m3

Ωi——內(nèi)圓柱轉(zhuǎn)速，rad/s