干雨新, 趙 寧

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

0 引 言

近年來(lái)，隨著計(jì)算機(jī)技術(shù)的迅猛發(fā)展，數(shù)值模擬方法成為了人們科學(xué)研究的重要熱點(diǎn)。在流體力學(xué)領(lǐng)域，計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)也成為了一門新興且熱門的學(xué)科。使用CFD方法進(jìn)行數(shù)值模擬，對(duì)計(jì)算區(qū)域的空間離散即網(wǎng)格生成技術(shù)是關(guān)鍵。網(wǎng)格質(zhì)量的好與壞，對(duì)計(jì)算結(jié)果的影響非常大，尤其是對(duì)于復(fù)雜外形的網(wǎng)格生成，即便對(duì)于一名熟練掌握網(wǎng)格生成技術(shù)的學(xué)者也是一件耗時(shí)耗力的事情。笛卡爾網(wǎng)格[1]應(yīng)運(yùn)而生，其相比于傳統(tǒng)的貼體網(wǎng)格，具有網(wǎng)格結(jié)構(gòu)簡(jiǎn)單、便于自動(dòng)化生成、人工干預(yù)少、實(shí)現(xiàn)自適應(yīng)容易等優(yōu)點(diǎn)，尤其是善于處理多物體和復(fù)雜幾何外形的網(wǎng)格生成問(wèn)題。對(duì)于N-S方程的求解，考慮流體的黏性作用，需要生成極薄的附面層網(wǎng)格。傳統(tǒng)的浸入邊界的各向同性笛卡爾網(wǎng)格方法，考慮到附面層網(wǎng)格的厚度，需要在附面層里布置大量的小尺度網(wǎng)格單元，這大大增加了網(wǎng)格數(shù)量(尤其對(duì)于三維問(wèn)題布置的大長(zhǎng)寬比六面體單元而言)。所以，在計(jì)算物體周圍生成貼體的附面層結(jié)構(gòu)網(wǎng)格，而在其他計(jì)算區(qū)域用笛卡爾網(wǎng)格填充的混合笛卡爾網(wǎng)格[2]成為了一種有效的解決問(wèn)題的手段。含有重疊區(qū)域的混合網(wǎng)格的難點(diǎn)在于混合網(wǎng)格的生成以及重疊區(qū)域的處理，Wang[3]使用重疊自適應(yīng)笛卡爾網(wǎng)格模擬了固定和運(yùn)動(dòng)邊界的流動(dòng)問(wèn)題；張來(lái)平等人[4]對(duì)于結(jié)構(gòu)非結(jié)構(gòu)和笛卡爾的混合動(dòng)網(wǎng)格技術(shù)進(jìn)行了深入的研究；肖天航等[5]則應(yīng)用變形重疊網(wǎng)格對(duì)非定常流動(dòng)進(jìn)行了數(shù)值模擬。

