高宜勝, 伍貽兆, 夏 健, 徐兆可

(南京航空航天大學(xué) 航空宇航學(xué)院, 江蘇 南京 210016)

基于物面弱邊界條件的并行全隱式流場(chǎng)求解方法

高宜勝, 伍貽兆*, 夏 健, 徐兆可

(南京航空航天大學(xué) 航空宇航學(xué)院, 江蘇 南京 210016)

提出了一種基于物面弱邊界條件的Navier-Stokes方程并行全隱式求解方法。首先，實(shí)現(xiàn)了非結(jié)構(gòu)網(wǎng)格格點(diǎn)格式有限體積法的粘性物面弱邊界條件。以此為基礎(chǔ)，推導(dǎo)了相應(yīng)的粘性Jacobian矩陣。相比傳統(tǒng)的物面強(qiáng)邊界條件，弱邊界條件的應(yīng)用簡(jiǎn)化了該推導(dǎo)過(guò)程。接著，為了解決Spalart-Allmaras (S-A)模型的收斂和保正性問(wèn)題，提出了基于物面弱邊界條件的S-A模型無(wú)條件保正性收斂方法，通過(guò)構(gòu)造特別的隱式矩陣同時(shí)保證S-A模型的收斂和湍流工作變量的正值。流場(chǎng)方程全隱式時(shí)間推進(jìn)和無(wú)條件保正性收斂方法所得到的線性方程組采用多色高斯-塞德?tīng)柕ㄟM(jìn)行求解，并由MPI實(shí)現(xiàn)并行化。通過(guò)若干二維和三維算例對(duì)本文所提出的方法進(jìn)行了測(cè)試，其結(jié)果表明：1) 隨著罰系數(shù)的增大，基于弱邊界條件的結(jié)果趨向于基于強(qiáng)邊界條件的結(jié)果；2) 對(duì)于力的計(jì)算，弱邊界條件自動(dòng)滿足通量一致關(guān)系而無(wú)需特別處理。

弱邊界條件；Navier-Stokes方程；全隱式方法；保正性方法； 湍流模型

0 引 言

邊界條件處理是CFD計(jì)算中極其重要的一個(gè)環(huán)節(jié)，直接關(guān)系到計(jì)算能否收斂、收斂速度快慢等方面。通常，邊界條件的處理可以分為強(qiáng)邊界條件和弱邊界條件。強(qiáng)邊界條件直接給定邊界上的數(shù)值，而弱邊界條件則是通過(guò)通量的方式間接引入邊界貢獻(xiàn)。在當(dāng)前的Navier-Stokes方程數(shù)值計(jì)算中，物面邊界為無(wú)滑移邊界，一般采用強(qiáng)邊界條件處理，即直接指定物面的速度。強(qiáng)邊界條件概念直觀，實(shí)現(xiàn)方便，是目前應(yīng)用最多的Navier-Stokes方程物面邊界條件處理方式。

近年來(lái)，物面弱邊界條件也逐漸得到了研究和應(yīng)用[1-4]。這是因?yàn)橄啾葌鹘y(tǒng)的強(qiáng)邊界條件，弱邊界條件具有以下優(yōu)勢(shì)：(1)在有限差分方法中，弱邊界條件配合分部求和算子(summation-by-parts operators)，能夠?qū)崿F(xiàn)數(shù)學(xué)上可證明的穩(wěn)定格式[5]，而基于強(qiáng)邊界條件尚無(wú)相應(yīng)的結(jié)果；(2)Eliasson等人利用非結(jié)構(gòu)網(wǎng)格求解器Edge研究了弱邊界條件和強(qiáng)邊界條件對(duì)于定常問(wèn)題收斂性的影響，發(fā)現(xiàn)弱邊界條件能夠改善收斂性，并提高多重網(wǎng)格的收斂速度[6]；(3)在實(shí)現(xiàn)流場(chǎng)全隱式求解或者推導(dǎo)離散伴隨方程的時(shí)候，需要推導(dǎo)Jacobian矩陣(或者其轉(zhuǎn)置形式)，由于強(qiáng)邊界條件下物面點(diǎn)的殘值和流場(chǎng)變量無(wú)關(guān)，因此Jacobian矩陣(或者其轉(zhuǎn)置矩陣)需要特別處理[7]；而在弱邊界條件下，可以直接推導(dǎo)相應(yīng)的Jacobian矩陣，簡(jiǎn)化了物面邊界的處理；(4)最近， Stück發(fā)現(xiàn)在強(qiáng)邊界條件下，傳統(tǒng)的力計(jì)算公式不滿足通量一致關(guān)系(flux consistency)，因此強(qiáng)邊界條件下計(jì)算力函數(shù)需要新的計(jì)算公式[8]。正因?yàn)榫哂猩鲜鰞?yōu)點(diǎn)，弱邊界條件正逐步得到應(yīng)用[9]。而目前國(guó)內(nèi)關(guān)于物面弱邊界條件的研究尚屬空白。

