焦子龍，姜利祥，黃建國(guó)，孫繼鵬，朱云飛

（北京衛(wèi)星環(huán)境工程研究所 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100094）

利用并行八叉樹(shù)劃分技術(shù)開(kāi)展復(fù)雜構(gòu)型航天器放氣污染分析研究

焦子龍，姜利祥，黃建國(guó)，孫繼鵬，朱云飛

（北京衛(wèi)星環(huán)境工程研究所 可靠性與環(huán)境工程技術(shù)重點(diǎn)實(shí)驗(yàn)室，北京 100094）

放氣污染對(duì)航天器造成的有害影響需要定量評(píng)估。放氣產(chǎn)物分子的輸運(yùn)過(guò)程可以采用 Monte-Carlo方法模擬，并用射線表示放氣分子輸運(yùn)軌跡。如此一來(lái)，即轉(zhuǎn)換為射線追蹤問(wèn)題。由于航天器構(gòu)型復(fù)雜，射線追蹤的計(jì)算很耗時(shí)，所以有必要研究其加速算法。文章將復(fù)雜構(gòu)型航天器表面用三角形非結(jié)構(gòu)網(wǎng)格表示，利用八叉樹(shù)劃分技術(shù)開(kāi)展分子運(yùn)動(dòng)軌跡與表面相交的加速模擬計(jì)算。同時(shí)，在基于共享內(nèi)存多核計(jì)算平臺(tái)，利用OpenMP API軟件實(shí)現(xiàn)了計(jì)算的并行化。該方法在對(duì)空間站的實(shí)例計(jì)算中顯著提高了運(yùn)算速度，并對(duì)空間站的污染情況進(jìn)行了初步評(píng)價(jià)分析。

航天器；放氣污染；模擬計(jì)算；復(fù)雜構(gòu)型；八叉樹(shù)；并行計(jì)算

0 引言

在航天器結(jié)構(gòu)中，廣泛應(yīng)用碳纖維/環(huán)氧復(fù)合材料、熱控涂層材料、隔熱材料、膠粘劑等[1-3]。航天器在軌運(yùn)行過(guò)程中，這些材料在真空、高低溫、粒子輻射等環(huán)境中將發(fā)生質(zhì)量損失，統(tǒng)稱為放氣[4]。這不僅影響材料本身的性能，還會(huì)對(duì)航天器造成污染。例如歐洲的Rosetta彗星探測(cè)器在利用四極質(zhì)譜計(jì)探測(cè)彗星周?chē)鷼怏w成分時(shí)，發(fā)現(xiàn)探測(cè)器材料放氣對(duì)探測(cè)結(jié)果產(chǎn)生了嚴(yán)重影響[5]。Cassini探測(cè)器在軌飛行中材料放氣使其窄視角相機(jī)受到嚴(yán)重污染，造成圖像出現(xiàn)光暈現(xiàn)象，不得不進(jìn)行了長(zhǎng)時(shí)間的加熱除氣[6]。因此，有必要對(duì)航天器用材料提出嚴(yán)格的放氣控制要求。

除需對(duì)材料放氣進(jìn)行嚴(yán)格控制外，還需要定量評(píng)估放氣帶來(lái)的影響。國(guó)內(nèi)外對(duì)放氣污染評(píng)估分析方法進(jìn)行了廣泛研究。例如，針對(duì)國(guó)際空間站的污染分析與評(píng)估，ESA和 NASA分別開(kāi)發(fā)了COMOVA和NASAN污染分析軟件[7]。這2款軟件均采用Monte-Carlo方法，并用射線來(lái)表征分子的輸運(yùn)過(guò)程，包括分子從放氣源表面出射，在航天器周?chē)\(yùn)動(dòng)、碰撞，在敏感表面被吸附等過(guò)程。

航天器構(gòu)型復(fù)雜，射線追蹤的計(jì)算很耗時(shí)，因此有必要研究其加速算法。本文將計(jì)算機(jī)圖形學(xué)中常用到的射線追蹤加速方法即并行八叉樹(shù)劃分技術(shù)引入復(fù)雜構(gòu)型航天器的放氣污染分析中，并在OpenMP API軟件平臺(tái)上開(kāi)展了模擬分析，獲得了很好的加速效果。

1 放氣污染計(jì)算

1.1 蒙特卡羅方法

