， ， ， ，

1 引 言

人工壓縮法是Chorin提出來(lái)的，通過(guò)改變方程的類(lèi)型來(lái)求解不可壓縮黏性流體流動(dòng)，具有形式簡(jiǎn)單、準(zhǔn)確有效且易于實(shí)施的特點(diǎn)。求解不可壓縮流動(dòng)問(wèn)題的關(guān)鍵是解決速度場(chǎng)與壓力場(chǎng)的耦合問(wèn)題。應(yīng)用原始變量方程求解非定常不可壓縮N-S方程，已有多種解決速度場(chǎng)與壓力場(chǎng)之間耦合的辦法。作為求解原始變量方程之一的人工壓縮法的基本思路是,通過(guò)在定常不可壓縮黏性流動(dòng)N-S方程組中添加一項(xiàng)時(shí)間導(dǎo)數(shù)項(xiàng)?/?t,將其轉(zhuǎn)化為一種虛擬的非定?？蓧嚎s黏性流動(dòng)N-S方程組，如果t→∞時(shí)，能保證?/?t→0,則非定?？蓧嚎s黏性流動(dòng)逼近定常不可壓縮黏性流動(dòng)。在研究工作[7-11]中，ACM已應(yīng)用于求解不可壓縮和弱可壓縮流動(dòng)問(wèn)題，并展示了ACM形式簡(jiǎn)單、準(zhǔn)確有效及易于實(shí)施的特點(diǎn)。

目前，在求解不可壓縮流動(dòng)問(wèn)題方面，離散控制方程常用的數(shù)值方法有有限差分法、有限體積法、有限元方法和譜方法等[12,13]，處理速度與壓力耦合問(wèn)題常用的方法有SIMPLE系列算法、渦量-流函數(shù)法、求解壓力泊松方程方法和人工壓縮法等。目前還沒(méi)有發(fā)現(xiàn)將譜方法與人工壓縮法進(jìn)行結(jié)合用于求解不可壓縮流動(dòng)問(wèn)題的相關(guān)文獻(xiàn)。因此，本文結(jié)合兩種方法的特點(diǎn)，開(kāi)發(fā)了SCM-ACM用于求解不可壓縮流動(dòng)問(wèn)題，并研究了SCM-ACM求解二維不可壓縮流動(dòng)問(wèn)題的收斂特性，為不壓縮流動(dòng)問(wèn)題的數(shù)值求解提供了新的選擇。

2 人工可壓縮格式的控制方程

本文選取流動(dòng)問(wèn)題的典型算例——二維方腔內(nèi)的頂蓋驅(qū)動(dòng)流為計(jì)算對(duì)象。方腔頂蓋驅(qū)動(dòng)流(簡(jiǎn)稱(chēng)方腔流)是一個(gè)非常經(jīng)典的流體力學(xué)問(wèn)題，對(duì)其研究也非常廣泛。盡管區(qū)域比較簡(jiǎn)單，但其包含了很復(fù)雜的流體流動(dòng)特性，如主渦、次級(jí)渦、剪切流及邊界層等。同時(shí)，方腔流也是驗(yàn)證新數(shù)值方法的經(jīng)典算例。

下面給出了求解頂蓋驅(qū)動(dòng)流所涉及到的不可壓縮格式的非穩(wěn)態(tài)控制方程。

(1)

(2)

(3)

式中 人工壓縮法中的參數(shù)c為可壓縮粘性流動(dòng)聲速，Re為雷諾數(shù)，U和V分別為x和y方向上的無(wú)量綱速度，P為無(wú)量綱壓力。

方腔內(nèi)為粘性不可壓縮流體，滿(mǎn)足如下無(wú)滑移邊界條件，

y=1處：U=1,V=0;

x=0,1,y=0處：U]0,V=0。

3 數(shù)值方法

本文采用配置點(diǎn)譜方法離散空間偏導(dǎo)數(shù)，其具體離散過(guò)程包括計(jì)算區(qū)域的轉(zhuǎn)換和系數(shù)的確定。

3.1.1 計(jì)算區(qū)域的轉(zhuǎn)換

采用譜方法進(jìn)行求解，必須要求所求解變量在[0,1]之間的規(guī)則區(qū)域內(nèi)，因此計(jì)算區(qū)域轉(zhuǎn)換是必要的步驟，在直角坐標(biāo)系內(nèi)，二維問(wèn)題的區(qū)域轉(zhuǎn)換如下。

(4)

(5)

