劉豫龍，徐明海，劉德緒，龔金海，莫 瓊

(1.中國石油大學(xué)儲(chǔ)運(yùn)與建筑工程學(xué)院，山東青島266580;2.中國石化勝利油田分公司勝利發(fā)電廠，山東東營257087;3.中國石化中原石油勘探局勘探設(shè)計(jì)研究院，河南濮陽457001)

不可壓縮流體流動(dòng)與傳熱數(shù)值模擬具有廣泛的工業(yè)應(yīng)用背景，Patankar［1］提出交錯(cuò)網(wǎng)格上SIMPLE(semiimplicit pressure linked equations)方法，成功解決了沒有壓力控制方程中的問題。然而隨著計(jì)算問題的復(fù)雜程度增加，同位網(wǎng)格方法逐漸顯示了其優(yōu)越性，Rhie等［2］提出動(dòng)量插值方法克服壓力速度失耦技巧后，同位網(wǎng)格方法在非結(jié)構(gòu)化網(wǎng)格和三維模擬中得到了廣泛應(yīng)用。Majumdar［3］發(fā)現(xiàn)了 Rhie＆Chow的動(dòng)量插值所得結(jié)果與松弛因子有關(guān)的問題并提出了改進(jìn)方案。Yu 等［4］和 Peric［5］觀察到其不僅與松弛因子有關(guān)，對(duì)于非穩(wěn)態(tài)問題還與時(shí)間步長有關(guān)，各自提出了與時(shí)間步長和松弛因子無關(guān)的界面速度插值格式。Date［6］提出了所有界面速度均采用算術(shù)平均、通過計(jì)算節(jié)點(diǎn)多維平均壓力以克服壓力鋸齒波的方法。插值過程中壓力梯度和系數(shù)應(yīng)作為乘積計(jì)算，而幾乎所有的動(dòng)量插值都分別插值到界面上，相當(dāng)于兩個(gè)數(shù)乘積的插值等于各自插值的乘積，顯然數(shù)學(xué)上這是不成立的。筆者針對(duì)動(dòng)量插值數(shù)學(xué)上的不足，提出一種簡單的修正格式，采用與原始的SIMPLE算法非常相似的推導(dǎo)過程，得到不可壓流體流動(dòng)的一種簡單的同位網(wǎng)格SIMPLE算法。

1 壓力速度耦合模型

不可壓縮流體流動(dòng)的數(shù)學(xué)模型為Navier-Stokes方程組:

在圖1所示的P控制容積上離散質(zhì)量和動(dòng)量方程，得

圖1 P控制容積示意圖Fig.1 Schematic diagram of control volume P

u、v、p布置在同一位置的弊端在于不能檢出壓力、速度和衰減掉壓力鋸齒分布，原因在于界面速度和壓力采用線性平均，質(zhì)量守恒方程中不出現(xiàn)本點(diǎn)的速度，而動(dòng)量方程中不出現(xiàn)本點(diǎn)的壓力。Rhie＆Chow提出了界面速度計(jì)算模型(以u(píng)速度分量為例):

式(7)實(shí)際上是人為地引入了界面壓力梯度，在壓力分布為線性的情況下，它實(shí)際上為速度的線性插值(此時(shí)，大括弧內(nèi)各項(xiàng)之和為零);壓力不是線性分布的情況下，相鄰單元的壓力梯度平均值不等于其共享界面處的壓力梯度，大括號(hào)中的項(xiàng)不為零，即 使 收 斂， 誤 差 也 不 能 消 除。Date［6］和Mujumdar［3］提出了整個(gè)動(dòng)量方程插值的形式:

無論Rhie＆Chow的原始動(dòng)量插值還是Peric的修正以及其他最近的修正格式，在數(shù)學(xué)中都有一個(gè)明顯的不足。把節(jié)點(diǎn)的壓力梯度項(xiàng)構(gòu)造到單元邊界上時(shí)，都采用了數(shù)學(xué)上不存在的公式，即兩個(gè)數(shù)積的插值等于各自插值的積?？刂迫莘eP與控制容積E之間的界面e動(dòng)量插值速度嚴(yán)格的數(shù)學(xué)關(guān)系為

假設(shè)離散動(dòng)量方程已經(jīng)求解，獲得了速度分布u*、v*，假設(shè)滿足質(zhì)量守恒的速度和壓力分別為

離散動(dòng)量方程為

同樣略去鄰點(diǎn)速度修正的影響，得

故

單元P和E之間的界面速度修正動(dòng)量插值為

同理可得

可以看出，式(19)與式(7)和(8)在數(shù)學(xué)上都有相同的不足，但隨著迭代過程的收斂，壓力修正值逐步趨于零，即壓力修正值的梯度趨于常數(shù)，因此式(19)逐步趨于正確。另外，隨著迭代次數(shù)增加，速度修正值也趨于零，這種數(shù)學(xué)上的不成立逐漸消失，到收斂時(shí)完全消失。而式(7)與(8)的結(jié)果，即使速度和壓力場收斂，也不能克服其數(shù)學(xué)上的不足，因此只能通過逐漸加密網(wǎng)格來消除這種缺陷，當(dāng)流場內(nèi)壓力分布明顯偏離線性分布時(shí)，會(huì)得出不合理的結(jié)果。

