二維Woods－Burnett方程的穩(wěn)定性分析

魏淑惠，包福兵，朱贈(zèng)好

（1.東北石油大學(xué) 數(shù)學(xué)科學(xué)與技術(shù)學(xué)院，黑龍江 大慶 163318；2.中國計(jì)量學(xué)院 流體檢測(cè)與仿真研究所，浙江 杭州 310018）

0 引 言

氣體的平均分子自由程與流動(dòng)特征尺度的比值定義為努森數(shù)，Kn，用來表征氣體的稀薄程度。努森數(shù)越大，稀薄程度越高。飛行器重新進(jìn)入大氣層時(shí)，當(dāng)?shù)氐钠骄肿幼杂沙瘫容^大，從而引起稀薄程度較高，這是稀薄氣體動(dòng)力學(xué)的傳統(tǒng)研究方向［1－2］。近年來，微機(jī)電系統(tǒng)（MEMS）和納機(jī)電系統(tǒng)（NEMS）迅猛發(fā)展，微／納機(jī)電系統(tǒng)中流動(dòng)特征尺度較小，從而使得努森數(shù)較大，流動(dòng)達(dá)到滑移甚至過渡流區(qū)。微納機(jī)電系統(tǒng)中的微納尺度流動(dòng)存在許多與宏觀流動(dòng)不一樣的特性，日益受到國內(nèi)外學(xué)者的關(guān)注［3－7］。但是，當(dāng)稀薄程度較高時(shí)（Kn＞0.1），流體應(yīng)力和應(yīng)變之間的線性本構(gòu)關(guān)系不再成立，因此，基于線性本構(gòu)關(guān)系的Navier－Stokes方程已不能夠準(zhǔn)確描述此時(shí)的流動(dòng)特性，必須引入新的方法［8］。

1935年，Burnett從Boltzmann方程出發(fā)，采用Chapman－Enskog展開，獲得了應(yīng)力張量偏離平衡態(tài)本構(gòu)關(guān)系的二階非線性近似，得到二階近似的Burnett方程以來，采用高階的展流體力學(xué)方程來研究稀薄氣體運(yùn)動(dòng)的嘗試就一直沒有停止。1939年，Chapman和Cowling把原始Burnett方程的隨體導(dǎo)數(shù)用Euler方程替換，得到了常規(guī)Burnett方程。1949年，Grad［9］通過采用矩方法得到13矩方程。1991年，Zhong等人［10］為了克服常規(guī)Burnett方程的穩(wěn)定性問題，通過采用線性穩(wěn)定性分析的方法，在常規(guī)Burnett方程里增加了部分線性三階項(xiàng)，得到增廣Burnett方程。近年來，Burnett方程在稀薄氣體流動(dòng)方面獲得了很大應(yīng)用［11］。但是，關(guān)于Burnett方程能否正確描述偏離平衡態(tài)的稀薄氣體流動(dòng)，一直有不同意見［12］。

1993年，Woods分析常規(guī)Burnett方程［13］，發(fā)現(xiàn)Burnett方程中的部分項(xiàng)不合理，經(jīng)過分析推導(dǎo)，得到了一種修正的Burnett方程，稱之為 Woods－Burnett方程。Woods－Burnett方程被認(rèn)為能夠較好地描述處于滑移過渡流區(qū)的氣體流動(dòng)和傳熱特性而受到很多學(xué)者的重視［13－16］。但是 Woods－Burnett方程在小擾動(dòng)下不穩(wěn)定［17］，數(shù)值模擬過程中當(dāng)網(wǎng)格較精細(xì)時(shí)計(jì)算就發(fā)散，這極大地限制了Woods－Burnett方程的應(yīng)用。

1988年，F(xiàn)iscko和Chapman發(fā)現(xiàn)當(dāng)計(jì)算網(wǎng)格尺寸比氣體平均分子自由程小時(shí)，擴(kuò)展流體力學(xué)方程就變得不穩(wěn)定［18］。這個(gè)穩(wěn)定性問題早在1982年就被Bobylev發(fā)現(xiàn)了，通過線性穩(wěn)定性分析，他指出線性化二階Burnett方程對(duì)小于特定波長(zhǎng)的擾動(dòng)是不穩(wěn)定的［19］。Balakrishnan和Agarwal［17］通過 線 性 穩(wěn) 定性，發(fā)現(xiàn)一維Woods－Burnett方程是不穩(wěn)定的。包福兵和林建忠［20］比較了 Woods－Burnett和各種類型的Burnett方程，同樣發(fā)現(xiàn) Woods－Burnett方程是不穩(wěn)定的。并且發(fā)現(xiàn)，Woods－Burnett方程的臨界Kn數(shù)為0.184。

