基于機(jī)器學(xué)習(xí)的非結(jié)構(gòu)網(wǎng)格陣面推進(jìn)生成技術(shù)初探1)

王年華 魯 鵬 常興華 張來平

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

?(西南科技大學(xué)信息工程學(xué)院，四川綿陽 621010)

??(重慶文理學(xué)院智能制造工程學(xué)院，重慶 402160)

??(國防科技創(chuàng)新研究院無人系統(tǒng)技術(shù)研究中心，北京 100071)

引言

網(wǎng)格生成是計(jì)算流體力學(xué)(computational flui dynamics,CFD)數(shù)值計(jì)算的第一步，張涵信院士將網(wǎng)格生成列為CFD 研究的5 個(gè)“M”之一[1-2]；在NASA 的研究報(bào)告《CFD Vision 2030 Study:A Path to Revolutionary Computational Aerosciences》[3]中，“幾何與網(wǎng)格生成”被列為未來六大重要研究領(lǐng)域之一，網(wǎng)格生成在CFD 數(shù)值模擬中的作用和重要性可見一斑.

在現(xiàn)代CFD 應(yīng)用過程中，自動生成復(fù)雜構(gòu)型的高質(zhì)量網(wǎng)格(包括自適應(yīng))依然是一個(gè)重大挑戰(zhàn)性問題.自動化程度和網(wǎng)格質(zhì)量是網(wǎng)格生成過程中最重要的兩個(gè)問題[4-5].據(jù)文獻(xiàn)統(tǒng)計(jì)，網(wǎng)格生成通常占據(jù)整個(gè)計(jì)算周期大約60%的人力時(shí)間，高度自動化的網(wǎng)格生成方法無疑可以很大程度節(jié)約CFD 計(jì)算周期內(nèi)的人工成本.

除此之外，網(wǎng)格質(zhì)量的好壞直接影響計(jì)算結(jié)果的精準(zhǔn)度，尤其是在復(fù)雜外形湍流數(shù)值模擬中，需要在流動參數(shù)梯度大的區(qū)域內(nèi)加密網(wǎng)格；在邊界層內(nèi)、激波附近、分離區(qū)內(nèi)也需要高質(zhì)量的網(wǎng)格.現(xiàn)階段借助商業(yè)軟件(如Gridgen,Pointwise 或ICEM 等)生成網(wǎng)格，其網(wǎng)格質(zhì)量高度依賴網(wǎng)格生成人員的經(jīng)驗(yàn)，導(dǎo)致不同條件下生成的網(wǎng)格質(zhì)量各不相同，針對同一問題得到的數(shù)值模擬結(jié)果也常常較為分散，這也是近年來舉辦的AIAA 阻力預(yù)測會議(drag prediction workshop,DPW)[6-7]、高升力預(yù)測會議(high lift prediction workshop,HiLiftPW)[8-9]和我國航空CFD 可信度研討會(aeronautical credibility workshop,AeCW)[10-11]等CFD 驗(yàn)證與確認(rèn)會議均提供了官方基準(zhǔn)網(wǎng)格的原因，目的就是避免因人工網(wǎng)格生成差異導(dǎo)致數(shù)值模擬結(jié)果無法在同等的條件下進(jìn)行比較.隨著CFD 應(yīng)用領(lǐng)域越來越廣泛，應(yīng)用案例越來越復(fù)雜，人們逐漸認(rèn)識到現(xiàn)有網(wǎng)格生成技術(shù)的不足仍然制約著復(fù)雜外形的數(shù)值模擬能力，仍需開展高度自動化的高質(zhì)量網(wǎng)格生成技術(shù)研究.

近年來，以深度學(xué)習(xí)和深度強(qiáng)化學(xué)習(xí)為代表的人工智能方法在各行各業(yè)得到廣泛應(yīng)用，也取得了舉世矚目的成功.2016 年3 月以來，基于深度強(qiáng)化學(xué)習(xí)的人工智能圍棋程序AlphaGo 相繼以4:1 和3:0 打敗了人類圍棋高手李世石和柯潔.以此為契機(jī)，人工智能和機(jī)器學(xué)習(xí)高調(diào)進(jìn)入大眾視野，并獲得極高關(guān)注度.隨著人工智能、大數(shù)據(jù)、超級計(jì)算機(jī)的發(fā)展，人工智能方法正在逐漸成為改變?nèi)祟愇磥砩鐣闹匾ぞ?

基于深度學(xué)習(xí)和深度強(qiáng)化學(xué)習(xí)的人工智能技術(shù)已經(jīng)成功應(yīng)用于工業(yè)社會的多個(gè)領(lǐng)域，如語音識別、機(jī)器翻譯、圖像識別、自動駕駛、智能推薦、搜索引擎等等.在流體力學(xué)專業(yè)領(lǐng)域，許多學(xué)者在人工智能方法與流體力學(xué)方法的結(jié)合領(lǐng)域也開展了許多探索性工作[12].比如傳統(tǒng)POD(proper orthogonal decomposition)/PCA(principal component analysis)分析方法可以由經(jīng)過訓(xùn)練的多層神經(jīng)網(wǎng)絡(luò)替代，用于數(shù)據(jù)降階模型建立[13-15]；結(jié)合深度神經(jīng)網(wǎng)絡(luò)和傳統(tǒng)PIV 的新型粒子圖像測速方法在精度、分辨率和計(jì)算效率上比傳統(tǒng)的相關(guān)分析法和光流法更具優(yōu)勢，并最終達(dá)到商用PIV 軟件的水平[16]；在流場可視化領(lǐng)域，卷積神經(jīng)網(wǎng)絡(luò)被用于漩渦等流場特征的識別與提取及原位可視化的數(shù)據(jù)壓縮[17]；在氣動外形優(yōu)化方面，通過對翼型氣動數(shù)據(jù)庫進(jìn)行學(xué)習(xí)，神經(jīng)網(wǎng)絡(luò)可用于翼型氣動力預(yù)測和反設(shè)計(jì)[18]；在湍流建模方面，采用人工神經(jīng)網(wǎng)絡(luò)對湍流黏性系數(shù)進(jìn)行建模，可代替?zhèn)鹘y(tǒng)的湍流模型，對雷諾平均NS(Reynolds averaged Navier-Stokes,RANS)方程進(jìn)行封閉[19-21].在魚類自主游動控制、無人機(jī)自主飛行控制、微型水下機(jī)器人控制等[22]領(lǐng)域，深度強(qiáng)化學(xué)習(xí)也展現(xiàn)出強(qiáng)大的應(yīng)用前景.

