林 超 李廷秋 林澤華 馬 麟 辛建建

(武漢理工大學(xué)交通學(xué)院 武漢 430063)

0 引 言

對流占優(yōu)流動問題是計算流體力學(xué)中的一個經(jīng)典問題，為此學(xué)者們發(fā)展了很多對流離散高精度格式，其中大多數(shù)格式基于結(jié)構(gòu)網(wǎng)格有限體積法，在界面處變量值的構(gòu)造涉及多個網(wǎng)格點的插值，而且通常與上游迎風(fēng)節(jié)點有關(guān)，這很容易引起數(shù)值發(fā)散.在科學(xué)與工程問題中，格式的單調(diào)性對于求解偏微分方程十分重要，因為單調(diào)格式能夠有效抑制數(shù)值振蕩，避免產(chǎn)生非物理解.學(xué)者們對此進行了深入研究，MUSCL格式(monotonic upstream-centered scheme for conservation laws)是一種具有二階精度的單調(diào)迎風(fēng)格式，為了保證數(shù)值求解的單調(diào)性，Van[1]采用了分段多項式計算界面通量，并在其中引入了通量限制器抑制數(shù)值振蕩；Harten[2]提出總殘差減少(total variable diminishing，TVD)格式求解雙曲線偏微分方程，通過引入r-ψ標(biāo)準(zhǔn)構(gòu)造通量限制器以確保格式的單調(diào)收斂性；Leonard[3]提出的歸一化變量表(normalized variable diagram，NVD)格式建立了界面f和迎風(fēng)節(jié)點U之間的聯(lián)系，確保了格式的單調(diào)性和有界性.這幾類格式所闡述的單調(diào)性在概念上相互聯(lián)系，但具體實現(xiàn)起來差別很大，尤其在非結(jié)構(gòu)網(wǎng)格中實現(xiàn)起來更加復(fù)雜.相比于結(jié)構(gòu)網(wǎng)格，非結(jié)構(gòu)網(wǎng)格在構(gòu)造高精度格式方面存在諸多困難，主要體現(xiàn)在以下兩個方面：

1) 上游迎風(fēng)網(wǎng)格單元重構(gòu) 高精度格式的構(gòu)造一般與上游節(jié)點有關(guān)，對于非結(jié)構(gòu)網(wǎng)格，控制體界面f的上游第二個節(jié)點U的位置和變量值均是未知的，因此，需要對上游迎風(fēng)網(wǎng)格單元進行重構(gòu).目前遠端重構(gòu)方法按插值方式可分為外插和內(nèi)插兩種算法，外插算法主要有Darwish′s算法、Bruner′s算法等；內(nèi)插算法主要有Hou′s算法、Li′s算法和FIFSAM算法等[4-5].

2) 單調(diào)收斂性 高階格式在處理間斷解問題時，例如激波、不連續(xù)界面等問題，間斷處極易出現(xiàn)數(shù)值振蕩甚至產(chǎn)生非物理解，隨著流動變化會繼續(xù)影響流場內(nèi)其他位置的精度甚至導(dǎo)致數(shù)值發(fā)散，如何確保格式的單調(diào)收斂性也一直是非結(jié)構(gòu)網(wǎng)格高精度格式的難點所在.目前抑制數(shù)值振蕩主要方法是通過引入通量限制器對界面插值多項式進行修正，本質(zhì)上是限制數(shù)值解的梯度以避免出現(xiàn)數(shù)值振蕩.

針對非結(jié)構(gòu)網(wǎng)格高精度格式中所面臨的問題，對于上游迎風(fēng)單元的重構(gòu)，文中采用一種線性插值方法重構(gòu)上游迎風(fēng)網(wǎng)格單元的空間位置，并根據(jù)頂點-單元中心間的距離加權(quán)計算上游迎風(fēng)單元處的變量值.相較于上述幾種遠端重構(gòu)算法這種方法同樣可以確保格式的高精度，而且處理方式更簡單有效；同時為了確保計算的單調(diào)收斂性，采用基于施主-受主通量近似思想的CICSAM(compressive interface capturing scheme for arbitrary meshes)格式離散對流方程.

1 數(shù)值方法

1.1 對流方程離散

1.1.1空間離散

對含間斷解或高梯度問題的研究，采用如下二維對流控制方程的守恒形式

(1)

式中：u，v分別為x，y方向上的輸運速度；φ為被輸運量.

采用有限體積法離散式(1)，取任意網(wǎng)格控制體CV(control volume)，控制體的體積為V，將上式寫成如下體積分形式

(2)

由高斯定理，體積分可以轉(zhuǎn)化為面積分形式

(3)

式中：Q為任意標(biāo)量場；V為控制體的體積；S為控制體的封閉單元界面面積；n為方向向量，根據(jù)高斯定理式(2)可以寫成面積分形式

(4)

式中：nx、ny分別為x,y方向的方向向量.將式(4)積分項移到右邊，即有

