楊關(guān)旬 盧志偉

摘要： 為了探討蓋驅(qū)動(dòng)腔流場(chǎng)的影響因素和分布規(guī)律，本文提出采用氣體動(dòng)理學(xué)格式對(duì)其流場(chǎng)開展研究。論文結(jié)合蓋驅(qū)動(dòng)腔的物理模型，采用有限體積法對(duì)Boltzmann方程開展了數(shù)值過程推導(dǎo)和數(shù)值計(jì)算。給出了數(shù)值計(jì)算過程中以虛網(wǎng)格技術(shù)實(shí)現(xiàn)的絕熱無滑移邊界條件，計(jì)算了不同雷諾數(shù)下蓋驅(qū)動(dòng)方腔的流場(chǎng)和流線特征，提取了兩條正交中心線上的速度分布并分析其衰減規(guī)律，得到了方腔流場(chǎng)中的渦量特征和流線變化，分析了方腔的渦量強(qiáng)度變化與流線變化成因。結(jié)果表明，采用氣體動(dòng)理學(xué)格式計(jì)算蓋驅(qū)動(dòng)腔流動(dòng)可以詳細(xì)地捕捉流場(chǎng)細(xì)節(jié)和分布規(guī)律，雷諾數(shù)的不同對(duì)流場(chǎng)分布有著顯著的影響，氣體動(dòng)理學(xué)格式可以作為研究蓋驅(qū)動(dòng)方腔的典型計(jì)算方法。

Abstract： In order to explore the influencing factors and distribution law of the flow field in the lid-driven cavity， this paper proposes to use the gas kinetic scheme to study its flow field. Based on the physical model of the lid-driven cavity， the paper uses the finite volume method to carry out the numerical process derivation and numerical calculation of the Boltzmann equation. The adiabatic non-slip boundary conditions realized by ghost cell technology in the numerical calculation process are given. The flow field and streamline characteristics of the square cavity lid-driven under different Reynolds numbers are calculated， and the velocities on two orthogonal center lines are extracted. Distribution and analysis of its attenuation law， the vorticity characteristics and streamline changes in the square cavity flow field are obtained， and the causes of the vorticity intensity changes and streamline changes in the square cavity are analyzed. The results show that using the gas kinetic scheme to calculate the flow in the lid-driven cavity can capture the details and distribution of the flow field in detail， and the difference in Reynolds number has a significant impact on the flow field distribution. The gas kinetic scheme can be used as a typical calculation method for studying the lid-driven cavity.

關(guān)鍵詞： 氣體動(dòng)理學(xué)格式;蓋驅(qū)動(dòng)腔;有限體積法;數(shù)值模擬

Key words： gas kinetic scheme;lid-driven cavity;finite volume method;numerical simulation

中圖分類號(hào)：TP271+.4 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 文獻(xiàn)標(biāo)識(shí)碼：A ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?文章編號(hào)：1674-957X（2022）03-0045-03

0 ?引言

氣體動(dòng)理學(xué)格式（Gas Kinetic Scheme，GKS）是20世紀(jì)90年代Xu Kun等人提出的以Boltzmann方程及其簡(jiǎn)化模型方程開展的離散空間直接建模數(shù)值方法[1][2]。相比于N-S方程，GKS適用于多尺度，跨流域流動(dòng)，數(shù)值計(jì)算量小，能有效模擬高超速飛行器飛行面臨的高雷諾數(shù)復(fù)雜湍流問題[3]。所以，GKS從建立之初就受到了廣泛的關(guān)注，更多的學(xué)者不斷地豐富和探討這一方法。

