陳浩，袁先旭，王田天，周丹，趙寧，唐志共，畢林,*

1. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點(diǎn)實驗室，綿陽 621000

2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

3. 中南大學(xué) 交通運(yùn)輸工程學(xué)院 軌道交通安全教育部重點(diǎn)實驗室，長沙 410075

4. 南京航空航天大學(xué) 航空學(xué)院，南京 210016

隨著計算機(jī)技術(shù)和計算流體力學(xué)學(xué)科的飛速發(fā)展，數(shù)值模擬在湍流研究中扮演著越來越重要的角色。然而，由于湍流流動具有非定常、多尺度、強(qiáng)擴(kuò)散與耗散等復(fù)雜特性，且目前面臨的幾何外形越來越復(fù)雜[1-2]，數(shù)值模擬要在認(rèn)識、解決、應(yīng)用湍流問題上發(fā)揮更重要的作用，還需攻克高質(zhì)量計算網(wǎng)格生成的難題。

目前，以多塊分區(qū)結(jié)構(gòu)網(wǎng)格[3-5]、三角形(四面體)貼體非結(jié)構(gòu)網(wǎng)格[6-7]等為基礎(chǔ)發(fā)展的湍流流動數(shù)值模擬方法，在提高湍流認(rèn)知水平、解決湍流工程應(yīng)用問題等方面發(fā)揮了重要作用。然而，這些網(wǎng)格類型對于湍流這樣的復(fù)雜流動，要么需要花費(fèi)大量的人力資源，自動化程度不高，要么網(wǎng)格分布難以滿足復(fù)雜湍流流動結(jié)構(gòu)動態(tài)、精細(xì)捕捉要求，自適應(yīng)能力不強(qiáng)，逐漸成為制約湍流精細(xì)化研究的瓶頸問題之一。

相對于傳統(tǒng)的分區(qū)結(jié)構(gòu)網(wǎng)格和貼體非結(jié)構(gòu)網(wǎng)格，自適應(yīng)笛卡爾網(wǎng)格[8-11]不依賴于物面網(wǎng)格直接生成空間網(wǎng)格，具有網(wǎng)格生成自動化程度高、復(fù)雜外形適應(yīng)性好、非定常/多尺度等流動結(jié)構(gòu)捕捉能力強(qiáng)等優(yōu)勢，適用于湍流等復(fù)雜流動問題的仿真模擬。因此自適應(yīng)笛卡爾網(wǎng)格方法具有重要的研究意義和廣闊的應(yīng)用前景。

笛卡爾網(wǎng)格也有其劣勢，主要是非貼體性帶來的物面邊界處理困難。對于非貼體的笛卡爾網(wǎng)格來說，網(wǎng)格與物面邊界相交，在處理黏性流動問題時面臨較大的局限性，尤其是對于邊界層內(nèi)大縱橫比流動，雷諾數(shù)越高，網(wǎng)格規(guī)模越大。因此，通過笛卡爾網(wǎng)格方法實現(xiàn)三維空間效應(yīng)強(qiáng)、旋渦尺度跨度大的非定常復(fù)雜流動問題的高精度求解具有一定難度。針對這種非貼體性，不同學(xué)者提出了不同的邊界處理方式，主要有：切割單元法和浸入邊界法。

切割單元法使笛卡爾網(wǎng)格具有貼體網(wǎng)格的特性，在處理時需要考慮網(wǎng)格被物面切割后的形狀和分布。該方法能使物面處滿足守恒性的要求，一般結(jié)合有限體積方法進(jìn)行計算，應(yīng)用較為廣泛[12-13]。商業(yè)軟件CART3D就是基于切割單元法進(jìn)行處理非貼體物面邊界的。中國空氣動力研究與發(fā)展中心的肖涵山等[14]對切割單元法進(jìn)行了研究和發(fā)展，并成功應(yīng)用于機(jī)彈分離問題的模擬仿真。雖然切割單元法在一些計算中取得了良好的結(jié)果，但是這種方法會產(chǎn)生形狀多樣、大小不一的網(wǎng)格單元，這樣不僅會限制時間步長，還會增加數(shù)值解的非物理震蕩，進(jìn)而增加求解難度，限制了其在黏性流動問題中的應(yīng)用。

浸入邊界法[15-16]則不需要對在物面附近的單元進(jìn)行特殊處理，保留了笛卡爾網(wǎng)格原有形狀，故此方法有一定的便捷性和靈活性。該方法最初由Peskin[16]在1972年提出，用于模擬人類心臟中的血液流動。這種方法的思想是：不直接考慮物面的存在，而是把物體作為邊界條件嵌入到流場中，對邊界的處理則是將邊界?；闪鲃涌刂品匠讨械牡牧υ错棥ohd[17]將浸入邊界法與B-Spline法結(jié)合來獲取高階精度，并且通過用時間離散導(dǎo)數(shù)去掉時間步長的限制，大大擴(kuò)展了浸入邊界法的使用范圍。Dadone和Grossman[18-20]在浸入邊界法的基礎(chǔ)上，提出了虛擬單元方法(Ghost Cell Method)，該方法是在Dadone提出的物面曲率修正技術(shù)(Curvature Corrected Symmetry Technique)的基礎(chǔ)上演化而來，最初應(yīng)用于貼體網(wǎng)格，之后拓展到笛卡爾網(wǎng)格，他成功地將這種方法應(yīng)用于二維和三維無黏流動問題的數(shù)值模擬，并在之后結(jié)合數(shù)值壁面函數(shù)，實現(xiàn)了對于高雷諾數(shù)可壓縮黏性流動問題的模擬。在上述工作的基礎(chǔ)上，Lee和Ruffin[21]將浸入邊界方法推廣到黏性問題的流場計算中，提出了一種新的壁面函數(shù)，并針對湍流問題進(jìn)行了模擬研究，取得了較好的效果。國內(nèi)在浸入邊界法和自適應(yīng)笛卡爾網(wǎng)格技術(shù)相結(jié)合的研究領(lǐng)域成果則相對較少，劉劍明[22]、胡偶[23]對自適應(yīng)笛卡爾網(wǎng)格技術(shù)和浸入邊界法進(jìn)行了發(fā)展，在黏性可壓縮流動問題中作了應(yīng)用研究，并擴(kuò)展到了高雷諾數(shù)黏性流動。陳浩等[24]將上述技術(shù)方法進(jìn)一步擴(kuò)展至超聲速和高超聲速黏性流動問題，在二維典型算例中得到了可靠的結(jié)果。上述研究所面臨的外形相對簡單，對于復(fù)雜構(gòu)型的擴(kuò)展應(yīng)用仍然有待開展。

為了保證求解精度，笛卡爾網(wǎng)格在自適應(yīng)加密過程會產(chǎn)生大量網(wǎng)格。尤其當(dāng)求解三維黏性流動問題時，為降低非物理振蕩、保證模擬精度，非貼體的笛卡爾網(wǎng)格會在邊界層區(qū)域進(jìn)行不計成本的加密，使得相同模擬精度下，笛卡爾網(wǎng)格明顯大于貼體類網(wǎng)格，甚至?xí)辛考壣系牟罹唷＞薮蟮木W(wǎng)格量會導(dǎo)致計算效率低下以及存儲量過大的難題。在應(yīng)對上述問題時，不同學(xué)者形成了兩種研究思路：

1) 在網(wǎng)格技術(shù)上做改進(jìn)，主要包括：混合笛卡爾網(wǎng)格技術(shù)、各向異性自適應(yīng)技術(shù)以及法向射線加密技術(shù)等。

