李 平，麻鐵昌，許香照，馬天寶

（北京理工大學(xué)機(jī)電學(xué)院，北京 100081）

爆炸力學(xué)是武器裝備研發(fā)，工業(yè)安全防護(hù)等問題的理論基礎(chǔ)。爆炸的特點(diǎn)為在極短的時(shí)間內(nèi)，爆炸場中各個(gè)物理量發(fā)生巨大的變化。爆炸過程的研究十分復(fù)雜，通常伴隨有化學(xué)變化、物理變化和核變化等。因此在理論研究和實(shí)驗(yàn)研究中都存在諸多困難。數(shù)值計(jì)算是研究爆炸與沖擊問題的重要手段之一。計(jì)算爆炸力學(xué)是以計(jì)算機(jī)為工具，探索爆炸力學(xué)的規(guī)律，豐富力學(xué)數(shù)據(jù)，解決爆炸力學(xué)問題。隨著計(jì)算機(jī)科學(xué)的飛速發(fā)展，數(shù)值模擬以其精度高，經(jīng)濟(jì)成本低，揭示規(guī)律清晰完整等諸多優(yōu)點(diǎn)越來越成為人們的研究熱點(diǎn)。爆炸力學(xué)的數(shù)值模擬按照其采用的坐標(biāo)形式可以分為歐拉方法和拉格朗日方法?；跉W拉坐標(biāo)的有限差分方法是爆炸與沖擊問題數(shù)值模擬的常用方法[1]，其特點(diǎn)為在固定的階梯型網(wǎng)格上進(jìn)行有限差分計(jì)算，網(wǎng)格位置和形狀不會(huì)隨著計(jì)算而改變，克服了拉格朗日方法中可能出現(xiàn)的網(wǎng)格畸變問題，因此能夠很好的模擬多物質(zhì)的大變形問題[2-3]。本文涉及的有限差分計(jì)算中使用的網(wǎng)格數(shù)據(jù)是一種階梯型笛卡爾網(wǎng)格，其特點(diǎn)是使用正六面體對整個(gè)區(qū)域進(jìn)行離散，不同物質(zhì)被離散為不同屬性的單元，不同屬性單元構(gòu)成的不同區(qū)域之間的交界面為階梯形[4-5]。相比于傳統(tǒng)貼體網(wǎng)格，階梯型有限差分網(wǎng)格的生成速度快，無需手動(dòng)劃分區(qū)域，生成過程自動(dòng)化程度高。階梯型有限差分網(wǎng)格是一種結(jié)構(gòu)網(wǎng)格，每個(gè)單元之間的連接關(guān)系固定并不隨著計(jì)算而變化，可以使用像所有結(jié)構(gòu)網(wǎng)格一樣的IJK 下標(biāo)來訪問網(wǎng)格單元和節(jié)點(diǎn)。它適用于所有的有限差分離散方法。