通常認(rèn)為，馬赫數(shù)低于0.3的流動(dòng)問(wèn)題，屬于不可壓流動(dòng)問(wèn)題。對(duì)于不可壓N-S方程的求解，直接采用可壓縮N-S方程求解的基于密度的方法由于速度過(guò)低，會(huì)導(dǎo)致方程的幾個(gè)特征值的量級(jí)差距過(guò)大，從而導(dǎo)致方程的剛性過(guò)大，使得求解收斂緩慢且結(jié)果錯(cuò)誤。而使用基于壓力的不可壓方法比如SIMPLE算法等則無(wú)法和可壓縮求解方法得到有效的統(tǒng)一。所以，對(duì)基于密度的可壓縮N-S方程求解算法進(jìn)行預(yù)處理，使方程的特征值處在同一個(gè)量級(jí)，從而使可壓縮N-S方程求解器能求解低速不可壓流場(chǎng)具備很大的意義。Turkel[6]經(jīng)過(guò)研究指出，對(duì)于低速流動(dòng)，預(yù)處理方法可以通過(guò)對(duì)時(shí)間導(dǎo)數(shù)的預(yù)處理達(dá)到加速計(jì)算收斂的效果，并且改變中心格式或迎風(fēng)格式的耗散來(lái)提高計(jì)算的精度。Weiss和Smith[7]則提出了基于N-S方程守恒變量的預(yù)處理矩陣，從而將基于密度的方法應(yīng)用到了低速黏性計(jì)算中。國(guó)內(nèi)也有很多學(xué)者對(duì)低速預(yù)處理方法進(jìn)行了研究[8-10]。本文使用基于混合笛卡爾網(wǎng)格的預(yù)處理方法，對(duì)NACA0012翼型在極低馬赫數(shù)下的定常繞流和動(dòng)態(tài)失速進(jìn)行了數(shù)值模擬，驗(yàn)證了所發(fā)展的低速預(yù)處理方法的正確性。耦合預(yù)處理技術(shù)的混合笛卡爾網(wǎng)格方法具備高效求解定常/非定常流動(dòng)的能力。

1 混合笛卡爾網(wǎng)格方法

1.1 混合笛卡爾網(wǎng)格生成

生成計(jì)算物體的貼體結(jié)構(gòu)網(wǎng)格，選定計(jì)算區(qū)域，并以貼體結(jié)構(gòu)網(wǎng)格的外邊界為幾何邊界，生成背景笛卡爾網(wǎng)格，兩套網(wǎng)格相互重疊，并在附近使用了加密的笛卡爾網(wǎng)格來(lái)盡量使得兩套網(wǎng)格系統(tǒng)之間的尺度過(guò)渡平緩。以NACA0012翼型為例，生成的混合笛卡爾網(wǎng)格如圖1所示。

圖1 NACA0012翼型的混合笛卡爾網(wǎng)格局部放大Fig.1 Local enlarged hybrid Cartesian grid of NACA0012 airfoil

在生成混合笛卡爾網(wǎng)格的過(guò)程中，本文首先在計(jì)算區(qū)域全場(chǎng)生成背景笛卡爾網(wǎng)格，然后根據(jù)貼體結(jié)構(gòu)網(wǎng)格的最外層(稱之為“交界面”)來(lái)進(jìn)行笛卡爾網(wǎng)格單元的分類，判斷其為“洞外單元”、“洞邊界單元”或“洞內(nèi)單元”，完成“挖洞”操作[11]。為了提高搜索效率，可以將貼體結(jié)構(gòu)網(wǎng)格存儲(chǔ)為二叉樹(shù)的數(shù)據(jù)結(jié)構(gòu)，使用了基于ADT技術(shù)[12]的單元搜尋方法，來(lái)快速找出洞內(nèi)單元。

1.2 貢獻(xiàn)單元查找

對(duì)于笛卡爾網(wǎng)格和貼體結(jié)構(gòu)網(wǎng)格的重疊區(qū)域，在進(jìn)行交界面上通量計(jì)算過(guò)程中，左單元在各自的網(wǎng)格系統(tǒng)里很容易確定，但右單元需要到另外一套網(wǎng)格系統(tǒng)里尋找，即尋找所謂的“貢獻(xiàn)單元”。本文采用Munikrishna和Liou給出的“貢獻(xiàn)單元”選取方法[13]?！柏暙I(xiàn)單元”的尋找過(guò)程中同樣使用了ADT技術(shù)來(lái)加速搜索過(guò)程。

1.3 非定常運(yùn)動(dòng)問(wèn)題的新現(xiàn)單元處理

在計(jì)算物面邊界運(yùn)動(dòng)的非定常問(wèn)題中，由于貼體網(wǎng)格隨著物面一起運(yùn)動(dòng)，所以原先被覆蓋的標(biāo)記為洞內(nèi)單元的笛卡爾網(wǎng)格會(huì)顯現(xiàn)出來(lái)，成為新現(xiàn)的網(wǎng)格單元，參與流體計(jì)算，所以對(duì)于新現(xiàn)單元的處理顯得尤為重要。

