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

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

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

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

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

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

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

相對(duì)于傳統(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)格生成自動(dòng)化程度高、復(fù)雜外形適應(yīng)性好、非定常/多尺度等流動(dòng)結(jié)構(gòu)捕捉能力強(qiáng)等優(yōu)勢，適用于湍流等復(fù)雜流動(dòng)問題的仿真模擬。因此自適應(yīng)笛卡爾網(wǎng)格方法具有重要的研究意義和廣闊的應(yīng)用前景。

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

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

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

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

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

其中，混合笛卡爾網(wǎng)技術(shù)相對(duì)使用較多，將貼體結(jié)構(gòu)網(wǎng)格和笛卡爾網(wǎng)格的優(yōu)勢相結(jié)合，通過在近壁使用貼體網(wǎng)格來避免笛卡爾網(wǎng)格邊界處理的困難。網(wǎng)格生成與空氣動(dò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ā)展了多個(gè)基于笛卡爾/結(jié)構(gòu)網(wǎng)格混合的知名網(wǎng)格技術(shù)求解器，例如，NASA著名的OVERFLOW[30]求解器就是使用笛卡爾背景網(wǎng)格耦合貼體結(jié)構(gòu)網(wǎng)格的形式。此外，由美國國防部支持開發(fā)的旋翼計(jì)算平臺(tái)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ù)對(duì)Wang的“黏性自適應(yīng)笛卡爾方法”進(jìn)行改進(jìn)，使之能夠處理包括縫隙、凹槽、凸起等結(jié)構(gòu)的更為復(fù)雜的幾何外形流動(dòng)。

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

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

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

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

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

國外研究人員將近壁模型和壁面函數(shù)方法在笛卡爾網(wǎng)格框架下進(jìn)行應(yīng)用，定義湍流壁面邊界條件，從而降低邊界層內(nèi)計(jì)算網(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ù)流動(dòng)問題。Capizzano[48]則采用雙層壁面模型修正浸入邊界網(wǎng)格交接面上湍流變量值的方法，進(jìn)行高雷諾數(shù)可壓縮流動(dòng)問題的模擬研究。Bernardini等[49]基于浸入邊界方法，通過壁面模型耦合SA模型，用于延遲脫體渦模擬(Delayed Detached Eddy Simulation, DDES)中，對(duì)于高雷諾數(shù)下壁面分離的湍流流動(dòng)開展了模擬研究，獲得了高精度、可靠的數(shù)據(jù)結(jié)果。這些學(xué)者的研究工作對(duì)基于笛卡爾網(wǎng)格技術(shù)進(jìn)行湍流問題的模擬起到了極大的推進(jìn)作用。

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

雖然國內(nèi)外學(xué)者對(duì)于黏性笛卡爾網(wǎng)格技術(shù)已經(jīng)開展了部分研究工作，促進(jìn)了笛卡爾網(wǎng)格方法在復(fù)雜黏性流動(dòng)問題中的應(yīng)用。但是由于自適應(yīng)笛卡爾網(wǎng)格技術(shù)模擬黏性流動(dòng)難度較大，且對(duì)于計(jì)算資源要求較高，國內(nèi)外對(duì)于黏性自適應(yīng)笛卡爾網(wǎng)格方法的研究相對(duì)其他網(wǎng)格技術(shù)而言仍然偏少，尤其是國內(nèi)，尚處于起步階段，公開的學(xué)術(shù)成果極少。因此，目前對(duì)于黏性自適應(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)實(shí)現(xiàn)黏性笛卡爾網(wǎng)格技術(shù)的工程實(shí)用化。本文聚焦國家數(shù)值風(fēng)洞工程，總結(jié)介紹了黏性笛卡爾網(wǎng)格方法相關(guān)課題的研究進(jìn)展，重點(diǎn)從笛卡爾網(wǎng)格高效生成技術(shù)、笛卡爾網(wǎng)格自適應(yīng)技術(shù)以及黏性物面邊界處理等3個(gè)方面展開。

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

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

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