原始的Rhie＆Chow的動(dòng)量插值另外一個(gè)問題是結(jié)果與松弛因子有關(guān)，即因其界面速度由動(dòng)量插值和速度線性插值兩部分構(gòu)成，所得結(jié)果依賴于壓縮因子?？梢钥闯觯?9)界面速度僅為節(jié)點(diǎn)速度的插值，而不是速度動(dòng)量插值和線性插值的組合，因此與松弛因子無關(guān)。這樣，克服了原始的動(dòng)量插值法的兩個(gè)缺陷。

無論動(dòng)量方程、質(zhì)量守恒方程還是其他求解變量方程離散中，用到的界面速度都采用線性插值，而不必計(jì)算假擬界面速度，也不必計(jì)算單元界面的壓力梯度，同時(shí)，也不像Date格式那樣，須計(jì)算節(jié)點(diǎn)的平均壓力，大大節(jié)省了計(jì)算時(shí)間和存儲(chǔ)空間需求，也簡化了計(jì)算編程過程，有利于向三維和非結(jié)構(gòu)化網(wǎng)格推廣。

2 算例

針對(duì)移動(dòng)頂蓋驅(qū)流場和方腔自然對(duì)流進(jìn)行了計(jì)算考核。

2.1 移動(dòng)頂蓋驅(qū)問題

移動(dòng)頂蓋驅(qū)的物理模型如圖2所示。方形區(qū)域，左右和底邊界為靜止的固體壁，頂部邊界水平和垂直速度分量分別為1和0，這樣在區(qū)域內(nèi)部會(huì)發(fā)生流體的旋轉(zhuǎn)流動(dòng)，形成幾個(gè)渦。采用多重網(wǎng)格方法求解渦量－流函數(shù)方程，得到了不同雷諾數(shù)的流場，給出了中心線上u、v速度分布表，可用于對(duì)比，同時(shí)也給出了幾個(gè)渦的位置分布情況。

圖2 移動(dòng)頂蓋驅(qū)示意圖Fig.2 Schematic diagram of top drive

圖3為移動(dòng)頂蓋驅(qū)中心線速度分布。由圖3可以看出，Ghia計(jì)算中采用渦量－流函數(shù)法(4階格式);本文計(jì)算時(shí)對(duì)流項(xiàng)采用二階迎風(fēng)，擴(kuò)散項(xiàng)采用中心差分格式，所得結(jié)果吻合較好。為了進(jìn)一步比較，驗(yàn)證了雷諾數(shù)Re=3 200的情況，結(jié)果如圖4所示。與Ghia的結(jié)果相比，吻合也較好。圖5給出了流函數(shù)等值線對(duì)比。由圖5可以看出，無論是渦的個(gè)數(shù)還是渦的位置，二者吻合都很好。

圖3 移動(dòng)頂蓋驅(qū)中心線速度分布(Re=400)Fig.3 Velocity distribution of top drive's centre line(Re=400)

圖4 Re=3200時(shí)移動(dòng)頂蓋驅(qū)中心線速度分布Fig.4 Velocity distribution of top drive's centre line with Re of 3200

圖5 流函數(shù)對(duì)比Fig.5 Comparison of stream function

2.2 方腔自然對(duì)流

方腔自然對(duì)流是數(shù)值傳熱學(xué)中經(jīng)常用到的程序和算法考核算例，其物理模型如圖6所示。頂部和底部邊界為絕熱邊界，左右兩側(cè)邊界為溫度分別是TH和TL的傳熱邊界，溫差引起流體的密度差別，導(dǎo)致其內(nèi)流體發(fā)生自然對(duì)流。

圖6 方腔自然對(duì)流示意圖Fig.6 Schematic diagram of square cavity natural convection

方腔自然對(duì)流主要考核指標(biāo)為冷邊界和熱邊界的平均對(duì)流換熱系數(shù)(努塞爾數(shù))。文獻(xiàn)2中給出了國內(nèi)外學(xué)者的計(jì)算結(jié)果。本文計(jì)算中采用81×81的網(wǎng)格進(jìn)行了計(jì)算對(duì)比，結(jié)果見表1。

表1的結(jié)果表明，雖然計(jì)算采用二階迎風(fēng)格式，計(jì)算節(jié)點(diǎn)數(shù)較少，努塞爾數(shù)平均值、最大和最小值與文獻(xiàn)結(jié)果吻合很好，努塞爾數(shù)最小值對(duì)應(yīng)的位置在高瑞利數(shù)Ra的情況下誤差大些，說明本方法是可行的。

