層流圓管潛射流生成蘑菇形渦結構特性數值研究*

陳云祥 陳科 尤云祥 胡天群

(上海交通大學海洋工程國家重點實驗室，上海 200240)

(2012年12月17日收到;2013年2月18日收到修改稿)

1 引言

在海洋遙感衛星圖像上經常會觀測到海面上存在一類形似蘑菇狀的特殊流動結構[1].研究表明，這種特殊流動結構實際上是一類大尺度相干渦結構在海洋衛星遙感渭像中的影像，而且源于海洋中的某種動量源效應[2,3].潛艇推進器和冷卻水排放的力學效應相當于有動量傳遞給周圍流體，形成射流動量尾跡[4].在海洋背景流體作用下射流動量會形成一種類似蘑菇狀的大尺度相干渦結構，稱為蘑菇形結構(mushroom-like pattern).這類結構的水平尺度可達1—200 km，垂向尺度可達10—100 m，衰減所需時間可以達到1—30 d，在海洋衛星遙感渭像中的影像為一對旋轉方向相反的偶極子對[1].因此，針對射流動量在背景流體作用下產生蘑菇形渦結構的形成機理及其表現特征等問題進行研究，對發展潛艇非聲遙感探測技術有重要軍事價值，同時相關問題也一直是國際流體力學界的前沿熱點研究領域之一[3,5].

對于同種流體中以射流出口動量作為源動力的動量射流，研究表明，其流態特性與雷諾數Re有關.當Re＜500時，射流屬于層流狀態;當Re＞500時，在距離水平噴口的某個位置處，射流變得不穩定，湍流渦出現，而且隨著Re的增大，層流區逐漸消失，射流逐漸發展為完全湍流，直至其錐形角近似為一個常數[6-8].在湍流射流中包含著剪切層不穩定性觸發的豐富擬序結構，包括渦環、單螺旋結構和雙螺旋結構等擬序結構等[9-12].當湍流射流動量作用在密度分層或密度均勻淺水流體中時，由于密度分層產生的浮力效應或淺水條件下流體上下邊界產生的淺水效應的作用，動量尾跡的垂向運動會受到抑制，從而形成準二維偶極子形大尺度相干渦結構[13-17].

層流射流涉及到層流失穩等非線性流動現象[18-20].在密度均勻黏性流體中，Reynolds[21]和McNaughton等[22]在實驗中發現，在背景靜止流體作用下，層流射流液體在維持有限長度后會匯集在射流前部，形成足靴狀或蘑菇狀.Reynolds[21]在實驗中進一步發現在層流射流上可以看到軸對稱縮脹和正弦波動，Gill[23]的研究表明這種軸對稱失穩是有限振幅擾動造成的.Schneider[24]從理論上指出，層流射流在有限范圍內將被背景靜止流體逐漸滯止，并形成環形回流流動，而Petrov[25]和陳遠等[26]則采用實驗方法研究了層流射流中環形回流流動的形成機理.

研究表明，蘑菇形渦結構實際上是層流射流在背景靜止流體作用下產生的一類環形回流流動結構[21,22,25,26].Afanasyev等[27]采用Stokes近似理論研究了蘑菇形渦結構的漸近理論解，基于該理論解所獲得的蘑菇形渦的形態結構、運動軌跡及流場特征等與實驗結果相符.Voropayev等[28]進一步采用該理論方法建立了兩個層流射流生成蘑菇形渦結構的頭部迎撞問題的漸近理論解，結果表明基于該理論解所獲兩個蘑菇形渦結構的形態變化特征與實驗結果相符.

總體上，采用實驗和漸近理論分析等方法，目前對蘑菇形渦結構的形成機理及相關特征已經有了一定的認識.但在層流圓管潛射流的演化過程中，蘑菇形渦結構的射流長度、螺旋形渦環半徑及其包絡外形長度等幾何特征參數隨時間是變化的，而且其變化特性不僅與射流雷諾數有關，而且還與射流時間等有關.有關這些問題的研究不僅在學術上是重要的，而且對解讀其在衛星遙感影像中的表現特征及其空間分布規律等也是重要的，但迄今對這些問題的認識尚不清楚.

