張曉瀟，張 武，褚學(xué)森，王良軍，楊廣文

(1. 上海大學(xué)，上海 201900；2. 清華大學(xué)，北京 100084；3. 中國船舶科學(xué)研究中心，無錫 214082)

0 引 言

格子Boltzmann方法（Lattice Boltzmann Method,LBM）[1]是計算流體力學(xué)中一種重要的數(shù)值方法。不同于傳統(tǒng)的宏觀方法和微觀模型方法，LBM結(jié)合了介觀動理學(xué)以及微觀模型的特點，擁有很多獨特的優(yōu)勢，如編程簡單、具有天生的并行性、物理背景清晰等。近年來，LBM在很多領(lǐng)域都得到了廣泛應(yīng)用，如多孔介質(zhì)流[2-3]、多相流[4]、燃燒[5]、湍流[6]問題等。LBM計算的模型也從方腔、圓柱等簡單幾何體變?yōu)轱w機(jī)機(jī)翼、汽車、動車等形狀十分復(fù)雜的幾何體，而精確快速的幾何前處理是實現(xiàn)LBM計算工程應(yīng)用的前提條件。

目前，CAD（Computational Aided Design）軟件在汽車、輪船等領(lǐng)域得到了廣泛應(yīng)用，在此基礎(chǔ)上，人們制定了很多通用的數(shù)據(jù)格式。流體計算中模型生成的文件是立體光刻（Stereolithographic, STL）數(shù)據(jù)類型，該通用數(shù)據(jù)文件格式是為了在不同的軟件、不同的用途之間進(jìn)行模型數(shù)據(jù)交換。這種文件格式在快速成型（Rapid Prototyping, RP）技術(shù)中有著廣泛的應(yīng)用，尤其是3D打印中。基于哈希表的拓?fù)渲貥?gòu)算法能夠快速的重建拓?fù)浣Y(jié)構(gòu)[7]、基于分層的優(yōu)化算法[8-9]能夠提高STL文件分層效率。文獻(xiàn)[10]主要對在3D打印過程中模型支撐架的設(shè)計提出算法，使得模型重建易于實現(xiàn)。文獻(xiàn)[11]對模型的表面進(jìn)行了擬合，使其更加平滑。文獻(xiàn)[12]基于自適應(yīng)層深度法向圖像針對模型重建問題提出了一個有效的方法。文獻(xiàn)[13]能夠很好地識別模型的輪廓，以及判定外輪廓以及內(nèi)輪廓。Michael等[14]以八叉樹為基本數(shù)據(jù)結(jié)構(gòu)，對GPU的特點提出了并行算法。文獻(xiàn)[15]對三維多邊形網(wǎng)格的26-鄰接體素模型提出了快速算法。

在利用LBM方法進(jìn)行流體力學(xué)計算的過程中，物面邊界的處理通常有三種方法：浸沒邊界法、插值曲面邊界、階梯狀近似的直接反彈法。三種方法中前兩種精度較高，但處理相對復(fù)雜，后一種精度稍低但處理簡單，數(shù)值穩(wěn)定性好，非常適合超大規(guī)模并行計算，在網(wǎng)格數(shù)量足夠的情況下同樣可以獲得準(zhǔn)確的結(jié)果。本文主要針對直接反彈法設(shè)計LBM的前處理算法。直接反彈法中需要用到流體計算區(qū)域內(nèi)的格點相對模型的位置關(guān)系，判斷計算流體點是在物體內(nèi)部、表面或者外部。LBM計算過程中物體外部點為流場點實施正常的碰撞遷移操作，內(nèi)部點為固體點不參與計算，表面點為物面邊界點實施反彈操作。

為了能夠根據(jù)STL幾何文件得到LBM計算所需的格點類型信息，完成LBM計算的前處理，本文提出了一種快速、高效、準(zhǔn)確的方法。該算法主要目標(biāo)是精簡計算量，提高計算效率，并能保證計算流體模型能夠很好地還原原始物體幾何特征。關(guān)鍵點包括：

（1） 讀取通用的物體模型文件即STL文件，利用STL文件存儲的信息生成物體原始模型；

（2）通過任意視角（視角指x、y、z三個方向）得到物體的邊界點，此邊界點是計算區(qū)域的格點，根據(jù)邊界點重構(gòu)幾何計算模型；

（3） 通過優(yōu)化判別方法，減少計算量，提高邊界點判定的效率；

（4）算法本身具有良好的并行性，為以后程序的并行提供了基礎(chǔ)。