本文提出了一種基于物面弱邊界條件的非結(jié)構(gòu)網(wǎng)格并行全隱式流場(chǎng)求解方法。首先，實(shí)現(xiàn)了非結(jié)構(gòu)網(wǎng)格格點(diǎn)格式有限體積法的物面弱邊界條件并推導(dǎo)了相應(yīng)的Jacobian矩陣。與強(qiáng)邊界條件相比，Jacobian矩陣可以直接從黏性通量表達(dá)式推導(dǎo)得到而無(wú)需特別處理。為了穩(wěn)定求解湍流問(wèn)題，改進(jìn)了Shende等人的Spalart-Allmaras(S-A)模型無(wú)條件保正性收斂方法(unconditionally positive-convergent implicit scheme， UPC)[10]，提出了弱邊界條件下的UPC方法，保證弱邊界條件下S-A模型的收斂性和湍流工作變量的正值。采用多色高斯-塞德?tīng)柕?multicolor Gauss-Seidel iteration，MCGS)求解由流場(chǎng)方程全隱式時(shí)間推進(jìn)和UPC方法所建立的線性方程組，并利用MPI實(shí)現(xiàn)并行計(jì)算。通過(guò)若干算例研究了弱邊界條件的計(jì)算結(jié)果，并通過(guò)與強(qiáng)邊界條件結(jié)果的對(duì)比，說(shuō)明了弱邊界條件所具有的優(yōu)點(diǎn)。

1 物面弱邊界條件的實(shí)現(xiàn)和Jacobian矩陣推導(dǎo)

1.1 物面弱邊界條件的具體實(shí)現(xiàn)

首先采用Eliasson等的方法[6]構(gòu)造物面弱邊界條件。對(duì)于Navier-Stokes方程的絕熱物面邊界，其對(duì)流通量的處理和Euler方程滑移物面邊界的處理一致，因此只考慮其黏性通量。物面邊界的黏性通量可以表示為:

(1)

[11]提供了多種計(jì)算式(1)中物面速度法向梯度?u/?n的表達(dá)式，這里選用其中的第一種，其表示為:

(2)

1.2 物面弱邊界條件的Jacobian矩陣推導(dǎo)

為了實(shí)現(xiàn)流場(chǎng)全隱式求解，需要推導(dǎo)相應(yīng)的Jacobian矩陣。在強(qiáng)邊界條件下，由于物面邊界點(diǎn)的速度直接設(shè)置為零，物面點(diǎn)的動(dòng)量方程不起作用，物面點(diǎn)的殘值和流場(chǎng)變量沒(méi)有直接聯(lián)系，因此需要采用類似有限元中Dirichlet 邊界條件的處理方式，將全局Jacobian矩陣中物面點(diǎn)的動(dòng)量方程所對(duì)應(yīng)的行消去，再將對(duì)應(yīng)的對(duì)角項(xiàng)設(shè)為1，并將對(duì)應(yīng)的右端項(xiàng)設(shè)為零。而在弱邊界條件下，物面無(wú)滑移邊界條件是利用式(1)計(jì)算黏性通量的方式引入的，因此其Jacobian矩陣可以直接通過(guò)式(1)的求導(dǎo)得到。這里給出二維情況下的推導(dǎo)，并只考慮對(duì)角項(xiàng)即物面點(diǎn)對(duì)自身黏性通量的貢獻(xiàn)。

設(shè)W為守恒變量，直接推導(dǎo)黏性通量對(duì)守恒變量的Jacobian矩陣?G0/?W較為復(fù)雜，所以先考慮黏性通量對(duì)原始變量的Jacobian矩陣?G0/?Q，其中Q=[ρuvp]T表示原始變量。設(shè)?u/?n=[?u/?n?v/?n]T,對(duì)式(2)求導(dǎo)，可以得到:

(3)

對(duì)式(1)求導(dǎo)，并代入式(3)的結(jié)果，有:

式中，nx，ny為n0的分量。根據(jù)式(4)可以得到:

2 基于弱邊界條件的Spalart-Allmaras模型無(wú)條件保正性收斂方法

湍流模型的求解是CFD計(jì)算中一個(gè)較為復(fù)雜的問(wèn)題。因?yàn)榫W(wǎng)格質(zhì)量不好、數(shù)值離散誤差等原因會(huì)造成計(jì)算過(guò)程中S-A模型的工作變量變?yōu)樨?fù)值，使得計(jì)算失敗。為了避免這個(gè)問(wèn)題，可以直接用max函數(shù)進(jìn)行截?cái)啵沁@樣往往會(huì)造成不連續(xù)而影響收斂。針對(duì)這個(gè)問(wèn)題，Mor-Yossef等人在兩方程模型上提出了無(wú)條件保正性收斂方法(UPC)，通過(guò)構(gòu)造特別的迭代矩陣使得每一步的工作變量都是非負(fù)值，并且能夠保證湍流模型計(jì)算的收斂[12-13]。隨后，Shende等人又將其發(fā)展到S-A模型上，實(shí)現(xiàn)了強(qiáng)邊界條件下S-A模型的穩(wěn)定計(jì)算[11]。本文改進(jìn)了該方法，提出了弱邊界條件下的UPC方法，用于湍流的穩(wěn)定求解。

(6)

對(duì)于每個(gè)控制體i，空間離散后的半離散形式S-A模型可以寫(xiě)為:

