微小擾動下中等長高比腔體內行進波對流的高精度數值模擬

王 濤， 葛永斌

（1.北方民族大學 預科教育學院，寧夏 銀川 750021；2.寧夏大學 數學計算機學院，寧夏 銀川 750021）

微小擾動下中等長高比腔體內行進波對流的高精度數值模擬

王 濤1， 葛永斌2

（1.北方民族大學 預科教育學院，寧夏 銀川 750021；2.寧夏大學 數學計算機學院，寧夏 銀川 750021）

文章采用高階緊致有限差分格式，通過二維流體力學擾動方程組的數值模擬，研究了具有較弱Soret效應下，附加一個微小的溫度擾動作為擾動源的中等長高比腔體內混合流體對流系統中時空結構的發展；它在經歷了瞬態的對傳波、調制對傳波，最終演變為穩定的定常行進波狀態；并進一步研究了在中等長高比腔體內混合流體對流中從調制對傳波向定常行進波的時空演化及其隨瑞利數的變化規律。

Soret效應；對傳波；調制對傳波；定常行進波；高階緊致格式

各種非平衡線性系統中的時空演化現象，是近年來科學研究的熱點問題之一。其中混合流體（如水和酒精，He3與 He4）Rayleigh－Benard對流為良好穩定性、分岔、復雜時空形態和湍流的非平衡非線性系統，該系統試驗簡單、易于控制，且描述系統的精確方程是已知的，便于理論分析。Rayleigh－Benard對流是當作用在上下平板間的溫度差超過一個臨界值（臨界瑞利數）時系統就會失穩，流體將從運動無結構的熱傳導狀態轉變為格子狀的靜止的對流渦卷。如果2種液體混合，就會形成混合流體對流系統，溫度場和濃度場的耦合效應會改變系統的穩定及分岔特性。垂直方向的溫度梯度可以誘導一個濃度梯度產生的現象，被稱作Soret效應，可用分離比φ來表示，φ的大小直接影響著分叉的形態與非線性特性。

就研究方法而言，早期主要以實驗研究和理論分析為主，實驗揭示了許多行進波對流斑圖［1－2］，近年來，數值模擬已成為重要的研究手段。文獻［3－6］對封閉腔體內的混合流體對流進行了研究，發現了對傳波（CPW）、雙局部行進波（DLTW）、具缺陷的混合流體行進波對流斑圖，探討了具有強Soret效應的混合流體Undulation行進波對流斑圖的動力學特性。文獻［7］對混合流體Rayleigh－Benard對流的線性穩定性進行了分析。文獻［8］給出了雙流體流動中脈沖擾動的時空演化。文獻［9］給出了水平流作用下的混合流體行進波對流的時空演化。

對于微小擾動的混合流體，文獻［10］討論了微小擾動下的線性階段到非線性階段的成長。另外，在數值模擬中均采用低精度格式進行計算，但要對系統結構成功地進行模擬，并得到較高精度的求解結果，對計算格式的求解效率及精度提出更高的要求。因此，采用高精度緊致格式來研究混合流體系統，探討混合流體對流運動的動力學特性及參數依賴性等是必要的。文獻［11－12］采用高精度格式數值模擬了長腔體（長高比Γ＝46）的矩形腔體內混合流體對流，確定出存在局部行進波和雙局部行進波的瑞利數的范圍。但目前采用高精度差分格式數值模擬長高比Γ＝12（即中等長高比）的矩形腔體內具較弱強度的Soret效應（即φ＝－0.11的對流系統）的報道不多。因此，本文基于對流系統的流體力學擾動方程組的高精度數值模擬，研究φ＝－0.11的具有較弱強度的Soret效應的對流系統中的中等長高比Γ＝12腔體內的混合流體對流情況，揭示了在腔體內混合流體對流中從調制對傳波向定常行進波的時空演化及其隨瑞利數R的變化規律。

1 數學模型

1.1 基本方程組

考慮一個底部加熱的四周為固壁的方腔內充滿混合流體的物理模型。如果上部平板的溫度保持不變，當下部平板的溫度升高到某個臨界數值時，在2個平板之間對流運動發生，如圖1所示，水平放置的矩形腔體，記Γ＝Lx／d，一薄層雙流體混合物封入矩形腔體之中，例如酒精與水的混合物的水平層，使流體處于一個均勻的重力場中，g＝－gez，g為重力加速度，方向向下，ez表示z方向的單位矢量。上下板之間存在一個正的溫度差ΔT＝Thot－Tcold，由瑞利數來表征。當R引發的浮力項超過極限時對流運動就開始發生。