其中，混合笛卡爾網(wǎng)技術(shù)相對使用較多，將貼體結(jié)構(gòu)網(wǎng)格和笛卡爾網(wǎng)格的優(yōu)勢相結(jié)合，通過在近壁使用貼體網(wǎng)格來避免笛卡爾網(wǎng)格邊界處理的困難。網(wǎng)格生成與空氣動力研究國際權(quán)威普林斯頓大學(xué)Baker[25]認(rèn)為混合網(wǎng)格技術(shù)是易用度和黏性模擬精度兼顧較好的方法。Charlton等[26]較早提出了內(nèi)部貼體網(wǎng)格和外部直角坐標(biāo)網(wǎng)格相結(jié)合的混合網(wǎng)格思想。之后，Delanaye[27]在Beger[28]與Aftosim等[29]的工作基礎(chǔ)上發(fā)展了貼體結(jié)構(gòu)與笛卡爾網(wǎng)格相混合的Navier-Stokes(N-S)求解器。目前，國外已發(fā)展了多個基于笛卡爾/結(jié)構(gòu)網(wǎng)格混合的知名網(wǎng)格技術(shù)求解器，例如，NASA著名的OVERFLOW[30]求解器就是使用笛卡爾背景網(wǎng)格耦合貼體結(jié)構(gòu)網(wǎng)格的形式。此外，由美國國防部支持開發(fā)的旋翼計算平臺HELIOS，是基于笛卡爾/貼體非結(jié)構(gòu)網(wǎng)格混合的形式。Wang和Chen[31]提出的 “黏性自適應(yīng)笛卡爾網(wǎng)格方法” 也是混合網(wǎng)格方法，其是在切割單元方法基礎(chǔ)上，通過投影方式構(gòu)造近壁貼體結(jié)構(gòu)網(wǎng)格。之后，F(xiàn)ujimoto等[32]繼續(xù)對Wang的“黏性自適應(yīng)笛卡爾方法”進(jìn)行改進(jìn)，使之能夠處理包括縫隙、凹槽、凸起等結(jié)構(gòu)的更為復(fù)雜的幾何外形流動。

相比于國外形成的較為成熟的網(wǎng)格求解技術(shù)與求解程序，國內(nèi)在基于笛卡爾網(wǎng)格的混合網(wǎng)格技術(shù)方面的研究相對有限，相關(guān)成果也較少。肖涵山等[33]在其發(fā)展自適應(yīng)笛卡爾網(wǎng)格技術(shù)基礎(chǔ)上，借鑒“投影”方法的思想，實現(xiàn)了混合笛卡爾網(wǎng)格的自動生成，同時構(gòu)建了Euler數(shù)值求解器，形成了飛行器氣動力快速預(yù)測和評估工具。張來平等[34]在混合網(wǎng)格方法上進(jìn)行了較為深入的研究，涉及到了混合笛卡爾網(wǎng)格方法的相關(guān)工作。他們針對不可壓非定常流動問題，在近壁面采用貼體結(jié)構(gòu)網(wǎng)格，遠(yuǎn)場采用笛卡爾網(wǎng)格，中間用非結(jié)構(gòu)網(wǎng)格連接過渡，建立了動態(tài)混合網(wǎng)格方法。田書玲等[35-37]發(fā)展了混合笛卡爾網(wǎng)格技術(shù)，并對網(wǎng)格生成、混合重疊網(wǎng)格裝配和動態(tài)重疊網(wǎng)格算法開展了深入研究。借鑒傳統(tǒng)重疊網(wǎng)格思想，胡偶[23]、沈志偉[38]在其所建立的自適應(yīng)混合笛卡爾網(wǎng)格方法基礎(chǔ)上，開展了非定常旋渦主導(dǎo)流動問題的數(shù)值模擬與流動機(jī)理分析。以上工作顯著促進(jìn)了國內(nèi)混合笛卡爾網(wǎng)格技術(shù)的發(fā)展，部分成果已應(yīng)用于相對復(fù)雜的流動問題研究。

各向異性笛卡爾網(wǎng)格技術(shù)是利用流場中很多區(qū)域是各向異性的特點(diǎn)，發(fā)展與之相匹配的網(wǎng)格加密方法，實現(xiàn)在保證網(wǎng)格質(zhì)量的同時明顯降低網(wǎng)格規(guī)模。傳統(tǒng)的笛卡爾網(wǎng)格方法多是采用各項同性加密的方法，事實上，流場中的很多區(qū)域是各向異性的，如附面層、激波間斷、剪切層等，在這些區(qū)域如果還采用各向同性的加密方法，勢必會造成網(wǎng)格點(diǎn)的大量浪費(fèi)。因此，很多學(xué)者將各向異性加密技術(shù)引入到笛卡爾網(wǎng)格方法當(dāng)中。Ham等[39]針對不可壓縮層流非定常流動問題，基于非結(jié)構(gòu)數(shù)據(jù)結(jié)構(gòu)，建立了各向異性笛卡爾網(wǎng)格加密和粗化方法，顯著減少了網(wǎng)格數(shù)目和內(nèi)存占用。在此基礎(chǔ)上，Keats和Lien[40]將該方推廣到可壓縮黏性流動問題的模擬當(dāng)中，經(jīng)驗證，他們所發(fā)展的方法在降低網(wǎng)格數(shù)目的同時能夠保證計算精度。Wang和Chen[31]提出的 “黏性自適應(yīng)笛卡爾網(wǎng)格方法”正是采用了基于2N叉樹數(shù)據(jù)結(jié)構(gòu)的各向異性流場解自適應(yīng)方法和“投影”方法相結(jié)合，實現(xiàn)了對于湍流黏性流動問題較高精度的模擬。Wu和Li[41]也對各向異性的網(wǎng)格加密技術(shù)進(jìn)行了研究和發(fā)展，他們最初是用各向異性網(wǎng)格方法捕捉黏性邊界層中的固有各向異性特性，在經(jīng)修正之后，在捕捉斜激波和正激波方面取得了較好的效果。Wu和Li[42]分別將各向同性和各向異性的網(wǎng)格加密技術(shù)應(yīng)用到笛卡爾網(wǎng)格方法中，進(jìn)行對比研究，驗證了所發(fā)展的各向異性解自適應(yīng)方法加密笛卡爾網(wǎng)格可以有效地捕捉流場中激波、接觸間斷等關(guān)鍵區(qū)域。Sang和Shi[43]基于全叉樹(Omni-Tree)數(shù)據(jù)結(jié)構(gòu)發(fā)展了自適應(yīng)笛卡爾網(wǎng)格方法，相對于傳統(tǒng)的僅支持各向同性加密的八叉樹數(shù)據(jù)結(jié)構(gòu)，該數(shù)據(jù)結(jié)構(gòu)同時支持各向同性和各向異性兩種加密方式，通過對民用飛機(jī)這種復(fù)雜外形開展的模擬仿真和對比研究，證明其所發(fā)展的各向異性笛卡爾網(wǎng)格技術(shù)能夠在一定程度上減少網(wǎng)格數(shù)目和提高計算效率。

法向加密技術(shù)則是借鑒貼體結(jié)構(gòu)網(wǎng)格大縱橫比的特點(diǎn)，通過沿壁面法向配置的射線，調(diào)整附面層內(nèi)的笛卡爾網(wǎng)格的粗細(xì)分布，減少附面層內(nèi)網(wǎng)格數(shù)目。這種方法最初由Ruffin等[44]在2012年提出，他們針對邊界層高精度模擬需要大量網(wǎng)格的問題，提出邊界層法向加密(Normal-Ray Refinement)的概念，并建立了物面法向射線的分布方式以及不同射線之間加密網(wǎng)格的數(shù)據(jù)交互方法，極大減少了網(wǎng)格點(diǎn)數(shù)目和內(nèi)存占用，節(jié)省了計算時間。但是這種方法需要人為地指定法向加密的位置，在非加密區(qū)其壁面處理精度較低，一般用于求解層流黏性流動問題，對于邊界層特征很復(fù)雜的湍流問題的高精度模擬仍然需要發(fā)展。

2) 引入通過簡化方法刻畫近壁流動的思想，降低壁面網(wǎng)格尺度要求，以達(dá)到減少網(wǎng)格量的目的，主要包括近壁模型(Wall Model)和壁面函數(shù)(Wall Function)方法。