深度學(xué)習(xí)方法[23]之所以能夠取得廣泛應(yīng)用，是因?yàn)槠浠谏顚由窠?jīng)網(wǎng)絡(luò)對復(fù)雜非線性關(guān)系進(jìn)行分層表示，通過大量樣本數(shù)據(jù)的訓(xùn)練，使其能夠掌握數(shù)據(jù)中的內(nèi)在規(guī)律.對于網(wǎng)格生成而言，經(jīng)歷了幾十年的發(fā)展和應(yīng)用，已經(jīng)積累了大量各種類型的網(wǎng)格數(shù)據(jù)，這些網(wǎng)格數(shù)據(jù)包含了網(wǎng)格生成人員對幾何、CFD 和流體物理知識的理解，如果采用大量現(xiàn)有網(wǎng)格數(shù)據(jù)進(jìn)行訓(xùn)練，基于人工神經(jīng)網(wǎng)絡(luò)的深度學(xué)習(xí)有望學(xué)習(xí)獲取網(wǎng)格中隱含的網(wǎng)格生成經(jīng)驗(yàn)和方法，形成魯棒、自動化、智能化的網(wǎng)格生成方法，將縮短CFD 數(shù)值模擬的計(jì)算周期，同時(shí)使網(wǎng)格生成人員得到一定程度地解放.

本文首先簡要介紹傳統(tǒng)非結(jié)構(gòu)/混合網(wǎng)格生成技術(shù)，并分析傳統(tǒng)方法的優(yōu)劣；簡要綜述人工智能方法在CFD 網(wǎng)格生成領(lǐng)域的研究進(jìn)展；隨后提出基于人工智能方法的網(wǎng)格生成技術(shù)有待解決的問題，最后介紹基于人工神經(jīng)網(wǎng)絡(luò)的陣面推進(jìn)法生成非結(jié)構(gòu)網(wǎng)格的初步探索，以期為人工智能方法在網(wǎng)格生成技術(shù)中的發(fā)展提供參考.

1 傳統(tǒng)非結(jié)構(gòu)/混合網(wǎng)格生成技術(shù)

在CFD 領(lǐng)域，最早得到發(fā)展與應(yīng)用的是結(jié)構(gòu)網(wǎng)格，結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)之間的連接關(guān)系存在隱含的順序，可以在幾何空間進(jìn)行維度分解，并可以通過各方向的指標(biāo)(i,j,k)增減直接得到對應(yīng)的連接關(guān)系，數(shù)據(jù)直接采用多維數(shù)組進(jìn)行存儲，如x(i,j,k)，y(i,j,k)，z(i,j,k)，結(jié)構(gòu)網(wǎng)格數(shù)據(jù)結(jié)構(gòu)和算法實(shí)現(xiàn)均簡單直接.但是隨著模擬的問題越來越復(fù)雜，結(jié)構(gòu)網(wǎng)格對復(fù)雜外形的適應(yīng)性不夠，導(dǎo)致復(fù)雜外形結(jié)構(gòu)網(wǎng)格生成困難，因而CFD 研究人員提出了各種解決辦法，如結(jié)構(gòu)分塊網(wǎng)格、對接或拼接網(wǎng)格、重疊網(wǎng)格等等，同時(shí)靈活性更高、適應(yīng)性更強(qiáng)的非結(jié)構(gòu)網(wǎng)格也得到了重視和迅速發(fā)展.

非結(jié)構(gòu)網(wǎng)格節(jié)點(diǎn)間的關(guān)聯(lián)不存在直接的順序關(guān)系，其節(jié)點(diǎn)和單元編號在空間是隨機(jī)分布的，節(jié)點(diǎn)和單元信息采用特定的數(shù)據(jù)結(jié)構(gòu)存儲，因此非結(jié)構(gòu)網(wǎng)格存在靈活性高的優(yōu)點(diǎn)，但也存在數(shù)據(jù)結(jié)構(gòu)復(fù)雜、存儲和計(jì)算效率低、黏性流動模擬較難達(dá)到高精度等缺點(diǎn).正因?yàn)榻Y(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格各有優(yōu)劣，結(jié)合兩者優(yōu)勢的混合網(wǎng)格便應(yīng)運(yùn)而生.針對復(fù)雜外形黏性流動數(shù)值模擬，在邊界層內(nèi)生成結(jié)構(gòu)或半結(jié)構(gòu)的各向異性四邊形網(wǎng)格(二維)或三棱柱/六面體(三維)，而在遠(yuǎn)離物面處填充各向同性的三角形或四面體網(wǎng)格，這樣既保證了邊界層內(nèi)黏性流動的模擬精度，又便于各向同性三角形/四面體的自動化生成，在保證計(jì)算精度的同時(shí)大大減少了網(wǎng)格生成工作量.

Baker 在總結(jié)了各種網(wǎng)格生成方法之后，認(rèn)為非結(jié)構(gòu)三棱柱/四面體混合網(wǎng)格在黏性流動模擬精度和易用性上達(dá)到了很好的平衡(圖1 所示)，是未來網(wǎng)格技術(shù)的發(fā)展趨勢[5].近年來，在遠(yuǎn)場網(wǎng)格生成時(shí)大量采用自適應(yīng)Cartesian 網(wǎng)格技術(shù)，“多網(wǎng)格”技術(shù)的綜合應(yīng)用無疑能更好地發(fā)揮各種網(wǎng)格技術(shù)的優(yōu)勢.

傳統(tǒng)的非結(jié)構(gòu)網(wǎng)格生成技術(shù)主要有[24]：陣面推進(jìn)法、Delaunay 方法和四叉樹/八叉樹方法，以及基于上述三類方法進(jìn)行改進(jìn)或推廣[25-26]，如非結(jié)構(gòu)四邊形網(wǎng)格的生成方法可以采用三角形網(wǎng)格的分解或者合并(間接法，如圖2)，或以拓展陣面推進(jìn)的思想，直接在計(jì)算域內(nèi)逐個(gè)生成[27]或者一排一排生成四邊形單元(也稱作鋪路法或?qū)油七M(jìn)法)等.除此之外，四邊形網(wǎng)格還可以通過幾何分解和子域網(wǎng)格生成兩個(gè)步驟直接生成[28].

