配置點譜方法-人工壓縮法(SCM-ACM)求解不可壓縮流體流動

， ， ， ，

1 引 言

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

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

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

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

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

(1)

(2)

(3)

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

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

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

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

3 數(shù)值方法

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

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

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

(4)

(5)

(6)

3.1.2 系數(shù)的確定

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

假設(shè)，

(7)

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

(8)

式中

寫成矩陣形式為

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

(9)

(10)

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

(11)

(12)

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

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

(13)

(14)

(15)

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

4 結(jié)果分析

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

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

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

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

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

節(jié)點數(shù)目10203040506070計算時間/s1350240032004800120001300020000

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

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

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

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

圖2 Re=100時方腔內(nèi)流場分布

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

圖3 Re=100時速度U在x=0.5中心線處的分布

Fig.3 Distribution ofUatx=0.5 with Re=100

圖4 Re=1000時方腔內(nèi)流場分布

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

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

圖5 Re=1000時速度U在x=0.5中心線處的分布

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)點，本文開發(fā)了SCM-ACM數(shù)值方法，并選用經(jīng)典的頂蓋驅(qū)動流為測試對象，檢驗了SCM-ACM求解不可壓縮流動問題的網(wǎng)格獨立性和代碼有效性，分析了SCM-ACM的收斂特點。研究結(jié)果顯示，本文開發(fā)的SCM-ACM能夠有效地用于求解不可壓縮流動問題，并具有指數(shù)收斂的特點，為不可壓縮粘性流動問題的求解提供了一種新的選擇。

:

[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] 李明昌，梁書秀，孫昭晨.異常天氣條件下潮位過程神經(jīng)網(wǎng)絡(luò)補遺預(yù)測方法的研究[J].計算力學(xué)學(xué)報，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].計算力學(xué)學(xué)報，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] 章爭榮，張湘?zhèn)?二維定常不可壓縮粘性流動N-S方程的數(shù)值流形方法[J].計算力學(xué)學(xué)報，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)流動的隱式最小二乘等幾何方法[J].計算力學(xué)學(xué)報，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] 梅 歡.譜元方法求解不可壓縮流體流動及流動線性穩(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.