本文采用線性小擾動(dòng)理論，經(jīng)過大量推導(dǎo)，在一維穩(wěn)定分析的基礎(chǔ)上，對(duì)二維常規(guī)Woods－Burnett方程進(jìn)行分析，成功獲得了二維Woods－Burnett方程的穩(wěn)定性方程，首次得到了二維Woods－Burnett方程的穩(wěn)定性曲線和臨界Kn數(shù)。

1 一維常規(guī) Woods－Burnett方程的線性穩(wěn)定性分析

在笛卡爾坐標(biāo)下，一維非定?？蓧嚎s粘性流動(dòng)的控制方程如下：

其中，σ11是粘性應(yīng)力張量xx方向的分量，q1是熱通量x方向的分量，p是理想氣體的壓強(qiáng)，滿足狀態(tài)方程：

當(dāng)考慮常規(guī) Woods－Burnett方程時(shí)［14］，粘性應(yīng)力張量σ11和熱通量q1的具體表達(dá)式如下：

其中ωi和θi是 Woods－Burnett方程中的系數(shù)。

遵從Bobylev的線性穩(wěn)定性分析方法［19］，假設(shè)擾動(dòng)前氣體處于穩(wěn)定狀態(tài)，氣體的速度為零，密度和溫度都為常數(shù)，可表示為：ρ＝ρ0＝常數(shù)，u＝u0＝0，T＝T0＝常數(shù) 。引入小擾動(dòng)量：

其中，V＝［ρu T］T。密度、速度、溫度、位移和時(shí)間等變量按照以下形式無量綱化［19］：

通過把擾動(dòng)方程（4）代入方程（1）－（3），同時(shí)采用方程（5）形式進(jìn)行無量綱化，忽略方程中小量相乘的非線性項(xiàng)之后，可以得到以下的無量綱控制方程（忽略了無量綱量上的上劃線）：

其中：

方程（6）－（8）中的擾動(dòng)量可以假設(shè)成如下形式：

其中，λ＝α＋iβ，α和β分別代表擾動(dòng)的增長(zhǎng)系數(shù)和擴(kuò)散系數(shù)，k是擾動(dòng)波數(shù)，V1是擾動(dòng)的振幅。

將方程（9）代入方程（6）－（8）可得如下的關(guān)于ρ1，T1和u1的齊次代數(shù)方程組：

為保證這個(gè)齊次方程組存在非奇異解，必須要求這個(gè)方程的系數(shù)行列式的值為零［19］，從而可以得到以下的穩(wěn)定性特征方程（對(duì)于Maxwell氣體）：

從方程（13）可以得到λ與k的關(guān)系式。隨著擾動(dòng)波數(shù)k的變化，擾動(dòng)的增長(zhǎng)系數(shù)和擴(kuò)散系數(shù)在復(fù)平面上的變化稱為穩(wěn)定性特征曲線。圖1給出了一維常規(guī)Woods－Burnett方程的穩(wěn)定性特征曲線。從方程（9）中λ的定義可知，當(dāng)擾動(dòng)增長(zhǎng)系數(shù)α＞0時(shí)，擾動(dòng)振幅會(huì)隨著時(shí)間的增長(zhǎng)而變大，從而方程不穩(wěn)定。從圖1可以看出，隨著波數(shù)k的變化，存在著正的增長(zhǎng)系數(shù)，說明一維常規(guī)Woods－Burnett方程是不穩(wěn)定的。圖2給出了擾動(dòng)的增長(zhǎng)系數(shù)隨波數(shù)的變化曲線。從圖中可以看出，當(dāng)波數(shù)大于0.907時(shí)，開始出現(xiàn)一個(gè)正的增長(zhǎng)系數(shù)，說明方程開始不穩(wěn)定。根據(jù)無量綱參數(shù)定義，波數(shù)與Kn數(shù)存在如下的關(guān)系：