本文認為對流以整齊的平行滾動形式出現，對流運動的斑圖隨上下板間的溫度差的變化而變化。在正交的x、z軸所在的平面上描述二維對流運動。

圖1 對流模型示意圖

對于該流動問題，可用二維流體動力學方程組來描述。假設坐標原點為底板與左側壁的交匯處，x軸向右為正，z軸向上為正。所有的幾何尺寸用流體層厚度d表示，時間t為d2／ν，擾動速度u，w為ν／d，擾動溫度θ為kν／（agd3），擾動濃度流束η為kνγ2D／（agd3），壓力p為ν2／d2無因次化，考慮了Soret效應的流體力學擾動方程組為：

在Boussinesq近似假設下，下標0表示相應物理量的平均值，T和C為溫度和濃度偏離參考值T0和C0小的偏差，質量密度的狀態方程為：

其中熱引起的體積膨脹系數α及濃度變化引起的體積膨脹系數β分別定義為：

α的作用是單調的，因為不論液體或氣體，溫度升高時密度都會減小，α始終大于0；而β的作用卻不同，當β＜0時，對流系統類似于純流體的情況為定常流動；當β＞0時，對流系統會出現行進波狀態、局部行進波及定常流動。

1.2 渦量－流函數形式的擾動方程組

流函數φ和渦量ω的定義為：

對（1）～（3）式進行化簡得：

其中，（6）式為渦量輸運方程，（7）式為流函數Poisson方程。方程（4）～（7）組成了在 Boussinesq假設下，二維非定常不可壓縮二成分混合流體的渦量－流函數形式的基本控制方程組。

1.3 邊界條件

為了求解方程組必須給出合理的邊界條件。所有壁面都是固定的，速度在壁面上為0，濃度流束在壁面上是不可穿透的，故

由于在上下邊界上溫度是固定的，溫度的擾動量可表示為：當z＝0，1時，θ＝0。

對于限定在矩形容器中的混合流體，容器是絕熱的，故側向邊界條件可表示為：

2 高精度格式

本文給出溫度方程、濃度方程和渦量方程的統一表達式為：

其中，f可表示ω、θ、η；q表示源項。

對于（8）式的求解，本文對時間方向采用三階R－K公式［13］離散，在空間方向采用四階精度的迎風格式［14］離散。

對流函數Poisson方程（7）的求解，采用四階緊致差分格式［15］，該格式為：

邊界的具體三階計算格式如下。

擾動溫度上下邊界為：

擾動溫度左右邊界為：

擾動濃度上下邊界為：

擾動濃度左右邊界為：

計算渦邊界值時需要渦邊界本身的值，故本文采用迭代求解法。

渦左邊界為：

同理可構造出其他3個邊界的迭代格式。

3 數值模擬及討論

給定初始擾動微小擾動波的波長為2倍的流體層厚度的平行滾動，其滾動的微小振幅包絡線為稍微偏離腔體中心的高斯分布函數，本文所模擬的是長高比Γ＝12的空腔，給定初始擾動波的振幅在x＝8處達到最大值。本文給出初始擾動波的數學表達式為：

計算采用均勻網格h＝Δx＝Δz＝1／10，時間步長Δt＝5×10－4，計算中采用的其他流體參數［4］分別取φ＝－0.11，Pr＝18，L＝0.015。為了研究方便，采用相對瑞利數r，即r＝R／Rco，其中，Rco為純流體對流時的臨界值（Rco＝1 708），R為混合流體對流時的瑞利數。

在求解描述二成分混合流體對流系統的渦量－流函數形式的擾動方程組時，采用內外兩重循環，在每個時間層上對渦量、擾動溫度、擾動濃度及擾動速度采用內循環步進求解。當前后2次迭代值之差小于給定的控制收斂準則，則內循環結束。在時間方向采用外循環，當推進到需要計算的時刻程序終止。

本文以對流臨界值r＝1.20處開始計算，在給定的微小振幅的初值下，將計算結果可視化，觀測到對傳波（Counter Propagating Wave，簡稱CPW）狀態圖。取腔體1／2處的溫度沿水平方向的分布形狀。當r＝1.20時，腔體1／2高度處擾動溫度場隨時間的演化如圖2所示，其中橫軸表示空腔長度，為12倍的空腔高度，縱軸表示時間的變化，各條曲線的時間間隔為dt＝1.25。