本文對(duì)于新現(xiàn)網(wǎng)格單元的具體處理方法可以參考文獻(xiàn)[14]，如圖2所示，貼體結(jié)構(gòu)網(wǎng)格從圖中黑色位置運(yùn)動(dòng)到紅色位置，則網(wǎng)格單元I為新現(xiàn)網(wǎng)格單元。對(duì)于新現(xiàn)網(wǎng)格單元I上流場(chǎng)信息的確定，主要做法是首先尋找新現(xiàn)網(wǎng)格單元I周圍的原洞邊界單元(即圖2中的網(wǎng)格單元1，2，3)，將這些單元作為插值的模板單元使用權(quán)系數(shù)為逆距離的插值公式計(jì)算新現(xiàn)網(wǎng)格單元上的流場(chǎng)信息，插值公式如下：

圖2 非定常計(jì)算新現(xiàn)單元的處理方法Fig.2 Method of new cell on unsteady calculation

2 控制方程與預(yù)處理方法

2.1 控制方程

積分形式的非定常N-S方程如下：

其中，

Vr為逆變速度V=V·n相對(duì)于運(yùn)動(dòng)網(wǎng)格的速度，其表達(dá)式為Vr=V-Vt=nxu+nyv+nzw-Vt,其中Vt為網(wǎng)格運(yùn)動(dòng)的速度。

2.2 預(yù)處理方法

首先將對(duì)于守恒變量W的N-S方程轉(zhuǎn)化成關(guān)于原始變量Q=(P,u,v,w,T)的形式：

預(yù)處理方法則是將雅克比矩陣M替換為預(yù)處理矩陣Γ，本文選用Weiss-Smith預(yù)處理矩陣，矩陣形式如下：

其中，

Ur=min(c,max(V,KV∞))(9)

其中，V∞為來(lái)流速度，常數(shù)K取值為1～4，對(duì)于黏性流動(dòng)，參考速度Ur還需要滿足：

其中，Δd為網(wǎng)格尺度。

引入預(yù)處理矩陣Γ后，N-S方程變?yōu)椋?/p>

則此時(shí)對(duì)流通量的雅克比矩陣的特征值為：

λ1,2,3=V,λ4,5=V′±c′

V′=V(1-a)

其中，V′和c′為偽速度和偽聲速。在低速條件下，參考速度Ur?c，偽速度V′和偽聲速c′都是V的量級(jí)，降低了方程的條件數(shù)，從而改善了方程的收斂性質(zhì)。當(dāng)參考速度Ur=c時(shí)，a=0，則偽速度和偽聲速等于真實(shí)的速度與聲速，方程的特征值自動(dòng)退回到預(yù)處理之前的形式，從而起到了關(guān)閉預(yù)處理的作用。

需要指出，使用預(yù)處理技術(shù)后，改變了N-S方程的特征系統(tǒng)，原本基于黎曼不變量的遠(yuǎn)場(chǎng)邊界條件將不再適用，需要根據(jù)新的特征值來(lái)計(jì)算新的黎曼不變量。為簡(jiǎn)化遠(yuǎn)場(chǎng)邊界的處理，本文采用Turkel等提出的邊界條件[15]，對(duì)于入流遠(yuǎn)場(chǎng)邊界：

ub=u∞,vb=v∞,wb=w∞,Tb=T∞,Pb=Pint(13)

對(duì)于出流遠(yuǎn)場(chǎng)邊界：

ub=uint,vb=vint,wb=wint,Tb=Tint,Pb=P∞(14)

其中，下標(biāo)b為遠(yuǎn)場(chǎng)邊界處的流場(chǎng)信息，下標(biāo)int為內(nèi)部單元的流場(chǎng)信息，下標(biāo)∞為來(lái)流的數(shù)值。

