展向周期變化槽道流中全局穩(wěn)定性的數(shù)值模擬

孔 瑋，羅紀(jì)生

(天津大學(xué)機(jī)械工程學(xué)院，天津 300072)

展向周期變化槽道流中全局穩(wěn)定性的數(shù)值模擬

孔 瑋，羅紀(jì)生

(天津大學(xué)機(jī)械工程學(xué)院，天津 300072)

用數(shù)值模擬的方法研究了基本流在展向周期變化的槽道流的全局穩(wěn)定性問題．首先用線性穩(wěn)定性理論求得的擾動特征函數(shù)作為初始值，計(jì)算擾動隨時(shí)間的演化．經(jīng)過充分的發(fā)展后，得到了全局穩(wěn)定性的特征解．給出了不同展向變化幅值下的計(jì)算結(jié)果和全局穩(wěn)定性下的中性曲線，分析了全局中性時(shí)雷諾數(shù)與展向不均勻程度的關(guān)系以及特征函數(shù)的性質(zhì)．

全局穩(wěn)定性；展向周期變化；槽道流；數(shù)值模擬

在大多數(shù)情況下，流動穩(wěn)定性理論中研究討論的是基本流只與一個(gè)坐標(biāo)有關(guān)的流動，如槽道流中，基本流只是壁面法向坐標(biāo)的函數(shù)．但一般來講，基本流可能還與其他的坐標(biāo)有關(guān)，這類流場的穩(wěn)定性分析就比較復(fù)雜了．因此，在流動穩(wěn)定性分析中，擾動的穩(wěn)定性問題通常用局部分析(local analysis)和全局分析(global analysis)兩種方法來研究．

局部的分析方法研究的是基本流只與一個(gè)坐標(biāo)有關(guān)或基于平行流假設(shè)下近似的只與一個(gè)坐標(biāo)有關(guān)的情況，即針對的是平行剪切流，通過求解 Orr-Sommerfeld方程(簡稱O-S方程)的特征值問題來進(jìn)行分析；而對于基本流是多個(gè)變量的函數(shù)，如果采用局部分析的方法，需要假定流場隨法向坐標(biāo)以外的其他變量是緩變的，從而轉(zhuǎn)化成參數(shù)形式的 O-S方程進(jìn)行求解，這樣在不同的參數(shù)下特征值都會不同，往往無法對整體穩(wěn)定性進(jìn)行準(zhǔn)確預(yù)測，所以通常需要采用全局的分析方法，按照全局理論，流場應(yīng)該在整體上具有一致的增長衰減性，即得到的是全局特征解．

全局穩(wěn)定性理論研究起始于 Huerre和Monkewitz[1]在 1990年的工作，2005年 Chomaz[2]對此問題進(jìn)行了回顧，在此期間也有許多研究工作有過一定程度上的進(jìn)展[3-8]．雖然流動穩(wěn)定性的全局分析與局部分析從概念上講是類似的，但局部分析研究的是常微分方程的特征值問題，而全局分析研究的是偏微分方程的特征值問題．因此，全局分析的計(jì)算量是巨大的，通過 QZ方法[9]直接求解這么大的特征值問題，計(jì)算效率非常低，有時(shí)甚至是不可能的．雖然可以用一些其他的方法進(jìn)行求解，但對計(jì)算技巧及計(jì)算量的要求都比較高．

由于通常研究的槽道流只和壁面法向坐標(biāo)有關(guān)，是一種典型的平行流，因此對于這種流動的穩(wěn)定性可以由局部線性穩(wěn)定性理論[10]進(jìn)行分析，并且早已有過研究且已經(jīng)成熟．但在工程或真實(shí)情況下流動是不可能完全平行的或只隨一個(gè)坐標(biāo)變化，比如Stokes層就是一種隨時(shí)間變化的流動，而研究此類流動更顯得具有實(shí)際意義．對于本文研究的展向周期變化的槽道流，則需要求解全局特征值問題，但是其計(jì)算量非常大．本文則通過數(shù)值模擬的方法研究了這種流動的全局穩(wěn)定性問題．