其中，近壁模型方法通過對湍流模型進(jìn)行修正，使其能夠?qū)τ诮陴ば杂绊憛^(qū)域進(jìn)行求解[45-46]。壁面函數(shù)方法則是考慮到邊界層存在解析解這一特征，通過構(gòu)造一維代數(shù)模型來代替黏性底層的直接求解，使得對壁面法向第1層網(wǎng)格的y+需求放寬至50～100，甚至200，大大降低計算資源的使用[21]。相比近壁模型方法，計算量相對較小，且更容易實現(xiàn)，在貼體類網(wǎng)格方法中得到了廣泛應(yīng)用。

國外研究人員將近壁模型和壁面函數(shù)方法在笛卡爾網(wǎng)格框架下進(jìn)行應(yīng)用，定義湍流壁面邊界條件，從而降低邊界層內(nèi)計算網(wǎng)格數(shù)量。Kalitzin和Iaccarino[47]提出了一種基于壁面律的浸入邊界方法處理高雷諾數(shù)湍流問題。Lee和Ruffin[21]通過壁面函數(shù)模型耦合k-ω湍流模型定義湍流壁面邊界條件，并用于高雷諾數(shù)問題的數(shù)值模擬，包括二維平板、翼型以及三維旋翼機(jī)的機(jī)身旋翼相互作用等。Dadone和Grossman[18-20]在其無黏虛擬單元方法的基礎(chǔ)上，提出在壁面處構(gòu)造數(shù)值壁面函數(shù)結(jié)合虛擬單元的方法，求解高雷諾數(shù)流動問題。Capizzano[48]則采用雙層壁面模型修正浸入邊界網(wǎng)格交接面上湍流變量值的方法，進(jìn)行高雷諾數(shù)可壓縮流動問題的模擬研究。Bernardini等[49]基于浸入邊界方法，通過壁面模型耦合SA模型，用于延遲脫體渦模擬(Delayed Detached Eddy Simulation, DDES)中，對于高雷諾數(shù)下壁面分離的湍流流動開展了模擬研究，獲得了高精度、可靠的數(shù)據(jù)結(jié)果。這些學(xué)者的研究工作對基于笛卡爾網(wǎng)格技術(shù)進(jìn)行湍流問題的模擬起到了極大的推進(jìn)作用。

國內(nèi)有關(guān)笛卡爾網(wǎng)格方法在高雷諾數(shù)黏性流動問題中的研究和應(yīng)用起步較晚，成果也較少。劉劍明[22]基于自適應(yīng)笛卡爾技術(shù)，發(fā)展了針對可壓縮Navier-Stokes方程的徑向基函數(shù)浸入邊界方法，結(jié)合Menter 剪切應(yīng)力輸運(yùn)(SST)兩方程湍流模型，對于高雷諾數(shù)下的圓柱和NACA0012翼型等簡單外形的繞流問題進(jìn)行了研究，取得了一定的成果。胡偶[23]針對SSTk-ω模型湍流模型首次構(gòu)造了一種壁面函數(shù)——虛擬單元方法(Wall Function-Ghost Cell Method, WF-GCM)，實現(xiàn)了湍流壁面邊界重構(gòu)在“非貼體”笛卡爾網(wǎng)格上的應(yīng)用。他的研究對象比較簡單，雷諾數(shù)量級也不是很高，對于湍流邊界層的流動特征研究較少，有進(jìn)一步發(fā)展的需求和空間。

雖然國內(nèi)外學(xué)者對于黏性笛卡爾網(wǎng)格技術(shù)已經(jīng)開展了部分研究工作，促進(jìn)了笛卡爾網(wǎng)格方法在復(fù)雜黏性流動問題中的應(yīng)用。但是由于自適應(yīng)笛卡爾網(wǎng)格技術(shù)模擬黏性流動難度較大，且對于計算資源要求較高，國內(nèi)外對于黏性自適應(yīng)笛卡爾網(wǎng)格方法的研究相對其他網(wǎng)格技術(shù)而言仍然偏少，尤其是國內(nèi)，尚處于起步階段，公開的學(xué)術(shù)成果極少。因此，目前對于黏性自適應(yīng)笛卡爾網(wǎng)格方法國內(nèi)外尚缺少系統(tǒng)性的綜述介紹，亟待相關(guān)工作進(jìn)行補(bǔ)充。

國家數(shù)值風(fēng)洞(NNW)工程關(guān)注黏性自適應(yīng)笛卡爾網(wǎng)格方法的發(fā)展?jié)摿蛻?yīng)用前景，資助了若干研究課題，力爭突破關(guān)鍵技術(shù)瓶頸，促進(jìn)實現(xiàn)黏性笛卡爾網(wǎng)格技術(shù)的工程實用化。本文聚焦國家數(shù)值風(fēng)洞工程，總結(jié)介紹了黏性笛卡爾網(wǎng)格方法相關(guān)課題的研究進(jìn)展，重點(diǎn)從笛卡爾網(wǎng)格高效生成技術(shù)、笛卡爾網(wǎng)格自適應(yīng)技術(shù)以及黏性物面邊界處理等3個方面展開。

1 網(wǎng)格高效生成技術(shù)

1.1 笛卡爾網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)是指用來存儲如位置、物理參數(shù)、鄰居關(guān)系等網(wǎng)格節(jié)點(diǎn)信息的一種結(jié)構(gòu)形式，它是進(jìn)行數(shù)值計算的前提和基礎(chǔ)。合適的網(wǎng)格數(shù)據(jù)結(jié)構(gòu)可以有效減小內(nèi)存的需求并提高訪問速度，為網(wǎng)格的高效生成提供支撐。

目前笛卡爾網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)主要分為以下3種：結(jié)構(gòu)式、非結(jié)構(gòu)式和叉樹式[48-52]。結(jié)構(gòu)網(wǎng)格由于其相鄰網(wǎng)格之間存在自然的聯(lián)系，只需要采用I、J、K的順序?qū)?shù)據(jù)進(jìn)行存儲即可，使用起來簡單且易于編程實現(xiàn)，但靈活性較差，自適應(yīng)困難。非結(jié)構(gòu)網(wǎng)格在存儲自身網(wǎng)格流場信息的基礎(chǔ)上，還要存儲其鄰居網(wǎng)格單元的編號及其與自身網(wǎng)格相鄰邊的相對關(guān)系信息，導(dǎo)致其數(shù)據(jù)結(jié)構(gòu)復(fù)雜，存儲量大，訪問效率相對較低。叉樹數(shù)據(jù)結(jié)構(gòu)是一種分層式的數(shù)據(jù)結(jié)構(gòu)，依據(jù)網(wǎng)格大小劃分出不同的層次，由父單元即上一層單元決定其網(wǎng)格的位置，同樣的，該單元也可以劃分成數(shù)個子單元，并通過指針來建立不同層級之間父子單元的訪問和聯(lián)系。對于自適應(yīng)笛卡爾網(wǎng)格而言，網(wǎng)格的加密過程實質(zhì)上是對空間沿著垂直坐標(biāo)方向進(jìn)行等分(一維為二等分，二維為四等分，三維為八等分)，因此傳統(tǒng)的叉樹數(shù)據(jù)結(jié)構(gòu)在自適應(yīng)加密和粗化、網(wǎng)格調(diào)用方面具有很大的優(yōu)勢，但是其所存在的問題也比較突出：查詢網(wǎng)格單元的鄰居單元所需要的搜索計算量大，例如最差情況需要從葉子節(jié)點(diǎn)回溯到根節(jié)點(diǎn)，然后從根節(jié)點(diǎn)的鄰居節(jié)點(diǎn)搜索到葉子節(jié)點(diǎn)；在網(wǎng)格進(jìn)行自適應(yīng)后需重新確定每個網(wǎng)格單元的鄰居關(guān)系；存儲網(wǎng)格信息所需要的存儲量較大。此外，叉樹數(shù)據(jù)結(jié)構(gòu)還存在一些無用的內(nèi)存分配，例如，叉樹數(shù)據(jù)結(jié)構(gòu)的每個單元都有子單元指針，但是這些指針對于葉子單元而言并不會被用到，這就會造成內(nèi)存的浪費(fèi)[28]。