圖2 四邊形網(wǎng)格的間接法生成示意圖[25]Fig.2 Sketch of quadrangle generation by the indirect method[25]

對于邊界層內(nèi)采用結(jié)構(gòu)或半結(jié)構(gòu)的四邊形(2D)、六面體或者三棱柱網(wǎng)格(3D)的非結(jié)構(gòu)混合網(wǎng)格，邊界層網(wǎng)格剖分可以采用層推進(jìn)[29]、求解偏微分方程法[30-31]及基于各向異性四面體聚合的三棱柱網(wǎng)格生成方法等[32-33]，如圖3 所示.

圖3 四面體聚合生成三棱柱網(wǎng)格過程[32]Fig.3 Merge of tetrahedra into a prism[32]

上述方法中，陣面推進(jìn)法流程清晰、對幾何邊界的適應(yīng)能力強(qiáng)、通過控制推進(jìn)步長進(jìn)而控制網(wǎng)格疏密分布，但是其相交性判斷較為繁瑣，網(wǎng)格質(zhì)量仍有優(yōu)化空間.Delaunay 方法能夠得到盡可能高質(zhì)量的三角形網(wǎng)格，效率也較陣面推進(jìn)法高，但是不能保證計(jì)算域邊界的完整性，需要設(shè)計(jì)合理的邊界恢復(fù)算法.依靠叉樹結(jié)構(gòu)生成非結(jié)構(gòu)網(wǎng)格效率較高，但是生成黏性邊界層內(nèi)的貼體網(wǎng)格比較困難.而生成邊界層網(wǎng)格的層推進(jìn)法在處理復(fù)雜外形幾何曲率劇烈變化處仍然難以避免網(wǎng)格相交性判斷.求解偏微分方程可以得到質(zhì)量較高的邊界層網(wǎng)格，但是在物面有凹凸、棱角等間斷信息時(shí)容易求解發(fā)散，導(dǎo)致網(wǎng)格生成失敗.

總之，上述幾種方法各有優(yōu)劣，這也使得發(fā)展新的網(wǎng)格生成算法成為迫切需求.2017 年AIAA 舉辦的第一屆幾何與網(wǎng)格生成研討會(Geometry and Mesh Generation Workshop,GMGW-1)就指出應(yīng)該鼓勵(lì)非傳統(tǒng)的網(wǎng)格生成方法參會[34]，針對未來CFD 的研究與應(yīng)用，結(jié)合當(dāng)前已有的網(wǎng)格生成方法或者基于全新思路，進(jìn)一步發(fā)展自動、高效、魯棒的網(wǎng)格生成技術(shù)仍是十分有價(jià)值的研究方向.

2 基于機(jī)器學(xué)習(xí)方法的非結(jié)構(gòu)網(wǎng)格生成技術(shù)

盡管機(jī)器學(xué)習(xí)方法在社會生活，甚至流體力學(xué)專業(yè)領(lǐng)域都得到了大量的應(yīng)用，但是這些方法在CFD 網(wǎng)格生成方面的應(yīng)用仍極為少見.

2.1 研究進(jìn)展

從20 世紀(jì)90 年代開始，隨著人工神經(jīng)網(wǎng)絡(luò)的發(fā)展，一些研究人員就逐漸嘗試將神經(jīng)網(wǎng)絡(luò)應(yīng)用于網(wǎng)格生成領(lǐng)域.這些工作主要可以歸為以下幾類：

(1)采用自組織特征映射神經(jīng)網(wǎng)絡(luò)進(jìn)行網(wǎng)格重建.比如Ahn 等[35]提出基于自組織特征射(selforganizing feature mapping,SOM)神經(jīng)網(wǎng)絡(luò)進(jìn)行有限元網(wǎng)格的生成，通過將權(quán)重設(shè)置為節(jié)點(diǎn)坐標(biāo)，自組織網(wǎng)絡(luò)可以通過調(diào)節(jié)權(quán)重自動將神經(jīng)網(wǎng)絡(luò)拓?fù)涮卣饔成涞綐颖緮?shù)據(jù)點(diǎn)，樣本數(shù)據(jù)自動形成含拓?fù)潢P(guān)系的非結(jié)構(gòu)網(wǎng)格；Alfonzetti 等[36-39]采用了一種能夠根據(jù)樣本分布特征自動增加或減少輸出層神經(jīng)元的增長型SOM 神經(jīng)網(wǎng)絡(luò)——Let-It-Grow(LIG)神經(jīng)網(wǎng)絡(luò)，根據(jù)密度分布函數(shù)和邊界信息自動逐步生長出有限元網(wǎng)格，同時(shí)能夠保持樣本的拓?fù)浜兔芏确植迹粎魏陱?qiáng)等[40]將SOM神經(jīng)網(wǎng)絡(luò)引入多重網(wǎng)格稀疏網(wǎng)格的生成，很好地保持了稀網(wǎng)格和密網(wǎng)格之間的拓?fù)浜兔芏确植缄P(guān)系.大連理工大學(xué)陳先華[41]對LIG 網(wǎng)絡(luò)進(jìn)行了改進(jìn)并提高了網(wǎng)格質(zhì)量；張偉等[42]基于SOM 神經(jīng)網(wǎng)絡(luò)對三維散點(diǎn)數(shù)據(jù)進(jìn)行精簡和三角形曲面網(wǎng)格進(jìn)行了重建；Jilani 等[43]基于SOM 神經(jīng)網(wǎng)絡(luò)生成的二維初始網(wǎng)格進(jìn)行有限元分析，并計(jì)算出應(yīng)力集中部位的網(wǎng)格自適應(yīng)參數(shù)，再次采用SOM 神經(jīng)網(wǎng)絡(luò)對初始網(wǎng)格進(jìn)行自適應(yīng).圖4 為SOM 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖.