階梯型有限差分網(wǎng)格生成方法主要包含兩種：一種是射線穿透法，這種方法即根據(jù)網(wǎng)格步長作出一系列射線直接與實(shí)體模型進(jìn)行求交計(jì)算[6]，這種方法的優(yōu)點(diǎn)是思路清晰，步驟簡單；但是由于空間中的射線與幾何模型的求交計(jì)算會(huì)消耗大量的時(shí)間，因此算法運(yùn)行效率不高。切片法的思路借鑒了快速成型加工領(lǐng)域的方法[7-8]，用一系列互相平行的平面截取幾何模型，得到一系列二維平面中的封閉輪廓線，進(jìn)而將三維空間中的幾何問題轉(zhuǎn)換成二維平面上的問題；此方法省去了部分幾何判斷、查詢占用的時(shí)間，網(wǎng)格生成效率相對較高。但是此種算法執(zhí)行過程中網(wǎng)格線的生成和介質(zhì)屬性映射兩個(gè)主要步驟耦合在一起，算法執(zhí)行的穩(wěn)定性和并行能力都有待提高。在三維有限差分計(jì)算中，網(wǎng)格尺寸越小，網(wǎng)格構(gòu)成的計(jì)算空間與實(shí)際物理空間越接近。因此為了保證數(shù)值模擬的精度，進(jìn)行三維有限差分計(jì)算通常需要至少千萬量級的網(wǎng)格[9]，如何高效的生成數(shù)量如此巨大的數(shù)值網(wǎng)格一直是研究熱點(diǎn)之一。MacGrillivray[10]使用高效的數(shù)據(jù)存儲(chǔ)方法和高效的空間線面相交判別方法實(shí)現(xiàn)了萬億量級的網(wǎng)格生成，網(wǎng)格生成過程內(nèi)存管理方法合理，使整個(gè)網(wǎng)格生成過程可以在小型工作站中完成。Berens 等[11]提出了非均勻笛卡爾網(wǎng)格的生成方法，詳細(xì)描述了根據(jù)計(jì)算需求而對網(wǎng)格尺寸的選取方法。以上都是串行的網(wǎng)格生成方法。并行計(jì)算是提高網(wǎng)格生成效率的最常用的方法。在計(jì)算集群中運(yùn)行的程序最常使用的是利用MPI(message passing interface)庫進(jìn)行并行計(jì)算，在小型工作站中常用的是OpenMP 庫實(shí)現(xiàn)共享內(nèi)存的并行計(jì)算，Ning 等[12]使用多線程技術(shù)實(shí)現(xiàn)了并行階梯型有限差分網(wǎng)格生成，但是其采用的是切片法生成網(wǎng)格，單個(gè)幾何體分別由單個(gè)線程執(zhí)行計(jì)算，并行化程度較低，對大型的網(wǎng)格數(shù)據(jù)生成效率不高。Ishida 等[13]提出采用OpenMP 并行生成自適應(yīng)笛卡爾網(wǎng)格。Foteinos 等[14]將并行計(jì)算應(yīng)用在“圖到網(wǎng)格”的網(wǎng)格生成算法中，采用分布式集群計(jì)算機(jī)對算法進(jìn)行加速。

最近，GPU （graphic processing unit）并行技術(shù)吸引了各個(gè)領(lǐng)域研究者的注意，不僅僅被應(yīng)用在圖形圖像相關(guān)算法，還被廣泛應(yīng)用于各種通用計(jì)算中。相對于CPU （central process unit）并行，GPU 對于大規(guī)模密集型數(shù)據(jù)的并行計(jì)算擁有巨大的優(yōu)勢。許多學(xué)者也對將GPU 應(yīng)用到網(wǎng)格生成領(lǐng)域做了相關(guān)研究。Qi 等[15]利用GPU 并行實(shí)現(xiàn)了二維Delaunay 非結(jié)構(gòu)網(wǎng)格的生成。Park 等[16]提出了基于GPU 生成3D 自適應(yīng)笛卡爾網(wǎng)格的方法。Schwarz 等[17]提出了一種方法使用GPU 來生成八叉樹結(jié)構(gòu)，并將其應(yīng)用于計(jì)算機(jī)圖形學(xué)中的像素化顯示。但是他們使用的基于三角形并行化的方法生成的單元數(shù)量較少，不適用于進(jìn)行有限差分計(jì)算。根據(jù)文獻(xiàn)調(diào)研情況，目前運(yùn)用GPU 并行技術(shù)生成三維笛卡爾網(wǎng)格的研究很少。本文中提出的并行算法能夠應(yīng)用GPU 并行技術(shù)在短時(shí)間內(nèi)生成大規(guī)模的階梯型有限差分網(wǎng)格，并且對多個(gè)擁有復(fù)雜外形的幾何體的大規(guī)模網(wǎng)格生成過程有很高的效率和準(zhǔn)確性。