(6)

3.1.2 系數(shù)的確定

本文所采用的是一種直接的配置點(diǎn)譜方法，控制方程中的參數(shù)和系數(shù)用離散的系數(shù)矩陣代替。

假設(shè)，

(7)

式中T(ξ)為Chebyshev多項(xiàng)式，ξ=cos(πj/N),j=0,1，…，N，相應(yīng)的變量u(ξ,t)對(duì)ξ的一階導(dǎo)數(shù)用Chebyshev多項(xiàng)式離散的表達(dá)式為

(8)

式中

寫(xiě)成矩陣形式為

其中一階導(dǎo)數(shù)矩陣系數(shù)(DN)i j為

(9)

(10)

由式(10)得出其一階導(dǎo)數(shù)和二階導(dǎo)數(shù)的形式分別為

(11)

(12)

式中d(2)i,j的表達(dá)式可由d(2)i,j=d(1)i,jd(1)i,j得出。

針對(duì)離散后的矩陣形式的代數(shù)方程采用隱格式求解方法進(jìn)行求解，其中對(duì)流項(xiàng)采用顯式的Adams-Bashforth格式，擴(kuò)散項(xiàng)和壓力項(xiàng)采用隱式格式，具體離散形式如下。

(13)

(14)

(15)

式中 Δt為時(shí)間步長(zhǎng)，本文取Δt=10-5,U為速度矢量。對(duì)離散后矩陣形式的N-S采用二步直接求解法[14]即可直接求解出速度分布。

4 結(jié)果分析

選取方腔內(nèi)的平均流函數(shù)作為網(wǎng)格獨(dú)立性測(cè)試的敏感變量，并在橫縱坐標(biāo)方向選取相同數(shù)目的網(wǎng)格節(jié)點(diǎn)(10～70)，雷諾數(shù)為100時(shí)，平均流函數(shù)隨網(wǎng)格節(jié)點(diǎn)的變化如圖1所示。

可以看出，網(wǎng)格節(jié)點(diǎn)數(shù)較少時(shí)(10～30)，平均流函數(shù)值變化較大；當(dāng)網(wǎng)格節(jié)點(diǎn)數(shù)超過(guò)30時(shí)，平均流函數(shù)值變化微小，所以對(duì)網(wǎng)格獨(dú)立性檢驗(yàn)的精度影響程度很小。除此之外，不同的網(wǎng)格節(jié)點(diǎn)數(shù)對(duì)應(yīng)的計(jì)算時(shí)間也有較大差別，網(wǎng)格節(jié)點(diǎn)數(shù)到70時(shí)所用的計(jì)算時(shí)間較長(zhǎng)，約為網(wǎng)格節(jié)點(diǎn)數(shù)為40時(shí)的 4倍，表1給出了不同節(jié)點(diǎn)數(shù)對(duì)應(yīng)的計(jì)算時(shí)間。因此，綜合考慮計(jì)算時(shí)間和精度，本文選取網(wǎng)格節(jié)點(diǎn)數(shù)為40。

由于缺少實(shí)驗(yàn)結(jié)果，選取數(shù)值計(jì)算結(jié)果對(duì)SCM-ACM求解二維方腔內(nèi)頂蓋驅(qū)動(dòng)流問(wèn)題進(jìn)行代碼有效性驗(yàn)證。本文分別對(duì)Re=100和 Re=1000時(shí)，方腔內(nèi)的流場(chǎng)分布和速度分布進(jìn)行了比較和驗(yàn)證。

表1 不同網(wǎng)格節(jié)點(diǎn)數(shù)(10～70)對(duì)應(yīng)的計(jì)算時(shí)間

Tab.1 Calculation time corresponding to the number of grid nodes (10～70)

節(jié)點(diǎn)數(shù)目10203040506070計(jì)算時(shí)間/s1350240032004800120001300020000

圖1 Re=100時(shí)平均流函數(shù)隨網(wǎng)格節(jié)點(diǎn)數(shù)的變化

Fig.1 Change of average stream function with the number of grid nodes with Re=100