粒子在軌道環(huán)境中傳輸分為直接流和背景散射返回流，對(duì)這些粒子流可以采用試驗(yàn)粒子蒙特卡羅方法（Test Particle Monte-Carlo, TPMC）進(jìn)行模擬。TPMC方法基于硬球碰撞（二體彈性碰撞）模型，采用統(tǒng)計(jì)方法對(duì)大量放氣分子的運(yùn)動(dòng)、碰撞和軌跡進(jìn)行模擬計(jì)算，以獲得相應(yīng)的散射分布和返回流比率，具體計(jì)算方法和步驟可參見(jiàn)文獻(xiàn)[8-10]。

1.2 八叉樹(shù)劃分算法

八叉樹(shù)（octa-tree）算法是一種常用的空間劃分技術(shù)，具有數(shù)據(jù)結(jié)構(gòu)簡(jiǎn)單、集合運(yùn)算高效、節(jié)點(diǎn)遍歷運(yùn)算易于實(shí)現(xiàn)等優(yōu)點(diǎn)，是射線追蹤普遍采用的加速方法，如圖1所示。

圖1 八叉樹(shù)空間劃分示意圖Fig.1 Sketch of octa-tree subdivision

其基本過(guò)程是：構(gòu)建包含所有空間對(duì)象的長(zhǎng)方體包圍盒作為根節(jié)點(diǎn)；再對(duì)該包圍盒進(jìn)行細(xì)分，即沿3個(gè)坐標(biāo)軸方向?qū)⑵浞譃?個(gè)大小相等的子立方體；接著檢查每個(gè)子立方體，如果子立方體內(nèi)空間對(duì)象的個(gè)數(shù)不大于給定的閾值，則停止分解該子立方體，否則，需要對(duì)該子立方體繼續(xù)進(jìn)行8等分，直到內(nèi)部空間對(duì)象數(shù)低于給定閾值。

采用八叉樹(shù)劃分后，當(dāng)需要求解射線與空間對(duì)象的交點(diǎn)時(shí)，不需要判斷射線與每個(gè)空間對(duì)象是否相交，而只需對(duì)與射線相交的包圍盒中的空間對(duì)象進(jìn)行判斷（即先對(duì)包圍盒進(jìn)行簡(jiǎn)單的求交運(yùn)算，然后略過(guò)與射線不相交的包圍盒），從而達(dá)到減少求交次數(shù)和計(jì)算時(shí)間的目的。

Agate等人[11]提出了一種自頂向下的八叉樹(shù)搜索算法，充分利用了八叉樹(shù)的空間結(jié)構(gòu)特性，與其他算法相比，計(jì)算速度更快。

射線可用參數(shù)方程表示，即

起點(diǎn)矢量為O=[Ox,Oy,Oz]T，方向矢量為D=[l,m,n]T。在判斷射線是否與坐標(biāo)軸平行的八叉樹(shù)子節(jié)點(diǎn)相交之前，可先求出射線與子節(jié)點(diǎn)邊界相交的t值：

其中x0和x1分別為子節(jié)點(diǎn)的x坐標(biāo)值的下邊界和上邊界。

類似地可以求出0yt、1yt、0zt、1zt。然后再求出tmin和tmax：

當(dāng)沿八叉樹(shù)從根節(jié)點(diǎn)向下搜索時(shí)，可以利用空間劃分的特性加速計(jì)算。八叉樹(shù)子節(jié)點(diǎn)劃分采用中間平面，因此其子節(jié)點(diǎn)的任一方向的t值恰好是對(duì)應(yīng)父節(jié)點(diǎn)t值的一半。然后，確定射線進(jìn)入節(jié)點(diǎn)的平面，再根據(jù)射線的參數(shù)確定后續(xù)的相交節(jié)點(diǎn)。如果該節(jié)點(diǎn)包含空間對(duì)象，則將其加入待求交點(diǎn)空間對(duì)象列表。遞歸調(diào)用搜索過(guò)程，直到射線穿出八叉樹(shù)為止。最后求出射線與前述空間對(duì)象列表每個(gè)元素相交的t值，則最小t值對(duì)應(yīng)的空間對(duì)象即為所求結(jié)果。

1.3 并行處理技術(shù)

