展向旋轉平板Couette湍流中的多態現象

夏振華

(浙江大學 航空航天學院工程力學系, 浙江 杭州 310027)

0 引 言

旋轉條件下的壁湍流問題是自然界和工業工程中常見的湍流現象之一。考慮地球自轉作用，那么地球上的壁湍流問題幾乎都可以視為旋轉壁湍流問題。而航空發動機的轉子、渦輪間的流動也是典型的旋轉壁湍流問題。為了設計性能更好的航空發動機等旋轉機械，我們急需了解旋轉條件下的壁湍流的統計規律，以及對應的流動結構隨旋轉強度的變化規律。另一方面，航空發動機等旋轉機械中關鍵湍流問題有很多，也非常復雜。為了更好地獨立出旋轉因素的影響，人們提出了一系列簡化的標準模型問題，其中展向旋轉平板Couette湍流(Spanwise rotating plane Couette flow, RPCF)就是其中之一[1-2]。

在旋轉系統下，由于科氏力的存在，基本的穩定性分析告訴我們，旋轉會根據其方向的選擇性使得剪切流動變得更穩定或者更不穩定。如果旋轉方向和平均流渦量同號(cyclonic rotation)，那么旋轉會使得流動更穩定；如果兩者異號(anticyclonic roation)，那么旋轉會讓流動變得更不穩定。在展向旋轉的槽道湍流中，由于其平均速度存在極值，平均流渦量會在槽道中變號，因此旋轉會在一部分槽道中使得流動變得不穩定；而在另一部分槽道中，其會使流動變得更穩定。而在RPCF中，其平均流渦量基本是相同符號的，因此旋轉在整個槽道保持相同的效應，即或全場讓流動更不穩定；或全場讓流動更穩定[3-5]。

圖1為RPCF流動的示意圖。相距為2h的兩塊平板以相對速度2Uw沿x方向運動，而整個系統沿展向(z方向)以常角速度Ωz旋轉。該系統有兩個控制參數：雷諾數Rew=Uwh/ν和旋轉數Ro=2Ωzh/Uw。

圖1 RPCF簡圖[5]Fig.1 Sketch of the flow geometry in RPCF[5]

在過去的幾十年中，很多研究組針對RPCF開展了系列的實驗和數值模擬研究。在實驗方面，Alfredsson教授和其團隊做出了一系列工作[6-9]。特別需要指出的是，Tsukahara、Tillmark和Alfredsson[7]在2010年時總結并擴充了之前的實驗工作，系統地研究了參數范圍在0

在數值研究上，Andersson教授及其合作者首先針對Rew=1300時的RPCF開展了直接數值模擬(Direct numerical simulation, DNS)研究[13-17]。他們研究了弱旋轉時的流動統計以及流動結構，中高旋轉數時的旋轉效應以及極強旋轉時的流動統計。他們發現在一定旋轉數范圍里，湍流和二次流隨著Ro增加而增強，但是超過一定的旋轉數范圍，其湍流和二次流則會隨著Ro增加而減弱。在極高旋轉數條件下，二次流消失，這個時候的湍流特征與無旋轉時的流動特征差別很大。Salewski和Eckhardt[18]用DNS研究了四個不同雷諾數(Rew= 650、1300、2600和5200)下，以壁面剪切為表征的動量流隨Ro的變化行為。他們考察了-0.04≤Ro≤0.6之間的算例，發現變化曲線在兩個較低雷諾數情況下只有一個極值點；而在兩個較高雷諾數下，變化曲線有兩個極值點，其中一個極值點位于Ro=0.02附近(較小旋轉數)。Gai等[5]也針對Rew=1300時的RPCF流動開展了一系列DNS模擬研究，他們分析了中心線處的平均剪切率、壁面摩擦雷諾數、體平均的總湍動能、體平均的二次流湍動能、體平均的剩余場湍動能等隨著旋轉數的變化規律。他們還分析了不同旋轉數下的流動結構的異同。