本文中，選用STL （stereolithography）文件來存儲(chǔ)幾何實(shí)體信息，并用其作為算法的輸入文件。在STL 文件中，幾何實(shí)體被離散為三角面片的形式，所有三角面片的三維坐標(biāo)和法向量都存貯其中[18]。在一些切片算法中，STL 文件需要過濾冗余幾何信息并重構(gòu)出拓?fù)湫畔ⅰτ谏渚€穿透法，STL 文件中的信息不再需要任何的處理，STL 文件中的三角形集合T 可以直接作為算法的輸入。階梯型有限差分網(wǎng)格生成算法的核心部分就是由T 作為輸入生成整個(gè)計(jì)算域大小的數(shù)據(jù)場F，F(xiàn) 中每一坐標(biāo)點(diǎn)的值為介質(zhì)標(biāo)志。

1 射線穿透法

1.1 計(jì)算域設(shè)置

（1）計(jì)算域的邊界線和計(jì)算域內(nèi)各個(gè)幾何體的AABB 包圍盒的邊線同時(shí)也是網(wǎng)格線。三維空間中幾何對象的AABB 包圍盒為包含該對象，且邊平行于坐標(biāo)軸的最小六面體。

（2）計(jì)算域內(nèi)同一維度方向上的網(wǎng)格尺寸s 相同。

條件1 和條件2 同時(shí)限定了整個(gè)計(jì)算域中網(wǎng)格尺寸大小相同。同時(shí)，條件1 消除了在幾何體AABB 包圍盒邊緣存在的狹窄網(wǎng)格單元。為了滿足上述兩個(gè)條件，網(wǎng)格單元尺寸需要根據(jù)輸入的尺寸進(jìn)行優(yōu)化。如下列方程分別為X、Y 和Z 方向上的網(wǎng)格尺寸函數(shù)：

式中：sopt為待優(yōu)化的網(wǎng)格尺寸，nx、ny和nz分別為X、Y 和Z 方向上的網(wǎng)格單元數(shù)，ΔXi為計(jì)算域被幾何體AABB 包圍盒邊界線分割出的子計(jì)算區(qū)域。

圖1 三維計(jì)算域，內(nèi)部幾何體和幾何體局部圖Fig. 1 Three-dimensional computational domain, the geometries in the domain and the details of the geometry

優(yōu)化方程(1)～(3)中的sopt，計(jì)算使得3 個(gè)方程同時(shí)取得最小值同時(shí)又不大于初始的網(wǎng)格尺寸s 的sopt，即為優(yōu)化后的網(wǎng)格單元尺寸。

獲得優(yōu)化的網(wǎng)格單元尺寸后，可以按照下列公式生成整個(gè)計(jì)算域中網(wǎng)格線：

1.2 射線與幾何體表面的相交判斷

對于任意一條射線穿過計(jì)算區(qū)域，可能與一些三角面片相交，也可能不相交。如圖2 中所示，由投影平面出發(fā)的射線穿過計(jì)算區(qū)域與幾何體相交，射線簇與幾何體的相交計(jì)算的本質(zhì)是空間中射線r 與三角面片的相交計(jì)算，V1、V2和V3為三角面片3 個(gè)頂點(diǎn)。為了計(jì)算射線與三角面片的交點(diǎn)，首先要判斷射線與哪些面片相交。為了減少機(jī)器誤差，這里采用一種沒有除法計(jì)算的判斷方法。對于一條射線和一個(gè)三角面片，通過下列方程計(jì)算射線與三角面片第i 條邊的相對位置Oi：

圖2 射線穿透法原理（以XY 平面為射線投射平面，Z 軸方向?yàn)樯渚€方向）Fig. 2 Principle of ray casting method (set the XY plane as the projection plane, and the Z dimension as the ray direction)