同時(shí)，由于預(yù)處理方法改變了方程的特征系統(tǒng)，對(duì)流通量的求解也要相應(yīng)的改變。本文對(duì)流通量求解使用HLLC格式[16]，進(jìn)行預(yù)處理后只需將波速計(jì)算中的特征值改為低速預(yù)處理后的特征值即可。

2.3 非定常流動(dòng)數(shù)值方法

對(duì)于非定常流動(dòng)，通常采用雙時(shí)間步對(duì)N-S方程的時(shí)間導(dǎo)數(shù)項(xiàng)進(jìn)行離散，本文采用雙時(shí)間步的LUSGS[17]預(yù)處理方法，物理時(shí)間采用二階向后差分離散，虛擬時(shí)間采用一階差分離散，離散后的方程為：

時(shí)間推進(jìn)采用隱式LUSGS方法，引入低速預(yù)處理后求解過(guò)程如下：

向前掃：

(MΓ-1ΔF*-λijΔW*)](16)

向后掃：

(MΓ-1ΔF-λijΔW)(17)

3 算例與分析

3.1 NACA0012翼型定常繞流

對(duì)NACA0012翼型在馬赫數(shù)為0.00001、0.0001、0.001、0.01、0.1、0.3、0.5和0.6這八個(gè)狀態(tài)進(jìn)行了數(shù)值計(jì)算，迎角α=3.59°，雷諾數(shù)Re=1.85×106，使用SSTk-ω湍流模型[18]，計(jì)算使用的混合笛卡爾網(wǎng)格如圖1所示。

該算例中，研究了不同馬赫數(shù)下預(yù)處理對(duì)計(jì)算收斂性的影響。從圖3可以看出，在馬赫數(shù)從0.01到0.3的低速范圍內(nèi)，預(yù)處理技術(shù)起到了非常明顯的加速收斂作用。在不使用預(yù)處理的情況下，隨著馬赫數(shù)的降低，收斂性能逐漸變差，當(dāng)馬赫數(shù)為0.01時(shí)，計(jì)算基本無(wú)法收斂，這是由于在馬赫數(shù)很低時(shí)，系統(tǒng)剛性明顯，方程組的耦合性很差，導(dǎo)致收斂性過(guò)差。而使用預(yù)處理技術(shù)后，系統(tǒng)剛性問(wèn)題得到有效改善，收斂性明顯提高，而且收斂的殘差下降過(guò)程基本和馬赫數(shù)無(wú)關(guān)。圖4給出了在三種極低馬赫數(shù)下的加入預(yù)處理技術(shù)后的殘差收斂曲線，同樣可以看出殘差的下降過(guò)程在收斂到“機(jī)器零”前和馬赫數(shù)無(wú)關(guān)。圖3和圖4表明，在低速情況下，使用預(yù)處理技術(shù)可以大大改善方程的剛性問(wèn)題，使收斂性得到大幅度提高。圖5為在馬赫數(shù)為0.5和0.6的亞、跨聲速范圍內(nèi)，預(yù)處理對(duì)流場(chǎng)計(jì)算的影響，此時(shí)由于馬赫數(shù)變高，預(yù)處理的影響變小，在高馬赫數(shù)情況下，方程具有良好的耦合性，預(yù)處理的影響也逐漸變小，但依然對(duì)計(jì)算收斂的加速有一定影響，因?yàn)榇藭r(shí)流場(chǎng)中有一部分的低速區(qū)，比如翼型附面層和駐點(diǎn)附近[19]。該算例表明，本文開(kāi)發(fā)的基于混合笛卡爾網(wǎng)格的預(yù)處理方法，將基于密度的可壓縮N-S方程求解方法拓展到了從低速(極低馬赫數(shù))到常規(guī)馬赫數(shù)的流場(chǎng)，使用預(yù)處理方法后可以達(dá)到加速收斂和提高計(jì)算精度的目的，這和Turkel的結(jié)論是一致的。

