基于二維非結(jié)構(gòu)網(wǎng)格的全隱式CICSAM格式*

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

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

0 引 言

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

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

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

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

1 數(shù)值方法

1.1 對(duì)流方程離散

1.1.1空間離散

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

(1)

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

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

(2)

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

(3)

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

(4)

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

(5)

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

1.1.2時(shí)間離散

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

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

(11)

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

圖1 有限體積法界面多點(diǎn)插值示意圖

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

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

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

對(duì)于任意頂點(diǎn)i，φi計(jì)算公式為

(12)

(13)

1.3 NVD格式

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

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

(14)

(15)

(16)

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

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

(17)

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

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

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

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

(18)

(19)

(20)

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

(21)

(22)

(23)

2 數(shù)值算例

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

(24)

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

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

圖5 網(wǎng)格劃分與計(jì)算域邊界條件設(shè)置(臺(tái)階流動(dòng))

2.2 驗(yàn)證算例

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

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

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

(25)

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

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

圖6 全隱式CICSAM格式的精度和 收斂性驗(yàn)證(臺(tái)階輪廓流動(dòng))

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

(26)

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

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

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

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

3 結(jié) 束 語(yǔ)

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