(5)

式中：Un為單元面上的法向速度；Nf為控制體界面數(shù)目；Sf為界面面積.

1.1.2時間離散

傳統(tǒng)的CICSAM格式采用具有二階精度的Crank-Nicolson格式離散時間項，這是一種半隱式離散方法.由于隱式方法的穩(wěn)定性更好，對于含間斷解、大梯度突變等對流占優(yōu)問題，為了確保數(shù)值解的單調(diào)收斂性，因此，本文采用穩(wěn)定性更好的全隱式離散方法.

(6)

時間項采用具有二階精度的Crank-Nicolson格式進行離散，由此式(6)可離散為如下形式

(7)

由于Pt+Δt是關(guān)于時間的函數(shù)，因此，可以對Pt+Δt進行泰勒級數(shù)展開，其在t時刻的泰勒級數(shù)一階展開可以寫成如下形式

(8)

將式(8)代入式(7)，可得

(9)

2) 全隱式離散 直接對下一時刻的對流通量Pt+Δt，以泰勒級數(shù)對時間t進行展開，有

(10)

同時，將式(10)代入式(6)，可推導(dǎo)為

(11)

式(9)～(10)中括號內(nèi)均為常數(shù)項，等號右邊項Pn中含有未知量φf，因此求解該線性方程組的關(guān)鍵在于求解界面值φf.前文提到界面值φf涉及多個控制體單元的節(jié)點插值，見圖1.通常與上游節(jié)點有關(guān)，下面將分別介紹NVD格式、上游迎風(fēng)節(jié)點的重構(gòu)以及非結(jié)構(gòu)網(wǎng)格CICSAM格式的構(gòu)造.

圖1 有限體積法界面多點插值示意圖

1.2 上游迎風(fēng)網(wǎng)格單元的重構(gòu)

對于任意形式的網(wǎng)格，構(gòu)造具有二階精度的對流離散格式都涉及控制體面f周圍的三個網(wǎng)格單元，包括兩個上游網(wǎng)格單元和一個下游網(wǎng)格單元(見圖1)，其中U，D，A分別表示上游迎風(fēng)網(wǎng)格單元、施主網(wǎng)格單元和受主網(wǎng)格單元.以二維三角形網(wǎng)格為例，見圖2a)，基于線性插值的思想本文將與面f相對的施主網(wǎng)格單元頂點(外節(jié)點)作為迎風(fēng)網(wǎng)格單元U，因此只需計算各單元節(jié)點處的變量值φi(i為節(jié)點編號)，見圖2b)，即可得到迎風(fēng)點的值φU.這種處理方法具有兩個優(yōu)勢：①相較于傳統(tǒng)的遠端重構(gòu)算法這種方法同樣可以確保格式的高精度，而且處理方式更簡單有效；②解決了NVD格式無法直接應(yīng)用于非結(jié)構(gòu)網(wǎng)格中的難點.

圖2 上游迎風(fēng)網(wǎng)格單元空間位置和變量值的重構(gòu)

對于任意頂點i，φi計算公式為

(12)

(13)

1.3 NVD格式

在結(jié)構(gòu)網(wǎng)格中，NVD格式是構(gòu)造許多現(xiàn)代高精度格式的基礎(chǔ)，例如CICSAM格式、HRIC格式、STACS格式和SMART格式等；而對于非結(jié)構(gòu)網(wǎng)格，由于上游迎風(fēng)網(wǎng)格單元的空間位置是未知的，因此一般很難直接采用NVD格式，為此Darwish等[6]發(fā)展了一種考慮空間位置變量的NVSF(normalized variable and space formulation)格式，與NVD格式相比，NVSF格式除了需要對標(biāo)量場φ進行歸一化處理(見式(9))，還需對各標(biāo)量場所在的空間位置x進行了歸一化處理.

而本文采用線性插值方法重構(gòu)上游迎風(fēng)點的位置，考慮到三角形網(wǎng)格和四面體網(wǎng)格的固有特性，將與面f相對的施主網(wǎng)格單元頂點(外節(jié)點)作為迎風(fēng)網(wǎng)格單元U，見圖2a)，這種處理方法消除了空間位置變量的影響，因此可以認(rèn)為迎風(fēng)單元的空間位置是已知的，因此以NVD格式為基礎(chǔ)構(gòu)造非結(jié)構(gòu)網(wǎng)格CICSAM格式依然是可行的，NVD格式的基本定義為

(14)

(15)

(16)

各單元中心和面上的歸一化變量分布情況見圖3.

圖3 原始變量和歸一化變量分布示意圖

(17)

式(17)所構(gòu)造的不等式關(guān)系可以表示為圖4a)的有界性區(qū)域.

圖4 CBC對流有界性準(zhǔn)則

1.4 非結(jié)構(gòu)網(wǎng)格CICSAM格式