1 STL文件格式

STL文件格式是由3D SYSTEMS公司在1988年制定的一個標(biāo)準(zhǔn)接口協(xié)議，是大多數(shù)CAD軟件通用的文件格式。STL文件由多個三角形面組合而成。其類型可以進(jìn)一步分為文本文件（ASCII格式）和二進(jìn)制文件（BINARY）。

（1）二進(jìn)制文件格式。二進(jìn)制文件利用固定的字符長度來給出三角面片的幾何信息。前80個字節(jié)存儲幾何相關(guān)的注釋信息。隨后存儲三角形的總數(shù)目，之后存儲每個三角形的法向量坐標(biāo)、根據(jù)右手法則排列的三角形頂點信息，兩兩三角形信息之間由兩個空格分開。該格式相對于文本文件格式更加簡潔，占用磁盤空間更少，后期對文件的操作相對更方便。

（2）文本文件格式。文本格式第一行以solid開始，最后一行以endsolid標(biāo)識。中間部分每7行存儲一個三角形的信息。首先存儲三角形法向量的三個坐標(biāo)，然后按照右手法則依次存儲三角形三個頂點的坐標(biāo)信息。這種格式存儲直觀，程序讀取方便。本文方法主要基于文本文件格式開展。

2 快速生成流體計算模型算法

計算模型重建的主要目的是得到計算區(qū)域內(nèi)各網(wǎng)格點的類型。格點類型可分為幾何體內(nèi)部點、幾何體外部點、幾何體邊界點。目前比較常用的方法有：

（1）數(shù)學(xué)表達(dá)式法。主要針對規(guī)則的幾何形狀，如球、立方體、橢球等能用數(shù)學(xué)表達(dá)式直接表示的幾何外形。這種方法可以根據(jù)表達(dá)式直接得到格點的類型，效率最高。但是它適用范圍很窄，無法較好解決現(xiàn)實中的復(fù)雜問題。

（2）直接法。例如文獻(xiàn)[16]中實現(xiàn)模型的重建過程都是利用的直接法。它首先找到一個比幾何體要大的包圍盒把幾何體包住，然后在這個包圍盒中遍歷所有的格點并對這些格點進(jìn)行判斷，包圍盒外部所有的格點都是幾何體外部點。在對包圍盒中格點的遍歷中，需要根據(jù)格點與所有的STL文件中三角形的位置關(guān)系來判斷這個格點是幾何內(nèi)部點或者是幾何外部點。遍歷完所有格點后，就可以得到計算區(qū)域內(nèi)幾何體的模型。

圖1（a）中，方形區(qū)域內(nèi)的格點為需要計算的格點范圍，圖1（b）展示的是圓球模型中對一個格點P1位置判定過程所有需要計算的三角形。為了顯示的直觀性，圖中給出的是圓球的一個投影面。其中，圖1（b）中P1點是需判定的格點，P0是模型內(nèi)部任意一點，因為圓球的對稱性，因此P0點選為球心。圖2是直接法判定格點類型的流程圖。

圖1 直接法示意圖Fig. 1 Schematic of the direct method

圖2 格點類型判斷流程圖Fig. 2 Flowchart for the grid type identification

直接法可以直接得到流場內(nèi)所有網(wǎng)格點的類型，但由于其遍歷次數(shù)的影響，在幾何體面網(wǎng)格復(fù)雜和計算區(qū)域較大的情況下，該方法的時間消耗很大，會占用流體計算過程的很大部分，影響了整個程序的效率。

（3）本文方法。本文提出了一種快速生成流體計算模型的算法。圖3（a）中陰影區(qū)域是計算格點區(qū)域。圖3（b）圓環(huán)中格點是需要計算的格點。只需要計算格點距離模型（內(nèi)圓）的距離，最近的格點即為所求邊界點。圖3（a）中，需要計算的格點均位于黑色圓環(huán)區(qū)域，該區(qū)域遠(yuǎn)小于直接法的方形區(qū)域。圖3（b）中的黑點是需要計算的格點，其在數(shù)量上也小于直接法的面三角形數(shù)。

圖3 本文提出的算法示意圖Fig. 3 Schematic of the algorithm proposed in this study

算法主要過程如圖4所示。算法主要步驟如下：

