應(yīng)用攝動(dòng)有限體積法求解Navier-Stokes方程組

范昌勝，郭 強(qiáng)

(1,陜西廣播電視大學(xué)工程管理系，西安 710119；2.西北工業(yè)大學(xué)理學(xué)院應(yīng)用數(shù)學(xué)系, 西安 710072)

有限體積法(Finite Volume Method，F(xiàn)VM)是在上世紀(jì)八十年代末發(fā)展起來的一種新數(shù)值方法。它吸收了有限差分法(Finite Difference Method，F(xiàn)DM)和有限元法(Finite Element Method，F(xiàn)EM)的優(yōu)點(diǎn)，是一種具有對(duì)區(qū)域剖分單元特征，能適應(yīng)復(fù)雜邊界，簡(jiǎn)單靈活的離散方法，同時(shí)具有間斷解的適應(yīng)性等優(yōu)點(diǎn)的數(shù)值方法。但是，F(xiàn)VM相對(duì)FDM和FEM的弱點(diǎn)是精度不高[1]。國(guó)內(nèi)外學(xué)者經(jīng)過多年的研究，雖然發(fā)展了一些實(shí)用的有限體積格式[2-3]，但這些格式的精度都不超過2階。當(dāng)求解復(fù)雜的或者局部變化比較劇烈的流動(dòng)時(shí)，不超過2階精度的格式并不能滿足工程計(jì)算的需要。這就需要構(gòu)造一些高精度的有限體積格式。在構(gòu)造高精度有限體積格式時(shí)，關(guān)鍵是構(gòu)造控制界面上的數(shù)值通量。數(shù)值通量構(gòu)造的好壞，直接關(guān)系到格式的精度、分辨率和計(jì)算效率。目前，大多數(shù)高精度有限體積格式都是基于多基點(diǎn)利用高精度插值多項(xiàng)式來構(gòu)造控制界面上的數(shù)值通量[4-5]。高智等人為克服上述缺點(diǎn)，提出了攝動(dòng)有限體積法(Perturbational Finite Volume Method，PFVM)[6-7]。PFVM在FVM的基礎(chǔ)上，利用對(duì)流項(xiàng)、擴(kuò)散項(xiàng)之間的相互關(guān)系來提高格式精度。具體做法是，對(duì)界面質(zhì)量通量的數(shù)值進(jìn)行攝動(dòng)，即將質(zhì)量通量表示成冪級(jí)數(shù)的形式。其中，冪級(jí)數(shù)系數(shù)是由方程本身的相互關(guān)系獲得的。王利業(yè)等人分別用同位和交錯(cuò)網(wǎng)格有限體積法，模擬了4:1平板收縮流動(dòng)過程，并將兩種方法所得模擬結(jié)果進(jìn)行了吻合分析。這說明同位網(wǎng)格有限體積法不僅算法簡(jiǎn)潔、實(shí)現(xiàn)方便，結(jié)果可靠,而且容易擴(kuò)展應(yīng)用于非結(jié)構(gòu)網(wǎng)格及高維問題的模擬[8]。