為了解決上述傳統(tǒng)叉樹數(shù)據(jù)結(jié)構(gòu)的問題，本文作者團(tuán)隊(duì)在笛卡爾網(wǎng)格存儲(chǔ)時(shí)采用了叉樹數(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)建快速訪問鄰居單元的“線程”，實(shí)現(xiàn)在不同檢索操作中數(shù)據(jù)信息的快速訪問，提供了高效、易并行的叉樹信息訪問路徑，同時(shí)節(jié)省數(shù)據(jù)節(jié)點(diǎn)，降低叉樹信息存儲(chǔ)的內(nèi)存占用,大大提高計(jì)算效率。

圖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%的時(shí)間。在內(nèi)存占用方面，采用傳統(tǒng)的叉樹數(shù)據(jù)結(jié)構(gòu)，要存儲(chǔ)完上述信息，每個(gè)單元需要17 words的內(nèi)存，而如果采用FTT方式，每個(gè)單元需要2 words的內(nèi)存，相對(duì)于傳統(tǒng)方式明顯降低[52]。

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

在生成網(wǎng)格的過程中，網(wǎng)格類型的判斷占據(jù)了大部分時(shí)間。為了高效判斷網(wǎng)格類型，Bonet和Peraire[54]采用ADT(Alternating Digital Tree)數(shù)據(jù)結(jié)構(gòu)來存儲(chǔ)物面網(wǎng)格，以縮小可能相交的物面單元檢測范圍，從而提高相交判斷的效率。在該方法中，每個(gè)ADT節(jié)點(diǎn)都可作為父節(jié)點(diǎn)或子節(jié)點(diǎn)對(duì)應(yīng)一個(gè)笛卡爾空間，父節(jié)點(diǎn)包含其所對(duì)應(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)隊(duì)采用其改進(jìn)版本KDT(K-Dimensional Tree)結(jié)構(gòu)，其更容易構(gòu)造平衡優(yōu)質(zhì)的二叉樹[56-57]。

KDT是一種多維的二叉樹結(jié)構(gòu)，KD樹的每一層都和某個(gè)維度相關(guān)，每個(gè)節(jié)點(diǎn)都是一個(gè)k維節(jié)點(diǎn)。KDT的叉樹在構(gòu)建時(shí)，只需要將當(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]，對(duì)于圖2(a) 中所示的物面網(wǎng)格，KDT叉樹數(shù)據(jù)結(jié)構(gòu)相比于ADT叉樹數(shù)據(jù)結(jié)構(gòu)更加平衡，其相應(yīng)的操作復(fù)雜度也明顯降低。對(duì)于物面網(wǎng)格較少的情況，宜采用ADT來存儲(chǔ)物面單元；當(dāng)物面網(wǎng)格量較多時(shí)，則宜采用KDT。這是因?yàn)?，KDT每進(jìn)行一次空間劃分，均需要進(jìn)行一次排序操作，當(dāng)網(wǎng)格量較大時(shí)，樹的深度將會(huì)比ADT小很多，此時(shí)這種排序的優(yōu)勢才會(huì)凸顯。

圖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)格相對(duì)于物面的位置關(guān)系，可將網(wǎng)格劃分為外部網(wǎng)格、內(nèi)部網(wǎng)格以及相交網(wǎng)格單元。實(shí)踐表明，相交判斷占據(jù)了笛卡爾網(wǎng)格生成的大量時(shí)間，并且是網(wǎng)格生成成敗的關(guān)鍵因素。

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