（1）為了縮小計算量，只計算由三角形生成包圍盒內(nèi)部的點（三角形包圍盒是指邊平行于坐標(biāo)軸的體積最小的長方體，長方體的頂點是計算格點，且三角形的三個頂點都在長方體的內(nèi)部）。為了進(jìn)一步減少計算量，繼續(xù)排除不滿足條件的格點，采取了以下處理方法：1）去除三角形對應(yīng)在物體內(nèi)部的格點：這里判斷的方法是，以靠近物體內(nèi)部的格點為終點，以三角形頂點為起點，構(gòu)成的向量與三角形的法向量方向相反；2）去除不是該三角形對應(yīng)的點：同時將格點和三角形沿著某一坐標(biāo)軸投影，如果格點的投影不在三角形的投影內(nèi)部，那么這個格點不是三角形對應(yīng)的格點。

（2）在篩選后的格點中，對于某個坐標(biāo)軸選取距離三角形最近的格點，例如選取x軸，則對于一個（y,z）坐標(biāo)選取距離x軸最近的格點，該距離可以直接用格點到三角形的投影距離，然后即可得到該三角形對應(yīng)的所有物體邊界點。

（3）對于每個三角形重復(fù)前面兩步，得到物體所有的邊界點。

（4）進(jìn)一步生成邊界點對，邊界點對是兩個邊界點構(gòu)成的點對，且這兩個點所確定的直線平行于坐標(biāo)軸。邊界點對中間的格點是物體內(nèi)部點。

（5）根據(jù)邊界點對可以確定所有的內(nèi)部點，即可完成物體模型的重建。

圖4 本文算法流程圖Fig. 4 Flowchat of the algorithm proposed in this study

3 算法應(yīng)用特性

通過測試，該算法能夠用于大部分模型文件讀取，并快速進(jìn)行模型重建。具體特點包括：

（1）不嚴(yán)格要求物體幾何封閉。對于模型STL文件中存在的平行于某一坐標(biāo)軸的三角形，該算法不需要這些三角形便可以完成模型的重建，即對STL文件的要求不是十分嚴(yán)格。

（2）能夠處理部分不規(guī)范網(wǎng)格文件。正常情況下網(wǎng)格文件只包含物體表面三角形，但是部分模型網(wǎng)格文件存在物體內(nèi)部三角形。而直接法有一個重要的假定就是網(wǎng)格三角形是模型內(nèi)部與外部的分界點，此時冗余的內(nèi)部三角形信息將會導(dǎo)致直接法無法正常判斷網(wǎng)格類型。而本文算法能夠較好地解決這類問題。

（3）能夠快速處理標(biāo)準(zhǔn)模型的重建。該算法是根據(jù)物體的面網(wǎng)格來確定物體的邊界格點，然后根據(jù)邊界格點重建物體模型。所以運行時間只與物體STL文件中三角形個數(shù)有關(guān)，在相同數(shù)量的面網(wǎng)格三角形下，運行時間基本不受物體幾何復(fù)雜度影響。

4 測試與結(jié)果比較

4.1 測試環(huán)境

測試使用“神威·太湖之光”計算機(jī)集群系統(tǒng)進(jìn)行測試，以“申威 26010”異構(gòu)眾核處理器為核心，訪存帶寬136.51 GB/s，主存容量32 GB，網(wǎng)絡(luò)接口雙向帶寬16 GB/s。采用sw5CC進(jìn)行程序編譯。

4.2 測試案例

測試中選取了3種比較有代表性的物體模型對直接法和本文算法進(jìn)行比較。三種模型分別為：圓球、NACA0012翼型、CHN-T1飛機(jī)標(biāo)模。圓球模型是一種簡單的模型，其形狀分布均勻、對稱、形狀簡單，網(wǎng)格量少。NACA0012翼型代表了形狀分布不均勻、表面三角形密集、網(wǎng)格量較大的模型，本文用到的是一種更為特殊的不封閉的NACA0012翼型模型，該模型只有型線面上有網(wǎng)格，側(cè)面為空。CHN-T1飛機(jī)標(biāo)模屬于復(fù)雜的模型，它結(jié)構(gòu)復(fù)雜、不同的部位的三角形有不同大小、網(wǎng)格量巨大。程序測試結(jié)果分別如下圖。圖5（a）、5（b）是圓球的面網(wǎng)格和經(jīng)過算法重建的模型，圖6（a）、6（b）是NACA0012的面網(wǎng)格和重建模型結(jié)果，圖7（a）、7（b）是CHN-T1飛機(jī)標(biāo)模的面網(wǎng)格和重建模型。由圖5、圖6、圖7中的結(jié)果可以看出，本文算法有很好的整體重建能力，能夠反應(yīng)原模型的復(fù)雜度，保留原有模型幾何特征。