研究現(xiàn)狀表明，氣體動(dòng)理學(xué)格式的研究更加深入，適用范圍廣泛，實(shí)用性更強(qiáng)。蓋驅(qū)動(dòng)方腔是計(jì)算流體力學(xué)（CFD）中驗(yàn)證計(jì)算方法的一個(gè)熱門問題。眾多的研究人員采用的不同的數(shù)值方法對(duì)方腔內(nèi)部的渦旋變化展開研究，Ghia等人[4]采用二維不可壓縮Navier-Stokes方程的渦量-流函數(shù)來研究耦合強(qiáng)隱式多重網(wǎng)格（CSI-MG）方法以確定高雷諾數(shù)下蓋驅(qū)動(dòng)方腔解中的有效性，此計(jì)算結(jié)果常被視為蓋驅(qū)動(dòng)方腔計(jì)算的基準(zhǔn)解。

目前眾多的研究者針對(duì)蓋驅(qū)動(dòng)方腔的研究主要還是以N-S方程為理論基礎(chǔ)，采用渦量-流函數(shù)法以有限差分的格式對(duì)蓋驅(qū)動(dòng)腔流動(dòng)進(jìn)行數(shù)值研究。本文運(yùn)用氣體動(dòng)理學(xué)格式對(duì)蓋驅(qū)動(dòng)方腔開展數(shù)值模擬，探討GKS數(shù)值模擬的數(shù)值過程。對(duì)比其他學(xué)者的計(jì)算結(jié)果，驗(yàn)證了氣體動(dòng)理學(xué)格式的適用性，計(jì)算結(jié)果的正確性。另外，對(duì)蓋驅(qū)動(dòng)方腔內(nèi)的渦旋展開了分析，揭示蓋驅(qū)動(dòng)方腔的流動(dòng)特征規(guī)律。

1 ?蓋驅(qū)動(dòng)方腔物理模型

本文研究的蓋驅(qū)動(dòng)方腔的幾何結(jié)構(gòu)如圖1所示，邊長(zhǎng)L=1.0m的二維方形腔體，方腔頂部有一速度為Ulid=1.0m·s-1的活動(dòng)平板沿水平方向平移，方腔內(nèi)密封氣體在平板的作用下產(chǎn)生運(yùn)動(dòng)，設(shè)腔體內(nèi)的氣體密度為？籽=1.0kg·m-3;運(yùn)動(dòng)粘度為ν=6.667×10-5m2·s-1。蓋驅(qū)動(dòng)腔的雷諾數(shù)Re=UlidL/v。

計(jì)算網(wǎng)格和邊界條件：蓋驅(qū)動(dòng)腔的計(jì)算過程采用均勻的笛卡爾直角網(wǎng)格來離散計(jì)算區(qū)域，邊界條件是流固邊界上的無滑移邊界條件，方腔流動(dòng)的左、右、下邊界均采用絕熱固壁無滑移邊界，上邊界采用絕熱邊界。GKS上不同的邊界條件一般采用虛網(wǎng)格技術(shù)來實(shí)現(xiàn)，一般假設(shè)流體邊界和固體邊界與網(wǎng)格界面重合，A1，A2表示流體區(qū)域內(nèi)的網(wǎng)格中心點(diǎn)，而A-1，A-2表示與A1，A2相對(duì)應(yīng)的虛網(wǎng)格中心點(diǎn)，可以得出A-1，A-2關(guān)于固體邊界分別與A1，A2對(duì)稱。為實(shí)現(xiàn)絕熱無滑移邊界條件，使得插值后邊界處的流體速度為0，即：

2 ?流場(chǎng)的數(shù)值分析

通過查閱文獻(xiàn)蓋驅(qū)動(dòng)方腔流場(chǎng)特征分布規(guī)律主要考慮雷諾數(shù)的變化關(guān)注初級(jí)渦旋和次級(jí)渦旋的變化。分別計(jì)算了雷諾數(shù)Re=100、400、1000、3200時(shí)方腔的流場(chǎng)分布情況。

2.1 正交中心線上速度分布