Oi表示射線與三角面片的一條邊的相對位置，如果Oi＞0，則射線r 在邊V1V2的左側(cè)，如果Oi＜0，則射線r 在邊V1V2的右側(cè)，Oi=0 時(shí)，射線與該邊相交。分別計(jì)算三角面片三條邊與射線的Oi值，如果Oi全部為正或全部為負(fù)值，則射線在3 條邊的同側(cè)，表明射線穿過了三角形，否則射線與三角形不相交。

通過一條穿過幾何體的射線將會(huì)計(jì)算出多個(gè)交點(diǎn)，將這些交點(diǎn)按照射線方向排序，得到交點(diǎn)集合P。在三維空間中，直線與三角面片的位置關(guān)系可分為7 類：穿過、錯(cuò)過、穿點(diǎn)、穿邊、平行、共面不相交和共面相交。很顯然，穿過和錯(cuò)過可以由上述方法判斷并計(jì)算交點(diǎn)。穿點(diǎn)是直線穿過三角面片任意一頂點(diǎn)，同時(shí)由于三角面片組成的是閉合幾何體，此直線一定也穿過其他三角面片的頂點(diǎn)。在這種情況下，在判別式中，O1、O2、O3中有兩個(gè)為0，另外一個(gè)不為0。交點(diǎn)計(jì)算公式仍然成立，但是會(huì)計(jì)算出多個(gè)相同的交點(diǎn)，在排序后所得的集合P 中，相同的交點(diǎn)只保留一個(gè)。對于穿邊的情況，依據(jù)相同的原理，O1、O2、O3中有一個(gè)為0，另外兩個(gè)同號。交點(diǎn)計(jì)算公式依然成立，在這個(gè)位置會(huì)有兩個(gè)交點(diǎn)，在P 中保留其中的一個(gè)。在直線與三角面片平行時(shí)，O1、O2、O3一定不會(huì)符號相同，因此會(huì)判定為錯(cuò)過，不計(jì)算交點(diǎn)。當(dāng)直線與三角面片共面不相交時(shí)，O1、O2、O3全部等于0，會(huì)判定為不相交，不計(jì)算交點(diǎn)。當(dāng)直線與三角面片共面并穿過面片時(shí)，說明這個(gè)面片為幾何物體的邊界，邊界上不會(huì)有網(wǎng)格點(diǎn)，因此判別式可以正確將其判斷為不相交。

綜上所述，在穿點(diǎn)與穿邊的情況中，會(huì)出現(xiàn)一項(xiàng)或兩項(xiàng)Oi為0，這種情況同樣應(yīng)視為相交，計(jì)算出的交點(diǎn)可能會(huì)在其相鄰三角面片中再次出現(xiàn)，使用上述原理僅保留一個(gè)交點(diǎn)即可。

2 GPU 并行網(wǎng)格生成策略

三維空間中幾何圖元的搜索、查找和相交計(jì)算是射線穿透法的主要耗時(shí)部分。同時(shí)，這些計(jì)算都可以轉(zhuǎn)化為某一簡單計(jì)算流程的疊加。因此，將這些耗時(shí)巨大的計(jì)算轉(zhuǎn)移到高度并行化的GPU 中并行執(zhí)行可以大大提高算法效率。GPU 采用單程序多數(shù)據(jù)（single program/multiple data，SPMD）模型，其目的就是通過其內(nèi)部大量的線程、線程束、線程塊和線程網(wǎng)格等并行層級來執(zhí)行大量的、比較簡單的計(jì)算任務(wù)。GPU 并行算法的效率主要由2 個(gè)部分決定：（1）單個(gè)線程中計(jì)算時(shí)間；（2）主機(jī)內(nèi)存與GPU 內(nèi)存數(shù)據(jù)交換時(shí)間。

下面，針對上述兩部分分別優(yōu)化并行算法。

2.1 幾何數(shù)據(jù)傳輸優(yōu)化