隨著計算機技術的發展，計算流體力學(computationalfluid dynamics，CFD)方法已成為研究射流演化和發展特性的重要手段之一.采用基于湍流模式的RANS方法以及大渦模擬方法等，在湍流射流形成機理及其大尺度擬序結構特性等方面，國內外學者已經開展大量研究工作[29-31]，但迄今還尚未見采用CFD方法對層流射流生成蘑菇形渦結構特性進行研究的文獻報道.有鑒于此，本文采用CFD方法，研究層流圓管潛射流在密度均勻黏性流體中的演化機理及其表現特征，分析蘑菇形渦結構的射流長度、螺旋形渦環半徑及其包絡外形長度等幾何特征參數隨時間的變化規律.

2 數值方法

本文考慮圓管射流在矩形容器中的演化問題.將一個圓形射流管置于盛水矩形容器中，其中心軸線位于一半水深處，一端位于矩形容器左側壁面處，如圖1所示.其中，射流管長為L0，直徑為D，矩形容器長為200D，寬和水深均為80D.建立直角坐標系，使坐標原點o位于圓管射流出口中心處，ox軸與圓管射流中心軸線重合且向右為正，oz軸垂直向上為正，oy軸垂直紙面向內為正.

圖1 計算區域側視圖及其坐標系

設射流流體與矩形容器中的背景流體為同種介質的流體，則圓管層流射流在矩形容器背景流體中的運動滿足如下不可壓Navier-Stokes方程：

其中，ν為運動黏性系數，p為動壓力，流體密度已被吸收入壓力中，u為流體運動的速度矢量.

為了研究層流射流在矩形容器背景流體中的演化問題，采用CFD方法對流場控制方程(1)進行數值模擬與分析.在數值模擬中，將射流管左端設置為入口邊界，入口速度為u0，在經過時間Tinj后停止射流;將矩形容器左側面也設置為入口邊界，入口速度為零;矩形容器底部和兩個側壁均設置為壁面邊界條件，自由面和右端出口均設置為自由出流條件;射流管外壁和內壁均設置為壁面邊界條件;以進口邊界條件作為流場計算的初始條件.此外，為保證數值模擬的穩定性，要求在射流管出口處的流動達到充分發展階段，因此其長度需要超過管內流動發展段的長度，將其設置為L0=10D[32].

計算網格均采用結構化六面體網格，將計算域劃分為三個子區域，分別為射流管內部區域、外部主體區域和外層區域.主體區域是一個環形結構的區域，其外徑為60D.采用等分圓周和等分半徑的方法獲得六面體網格，外層區域采用分塊方法將其劃分為六面體網格.射流管內部區域網格劃分足夠密集，而外層區域幾乎沒有射流流體到達，用Nx×Nr×Nθ表示環形主體區域的網格數，其中Nx為沿射流軸向的網格劃分數，Nr和Nθ分別為徑向和圓周方向網格劃分數.在本文中，采用的網格為300×120×48.在圖2中，給出了含射流管區域網格的橫視圖，圖2(a)中白點為射流管的內部區域，其網格劃分如圖2(b)所示.

圖2 含射流管區域網格劃分橫視圖

使用商業軟件FLUENT對所述問題進行數值模擬與分析，在求解中采用有限體積法離散動量方程和連續性方程，對流項離散采用三階精度QUICK(quadratic upstream interpolation for convective kinetics)格式，擴散項離散采用二階中心差分格式，壓力速度耦合迭代采用PISO(pressure implicit with splitting of operators)算法.計算初始時間步長為Δt=5×10-4s，計算過程中根據每個時間步長的收斂情況逐漸增加時間步長以縮短計算時間.

由于圓管射流的管徑遠小于矩形水箱寬度和水深，因此可將其等效為一個點動量源，其動量流量為[33]

3 結果與分析

3.1 蘑菇形渦形成機理

由圖3(a)和圖4(a)可知，當射流啟動后，在管內射流動量的作用下，射流液體從水平管出口處流出，并在背景流體中沿縱向以直線軌跡運動，形成軸對稱的柱狀主干流動，其形態結構如圖5(a)所示，這一階段稱為射流啟動階段.由圖3(b)和圖4(b)可知，當射流液體在主干流前方匯集的過程，由于背景流體的阻滯作用，主干流開始出現橫向輸運現象，并在其頭部逐漸形成一個類似半球形的圓突結構，其形態結構如圖5(b)所示，這一階段稱為射流過渡階段.