圖2給出了不同雷諾數(shù)下，方腔的縱向中心線上U/Ulid和橫向中心線上V/Vlid速度分布曲線。當(dāng)Re=100時(shí)，橫向中心線上V/Vlid的速度關(guān)于方腔縱向中心線反向?qū)ΨQ;方腔的縱向中心線上U/Ulid的速度曲線斜率整體變化范圍較大，且大致在水平中心線以下位置速度陡降，雷諾數(shù)相對(duì)較小，流體的運(yùn)動(dòng)粘度大，方腔頂蓋的平板驅(qū)動(dòng)對(duì)方腔底部產(chǎn)生的作用較小，腔內(nèi)底部的氣體流動(dòng)相對(duì)緩慢。結(jié)合兩條正交中心線的速度相對(duì)比分析，可以預(yù)測(cè)，當(dāng)Re=100時(shí)，方腔底部的的流場(chǎng)流態(tài)大致關(guān)于縱向中心線對(duì)稱。當(dāng)Re=400，1000，3200時(shí)，兩條正交中心線上的速度曲線斜率在遠(yuǎn)離壁面的中間位置變化范圍很小，當(dāng)Re=1000時(shí)，這種現(xiàn)象尤其明顯，斜率幾乎不變。這表明在流場(chǎng)的內(nèi)部，速度恒定，流體呈均勻運(yùn)動(dòng)。另外，靠近方腔壁面邊界時(shí)，速度斜率反生劇烈變化的的先后次序?qū)?yīng)的雷諾數(shù)依次是Re=100，400，1000，3200，且速度斜率發(fā)生劇烈變化與到壁面的距離相關(guān)，雷諾數(shù)越大，越接近壁面，速度斜率才發(fā)生劇烈變化，這表明隨著雷諾數(shù)的增大，壁面處邊界層的厚度減小，方腔中的大多數(shù)氣體繞初級(jí)渦旋做圓周運(yùn)動(dòng)。

2.2 渦量分布和流線特征

圖3給出了不同的雷諾數(shù)下蓋驅(qū)動(dòng)腔的渦量等值線。當(dāng)Re=100時(shí)，由渦量的等值線可以得出，方腔的頂蓋附近，渦量的強(qiáng)度較大，自上而下，渦量強(qiáng)度依次衰減，腔內(nèi)的渦量等值線向右偏移。方腔的底部壁面和左側(cè)壁面渦量強(qiáng)度較小，右壁面的渦量強(qiáng)度強(qiáng)于左壁面。

當(dāng)Re=1000時(shí)，方腔內(nèi)的總體渦量強(qiáng)度與Re=100相比明顯增大，方腔的左、右底角渦量明顯增強(qiáng)。壁面對(duì)腔內(nèi)氣體的剪切力作用，方腔的左右壁面附近渦量強(qiáng)度較大。壁面到腔體內(nèi)部的過渡位置形成了一條兩邊渦量強(qiáng)度大，內(nèi)部渦量強(qiáng)度小的環(huán)形渦量等值帶，這是壁面剪切力和腔體內(nèi)部氣體做圓周運(yùn)動(dòng)共同作用的結(jié)果。在方腔的內(nèi)部，渦量強(qiáng)度差別變化不大，左右底角的渦量強(qiáng)度與方腔內(nèi)部的渦量強(qiáng)度大致相等。Re=3200時(shí)，與Re=1000相比渦量強(qiáng)度除壁面附近外變化不大，左右底角的渦量強(qiáng)度增大，且右底角的渦量強(qiáng)度變化明顯強(qiáng)于左底角。渦量等值線在方腔中心形成了完成的渦量流線和恒定的渦量強(qiáng)度。腔體內(nèi)部作圓周運(yùn)動(dòng)形成的渦旋增大，似有向壁面擴(kuò)大的趨勢(shì)。這與之前從正交中心線上特征速度分布曲線得到的現(xiàn)象是一致的。