從而可以求得臨界Kn數(shù)為0.184。當(dāng)Knudsen數(shù)大于這個(gè)值時(shí)一維常規(guī) Woods－Burnett方程就會(huì)發(fā)散［19］。

圖1 一維常規(guī)Woods－Burnett方程的穩(wěn)定性曲線Fig.1 Stability curve of 1－D conventional Woods－Burnett equations

圖2 一維常規(guī)Woods－Burnett方程的增長(zhǎng)系數(shù)隨波數(shù)的變化Fig.2 Variations of attenuation coefficient with wave frequency for 1－D conventional Woods－Burnett equations

2 二維 Woods－Burnett方程的穩(wěn)定性分析

在笛卡爾坐標(biāo)下，二維非定?？蓧嚎s粘性流動(dòng)的控制方程如下：

采用一維穩(wěn)定性分析類似的方法，引入小擾動(dòng)，代入方程無量綱化之后的二維 Woods－Burnett方程（14）－（17），并去掉非線性項(xiàng)之后，可以得到如下的控制方程：

與一維分析類似，引入如下形式的小擾動(dòng)：

其中，V＝［ρu T］T。將方程（22）形式的小擾動(dòng)代入方程（18）－（21），可得關(guān)于ρ1，u1，v1和T1的齊次代數(shù)方程組，從而可以得到以下的穩(wěn)定性特征方程：

根據(jù)方程（23），可以得到圖3所示的二維Woods－Burnett方程的穩(wěn)定性曲線。從圖3可以看出，二維常規(guī) Woods－Burnett方程的穩(wěn)定性分析中，隨著波數(shù)的變化，有部分增長(zhǎng)系數(shù)大于零，這說明二維常規(guī)Woods－Burnett方程是不穩(wěn)定的。圖4給出了增長(zhǎng)系數(shù)隨著波數(shù)的變化曲線。從圖中可以看出，當(dāng)波數(shù)達(dá)到0.642時(shí)，有一個(gè)增長(zhǎng)系數(shù)變得大于零，這說明當(dāng)波數(shù)大于0.642時(shí)，方程就變得不穩(wěn)定。從而說明二維常規(guī) Woods－Burnett方程的臨界波數(shù)為0.642。而根據(jù)波數(shù)與Kn數(shù)的關(guān)系，可以得到二維常規(guī)Woods－Burnett方程的臨界Kn數(shù)為0.130。這個(gè)值小于一維常規(guī) Woods－Burnett方程的臨界Knudsen數(shù)為0.184。說明，二維常規(guī) Woods－Burnett方程更加不穩(wěn)定。

圖3 二維常規(guī)Woods－Burnett方程的穩(wěn)定性曲線Fig.3 Stability curve of 2－D conventional Woods－Burnett equations

圖4 二維常規(guī)Woods－Burnett方程的增長(zhǎng)系數(shù)隨波數(shù)的變化Fig.4 Variations of attenuation coefficient with wave frequency for 2－D conventional Woods－Burnett equations

3 結(jié) 論

Woods－Burnett方程是Boltzmann方程的二階近似，能夠描述輕微偏離熱力學(xué)平衡時(shí)的稀薄氣體流動(dòng)。但是Woods－Burnett方程在小擾動(dòng)下不穩(wěn)定，這是限制Woods－Burnett方程廣泛應(yīng)用的一個(gè)重要原因。本文通過線性小擾動(dòng)理論，在一維文定性分析的基礎(chǔ)上，首次得到了二維Woods－Burnett方程的穩(wěn)定性特征方程，并把穩(wěn)定性方程的解表示在復(fù)平面上，得到了二維穩(wěn)定性曲線.二維常規(guī)Woods－Burnett方程的臨界 Knudsen數(shù)為0.130，而一維 Woods－Burnett方程的臨界Knudsen數(shù)為0.184，這說明二維Woods－Burnett方程比一維 Woods－Burnett方程更加不穩(wěn)定。