隨著計(jì)算機(jī)技術(shù)的發(fā)展，為了解決單核心CPU處理能力難以繼續(xù)提高的瓶頸，多核心CPU技術(shù)得到了飛速發(fā)展和廣泛應(yīng)用。因此模擬計(jì)算中應(yīng)充分利用多核CPU的計(jì)算加速能力。在確定計(jì)算硬件條件后，接下來(lái)的主要工作是計(jì)算程序的并行化編程。

多核CPU屬于共享存儲(chǔ)的對(duì)稱多處理器[12]。OpenMP是為共享存儲(chǔ)環(huán)境編寫(xiě)并行程序而設(shè)計(jì)的一個(gè)應(yīng)用編程接口，支持Fortran、C/C++語(yǔ)言。OpenMP標(biāo)準(zhǔn)中包括一套編譯指導(dǎo)語(yǔ)句和一個(gè)支持函數(shù)庫(kù)。在并行執(zhí)行程序時(shí)，它采用 Fork/Join方式。其基本流程是：程序開(kāi)始時(shí)只有一個(gè)主線程，即程序中的串行部分都由主線程執(zhí)行，而并行部分則通過(guò)派生（Fork）其他線程來(lái)執(zhí)行。待并行執(zhí)行程序結(jié)束后（Join）才能執(zhí)行后面非并行部分代碼，如圖2所示。

分析TPMC算法計(jì)算過(guò)程，可以在網(wǎng)格循環(huán)部分和試驗(yàn)粒子循環(huán)部分進(jìn)行并行。若在試驗(yàn)粒子循環(huán)部分采用并行，則需要進(jìn)行多次的線程派生注銷(xiāo)操作，效率低。因此，需要在網(wǎng)格循環(huán)級(jí)別進(jìn)行并行，即在C/C++語(yǔ)言中具體實(shí)現(xiàn)方法是需要在循環(huán)前加入“#pragma omp parallel for”編譯指導(dǎo)語(yǔ)句。

圖2 共享存儲(chǔ)并行程序的并行執(zhí)行流程Fig.2 Flowchart of parallel program on shared memory machine

2 計(jì)算結(jié)果

2.1 圓盤(pán)返回流

對(duì)文獻(xiàn)[10]中介紹的圓盤(pán)返回流進(jìn)行計(jì)算，顯示了八叉樹(shù)劃分算法的加速效果。算例的模型如圖3所示。計(jì)算結(jié)果為返回流比率，即與背景分子碰撞后返回放氣表面的分子數(shù)與總放氣分子數(shù)的比。算例的輸入?yún)?shù)見(jiàn)表1。

圖3 模型示意圖Fig.3 Model sketch

表1 算例輸入?yún)?shù)Table 1 Input parameters for test case

對(duì)粗網(wǎng)格（90個(gè)三角形網(wǎng)格）和細(xì)網(wǎng)格（1128個(gè)三角形網(wǎng)格）2種情況進(jìn)行了比較計(jì)算。網(wǎng)格劃分如圖4所示。

圖4 兩種網(wǎng)格劃分示意圖Fig.4 Comparison between coarse meshes and fine meshes

計(jì)算機(jī)硬件配置為Intel Core2 Duo E7500雙核CPU、內(nèi)存4 GB、操作系統(tǒng)Windows XP。文獻(xiàn)[10]中的計(jì)算結(jié)果為 1.367 7×10-6，本文的計(jì)算結(jié)果為1.326 2×10-6，二者偏差為3%，符合度較好。加速效果如表2所示。

表2 計(jì)算時(shí)間對(duì)比Table 2 Comparison of calculation time

可以看出，在網(wǎng)格數(shù)較多時(shí)，八叉樹(shù)的加速效果極為顯著。并且，加速系數(shù)可用ncell/(klog8ncell)近似表示，其中k為八叉樹(shù)子節(jié)點(diǎn)中空間單元的最大個(gè)數(shù)，本文中設(shè)置為 8，ncell表示總的空間單元個(gè)數(shù)。

并行加速比定義為S=T串/T并，即串行執(zhí)行時(shí)間與并行執(zhí)行時(shí)間的比值，當(dāng)并行加速比等于計(jì)算內(nèi)核數(shù)時(shí)，稱為線性加速。本文中加速比接近計(jì)算內(nèi)核數(shù)2，說(shuō)明算法具有較好的加速特性。當(dāng)計(jì)算任務(wù)復(fù)雜，需要較短時(shí)間獲得結(jié)果時(shí)，可通過(guò)簡(jiǎn)單增加線程數(shù)目實(shí)現(xiàn)。