圖4 SOM 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)示意圖[44]Fig.4 Structure of SOM neural network[44]

SOM 神經(jīng)網(wǎng)絡(luò)的輸出層(也稱作競爭層)具有網(wǎng)格拓?fù)?如二維的三角形拓?fù)浜退倪呅瓮負(fù)?，通過樣本數(shù)據(jù)的學(xué)習(xí)，能將樣本數(shù)據(jù)在輸出層映射成一維或多維網(wǎng)格；網(wǎng)絡(luò)通過對輸入樣本的反復(fù)學(xué)習(xí)可以使權(quán)重向量空間與輸入樣本的概率分布趨于一致，即概率保持性.其生成網(wǎng)格的拓?fù)溆缮窠?jīng)網(wǎng)絡(luò)輸出層的拓?fù)浣Y(jié)構(gòu)決定，因而比較適合拓?fù)浜唵蔚母飨蛲匀切位蛘咚倪呅尉W(wǎng)格的生成，而較難處理具有復(fù)雜拓?fù)浣Y(jié)構(gòu)的混合網(wǎng)格.同時(shí)SOM 神經(jīng)網(wǎng)絡(luò)還需要輸入具有合理密度分布的樣本點(diǎn)，這對各向異性網(wǎng)格也比較困難，因而也難以生成計(jì)算流體力學(xué)中常用各向異性網(wǎng)格.對于三維問題，尤其是混合網(wǎng)格，SOM神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)也會十分復(fù)雜.

(2)基于多層神經(jīng)網(wǎng)絡(luò)進(jìn)行網(wǎng)格密度預(yù)測.如Lowther 等[45]將多層神經(jīng)網(wǎng)絡(luò)應(yīng)用于有限元網(wǎng)格自適應(yīng)中的網(wǎng)格密度預(yù)測，通過訓(xùn)練神經(jīng)網(wǎng)絡(luò)建立起幾何外形或流場特征與網(wǎng)格密度分布之間的關(guān)系，可以為網(wǎng)格節(jié)點(diǎn)分布和自適應(yīng)提供更好的依據(jù)，提高網(wǎng)格自適應(yīng)的效率和計(jì)算結(jié)果的精準(zhǔn)度.Chedid 等[46]也通過前向神經(jīng)網(wǎng)絡(luò)建立幾何特征與網(wǎng)格密度之間的關(guān)系，為自組織映射神經(jīng)網(wǎng)絡(luò)提供網(wǎng)格密度分布作為樣本輸入.Zhang 等[47]通過在不同密度的均勻網(wǎng)格上進(jìn)行誤差評估，建立起邊界特征、流場特征與網(wǎng)格密度分布的關(guān)系，并采用人工神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)這種關(guān)系，實(shí)現(xiàn)網(wǎng)格密度分布的預(yù)測.

(3)采用人工神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)網(wǎng)格生成規(guī)則，直接生成新網(wǎng)格點(diǎn)或單元.Yao 等[48]最早于2005 年基于人工神經(jīng)網(wǎng)絡(luò)的單元提取法(element extraction method)直接生成網(wǎng)格，采用反向傳播算法訓(xùn)練深度神經(jīng)網(wǎng)絡(luò)進(jìn)行二維有限元網(wǎng)格剖分研究，通過神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)四邊形網(wǎng)格生成規(guī)則，訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)可以替代原先的復(fù)雜繁瑣的條件判斷，在二維情況下取得了較好的生成效果.該方法中神經(jīng)網(wǎng)絡(luò)基本取代了傳統(tǒng)的網(wǎng)格生成規(guī)則，發(fā)揮了人工神經(jīng)網(wǎng)絡(luò)非線性表達(dá)能力較強(qiáng)的優(yōu)勢.但是該方法基于所謂的單元提取法，并非網(wǎng)格生成領(lǐng)域常用的方法，自動化程度有待提高，擴(kuò)展到三維、混合網(wǎng)格和各向異性網(wǎng)格時(shí)也會存在較大難度.同時(shí)，該方法中神經(jīng)網(wǎng)絡(luò)的訓(xùn)練并非自動提取大量樣本數(shù)據(jù)進(jìn)行訓(xùn)練，而是手動提取少量樣本進(jìn)行反復(fù)訓(xùn)練.

除此之外，在CAD 模型修補(bǔ)簡化、網(wǎng)格質(zhì)量檢測與優(yōu)化等前處理步驟中，也有一些研究人員嘗試應(yīng)用機(jī)器學(xué)習(xí)技術(shù).例如，Danglade 等[49]采用機(jī)器學(xué)習(xí)方法對CAD 數(shù)模進(jìn)行特征簡化；Stadler 等[50]將人工神經(jīng)網(wǎng)絡(luò)應(yīng)用于網(wǎng)格變形，并與基于徑向基函數(shù)(RBF)和反距離權(quán)插值的變形方法進(jìn)行了比較；高翔[51]也提出了基于支持向量機(jī)的并行網(wǎng)格變形算法；黃尚坤[52]發(fā)展了基于卷積神經(jīng)網(wǎng)絡(luò)的網(wǎng)格質(zhì)量判別技術(shù).