圖3 當Re=160和時，蘑菇形渦生成演化過程三維結構的數值結果(觀察視角：方位角143°，仰角110°) (a)t?=0.045;(b)t? =0.089;(c)t? =0.22;(d)t? =2.2;(e)t? =4.5;(f)t? =17.8

圖4 當Re=160和時，蘑菇形渦生成演化過程二維俯視結構的數值結果 (a)t?=0.045;(b)t?=0.089;(c)t? =0.22;(d)t? =2.2;(e)t? =4.5;(f)t? =17.8

圖5 射流演化形態結構示意圖

在過渡階段之后的階段稱為射流發展階段，這個階段包括三個過程.第一個過程如圖3(c)和圖4(c)所示，在射流液體被卷吸進入圓突結構的過程中，圓突結構與背景流體形成一個剪切層，并在剪切層作用下發生彎曲，其形態結構如圖5(c)所示.第二個過程如圖3(d)和圖4(d)所示，在圓突結構彎曲的過程中，會逐漸轉變為與射流方向相反的回流，其中回流液體緊貼在主干流周圍，被包裹在圓突結構內部，其形態結構如圖5(d)所示.第三個過程如圖3(e)和圖4(e)所示，隨著射流液體不斷被卷吸進入圓突結構內部，回流逐漸遠離主干流，使得圓突結構不斷“膨脹”，并在剪切層的作用下產生二次回流，其形態結構如圖5(e)所示.

在射流發展階段，射流形態在橫截面上呈現為螺旋渦結構，而其三維結構則表現為螺旋渦環形式.由于螺旋渦環的包絡外形猶如張開的蘑菇菌傘，而射流主干流形似菌柄，因此這是一種典形的蘑菇形渦結構.當蘑菇形渦結構形成后，在射流液體的不斷注入下隨時間穩定增長.當射流停止后，主干流動逐漸消失，但殘余射流液體繼續進入蘑菇形渦結構內部，直到射流液體全部被吸入其內部，最終形成“無柄”的蘑菇形渦結構，其形態結構如圖5(e)所示.

由圖3(f)和圖4(f)可知，當射流液體全部被吸入蘑菇形渦結構內部后，主干流動完全消失，渦環逐漸脫離蘑菇菌傘，其形態結構如圖5(f)所示.由于螺旋渦環得不到射流提供的動量，其發展受到抑制，蘑菇形渦結構開始逐漸衰退，最終在背景流體中侮全耗散，這個階段稱為射流衰退階段.

3.2 蘑菇形渦的幾何結構參數特性

結果表明，在圓管潛射流的三個不同階段，其演化特征有著較大的差異.為了對三個不同演化階段中蘑菇形渦幾何結構的變化規律進行定量化分析，定義如下三個特征參數，分別稱為射流長度L、螺旋形渦環半徑R和螺旋渦環包絡外形長度d，如圖6所示.其中，射流長度為射流頭部與射流噴嘴出口O之間的距離，螺旋渦環包絡外形長度為射流頭部與回流所在位置的縱向間距.O’為射流虛擬原點，T，G和M分別為蘑菇形渦結構的二次回流點、幾何中心和動量源作用中心，具體討論見下一節.

圖7給出了L?，R?和d?隨無量綱時間t?變化特性的數值結果，其中S?表示L?，R?或d?.由圖可知，隨著無量綱時間t?的變化，S?與t?具有不同的對應關系，整個曲線有兩個明顯的轉折點，這兩個轉折點正好對應于射流演化的三個階段，即啟動階段、發展階段與衰退階段.

圖6 蘑菇形渦幾何結構參數定義示意圖

圖7 當Re=160和時，L*、R*和d*隨無量綱時間t*的變化特性

在射流啟動階段，R?近似為一個常數，雖然L?和d?均隨無因量綱時間t?近似線性變化，但其斜率是不同的.另一方面，在射流過渡階段，d?近似為一個常數，雖然L?和R?均隨無因量綱時間t?近似線性變化，但其斜率也是不同的.這意味著射流啟動和過渡階段的R?，L?和d?都不是自相似或自保持的.但當流動進入發展階段后，如圖7所示，在對數坐標系下，R?，L?和 d?隨無量綱時間 t?1/2的變化規律近似為一條直線，這意味著在射流發展階段中蘑菇形渦結構是自相似或自保持的.利用圖7中的數據，采用數據擬合方法可得