下面以一個(gè)單獨(dú)幾何體G 為例說明數(shù)據(jù)傳輸過程。作為輸入，STL 中的三角面片數(shù)據(jù)被存儲(chǔ)到內(nèi)存序列中，以數(shù)組T 表示。幾何體G 的AABB 包圍盒被一系列水平平面二次劃分為M 個(gè)包圍盒子區(qū)域。隨后，將三角形序列T 按照其在空間中所屬的包圍盒子區(qū)域分解為M 個(gè)子三角形序列，將一個(gè)包圍盒子區(qū)域的信息SubAABB 和與其對應(yīng)的子三角形序列Tsub作為一個(gè)批次（batch），則輸入數(shù)據(jù)被劃分為M 個(gè)批次，每個(gè)批次都包含了處理每個(gè)批次中所有射線計(jì)算所需的全部數(shù)據(jù)。隨后，將M 個(gè)批次數(shù)據(jù)逐個(gè)傳入GPU 內(nèi)存中執(zhí)行并行計(jì)算。每一個(gè)線程執(zhí)行計(jì)算后的輸出數(shù)據(jù)為一條射線r 方向上的網(wǎng)格個(gè)數(shù)序列Cr，Cr中元素的數(shù)量表示射線r 方向上的總網(wǎng)格數(shù)被幾何體分割出的段數(shù)，每段的網(wǎng)格數(shù)由Cr中一個(gè)元素表示。幾何體G 的AABB 包圍盒內(nèi)全部射線對應(yīng)的網(wǎng)格個(gè)數(shù)序列的集合可以表示為C={Cr1,Cr2,···,Cri,···,Crn}，C 中包含的數(shù)據(jù)即可表示幾何體G 的網(wǎng)格數(shù)據(jù)。對于計(jì)算域中的每一個(gè)幾何體，重復(fù)上述過程，將每個(gè)幾何體計(jì)算所得的網(wǎng)格個(gè)數(shù)序列集合C 合并在一起，即可得到整個(gè)計(jì)算域中的階梯型有限差分網(wǎng)格劃分結(jié)果。上述數(shù)據(jù)轉(zhuǎn)換與計(jì)算過程如圖3 所示。圖中左上部分表示待計(jì)算區(qū)域及計(jì)算區(qū)域中的幾何體，其中幾何體以三角面片形式被表示。右上圖表示幾何體的AABB 包圍盒構(gòu)成的子計(jì)算區(qū)域。圖中下半部分表示幾何數(shù)據(jù)及邊界條件等被組織為數(shù)據(jù)批次，并逐批次傳如GPU 進(jìn)行計(jì)算的并行計(jì)算流程。

圖3 三維幾何數(shù)據(jù)轉(zhuǎn)換與并行計(jì)算過程圖Fig. 3 Diagram of three-dimensional geometric data conversion and parallel computing process

2.2 線程內(nèi)射線生成

GPU 并行計(jì)算開始執(zhí)行之后，大量的線程會(huì)同時(shí)執(zhí)行相同的指令，這些指令的設(shè)計(jì)直接影響并行算法的并行執(zhí)行效率。對于域中的每個(gè)幾何體，網(wǎng)格線和射線的分布可以通過AABB 包圍盒及包圍盒子區(qū)域的對角點(diǎn)坐標(biāo)(x1, y1)和(x2, y2)與網(wǎng)格尺寸sx和sy來確定。當(dāng)GPU 執(zhí)行某個(gè)計(jì)算任務(wù)時(shí)，通常需要將輸入數(shù)據(jù)從主機(jī)內(nèi)存復(fù)制到GPU 中的設(shè)備內(nèi)存。為了避免主機(jī)和設(shè)備之間數(shù)據(jù)傳輸?shù)臅r(shí)間消耗和GPU 內(nèi)存中的容量消耗，可以在GPU 中的每個(gè)線程中自動(dòng)生成網(wǎng)格線和射線的起點(diǎn)坐標(biāo)。