在求解非定常對流方程時，庫朗數(shù)是影響數(shù)值解精度的一個重要因素.為此Leonard等[8-9]引入庫朗數(shù)CFL構(gòu)建了顯式CBC有界性準(zhǔn)則，見圖4b)，圖中特征線構(gòu)成的陰影區(qū)域為有界區(qū)域，這種限制準(zhǔn)則確保了數(shù)值解的有界性.在此基礎(chǔ)上，CICSAM格式以色散格式Hyper-C和耗散格式ULTIMATE-QUICKEST(UQ)為基本差分格式，通過權(quán)重函數(shù)f(θf)構(gòu)造將兩種格式進行加權(quán)組合式(18).權(quán)重函數(shù)f(θf)是關(guān)于θf的函數(shù)式(19)，θf定義為界面外法向與速度矢量之間的夾角式(20).

(18)

(19)

(20)

式中：dDA為施主網(wǎng)格中心和受主網(wǎng)格中心之間的向量；φD為施主網(wǎng)格中心處的梯度值.

(21)

(22)

(23)

2 數(shù)值算例

為了檢驗本文全隱式CICSAM格式的數(shù)值精度、單調(diào)收斂性以及上游節(jié)點重構(gòu)方法的可行性，選擇兩個經(jīng)典純對流問題進行數(shù)值驗證，分別是臺階輪廓(Step-Profile)流動、半橢圓輪廓(Semi-Ellipse)流動，基本控制方程見式(1)，每個算例均給定對應(yīng)的速度場.計算數(shù)值解的絕對誤差，采用誤差的L2范數(shù)監(jiān)測數(shù)值收斂性能，為

(24)

2.1 計算域和網(wǎng)格劃分

兩個算例的計算域和和網(wǎng)格均相同，見圖5a)，計算域為1×1的正方形區(qū)域，區(qū)域內(nèi)劃分22 464個三角形網(wǎng)格單元.

圖5 網(wǎng)格劃分與計算域邊界條件設(shè)置(臺階流動)

2.2 驗證算例

1) 臺階輪廓流動 臺階輪廓流動在數(shù)學(xué)物理上呈現(xiàn)為含間斷函數(shù)的特性，存在極端邊界梯度突變，隸屬于間斷解、非連續(xù)問題，例如空氣動力學(xué)的激波問題、水動力學(xué)的自由液面問題.由于臺階流動在間斷處變量的梯度發(fā)生突變，因此該算例可以用來檢驗算法描述突變界面的能力.

φ(0,x)=0 0≤x≤1

φ(0,y)=1 0≤y≤1

(25)

出口處的邊界條件設(shè)置為Neumann邊界條件，見圖5b).

圖6a)為x=0.7截面上不同格式數(shù)值解的精度比較，三條曲線依次表示全隱式CICSAM格式、傳統(tǒng)CICSAM格式(半隱式)的數(shù)值解以及臺階輪廓流動的解析解，可以看出本文構(gòu)造的全隱式CICSAM格式精度比傳統(tǒng)CICSAM格式更好。

圖6 全隱式CICSAM格式的精度和 收斂性驗證(臺階輪廓流動)

2) 半橢圓輪廓流動 對于半橢圓輪廓流動，橢圓輪廓整體過渡平滑，但在兩端仍然存在變量的梯度突變，因此該算例可以檢驗算法同時描述尖銳界面和平滑曲面的能力.本算例采用與臺階流動算例相同的恒定速度場、初始條件和離散方案，在入口處的邊界條件設(shè)置為

(26)

出口處的邊界條件設(shè)置為Neumann邊界條件，見圖7.

圖7 半橢圓輪廓流動計算域邊界條件設(shè)置

圖8a)為x=0.7截面上不同格式數(shù)值解的精度比較，圖8b)為的L2范數(shù)殘差收斂曲線，從兩個圖中的結(jié)果可以得到與臺階流動算例相同的結(jié)論.

圖8 全隱式CICSAM格式的精度和收斂性驗證(Semi-Ellipse輪廓流動)

3 結(jié) 束 語

在NVD格式和顯式CBC對流有界性準(zhǔn)則的基礎(chǔ)上，以色散格式Hyper-C和耗散格式UQ為基本差分格式構(gòu)造了全隱式非結(jié)構(gòu)網(wǎng)格CICSAM格式，采用線性插值的方法重構(gòu)上游迎風(fēng)點的位置，這種處理方法不僅簡化了上游迎風(fēng)節(jié)點的重構(gòu)，同時也實現(xiàn)了NVD格式在非結(jié)構(gòu)網(wǎng)格中的應(yīng)用，這種方法也易于拓展到三維四面體網(wǎng)格.數(shù)值模擬結(jié)果表明，在此基礎(chǔ)上構(gòu)造的CICSAM格式實現(xiàn)了對梯度突變的間斷界面精確捕捉，相比傳統(tǒng)CICSAM格式其具有更好的數(shù)值收斂性能.