針對(duì)上述傳統(tǒng)射線相交方法存在的問題，本文作者團(tuán)隊(duì)提出了一種新的魯棒的判斷算法，與射線相交方法相配合，可以準(zhǔn)確識(shí)別射線經(jīng)過外形拐點(diǎn)的不同情形，稱為奇異性判斷算法[52]。這種方法通過判斷相應(yīng)的線段或面與射線的相對(duì)位置從而更加準(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對(duì)應(yīng)的向量l1(對(duì)應(yīng)外法向n1)和射線的反向矢量以及PB對(duì)應(yīng)的向量l2(對(duì)應(yīng)外法向n2)與射線R的夾角分別為α和β。

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

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

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

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

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

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

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

1.4 并行加速技術(shù)

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

本文作者團(tuán)隊(duì)針對(duì)自適應(yīng)笛卡爾網(wǎng)格下高性能并行計(jì)算技術(shù)的應(yīng)用開展了深入研究。在流場計(jì)算時(shí)基于空間填充曲線(Z型排序是Z曲線，如圖7所示[61]，逆時(shí)針排序是Hilbert曲線)按照排序規(guī)則將所有子節(jié)點(diǎn)連接起來形成一維鏈表型數(shù)據(jù)結(jié)構(gòu)；然后，采用網(wǎng)格分區(qū)軟件METIS的思想，直接對(duì)一維鏈表式數(shù)據(jù)進(jìn)行分區(qū)。此外，在流場計(jì)算過程中，隨著解自適應(yīng)的進(jìn)行，該方法還便于實(shí)現(xiàn)動(dòng)態(tài)的網(wǎng)格分區(qū)，很好地保障計(jì)算負(fù)載平衡。同時(shí)，前述網(wǎng)格自適應(yīng)加密時(shí)遵循2:1加密準(zhǔn)則，即加密的網(wǎng)格單元和其周圍單元的加密層級(jí)最多不超過1，結(jié)合網(wǎng)格單元快速ADT搜索算法，可開發(fā)最小化的分區(qū)交接面信息交換算法，實(shí)現(xiàn)分區(qū)網(wǎng)格之間信息的快速傳遞。針對(duì)三維等熵渦問題對(duì)動(dòng)態(tài)分區(qū)效果、負(fù)載平衡以及并行效率等方面展開了技術(shù)驗(yàn)證，如圖8所示?？梢钥闯鲭S著渦的運(yùn)動(dòng)，網(wǎng)格進(jìn)行了自適應(yīng)并實(shí)現(xiàn)了自動(dòng)分區(qū)，與此同時(shí)也較好地保持了物理區(qū)域的連續(xù)性，能夠有效地減少進(jìn)程之間的通訊開銷，進(jìn)而提升計(jì)算效率。隨著進(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)動(dòng)問題網(wǎng)格分區(qū)及計(jì)算時(shí)間占比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)確識(shí)別模型的幾何信息和流場的特征結(jié)構(gòu)信息，需要建立相應(yīng)的加密判據(jù)。對(duì)于基于幾何特征信息進(jìn)行的加密，一般稱為幾何自適應(yīng)。其基本思路是：通過射線相交或染色算法判斷出相交網(wǎng)格后，對(duì)于與物面相交的單元進(jìn)行標(biāo)記和初步加密一定的次數(shù)。為了保證物面與流場之間的網(wǎng)格過渡均勻，一般對(duì)過渡層網(wǎng)格也進(jìn)行加密。之后，通過識(shí)別幾何的曲率、厚度等特征信息，對(duì)精細(xì)結(jié)構(gòu)進(jìn)行進(jìn)一步加密，以保證模型信息的高保真度。最后，對(duì)于網(wǎng)格質(zhì)量進(jìn)行檢測和優(yōu)化，即對(duì)于加密不合理、影響叉樹結(jié)構(gòu)平衡性、質(zhì)量不高等區(qū)域進(jìn)行排查和修正，直至達(dá)到需要的高質(zhì)量網(wǎng)格。例如，對(duì)需要加密的網(wǎng)格進(jìn)行合適的細(xì)化、對(duì)相鄰網(wǎng)格層次差大于1的網(wǎng)格進(jìn)行降層等?；跇?gòu)建的幾何自適應(yīng)方法，針對(duì)多種復(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)，可以提高對(duì)流動(dòng)特征結(jié)構(gòu)的分辨率，降低數(shù)值耗散，提高流動(dòng)模擬的精度，同時(shí)也能夠使網(wǎng)格分布更加合理。為了精細(xì)捕捉流動(dòng)中出現(xiàn)的渦旋、激波、剪切層等特征結(jié)構(gòu)，需要構(gòu)建適用性好、精度可靠的判據(jù)。為此，本文作者團(tuán)隊(duì)發(fā)展了一種綜合判據(jù)，即速度散度和旋度相結(jié)合的判斷方法[52]。該判據(jù)既可以通過速度的散度捕捉激波，又可以基于速度的旋度捕捉渦旋和剪切層等特征，相對(duì)于單一判據(jù)具有明顯的優(yōu)勢。針對(duì)典型算例測試了構(gòu)建的自適應(yīng)方法，得到了較好的效果。流場自適應(yīng)結(jié)果如圖10 所示。

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

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

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

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