1 基本流、控制方程與數(shù)值方法

1.1 基本流

本文研究的是基本流展向周期變化的槽道流，即流動除了流向的壓力梯度外，還由展向周期變化的體積力引起．

槽道流的坐標(biāo)如圖 1所示．其中 x、y、z分別為流向、法向和展向，h為半槽道寬，Um為最大基本流速度．

用 Um作為特征速度，半槽寬h作為特征長度，密度ρ作為特征速度，h /Um作為特征時(shí)間尺度，作為特征壓力尺度，對流體力學(xué)的不可壓縮原始Navier-Stokes方程進(jìn)行無量綱化．無量綱化后的Navier-Stokes方程對應(yīng)的層流方程為

式中：ui是速度；是坐標(biāo)(x,y,z)；p是壓力；Re = Umh/ν是雷諾數(shù)；fxi是體積力項(xiàng)．

設(shè) U1、U2、U3分別為流向、法向和展向的基本流速度．考慮在流向加入沿展向周期性變化的體積力及流向的壓力梯度，可將式(1)化簡為

式中：b是基本流展向變化部分的的幅值；β0為展向基本波數(shù)．

當(dāng)b＝0時(shí)，基本流為通常槽道流的速度剖面；當(dāng)b≠0時(shí)，為展向周期變化的槽道流．圖 2給出了b＝0.2時(shí)基本流在 yz平面上的等值線圖，可以看出，基本流速度 U1在展向是周期性變化的，圖中給出了沿展向2個(gè)周期的結(jié)果．

圖2 b＝0.2時(shí)基本流速度 U1在yz平面的等值線圖Fig.2 Contours of base flow U1in the yz plane when b＝0.2

1.2 控制方程

將基本流寫成只含有 y的部分和含有展向變化的部分，即

把與 Uz相關(guān)的項(xiàng)寫到方程右邊，得到數(shù)值求解的擾動方程(為求簡潔，略去擾動量的上標(biāo))為

其中

1.3 數(shù)值方法

本文在流向和展向采取傅里葉譜方法進(jìn)行計(jì)算，其中展向網(wǎng)格數(shù)依據(jù)計(jì)算具體情況而定；在法向使用Chebyshev配置點(diǎn)法，網(wǎng)格數(shù)為 129，時(shí)間上采取二階精度的顯隱差分格式[11]進(jìn)行計(jì)算，時(shí)間步長依據(jù)計(jì)算所取不同展向網(wǎng)格數(shù)而定．

將式(4)中的各函數(shù)在流向和展向展開為Fourier級數(shù)的形式，再在譜空間求解．

式中 ωr是由 O-S方程求解出的特征值實(shí)部；α 為流向波數(shù)．流向和展向計(jì)算域長度分別是

式(4)變?yōu)?/p>

其中

將方程(6)寫成算子的形式，即

式中：L(u)是與 Uz無關(guān)的項(xiàng)，采用時(shí)間隱格式求解；為與 Uz相關(guān)的項(xiàng)，采用顯格式進(jìn)行求解．于是在時(shí)間方向離散后的算子形式可寫成

時(shí)間方向離散后，式(6)可寫成

其中

聯(lián)立方程組(7)可解得η、vn、pn，求出 vn后，再代入式(7)中第 4式求出ξ，再由出了0～8階譜的結(jié)果)．從圖中可以看出，擾動在經(jīng)過前期的調(diào)整后，隨時(shí)間是呈現(xiàn)指數(shù)增長的，并且從求得

2 計(jì)算結(jié)果及分析

圖4 擾動速度 u各階展向譜沿法向最大幅值隨時(shí)間的演化Fig.4 Variation of maximum amplitude of disturbance velocity u with time

圖3 一般槽道流的中性曲線Fig.3 Neutral curve of common channel flow