為了解決上述傳統(tǒng)叉樹數(shù)據(jù)結(jié)構(gòu)的問題，本文作者團(tuán)隊在笛卡爾網(wǎng)格存儲時采用了叉樹數(shù)據(jù)結(jié)構(gòu)的改進(jìn)形式FTT數(shù)據(jù)結(jié)構(gòu)[51-53]，如圖1[52]所示。它將叉樹數(shù)據(jù)結(jié)構(gòu)中葉子節(jié)點(diǎn)閑置的指針利用起來，指向鄰居單元的父單元，從而構(gòu)建快速訪問鄰居單元的“線程”，實現(xiàn)在不同檢索操作中數(shù)據(jù)信息的快速訪問，提供了高效、易并行的叉樹信息訪問路徑，同時節(jié)省數(shù)據(jù)節(jié)點(diǎn)，降低叉樹信息存儲的內(nèi)存占用,大大提高計算效率。

圖1 叉樹數(shù)據(jù)結(jié)構(gòu) vs FTT數(shù)據(jù)結(jié)構(gòu)[52]Fig.1 Fork tree data structure vs FTT data structure[52]

總體而言，采用FTT數(shù)據(jù)結(jié)構(gòu)進(jìn)行鄰居單元的查找，平均遍歷的數(shù)據(jù)層數(shù)是傳統(tǒng)的鄰居查找方式的1/4，這樣在搜索鄰居上可以節(jié)約大約75%的時間。在內(nèi)存占用方面，采用傳統(tǒng)的叉樹數(shù)據(jù)結(jié)構(gòu)，要存儲完上述信息，每個單元需要17 words的內(nèi)存，而如果采用FTT方式，每個單元需要2 words的內(nèi)存，相對于傳統(tǒng)方式明顯降低[52]。

1.2 物面網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

在生成網(wǎng)格的過程中，網(wǎng)格類型的判斷占據(jù)了大部分時間。為了高效判斷網(wǎng)格類型，Bonet和Peraire[54]采用ADT(Alternating Digital Tree)數(shù)據(jù)結(jié)構(gòu)來存儲物面網(wǎng)格，以縮小可能相交的物面單元檢測范圍，從而提高相交判斷的效率。在該方法中，每個ADT節(jié)點(diǎn)都可作為父節(jié)點(diǎn)或子節(jié)點(diǎn)對應(yīng)一個笛卡爾空間，父節(jié)點(diǎn)包含其所對應(yīng)所有子節(jié)點(diǎn)的笛卡爾空間[55]。ADT數(shù)據(jù)結(jié)構(gòu)利用以上的層次包含關(guān)系，只要父節(jié)點(diǎn)不與目標(biāo)空間相交，則該父節(jié)點(diǎn)所包含的子節(jié)點(diǎn)均將在查詢中被忽略，從而大大減少查詢次數(shù)，提高了搜索效率。然而傳統(tǒng)的ADT所構(gòu)建的二叉樹的平衡性往往受物面單元的分布密度影響，導(dǎo)致查詢效率不穩(wěn)定。為此，本文作者團(tuán)隊采用其改進(jìn)版本KDT(K-Dimensional Tree)結(jié)構(gòu)，其更容易構(gòu)造平衡優(yōu)質(zhì)的二叉樹[56-57]。

KDT是一種多維的二叉樹結(jié)構(gòu)，KD樹的每一層都和某個維度相關(guān)，每個節(jié)點(diǎn)都是一個k維節(jié)點(diǎn)。KDT的叉樹在構(gòu)建時，只需要將當(dāng)前層的相關(guān)點(diǎn)進(jìn)行排序后取其中值點(diǎn)，隨后確定經(jīng)過該點(diǎn)與相關(guān)維度垂直的平面,用來剖分空間的搜索區(qū)域，與ADT數(shù)據(jù)結(jié)構(gòu)在區(qū)域等分點(diǎn)進(jìn)行劃分不同，KDT可以保證多層叉樹的平衡性并減少樹的深度，從而便于查找單元信息。如圖2所示[51]，對于圖2(a) 中所示的物面網(wǎng)格，KDT叉樹數(shù)據(jù)結(jié)構(gòu)相比于ADT叉樹數(shù)據(jù)結(jié)構(gòu)更加平衡，其相應(yīng)的操作復(fù)雜度也明顯降低。對于物面網(wǎng)格較少的情況，宜采用ADT來存儲物面單元；當(dāng)物面網(wǎng)格量較多時，則宜采用KDT。這是因為，KDT每進(jìn)行一次空間劃分，均需要進(jìn)行一次排序操作，當(dāng)網(wǎng)格量較大時，樹的深度將會比ADT小很多，此時這種排序的優(yōu)勢才會凸顯。

圖2 隨機(jī)分布的網(wǎng)格點(diǎn)和ADT、KDT數(shù)據(jù)結(jié)構(gòu) [51]Fig.2 Randomly distributed grids and ADT、KDT data structure [51]

1.3 相交判斷方法

笛卡爾網(wǎng)格的相交判斷是開展自適應(yīng)和浸入邊界處理的前提?；诰W(wǎng)格相對于物面的位置關(guān)系，可將網(wǎng)格劃分為外部網(wǎng)格、內(nèi)部網(wǎng)格以及相交網(wǎng)格單元。實踐表明，相交判斷占據(jù)了笛卡爾網(wǎng)格生成的大量時間，并且是網(wǎng)格生成成敗的關(guān)鍵因素。

Aftosmis等[58]以及Dadone和Grossman[20]采用簡單且有效的射線相交方法進(jìn)行網(wǎng)格點(diǎn)的類型判斷。該方法通過計算由目標(biāo)點(diǎn)生成的射線與物面的交點(diǎn)個數(shù)，來判斷網(wǎng)格點(diǎn)與物面的相對位置：相交數(shù)目為奇數(shù)，網(wǎng)格點(diǎn)在物面內(nèi)部；相交數(shù)目為偶數(shù)，網(wǎng)格點(diǎn)在物面外部。然而在使用射線相交方法時存在局限性，即可能存在射線恰好經(jīng)過外形尖點(diǎn)的情況，導(dǎo)致判斷的結(jié)果出現(xiàn)錯誤。

針對上述傳統(tǒng)射線相交方法存在的問題，本文作者團(tuán)隊提出了一種新的魯棒的判斷算法，與射線相交方法相配合，可以準(zhǔn)確識別射線經(jīng)過外形拐點(diǎn)的不同情形，稱為奇異性判斷算法[52]。這種方法通過判斷相應(yīng)的線段或面與射線的相對位置從而更加準(zhǔn)確地判斷目標(biāo)網(wǎng)格點(diǎn)的類型。

圖3[52]顯示的是該方法在二維情況中的應(yīng)用，在圖3(a)中，射線R以待判斷的目標(biāo)點(diǎn)T為起始點(diǎn)，經(jīng)過離散線段之間的尖點(diǎn)P，兩條以P為端點(diǎn)的線段分別為PA和PB。PL和PR均在射線的同一側(cè)，因此目標(biāo)點(diǎn)為物面外的點(diǎn)。在圖3(b) 中，尖點(diǎn)處的線段PA和PB在射線R的不同側(cè)，則目標(biāo)點(diǎn)為物面內(nèi)的點(diǎn)。其中，r1為以P為起點(diǎn)，方向與R相同的矢量，r2則為方向相反的矢量。PA對應(yīng)的向量l1(對應(yīng)外法向n1)和射線的反向矢量以及PB對應(yīng)的向量l2(對應(yīng)外法向n2)與射線R的夾角分別為α和β。

圖3 射線經(jīng)過拐點(diǎn)[52]Fig.3 Ray passing inflection point [52]