式中，Ri表示控制體i的對(duì)流項(xiàng)、擴(kuò)散項(xiàng)和反擴(kuò)散項(xiàng)綜合起來(lái)的殘值，Ωi表示控制體i的體積，Si表示由生成項(xiàng)和消耗項(xiàng)組成的源項(xiàng)。采用一階后向Euler差分離散時(shí)間項(xiàng)，則得到全局線性方程組:

(a)M是M矩陣[15]；

(b) -Rn+ΩSn+Mνn是非負(fù)矢量，即每個(gè)分量都是非負(fù)值。

參考文獻(xiàn)[13]證明了如果矩陣M滿足上述2個(gè)條件，則全隱式時(shí)間推進(jìn)

(9)

具有無(wú)條件收斂性和保正性，即ν始終為非負(fù)矢量。

基于這個(gè)結(jié)果，需要構(gòu)造M矩陣并且使其同時(shí)滿足條件(b)。參考文獻(xiàn)[15]對(duì)于構(gòu)造M矩陣提供了一個(gè)充分條件：如果一個(gè)矩陣滿足下面3個(gè)要求，則該矩陣為M矩陣：(1)矩陣的所有對(duì)角線元素大于零；(2)矩陣的所有非對(duì)角線元素非正；(3)對(duì)角占優(yōu)，即每行的對(duì)角線元素大于該行所有非對(duì)角元素的絕對(duì)值之和。根據(jù)這個(gè)充分條件，下面分別考慮式(6)中對(duì)流項(xiàng)、擴(kuò)散項(xiàng)、反擴(kuò)散項(xiàng)、源項(xiàng)的近似Jacobian矩陣構(gòu)造。

令Rinv表示對(duì)流項(xiàng)的殘值，由于對(duì)流項(xiàng)采用一階逆風(fēng)格式進(jìn)行離散，那么得到

(10)

(11)

u⊥i,j=ui,jnx+vi,jny+wi,jnz

(12)

(13)

這里，J表示控制體i和控制體j的交界面，SJ表示交界面面積，ui,j，vi,j，wi,j表示控制體i/j的速度分量，nx、ny、nz為單位法向矢量。根據(jù)式(10)～(13)，可以得到:

(14)

(15)

即對(duì)角線元素大于零，非對(duì)角線元素非正，但無(wú)法滿足對(duì)角占優(yōu)的要求。為了使得對(duì)角占優(yōu)，令:

(16)

(17)

則式(14)可以寫(xiě)為:

(18)

那么

(19)

(20)

因此對(duì)流項(xiàng)的近似Jacobian矩陣滿足上述的充分條件的要求。

(21)

(22)

(23)

這里，lij表示邊ij的長(zhǎng)度，rij表示邊ij的單位邊矢量，nJ控制體界面J的單位法向矢量。因此有:

Rvis=

(24)

(25)式(25)表示擴(kuò)散項(xiàng)的近似Jacobian矩陣?？紤]到整體Jacobian矩陣中還包括對(duì)流項(xiàng)和反擴(kuò)散項(xiàng)的矩陣構(gòu)造，所以擴(kuò)散項(xiàng)的矩陣構(gòu)造只要求對(duì)角項(xiàng)大于零而非對(duì)角項(xiàng)小于零，而不要求對(duì)角占優(yōu)。

對(duì)于反擴(kuò)散項(xiàng)，其離散形式為:

(26)

(27)

采用類似耗散項(xiàng)的處理，有

(28)

考慮到

(29)

令

(30)

這里反耗散項(xiàng)的矩陣構(gòu)造本身也不滿足M矩陣的充分條件，但是綜合對(duì)流項(xiàng)、耗散項(xiàng)和反耗散項(xiàng)的矩陣構(gòu)造以后仍然滿足M矩陣的充分條件。

接著考慮邊界點(diǎn)對(duì)近似Jacobian矩陣的貢獻(xiàn)。對(duì)于遠(yuǎn)場(chǎng)邊界，由于遠(yuǎn)場(chǎng)工作變量梯度較小，所以其耗散項(xiàng)和反耗散項(xiàng)較小，忽略它們的貢獻(xiàn)，只考慮遠(yuǎn)場(chǎng)邊界對(duì)流項(xiàng)的貢獻(xiàn)。遠(yuǎn)場(chǎng)邊界的對(duì)流項(xiàng)殘值可表示為

Rfar=F⊥J,farSfar

(31)

(32)

u⊥,far=ufarnx+vfarny+wfarnz

(33)

(34)

這里，ufar，vfar，wfar為遠(yuǎn)場(chǎng)邊界的速度，由流場(chǎng)的特征邊界條件得到，Sfar為遠(yuǎn)場(chǎng)邊界面積。相應(yīng)的矩陣構(gòu)造為:

(35)

對(duì)于對(duì)稱邊界，由于對(duì)稱邊界的法向速度為零，并且工作變量的法向梯度為零，因此沒(méi)有通量貢獻(xiàn)，也沒(méi)有Jacobian矩陣貢獻(xiàn)。

對(duì)于物面邊界，由于物面上工作變量為零，因此沒(méi)有對(duì)流項(xiàng)和反擴(kuò)散項(xiàng)的貢獻(xiàn)，只有擴(kuò)散項(xiàng)的貢獻(xiàn)。弱邊界條件可以按照類似式(2)的方式計(jì)算工作變量的法向梯度:

(36)

這里，φ為S-A模型弱邊界條件的罰系數(shù)。那么物面邊界的擴(kuò)散項(xiàng)對(duì)殘值的貢獻(xiàn)為:

(37)

式中，Swall為物面邊界面積。相應(yīng)的矩陣構(gòu)造為:

對(duì)于源項(xiàng)，將其表示為:

(39)

對(duì)應(yīng)的矩陣構(gòu)造為:

(40)

這里采用類似參考文獻(xiàn)[14]中源項(xiàng)的處理方法，用max函數(shù)限制源項(xiàng)Jacobian矩陣為非負(fù)值。

綜合對(duì)流項(xiàng)、擴(kuò)散項(xiàng)、反擴(kuò)散項(xiàng)和源項(xiàng)的矩陣構(gòu)造，得到:

=-Ri+ΩiSi

(41)

式(41)中的近似矩陣構(gòu)造滿足了UPC的條件(a)，但由于實(shí)際求解的是非線性方程組，每個(gè)偽時(shí)間步迭代中式(41)的右端項(xiàng)都會(huì)改變，所以條件(b)不一定能滿足。令

-Ri+ΩiSi

(42)

如果Φ<0，令:

(43)系數(shù)1.5考慮了計(jì)算機(jī)采用浮點(diǎn)數(shù)表示實(shí)數(shù)而帶來(lái)數(shù)值精度影響。將式(43)加入式(41)的對(duì)角項(xiàng)，得到:

=-Ri+ΩiSi

(44)

需要說(shuō)明的是，由于UPC方法采用了特別的矩陣構(gòu)造，為了保證UPC特性，在式(44)的求解過(guò)程中必須保證矩陣不發(fā)生改變，所以不能使用LU-SGS這類引入矩陣分解誤差的求解方法。為此，本文采用多色高斯-塞德?tīng)柕ㄇ蠼饩€性方程組，這是因?yàn)楦咚?塞德?tīng)柕芯仃嚤３植蛔?，因此保證了UPC特性。其具體過(guò)程在下節(jié)詳細(xì)敘述。

3 多色高斯-塞德?tīng)柕?/h2>

流場(chǎng)的全隱式時(shí)間推進(jìn)和上一節(jié)介紹的UPC方法都需要在每個(gè)偽時(shí)間步中求解大型稀疏線性方程組。同時(shí)考慮到上文所述的UPC求解不能引入矩陣分解誤差的要求，因此采用多色高斯-塞德?tīng)柕?MCGS)[17]求解線性方程組，并利用MPI實(shí)現(xiàn)了該方法的分布式并行計(jì)算。

待求解的線性方程組可以表示為:

Ax=b

(45)

這里，A表示(大型稀疏)矩陣，x為解向量，b為已知右端列向量。傳統(tǒng)的高斯-塞德?tīng)柕梢詫?xiě)為:

(i=1,2,…,n)

(46)

式中，i表示點(diǎn)的編號(hào)，對(duì)應(yīng)矩陣的第i行，上標(biāo)m表示迭代計(jì)數(shù)，aii為矩陣第i行的對(duì)角項(xiàng)，bi為第i行對(duì)應(yīng)的右端項(xiàng)，Adj(i)表示點(diǎn)i的相鄰點(diǎn)，aij為矩陣第i行的非對(duì)角項(xiàng)，n為點(diǎn)數(shù)，即解向量的維數(shù)。為了實(shí)現(xiàn)并行求解，首先對(duì)網(wǎng)格進(jìn)行染色。非結(jié)構(gòu)網(wǎng)格點(diǎn)周圍存在數(shù)量不一的相鄰點(diǎn)，需要多種顏色進(jìn)行染色，來(lái)保證每個(gè)點(diǎn)和周圍相鄰點(diǎn)的顏色不同，如圖2所示的網(wǎng)格就需要6種顏色進(jìn)行染色。本文采用最常用的貪心染色算法[18]，其流程如圖3所示。數(shù)值測(cè)試顯示，非結(jié)構(gòu)網(wǎng)格染色所需的顏色數(shù)量大致在5～13之間。染色完成后，染成同一種顏色的點(diǎn)相互之間不相鄰，沒(méi)有數(shù)據(jù)依賴，因此可以按照顏色順序進(jìn)行并行更新。整個(gè)并行MCGS迭代法的計(jì)算過(guò)程如圖4所示。圖中在Loop3循環(huán)完成以后，采用MPI非阻塞通信更新分區(qū)邊界上相應(yīng)顏色的ghost點(diǎn)的值，就實(shí)現(xiàn)了并行求解。

4 計(jì)算結(jié)果

4.1 NACA0012翼型的層流流動(dòng)