PFVM與FVM都是包括積分近似和插值(重構(gòu))近似的二級(jí)近似方法。與FVM 相比，PFVM由于對(duì)質(zhì)量通量進(jìn)行攝動(dòng)而提高了插值近似的精度。在計(jì)算中，一般需要對(duì)冪級(jí)數(shù)采取截?cái)嗵幚?。假設(shè)取冪級(jí)數(shù)的前N項(xiàng)作為近似，那么，以迎風(fēng)有限體積方法(Upwind Finite Volume Method ，UFVM)為基礎(chǔ)的迎風(fēng)攝動(dòng)有限體積方法(Upwind Perturbational Finite Volume Method ，UPFVM)就是N+ 1 階插值近似、二階積分近似的混合格式；以中心有限體積方法(Center Finite Volume Method ，CFVM)為基礎(chǔ)的中心攝動(dòng)有限體積法(Center Perturbational Finite Volume Method ，CPFVM)則是具有2N+2階插值近似、二階積分近似的混合格式。劉百倉介紹了一種基于非結(jié)構(gòu)網(wǎng)格的二階迎風(fēng)有限體積離散格式，計(jì)算結(jié)果表明該二階迎風(fēng)有限體積算法有很好的收斂性和穩(wěn)定性，且松弛因子對(duì)計(jì)算結(jié)果影響很小[9]。近年來發(fā)展起來一類新的無網(wǎng)格數(shù)值方法，該類方法基于點(diǎn)近似，不需要將節(jié)點(diǎn)連成單元，克服了網(wǎng)格類方法遇到的一些困難。由于無網(wǎng)格近似函數(shù)不是多項(xiàng)式，積分不能精確計(jì)算，一般需要用到較高階的Gauss積分，從而增加了計(jì)算量。這樣就需要包括更多的控制單元，導(dǎo)致求解的代數(shù)方程非常復(fù)雜，邊界條件的處理也更加困難[10]。

采用PFVM和FVM計(jì)算方腔頂蓋驅(qū)動(dòng)流。數(shù)值結(jié)果表明，PFVM比一階迎風(fēng)、二階中心FVM的精度和分辨率高、穩(wěn)定性好。與此同時(shí)，應(yīng)用PFVM計(jì)算了4:1平板收縮流，給出了不同雷諾數(shù)下的流線圖、對(duì)稱軸的速度圖以及收縮區(qū)下游管道的速度截面圖，詳細(xì)分析了雷諾數(shù)對(duì)流場(chǎng)的影響，以及收縮區(qū)下游管道長(zhǎng)度對(duì)流場(chǎng)充分發(fā)展的影響。

1 Navier-Stokes方程組的攝動(dòng)有限體積格式

假設(shè)流體為二維定常不可壓縮流的，則控制方程的積分形式如下：

采用基于同位網(wǎng)格上的PFVM離散格式求解上述方程，用動(dòng)量插值解決速度和壓力的失耦問題。具體計(jì)算采用SIMPLEC算法。其中UPFVM離散格式如下式(CPFVM離散格式與之類似)。

(2a)

(2b)

(2c)

2 數(shù)值結(jié)果和討論

2.1 方腔頂蓋驅(qū)動(dòng)流

為了應(yīng)用PFVM求解復(fù)雜的4:1平板收縮流動(dòng)，本文首先用該方法求解具有明確提法和公認(rèn)基準(zhǔn)解的方腔頂蓋驅(qū)動(dòng)流，以檢驗(yàn)算法的有效性和程序的正確性。其中方腔長(zhǎng)寬各1個(gè)單位，方腔上邊界以恒定速度u=1向右運(yùn)動(dòng)，其它邊界上的速度值為零，整個(gè)腔內(nèi)流體在粘性作用下隨之流動(dòng)。

圖1給出雷諾數(shù)Re=1 000和5 000時(shí)，用PFVM計(jì)算二維方腔頂蓋驅(qū)動(dòng)流的流線圖。圖2為Re=1 000和5 000時(shí)FVM和PFVM計(jì)算的方腔垂直中心線上u速度線和水平中心線上v速度線。在計(jì)算時(shí)，UPFVM和CPFVM的重構(gòu)近似精度都取為8階。

從圖1中發(fā)現(xiàn)，當(dāng)雷諾數(shù)較小時(shí)，左下角和右下角出現(xiàn)明顯的一級(jí)渦。隨著雷諾數(shù)的增大，右下角出現(xiàn)明顯的二級(jí)渦，左上角出現(xiàn)規(guī)則的旋渦。并且原始渦的渦心位置越來越接近方腔的幾何中心。在不同雷諾數(shù)下方腔中各級(jí)渦的渦心位置以及渦心處的流函數(shù)值見表1.這與文獻(xiàn)[11]中的結(jié)果相吻合。

圖1 流線圖
Fig.1 Streamline diagram表1 8階UPFVM計(jì)算的渦心位置和流函數(shù)
Tab.1 8-order UPFVM calculation of the vortex center position and stream function