在各向異性自適應(yīng)加密過程中，會(huì)產(chǎn)生方向、大小不同的相鄰網(wǎng)格，因此鄰居單元查詢算法比各向同性更為復(fù)雜。為便于顯示，以二維為例說明可能出現(xiàn)的鄰居情況，如圖12所示。為了防止相鄰單元層級(jí)差過大導(dǎo)致解不穩(wěn)定的情況，本文作者團(tuán)隊(duì)采用相鄰單元層級(jí)差不大于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)行各向異性加密，若單元對(duì)應(yīng)的物面斜率在i方向上分量和臨近單元在相同方向上的斜率分量相差大于設(shè)定閾值，則此單元在i方向上進(jìn)行加密。至于流場解自適應(yīng)判據(jù)，若繼續(xù)采用各向同性笛卡爾網(wǎng)格解自適應(yīng)判據(jù)會(huì)出現(xiàn)失效的情況。因此，采用基于速度旋度和散度以及壓力梯度的綜合判據(jù)，以實(shí)現(xiàn)不同流場結(jié)構(gòu)(激波、漩渦等)的細(xì)致捕捉。

在進(jìn)行流場自適應(yīng)之前，首先需要解決插值模板點(diǎn)的問題。由于各向異性網(wǎng)格鄰居排列組合形式多樣，為了提高結(jié)果精度，插值模板點(diǎn)的選取需要對(duì)特定情況進(jìn)行特殊處理。以圖13為例，當(dāng)網(wǎng)格i為Type=1的Case 1情況時(shí)，i+1模板點(diǎn)的流場值記為其所在單元的值，而Case 2、Case 3以及Case4的i+1模板點(diǎn)的流場值需要周圍網(wǎng)格單元插值得到。圖14為球柱流場解自適應(yīng)結(jié)果，圖例中的ρ代表無量綱的密度分布。從圖中可以看出，各向異性笛卡爾網(wǎng)格可較好捕捉激波和尾流渦結(jié)構(gòu)，同時(shí)，在激波附近，網(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 法向射線加密

針對(duì)傳統(tǒng)自適應(yīng)笛卡爾網(wǎng)格方法在面對(duì)高雷諾數(shù)黏性流動(dò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ù)黏性流動(dòng)的笛卡爾方法。實(shí)際上，在邊界層法向上布置高分辨率網(wǎng)格相對(duì)于沿壁面方向布置同等密度的網(wǎng)格會(huì)擁有更高的計(jì)算收益。法向射線加密技術(shù)則是基于該思想，在法向射線上布置足夠細(xì)密的網(wǎng)格單元，在沿壁面方向就可以使用比較粗的網(wǎng)格。本文作者團(tuán)隊(duì)對(duì)于該方法進(jìn)行了發(fā)展和實(shí)現(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ū)域劃分、射線間信息交互、計(jì)算插值與修正。射線區(qū)域的劃分如圖16所示，將整個(gè)計(jì)算域分為4個(gè)區(qū)域。4個(gè)區(qū)域中，藍(lán)色和灰色區(qū)域并不參與流場求解。藍(lán)色網(wǎng)格的流場值通過射線插值得到，通過應(yīng)用射線間信息交互技術(shù)，構(gòu)建了基于指針方法的相鄰射線數(shù)據(jù)交互方法，并根據(jù)這些流場信息可靠度高的射線單元進(jìn)行流場計(jì)算，射線間的粗網(wǎng)格的流場信息基于插值方法獲得，實(shí)現(xiàn)了在較少法向射線的情況下高精度模擬邊界層流動(dòng)，大幅度降低了網(wǎng)格規(guī)模。

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

3 黏性物面邊界處理

由于笛卡爾網(wǎng)格的非貼體特性，與物面相交會(huì)產(chǎn)生大小不同、形狀各異的切割單元。與傳統(tǒng)的貼體結(jié)構(gòu)網(wǎng)格相比，這種單元形式不利于邊界層流動(dòng)特征的捕捉，容易引起邊界層流動(dòng)數(shù)據(jù)的非物理振蕩。為了減小這種數(shù)值誤差，若通過不計(jì)成本的加密方式，還會(huì)帶來網(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ù)值方法對(duì)網(wǎng)格尺度的依賴、減小邊界層內(nèi)計(jì)算網(wǎng)格數(shù)量。

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

基于重疊網(wǎng)格開展數(shù)值計(jì)算需要網(wǎng)格間的數(shù)據(jù)傳遞和交互，在由粗網(wǎng)格貢獻(xiàn)單元向細(xì)網(wǎng)格插值(或者雙方長細(xì)比相差很大)時(shí)，都會(huì)引入較大數(shù)值誤差。因此貢獻(xiàn)單元與插值點(diǎn)所在網(wǎng)格單元應(yīng)盡量保持一致的形狀和尺度?；诘芽柧W(wǎng)格的自適應(yīng)優(yōu)勢，本文作者團(tuán)隊(duì)發(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ù)之一。對(duì)笛卡爾背景網(wǎng)格和近壁貼體網(wǎng)格的重疊邊界貢獻(xiàn)單元搜索過程中，可結(jié)合笛卡爾網(wǎng)格正交規(guī)則的幾何特征這一優(yōu)勢，實(shí)現(xiàn)搜索的高效、準(zhǔn)確以及存儲(chǔ)的低需求。笛卡爾網(wǎng)格貢獻(xiàn)單元搜索如圖18所示，以目標(biāo)笛卡爾網(wǎng)格格心為中心，以網(wǎng)格的空間特征尺寸為參考，構(gòu)建空間大小合適的包圍盒，其中P1和P2表示包圍盒兩對(duì)角頂點(diǎn)；以空間坐標(biāo)為衡量尺度，搜索該包圍盒內(nèi)所含有的結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)；統(tǒng)計(jì)包圍盒內(nèi)結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)數(shù)目，計(jì)算與目標(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)格流動(dò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ù)。同時(shí)開展了真實(shí)管道列車的算例驗(yàn)證工作，初步考核了所發(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è)的臺(tái)前幕后培養(yǎng)專門人才，所有專業(yè)人才培養(yǎng)與課程設(shè)計(jì)都與傳媒行業(yè)緊密結(jié)合。由于兩館原來服務(wù)的讀者對(duì)象不同，藏書重點(diǎn)和收藏范圍也不一樣。筆者所在圖書館在這方面所做的規(guī)劃和實(shí)踐過程如下：

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

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

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

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

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

本文作者團(tuán)隊(duì)發(fā)展了一種適用于邊界層全域的壁面函數(shù)。針對(duì)T3B平板邊界層算例，如圖24所示，其中New_wf為基于新型壁面函數(shù)方法得>到的計(jì)算結(jié)果，F(xiàn)ull_grd為全網(wǎng)格下的摩阻計(jì)算結(jié)果，Exp代表實(shí)驗(yàn)參照結(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é)、中國空氣動(dòng)力研究與發(fā)展中心等單位在笛卡爾網(wǎng)格技術(shù)優(yōu)化、非貼體笛卡爾網(wǎng)格物面邊界處理等方面深入開展的工作，主要得出了以下結(jié)論：

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

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

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

本文對(duì)黏性自適應(yīng)笛卡爾網(wǎng)格方法的國內(nèi)外研究進(jìn)展進(jìn)行了簡要梳理，重點(diǎn)對(duì)國家數(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ù)雜流動(dòng)模擬能力提供技術(shù)儲(chǔ)備。