圖2給出了Re=100時(shí)本文的計(jì)算結(jié)果和參考文獻(xiàn)結(jié)果的比較?？梢钥闯?，本文的計(jì)算結(jié)果與Li等[14]的計(jì)算結(jié)果吻合較好，而與Rahman等[15]的結(jié)果則有輕微差異。與圖2(a，b)相比，Rahman等[15](圖2(c))計(jì)算的流場(chǎng)中心漩渦偏右上。為了更細(xì)致地檢驗(yàn)結(jié)果的準(zhǔn)確性，圖3給出了方腔中線(xiàn)x=0.5處速度U沿y方向的分布。流場(chǎng)分布(圖2)能從宏觀(guān)上檢驗(yàn)流動(dòng)分布的有效性，而速度U的分布則能更細(xì)致地檢驗(yàn)本文計(jì)算結(jié)果的準(zhǔn)確性。從中線(xiàn)x=0.5處速度U的分布可以看出，本文的計(jì)算結(jié)果與參考文獻(xiàn)的結(jié)果幾乎重合，說(shuō)明本文的計(jì)算結(jié)果具有很好的準(zhǔn)確度。

為了更好地檢驗(yàn)SCM-ACM代碼的有效性，本文給出了Re=1000時(shí)流場(chǎng)分布和速度分布的驗(yàn)證，分別如圖4和圖5所示。在流場(chǎng)宏觀(guān)分布為了檢測(cè)SCM-ACM的收斂特性，本文給出了收斂誤差隨迭代次數(shù)的變化趨勢(shì)。收斂誤差定義為后一時(shí)層所有節(jié)點(diǎn)速度值(或壓力)與前一時(shí)層所有節(jié)點(diǎn)速度值(或壓力)的最大差值，如速度U的誤差為ε(U)=max(Un+1-Un)。

圖2 Re=100時(shí)方腔內(nèi)流場(chǎng)分布

Fig.2 Flow field distribution in a square cavity with Re=100上，Re=1000時(shí)本文的計(jì)算結(jié)果(圖4(a))與參考文獻(xiàn)[16,17]的結(jié)果吻合很好。中線(xiàn)x=0.5處速度U的分布結(jié)果顯示，本文的計(jì)算結(jié)果與Erturk等[18]以及Li等[14]的結(jié)果幾乎重合，說(shuō)明在Re=1000時(shí)，本文開(kāi)發(fā)的SCM-ACM仍具有很好的計(jì)算精確度。

圖3 Re=100時(shí)速度U在x=0.5中心線(xiàn)處的分布

Fig.3 Distribution ofUatx=0.5 with Re=100

圖4 Re=1000時(shí)方腔內(nèi)流場(chǎng)分布

Fig.4 Flow field distribution in a square cavity with Re=1000

圖6給出了Re=100和1000時(shí)，速度分量U，V和壓力P的收斂誤差隨迭代次數(shù)的變化，其中縱坐標(biāo)為對(duì)數(shù)坐標(biāo)。可以看出，速度和壓力的誤差隨迭代次數(shù)的增加呈指數(shù)下降(迭代初期除外)，這種收斂特性繼承自譜方法。在迭代初期誤差的快速變化是由U，V和P給定的初值所引起的。

圖5 Re=1000時(shí)速度U在x=0.5中心線(xiàn)處的分布

Fig.5 Distribution ofUatx=0.5 with Re=1000

圖6 速度分量U，V和壓力P的誤差隨迭代次數(shù)的變化

Fig.6 Errors of the velocity componentU，Vand pressurePagainst calculation steps with Re=100 and Re=1000

5 結(jié) 論

基于SCM和ACM各自的優(yōu)點(diǎn)，本文開(kāi)發(fā)了SCM-ACM數(shù)值方法，并選用經(jīng)典的頂蓋驅(qū)動(dòng)流為測(cè)試對(duì)象，檢驗(yàn)了SCM-ACM求解不可壓縮流動(dòng)問(wèn)題的網(wǎng)格獨(dú)立性和代碼有效性，分析了SCM-ACM的收斂特點(diǎn)。研究結(jié)果顯示，本文開(kāi)發(fā)的SCM-ACM能夠有效地用于求解不可壓縮流動(dòng)問(wèn)題，并具有指數(shù)收斂的特點(diǎn)，為不可壓縮粘性流動(dòng)問(wèn)題的求解提供了一種新的選擇。

:

[1] Benwen L I，Zhang W.Direct solvers for discrete ordinates equations based on Chebyshev collocation spectral method[J].JournaloftheChemicalIndustry&EngineeringSocietyofChina，2010，61(2):296-301.

[2] Shen J，Yang X F.An efficient moving mesh spectral method for the phase -field model of two -phase flows[J].JournalofComputationalPhyscis，2009，228(8):2978-2992.