由圖2a可看出，存在一個瞬態的CPW狀態，在CPW中心部的駐波在改變振動相位180°時，其兩側各產生一個與原駐波同相位的行進波。文獻［2］在φ＝－0.12條件下的實驗也說明在CPW狀態，對流振幅在不斷成長，始終沒有達到一個穩定的狀態。

對流系統進入CPW狀態后，出現了左行進波和右行進波競爭的現象，在該參數條件下，這種CPW狀態是瞬態的，是極其不穩定的，在競爭階段的開始，由于駐波到左右壁面的距離不等，造成了左右行進波到達左右壁面的時間不同，在左右行進波到達左右壁面后，反射回來的時間也不同，這樣隨著時間的推進，兩列波進行疊合，因而左右行進波的振幅隨時間交替增長，導致了CPW駐波中心位置的移動。CPWS（Counterpropagating Wave Source）包括一個反傳播波和一個駐波，開始向一邊移動，向右傳播，CPWS的運動引起了CPW振幅的調制。

t＝37.5～57.5時刻腔體1／2處的溫度場的時間演化圖如圖2b所示。圖2b表明，左行進波的振幅明顯大于向右傳播的行進波的振幅，并且CPWS向右移動，從虛線的位置可以看出駐波已移至x＝9處。

t＝137.5～157.5時刻腔體1／2處的溫度場的時間演化圖如圖2c所示。圖2c表明，右行進波的振幅明顯大于向左傳播的行進波的振幅，并且CPWS向左移動，從虛線的位置可以看出駐波已移至x＝3處。圖2b和圖2c為2種類型的調制對傳行進波（Modulated Counterpropagating Wave，簡稱 MCPW）。

在圖2c的基礎上，經過長時間的計算后取值，其溫度場的變化如圖2d所示。由圖2d可以看出，擾動溫度場不向左或向右傳播了，而是在原有的位置連續地滾動，隨時間不再變化。在這種狀態下，對流的渦卷不再向某個方向移動，在固定的位置上不斷地旋轉，行進波向左或向右的傳播速度為0?？傊?，獲得了很好的定常狀態。

圖2 r＝1.20時腔體1／2高度處擾動溫度場隨時間演化

本文對φ＝－0.11，Pr＝18，L＝0.015下，Γ＝12的腔體內酒精與水的混合流體在實驗可控制的溫度下，進行了數值模擬研究，得到與實驗結果［4］相吻合的結論。

通過數值模擬的研究發現，瑞利數對腔體內流體的對流影響非常大。通過計算得出，只有R＝1.20×1 708時出現CPW 狀態，系統發展并穩定在SOC（stationary overturning convection）狀態。當R≥1.20×1 708時，計算了多個R值的對流情況，發現在這些瑞利數下，系統的對流振幅增大到其SOC狀態時的水平，直接進入SOC狀態，即R＝1.20×1 708為臨界值時，當R超過該臨界值時，經過長時間的計算取值，行進波過渡成定常行進波對流狀態。并且R越大，進入SOC狀態的時間越早。

為了表征行進波對流特征，本文采用垂直擾動流速最大值這個特性數，該特性數｜w｜max可直接表示對流振幅的大小。R＝1 708×1.50和R＝1 708×5.0下對流振幅隨時間的發展情況如圖3所示。

不同瑞利數下MCPW狀態轉變為SOC狀態時對應的時間見表1所列。

從圖3可看出，隨著時間的推進，對流振幅穩定在0.497和1.379左右，即均穩定在一個固定的值，對流達到SOC狀態。

圖3 垂直流速最大值隨時間的變化

表1 不同瑞利數時MCPW轉變為SOC狀態時的時間

在R≥1.20×1 708時，本文計算了多個R值的對流情況，發現在瑞利數下，系統的對流振幅增大到其SOC狀態時的水平，且R越大，對流振幅隨著瑞利數R的增大而增大。不同瑞利數下到達SOC狀態時相對應的對流振幅見表2所列。

表2 不同瑞利數時達到SOC狀態時對應的對流振幅

4 結束語

本文數值模擬了具有分離比為－0.11時混合流體對流系統的行進波對流情形，討論了附加一個微小的溫度擾動作為擾動源，小擾動成長發展，它在經歷了瞬態的對傳波，調制對傳波，最終演變成為穩定的定常行進波狀態，在與文獻相同的參數下，模擬得到的結果與文獻結果吻合得很好，表明本文所采用的方法適用于復雜流動問題的計算。另外，在本文參數下，隨著瑞利數R的增大，系統進入定常行進波的時間隨著瑞利數R的增大而減小，在該參數下，對流振幅整體上隨著瑞利數R的增大而增大。