Re=100Re=1000Re=5000原 始 渦流 函 數(shù)渦心位置?0.103420(0.6172,0.7305)?0.115408(0.5302,0.5603)?0.116057(0.5188, 0.5313)左下一級(jí)渦流 函 數(shù)渦心位置1.735×10?6(0.0313,0.0351)2.172×10?4(0.0829,0.0729)1.37×10?3(0.0737,0.1364)左下二級(jí)渦流 函 數(shù)渦心位置——————?6.0×10?8(0.0117,0.0078)右下一級(jí)渦流 函 數(shù)渦心位置1.248×10?5(0.9453,0.0664)1.786×10?3(0.8669,0.1131)3.0214×10?3(0.8229,0.0737)右下二級(jí)渦流 函 數(shù)渦心位置———?9.317×10?8(0.9922,0.0117)?1.98×10?6(0.9828,0.0204)

從圖2中可以發(fā)現(xiàn)，在低雷諾數(shù)(Re=1 000)時(shí)，一階UFVM的計(jì)算結(jié)果和基準(zhǔn)解[11]的差距較大，而UPFVM和CPFVM的計(jì)算結(jié)果和基準(zhǔn)解完全吻合。隨著雷諾數(shù)增大，流場(chǎng)的變化越來越劇烈，控制方程的非線性也越來越強(qiáng)。在較高雷諾數(shù)(Re=5 000)時(shí)，一階UFVM計(jì)算結(jié)果與基準(zhǔn)解的差距仍是最大的。但此時(shí)UPFVM以及CPFVM的計(jì)算結(jié)果與基準(zhǔn)解相吻合。同時(shí)，與CPFVM的結(jié)果相比，UPFVM的計(jì)算結(jié)果更好。這與文獻(xiàn)[6]中的結(jié)論一致，它充分說明了PFVM具有高精度高分辨率的特性。

3.2 4:1平板收縮流

4:1平板收縮流能夠反映流體經(jīng)復(fù)雜剪切和拉伸變形后的各種特性，并且在聚合物加工過程中的擠出、注塑成型、石油的輸送以及流體機(jī)械等領(lǐng)域有廣泛的應(yīng)用[12-15]。所以本文選用具有高精度高分辨率特性的PFVM來數(shù)值模擬4:1平板收縮流的流動(dòng)過程以及各種物理現(xiàn)象。4:1平板收縮流的模擬區(qū)域如圖3所示。網(wǎng)格劃分[16]在計(jì)算中取L1=5H2，L2=12H2，H1/H2=4，其中邊界條件設(shè)定如下：

(1)入口處的流速為：

(2)平板固壁采用無滑移條件，

即u=0，v=0.

圖2 中軸速度
Fig.2 Shaft speed

圖3 4:1平板收縮流模擬區(qū)域示意圖
Fig.3 Diagram of 4: 1 flat contraction flow simulation area

圖4是不同雷諾數(shù)下，由PFVM計(jì)算的流線圖。就流場(chǎng)流動(dòng)而言，在靠近壁面且遠(yuǎn)離拐角的區(qū)域，流動(dòng)以剪切為主；在接近中軸線的位置，主要是拉伸流動(dòng)；而在上游拐角區(qū)域，流動(dòng)更為復(fù)雜，以旋轉(zhuǎn)為主；在收縮入口處，流場(chǎng)突然收縮，導(dǎo)致下游的流速增大。由圖中發(fā)現(xiàn)，隨著雷諾數(shù)的增大，回流區(qū)的主渦逐漸由豎直向水平旋轉(zhuǎn)。即：在雷諾數(shù)較小時(shí)，主渦在y方向的再附著長(zhǎng)度大于在x方向的再附著長(zhǎng)度，但隨著雷諾數(shù)的增大，主渦在y方向的再附著長(zhǎng)度小于x方向的再附著長(zhǎng)度。另外，主渦整體大小隨著雷諾數(shù)增大而呈先減小后增大的趨勢(shì)，并且當(dāng)雷諾數(shù)增大到一定值時(shí)，在收縮區(qū)的下游出現(xiàn)一個(gè)扁長(zhǎng)的次級(jí)渦，此渦在x、y方向的再附著長(zhǎng)度隨著雷諾數(shù)增大而增大，且在x方向增長(zhǎng)速度比在y方向增長(zhǎng)速度快的多。