圖3 預(yù)處理對(duì)NACA0012翼型低速計(jì)算收斂性的影響Fig.3 Effect of precondition on the convergence of NACA0012 airfoil at low speed

圖4 預(yù)處理對(duì)NACA0012翼型極低馬赫數(shù)計(jì)算收斂性的影響Fig.4 Effect of precondition on the convergence of NACA0012 airfoil at very low Mach number

圖5 預(yù)處理對(duì)NACA0012翼型亞、跨聲速計(jì)算收斂性的影響Fig.5 Effect of precondition on the convergence of NACA0012 airfoil at subsonic and transonic simulations

3.2 NACA0012翼型的低速俯仰震蕩

在翼型的低速俯仰震蕩問(wèn)題中，尤其在迎角超過(guò)翼型的失速迎角后，翼型的背風(fēng)面流動(dòng)會(huì)發(fā)生分離從而動(dòng)態(tài)失速，此時(shí)翼型背風(fēng)面的流動(dòng)情況及其復(fù)雜，不少學(xué)者對(duì)該問(wèn)題做了數(shù)值計(jì)算[20]與實(shí)驗(yàn)進(jìn)行研究，選用Lee等[21]做的典型的低速實(shí)驗(yàn)狀態(tài)來(lái)驗(yàn)證本文的非定常低速預(yù)處理方法。該狀態(tài)下來(lái)流馬赫數(shù)為Ma=0.0385，雷諾數(shù)為Re=1.35×105，使用SSTk-ω兩方程湍流模型，迎角變化規(guī)律為α(t)=αm+α0sin(ωt)，其中αm=10°，α0=15°，縮減頻率為k=ωt/(2V∞)=0.1。網(wǎng)格與前文保持一致，貼體結(jié)構(gòu)網(wǎng)格部分隨著翼型一起做俯仰震蕩，交界面新現(xiàn)單元的處理在2.3節(jié)已經(jīng)做了闡述。

使用混合笛卡爾網(wǎng)格的非定常預(yù)處理方法，對(duì)NACA0012翼型在上述狀態(tài)下進(jìn)行了6個(gè)周期的俯仰震蕩計(jì)算，圖6為非定常計(jì)算雙時(shí)間步中第一個(gè)物理時(shí)間步下的收斂曲線，從曲線中可以看出預(yù)處理技術(shù)可以有效改善非定常計(jì)算偽時(shí)間迭代上的系統(tǒng)剛性，起到加速收斂的效果。

圖7給出了翼型俯仰震蕩過(guò)程中幾個(gè)迎角下的預(yù)處理結(jié)果和Lee等的實(shí)驗(yàn)[21]觀察到的流場(chǎng)的比較，從中可以看出，一開(kāi)始流動(dòng)附著在翼型表面，隨著迎角的增加，翼型上表面前緣逐漸出現(xiàn)前緣渦(LEV)，然后隨著迎角繼續(xù)增大，前緣渦也逐漸變大并往翼型

圖6 預(yù)處理對(duì)NACA0012翼型低馬赫數(shù)非定常計(jì)算收斂性的影響Fig.6 Effect of precondition on the convergence of NACA0012 airfoil′s unsteady calculation at low Mach number

圖7 NACA0012翼型低速動(dòng)態(tài)失速過(guò)程中不同時(shí)刻的流場(chǎng)Fig.7 NACA0012 airfoil dynamic stall at low speed and different time

上表面后緣移動(dòng)，升力也逐漸增大。隨著迎角增加到24.7°，前緣渦破裂，翼型上表面流動(dòng)完全分離，升力迅速下降。接著迎角逐漸變小，翼型上表面的流動(dòng)也隨著迎角變小而逐漸慢慢附著到翼型表面上，直到最后完全附著。