近些年，本文作者與其合作者還發現了RPCF中存在多湍態現象[19-22]。在經典的湍流理論中，人們都假設雷諾數足夠高的時候，湍流是各態遍歷的，也就是說所有的湍流態之間都是等價的，或者說湍流態是唯一的。在這個假設下，雖然瞬時場跟初始條件有關，但是湍流的統計規律是與其初始流場無關的[23]。湍流的唯一統計態假設是經典湍流理論的基石，也是很多湍流結論和假設能夠從低雷諾數外推到更高雷諾數的原因之一[24]。然而，近些年來，一些實驗和數值模擬的結果顯示，部分湍流中可能存在不唯一的湍流態[25]。在多態條件下，相同的控制參數，可能會出現不同的湍流統計和流動結構。例如，在Taylor-Couette流動中，Huismann等[24]就通過實驗發現在Re為106量級下，通過不同途徑達到相同的內外筒轉速比a時，流動的Taylor渦對數和扭矩等結果在一定的參數范圍里(0.17

本文將就RPCF中的多態現象的研究現狀進行綜述介紹。

1 RPCF的控制方程和數值解法

1.1 控制方程

RPCF的幾何構型如圖1所示。在通常的設置中，我們假設流向為x，法向為y，展向為z(或者x1、x2、x3)，其對應的速度為u、v、w(或者u1、u2、u3)。用Uw，h無量綱化之后，系統的控制方程為：

(1)

(2)

這里p是包含了離心力的有效壓力(effective pressure)。其中雷諾數Rew和旋轉數Ro的定義前文已給出。

1.2 數值解法介紹

在實際計算中，為了保證程序的精度，我們采用高精度Fourier-Chebyshev偽譜法對空間項進行離散。具體而言，在流向和展向周期方向，我們采用Fourier級數展開；而在法向采用Chebyshev多項式展開。在時間推進上，我們采用了半隱式離散格式，即線性項采用2階精度的Crank-Nicolson根式而非線性項采用的是2階的Adams-Bashforth格式。此外，我們用3/2法則消除混淆誤差。值得說明的是，我們的方法比Bech和Andersson[14-15]和Barri和Andersson[17]的數值解法精度更高。Bech和Andersson用的是2階精度的有限差分程序ECCLES[14-15]，而Barri和Andersson則用的是2階精度的有限體積程序MGLET[17]。

在實際計算中，我們在大小為Lx×Ly×Lz的計算域上開展計算，離散的網格為Nx×Ny×Nz。具體的計算網格，在后續會做介紹。

1.3 雷諾分解和三重分解

(3)

由于在RPCF中存在二次流，因此我們在分析的時候，還可以開展三重分解[5,14,19]：

(4)

(5)

(6)

即湍動能和雷諾應力都可以分解為由二次流貢獻的部分和剩余場貢獻的部分。

2 結果與討論

本部分將著重介紹相關的結果與討論部分。

2.1 旋轉數對湍流統計和流動結構的影響

Gai等人[5]針對RPCF開展了一系列的DNS模擬，基本控制參數如表1所示。其主要研究目標是旋轉數對湍流統計和流動結構的影響。在分析該問題時，研究者們一般重點關注兩個量：平均速度在中心線處的法向梯度Ψ和壁面摩擦雷諾數Reτ。

表1 計算參數Table 1 Computational parameters

圖2(a、b)分別給出了計算中20個算例的中心線速度的法向梯度的平均值及其標準差以及壁面摩擦雷諾數Reτ隨旋轉數Ro的變化規律。可以看到，在Ro從0增加到0.02附近時，Ψ逐漸減小；而當Ro繼續從0.02增加到0.9時，Ψ則隨著Ro增大而增加。非常有意思的是，在Ro=0.02時，Ψ是負值，即在該旋轉數條件下，平均流動存在局部回流現象。Salewski和Eckhardt[18]在Ro=0.04時也發現了平均速度的局部回流現象。不過需要說明的是，他們所采用的流向計算域過小，對流動結構和統計有一定的影響(在下文，我們將討論這種差別)。而Kawata和Alfredsson[10]他們也通過實驗測到過這種平均速度局部回流現象。相反地，壁面摩擦雷諾數Reτ隨旋轉數Ro的變化則呈現先增加后減小的趨勢。其最大值約在Ro=0.15附近。在該問題中，Rew保持固定，因此Reτ可以表征壁面摩擦的大小。由圖可以看出，在00.611時，旋轉可以起到減阻的作用。

(b)

圖3給出了體平均的湍動能[k]y以及二次流和剩余場貢獻的部分[ks]y和[k″]y隨旋轉數Ro的變化曲線。其中，算子“[]y”定義如下

(7)

圖3 體平均的湍動能及其分量隨旋轉數Ro的變化規律。相關量都用Ro=0時的無量綱化[5]Fig.3 Volume-averaged kinetic energy and its two shears (the contributions from the secondary flows and the residual field) as a function of Ro. All quantities are normalized by at Ro=0[5]

由于湍動能等項已經是時間、流向和展向平均之后的量，因此其法向平均可以認為是體平均的結果。由圖中可以看出，總湍動能隨著Ro是先增加后減小的，在Ro=0.25附近達到最大值。即當Ro<0.25時，旋轉會增強湍流；而當Ro>0.25時，旋轉相反會抑制湍流。二次流部分的湍動能的結果同總湍動能基本保持類似的趨勢，在Ro>0.6時，二次流基本消失。剩余場部分的湍動能則體現出同總動能相反的趨勢。它是先減小，后增大，最后又減小的過程。在Ro=0.1附近其達到極小值，而在Ro=0.5時，其已經非常接近總湍動能的值。

圖4是六個不同旋轉數下二次流的統計圖。該圖中，橫截面速度向量由(vs,ws)表征；流向速度分量us則由等值面表征。從圖中可以看出，在Ro<0.25時隨著Ro增加，二次流逐漸增強；而當Ro>0.25時，二次流則隨著Ro增加而減弱，在Ro=0.7時已經幾乎無法看到明顯的二次流了。這和前面的圖3的結論基本一致。值得注意的是，當Ro=0.25和0.32時，二次流的大渦兩側存在次生的二次流。其中Ro=0.25時，次生二次流強于Ro=0.32時的結果。

(a) Ro=0

(b) Ro=0.01

(c) Ro=0.02

(d) Ro=0.25

(e) Ro=0.32

(f) Ro=0.7

2.2 弱旋轉條件下的流動結構和統計

在上一節中，我們提到Rew=1300、Ro=0.02時平均流在槽道中心位置有局部回流現象。本部分將介紹在弱旋轉Ro=0.02時的流動結構和統計分析。在Huang等人[20]的文章中，他們開展了一系列DNS研究，分析了網格分辨率、流向計算域、展向計算域和計算時間的影響。在計算結果中，他們發現，在Ro=0.02時，流動結構可能會隨著計算時長而改變，不過這種計算結果還跟流向計算域有關。當流向計算域較小，例如Lx=4π時，展向Lz=4π的槽道中二次流的渦對(Roll-cells)數恒定為3對；而當Lx≥8π時，展向Lz=4π的槽道中會先出現3對渦對而后流動會演變為2對渦對。在流動處在3對渦對時，平均速度在中心區存在局部回流；而當流動處于2對渦對時，該局部回流消失。如果Lx=8π，Lz=6π時，該流態轉換依舊存在，流動會從4對渦對轉換成3對渦對。圖5給出了該計算域下x=0位置處的中心線瞬時流向速度在展向和時間上的變化規律。可以看出，在長達5000個無量綱時間中，流動從層流加擾動初始解首先演化出4對流向渦對，該4對渦對的狀態持續了大約1800無量綱時間(200～2000)。之后該流動的渦對開始出現合并，最后在3000個無量綱時間的時候，流動變成3對渦對結構。值得說明的是，我們嘗試過計算超過10 000個無量綱時間，流動后期依舊保持在3對渦對的狀態。

圖5 在x=0處中心線流向瞬時速度信號u(0,0,z,t)隨時間的變化關系[20]Fig.5 Contours of u(0,0,z,t) as time evolves[20]

圖6給出了4對渦和3對渦結構下，流場的湍流統計規律。可以看出，兩者之間是存在明顯差別的。特別的，在4對渦時，平均速度在中心線處存在局部回流(從圖6(d)中雷諾剪應力大于1可以推出來)；而當流動變成3對渦的時候，該平均速度的局部回流消失，雷諾剪應力的最大值小于1。

(a) 平均速度

和

(d) -〈u′v′〉+

需要說明的是，在前人的工作中，Lz=4π的槽道寬度中只發現了3對渦的狀態[5,14-15]。我們認為這主要原因在于過去的計算時間不夠長。在過去的計算中，最多只算了幾百無量綱時間[5]，圖5顯示了計算時間如果不夠長，容易忽略掉后面穩定的流態。因此，我們建議算湍流問題時，要用更長的計算時間。

2.3 多態現象及其統計結果、魯棒性分析

事實上，RPCF中的多湍流態現象在文獻[19]中已有論述。在該文章中，本文作者和合作者針對Rew=1300，Ro=0.2的RPCF流動開展了DNS模擬。由于在文獻[5]的基礎上，我們知道隨著Ro從0.02增加到0.2時，二次流是會變得更強更穩定的；而在前期研究Ro=0.02時，我們也發現了Ro=0.02時存在流態轉換過程。因此，在設計計算的時候，我們選取了兩個不同的初始場，一個初始場是Ro=0.02時較前期的瞬時場，此時該流場擁有3對渦對結構(Lz=4π)；另一個初始場則選取Ro=0.02時后期的瞬時場，該流場擁有2對渦對結構。在計算中，兩個初始場分別按照相同的程序，在相同的控制參數，相同的計算域和相同的網格分辨率條件下進行時間推進。經過超過3000個無量綱時間的計算，我們發現流動的最終統計結果和流動結構不一樣。另外，我們還在流向和展向進行了網格加密，結果依舊出現多態現象。

圖7給出了穩態條件下，兩種不同流態下的二次流的橫流速度矢量圖和流向分量的等值線圖。可以看出，兩組不同的初始條件會等到不同的流動結構：一個會得到2對渦對結構；另一個則會保持3對渦對的狀態。

(a) 2對渦對的結果

(b) 3對渦對的結果

兩種不同流態下，其平均速度存在差異。為了說明問題，表2中給出了兩組不同網格下，2對渦和3對渦的一些統計結果。從表中的數據可以看出：相同流態下，兩組不同的網格下壁面摩擦雷諾數Reτ相差很小，誤差都小于1‰。而不同流態之間的Reτ之間的差別則明顯增大，差別約為6%。類似地，平均速度在中心線處的法向梯度也呈現類似的行為，只是誤差相對較大。在3對渦流態下，兩組不同網格的計算給出的值約為0.2；而在2對渦流態下對應的值減小到約0.14。從穩定性分析的角度看，前者在中心區趨向于中性穩定；而后者則處于穩定區域，這也可以解釋為何2對渦條件下湍動能在中心區會相對較小。

表2 不同流態下的湍流特征統計[19]Table 2 Turbulent statistics at different states[19]

圖8給出了2對渦和3對渦時的湍流二階統計結果。可以看出，兩種不同流態之間的統計存在非常明顯的差異，而同一種流態在不同網格下的統計差異幾乎可以忽略。在靠近壁面的區域，3對渦條件下的流向和展向速度脈動更強；同時，其對應的法向速度脈動和雷諾剪應力則在槽道中心區更強。

為了驗證多態現象的存在是物理的，而不是計算域造成的虛假結果，我們在文獻[19]中擴大了流向和展向的計算域，計算結果基本還能保持多態的結果。為了進一步驗證多態現象的魯棒性，我們還考察了多態現象對初始流場的依賴性。我們構造了一個新的流場ui：

(8)

更進一步，我們還模擬實驗的方式[25]，在計算中逐漸增加旋轉數或者逐漸減小旋轉數，進而來研究多態現象的存在區間。在文獻[21]中，Huang等人開展了兩組DNS模擬。在相同的計算域和計算網格下，兩組計算的初始場都是層流解加擾動的無散場。其中一組的旋轉數Ro逐漸從0.01增加到0.02，0.03，0.04，0.05，0.1，0.2，0.3，0.5，0.6；而另一組則按照相同的Ro從0.6逐漸減小到0.01。在每個Ro上，都計算1250個無量綱時間，后一個算例(Ro)是以前一個算例的最后一個時刻的瞬時流場作為初始場。

(a) 〈u′u′〉

(b) 〈v′v′〉

(c) 〈w′w′〉

(d) -〈u′v′〉

圖9 不同參數α下初場的最終流態示意圖[22]Fig.9 A sketch of final states with different initial conditions defined by α[22]

圖10給出了法向瞬時速度的前四個模態的幅值隨計算時間的變化曲線。其中，模態幅值定義如下：

(9)

其中，Fz(f(z))表征的是函數f(z)的離散傅里葉變換，“〈〉x”表征流向平均。通常情況下，幅值最大的kz值可以視為渦對數。由圖10可以看出，在0.03≤Ro≤0.3范圍里，隨著旋轉數Ro逐漸增加，kz=2的模態逐漸占主導；而當Ro逐漸減小時，kz=3的模態逐漸占主導。

(a) 隨著Ro增加的結果

(b) 隨著Ro減小的結果

圖11給出了在Ro在0.01到0.5之間的流向瞬時速度在x=0、y=0處的時間歷程。其中，圖11(a)給出的是隨著Ro增加的情況，圖11(b)給出的是隨著Ro減小時的情況。可以看出在0.03≤Ro≤0.3范圍里，隨著Ro增加時，流場呈現2對渦；而隨著Ro減小時，流場則有3對渦。圖11的結果與Huisman等人[25]在Taylor-Couette流動的實驗中的結果非常相似。

(a) Ro逐漸增加

(b) Ro逐漸減小

我們還分析了壁面摩擦雷諾數Reτ、體平均的湍動能以及體平均的二次流湍動能隨著Ro增加和減小的結果。圖12中給出了體平均的湍動能隨著Ro增加和減小時的變化規律。由圖可以看出，在0.03≤Ro≤0.3范圍里，兩組計算出來的體平均湍動能在相同的Ro下存在明顯差別，3對渦時(Ro減小方向)的湍動能明顯大于2對渦時的湍動能。在文章中，我們還發現Reτ和體平均的二次流的湍動能都存在類似于體平均湍動能一樣的遲滯環[21]。

需要說明的是，在這兩組計算中，10個Ro的具體數值并非預先有意設計的，而是在計算中隨機選取的。計算時間間隔1250個無量綱時間也不是特意選取的。在所有計算結果出來之前，我們也不能保證我們能觀測到遲滯環現象。最終結果我們發現了遲滯環，這至少說明了多態現象可以在很大的Ro范圍里存在。我們承認如果改變Ro的具體值選取或者改變時間間隔，結果可能不一樣。事實上，在Huisman等人的實驗研究中，他們也是在一些情況下遇到了多態現象[25]。

圖12 體平均的湍動能隨著Ro增加和減小的變化規律[21]Fig.12 Changes of the volume-averaged turbulent kinetic energy as Ro increases or decreases[21]

2.4 多態條件下的能量流動分析

在Taylor-Couette湍流中多態現象的工作中，Huisman等人曾經猜測大尺度相干結構的選擇性作用在多態現象中占據重要地位[24]。在RPCF中，二次流就是這種大尺度的流動結構。為了驗證Huisman等人的假設，我們研究了二次流和剩余場之間的能量傳輸過程。參考Bech和Andersson[14]和Cai等人[26]的工作，我們將速度場分成四部分：

(u,v,w)=(〈u〉,〈v〉,〈w〉)+(us,0,0)+

(0,vs,ws)+(u″,v″,w″)

(10)

(11)

(12)

(13)

(a) 3對渦的流態

(b) 2對渦的流態

3 結論與展望

我們通過一系列直接數值模擬在Rew=1300的RPCF中發現了多湍流態現象。在不同的湍流態下其流動結構和湍流統計差別明顯。該多態現象比較魯棒，可以在非常大的Ro范圍里被觀察到。

我們還分析了不同湍流態下，能量在平均場、二次流場和剩余場之間的傳遞關系。通過數據分析我們發現，3對渦時(Lz=4π)，二次流較強，此時二次流相當于能量傳輸的一個橋梁：平均場將大部分能量傳輸給二次流，再由二次流傳給剩余場；而當流動處在2對渦的流態時，此時二次流相對較弱，平均流雖然傳輸了很多能量給二次流，但是大部分都用于維持二次流本身，只有少部分能量能傳輸給剩余場。而剩余場的能量主要是通過平均場直接傳輸過來的。這種能量傳輸的不同表現，側面驗證了大尺度結構在多態現象中占據重要地位。

經典的湍流理論假設湍流是各態遍歷的，是跟初始場無關的。這是湍流建模和模擬的基礎之一。多態現象的出現將對湍流建模和模擬提出新的挑戰。如果在模型的構建(生成項的大小都不一樣)和驗證中，將多態現象考慮進去是值得研究的課題。此外，多湍態現象在工程問題中是否存在？如果存在，它會對工程設計和安全有什么影響？這些也是非常值得研究的問題。