［1］Heinrichs R，Ahlers G，Cannel D S.Traveling waves and spatial variation in the convection of a binary mixture［J］.Physical Review A，1987，35（6）：2761－2764.

［2］Moses E，Fineberg J，Steinberg V.Multistability and confined traveling wave patterns in a convecting binary mixture［J］.Physical Review A，1987，35（6）：2757－2760.

［3］Ning L Z，Harada Y，Yahata H.Formation process of the traveling wave state with a defect in binary fluid convection［J］.Progress of Theoretical Physics，1997，98 （3）：551－566.

［4］Ning L Z，Qi X，Yuan Z，et al.A counter propagating wave states with a periodically horizontal motion of defects［J］.Journal of Hydrodynamics，2008，20（5）：567－573.

［5］寧利中，周 洋，王思怡，等.Pioseuille－Rayleigh－Benard流動中的局部行進波 ［J］.水動學研究與進展：A輯，2010，25（3）：299－306.

［6］寧利中，齊 昕，王思怡，等.具有大負分離比的混合流體局部行進波對流 ［J］.力學季刊，2010，31（2）：194－200.

［7］李國棟，田飛飛，寧利中，等.混合流體Rayleigh－Benard對流的線性穩定性分析 ［J］.水動學研究與進展：A輯，2010，25（4）：446－452.

［8］胡 軍，尹協遠.雙流體 Poiseuille－Rayleigh－Benard流動中脈沖擾動的時空演化 ［J］.中國科學技術大學學報，2007，37（10）：1267－1270.

［9］Li G D.Traveling－wave patterns in binary mixtures convection with through flow ［J］.Prog Phys，2001，106（2）：293－313.

［10］周 倩，寧利中，張淑蕓，等.混合流體對流中擾動的成長［J］.電網與清潔能源，2011，27（2）：82－85.

［11］王 濤，田振夫，葛永斌.長腔體內混合流體行進波對流的高精度數值模擬［J］.水動學研究與進展：A輯，2011，26（1）：41－47.

［12］王 濤，葛永斌.雙局部行進波對流的高精度數值模擬［J］.合 肥 工 業 大 學 學 報：自 然 科 學 版，2012，35（6）：842－847.

［13］Shu C W，Osher S.Efficient implementation of essentially nonoscillatory shock capturing schemes［J］.Journal of Computational Physics，1988，77：439－471.

［14］Tian Z F，Li Y A.Numerical solution of the incompressible navier－strokes equation with a three－point fourth－order upwind compact difference schemes［C］／／Proceedings of 4th International Conference on Nonlinear Mechanics，Shanghai，China，2002：942－946.

［15］Gupta M M.A fourth－order Poisson solver［J］.Journal of Computational Physics，1985，55（1）：166－172.

High accuracy numerical simulation of traveling wave convection in an intermediate－aspect－ratio rectangular cell with small perturbation

WANG Tao1， GE Yong－bin2
（1.College of Preparatory Education，Beifang University of Nationalities，Yinchuan 750021，China；2.College of Mathematics and Computer Science，Ningxia University，Yinchuan 750021，China）

By using the high order compact finite difference scheme，and through the numerical simulations of two－dimensional hydrodynamic perturbation equations，the evolution of spatio－temporal patterns of the traveling wave convection in fluid mixtures in an intermediated－aspect－ratio rectangular cell with small temperature perturbation under the weakly Soret effect is performed.The result shows that counterpropagating waves are generated and modulated until the system evolves to spatially modulated counterpropagating and finally to a stable stationary overturning convection.Based on this，the Rayleigh effect on the course of spatio－temporal evolution from modulated counterpropagating wave to stationary overturning convection is discussed.

Soret effect；counterpropagating wave；modulated counterpropagating wave；stationary traveling wave；high order compact scheme

O357.1

A

1003－5060（2012）11－1569－05

10.3969／j.issn.1003－5060.2012.11.030

2012－04－10；

2012－07－28

國家自然科學基金資助項目（11061025）；教育部科學技術研究重點資助項目（210239）和北方民族大學科學研究基金資助項目（2010Y038）

王 濤（1978－），女，寧夏銀川人，北方民族大學講師；

葛永斌（1975－），男，寧夏青銅峽人，博士，寧夏大學教授，碩士生導師.

（責任編輯 閆杏麗）