2.1 展向不均勻性對穩(wěn)定性的影響

為了了解展向幅值 b對穩(wěn)定性的影響，下面以α＝1.02、Re＝8,000的情況為例來進(jìn)行說明．

圖4是b＝0.2時(shí)擾動速度u各階展向譜沿法向最大幅值隨時(shí)間的演化及其在計(jì)算前期的演化及擾動速度u的幅值取自然對數(shù)后隨時(shí)間的演化(其中給取對數(shù)后各階譜擾動速度幅值隨時(shí)間的變化曲線可以看出，各階譜互相平行，這證明它們的增長率是相同的，即各階譜擾動的特征函數(shù)都具有一致的增長率．同時(shí)為了驗(yàn)證特征函數(shù)的形狀是否不發(fā)生改變，用 0階譜沿法向的最大幅值為標(biāo)準(zhǔn)對各階譜進(jìn)行歸一化，得到不同時(shí)刻擾動速度u各階展向譜的特征函數(shù)的比較，見圖5(其中給出了0～5階譜的結(jié)果)．從圖5中可以看到，歸一化后的各階譜擾動特征函數(shù)的形狀不隨時(shí)間發(fā)生改變，與t＝500的結(jié)果重合，說明得到了定常解．其他兩個(gè)速度分量 v和 w的情況類似，在此不再給出．由此表明，擾動在整體上具有一致的增長，得到了擾動全局增長的特征解．

圖5 不同時(shí)刻以0階譜歸一化后的擾動速度u各階展向譜特征函數(shù)Fig.5 Eigenfuctions of disturbance velocity u normalized with Fourier mode n＝0 at different time

圖6 展向譜n＝0時(shí)擾動速度u及其增長率隨時(shí)間的演化Fig.6 Variation of maximum amplitude and growth rate of disturbance velocity u(n＝0) with time

下面討論展向幅值 b對全局穩(wěn)定性的影響．圖6(a)給出的是 b分別為 0.2、0.324,95、0.4時(shí)擾動速度 u的 0階展向譜幅值隨時(shí)間的演化．從圖中可以看出，隨著b的增大，擾動是逐漸趨向于衰減的，即b對擾動起著穩(wěn)定的作用．圖 6(b)給出了擾動的增長率隨時(shí)間的變化，即先將擾動幅值對時(shí)間求導(dǎo)數(shù)，然后再除以其本身，得到了增長率的大?。畯膱D中可以看到，b＝0.2時(shí)擾動的增長率大概是 0.001,3左右，而當(dāng) b＝0.4時(shí)增長率為-0.000,92，變?yōu)樨?fù)值，而當(dāng)b＝0.324,95時(shí)增長率大約為 0，即對應(yīng)全局的中性解. 由此可知，在同一雷諾數(shù)下，隨著b的增大，擾動波的增長率是逐漸減小直至變?yōu)樗p，一般總可以找到一個(gè)使得增長率基本為0即全局中性下的b值，統(tǒng)一稱其為b0．

通過計(jì)算，在不同的雷諾數(shù)下找到了相應(yīng)的 b0，圖7給出的是在不同α下Re與b0的關(guān)系曲線．從圖中可以看出，當(dāng)雷諾數(shù)處于一般槽道流中性曲線下分支附近時(shí)，b0隨著雷諾數(shù)的增大而增大；而當(dāng)雷諾數(shù)處于中性曲線上分支附近時(shí)，b0是隨著雷諾數(shù)的增大而減小的，并且隨著α的增大，曲線的范圍是逐漸縮小的．

圖7 b0隨Re的變化曲線Fig.7 Variation of b0with Re

2.2 全局中性情況下計(jì)算結(jié)果分析

以下分別給出了擾動既不增長也不衰減即中性的情況下，一些物理量的計(jì)算結(jié)果及其分析．

2.2.1 擾動速度u各階譜沿法向的分布

圖8 歸一化后擾動速度的特征分布Fig.8 Eigenfunctions of disturbance velocity u normalized with Fourier mode n＝0