2.2 空間站污染仿真分析

對(duì)特定空間站構(gòu)型下的污染分布狀況進(jìn)行了初步分析?？臻g站構(gòu)型取自文獻(xiàn)[13]，共劃分24 166個(gè)三角形網(wǎng)格，網(wǎng)格劃分情況如圖5所示。可以看出，在表面曲率變化大的部位進(jìn)行了自動(dòng)加密，使得網(wǎng)格能夠較精確地貼合表面。

圖5 空間站網(wǎng)格劃分Fig.5 Meshing of a space station

已有試驗(yàn)證明太陽(yáng)電池板放氣較為嚴(yán)重[14]，因此污染源選擇其中一塊面積較大的太陽(yáng)電池板。軌道大氣來(lái)流速度方向?yàn)?z，其他參數(shù)與2.1節(jié)中算例一致。該計(jì)算采用雙線程,耗時(shí) 387.822 s。計(jì)算結(jié)果如圖6所示，顯示的是從污染源釋放的分子直接運(yùn)動(dòng)或經(jīng)過(guò)與大氣分子碰撞后運(yùn)動(dòng)至其他表面的情況，即單位面積表面網(wǎng)格上的入射分子數(shù)。結(jié)果表明，污染分布與污染源及敏感表面之間的視角因子成正比。因此，與污染源距離越近，污染越嚴(yán)重。圖6(a)中正對(duì)污染源的位置所受污染最為嚴(yán)重。圖6(b)中，背對(duì)污染源的表面，未受到污染影響。從圖6中還可以發(fā)現(xiàn)由于表面遮擋造成的污染分布不均勻。污染源表面以及與其處于同平面的另一個(gè)太陽(yáng)電池板也有分子沉積，這是由于污染分子與軌道大氣分子碰撞后返回而形成的。

圖6 污染分布計(jì)算結(jié)果Fig.6 The calculated contamination distribution

3 結(jié)束語(yǔ)

航天器的放氣污染需要定量評(píng)估其影響。采用Monte-Carlo方法進(jìn)行模擬，用射線表征放氣分子的輸運(yùn)過(guò)程，可將污染輸運(yùn)過(guò)程轉(zhuǎn)化為射線追蹤問(wèn)題。由于航天器構(gòu)型復(fù)雜，射線追蹤的計(jì)算很耗時(shí)。本文將復(fù)雜構(gòu)型航天器表面用三角形非結(jié)構(gòu)網(wǎng)格表示，用八叉樹(shù)劃分算法開(kāi)展分子運(yùn)動(dòng)軌跡與表面相交的加速計(jì)算。同時(shí)，利用基于共享內(nèi)存多核計(jì)算平臺(tái)的OpenMP軟件實(shí)現(xiàn)了計(jì)算的并行化。計(jì)算實(shí)例證明該算法可大大地提高計(jì)算速度，且與串行計(jì)算結(jié)果一致。此外，上述算法還可應(yīng)用于輻射傳熱角系數(shù)計(jì)算、空間粒子輻射影響分析等。

（

）

[1]陳烈民.航天器結(jié)構(gòu)與機(jī)構(gòu)[M].北京: 中國(guó)科學(xué)技術(shù)出版社, 2005: 63-70

[2]侯增祺, 胡金剛.航天器熱控制技術(shù): 原理及其應(yīng)用[M].北京: 中國(guó)科學(xué)技術(shù)出版社, 2007: 129-140

[3]陳樹(shù)海.衛(wèi)星結(jié)構(gòu)中的非金屬材料[J].上海航天, 2004, 21(3): 39-44 CHEN S H.Nonmetal materials in satellite structure[J].Aerospace Shanghai, 2004, 21(3): 39-44

[4]MILLER S K R, BANKS B.Degradation of spacecraft materials in the space environment[J].MRS Bulletin, 2010, 35(1): 20-24

[5]SCHLAPPI B, ALTWEGG K, BALSIGER H, et al.Characterization of the gaseous spacecraft environment of Rosetta by ROSINA[C]//3rdAIAA Atmospheric Space Environments Conference.Honolulu, HI, United States: American Institute for Aeronautics and Astronautics, 2011: 1-8