針對三維外形，可通過切割平面方法建立目標(biāo)射線與相交線的相對位置關(guān)系，從而將三維問題轉(zhuǎn)化為二維問題。以圖4[29]為例，在射線R與切割線PA和PB組成的切割平面內(nèi)，即可按照前面的二維問題處理方式來識別射線與離散三角形的相對位置關(guān)系。其中，nPA和nPB分別為切割線PA和PB對應(yīng)的平面內(nèi)的外法向，R為過目標(biāo)點(diǎn)T的射線。

圖4 通過切割面將三維問題轉(zhuǎn)化為二維問題[29]Fig.4 Conversion of 3D problem to 2D problem by cutting plane passing ray[29]

對于網(wǎng)格在物面內(nèi)外的判斷，Aftosmis等提出了染色算法[58]，以快速確定與物面不相交的網(wǎng)格單元類型。這種算法通過查詢目標(biāo)網(wǎng)格單元相鄰非物面相交單元的類型，直接給予其相同的網(wǎng)格類型，目標(biāo)單元就像被鄰居單元“染色”，故稱為染色算法。肖涵山等[33]提出了一種與之相類似的“推進(jìn)查找”方法，以判斷網(wǎng)格單元是否在物面內(nèi)部。該方法首先將與物面相交的單元刪除，這樣便切斷了外部單元和物面內(nèi)部單元的聯(lián)系，然后將流場外邊界上的單元標(biāo)記為外部單元，將與之相連的也標(biāo)記為外部單元并進(jìn)行推進(jìn)，從而快速準(zhǔn)確地判定出內(nèi)外單元。

本文作者團(tuán)隊基于上述方法思想，結(jié)合FTT數(shù)據(jù)結(jié)構(gòu)，發(fā)展了適用于自適應(yīng)笛卡爾網(wǎng)格的染色算法[52]。如圖5所示，首先在外部和內(nèi)部分別找一種子單元，并標(biāo)記為OutCell和InCell，接下來對其鄰居單元進(jìn)行判斷，只要是不與物面相交的鄰居單元便標(biāo)記為與種子單元相同的網(wǎng)格類型，從而快速確定網(wǎng)格單元是在物面內(nèi)部還是外部。從表1中可以看出，使用改進(jìn)后的染色算法相比于傳統(tǒng)的判別方式，所用的時間大大減少。

表1 染色算法 vs 傳統(tǒng)判別方法

由于空間笛卡爾網(wǎng)格單元都是標(biāo)準(zhǔn)的AABB(Axis-Aligned Bounding Box)[59]，因此另一種笛卡爾網(wǎng)格單元與三角形的相交判斷方法可歸結(jié)為AABB與三角形的相交判定，如圖6所示。該種方法基于獨(dú)立軸理論，此理論可表述為當(dāng)兩個凸多面體A和B滿足以下兩個條件之一則可判定二者不相交：①A和B可以被平行于A或B的任意一個面的法向的軸分離；②A和B可以沿著由A與B的邊叉乘形成的軸分離。在具體實現(xiàn)過程中，需要對13根軸進(jìn)行測試：前3根軸為AABB的面法向軸(e1,e2,e3)，用來判定三角形的AABB與笛卡爾網(wǎng)格的AABB是否有重疊，如果沒有，則不相交；第4根軸為三角形的面法向軸(n)，用來判定AABB是否與三角形所在平面相交；最后9根軸為AABB的3個面法向軸分別與三角形的三條邊叉乘形成的軸。以上13根軸的測試均通過，則表明沒有分離軸存在，即AABB與三角形不相交；否則，只要有分離軸存在，即可結(jié)束測試，AABB與三角形相交。通過這種方法可以對待判斷網(wǎng)格單元進(jìn)行快速篩選，減小一定的計算量。

圖6 AABB與三角形相交判斷示意圖Fig.6 Diagram of intersection judgment between AABB and triangle

1.4 并行加速技術(shù)

笛卡爾網(wǎng)格在復(fù)雜構(gòu)型/流動問題中應(yīng)用時存在網(wǎng)格規(guī)模過大的難題，為了保證網(wǎng)格生成的效率，有必要開展大規(guī)模并行計算技術(shù)的應(yīng)用。事實上，高可擴(kuò)展性、高負(fù)載平衡的并行流場算法和實現(xiàn)仍是目前國際上的一個難點(diǎn)，是CFD模擬中提高計算效率的重要環(huán)節(jié)[60]。若在計算分區(qū)過程中無法做到負(fù)載平衡，則并行效率會急劇下降。在非結(jié)構(gòu)網(wǎng)格和混合網(wǎng)格并行計算方法中廣泛使用的網(wǎng)格分區(qū)軟件如METIS和Zoltan等，為了保持分區(qū)的負(fù)載平衡，每次笛卡爾網(wǎng)格進(jìn)行解自適應(yīng)后都需要重新對網(wǎng)格進(jìn)行分區(qū)，這在處理非定常流動問題時，由于笛卡爾網(wǎng)格需要頻繁地進(jìn)行自適應(yīng)，會在分區(qū)上耗費(fèi)大量的時間。針對以上問題和現(xiàn)狀，當(dāng)下需發(fā)展高效的適用于自適應(yīng)笛卡爾網(wǎng)格的并行策略。

本文作者團(tuán)隊針對自適應(yīng)笛卡爾網(wǎng)格下高性能并行計算技術(shù)的應(yīng)用開展了深入研究。在流場計算時基于空間填充曲線(Z型排序是Z曲線，如圖7所示[61]，逆時針排序是Hilbert曲線)按照排序規(guī)則將所有子節(jié)點(diǎn)連接起來形成一維鏈表型數(shù)據(jù)結(jié)構(gòu)；然后，采用網(wǎng)格分區(qū)軟件METIS的思想，直接對一維鏈表式數(shù)據(jù)進(jìn)行分區(qū)。此外，在流場計算過程中，隨著解自適應(yīng)的進(jìn)行，該方法還便于實現(xiàn)動態(tài)的網(wǎng)格分區(qū)，很好地保障計算負(fù)載平衡。同時，前述網(wǎng)格自適應(yīng)加密時遵循2:1加密準(zhǔn)則，即加密的網(wǎng)格單元和其周圍單元的加密層級最多不超過1，結(jié)合網(wǎng)格單元快速ADT搜索算法，可開發(fā)最小化的分區(qū)交接面信息交換算法，實現(xiàn)分區(qū)網(wǎng)格之間信息的快速傳遞。針對三維等熵渦問題對動態(tài)分區(qū)效果、負(fù)載平衡以及并行效率等方面展開了技術(shù)驗證，如圖8所示?？梢钥闯鲭S著渦的運(yùn)動，網(wǎng)格進(jìn)行了自適應(yīng)并實現(xiàn)了自動分區(qū)，與此同時也較好地保持了物理區(qū)域的連續(xù)性，能夠有效地減少進(jìn)程之間的通訊開銷，進(jìn)而提升計算效率。隨著進(jìn)程數(shù)的增加，并行效率保持情況良好，擁有可行的并行規(guī)模擴(kuò)展前景。

圖7 排序方案和Z型空間填充曲線示意圖[61]Fig.7 Schematic diagram of sorting scheme and Z-type space filling curves[61]

圖8 三維渦運(yùn)動問題網(wǎng)格分區(qū)及計算時間占比Fig.8 Grid partition and calculation time of 3D vortex motion problem

2 網(wǎng)格自適應(yīng)技術(shù)

2.1 各向同性自適應(yīng)

上述工作為構(gòu)建面向任意外形的自適應(yīng)笛卡爾網(wǎng)格生成方法奠定了基礎(chǔ)。接下來根據(jù)幾何和流場特征進(jìn)行網(wǎng)格的自適應(yīng)加密或粗化。