首先，圖 8給出了α＝1.02、b＝0、Re＝5,772.3時(shí)擾動速度的特征分布，這是通常槽道流中性曲線臨界雷諾數(shù)對應(yīng)的速度特征函數(shù)的結(jié)果．

圖9 擾動速度 u各階展向譜的歸一化特征分布(Re＝5,800，8,000)Fig.9 Eigenfunctions of disturbance velocity u normalized with Fourier mode n＝0(Re＝5,800，8,000)

圖9和圖10分別給出的是α＝1.02時(shí)中性曲線下支界和上支界附近的不同雷諾數(shù)在對應(yīng) b0下擾動速度u展向譜歸一化后特征分布的比較結(jié)果，其中下支界附近給出的雷諾數(shù)分別為 5,800和 8,000，上支界附近給出的雷諾數(shù)分別為 22,000和 26,000，給出了展向譜 n＝1～4的結(jié)果．從比較結(jié)果上看，在中性情況下不同雷諾數(shù)下各階諧波沿法向的分布基本是類似的．其他兩個(gè)擾動速度分量 v和 w也具有相同的性質(zhì)．

圖10 擾動速度u的各階展向譜的歸一化特征分布(Re＝22,000，26,000)Fig.10 Eigenfunctions of disturbance velocity u normalized with Fourier mode n=0(Re＝22,000，26,000)

2.2.2 擾動速度u在xz平面上的等值線圖

為了比較不同雷諾數(shù)下擾動在 xz平面的分布，將譜空間的擾動變換到物理空間中，并在擾動沿法向最大值處給出 xz平面上 u的等值線圖．圖 11和圖12分別給出的是α＝1.02時(shí)，中性曲線下分支和上分支附近雷諾數(shù)的計(jì)算結(jié)果．由圖中可以看出，擾動在展向呈現(xiàn)明顯的周期性變化，在中性曲線下分支附近，當(dāng)雷諾數(shù)增大時(shí)，展向的不均勻性增強(qiáng)；在中性曲線上分支附近，當(dāng)雷諾數(shù)減小時(shí)，展向的不均勻性增強(qiáng)．

圖11 擾動速度 u在 xz平面上的等值線圖(中性曲線下分支附近)Fig.11 Contours of disturbance velocity u in the xz plane (near the lower branch of the neutral curve)

圖12 擾動速度 u在 xz平面上的等值線圖(中性曲線上分支附近)Fig.12 Contours of disturbance velocity u in the xz plane (near the upper branch of the neutral curve)

2.2.3 全局穩(wěn)定性下的中性曲線

在不同流向波數(shù)α下進(jìn)行計(jì)算，得到了如圖 13的b0在Re-α面上的等值線圖，即全局穩(wěn)定性下的中性曲線，本文計(jì)算的流向波數(shù)范圍是(0.95，1.09)．從圖中可以看到，最外層的線對應(yīng)的是 b0＝0的情況，即一般槽道流的中性曲線，并且隨著b0的增大，中性曲線逐漸向下移動．

圖13 b0在Re-α 平面上的等值線圖Fig.13 Contours of b0in the Re-α plane

3 結(jié)論及展望

本文通過數(shù)值模擬的方法對基本流在展向周期變化的槽道流的全局穩(wěn)定性問題進(jìn)行了研究，得出以下主要結(jié)論．

(1) 全局穩(wěn)定性下的擾動呈指數(shù)演化的形式．

(2) 展向幅值b對擾動起著穩(wěn)定性作用，即隨著b的增大，擾動的增長率減小，流動變得越來越穩(wěn)定.

(3) 全局中性情況下各階諧波沿法向的分布類似．

(4) 隨著 b0的增大，全局穩(wěn)定性下的中性曲線逐漸向下移動．