假設(shè)幾何模型G 的AABB 包圍盒可以表示為兩個(gè)對角點(diǎn)坐標(biāo)(x1, y1)和(x2, y2)。X 與Y 方向上的網(wǎng)格尺寸可以分別表示為sx和sy。在GPU 的每個(gè)線程中，當(dāng)前線程開始計(jì)算時(shí)，可以通過下列公式生成射線的起點(diǎn)坐標(biāo)(xr, yr)：

式中：Nindex和nx分別表示GPU 中當(dāng)前線程的索引序號和當(dāng)前包圍盒子區(qū)域中X 坐標(biāo)方向的網(wǎng)格數(shù)。Nindex和nx分別可以由下式計(jì)算：

生成射線起點(diǎn)坐標(biāo)以后，應(yīng)用1.2 節(jié)中所述方法判斷射線與哪些三角形相交并計(jì)算出交點(diǎn)坐標(biāo)。將各個(gè)交點(diǎn)Z 坐標(biāo)連同計(jì)算域Z 方向的最大、最小值按照由小到大的次序排列，可得到Z 坐標(biāo)序列(0＜i＜2n)。由于計(jì)算域中幾何體都是封閉的，交點(diǎn)個(gè)數(shù)必定為偶數(shù)，以2n 表示。則對于射線r 方向上所有的網(wǎng)格單元，第i 個(gè)網(wǎng)格單元的屬性可以按照下列公式賦值：

式中：sz為計(jì)算域Z 方向上的網(wǎng)格尺寸。這一步驟通常稱為介質(zhì)屬性映射。本節(jié)所述計(jì)算全部由單個(gè)線程獨(dú)立完成，計(jì)算流程如圖4 所示。

圖4 單個(gè)GPU 線程內(nèi)數(shù)據(jù)計(jì)算流程Fig. 4 Data computing process in a single GPU thread

3 數(shù)值試驗(yàn)結(jié)果和討論

應(yīng)用文中提出的基于射線穿透法的GPU 并行網(wǎng)格生成方法，使用Visual C++和Nvida CUDA[19]編制基于GPU 的并行階梯型有限差分網(wǎng)格生成程序并對程序進(jìn)行性能測試。所用測試服務(wù)器配置如下：CPU 型號，Intel Xeon E5-2650 v2（2.6 GHz）；GPU 型號，Nvidia Quadro K2200；顯存容量，4 GB。

3.1 網(wǎng)格生成結(jié)果與可視化

如圖5(a)所示橋梁模型包含3 種材質(zhì)，由多個(gè)幾何體組合而成。模型總計(jì)包含24 184 個(gè)三角面片，被劃分為1.15×109個(gè)網(wǎng)格單元?？傆?jì)用時(shí)21.45 s。圖5(b)為網(wǎng)格圖，由放大后的局部網(wǎng)格圖可以清晰地看出，三角面片構(gòu)成的幾何模型被轉(zhuǎn)變?yōu)橛刹煌馁|(zhì)的六面體單元構(gòu)成的階梯型有限差分網(wǎng)格。