首先考慮二維層流流動(dòng)的結(jié)果。計(jì)算網(wǎng)格如圖5所示，物面附近采用四邊形單元而其它部分采用三角形單元，網(wǎng)格單元數(shù)為13156，網(wǎng)格點(diǎn)數(shù)為8741。計(jì)算狀態(tài)為馬赫數(shù)Ma=0.5，迎角α=1°，Reynolds數(shù)Re=5000。對(duì)流通量采用JST (Jameson-Schmidt-Turkel)格式計(jì)算[19]。而在黏性通量的計(jì)算中，先由Green-Gauss公式得到原始變量的梯度，再采用基于邊修正的方法計(jì)算控制體邊界的梯度[20]。前10步CFL數(shù)從10線性增長(zhǎng)至10000，然后維持為10000。MCGS迭代次數(shù)為15次。計(jì)算平臺(tái)為MacBook Pro 15(2013 Later)，編譯器為GFortran 4.9，采用串行方式運(yùn)行。

圖6給出了弱邊界條件下物面罰系數(shù)σ分別取1、10、100時(shí)的翼型表面壓強(qiáng)系數(shù)Cp的分布情況，同時(shí)為了與強(qiáng)邊界條件的結(jié)果進(jìn)行對(duì)比，還給出了開(kāi)源CFD代碼SU2[21]在相同狀態(tài)下采用強(qiáng)邊界條件計(jì)算的結(jié)果(需要說(shuō)明的是，SU2的對(duì)流通量計(jì)算同樣采用JST格式，但黏性通量計(jì)算方法和邊界條件處理與本文方法不同)?？梢钥吹?條曲線在除了后緣的位置以外幾乎完全重合。圖7給出了后緣部分的放大圖，可以看到隨著罰系數(shù)的增大，弱邊界條件的結(jié)果將趨向于強(qiáng)邊界條件的結(jié)果。

圖8給出了物面罰系數(shù)σ分別取1、10、100時(shí)以及SU2計(jì)算的翼型表面摩擦力系數(shù)Cf的分布情況，可以看出4條曲線幾乎完全重合，而圖9給出了后緣部分的放大圖，4條曲線仍然重合得很好。表明采用弱邊界條件可以得到和強(qiáng)邊界條件一致的結(jié)果。

圖10顯示了物面罰系數(shù)σ分別取1、10、100時(shí)翼型上表面一點(diǎn)以及附近點(diǎn)的速度矢量和速度大小。從圖中可以看出，由于使用了弱邊界條件，壁面點(diǎn)的速度不再為零，但是其數(shù)值大小比其法向點(diǎn)小的多，而且隨著罰系數(shù)的增大將更加趨近于零，因此其對(duì)流場(chǎng)計(jì)算的影響可以忽略，表明了弱邊界條件可用于黏性流動(dòng)計(jì)算。

圖11給出了物面罰系數(shù)σ分別取1、10、100時(shí)密度的收斂歷史，可以看到不同的罰系數(shù)對(duì)于收斂影響很小，為了更接近實(shí)際的無(wú)滑移條件，計(jì)算中可以采用較大的罰系數(shù)。

4.2 RAE2822翼型的湍流流動(dòng)

本算例說(shuō)明本文方法用于二維湍流計(jì)算的情況。計(jì)算網(wǎng)格如圖 12所示，物面附近采用四邊形單元捕捉湍流邊界層，而其余部分采用三角形單元，網(wǎng)格單元數(shù)為22842，網(wǎng)格點(diǎn)數(shù)為13937。計(jì)算狀態(tài)為馬赫數(shù)Ma=0.729，迎角α=2.31°，Reynolds數(shù)Re=6.5×106。對(duì)流通量和黏性通量計(jì)算方法同上。流場(chǎng)方程和湍流模型方程的CFL數(shù)在前10步從10線性增長(zhǎng)至1000，然后維持為1000。MCGS迭代次數(shù)均為15次。物面罰系數(shù)均取100。計(jì)算平臺(tái)和上例一樣，同樣采用串行計(jì)算。

圖13給出了本文采用弱邊界條件和SU2采用強(qiáng)邊界條件計(jì)算的表面壓強(qiáng)系數(shù)以及和實(shí)驗(yàn)數(shù)據(jù)[22]的對(duì)比?？梢钥吹剑疚牡挠?jì)算結(jié)果和SU2的結(jié)果幾乎重合，與實(shí)驗(yàn)數(shù)據(jù)也吻合較好，表明本文基于弱邊界條件的求解體系能夠有效地計(jì)算湍流流動(dòng)，得到與傳統(tǒng)強(qiáng)邊界條件一致的結(jié)果。

下面利用該算例來(lái)說(shuō)明不同邊界條件處理對(duì)力計(jì)算的影響。對(duì)于力的計(jì)算，可以采用基于物面或者基于遠(yuǎn)場(chǎng)的方法，如圖14所示，左圖表示對(duì)整個(gè)物面積分計(jì)算力，右圖表示對(duì)整個(gè)遠(yuǎn)場(chǎng)積分計(jì)算力。其標(biāo)準(zhǔn)計(jì)算公式可以表示為[8]：