本文只給出了流向波數(shù)α范圍為(0.95，1.09)的中性曲線，對于其他波數(shù)可以進(jìn)一步計(jì)算，并且本文只針對展向基本波數(shù)為1.0的情況進(jìn)行了研究，對于其他波數(shù)可以進(jìn)行更深入的探討．

［1］ Huerre P，Monkewitz P. Local and global instabilities in spatially developing flows[J]. Annu Rev Fluid Mech，1990，22：473-537.

［2］ Chomaz J. Global instabilities in spatially developing flows：Non-normality and nonlinearity[J]. Annu Rev Fluid Mech，2005，37(2)：357-392.

［3］ Hammond D A，Redekopp L G. Local and global insta

bility properties of separation bubbles[J]. Eur J Mech-B/Fluids，1998，17(2)：145-164.

［4］ Zulauf M，Hart J，Leben R，et al. Local instability in a periodically forced sliced cylinder[J]. Dynamics of Atmospheres and Oceans，1996，25(2)：87-107.

［5］ Monkewitz P A，Huerre P，Chomaz J M. Global linear stability analysis of weakly non-parallel shear flows[J]. Fluid Mech，1993，251(1)：1-20.

［6］ Crouch J，Garbaruk A，Magidov D. Predicting the onset of flow unsteadiness based on global instability[J]. Journal of Computional Physics，2007，10(2)：924-940.

［7］ Tezuka A，Suzuki K. The global stability analysis of steady flow around a blunted cyclinder[C]//Proceedings of the 4th Asia Workshop on Computational Fluid Dynamics. Tokyo，Japan，2004：89-94.

［8］ Couairon A，Chomaz J. Global instability in fully nonlinear systems[J]. Physical Review Letters，1996，77(4)：4015-4018.

［9］ Dongarra J，Straughan B，Walker D. Chebyshev tau-qz algorithm method for calculating spectra of hydrodynamic stability problems[J]. Applied Numerical Mathematics，1996，22(4)：399-434.

［10］ 周 恒，趙耕夫. 流動穩(wěn)定性[M]. 北京：國防工業(yè)出版社，2004.

Zhou Heng，Zhao Gengfu. Hydrodynamic Stability[M]. Beijing：National Defence Industry Press，2004(in Chinese).

［11］ 肖紅林. 槽道湍流的大渦模擬[D]. 天津：天津大學(xué)機(jī)械工程學(xué)院，2004.

Xiao Honglin. Large Eddy Simulation of Turbulent Channel Flow[D]. Tianjin：School of Mechanical Engineering，Tianjin University，2004(in Chinese).

(責(zé)任編輯：樊素英)

Numerical Simulation of Global Stability Investigation on Channel Flow with Periodical Spanwise Variation

Kong Wei，Luo Jisheng
(School of Mechanical Engineering，Tianjin University，Tianjin 300072，China)

The global stability of the channel flow with periodical spanwise variation is investigated using numerical simulation in this paper. At first，the eigenfunction of disturbance obtained by linear stability theory is given as the initial value to calculate the time development of the disturbance. Then，the eigenvalue of global stability is obtained after fully development. The results of different spanwise variation amplitudes and neutral curves in global stability are given，and detailed analysis is made on the variation of b0with the Reynolds number and the characteristics of eigenfuctions.

global stability；periodical spanwise variation；channel flow；numerical simulation

O357.1

A

0493-2137(2015)03-0246-09

10.11784/tdxbz201307029

2013-07-11；

2013-10-31.

國家自然科學(xué)基金重點(diǎn)資助項(xiàng)目(10632050)；國家自然科學(xué)基金青年科學(xué)基金資助項(xiàng)目(11202147)；高等學(xué)校博士學(xué)科點(diǎn)專項(xiàng)科研基金(新教師類)資助項(xiàng)目(20120032120007).

孔 瑋（1986— ），男，博士研究生，kongwei203@163.com.

羅紀(jì)生，jsluo@tju.edu.cn.

時(shí)間：2013-11-22.

http://www.cnki.net/kcms/detail/12.1127.N.20131122.1101.006.html．