除此之外，最近10 年機(jī)器學(xué)習(xí)在網(wǎng)格生成領(lǐng)域的研究工作較為罕見，這與當(dāng)下機(jī)器學(xué)習(xí)在流體力學(xué)領(lǐng)域甚至人工智能領(lǐng)域的研究熱潮并不相稱.通過前期調(diào)研和初步研究，發(fā)現(xiàn)在將機(jī)器學(xué)習(xí)方法推廣應(yīng)用于CFD 網(wǎng)格生成時(shí)，仍存在一些尚未解決的關(guān)鍵問題，這也是機(jī)器學(xué)習(xí)在網(wǎng)格生成領(lǐng)域的研究工作罕見的原因.

2.2 有待解決的問題

(1)如何建立通用的基于機(jī)器學(xué)習(xí)的計(jì)算流體力學(xué)網(wǎng)格生成方法? 自組織神經(jīng)網(wǎng)絡(luò)不適合混合網(wǎng)格和各向異性網(wǎng)格的生成，因而只適用于簡單各向同性網(wǎng)格生成.傳統(tǒng)意義上，CFD 網(wǎng)格大多還是基于陣面推進(jìn)法、Delaunay 方法、四叉樹/八叉樹方法、層推進(jìn)法和求解偏微分方程等方法生成，因此在設(shè)計(jì)基于機(jī)器學(xué)習(xí)的網(wǎng)格生成方法時(shí)，必須充分考慮到CFD 網(wǎng)格特征及其傳統(tǒng)生成方法的特點(diǎn)和優(yōu)劣.

(2)如何在現(xiàn)有網(wǎng)格數(shù)據(jù)中自動識別和提取網(wǎng)格特征并構(gòu)造訓(xùn)練樣本數(shù)據(jù)集？網(wǎng)格特征往往與幾何特征和流動特征相關(guān)聯(lián)，如網(wǎng)格類型和網(wǎng)格疏密等特征一般由物面幾何曲率、流動物理量梯度等決定，這些關(guān)聯(lián)關(guān)系也是網(wǎng)格生成經(jīng)驗(yàn)的具體體現(xiàn).要采用神經(jīng)網(wǎng)絡(luò)建立起網(wǎng)格類型、網(wǎng)格疏密分布與物面幾何曲率、流動物理量梯度等之間的關(guān)系，就需要設(shè)計(jì)恰當(dāng)?shù)臉颖緮?shù)據(jù)格式，確保這些關(guān)聯(lián)關(guān)系和網(wǎng)格生成經(jīng)驗(yàn)仍然能保留；此外，手工提取大量訓(xùn)練數(shù)據(jù)是不現(xiàn)實(shí)，需要實(shí)現(xiàn)樣本數(shù)據(jù)的自動提取.

(3)如何通過深度學(xué)習(xí)訓(xùn)練提取大量網(wǎng)格數(shù)據(jù)中的特征以及如何使用深度學(xué)習(xí)生成網(wǎng)格？深度學(xué)習(xí)具有較強(qiáng)的非線性擬合能力，而神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)和參數(shù)對學(xué)習(xí)效率和學(xué)習(xí)效果影響較大，需要進(jìn)一步探索優(yōu)化神經(jīng)網(wǎng)絡(luò)參數(shù)，使得訓(xùn)練更快更準(zhǔn)確.此外，參考文獻(xiàn)中的方法，可采用深度學(xué)習(xí)進(jìn)行網(wǎng)格密度分布的學(xué)習(xí)訓(xùn)練.但是在訓(xùn)練完成后，如何將訓(xùn)練好的神經(jīng)網(wǎng)絡(luò)用于復(fù)雜外形的網(wǎng)格生成，是否需要結(jié)合傳統(tǒng)方法進(jìn)行，如結(jié)合陣面推進(jìn)法、Delaunay 方法、四叉樹/八叉樹方法，或是基于全新的方法等.

(4)在CAD 模型修補(bǔ)簡化和網(wǎng)格質(zhì)量檢測與優(yōu)化等前處理步驟中，探索應(yīng)用機(jī)器學(xué)習(xí)技術(shù)，解決傳統(tǒng)方法存在的弊端.

3 基于人工神經(jīng)網(wǎng)絡(luò)的非結(jié)構(gòu)網(wǎng)格陣面推進(jìn)生成技術(shù)初探

由于陣面推進(jìn)法具備較好的通用性，能夠較為方便地推廣到混合網(wǎng)格、各向異性網(wǎng)格和三維情況，這是其他傳統(tǒng)網(wǎng)格生成方法不具備的優(yōu)勢.通過結(jié)合傳統(tǒng)的陣面推進(jìn)法，本文在基于人工神經(jīng)網(wǎng)絡(luò)的非結(jié)構(gòu)網(wǎng)格生成技術(shù)方面進(jìn)行了初步探索.主要包括3 方面內(nèi)容：(1)網(wǎng)格訓(xùn)練樣本數(shù)據(jù)格式的設(shè)計(jì)與數(shù)據(jù)集的自動提取；(2)適用于網(wǎng)格數(shù)據(jù)學(xué)習(xí)的人工神經(jīng)網(wǎng)絡(luò)設(shè)計(jì)與訓(xùn)練；(3)基于人工神經(jīng)網(wǎng)絡(luò)的三角形網(wǎng)格的陣面推進(jìn)生成.

3.1 訓(xùn)練樣本的設(shè)計(jì)與自動提取

本文基于陣面推進(jìn)法進(jìn)行網(wǎng)格生成，因此在設(shè)計(jì)訓(xùn)練樣本時(shí)就需要考慮到陣面推進(jìn)的特點(diǎn).