從機理上看，層流射流的發展是依靠射流提供動量的不斷作用來維持的.當射流動量在作用一段時間消失后，射流的發展將必然無法繼續維持，射流結構不可避免地會進入衰退階段.由圖7可知，當射流演化進入衰退階段后，L?，R?和d?隨無量綱時間t?的變化規律發生顯著變化，而且并不存在一個明顯的過渡期.其中，L?和R?近似與t?1/5相關，而d?近似為一個常數，即此時螺旋渦環包絡外形長度已不再隨時間發生變化，這意味著在衰退階段中蘑菇形渦結構已不在具有自相似或自保持性.

下面進一步考察不同雷諾數和射流時間的組合對蘑菇形渦結構參數特性的影響，結果如圖8所示.由圖可知，對發展階段的蘑菇形渦結構，在雷諾數和無量綱射流時間的不同組合下，其L?，R?和d?隨t?的變化規律仍然滿足(5)式中的正比關系，即雷諾數和無量綱射流時間的不同組合并不影響L?，R?和d?隨t?的變化規律，其影響主要表現在啟動和衰退階段.這意味著在發展階段中雷諾數和無量綱射流時間的不同組合并不影響蘑菇形渦結構的自相似或自保持性.

下面對數值結果進行簡化理論分析與解釋.在圓管潛射流的發展階段，可以將其視為一個軸對稱局部點動量源在一個無界同質流體中的演化問題.由動量平衡及質量守恒條件可得[27]

其中，J=(J,0,0)為射流動量矢量，δ(x)為狄拉克函數，H(t)為單位躍階函數，即當t＜0時，H(t)=0;當t=0時，H(t)=0.5;當t＞0時，H(t)=1.

圖8 不同雷諾數和射流時間組合下，L?，R?和d?隨無量綱時間t?的變化特性

對處于發展階段的射流，可將其演化視作一個遠場問題，即在坐標原點處的奇性可以忽略.忽略掉動量方程中的非線性部分，在線性近似下，可得(6)式的時間依賴線性解為[27]

上式為歐拉觀點下的速度場解，若以拉格朗日的觀點來考量，從射流啟動瞬間，設有一個流體質點從射流的水平管出口流出，沿ox軸正方向沿著射流中心軸線運動，并且處于射流的最前端.那么，此流體質點隨時間所經過的位移即射流最前端的位置L，設其運動速度為uhp，此即射流頭部的運動速度.由(7)可得

將uhp=dL/dt代入(8)式可得

對蘑菇形渦結構，其體積V的增長是由射流流體的進入所引起的，參考液滴模形[6]，V隨時間的變化率可表示為

其中，uB為蘑菇形渦結構的整體水平移動速度，um為射流中心軸線上縱向速度的最大值.由文獻[7]可知，um和uB皆與t1/2成正比關系，而x也與t1/2成正比關系，故(10)式可改寫為

其中，CB為一個常數.

由(11)式可得

根據蘑菇形渦結構的空間結構特征，其體積可表示為V=V(R,d)，將啟動階段結束時蘑菇形渦結構的體積作為其初值，即

其中，R0和d0分別為啟動階段結束時蘑菇形渦結構的半徑及其包絡外形長度.

根據發展階段蘑菇形渦結構在其演化過程中的自相似性，結合(13)式可得，V∝R2和V∝d，再由(12)式可得

由此可見，在發展階段，蘑菇形渦結構的無量綱射流長度L?、螺旋形渦環半徑R?和螺旋渦環包絡外形長度d?均與t?1/2成正比關系，理論預測結果與數值結果一致.

下面考慮網格適應性問題.為此，將300×120×48稱為網格1，同時采用另外兩種網格，分別為，網格2：200×120×64和網格3：200×120×48.在三套網格下，L?，R?和d?的數值結果如圖9所示.由圖可見，L?和R?在這三套網格下的差異均較小，特別是在網格1和2下基本重合.d?在網格1和2基本重合，但在啟動階段網格3與網格1及2的結果有較大差異.由此可見，采用網格1進行計算是合理可靠的.