圖5 圓球模型Fig. 5 Sphere model

圖6 NACA0012翼型Fig. 6 NACA0012 airfoil

圖7 CHN-T1飛機(jī)標(biāo)模Fig. 7 CHN-T1 aircraft model

4.3 測試結(jié)果比較

在神威平臺上的測試結(jié)果如表1所示。兩個算法除了判斷處理部分步驟不同，其他都相同。從表1中可以看出，本文算法相對于直接法有著明顯的速度提升。在圓球、NACA0012翼型測試中，直接法和本文算法均采用單核測試，判斷時間由525.4 s和5 895.5 s降到了0.17 s和0.4 s。從結(jié)果可以看出直接法隨著幾何復(fù)雜度增加，計算耗時急劇增加，而本文方法耗時并未顯著增加。在對CHN-T1模型建模過程中直接法采用單核已無法完成，通過120核的并行生成判斷花了40 056 s（11.13 h），而采用本文算法即使采用1核計算判斷過程也只耗時4.5 s，加文件讀寫，前處理總耗時約20 s。測試結(jié)果表明本文算法極大提高了LBM針對復(fù)雜幾何的前處理效率。

表1 本文算法和直接法的測試結(jié)果Table 1 Testing results for the present algorithm and the direct method

5 應(yīng)用案例

針對不同算例，基于本文的前處理算法開展網(wǎng)格生成，并采用SWLBM[17]軟件開展了相應(yīng)流場計算。圖8為針對圓球模型開展以直徑為特征長度Re=3.3×106工況下計算獲得的Q渦結(jié)構(gòu)與軸向切面速度場，可以發(fā)現(xiàn)基于本文方法構(gòu)建的網(wǎng)格可以獲得圓球繞流的復(fù)雜渦系結(jié)構(gòu)。

圖8 圓球繞流Q渦結(jié)構(gòu)與速度云圖Fig. 8 Vortical structures based on Q-criterion and the velocity contour for flow around a sphere

圖9為本文算法所建模型計算獲得的在8°攻角下NACA0012翼型Re= 1 000工況下的繞流場速度云圖。圖10為本文計算所得翼型表面壓力系數(shù)與文獻(xiàn)[18]的結(jié)果對比，由圖可見兩者基本一致，表明基于本文所建立的前處理算法獲得的網(wǎng)格模擬,可以獲得同樣準(zhǔn)確的結(jié)果。

圖9 NACA0012翼型繞流速度云圖Fig. 9 Velocity contour for flow around the NACA0012 airfoil

圖10 NACA0012表面壓力系數(shù)Fig. 10 Pressure coefficient along NACA0012 airfoil

圖11為基于本文算法前處理獲得CHN-T1標(biāo)模2 mm分辨率笛卡爾網(wǎng)格開展4°攻角Re= 3.3×106、Ma= 0.2工況下，流場計算獲得的不同視角瞬態(tài)Q渦結(jié)構(gòu)。圖12為同狀態(tài)下不同特征截面上的速度場云圖?？梢钥吹交诒疚那疤幚矸椒ǐ@得的網(wǎng)格開展的流場計算，可以獲得從機(jī)翼脫落的渦結(jié)構(gòu)，垂尾形成的卡門渦街也清晰可見。

圖11 CHN-T1模型不同視角的Q渦結(jié)構(gòu)Fig. 11 Vortical structures based on Q-criterion for flow around the CHN-T1 model from different perspectives

圖12 CHN-T1模型不同特征面速度場Fig. 12 Velocity contours in different characteristic planes of the CHN-T1 model

以上三個算例結(jié)果，驗證了本文前處理方法在LBM計算中的適用性。

6 結(jié)束語

本文針對LBM計算中流場網(wǎng)格類型判斷需求，建立了一套基于STL幾何文件的快速前處理方法。與常用的直接法比較，克服了其隨著幾何復(fù)雜度增加耗時急劇增加問題，并且對復(fù)雜幾何具有同樣高效。該算法從模型的面網(wǎng)格三角形出發(fā)，找出面網(wǎng)格對應(yīng)的邊界格點，進(jìn)一步由邊界格點完成對模型的重構(gòu)。基于本算法的LBM前處理技術(shù)已經(jīng)在SWLBM軟件中獲得很好的應(yīng)用。在之后的工作中，還將從并行以及針對更為特殊的模型等方面考慮，進(jìn)一步擴(kuò)大算法的適用性，提高算法的模型重建效率。