這里dk為單位矢量，對(duì)于阻力取來(lái)流方向，對(duì)于升力取與來(lái)流垂直方向，Γ0根據(jù)不同的力計(jì)算方法表示遠(yuǎn)場(chǎng)邊界或者壁面邊界，可見(jiàn)圖14，δkl為克羅內(nèi)克符號(hào)，nl為單位外法向矢量，ΔΓ為表面單元面積。由于流場(chǎng)中沒(méi)有額外的力源，理論上基于物面或者基于遠(yuǎn)場(chǎng)這兩種方法計(jì)算得到的力應(yīng)該是完全一樣的。但最近Stuck[8]發(fā)現(xiàn)，在強(qiáng)邊界條件下，兩種方法使用式(47)計(jì)算出來(lái)的力是不一樣的，存在著通量不一致(flux inconsistency)的問(wèn)題。這是因?yàn)閺?qiáng)邊界條件直接給定物面點(diǎn)的速度，動(dòng)量方程不起作用，使得物面處沒(méi)有嚴(yán)格考慮通量項(xiàng)的貢獻(xiàn)，而在遠(yuǎn)場(chǎng)邊界，特征邊界條件是一種弱邊界條件處理，是通過(guò)通量的方式引入的，因此遠(yuǎn)場(chǎng)處的所有通量項(xiàng)對(duì)于力的計(jì)算都是有貢獻(xiàn)的，這樣造成了整個(gè)計(jì)算區(qū)域的通量貢獻(xiàn)不守恒。而采用了弱邊界條件以后，物面邊界條件也是以通量的方式引入的，同時(shí)整個(gè)流場(chǎng)采用守恒形式有限體積法離散，因此整個(gè)流場(chǎng)區(qū)域保持通量守恒，自動(dòng)滿足了通量一致關(guān)系。表1給出了強(qiáng)邊界條件和弱邊界條件下由基于物面和基于遠(yuǎn)場(chǎng)方法所得到的升力系數(shù)和阻力系數(shù)的值?？梢钥吹?，在強(qiáng)邊界條件下，基于物面和基于遠(yuǎn)場(chǎng)的方法計(jì)算得到的升力系數(shù)和阻力系數(shù)是不同的，表現(xiàn)出了通量不一致的問(wèn)題，這與Stuck的結(jié)論一致。為了使得強(qiáng)邊界條件下兩種方法計(jì)算的力相同，必須采用參考文獻(xiàn)[8]提出的公式而不能使用式(47)。而在弱邊界條件下，基于物面和基于遠(yuǎn)場(chǎng)方法計(jì)算得到的升力系數(shù)和阻力系數(shù)在機(jī)器精度的范圍內(nèi)一致，自動(dòng)滿足通量一致關(guān)系，無(wú)需任何特別處理，表明弱邊界條件相比強(qiáng)邊界條件在計(jì)算力的時(shí)候更加方便可靠。

表1 RAE2822翼型湍流流動(dòng)強(qiáng)邊界條件和弱邊界條件下計(jì)算得到的力Table 1 Force functionals obtained by strong and weak imposedboundary conditions for RAE2822 airfoil turbulent flow

4.3 ONERA M6機(jī)翼的湍流流動(dòng)

接下來(lái)研究本文方法應(yīng)用于三維湍流流場(chǎng)求解的情況。計(jì)算采用非結(jié)構(gòu)混合網(wǎng)格，如圖15所示，邊界層內(nèi)采用三棱柱單元，遠(yuǎn)場(chǎng)采用四面體單元，過(guò)渡部分采用金字塔單元。網(wǎng)格單元數(shù)為915923，網(wǎng)格點(diǎn)數(shù)為2323893。計(jì)算狀態(tài)為馬赫數(shù)Ma=0.8395，迎角α=3.06°，Reynolds數(shù)Re=11.72×106。對(duì)流通量和黏性通量計(jì)算方法同上。流場(chǎng)方程和湍流模型方程的CFL數(shù)在前10步從10線性增長(zhǎng)至500，然后維持為500。MCGS迭代次數(shù)均為15次。物面罰系數(shù)均取10000。計(jì)算環(huán)境為天津超級(jí)計(jì)算中心天河1A的2個(gè)計(jì)算節(jié)點(diǎn)。每個(gè)節(jié)點(diǎn)擁有2個(gè)Intel Xeon X5670處理器，每個(gè)處理器擁有6個(gè)核心，主頻為2.93GHz，操作系統(tǒng)為RHEL，編譯器為Intel Fortran 11.1，MPI環(huán)境為天河自主實(shí)現(xiàn)的MPI。計(jì)算采用24核并行計(jì)算。

圖16給出了本文和SU2計(jì)算的機(jī)翼展向44%和65%這兩個(gè)位置的壓強(qiáng)系數(shù)分布以及同實(shí)驗(yàn)數(shù)據(jù)[23]的對(duì)比。本文的計(jì)算結(jié)果與SU2的結(jié)果非常接近，因此在各展向位置的壓強(qiáng)系數(shù)分布曲線幾乎完全重合，并且和實(shí)驗(yàn)數(shù)據(jù)也吻合得很好，說(shuō)明本文的弱邊界條件適用于三維湍流流動(dòng)的計(jì)算。

表2給出了該算例中采用強(qiáng)邊界條件和弱邊界條件時(shí)由基于物面和基于遠(yuǎn)場(chǎng)方法所得到的升力系數(shù)和阻力系數(shù)的值?？梢钥吹剑投S情況類似，強(qiáng)邊界條件下同樣出現(xiàn)了通量不一致的問(wèn)題，而弱邊界條件計(jì)算得到的力系數(shù)一致，自動(dòng)滿足通量一致關(guān)系。