圖9 在三套不同網格下，L?，R?和d?的數值結果

3.3 速度與渦量場變化特性

在湍流射流的情況，保持射流出口速度不變的區域稱為射流核心區，沿著水平管縱向從出口至核心區末端的區域稱為射流的初始段，在初始段下游區域絕大部分為充分發展的湍流混摻區，稱為射流主體段.湍流射流的一個重要特點是：在射流主體段，各斷面縱向速度分布是自相似的或自保持的[34].

圖10 當 Re=160，和t?=3.6時，不同斷面上縱向速度分布數值結果

將在蘑菇形渦結構內部位于射流中軸線上速度的局部極大值點稱為動量源作用中心，設G為蘑菇形渦結構的幾何中心，而uM和uG分別為動量源作用中心和幾何中心沿中心軸線的縱向速度.記u?=uM/u0或u?=uG/u0，則在蘑菇形渦結構的發展階段u?隨時間變化特性的數值結果如圖12所示.由圖可知，蘑菇形渦結構的動量源作用中心M和幾何中心G沿中心軸線的縱向速度幾乎是相同的，這意味著可以將幾何中心近似視為蘑菇形渦結構的動量源作用中心.利用圖12中的數據，采用數據擬合方法可得

圖11 當Re=160和時，斷面中心軸線上縱向速度變化特性的數值結果

圖12 當Re=160和時，動量源作用中心和幾何中心沿中心軸線的縱向速度隨時間變化特性的數值結果

最后考察渦量-流函數特性，其中流函數ψ控制方程為?2ψ=-ωz′.設ud和vd分別為蘑菇形渦結構整體移動速度的縱向和橫向分量，將流函數ψ修正為 ψ′=ψ-udy′+vdx′.在圖 14中，給出了在三個不同時刻下Mx′y′平面上渦量-流函數散點圖.由圖可知，所得散點圖由兩個中心對稱的分支組成，由于在動量源作用中心后方區域渦量ωz′存在紐結現象，因此各分支的散點并不集中在某個曲線附近，而是形成一個形似“鐮刀”形的區域，這與準二維偶極子渦的渦量-流函數散點圖特性不同[36,37].將上分支右側邊緣與下分支左側邊緣連成一條曲線，由圖可知，隨著演化時間的增大該曲線逐漸近似成為一條直線.

圖13 當 Re=160時，在t?=4.5 時刻渦量 ωz′變化特性的數值結果

圖14 當Re=160和時，Mx′y′平面上渦量 -流函數散點圖

4 結論

采用基于不可壓Navier-Stokes方程的CFD方法，研究了在不同雷諾數和射流時間組合下蘑菇形渦結構的形成機理及其表現特性，主要結果如下：

蘑菇形渦結構是層流射流在背景流體阻滯作用下產生的一類環形回流現象，其形成與演化過程可分為三個不同的階段，分別為啟動階段、發展階段和衰退階段.在啟動階段，回流液體緊貼在主干流動周圍，L?與d?均隨t?線性變化，而R?則近似為一個常數.在發展階段，回流逐漸遠離主干流并產生二次回流，蘑菇形渦結構具有自相似性，L?，R?和d?均與t?1/2成正比關系，而且雷諾數和無量綱射流時間不影響該正比關系.在衰退階段，L?和R?與t?1/5相關，而d?則近似為一個常數.

在發展階段，射流軸線上的縱向速度存在局部極小和極大值點，局部極小值點與蘑菇形渦結構的二次回流點位置重合，局部極大值點與蘑菇形渦結構的幾何中心重合，該點亦為蘑菇形渦結構的動量源作用中心.在過射流中心軸線的平面上，垂向渦量關于該中心軸線反對稱，而且在動量源作用中心后方垂直于射流軸線的直線上垂向渦量分布出現紐結現象.渦量-流函數散點圖由兩個中心對稱的分支組成，但各分支的散點并不集中在某個曲線附近，而是隨著演化時間的增大上分支右側邊緣與下分支左側邊緣連線逐漸近似成為一條直線.

[1]Ivanov A Y，Ginsburg A I 2002 Proc.Indan Acad.Sci.(Earth Planet.Sci.)111 281