表1 方腔自然對(duì)流典型數(shù)據(jù)與文獻(xiàn)結(jié)果對(duì)比(81×81)Table 1 Comparison between typical data and literature data in square cavity natural convection(81×81)

為了考察計(jì)算結(jié)果與松弛因子的關(guān)系，在Ra=105時(shí)，速度松弛因子分別用αu=αv=0.25和αu=αv=0.65進(jìn)行了計(jì)算，計(jì)算中壓力松弛因子取0.675，溫度松弛因子取0.8，網(wǎng)格采用81×81，比較所得冷面和熱面努塞爾數(shù)的變化。當(dāng)αu、αv分別為0.25和0.65時(shí)，Nu均為4.536 4，這表明所得結(jié)果與速度松弛因子無關(guān)。

3 結(jié)束語

采用與原始的SIMPLE算法相似的推導(dǎo)過程，克服了同位網(wǎng)格上用動(dòng)量插值計(jì)算界面速度的數(shù)學(xué)上的缺陷，得到了不可壓流體流動(dòng)的一種非常簡單的同位網(wǎng)格SIMPLE算法。用有限數(shù)值算例證明了本文提出的同位網(wǎng)格算法的有效性。所有界面速度均采用算術(shù)平均值，具有編程簡單、內(nèi)存占用少的優(yōu)點(diǎn)，同時(shí)有效地克服了計(jì)算結(jié)果與松弛因子有關(guān)的缺點(diǎn)。

［1］ PATANKAR S V.Numerical heat transfer and fluid flow［M］.New York:Hemisphere，1980.

［2］ RHIE C M，CHOW W L.A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation ［J］.AIAA Journal，1983，21:1525-1532.

［3］ MAJUMDAR S.Role of under relaxation in momentum interpolation for calculation of flow with non-staggered grids［J］.Numerical Heat Transfer，1988，13:125-132.

［4］ YU B，TAO W Q，WEI J J.Discussion on momentum interpolation method for colocated grids of incompressible flow ［J］.Numerical Heat Transfer:Part B，2002，42:141-166.

［5］ PERIC M.A finite volume method for the prediction of three-dimensional fluid flow in complex ducts［D］.London:Imperial College，University of London，1985.

［6］ DATE A W.Solution of Navier-Stokes equations on nonstaggered grid［J］.International Journal for Heat and Mass Transfer，1993，36:1913-1922.

［7］ MILLER T F，SCHMIDT F W.Use of a pressure-weighted interpolation method for the solution of the incompressible Navier-Stokes equations on a non-staggered grid system ［J］.Numerical Heat Transfer，1988，14:213-233.

［8］ TAO Wen-quan.Numerical heat transfer［M］.2nd ed.Xi'an:Xi'an Jiaotong University Press，2001:202-245.

［9］ GHIA U，GHIA K N，SHIN C T.High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method.［J］.Journal of Computational Physics，1982，48，387-411.

［10］ 王振興，徐明海.不同多邊形網(wǎng)格計(jì)算性能的比較［J］.工程熱物理學(xué)報(bào)，2009，30(1):93-95.

WANG Zhen-xing，XU Ming-hai.The computing characteristics of polygonal meshes for solving laminar Navier-Stokes equations［J］.Journal of Engineering Thermophysics，2009，30(1):93-95.

［11］ 徐明海，陶文銓.一種求解不可壓N-S方程的非結(jié)構(gòu)化網(wǎng)格方法［J］.西安交通大學(xué)學(xué)報(bào)，2005，39(1):83-86.

XU Ming-hai，TAO Wen-quan.Method on unstructured grid for incompressible flows［J］.Journal of Xi'an Jiaotong University，2005，39(1):83-86.

［12］ 徐明海，陶文銓.非均質(zhì)各向異性材料導(dǎo)熱類問題的數(shù)值求解方法［J］.石油大學(xué)學(xué)報(bào):自然科學(xué)版，1999，6(2):59-62.

XU Ming-hai，TAO Wen-quan.Inhomogeneous and anisotropic materials heat conduction numerical method［J］.Journal of the University of Petroleum，China(E-dition of Natural Science)，1999，6(2):59-62.

［13］ 徐明海，陶文銓.二維各向異性非結(jié)構(gòu)化網(wǎng)格的自動(dòng)生成與應(yīng)用［J］.西安交通大學(xué)學(xué)報(bào)，2002，36(3):221-225.

XU Ming-hai，TAO Wen-quan.Generation and application of 2D anisotropic unstructured grid［J］.Journal of Xi'an Jiaotong University，2002，36(3):221-225.

［14］ 修榮榮，徐明海，黃善波，等.網(wǎng)格單元形狀對(duì)數(shù)值計(jì)算的影響［J］.工程熱物理學(xué)報(bào)，2004，25(2):317-319.

XIU Rong-rong，XU Ming-hai，HUANG Shan-bo，et al.The effects of grid shape on the solution of N-S equations［J］.Journal of Engineering Thermophysics，2004，25(2):317-319.