圖8和圖9為NACA0012翼型俯仰震蕩升力系數(shù)和力矩系數(shù)的遲滯曲線，從圖中可以看出，數(shù)值模擬結(jié)果和實(shí)驗(yàn)結(jié)果雖然稍有差別，但較為準(zhǔn)確的模擬出了翼型動(dòng)態(tài)失速的規(guī)律。在失速前的上仰階段，翼型氣動(dòng)力和實(shí)驗(yàn)結(jié)果符合較好，但隨著迎角增大，翼型失速后的升力系數(shù)和力矩系數(shù)均與實(shí)驗(yàn)值有所偏差，但和未加預(yù)處理的結(jié)果相比，和實(shí)驗(yàn)值的誤差大大降低。不進(jìn)行預(yù)處理的計(jì)算結(jié)果無(wú)論在數(shù)值上還是在變化趨勢(shì)上和實(shí)驗(yàn)結(jié)果相去甚遠(yuǎn)，預(yù)處理的結(jié)果雖然和實(shí)驗(yàn)結(jié)果稍有偏差，但依然能較為準(zhǔn)確的捕捉到動(dòng)態(tài)失速的規(guī)律。圖8和圖9中紅色點(diǎn)劃線為變形網(wǎng)格方法[22]的計(jì)算結(jié)果，可以看到，其模擬結(jié)果也和實(shí)驗(yàn)結(jié)果存在偏差，而且從圖中可以看出本文用混合笛卡爾網(wǎng)格預(yù)處理的模擬結(jié)果和實(shí)驗(yàn)更為接近。

圖8 升力系數(shù)隨迎角的遲滯曲線Fig.8 Lift coefficient hysteresis curves

圖9 力矩系數(shù)隨迎角的遲滯曲線Fig.9 Moment coefficient hysteresis curves

4 結(jié) 論

本文使用混合笛卡爾網(wǎng)格方法，對(duì)NACA0012翼型的低速定常繞流和非定常俯仰震蕩進(jìn)行了數(shù)值模擬研究，得到結(jié)論如下：

1) 使用混合笛卡爾網(wǎng)格的隱式雙時(shí)間步預(yù)處理方法對(duì)翼型的定常/非定常數(shù)值模擬是可行的，使用預(yù)處理方法后可大大加快翼型在低馬赫數(shù)下數(shù)值計(jì)算的收斂速度，并提高計(jì)算的準(zhǔn)確性。

2) 使用預(yù)處理方法后，隨著馬赫數(shù)的升高，預(yù)處理的影響在逐漸去除，但對(duì)收斂依然有著一定的加速作用，原因是翼型附面層或駐點(diǎn)附近仍然存在著一部分的低速區(qū)。

3) 使用混合笛卡爾網(wǎng)格方法可以較為簡(jiǎn)便的處理邊界運(yùn)動(dòng)的非定常問(wèn)題，并且避免了在非定常運(yùn)動(dòng)中的網(wǎng)格單元變形和運(yùn)動(dòng)導(dǎo)致的單元負(fù)體積產(chǎn)生，同時(shí)也避免了因網(wǎng)格運(yùn)動(dòng)項(xiàng)的引入而產(chǎn)生的幾何守恒率問(wèn)題[23]，只需要在交界面進(jìn)行新現(xiàn)單元的處理即可。

4) 本文將混合笛卡爾網(wǎng)格方法和預(yù)處理技術(shù)相結(jié)合，開(kāi)發(fā)了一套基于混合笛卡爾網(wǎng)格的從低速(極低馬赫數(shù))到常規(guī)馬赫數(shù)的流場(chǎng)求解器，將混合笛卡爾網(wǎng)格方法用于不可壓?jiǎn)栴}的計(jì)算，結(jié)果和實(shí)驗(yàn)值吻合良好，顯示了基于混合笛卡爾網(wǎng)格的預(yù)處理方法對(duì)從低速到常規(guī)馬赫數(shù)流動(dòng)問(wèn)題的求解是有效的。