如圖5 所示，針對三角形網(wǎng)格，訓(xùn)練樣本考慮基準(zhǔn)邊(B1,B2)及其左鄰點(diǎn)(L1)和右鄰點(diǎn)(R1)，其目標(biāo)點(diǎn)(TP)與基準(zhǔn)邊構(gòu)成目標(biāo)三角形.因此可將訓(xùn)練樣本設(shè)計(jì)為輸入L1,B1,B2,R14 個(gè)點(diǎn)的坐標(biāo)，輸出目標(biāo)點(diǎn)TP的坐標(biāo)和生成模式，寫成數(shù)組形式為

圖5 三角形網(wǎng)格訓(xùn)練樣本Fig.5 Training sample of triangular grids

在原始陣面推進(jìn)法中需要進(jìn)行大量相交性判斷，新生成的邊不僅需要跟活躍陣面進(jìn)行相交性判斷，還要與已經(jīng)生成的非活躍網(wǎng)格邊進(jìn)行判斷，操作十分費(fèi)時(shí).為減少相交性判斷次數(shù)，參考文獻(xiàn)[48]的做法，本文引入網(wǎng)格生成模式(mode).其內(nèi)涵是：假設(shè)訓(xùn)練樣本的4 個(gè)模板點(diǎn)分別為L1，B1，B2和R1，其構(gòu)成的幾何關(guān)系可以根據(jù)其連線的夾角的大小人為分成3 類，即目標(biāo)點(diǎn)選擇生成的新點(diǎn)為模式1；目標(biāo)點(diǎn)不選擇新點(diǎn)而直接選擇L1點(diǎn)為模式2；目標(biāo)點(diǎn)直接選擇模板R1點(diǎn)為模式3，如圖6 所示是一種模式分類實(shí)例.

圖6 三角形網(wǎng)格生成模式Fig.6 Generation mode of triangular grids

采用人工神經(jīng)網(wǎng)絡(luò)進(jìn)行模式預(yù)測，對于模式2 和模式3，由于模式分類直接指定了目標(biāo)點(diǎn)，可省去候選點(diǎn)相交性判斷等篩選過程；直接與模板現(xiàn)有點(diǎn)相連也不會出現(xiàn)與非活躍邊相交的情況，無需與非活躍邊進(jìn)行相交性判斷，可以減少相交性判斷；比如，如圖7(a)所示，采用原始陣面推進(jìn)法時(shí)，新生成的邊必須與非活躍邊(邊1-2)進(jìn)行相交性判斷，否則陣面4-5 在生成時(shí)會選擇3 點(diǎn)，生成不符合要求的4-5-3 單元，而若是采用人工神經(jīng)網(wǎng)絡(luò)提前進(jìn)行生成模式預(yù)測，則神經(jīng)網(wǎng)絡(luò)能夠提前預(yù)測出生成模式2，從而減少了大量的與非活躍邊的相交性判斷，直接與現(xiàn)有點(diǎn)2 相連，生成4-5-2 單元，如圖7(b)所示.

圖7 三角形網(wǎng)格生成模式預(yù)測實(shí)例Fig.7 Example of generation mode on triangular grids

為消除不同位置訓(xùn)練樣本由于坐標(biāo)絕對位置差異帶來的影響，提升訓(xùn)練效果，對樣本數(shù)據(jù)進(jìn)行歸一化處理.具體來說，通過平移、縮放、旋轉(zhuǎn)等操作將B1，B2變換到(0,0)和(1,0)，而L1，R1，TP通過同樣的變換操作，變換到新的位置，分別記為因此歸一化之后的樣本為

此外，采用one-hot 向量將生成模式轉(zhuǎn)換為機(jī)器學(xué)習(xí)更易于利用的形式

為自動生成大量訓(xùn)練樣本，可以采用商業(yè)軟件或其它現(xiàn)有軟件生成基準(zhǔn)網(wǎng)格，對網(wǎng)格中每個(gè)邊進(jìn)行遍歷，以該邊作為基準(zhǔn)邊，自動提取其訓(xùn)練樣本.

對于如圖8 所示的同一基準(zhǔn)邊(node14,node15)、同一目標(biāo)點(diǎn)(node11)，但是基準(zhǔn)邊鄰近點(diǎn)不同的情況(鄰近點(diǎn)為node10 和node17，或node8 和node17，或node8 和node11，或node11 和node17)，本文將其視作不同的訓(xùn)練樣本，因此樣本數(shù)量將是網(wǎng)格邊數(shù)量的數(shù)倍.

圖8 三角形網(wǎng)格訓(xùn)練樣本提取實(shí)例Fig.8 Examples of sample extraction on triangular grids

如圖9 所示，對各向同性三角形網(wǎng)格進(jìn)行訓(xùn)練樣本自動提取，對于含有12 582 個(gè)邊的二維NACA0012 翼型三角形網(wǎng)格，經(jīng)過自動提取，可得到293 746 個(gè)樣本，其輸入和輸出分別為293 746×8和293 746×5 的兩個(gè)數(shù)組.

圖9 三角形基準(zhǔn)網(wǎng)格Fig.9 Triangular base grids for sample extraction

3.2 適用于網(wǎng)格數(shù)據(jù)學(xué)習(xí)的ANN 設(shè)計(jì)與訓(xùn)練

本文基于Matlab 神經(jīng)網(wǎng)絡(luò)訓(xùn)練工具設(shè)計(jì)全連接的人工神經(jīng)網(wǎng)絡(luò)，網(wǎng)絡(luò)含有1 個(gè)輸入層，1 個(gè)隱藏層，1 個(gè)輸出層.輸入層含有8 個(gè)神經(jīng)元(分別輸入4 個(gè)模板點(diǎn)的坐標(biāo))，輸出層含有5 個(gè)神經(jīng)元(分別輸出1 個(gè)新點(diǎn)坐標(biāo)和1 個(gè)生成模式onehot 向量)，隱藏層的神經(jīng)元數(shù)量為20 個(gè)，激活函數(shù)采用Sigmoid 函數(shù)，損失函數(shù)為均方誤差函數(shù)，訓(xùn)練方法采用Levenberg-Marquardt 反向傳播方法，該方法比常規(guī)的梯度下降反向傳播算法訓(xùn)練效率更高，具體可以參考文獻(xiàn)[53].神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)如圖10所示.