[3] Chen Y Y，Li B W，Zhang J K.Spectral collocation method for natural convection in a square porous cavity with local thermal equilibrium and non-equilibrium models[J].InternationalJournalofHeatandMassTransfer，2016，96:84-96.

[4] Sun Y S，Ma J，Li B W.Spectral collocation method for convective-radiative transfer of a moving rod with variable thermal conductivity[J].InternationalJournalofThermalSciences，2015，90:187-196.

[5] Zhang C M,Song X T,Treeviriyanupab P，et al.Delayed error verification in quantum key distribution[J].ChineseScienceBulletin，2014，59(23):2825-2828.

[6] 李明昌，梁書(shū)秀，孫昭晨.異常天氣條件下潮位過(guò)程神經(jīng)網(wǎng)絡(luò)補(bǔ)遺預(yù)測(cè)方法的研究[J].計(jì)算力學(xué)學(xué)報(bào)，2008，25(3):368-372.(LI Ming-chang，LIANG Shu-xiu，SUN Zhao-chen.Research on tide supplemental prediction using ANN methods under unusual weather[J].ChineseJournalofComputationalMechanics，2008，25(3):368-372.(in Chinese))

[7] 張廼龍，郭小明.多尺度模擬與計(jì)算研究進(jìn)展[J].計(jì)算力學(xué)學(xué)報(bào)，2010，28(s1):1-5.(ZHANG Nai-long，GUO Xiao -ming.Review on multiscale modeling and computation[J].ChineseJournalofComputationalMechanics，2010，28(s1):1-5.(in Chinese))

[8] Tanno I，Morinishi K，Satofuka N,et al.Calculation by artificial compressibility method and virtual flux method on GPU[J].Computers&Fluids，2011，45(1):162-167.

[9] Asinari P，Ohwada T，Chiavazzo E，et al.Link-wise artificial compressibility method[J].JournalofComputationalPhyscis，2012，231(15):5109-5143.

[10] Silva G，Semiao V.First- and second-order forcing expansions in a lattice Boltzmann method reproducing isothermal hydrodynamics in artificial compressibility form[J].JournalofFluidMechanics，2012，698:282-303.

[11] Clausen J R.Entropically damped form of artificial compressibility for explicit simulation of incompressible flow[J].PhysicalReviewE，2013，87(1):013309.

[12] 章?tīng)?zhēng)榮，張湘?zhèn)?二維定常不可壓縮粘性流動(dòng)N-S方程的數(shù)值流形方法[J].計(jì)算力學(xué)學(xué)報(bào)，2010，27(3):415-421.(ZHANG Zheng-rong，ZHANG Xiang-wei.Numerical manifold method for steady incompressible viscous 2D flow Navier-Stokes equations[J].ChineseJournalofComputationalMechanics，2010，27(3):415-421.(in Chinese))

[13] 陳德祥，劉 帥，徐自力,等.非穩(wěn)態(tài)流動(dòng)的隱式最小二乘等幾何方法[J].計(jì)算力學(xué)學(xué)報(bào)，2015，32(5):639-643.(CHEN De -xiang，LIU Shuai，XU Zi-li,et al.Implicit least squares isogeometric method for unsteady flow[J].ChineseJournalofComputationalMechanics，2015，32(5):639-643.(in Chinese))

[14] Li Z Q，Wood R.Accuracy analysis of an adaptive mesh refinement method using benchmarks of 2D steady incompressible lid-driven cavity flows and coarser meshes[J].JournalofComputationalandAppliedMathematics，2015，275:262-271.

[15] Rahman M M，Alim M A，Sarker M M A.Numerical study on the conjugate effect of joule heating and magnato -hydrodynamics mixed convection in an obstructed lid-driven square cavity[J].InternationalCommunicationsinHeatandMassTransfer，2010，37(5):524-534.

[16] 梅 歡.譜元方法求解不可壓縮流體流動(dòng)及流動(dòng)線(xiàn)性穩(wěn)定性分析[D].重慶大學(xué)，2012.(MEI Huan.Spectral Element Methods for Incompressible Fluid Flow and Linear Stability Analysis of Flow[D].Chongqing University，2012.(in Chinese))

[17] Ghia U，Ghia K N，Shin C T.High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method[J].JournalofComputationalPhyscis，1982，48(3):387-411.

[18] Erturk E，Corke T C，G?k??l G.Numerical solutions of 2D steady incompressible driven cavity flow at high Reynolds numbers[J].InternationalJournalforNumericalMethodsinFluids，2005，48(7):747-774.