圖4出了不同的雷諾數(shù)下的流場(chǎng)流線分布，流場(chǎng)的渦量變化是流場(chǎng)流線分布的內(nèi)部歸因，流場(chǎng)流線分布是渦量強(qiáng)度變化的外部具體表現(xiàn)。當(dāng)Re=100時(shí)，頂蓋附近流線略微擁擠，在那里出現(xiàn)最大的速度，頂蓋的平移帶動(dòng)氣體產(chǎn)生了剪切力，因此頂蓋附近流場(chǎng)的渦量強(qiáng)度大。右壁面產(chǎn)生剪切力強(qiáng)于左壁面，流場(chǎng)的初級(jí)渦旋在頂蓋附近且偏右。左右底角的渦量強(qiáng)度并無明顯變化，所以在方腔的底部附近流線基本呈對(duì)稱分布。

隨著雷諾數(shù)的增大，頂蓋驅(qū)動(dòng)產(chǎn)生的初級(jí)渦旋流線向方腔的中心附近靠近，初級(jí)渦旋流線的左上角向內(nèi)凹陷。其次，方腔的左右兩底角產(chǎn)生了次級(jí)渦旋。方腔內(nèi)流場(chǎng)的流線除左右底角外分布地更加均勻，這表明腔內(nèi)大部分氣體運(yùn)動(dòng)趨于平穩(wěn)。當(dāng)Re=3200時(shí)，流場(chǎng)分布與Re=1000相比方腔的初級(jí)渦旋流線并無明顯變化。流場(chǎng)的左右底角次級(jí)渦旋隨著渦量強(qiáng)度增大流線分布更加密集，左右底角流線分布情況與渦量強(qiáng)度變化是一致的。除此之外方腔的左上角也產(chǎn)生了次級(jí)渦旋流線。

3 ?結(jié)論

本文采用氣體動(dòng)理學(xué)格式以有限體積法的方式對(duì)蓋驅(qū)動(dòng)方腔進(jìn)行了數(shù)值計(jì)算，探討和研究了蓋驅(qū)動(dòng)腔中不同的雷諾數(shù)對(duì)流場(chǎng)的渦旋變化和流線特征的影響規(guī)律。主要結(jié)論如下：雷諾數(shù)的改變不僅改變了流場(chǎng)的運(yùn)動(dòng)粘度，直接影響是雷諾數(shù)增大使運(yùn)動(dòng)的頂蓋產(chǎn)生的壁面剪切力增強(qiáng)，頂蓋帶動(dòng)腔內(nèi)氣體運(yùn)動(dòng)的對(duì)稱性遭到了破壞，流場(chǎng)運(yùn)動(dòng)的渦量等值線和流線向右傾斜;隨著雷諾數(shù)的增大，方腔的左右底角產(chǎn)生次級(jí)渦旋，并且渦旋的渦量強(qiáng)度是隨著雷諾數(shù)的增大而增強(qiáng)的。運(yùn)用氣體動(dòng)理學(xué)格式，建立基本的數(shù)學(xué)模型過程，詳細(xì)地闡述了虛網(wǎng)格技術(shù)邊界條件的構(gòu)造，揭示了GKS在流體數(shù)值計(jì)算的一般過程。蓋驅(qū)動(dòng)方腔可以作為研究氣體動(dòng)理學(xué)格式一個(gè)很好的切入點(diǎn)，同時(shí)氣體動(dòng)理學(xué)格式為進(jìn)一步研究蓋驅(qū)動(dòng)腔提供了新的理論指導(dǎo)。

參考文獻(xiàn)：

[1]Xu K. Numerical Navier–Stokes solutions from gas-kinetic theory[J]. Comput. Phys.， 1994， 114（1）： 9-17.

[2]Xu K. A gas-kinetic BGK scheme for the Navier–Stokes equations and its connection with artificial dissipation and Godunov method[J]. Comput. Phys.， 2001， 171（1）： 289-335.

[3]李詩一，張潮，譚爽，等.氣體動(dòng)理學(xué)格式及其在再入問題中的應(yīng)用[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2018，36（05）：885-890.

[4]Ghia U，Ghia K N，Shin C T. High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J]. Academic Press，1982，48（3）：387-411.