網(wǎng)格訓(xùn)練樣本按80%,10%和10%的比例隨機(jī)劃分為訓(xùn)練集、測試集和驗(yàn)證集，圖11 給出了在三角形網(wǎng)格訓(xùn)練集上的Loss 值及在驗(yàn)證集上的預(yù)測精度收斂歷程.結(jié)果顯示：經(jīng)過500 多次迭代后，Loss 值及預(yù)測精度均下降到了0.013 左右.

圖10 基于Matlab 的人工神經(jīng)網(wǎng)絡(luò)訓(xùn)練工具Fig.10 Artificia neural network training tool based on Matlab

圖11 三角形網(wǎng)格訓(xùn)練Loss值和精度收斂歷程Fig.11 Convergence of loss and accuracy on triangular grids

3.3 基于人工神經(jīng)網(wǎng)絡(luò)的非結(jié)構(gòu)網(wǎng)格陣面推進(jìn)生成

傳統(tǒng)陣面推進(jìn)法(advancing front method,AFT method)的基本思想是首先將計(jì)算域的邊界劃分為小的陣元，如二維情況下的線段、三維情況下的表面三角形.由此構(gòu)成初始陣面，然后選定某一陣元，將某一在計(jì)算域中新插入的網(wǎng)格節(jié)點(diǎn)或原陣面上已經(jīng)存在的點(diǎn)與該陣元相連構(gòu)成基本單元(二維時(shí)為三角形，三維時(shí)為四面體).初始陣面不斷向計(jì)算域中推進(jìn)，逐步填充整個(gè)計(jì)算域，如圖12 所示，具體算法流程可參考文獻(xiàn)[24]，本文不再贅述.

陣面推進(jìn)法中需要進(jìn)行大量候選點(diǎn)的相交性判斷等篩選過程，繁瑣耗時(shí).經(jīng)過訓(xùn)練的神經(jīng)網(wǎng)絡(luò)能夠根據(jù)輸入模板點(diǎn)，輸出生成模式和新點(diǎn)坐標(biāo)，結(jié)合生成模式可以省去繁瑣耗時(shí)的相交性判斷操作.因此本文結(jié)合陣面推進(jìn)的思路，逐步生成非結(jié)構(gòu)三角形網(wǎng)格，新方法流程如圖13 所示，本文稱之為基于人工神經(jīng)網(wǎng)絡(luò)的陣面推進(jìn)方法，簡稱ANN-AFT方法，具體方法流程如下：

(1)將幾何邊界離散成陣元；

(2)從最小陣面出發(fā)，自動選擇陣面左右相鄰的網(wǎng)格點(diǎn)作為模板點(diǎn)，歸一化處理后作為神經(jīng)網(wǎng)絡(luò)輸入，人工神經(jīng)網(wǎng)絡(luò)輸出生成模式和新點(diǎn)的歸一化坐標(biāo)；

圖13 基于ANN 的陣面推進(jìn)網(wǎng)格生成流程圖Fig.13 Flow chart of ANN-based advancing front method

(3)根據(jù)生成模式和新點(diǎn)坐標(biāo)生成新網(wǎng)格單元，具體來說：首先根據(jù)不同生成模式，確定是否使用新點(diǎn).若為模式1，則使用新點(diǎn)，根據(jù)局部網(wǎng)格步長控制確定新網(wǎng)格點(diǎn)在計(jì)算域中的絕對坐標(biāo)(“反歸一化”)，按照傳統(tǒng)陣面推進(jìn)法，將新點(diǎn)與其周圍一定范圍內(nèi)的陣面點(diǎn)搜集起來組成候選點(diǎn)集合，按照候選點(diǎn)構(gòu)成的網(wǎng)格質(zhì)量高低，依次進(jìn)行相交性判斷等篩選；若為模式2 或模式3，則直接選擇相應(yīng)的模板點(diǎn)，跳過步長控制和候選點(diǎn)相交性判斷；

(4)判斷新單元是否合適，比如單元網(wǎng)格質(zhì)量是否滿足要求等，如合適，則更新數(shù)據(jù)結(jié)構(gòu)，即更新點(diǎn)面體之間的關(guān)聯(lián)關(guān)系、標(biāo)記面是否活躍等；

(5)回到步驟2，直至所有面變成非活躍面，整個(gè)計(jì)算域被網(wǎng)格填滿.

通常，陣面推進(jìn)法生成的初始網(wǎng)格要經(jīng)過優(yōu)化才能得到高質(zhì)量的網(wǎng)格.在本文的算例中，引入Delaunay 準(zhǔn)則對質(zhì)量較差的單元進(jìn)行對角線變換，如圖14，同時(shí)引入彈簧松弛法對節(jié)點(diǎn)坐標(biāo)進(jìn)行優(yōu)化[27]，如圖15，以提升網(wǎng)格質(zhì)量.

陣面推進(jìn)法中必須給定網(wǎng)格步長分布，這里采用背景網(wǎng)格方法人為指定網(wǎng)格步長，具體方法不再贅述.

圖14 Delaunay 準(zhǔn)則進(jìn)行對角線變換Fig.14 Diagonal line exchange based on Delaunay criterion

圖15 彈簧法網(wǎng)格變形優(yōu)化Fig.15 Grid quality optimization by spring method

根據(jù)上述思路，本文生成了幾個(gè)簡單外形(矩形、圓柱、NACA0012 翼型和三段翼型)的非結(jié)構(gòu)三角形網(wǎng)格.如圖16 所示，初始網(wǎng)格中紅色邊為神經(jīng)網(wǎng)絡(luò)生成，藍(lán)色邊為選擇現(xiàn)有點(diǎn)生成的邊.

圖16 ANN-AFT 生成三角形網(wǎng)格Fig.16 Triangular grids generation based on ANN-AFT method