[6]HAEMMERLE V R, GERHARD J H.Cassini camera contamination anomaly: experiences and lessons learned: AIAA 2006-5834[R]

[7]ROUSSEL J F, SOARES C E, SCHMIDL W D.Development and validation of a new return flux model for International Space Station[C]//Optical System Contamination: Effects, Measurements, and Control VII, Proceedings of SPIE, 2002, 4774: 172-183

[8]焦子龍, 龐賀偉.航天器分子污染返回流計(jì)算方法[J].宇航學(xué)報(bào), 2009, 30(2): 415-421 JIAO Z L, PANG H W.Simulation methods for return flux of molecular contamination[J].Journal of Astronautics, 2009, 30(2): 415-421

[9]焦子龍.航天器分子污染計(jì)算方法綜述[J].航天器環(huán)境工程, 2010, 27(2): 181-189 JIAO Z L.Review on simulation method for spacecraft molecular contamination[J].Spacecraft Environment Engineering, 2010, 27(2): 181-189

[10]FAN C, GEE C, FONG M C.Monte Carlo simulation of molecular flux on simple spacecraft surfaces due to selfand ambient-scatter of outgassing molecules: AIAA 1993-2867[R]

[11]AGATE M, GRIMSDALE R L, LISTER P F.The HERO algorithm for ray tracing octa-trees[M]//GRIMSDALER L, STRASSER W.Advances in computer graphics hardware IV.New York: Springer-Verlag, 1991: 23-27

[12]周偉明.多核計(jì)算與程序設(shè)計(jì)[M].武漢: 華中科技大學(xué)出版社, 2009: 5-10

[13]翟睿瓊, 姜海富, 田東波, 等.空間站原子氧環(huán)境仿真研究[J].裝備環(huán)境工程, 2014, 11(3): 35-39 ZHAI R Q, JIANG H F, TIAN D B.et al.Simulation of atomic oxygen interaction with space station[J].Equipment Environmental Engineering, 2014, 11(3): 35-39

[14]焦子龍, 龐賀偉, 易忠, 等.航天器真空熱試驗(yàn)污染物成分分析[J].航天器環(huán)境工程, 2010, 27(6): 711-714 JIAO Z L, PANG H W, YI Z.et al.Analysis of contaminant components of spacecraft thermal vacuum test[J].Spacecraft Environment Engineering, 2010, 27(6): 711-714

（編輯：肖福根）

Parallel and octa-tree-based acceleration of outgassing contamination simulation for spacecraft of complex geometries

JIAO Zilong, JIANG Lixiang, HUANG Jianguo, SUN Jipeng, ZHU Yunfei
(Science and Technology on Reliability and Environmental Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering, Beijing 100094, China)

The organic materials used in spacecraft may outgas under the space environment, including, mainly, the vacuum, the high & low temperature, and the influence of the ionization irradiation.The outgassing product has harmful effects on the spacecraft itself.So a quantitative assessment of the contamination effects is necessary.The molecular transport process is modeled by the Monte-Carlo method, with the rays representing the molecule trajectories, to simulate the motion process of a large number of molecules.Thus the transport process is transformed into a ray tracing problem.Due to the complexity of the spacecraft configurations, the ray tracing computation might be very time consuming, therefore, it is extremely necessary to have an accelerating algorithm.In this paper, the complex geometry of the spacecraft surface is modeled with triangular unstructured grids, the octa-tree space subdivision technology is used to accelerate the molecular motion and surface intersection calculations.At the same time, based on the shared memory multi-core computing platform, the OpenMP API is used to implement parallel computing for the simulation.The test case shows a significant acceleration in the calculation, and the developed method is applied to the contamination analysis of a space station.

spacecraft; outgassing contamination; simulation computation; complex geometries; octa-tree; parallel computing

X502

：A

：1673-1379(2016)06-0601-05

10.3969/j.issn.1673-1379.2016.06.005

焦子龍（1979—），男，博士學(xué)位，高級(jí)工程師，從事航天器特殊空間環(huán)境效應(yīng)分析、試驗(yàn)及防護(hù)技術(shù)研究。E-mail: novanova2009@163.com。

2016-04-01；

：2016-11-10