由于初始網(wǎng)格一般尺寸較大，無法準(zhǔn)確識別模型的幾何信息和流場的特征結(jié)構(gòu)信息，需要建立相應(yīng)的加密判據(jù)。對于基于幾何特征信息進(jìn)行的加密，一般稱為幾何自適應(yīng)。其基本思路是：通過射線相交或染色算法判斷出相交網(wǎng)格后，對于與物面相交的單元進(jìn)行標(biāo)記和初步加密一定的次數(shù)。為了保證物面與流場之間的網(wǎng)格過渡均勻，一般對過渡層網(wǎng)格也進(jìn)行加密。之后，通過識別幾何的曲率、厚度等特征信息，對精細(xì)結(jié)構(gòu)進(jìn)行進(jìn)一步加密，以保證模型信息的高保真度。最后，對于網(wǎng)格質(zhì)量進(jìn)行檢測和優(yōu)化，即對于加密不合理、影響叉樹結(jié)構(gòu)平衡性、質(zhì)量不高等區(qū)域進(jìn)行排查和修正，直至達(dá)到需要的高質(zhì)量網(wǎng)格。例如，對需要加密的網(wǎng)格進(jìn)行合適的細(xì)化、對相鄰網(wǎng)格層次差大于1的網(wǎng)格進(jìn)行降層等?；跇?gòu)建的幾何自適應(yīng)方法，針對多種復(fù)雜構(gòu)型開展了算例測試，以導(dǎo)彈和翼身組合體為例，如圖9所示[52]。

圖9 幾何自適應(yīng)網(wǎng)格[52]Fig.9 Geometrical adaptive mesh[52]

笛卡爾網(wǎng)格在流場結(jié)構(gòu)自適應(yīng)方面具有天然優(yōu)勢。通過流場自適應(yīng)，可以提高對流動特征結(jié)構(gòu)的分辨率，降低數(shù)值耗散，提高流動模擬的精度，同時也能夠使網(wǎng)格分布更加合理。為了精細(xì)捕捉流動中出現(xiàn)的渦旋、激波、剪切層等特征結(jié)構(gòu)，需要構(gòu)建適用性好、精度可靠的判據(jù)。為此，本文作者團(tuán)隊發(fā)展了一種綜合判據(jù)，即速度散度和旋度相結(jié)合的判斷方法[52]。該判據(jù)既可以通過速度的散度捕捉激波，又可以基于速度的旋度捕捉渦旋和剪切層等特征，相對于單一判據(jù)具有明顯的優(yōu)勢。針對典型算例測試了構(gòu)建的自適應(yīng)方法，得到了較好的效果。流場自適應(yīng)結(jié)果如圖10 所示。

圖10 超聲速流動自適應(yīng)網(wǎng)格Fig.10 Adaptive mesh based on supersonic flow

2.2 各向異性自適應(yīng)

各向異性自適應(yīng)即笛卡爾網(wǎng)格在加密過程中會存在多種加密方式，從而生成形狀、大小不一的子單元。以二維各向異性為例，一個父單元有4種加密方式，分別為在x、y兩個方向同時加密,在x、y其中一個方向加密以及在x、y兩個方向均不加密。加密示意圖如圖11所示。

圖11 二維各向異性笛卡爾網(wǎng)格加密形式Fig.11 2D anisotropic Cartesian mesh refinement form

在各向異性自適應(yīng)加密過程中，會產(chǎn)生方向、大小不同的相鄰網(wǎng)格，因此鄰居單元查詢算法比各向同性更為復(fù)雜。為便于顯示，以二維為例說明可能出現(xiàn)的鄰居情況，如圖12所示。為了防止相鄰單元層級差過大導(dǎo)致解不穩(wěn)定的情況，本文作者團(tuán)隊采用相鄰單元層級差不大于1的加密/粗化原則，從根本上解決各向異性笛卡爾網(wǎng)格鄰居查詢和網(wǎng)格生成質(zhì)量問題。

圖12 各向異性笛卡爾網(wǎng)格鄰居分布示意圖Fig.12 Schematic diagram of anisotropic Cartesian mesh neighborhood distribution

在網(wǎng)格自適應(yīng)過程中，先依據(jù)給定的網(wǎng)格尺寸構(gòu)造初始粗網(wǎng)格，以粗網(wǎng)格為基礎(chǔ)根據(jù)需要進(jìn)行加密。為了減少加密邊界的復(fù)雜程度，首先將網(wǎng)格單元進(jìn)行一定程度的各向同性加密；然后在物面曲率變化較大的區(qū)域進(jìn)行各向異性加密，若單元對應(yīng)的物面斜率在i方向上分量和臨近單元在相同方向上的斜率分量相差大于設(shè)定閾值，則此單元在i方向上進(jìn)行加密。至于流場解自適應(yīng)判據(jù)，若繼續(xù)采用各向同性笛卡爾網(wǎng)格解自適應(yīng)判據(jù)會出現(xiàn)失效的情況。因此，采用基于速度旋度和散度以及壓力梯度的綜合判據(jù)，以實現(xiàn)不同流場結(jié)構(gòu)(激波、漩渦等)的細(xì)致捕捉。

在進(jìn)行流場自適應(yīng)之前，首先需要解決插值模板點(diǎn)的問題。由于各向異性網(wǎng)格鄰居排列組合形式多樣，為了提高結(jié)果精度，插值模板點(diǎn)的選取需要對特定情況進(jìn)行特殊處理。以圖13為例，當(dāng)網(wǎng)格i為Type=1的Case 1情況時，i+1模板點(diǎn)的流場值記為其所在單元的值，而Case 2、Case 3以及Case4的i+1模板點(diǎn)的流場值需要周圍網(wǎng)格單元插值得到。圖14為球柱流場解自適應(yīng)結(jié)果，圖例中的ρ代表無量綱的密度分布。從圖中可以看出，各向異性笛卡爾網(wǎng)格可較好捕捉激波和尾流渦結(jié)構(gòu)，同時，在激波附近，網(wǎng)格分布合理，過渡均勻。網(wǎng)格數(shù)目較各向同性笛卡爾網(wǎng)格可降低15%～30%，內(nèi)存占用下降10%～25%，網(wǎng)格生成效率提高20%～30%。

圖13 部分插值模板點(diǎn)選取示意圖Fig.13 Schematic diagram of partial interpolation template points selection

圖14 各向異性流場解自適應(yīng)結(jié)果Fig.14 Adaptive results of anisotropic flow field solutions

2.3 法向射線加密

針對傳統(tǒng)自適應(yīng)笛卡爾網(wǎng)格方法在面對高雷諾數(shù)黏性流動問題時，存在網(wǎng)格數(shù)目過大的問題，Ruffin等[44]借鑒貼體結(jié)構(gòu)網(wǎng)格大縱橫比的特點(diǎn)，通過在壁面配置射線來確定加密區(qū)域，射線間則用粗網(wǎng)格過渡，構(gòu)建了一種擬結(jié)構(gòu)化的網(wǎng)格生成方法，即法向射線加密技術(shù)，提供了一種在可接受的網(wǎng)格數(shù)量下模擬高雷諾數(shù)黏性流動的笛卡爾方法。實際上，在邊界層法向上布置高分辨率網(wǎng)格相對于沿壁面方向布置同等密度的網(wǎng)格會擁有更高的計算收益。法向射線加密技術(shù)則是基于該思想，在法向射線上布置足夠細(xì)密的網(wǎng)格單元，在沿壁面方向就可以使用比較粗的網(wǎng)格。本文作者團(tuán)隊對于該方法進(jìn)行了發(fā)展和實現(xiàn)，如圖15所示，采用法向射線技術(shù)相應(yīng)的網(wǎng)格數(shù)目比傳統(tǒng)在邊界層均勻加密的笛卡爾網(wǎng)格數(shù)目少很多，其網(wǎng)格總體分布與與貼體結(jié)構(gòu)網(wǎng)格類似，利于在可接受的網(wǎng)格規(guī)模下捕捉邊界層特征。

圖15 自適應(yīng)下笛卡爾網(wǎng)格和法向射線比較Fig.15 Adaptive Cartesian grid and normal ray comparison