為了驗(yàn)證并行算法的網(wǎng)格生成效率，將上述模型分別離散為不同網(wǎng)格規(guī)模的階梯型有限差分網(wǎng)格，并對比網(wǎng)格生成時(shí)間與傳統(tǒng)串行射線穿透法的網(wǎng)格生成時(shí)間。如圖6 所示，橋梁模型分別被離散為1.44×108、1.15×109和9.2×109個(gè)網(wǎng)格單元?？梢钥闯?，對同一計(jì)算模型應(yīng)用本文提出的GPU 并行算法的網(wǎng)格生成效率遠(yuǎn)高于傳統(tǒng)串行CPU 算法的執(zhí)行效率。傳統(tǒng)串行算法在生成網(wǎng)格規(guī)模達(dá)到1×1010數(shù)量級的時(shí)候耗費(fèi)的時(shí)間超過2 000 s，這無疑大大影響了建模與計(jì)算時(shí)的靈活性，在很多需要多次調(diào)整計(jì)算域參數(shù)的數(shù)值模擬中是難以接受的。另外，表示這種數(shù)量級的網(wǎng)格數(shù)據(jù)通?？梢赃_(dá)到4 GB 以上，本文中提出的分批次的數(shù)據(jù)處理方法使得算法能夠處理的數(shù)據(jù)規(guī)模不依賴于GPU 內(nèi)存大小，使程序在常見的4 GB 顯存容量的GPU 中可以高效執(zhí)行。CPU 的硬件發(fā)展已趨于穩(wěn)定，短時(shí)間內(nèi)難以有巨大的提升。本文中提出的并行算法可以解決這些瓶頸，使得在擁有一顆普通GPU 的PC 機(jī)上就足以進(jìn)行1×1010數(shù)量級的網(wǎng)格生成。對于不同的初始幾何模型，網(wǎng)格生成算法的效率應(yīng)該不受模型的復(fù)雜程度的影響。表1 中展示了3 種不同的初始幾何模型分別包含19 202、78 354 和95 062 個(gè)三角面片，分別生成相同規(guī)模（1×109個(gè)網(wǎng)格單元）的三維笛卡爾網(wǎng)格。表中分別給出了CPU 串行算法和GPU 并行算法的生成時(shí)間。可以看出，對于復(fù)雜程度不同的幾何模型生成相同規(guī)模的網(wǎng)格，本文提出的并行算法的效率是傳統(tǒng)串行算法的8～11 倍。

圖5 橋梁模型及其階梯型有限差分網(wǎng)格生成結(jié)果圖及細(xì)節(jié)放大圖Fig. 5 A bridge model, the finite difference mesh generated and the details of the mesh

圖6 本文提出的并行算法與傳統(tǒng)算法的網(wǎng)格生成時(shí)間比較折線圖Fig. 6 Mesh generation time comparison curve between the proposed parallel algorithm and the traditional algorithm

表1 不同模型生成相同數(shù)量網(wǎng)格單元(1×109)的執(zhí)行時(shí)間Table 1 Generation times of different models with the cell number of 1×109

綜上所述，本文中提出的并行算法的網(wǎng)格生成效率是傳統(tǒng)串行算法的8～11 倍。隨著網(wǎng)格數(shù)量級和模型復(fù)雜度的增高，并行算法節(jié)省的網(wǎng)格生成時(shí)間越來越多。

3.2 數(shù)值模擬結(jié)果

工廠廠房內(nèi)的爆炸是工業(yè)生產(chǎn)中危害巨大的一種安全事故。廠房合理的結(jié)構(gòu)設(shè)計(jì)是降低廠房爆炸造成的損失的有效手段。在這一節(jié)，通過文中提出的并行階梯型有限差分網(wǎng)格生成方法，對某工廠廠房進(jìn)行網(wǎng)格劃分，并對廠房內(nèi)爆炸問題進(jìn)行數(shù)值模擬研究。如圖7 所示，計(jì)算域中包括墻壁、立柱、屋頂、炸藥和空氣5 種介質(zhì)。廠房為一座兩層建筑，上下兩層有樓梯連接，64 kg TNT 炸藥位于廠房一層中間位置。為了應(yīng)用有限差分法進(jìn)行數(shù)值模擬，將計(jì)算域離散為三維階梯型有限差分網(wǎng)格，如圖8 所示。為了清晰顯示網(wǎng)格生成結(jié)果，空氣網(wǎng)格被隱藏掉，圖8 中的炸藥，墻壁、立柱和屋頂4 種介質(zhì)分別以4 種顏色顯示。網(wǎng)格生成單元總數(shù)為1.5×108個(gè)。并行生成消耗總時(shí)間為2.47 s。

圖8 某廠房三維階梯型有限差分網(wǎng)格Fig. 8 Three-dimensional finite difference mesh of a factory model