對圓柱、NACA0012 翼型和三段翼型網(wǎng)格進(jìn)行質(zhì)量檢測，考察三角形質(zhì)量參數(shù)α，及相鄰單元的面積比γ.其中三角形質(zhì)量參數(shù)α 的定義為

α 越接近1，網(wǎng)格質(zhì)量越好，正三角形α=1，表1 給出了網(wǎng)格質(zhì)量檢測結(jié)果，結(jié)果顯示：本文方法優(yōu)化生成的網(wǎng)格質(zhì)量最小α 和最小γ(表中記為min.)均比某知名商業(yè)軟件略高，而本文生成的網(wǎng)格(表中present 值)平均α 和平均γ(表中記為avg.)均與商業(yè)軟件(表中ref.值)接近，證明本文生成的三角形網(wǎng)格質(zhì)量較高.

表1 二維三角形網(wǎng)格質(zhì)量檢測Table 1 Quality statistics of 2D triangular grids

此外，還對原始陣面推進(jìn)法和基于人工神經(jīng)網(wǎng)絡(luò)的陣面推進(jìn)法耗時(shí)進(jìn)行比較，如表2 所示.與原始的陣面推進(jìn)相比，ANN-AFT 相交性判斷減少30%以上，推進(jìn)耗時(shí)節(jié)約30%，新方法具有更高的效率.預(yù)計(jì)在三維情況下，原始方法的相交性判斷量更大，ANN-AFT 方法能節(jié)約更多時(shí)間.

進(jìn)一步的，采用合并法生成非結(jié)構(gòu)三角形/四邊形混合網(wǎng)格，根據(jù)與相鄰單元合并后的新單元質(zhì)量β 的高低，選擇質(zhì)量最高的單元進(jìn)行合并，生成新的四邊形單元.四邊形單元的質(zhì)量β 定義為

式中，α 為四邊形沿兩條對角線分成的4 個(gè)三角形對應(yīng)的質(zhì)量參數(shù)，且滿足α1≥α2≥α3≥α4.凸四邊形0 <β≤1，越接近1 四邊形質(zhì)量越好，凹四邊形β<0.

圖17 和表3 給出了商業(yè)軟件生成的混合網(wǎng)格和本文生成的混合網(wǎng)格的對比結(jié)果.由表3 中數(shù)據(jù)可見，本文生成的混合網(wǎng)格的平均面積比γ 和平均質(zhì)量β 均與商業(yè)軟件生成的混合網(wǎng)格(表中ref.值)接近，證明本文的混合網(wǎng)格質(zhì)量達(dá)到了較好的水平.

圖17 ANN-AFT 生成三角形/四邊形混合網(wǎng)格Fig.17 Hybrid grids generation based on ANN-AFT method

表3 二維三角形/四邊形混合網(wǎng)格質(zhì)量檢測Table 3 Quality statistics of 2D hybrid grids

4 總結(jié)與展望

本文首先簡要回顧了傳統(tǒng)的非結(jié)構(gòu)/混合網(wǎng)格生成方法，隨后綜述了機(jī)器學(xué)習(xí)方法在非結(jié)構(gòu)網(wǎng)格生成領(lǐng)域的應(yīng)用，提出了在非結(jié)構(gòu)/混合網(wǎng)格生成中應(yīng)用機(jī)器學(xué)習(xí)方法存在的4 個(gè)方面的關(guān)鍵問題，最后，結(jié)合人工神經(jīng)網(wǎng)絡(luò)和陣面推進(jìn)法進(jìn)行了非結(jié)構(gòu)/混合網(wǎng)格生成的初步探索，主要包括：

(1)設(shè)計(jì)了非結(jié)構(gòu)網(wǎng)格樣本數(shù)據(jù)格式并實(shí)現(xiàn)了樣本數(shù)據(jù)集的自動提取；

(2)初步發(fā)展了一種基于人工神經(jīng)網(wǎng)絡(luò)的二維非結(jié)構(gòu)網(wǎng)格陣面推進(jìn)生成方法；

(3)采用新發(fā)展的方法生成了幾個(gè)典型二維各向同性非結(jié)構(gòu)混合網(wǎng)格并測試了網(wǎng)格質(zhì)量和生成耗時(shí).

結(jié)果顯示本文方法生成的網(wǎng)格質(zhì)量可以達(dá)到商業(yè)軟件的水平，且由于減少了相交性判斷，生成效率較傳統(tǒng)陣面推進(jìn)法更高，顯示人工智能方法在非結(jié)構(gòu)網(wǎng)格生成中有較好的應(yīng)用前景.

下一步將持續(xù)開展更加深入的研究，比如將本文方法推廣到三維情況.在三維情況下，陣面是三角形面，網(wǎng)格模板是多個(gè)面的集合，比二維更復(fù)雜；同時(shí)，三維時(shí)生成模式預(yù)計(jì)不止3 種，建立合理完整的生成模式分類也比二維更難，值得進(jìn)一步探索.

除此之外，結(jié)合人工神經(jīng)網(wǎng)絡(luò)與層推進(jìn)方法發(fā)展各向異性網(wǎng)格生成方法，并進(jìn)一步生成各向異性混合網(wǎng)格；也可考慮改進(jìn)樣本數(shù)據(jù)格式，引入網(wǎng)格夾角、到壁面的最小距離、增加相鄰點(diǎn)等參數(shù)；考察神經(jīng)網(wǎng)絡(luò)參數(shù)對訓(xùn)練結(jié)果的影響；通過改進(jìn)現(xiàn)有的陣面推進(jìn)法，并行推進(jìn)生成網(wǎng)格，進(jìn)一步提高網(wǎng)格生成效率；訓(xùn)練人工神經(jīng)網(wǎng)絡(luò)進(jìn)行網(wǎng)格密度分布預(yù)測；在CAD 數(shù)模修補(bǔ)與簡化等方面引入機(jī)器學(xué)習(xí)方法等.