法向射線加密技術(shù)分為3部分，分別是射線區(qū)域劃分、射線間信息交互、計算插值與修正。射線區(qū)域的劃分如圖16所示，將整個計算域分為4個區(qū)域。4個區(qū)域中，藍(lán)色和灰色區(qū)域并不參與流場求解。藍(lán)色網(wǎng)格的流場值通過射線插值得到，通過應(yīng)用射線間信息交互技術(shù)，構(gòu)建了基于指針方法的相鄰射線數(shù)據(jù)交互方法，并根據(jù)這些流場信息可靠度高的射線單元進(jìn)行流場計算，射線間的粗網(wǎng)格的流場信息基于插值方法獲得，實現(xiàn)了在較少法向射線的情況下高精度模擬邊界層流動，大幅度降低了網(wǎng)格規(guī)模。

圖16 法向射線加密技術(shù)的計算區(qū)域劃分Fig.16 Region division of normal ray refinement technology

3 黏性物面邊界處理

由于笛卡爾網(wǎng)格的非貼體特性，與物面相交會產(chǎn)生大小不同、形狀各異的切割單元。與傳統(tǒng)的貼體結(jié)構(gòu)網(wǎng)格相比，這種單元形式不利于邊界層流動特征的捕捉，容易引起邊界層流動數(shù)據(jù)的非物理振蕩。為了減小這種數(shù)值誤差，若通過不計成本的加密方式，還會帶來網(wǎng)格規(guī)模過大的問題。因此，發(fā)展了兩種處理方法。

一種方法是采用重疊網(wǎng)格技術(shù)[62]，將貼體結(jié)構(gòu)網(wǎng)格和笛卡爾網(wǎng)格的優(yōu)勢相結(jié)合，在物面邊界附近提前生成貼體的結(jié)構(gòu)或非結(jié)構(gòu)網(wǎng)格，而在流場其他區(qū)域生成笛卡爾網(wǎng)格。通過兩套網(wǎng)格的銜接，將笛卡爾網(wǎng)格的非貼體物面邊界問題轉(zhuǎn)移到兩套網(wǎng)格的交界面處。

另一種方法是在笛卡爾網(wǎng)格框架下引入通過邊界層近似理論得到的壁面函數(shù)方法，以降低數(shù)值方法對網(wǎng)格尺度的依賴、減小邊界層內(nèi)計算網(wǎng)格數(shù)量。

3.1 重疊網(wǎng)格技術(shù)

基于重疊網(wǎng)格開展數(shù)值計算需要網(wǎng)格間的數(shù)據(jù)傳遞和交互，在由粗網(wǎng)格貢獻(xiàn)單元向細(xì)網(wǎng)格插值(或者雙方長細(xì)比相差很大)時，都會引入較大數(shù)值誤差。因此貢獻(xiàn)單元與插值點(diǎn)所在網(wǎng)格單元應(yīng)盡量保持一致的形狀和尺度?；诘芽柧W(wǎng)格的自適應(yīng)優(yōu)勢，本文作者團(tuán)隊發(fā)展了重疊區(qū)尺度匹配的笛卡爾網(wǎng)格自適應(yīng)技術(shù)，如圖17所示。在初始笛卡爾網(wǎng)格的基礎(chǔ)上，依據(jù)背景網(wǎng)格與貼體網(wǎng)格的尺寸特征進(jìn)行自適應(yīng)加密處理。通過自適應(yīng)匹配，保證邊界附近網(wǎng)格尺度的一致性，使流場在混合網(wǎng)格邊界處光滑過渡，避免由于背景直角坐標(biāo)網(wǎng)格與貼體結(jié)構(gòu)網(wǎng)格在邊界處的網(wǎng)格尺度相差太大而引起數(shù)值耗散。

圖17 重疊網(wǎng)格交界面自適應(yīng)笛卡爾網(wǎng)格Fig.17 Adaptive Cartesian mesh at hybrid mesh interfaces

貢獻(xiàn)單元搜索[63]即“尋點(diǎn)”，是確定網(wǎng)格點(diǎn)在其他網(wǎng)格中的位置、查找網(wǎng)格貢獻(xiàn)單元的統(tǒng)稱。在網(wǎng)格的生成過程中需要進(jìn)行大量的“尋點(diǎn)”操作，準(zhǔn)確、高效的貢獻(xiàn)單元搜索方法是重疊網(wǎng)格的關(guān)鍵技術(shù)之一。對笛卡爾背景網(wǎng)格和近壁貼體網(wǎng)格的重疊邊界貢獻(xiàn)單元搜索過程中，可結(jié)合笛卡爾網(wǎng)格正交規(guī)則的幾何特征這一優(yōu)勢，實現(xiàn)搜索的高效、準(zhǔn)確以及存儲的低需求。笛卡爾網(wǎng)格貢獻(xiàn)單元搜索如圖18所示，以目標(biāo)笛卡爾網(wǎng)格格心為中心，以網(wǎng)格的空間特征尺寸為參考，構(gòu)建空間大小合適的包圍盒，其中P1和P2表示包圍盒兩對角頂點(diǎn)；以空間坐標(biāo)為衡量尺度，搜索該包圍盒內(nèi)所含有的結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)；統(tǒng)計包圍盒內(nèi)結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)數(shù)目，計算與目標(biāo)笛卡爾網(wǎng)格格心的距離。若網(wǎng)格節(jié)點(diǎn)數(shù)滿足貢獻(xiàn)單元所需數(shù)目，依據(jù)距離參數(shù)選取所需結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)；若網(wǎng)格節(jié)點(diǎn)數(shù)少于貢獻(xiàn)單元所需數(shù)目，則繼續(xù)擴(kuò)大包圍盒范圍，直到可選貢獻(xiàn)單元數(shù)目滿足條件為止。以包圍盒中所選結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)為貢獻(xiàn)單元，結(jié)合相應(yīng)插值算法，獲得笛卡爾網(wǎng)格流動參數(shù)。

圖18 笛卡爾網(wǎng)格貢獻(xiàn)單元搜索Fig.18 Cartesian grids donor cell search

通過上述發(fā)展的重疊區(qū)網(wǎng)格匹配技術(shù)和數(shù)據(jù)交互技術(shù)，與附面層貼體結(jié)構(gòu)網(wǎng)格生成方法相結(jié)合[64]，構(gòu)建了基于近壁貼體結(jié)構(gòu)網(wǎng)格/笛卡爾網(wǎng)格的重疊網(wǎng)格技術(shù)。同時開展了真實管道列車的算例驗證工作，初步考核了所發(fā)展的重疊笛卡爾網(wǎng)格方法的有效性和可靠性(圖19)。

圖19 光順外形管道列車Fig.19 Smooth shape pipe train

3.2 壁面函數(shù)方法

傳統(tǒng)的湍流壁面函數(shù)基于貼體網(wǎng)格提出，求解過程中直接應(yīng)用壁面邊界條件，這在非貼體網(wǎng)格下難以滿足。非貼體網(wǎng)格的單元與物面的不一致性，將傳統(tǒng)壁面函數(shù)直接推廣至非貼體網(wǎng)格存在一定困難。

2002年2月，中國礦業(yè)大學(xué)北京東校區(qū)（其圖書館以下簡稱“北館”）并入中國傳媒大學(xué)。經(jīng)過多年建設(shè)與深耕，中國傳媒大學(xué)以新聞傳播學(xué)、戲劇與影視學(xué)、信息與通信工程為龍頭，逐漸形成文、工、管、經(jīng)、法、理多學(xué)科協(xié)調(diào)發(fā)展的新聞傳媒特色高校，為傳媒行業(yè)的臺前幕后培養(yǎng)專門人才，所有專業(yè)人才培養(yǎng)與課程設(shè)計都與傳媒行業(yè)緊密結(jié)合。由于兩館原來服務(wù)的讀者對象不同，藏書重點(diǎn)和收藏范圍也不一樣。筆者所在圖書館在這方面所做的規(guī)劃和實踐過程如下：