[2]Fedorov K N，Ginsburg A I 1989 Mesoscale/Synoptic Coherent structures in GeophysicalTurbulence 50 1

[3]Van Heijst G J F，Clercx H J H 2009 Annu.Rev.Fluid Mech.41 143

[4]Meunier P，Spedding G R 2006 J.Fluid Mech.552 229

[5]WeiG，DaiS Q 2006 Advances in Mechanics 36 111(in Chinese)[魏崗，戴世強2006力學進展36 111]

[6]Abramovich S，Solan A 1973 J.Fluid Mech.59 791

[7]Sozou G，Pickering W M 1977J.Fluid Mech.80 673

[8]Sozou G 1979 J.Fluid Mech.91 541

[9]Fan Q L，Wang X L，Zhang H Q，GuoY C，Lin W Y 2002 Advances in Mechanics 32 109(in Chinese)[范全林，王希麟，張會強，郭印誠，林文漪2002力學進展32 109]

[10]Yule A J 1978 J.Fluid Mech.89 413

[11]TsoJ，Hussain F 1989 J.Fluid Mech.203 425

[12]Fiedler H E 1988 Progress in Aerospace Science 25 231

[13]Sous D，Bonneton N，Sommeria J 2004 Phys.Fluids 16 2886

[14]Sous D，Bonneton N，Sommeria J 2005 European Journalof Mechanics B/Fluids 24 19

[15]Praud O，Fincham A M 2005 J.Fluid Mech.544 1

[16]Meunier P，Diamessis P J，Spedding G R 2006 Phys.Fluids 18 1

[17]Chen K，You Y X，Hu T Q，Zhu M H，Wang X Q 2011 Acta Phys.Sin.60 024702(in Chinese)[陳科，尤云祥，胡天群，朱敏慧，王小青2011物理學報60 024702]

[18]Danaila I，Dusek J，Anselmet F 1997 Phys.Fluids 9 3323

[19]O’NeilP，Soria J，Honnery D 2004 Exp.Fluids 36 473

[20]Kwon S J，SeoI W 2005 Exp.Fluids 38 801

[21]Reynolds A J 1962 J.Fluid Mech.14 552

[22]McNaughton K J，Sincair C G 1966 J.Fluid Mech.25 367

[23]GillA E 1965 J.Fluid Mech.14 557

[24]Schneider W 1985 J.Fluid Mech.154 91

[25]Petrov P A 1973 Fluid Dynamics 8 190

[26]Chen Y，LiD H，E X Q 1991 Acta Mechanica Sinica 23 721(in Chinese)[陳遠，李東輝，鄂學金1991力學學報23 721]

[27]Afanasyev Y D，KorabelV N 2004 J.Fluid Mech.16 3850

[28]Voropayev S I，Afanasyev Y D，KorabelV N，Filippov I A 2003 Phys.Fluids 15 3420

[29]Qu Y P，Chen S Y，Wang X P，Teng S G 2008 Journalof Engineering Thermophysics 29 957[曲延鵬，陳頌英，王小鵬，騰書格2008工程熱物理學報29 957]

[30]ZhaoP H，Ye T H，Zhu H M，Chen J，Chen Y L 2012 Journalof Engineering Thermophysics 33 529[趙平輝，葉桃紅，朱昊明，陳靖，陳義良2012工程熱物理學報33 529]

[31]Du C，Xu M Y，MiJ C 2010 Acta Phys.Sin.59 6331(in Chinese)[杜誠，徐敏義，米建春2010物理學報59 6331]

[32]Durst F，Ray S，UnsalB，BayoumiO A 2005 J.Fluids Eng.127 1154

[33]Lee J H W，Chu V H 2003 Turbulent Jets and Plumes：A Lagrangian Approach(Dordrecht：Kluwer Academic Publisher)p6

[34]Schlichting H 1979 Boundary Layer Theory(7th Edn.)(New York：McGraw-Hill)230

[35]AntonioR，AntonioL S 2002 Phys.Fluids 14 1821

[36]Beckers M，Clercx H J H，vanHeijst G J F 2002 Phys.Fluids 14 704

[37]MeleshkoV V，van Heijst G J F 1994 J.Fluid Mech.272 157