圖5是不同雷諾數(shù)下中軸線上u速度隨軸向距離變化的對(duì)比圖。從圖5(a)中可以看出，當(dāng)L2=12H2時(shí)，雷諾數(shù)Re<100時(shí)的流場(chǎng)已充分發(fā)展。但當(dāng)雷諾數(shù)Re≥100時(shí)，流場(chǎng)還沒有完全充分發(fā)展，這就需要將下游收縮管道長(zhǎng)度延長(zhǎng)。從圖5(b)中發(fā)現(xiàn)，雖然將下游管道延長(zhǎng)到50H2，也只有雷諾數(shù)Re=100時(shí)流場(chǎng)充分發(fā)展。雷諾數(shù)Re=400時(shí)流場(chǎng)依舊沒有充分發(fā)展，直到將下游管道延長(zhǎng)到75H2時(shí)，流場(chǎng)才得到充分發(fā)展?？偟膩碚f，隨著雷諾數(shù)的增大，從收縮處到流場(chǎng)充分發(fā)展處的距離變大，并且這一變化相當(dāng)劇烈。

圖4 流線圖
Fig.4 flow chart

圖5 中軸速度圖
Fig.5 diagram of axis speed

圖6是下游收縮管道不同x處截面的速度分布圖。從圖可以看出，當(dāng)雷諾數(shù)較小時(shí)，收縮處u速度和出口處u速度相差不大，也就是說，流場(chǎng)在收縮處經(jīng)復(fù)雜的剪切和拉伸變形后，能夠在收縮管道下游很短的距離內(nèi)就達(dá)到充分發(fā)展；當(dāng)雷諾數(shù)較大時(shí)，收縮處u速度隨y的變化非常小，只是在上板壁處有很大的梯度，并且和出口處的u速度相差也很大。同時(shí)，在收縮管道下游達(dá)到充分發(fā)展的距離也越來越大。也就是說，雷諾數(shù)較大時(shí)收縮處對(duì)流動(dòng)的剪切和拉伸變形使流場(chǎng)發(fā)生劇烈變化，且這一變化的影響區(qū)域很大，以至于在收縮區(qū)域下游的很長(zhǎng)距離內(nèi)流場(chǎng)都沒有達(dá)到充分發(fā)展。

圖6下游收縮管道不同橫截面u速度分布
Fig.6 The velocity distribution of different cross-sections of the downstream contraction pipelines

3 結(jié)論

采用PFVM和FVM計(jì)算方腔頂蓋驅(qū)動(dòng)流，并且將PFVM推廣到4:1平板收縮流，數(shù)值結(jié)果表明：

(1)PFVM既保留了一階迎風(fēng)、二階中心格式基點(diǎn)少的優(yōu)點(diǎn)，又具有高精度和高分辨率的特性。

(2)在4:1平板收縮流中，隨著雷諾數(shù)的增大，收縮區(qū)域上游拐角的主渦呈先減小后增大的趨勢(shì)。并且渦逐漸由豎直方向向水平方向旋轉(zhuǎn)。

(3)在4:1平板收縮流中，當(dāng)雷諾數(shù)達(dá)到一定值后，收縮區(qū)下游會(huì)出現(xiàn)一個(gè)扁長(zhǎng)的渦。此渦隨著雷諾數(shù)的增大而增大，并且會(huì)變的越來越“扁”。

(4)流動(dòng)達(dá)到充分發(fā)展所需的收縮管道長(zhǎng)度，隨著雷諾數(shù)的增大而逐漸變長(zhǎng)，并且變化非常劇烈。