本文作者團(tuán)隊基于Spalding[65]提出的壁面函數(shù)模型，對該模型在笛卡爾網(wǎng)格下進(jìn)行適用性修正。該模型格式較簡單，且經(jīng)過很多學(xué)者的考核驗證，具有較高的可信度。針對可壓縮湍流黏性流動問題的計算，為解決邊界層非線性結(jié)構(gòu)特征，引入虛擬點(diǎn)相對物面的參考點(diǎn)來取代無黏情況下的對稱點(diǎn)。同時基于壁面函數(shù)修正湍流壁面邊界條件，建立了面向笛卡爾網(wǎng)格的壁面函數(shù)方法，并開展了圓柱(圖20，Cp為圓柱表面壓力系數(shù)分布，θ為沿圓柱周向角度值)、翼型(圖21)等典型二維算例的考核驗證工作。通過計算結(jié)果對比，證明了采用壁面函數(shù)方法可以有效降低壁面網(wǎng)格尺度要求，并在網(wǎng)格較粗時達(dá)到良好的模擬精度。

圖20 二維圓柱算例Fig.20 Two dimensional cylindrical calculation case

圖21 NACA0012翼型算例Fig.21 NACA0012 airfoil calculation case

胡偶在其前期的工作基礎(chǔ)上，繼續(xù)對基于SSTk-ω湍流模型的Spalding壁面函數(shù)模型進(jìn)行發(fā)展[23]。為了使參考點(diǎn)在壁面附近均勻分布以及提高壁面附近的湍流模擬精度和數(shù)值插值精度，提出了一種探針型浸入邊界法(Ghost Cell Probe Immersed Boundary Method), 如圖22所示，再結(jié)合針對SA和SST湍流模型的新型常微分化(ODE)壁面模型定義黏性壁面條件，形成耦合虛擬單元-壁面模型的湍流壁面邊界條件重構(gòu)方法。此外，除常用的Spalding等提出的壁面函數(shù)外，在現(xiàn)有浸入邊界-虛擬單元方法體系下，耦合了Knopp[66]提出的自適應(yīng)壁面函數(shù)模型(Adaptive Wall Function)，能夠根據(jù)不同湍流模型自適應(yīng)地調(diào)整經(jīng)驗公式，擴(kuò)展了壁面函數(shù)模型的適用范圍。針對三維問題，基于自適應(yīng)混合笛卡爾網(wǎng)格體系進(jìn)行了壁面函數(shù)理論研究，并通過典型的ONERA M6翼型繞流問題初步驗證了理論方法，如圖23所示，為后續(xù)移植到三維自適應(yīng)笛卡爾網(wǎng)格體系做了技術(shù)積累。

圖22 探針型浸入邊界法示意圖Fig.22 Schematic diagram of probe immersed boundary method

圖23 ONERA M6 網(wǎng)格及繞流模擬結(jié)果Fig.23 ONERA M6 grid and flow simulation results

目前，針對非貼體網(wǎng)格下的邊界條件的主要處理辦法是，通過線性重構(gòu)或插值方式獲得壁面邊界的流動變量。然而在湍流邊界層內(nèi)，不合理的線性重構(gòu)得到的邊界條件難以滿足湍流情況下的邊界層分布特性，容易誘發(fā)非物理的數(shù)值求解振蕩。另一方面，常見的基于貼體網(wǎng)格發(fā)展的湍流模型(如SA、SST等)要求積分到壁面，對于湍流邊界層內(nèi)的網(wǎng)格具有較高的精度要求，并且在壁面邊界及近壁區(qū)存在相匹配的近壁處理手段(Near-wall Treatment)，例如一方程SA模型使用壁面距離，k-ω模型使用Damping Function和Yap修正項等。在非貼體笛卡爾網(wǎng)格條件下，這也帶來了匹配性問題，即在笛卡爾網(wǎng)格下特殊構(gòu)造的壁面函數(shù)在參數(shù)形態(tài)上難以無縫對接原有成熟的湍流模型。這種匹配問題會帶來計算誤差，污染邊界層求解精度。因此非貼體笛卡爾網(wǎng)格下魯棒、高精度的湍流壁面函數(shù)是關(guān)鍵問題，有待進(jìn)一步深入研究。

本文作者團(tuán)隊發(fā)展了一種適用于邊界層全域的壁面函數(shù)。針對T3B平板邊界層算例，如圖24所示，其中New_wf為基于新型壁面函數(shù)方法得>到的計算結(jié)果，F(xiàn)ull_grd為全網(wǎng)格下的摩阻計算結(jié)果，Exp代表實驗參照結(jié)果?？梢钥闯?，基于SST模型，第1層網(wǎng)格高度在y+=50處；預(yù)測的摩阻系數(shù)Cf與y+(1)=0.01的全網(wǎng)格摩阻相當(dāng)，具有較高的精度。后續(xù)將開展其在笛卡爾網(wǎng)格框架下的適用性應(yīng)用研究。

圖24 沿平板邊界層流向摩擦阻力分布Fig.24 Friction drag distribution along flow direction of plate boundary layer

4 結(jié) 論

本文總結(jié)了國家數(shù)值風(fēng)洞工程黏性自適應(yīng)笛卡爾網(wǎng)格方法的研究進(jìn)展，涵蓋了網(wǎng)格生成、網(wǎng)格自適應(yīng)方法和物面邊界處理技術(shù)等內(nèi)容，涉及南京航空航天大學(xué)、中南大學(xué)、中國空氣動力研究與發(fā)展中心等單位在笛卡爾網(wǎng)格技術(shù)優(yōu)化、非貼體笛卡爾網(wǎng)格物面邊界處理等方面深入開展的工作，主要得出了以下結(jié)論：

1) 網(wǎng)格生成技術(shù)方面，分別在網(wǎng)格數(shù)據(jù)結(jié)構(gòu)、物面單元檢索技術(shù)、網(wǎng)格單元類型判斷方法、高性能并行計算技術(shù)等開展了應(yīng)用或改進(jìn)工作，提高了三維自適應(yīng)笛卡爾網(wǎng)格生成的效率和魯棒性，為開展大規(guī)模網(wǎng)格量下復(fù)雜構(gòu)型黏性流動問題的模擬研究奠定了基礎(chǔ)。

2) 網(wǎng)格自適應(yīng)技術(shù)方面，分別對各向同性自適應(yīng)方法、各向異性自適應(yīng)方法、法向射線加密方法進(jìn)行了實現(xiàn)，構(gòu)建了高效可靠的自適應(yīng)判據(jù)，可優(yōu)化網(wǎng)格分布，并有效降低網(wǎng)格規(guī)模，為復(fù)雜外形/流動問題的高效、高精度模擬豐富了技術(shù)手段。

3) 物面邊界處理方面，分別對貼體式的重疊笛卡爾網(wǎng)格方法、非貼體式的浸入邊界方法進(jìn)行了研究與發(fā)展，提供了黏性物面邊界模擬的不同技術(shù)手段，通過壁面函數(shù)方法的引入和適用性改進(jìn)，初步形成了可有效降低湍流模擬的網(wǎng)格規(guī)模、保證非貼體黏性物面模擬精度的技術(shù)方法，為基于自適應(yīng)笛卡爾網(wǎng)格方法開展高雷諾數(shù)黏性流動高效、高精度的模擬提供了有力的技術(shù)支撐。

本文對黏性自適應(yīng)笛卡爾網(wǎng)格方法的國內(nèi)外研究進(jìn)展進(jìn)行了簡要梳理，重點(diǎn)對國家數(shù)值風(fēng)洞工程黏性笛卡爾網(wǎng)格技術(shù)課題進(jìn)展進(jìn)行了階段性總結(jié)介紹。在后續(xù)工作中，將在國家數(shù)值風(fēng)洞工程的支撐下，直面黏性笛卡爾網(wǎng)格技術(shù)的難點(diǎn)和痛點(diǎn)，繼續(xù)開展笛卡爾網(wǎng)格方法基礎(chǔ)理論的深化研究和工程推廣應(yīng)用，為國家數(shù)值風(fēng)洞工程豐富復(fù)雜構(gòu)型和復(fù)雜流動模擬能力提供技術(shù)儲備。