表2 ONERA M6機(jī)翼湍流流動(dòng)強(qiáng)邊界條件和弱邊界條件下計(jì)算得到的力Table 2 Force functionals obtained by strong and weak imposed boundary conditions for ONERA M6 wing turbulent flow

圖 17對(duì)本文方法的并行加速比和理想加速比進(jìn)行了對(duì)比?？梢钥闯霎?dāng)CPU核心數(shù)較少時(shí)，并行效率較高，而當(dāng)CPU核心數(shù)達(dá)到24核時(shí)并行效率約為77%。因?yàn)楸疚牟捎玫腗CGS迭代，在每次循環(huán)一種顏色時(shí)都需要通信一次，通信量較大，如果需要進(jìn)一步改進(jìn)并行效率，后續(xù)將考慮通信和計(jì)算重疊的方法或者采用通信量更小的迭代方法。

5 結(jié) 論

本文基于近年來(lái)逐步發(fā)展起來(lái)的物面弱邊界條件，建立了并行全隱式流場(chǎng)求解方法。在非結(jié)構(gòu)網(wǎng)格格點(diǎn)格式有限體積法中實(shí)現(xiàn)了物面弱邊界條件，并推導(dǎo)了相應(yīng)的Jacobian矩陣。其Jacobian矩陣可以直接由黏性通量表達(dá)式推導(dǎo)得到，簡(jiǎn)化了流場(chǎng)全隱式計(jì)算的推導(dǎo)處理。在此基礎(chǔ)上，提出了弱邊界條件下S-A模型的無(wú)條件保正性收斂方法，能夠保證弱邊界條件下S-A模型求解的收斂性和湍流工作變量的正值，實(shí)現(xiàn)了湍流流動(dòng)的穩(wěn)定求解。采用多色高斯-塞德?tīng)柕ㄇ蠼馊[式時(shí)間推進(jìn)和S-A模型無(wú)條件保正性收斂方法所得到的線性方程組，并利用MPI實(shí)現(xiàn)并行計(jì)算。對(duì)二維、三維條件下的若干層流和湍流問(wèn)題進(jìn)行了計(jì)算分析，結(jié)果表明：(1)隨著物面罰系數(shù)的增大，弱邊界條件下的結(jié)果趨向于強(qiáng)邊界條件的結(jié)果；(2)對(duì)于力的計(jì)算，弱邊界條件保證了整個(gè)流場(chǎng)的通量守恒，自動(dòng)滿足通量一致關(guān)系，無(wú)需特別的處理，因而比傳統(tǒng)的強(qiáng)邊界條件更加方便可靠。

后續(xù)研究工作將深入以下幾個(gè)方面：(1)力計(jì)算的通量一致問(wèn)題與對(duì)偶一致問(wèn)題(dual consistency)[24]密切相關(guān)，近年來(lái)正得到越來(lái)越多的重視，對(duì)于計(jì)算結(jié)果的影響需要進(jìn)一步研究；(2)本文目前的并行效率在CPU數(shù)量增多以后下降較為明顯，需要研究提高并行效率的策略，比如通信和計(jì)算重疊的方法，以便用于更大規(guī)模的問(wèn)題。

參 考 文 獻(xiàn)：

[1]Nordstrom J, Eriksson S, Eliasson P. Weak and strong wall boundary procedures and convergence to steady-state of the Navier-Stokes equations[J]. Journal of Computational Physics, 2012, 231(14):4867-4884.

[2]Bazilevs Y, Michler C, Calo V M, et al. WeakDirichlet boundary conditions for wall-bounded turbulent flows[J]. Computer Methods in Applied Mechanics and Engineering, 2007, 196(49-52):4853-4862.

[3]Bazilevs Y, Hughes T J R. Weak imposition of Dirichlet boundary conditions in fluid mechanics[J]. Computers & Fluids, 2007, 36(1):12-26.

[4]Chandrashekar P. Finite volume discretization of heat equation and compressible Navier-Stokes equations with weak Dirichlet boundary condition on triangular grids[J]. International Journal of Advances in Engineering Sciences and Applied Mathematics, 2016, 8(3):174-193.

[5]Fisher T C, Carpenter M H. High-order entropy stable finite difference schemes for nonlinear conservations laws: Finite domains[J]. Journal of Computational Physics, 2013, 252: 518-557.

[6]Eliasson P, Eriksson S, Nordstrom J. The influence of weak and strong solid wall boundary conditions on the convergence to steady-state of the Navier-Stokes equations. AIAA 2009-3551[R]. Reston: AIAA, 2009.

[7]Giles M B, Duta M C, Muller J D, et al. Algorithm developments for discrete adjoint methods[J]. AIAA Journal, 2003, 41(2):198-205.

[8]Stuck A. Anadjoint view on flux consistency and strong wall boundary conditions to the Navier-Stokes equations[J]. Journal of Computational Physics, 2015, 301: 247-264.

[9]Lakshminarayan V, Sitaraman J, Roget B, et al. Development and validation of a multi-strand solver for complex aerodynamic flows. AIAA 2016-1581[R]. Reston: AIAA, 2016.

[10]Shende N, Mor-Yossef Y. Robust implementation of the Spalart-Allmaras turbulence model for unstructured grid[C]//European Conference on Computational Fluid Dynamics 2010, 2010.

[11]Eliasson P, Nordstrom J. The influence of viscous operator and wall boundary conditions on the accuracy of the Navier-Stokes equations. AIAA 2013-2956[R]. Reston: AIAA, 2013.

[12]Moryossef Y, Levy Y. Unconditionally positive implicit procedure for two-equation turbulence models: Appliciation tok-ωturbulence models[J]. Journal of Computational Physics, 2006, 220(1):88-108.

[13Mor-Yossef Y, Levy Y. The unconditionally positive-convergent implicit time integration scheme for two-equation turbulence models: Revisited[J]. Computers & Fluids, 2007, 38(10):1984-1994.

[14]Spalart P, Allmaras S. A one-equation turbulence model for aerodynamic flows. AIAA 92-0439[R]. Reston: AIAA, 1992.

[15]Berman A, Plemmons R. Nonnegative matrices in the mathematical sciences[M]. New York: Academic Press, 1979.

[16]Svard M, Nordstrom J. Stability of finite volume approximations for the Laplacian operator on quadrilateral and triangular grids[J]. Applied Numerical Mathematics, 2004, 51(1):101-125.

[17]Sato Y, Hino T, Ohashi K. Parallelization of an unstructured Navier-Stokes solver using a multi-color ordering method for OpenMP[J]. Computers & Fluids, 2013, 88(1):496-509.

[18]Saad Y. Iterative methods for sparse linear systems, second edition[M]. Philadelphia: Society for Industrial and Applied Mathematics, 2003.

[19]Jameson A, Schmidt W, Turkel E. Numerical solution of the Euler equations by finite volume methods using Runge-Kutta time-stepping schemes. AIAA 1981-1259[R]. Reston: AIAA, 1981.

[20]Haselbacher A, Blazek J.Accurate and efficient discretization of Navier-Stokes equations on mixed grids[J]. AIAA Journal, 2000, 38(11):2094-2102.

[21]Economon T D, Palacios F, Copeland S R, et al. SU2: an open-source suite for multiphysics simulation and design[J]. AIAA Journal, 2016, 54(3):828-846.

[22]Cook P, Mcdonald M, Firmin M. Aerofoil RAE 2822-pressure distributions, and boundary layer and wake measurements. AGARD Report AR 138[R]. AGARD, 1979.

[23]Schmitt V,Charpin F. Pressure distribution on the ONERA-M6-Wing at transonic Mach numbers. AGARD Report AR 138[R]. AGARD, 1979.

[24]Hicken J E, Zingg D W. Dual consistency and functional accuracy: a finite-difference perspective[J]. Journal of Computational Physics, 2014, 256:161-182.

A parallel full implicit flow solver based on weakly imposed wall boundary condition

Gao Yisheng, Wu Yizhao*, Xia Jian, Xu Zhaoke

(NanjingUniversityofAeronautics&Astronautics,Nanjing210016,China)

A parallel full implicit Navier-Stokes solver based on weakly imposed wall boundary condition is proposed. First, a weak boundary procedure for viscous solid wall is implemented on the unstructured node-centered finite volume method. Upon this procedure, the corresponding viscous Jacobian is derived. Compared to the traditional strongly imposed wall boundary condition, the application of weakly imposed boundary condition eases the derivation. To address the issues of convergence and positivity preserving difficulties with regard to the Spalart-Allmaras (S-A) model, an unconditionally positive-convergent (UPC) implicit scheme for the S-A model using weakly imposed wall boundary condition is presented, which is based on the construction of a special implicit matrix to ensure both the convergence of the S-A model and the positivity of the turbulence working variables. The linear equations resulted from the full implicit time integration and the UPC scheme are solved using multicolor Gauss-Seidel iteration (MCGS) and MPI is employed for parallelization. The proposed approach is tested by simulating several two-dimensional and three-dimensional cases. Results from the numerical simulations demonstrate that: 1) with larger penalty strength parameter, the results of weak boundary conditions tend to those of strong boundary conditions; 2) the weak boundary conditions satisfy flux consistency for force computation automatically, and therefore no special treatment is required.

weakly imposed boundary condition; Navier-Stokes equations; full implicit method; positive scheme; turbulence model

0258-1825(2017)01-0046-11

2016-10-10;

2016-12-04

江蘇高校優(yōu)勢(shì)學(xué)科建設(shè)工程資助項(xiàng)目(PAPD)

高宜勝(1984-)，男，福建福州人，博士研究生，研究方向：計(jì)算流體力學(xué). E-mail：gaoyisheng@nuaa.edu.cn

伍貽兆(1945-)，男，江蘇六合人，教授，研究方向：計(jì)算流體力學(xué). E-mail：wyzao@nuaa.edu.cn.

高宜勝, 伍貽兆, 夏健, 等. 基于物面弱邊界條件的并行全隱式流場(chǎng)求解方法[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2017, 35(1): 46-56.

10.7638/kqdlxxb-2016.0121 Gao Y S, Wu Y Z, Xia J. Key technique and aerodynamic configuration characteristic of UCAV with command of the air[J]. Acta Aerodynamica Sinica, 2017, 35(1): 46-56.

V211.3

A doi: 10.7638/kqdlxxb-2016.0121