二維穩(wěn)定性的發(fā)現(xiàn)對(duì)進(jìn)一步研究 Woods－Burnett方程具有重要意義，使得我們能夠確定二維Woods－Burnett方程在什么情況下會(huì)變得不穩(wěn)定，二維Woods－Burnett方程能夠應(yīng)用于哪些范圍內(nèi)的流場(chǎng)。如何增加 Woods－Burnett方程的穩(wěn)定性以及把各類Burnett方程應(yīng)用到各種處于稍微偏離平衡態(tài)的連續(xù)流或近連續(xù)滑移流動(dòng)中是今后研究的重點(diǎn)。

［1］沈青.稀薄氣體動(dòng)力學(xué)［M］.北京：國防工業(yè)出版社，2003.

［2］李志輝，張涵信.稀薄流到連續(xù)流的氣體運(yùn)動(dòng)論統(tǒng)一算法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2003，21（3）：255－266.

［3］林建忠，包福兵，王瑞金，等.微納流動(dòng)理論及應(yīng)用［M］.北京：科學(xué)出版社，2010.

［4］GAD－EL－HAK M.Review：flow physics in MEMS［J］.Mecan Indust，2001，2：313－341.

［5］KARNIADAKIS G，BESKOK A，ALURU N.Microflows and Nanoflows，fundamentals and Simulation［M］：Springer，2005.

［6］秦豐華，孫德軍，尹博遠(yuǎn).長(zhǎng)微管道氣體流動(dòng)的近似解析解［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2002，20（1）：49－56.

［7］沈青，蔣建政.有限長(zhǎng)度的過渡領(lǐng)域微槽道流動(dòng)——分析預(yù)測(cè)和試驗(yàn)的比較［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2008，26（2）：257－262.

［8］HO C M，TAI Y C.Micro－Electro－Mechanical－Systems（MEMS）and fluid Flows［J］.Annu.Rev.Fluid Mech.，1998，30：579－612.

［9］GRAD H.On the kinetic theory of rarefied gases［J］.Commun.Pure Appl.Math.，1949，2：331－407.

［10］ZHONG X L，MACCORMACK R W，CHAPMAN D R.Stabilization of the Burnett equations and application to hypersonic flows［J］.AIAA J.，1993，31（6）：1036－1043.

［11］AGARWAL R K，YUN K Y，BALAKRISHNAN R.Beyond Navier－Stokes：Burnett equations for flows in the continuum－transition regime［J］.Phys.Fluids，2001，13（10）：3061－3085.

［12］COMEAUX K A，CHAPMAN D R，MACCORMACK R W.An analysis of the Burnett equations based in the second law of thermodynamics［R］.AIAA，1995，95－0415.

［13］WOODS L C.An Introduction to the Kinetic Theory of Gases and Magnetoplasmas［M］，Oxford：Oxford University Press，1993.

［14］REESE J M，WOODS L C，THIVET F J P，et al.A second－order description of shock structure［J］.J.Comput Phys.，1995，117（2）：240－250.

［15］REESE J M，GALLIS M A，LOCKERBY D A.New directions in fluid dynamics：non－equilibrium aerodynamic and microsystem flows［J］.Philos.T.Roy.Soc.A.，2003，361（1813）：2967－2988.

［16］LOCKERBY D A，REESE J M，GALLIS M A.The usefulness of higher－order constitutive relations for describing the Knudsen layer［J］.Phys.Fluids，2005，17（10）：100609－1－6.

［17］BALAKRISHNAN R，AGARWAL R K.A comparative study of several higher－order kinetic formulations beyond Navier－Stokes for computing the shock structure［A］.In 37th AIAA Aerospace Sciences Meeting and Exhibit［C］，1999.Reno，NV.

［18］FISCKO K A，CHAPMAN D R.Comparison of Bur－nett，super－Burnett and Monte Carlo solutions for hypersonic shock structure［A］.in Proceedings of the 16th International Symposium on Rarefied Gas Dynamics［C］，1988.

［19］BOBYLEV A V.The Chapman－Enskog and Grad methods for solving the boltzmann equation［J］.Soviet Physics，1982，27（1）：29－31.

［20］BAO F B，LIN J Z.Linear stability analysis for various forms of one－dimensional Burnett equations［J］.Int.J.Nonlinear Sci.，2005，6（3）：295－303.