圖7 某廠房三維幾何模型Fig. 7 A three-dimensional factory model

應(yīng)用北京理工大學(xué)爆炸與科學(xué)國家重點(diǎn)實(shí)驗(yàn)室自主開發(fā)的三維多介質(zhì)流體動(dòng)力學(xué)仿真軟件PMMIC-3D[20]對上述網(wǎng)格數(shù)據(jù)進(jìn)行數(shù)值模擬，得到計(jì)算結(jié)果如圖9 所示。圖9（a）中左側(cè)為爆炸產(chǎn)生后3 個(gè)有代表性的時(shí)刻的三維可視化結(jié)果，右側(cè)為與左側(cè)對應(yīng)的時(shí)刻的二維剖面可視化結(jié)果。從圖9（a）可以清晰地看到爆炸在廠房一層發(fā)生后，形成沖擊波（圖中以白色表示）并膨脹擴(kuò)散到二層的過程。圖9（b）為廠房二層2 維剖面圖，以壓力的變化為可視化屬性顯示了在4 個(gè)典型的時(shí)刻沖擊波在二層傳播過程。

圖9 廠房爆炸數(shù)值模擬結(jié)果圖Fig. 9 Numerical simulation results of factory explosion

為了更進(jìn)一步分析數(shù)值模擬的準(zhǔn)確性，選取了廠房二層距樓梯口10 m 處的點(diǎn)作為關(guān)鍵點(diǎn)M 進(jìn)行超壓測試。記錄點(diǎn)M 超壓隨時(shí)間的變化如圖10 所示?？梢钥吹剑c(diǎn)M 處初始超壓為0，隨后沖擊波傳入二樓并傳到點(diǎn)M 處，壓力值開始上升達(dá)到第一個(gè)峰值。沖擊波傳播到墻壁后反射再次到達(dá)點(diǎn)M，壓力值達(dá)到第二個(gè)峰值。隨后，隨著不斷的反射，M 點(diǎn)超壓值逐漸減小直至衰減為0?？梢钥闯鰯?shù)值模擬與沖擊波傳播反射理論和實(shí)際經(jīng)驗(yàn)相吻合，說明計(jì)算域的網(wǎng)格剖分能夠滿足大規(guī)模有限差分計(jì)算的需要。

圖10 廠房二層距離樓梯口10 m 處關(guān)鍵點(diǎn)的超壓變化曲線Fig. 10 Change of overpressure with time at a key point which is 10 m from stairway entrance on the second floor

4 結(jié) 論

本文提出了一種基于傳統(tǒng)射線穿透法的GPU 并行階梯型有限差分網(wǎng)格生成方法。在這種并行方法中，提出了一種分批次的數(shù)據(jù)傳輸策略，使得算法能夠處理的數(shù)據(jù)規(guī)模不依賴于GPU 內(nèi)存大小，打破了硬件對網(wǎng)格劃分規(guī)模的限制，平衡了數(shù)據(jù)傳輸效率和網(wǎng)格生成規(guī)模之間的關(guān)系。為了減少數(shù)據(jù)傳輸量，本文提出的并行算法可以由GPU 線程相互獨(dú)立的生成射線起點(diǎn)坐標(biāo)，射線相交計(jì)算在GPU 的每個(gè)線程中獨(dú)自計(jì)算，進(jìn)一步提高了并行算法的執(zhí)行效率，通過數(shù)值試驗(yàn)的對比可以看出，并行算法的執(zhí)行效率是傳統(tǒng)射線穿透法執(zhí)行效率的8～11 倍，并且隨著計(jì)算規(guī)模的提升，并行算法的加速比有上升趨勢。最后，通過有限差分計(jì)算實(shí)例驗(yàn)證了應(yīng)用并行算法生成的階梯型有限差分網(wǎng)格能夠滿足基于有限差分的數(shù)值模擬需求，得到了與理論和實(shí)驗(yàn)一致的數(shù)值模擬結(jié)果。