類天宮飛行器軌道衰降過程空氣動力特性一體化建模并行優化設計

張子彬，李志輝*，白智勇，彭傲平

（1.中國空氣動力研究與發展中心超高速空氣動力研究所，綿陽621000；2.北京航空航天大學北京前沿創新中心國家計算流力學實驗室，北京100191）

1 引言

航天器軌道計算與飛行航跡數值預報的準確可靠直接依賴于軌道動力學方程空氣動力項求解的快速高效［1］，研究建立適于大型航天器在軌與服役期滿軌道衰降過程全飛行流域高超聲速空氣動力特性一體化模擬平臺［2］與融合修正沿軌道空氣動力一體化快速計算方法是關鍵基礎［3］。當航天器運動在200 km以上時，飛行繞流已屬于高稀薄自由分子流，此時氣體分子平均自由程將超過200 m，Knudsen數將大于10。對于這種情況，傳統觀點認為航天器飛行阻力、升力等空氣動力將不再隨飛行高度發生變化，因此應該采用固定的氣動力系數參與飛行動力學軌道控制的［3］。但是，中國天宮飛行器在采用固定阻力系數2.2開展低軌道飛行試驗時，發現有時飛行偏差越來越大，須使用RCS姿控系統強制控回原標稱軌道，造成燃料消耗，甚至難以準確有效控制［4］。因此，提出這樣的問題：天宮飛行器低軌道飛行所受大氣阻尼系數是否會變化？為了解決這一問題，李志輝等［4］使用修正Boettcher/Legge非對稱橋函數，發展了適于高超聲速稀薄過渡流區氣動力計算的當地化橋公式法，建立了專門針對復雜巨大的天宮飛行器跨越大氣層不同高度、不同馬赫數、不同迎角與側滑角空氣動力特性快速計算方法與程序軟件，以減小數值計算工作量，并融合彈（軌）道飛行動力學計算推進，得到了具有實用價值的結果。實際應用發現，在空氣動力融合軌道飛行航跡數值預報過程中，該類算法程序的運算速度與能力效率仍有待進一步提升［3］。為此，本文研究基于CUDA架構的GPU異構高性能并行計算技術，開展大型航天器軌道衰降過程空氣動力特性一體化建模并行優化設計與應用嘗試。

隨著科學技術的迅猛發展，高性能計算已成為科學技術發展和重大工程設計中具有戰略意義的研究手段［5］。早期計算機性能的提升主要依賴CPU（中央處理器）主頻的提升，但由于CPU的功耗與頻率的三次方近似成正比，無限制地提升頻率已不可能。到2003年，CPU頻率的提升接近停止［6］。為了應對這一情況，人們往往采用基于CPU的并行計算方法來解決這一問題（例如消息傳遞接口，簡稱MPI），但由于單個的CPU最多只能集成十幾、幾十核，導致基于CPU的并行計算能力受到很大限制。

2007 年，隨著基于英偉達CUDA架構圖形處理單元（GPU）硬件的出現，集成核心的能力比CPU要強大得多（例如，英偉達Tesla架構的K20X芯片包含2688個核心），這使得GPU的運算能力比CPU更勝一籌［7］。除此之外，CUDA架構不但包含英偉達GPU硬件，而且包含一個軟件編程環境，在這樣的編程模型下，CPU負責整個任務的流程控制和任務的串行計算，而GPU采用并行算法完成任務中計算密集部分的計算，人性化的編程環境有力助推了GPU的廣泛應用。

時至今日，GPU和并行技術的使用已成為當今數值計算的潮流和趨勢，并被廣泛應用于包括計算流體力學在內的許多領域。同樣地，它們也為求解稀薄氣體流動問題提供了一條行之有效的途徑。特別是使用NVIDIA的Tesla GPU在計算流體動力學中的應用加速明顯，針對不可壓縮的納維—斯托克斯（Navier-Stokes）方程數值求解模型的并行實現，在1024×32×1024的規模下使用4個Tesla C870 s與使用AMD Opteron 2.4 GHz的性能相比約提高100倍；針對格子Boltzmann方法在格子大小128×128規模下，使用Tesla 8-series GPU相比于使用Intel Xeon的性能約提高了130倍［8］。可見，GPU和CPU融合的并行架構是未來加速通用計算必然趨勢。

本文基于英偉達CUDA架構的GPU，以微型計算機為平臺，根據近地軌道類天宮飛行器空氣動力特性當地化橋函數關聯快速算法的運算流程及對應程序軟件的編寫特點，提出相應并行優化方案，并付諸實施，以提升程序計算效率，減少程序運行時間。并以類天宮一號目標飛行器無控飛行軌道衰降過程空氣動力特性一體化建模作為算例，給出其在低軌道飛行過程相應高度飛行繞流狀態的氣動特性，檢驗檢測本文并行計算模型方案的并行加速有效性與高效率，以此為基礎，逐步開拓GPU并行計算在融合彈（軌）道數值預報關于空氣動力計算的快速高效與準確性。

2 空氣動力特性當地化計算方法

融合航天器在軌飛行與離軌再入飛行航跡彈（軌）道數值預報的跨流域高超聲速空氣動力計算所用的當地化計算方法，是基于自由分子流與連續流當地橋函數理論關聯的橋公式法［9-10］。自由分子流和牛頓連續流作為航天器再入飛行過程2個對應極限流態，其中稀薄過渡流區域的氣動力系數可以用當地氣動力的解析來表述，當地面元氣動力系數采用基于自由分子流與連續流的當地化橋函數光滑搭接［11］。這是一種加權函數，依賴于獨立的關聯或相似參數（如克努森數Kn），通常由實驗結果與數值計算分析確定。橋公式法根據自由分子流和牛頓連續流不同的物理定律，采用兩個不同的數學函數FCont和FFM來描述，而在中間區域給出一個新的函數Fb，即橋函數。由此，對于給定的入射角θ，當地面元上經歸一化壓力和摩擦力系數表達為式（1）［1］：

式中，Fb，p和Fb，τ分別是壓力和摩擦力橋函數，依賴于獨立參數Kn、TW/T0和θ。

對過渡流區域，當地橋公式［12］可以提供依賴于獨立參數的歸一化壓力系數和摩擦力系數。當地橋函數的表達式有許多種，本文研究發展修正的Boettcher/Legge非對稱橋函數表達式，其中，非對稱的壓力橋函數可表示為式（2）：

式中，Knn，p、ΔKnp1、ΔKnp2根據飛行器外形及繞流特點為可調關聯參數。

非對稱摩擦力橋函數可表示為式（3）：

發展可分段描述的非對稱壓力與摩擦力系數關聯橋函數快速計算模型，需要根據飛行器外形及繞流特點調試確定Knn，p、ΔKnp1、ΔKnp2及Knn，τ、ΔKnτ1、ΔKnτ26個關聯參數，由于天宮飛行器的大尺度復雜結構與在軌飛行跨流域多尺度非平衡繞流特點，需要依靠高性能大規模并行計算機開展求解Boltzmann模型方程氣體動理論統一算法等跨流域空氣動力學數值計算典型飛行狀態模擬結果［2］，最優擬合計算修正確定，使當地化工程計算結果最大程度地滿足實際計算精度［9］。

關于飛行器物形處理，選取一種較為通用和近似的面元非結構網格物形表征處理方法［4］，將飛行器物面劃分為若干個面元，利用四邊形或三角形來逼近復雜外形，圖1繪出類似表征天宮飛行器兩艙體的大尺度圓柱體計算模型表面三角非結構網格布局，可計算出每個面元的中心點坐標、中心點處的法向及面元面積［13］。由于類天宮飛行器屬于復雜非規則巨大板艙組合體，為了解決類天宮飛行器及附件復雜物形表征的困難，引入分區網格與非結構網格生成技術，進行復雜物形表征處理，以便計算面元所受氣動力，將所有面元上的氣動力加起來就可以得到整個飛行器的氣動力。

圖1 類天宮飛行器兩艙體圓柱體外形表面非結構三角形面元網格Fig.1 Unstructured triangular surface elementmesh on cylinder surface for two-capsule body of Tiangong type spacecraft

3 并行優化方法

傳統的空氣動力特性當地化計算程序沒有采用并行方案，因此只能使用其中1個核心和1個線程，且每次只能執行1次循環，造成了計算資源的閑置與浪費。為了滿足氣動融合軌道飛行航跡數值預報對空氣動力特性實時計算快速高效要求，需要發展多核處理單元的并行優化方法［14］。

為此，對類天宮飛行器空氣動力特性當地化算法的運算流程及對應程序代碼進行整體分析發現：

1）程序熱點代碼較為集中。原有程序主要執行的代碼集中在算法的求解部分，具體表現為多個連續的3層嵌套循環。其中最內層的循環次數最多，達到十萬級別，適合數百個以上核心同時計算，有利于發揮GPU設備的性能優勢。

2）數據依賴性較弱。數據依賴是指下一條對數據操作的指令必須等待上一條操作該數據的指令完成。原程序各個循環間保持相互獨立，不存在數據依賴關系，有利于發揮GPU設備的性能優勢。

3）循環內外數據傳輸需求較少。長期以來，GPU與CPU之間數據傳輸較慢，一直是GPU在應用過程中的詬病。具體來說，由于GPU與CPU之間的數據傳輸速度遠低于GPU、CPU芯片與各自寄存器的傳輸速度，甚至不及CPU與內存、GPU與顯存之間的傳輸速度，因此造成引入GPU設備被迫增加的數據傳輸時間，會小于引入GPU設備所節省時間的情況，即損耗大于收益。而本文所涉及的程序由于循環內外數據傳輸需求較少，可以避免這一情況，因此有利于發揮GPU設備的性能優勢。

由上述特點可見，原有程序擁有改寫為并行程序的潛力，并可以較好地發揮GPU設備的性能優勢。因此，針對原有程序的編寫特點，我們引入CUDA架構的GPU設備，根據調試過程中出現的問題，分別采用系統級別、算法級別以及語句級別3個層次的優化手段對其進行優化。

3.1 系統級別優化

系統級別手段優化主要考慮程序性能的控制因素。對本文要解決的問題，首先分析原程序是否存在以下3個限制因素：①處理器利用率過高還是過低；②存儲器帶寬利用是否不足；③阻塞處理器運算的其他因素。

CPU、GPU及存儲器連接關系如圖2所示，測試發現，GPU與CPU以及GPU與顯存的數據傳輸速度比較緩慢，因此對這兩部分著重進行了優化設計如下：

1）對齊訪問。由于在GPU顯存數據傳輸時，可以通過32、64或128 Byte的作業來訪問，這些作業都對齊到它們的尺寸。因此未對齊的訪問會增加訪問次數，浪費傳輸資源。就本程序而言，由于算法使用GPU運算時有很多零散的參數需要傳入，因此將其整合為一個數組，進行一次性傳輸，方便訪問數據的對齊。圖3為本文設計使用的GPU內存對齊訪問格式。

圖2 CPU、GPU及存儲器連接關系示意圖Fig.2 Connection diagram of CPU，GPU and correspondingmemory

圖3 基于CUDA架構的GPU內存對齊訪問Fig.3 Aligned access of GPU memory based on CUDA

2）使用數據傳輸函數。由于GPU設備傳統的數據傳輸方式效率較低，因此高版本CUDA架構的GPU設備提供了豐富的數據傳輸函數，可以完成CPU到GPU或者GPU到另一臺GPU設備的數據傳輸，通過測試，可選用cudaMemcpy函數作為傳輸接口，具體格式為：cudaMemcpy（dst，src，count，kdir）。其中，dst、src分別為目標變量和源變量，count為拷貝的數據長度，kdir為可省略參數。

3.2 算法級別優化

算法級別的優化通常涉及程序部分，而這部分包含一個或多個函數，甚至只有一段語句。算法級別優化主要涉及算法實現時要考慮的問題，通常算法要考慮數據的組織、算法實現的策略。

一般串行算法求解的程序部分表現為多個連續的三層嵌套循環，且循環主要集中在最內層，達到十萬級別。以NVIDIA公司的GTX750Ti圖形處理單元為例，它擁有640個計算核心，理論上可支持131 072個線程同時計算［14］。從線程數來說，每個循環都可以交給一個指定的線程進行計算，具有較好的并行性；從核心數來講，十萬余個循環可以拆分給多個GPU芯片共同完成，具有較好的可擴展性。

因此，本文將算法求解部分的最內層循環經過展開與合并處理后，整理為一個循環，再整體移植入核函數，從而利用GPU強大的并行計算能力提升運算效率。具體代碼如下所述。

原有串行程序代碼為：

do i=1，ntr！ntr=121822

……！數組的循環運算A

end do

……！數組的共同運算B

do i=1，ntr！

……！數組的循環運算C

end do

……！數組的共同運算D

經過優化的程序代碼為：

attributes subroutine CUDA（input，output）

！GPU設備核函數開頭

do i=1，ntr！ntr=121822

……！運算A+B+C+D

end do

end subroutine

3.3 語句級別優化

語句級別的優化包括對函數、循環、指令等語句細節的優化。

CUDA架構的GPU設備顯式地將存儲器分成4種不同的類型：全局存儲器（globalmemory）、常量存儲器（constantmemory）、共享存儲器（share memory）或稱局部存儲器（localmemory）、私有存儲器（private memory）。雖然其中除全局存儲器有較大容量之外（如tesla K20全局內存為5 GB），其余3種存儲器的存儲容量只有64KB左右［14］，但是這4種存儲器的讀取速度卻有很大差異：全局存儲器是最慢的，其延遲為幾百個時鐘；常量存儲器次之，延遲為幾十個周期；共享存儲器的延遲在幾個至幾十個時鐘周期之內；而私有存儲器通常就是硬件的寄存器，它是GPU的存儲器中讀取最快的，幾乎不用耗費時間［14］。而相比較全局變量耗費的幾百個周期，1個時鐘周期可以做4次整形加法、2次浮點乘加，3個周期可以做一次乘法，十幾個周期可以做一次乘法［15］。因此如果能夠合理利用常量存儲器、共享存儲器和私有存儲器，盡量減少讀取全局存儲器所產生的時間耗費，那么所產生的收益將是相當可觀的。

以下列出了私有存儲器、共享存儲器和常量存儲器的聲明使用方式：

integer：：t，tr！ 使用私有存儲器

real，shared：：s（64） ！ 使用共享存儲器

real，parameter：：PI=3.14159！ 使用常量存儲器

4 算法驗證與結果分析

將所建立類天宮飛行器空氣動力特性當地化串行算法及對應程序按提出的并行計算方案進行CPU+GPU移植優化后，對類天宮一號目標飛行器外形進行運算，分析優化前后程序運行性能。

4.1 計算模型與條件

實施對于天宮一號目標飛行器低地球軌道飛行過程氣動力系數一體化計算，采用基于三角非結構網格表征物面，繪出計算模型表面三角非結構網格布局如圖4所示。

圖4 天宮飛行器計算模型表面非結構網格布局Fig.4 Unstructured triangular surface elementmesh on the com puting model of Tiangong spacecraft

關于極不規則大型復雜結構航天器氣動特性計算，參考面積選取可分為2部分考慮：第一部分是本體迎風面積，即飛行器本體的底面積；第二部分是帆板的迎風面積，由于太陽電池帆板需要圍繞中心軸轉動，于是帆板迎風面積采用等效面積計算；總的參考面積為上述兩部分面積之和［4］。

本算例采用多核CPU+GPU架構計算機進行運算，CPU及GPU性能參數見表1所列。對比試驗中，串行程序僅由CPU獨立完成，共調用1個核心和一個線程。并行程序的算法求解運算由GPU完成，調用640個核心和131 072個線程，其余步驟如數據載入、結果收集及試驗監測等過程仍由CPU完成。

表1 試驗設備性能列表Table 1 Performance list of test equipments

4.2 GPU并行算法與CPU串行計算結果驗證

圖5與圖6分別繪出天宮飛行器角α為22°，阻力系數與升力系數隨側滑角β變化趨勢其中藍色曲線表示CPU串行程序計算結果，而紅色方格代表并行程序計算結果。由圖看出，兩者完全吻合，本文經過對156次計算結果進行統計，除部分個別數值末位略有不同且結果僅差0.000 01，誤差緣由可能是CPU和GPU芯片取舍原則不同所致，其余數據完全一致，這種GPU與CPU計算的高精度一致性，進一步證實本文GPU并行算法實現的正確性與高效率。

圖5 GPU加速前后阻力系數隨側滑角變化計算驗證Fig.5 Computing verification of drag variation w ith sideslip angle before and after GPU acceleration

圖6 GPU加速前后升力系數隨側滑角變化計算驗證Fig.6 Com puting verification of lift variation w ith sideslip angle before and after GPU acceleration

表2列出原串行程序CPU計算與本文GPU并行算法程序及優化后的GPU并行程序在求解過程中，單步求解各模塊程序段的平均消耗時間，其中GPU運算是在GPU中發生的過程，數據歸約是在CPU中各節點進行。由表2看出：①GPU并行程序雖然增加了數據輸入、輸出的過程，但總的運行時間仍比CPU串行程序計算時間有了明顯減?。虎诮涍^并行優化后，數據在CPU與GPU之間傳輸的時間也有了明顯減少。

表3給出了CPU串行、并行及引入GPU設備并優化后的程序平均消耗時間，從表中可以看出，受制于CPU線程與計算核心數影響（CPU i7-4770 K只能使用8個線程與4個處理單元），CPU并行程序在單步求解時間方面的加速比只能達到7.0左右，而無法超過8，而GPU設備的相應加速比可以達到13.0，受此影響，CPU并行程序與優化后的GPU并行程序的整體加速比仍有一定差距，雖然近年來CPU的最大線程數和性能有所提升，但GPU設備的并行能力也在大幅度增長，因此仍然沒有實質性的區別，這樣的案例進一步說明了采用GPU并行的重要性和研究價值。

表2 算法單步求解各階段平均消耗時間Table 2 Average runtime per step in solving process/ms

表3 CPU串行、并行及GPU優化后程序平均消耗時間Table 3 The average runtime time of CPU serial program，parallel program and GPU optim ized program

4.3 性能分析

根據常用并行算法性能衡量手段，采用加速比對計算試驗進行分析。加速比的定義為順序執行時間與并行執行時間比值。由表3可知，CPU串行程序消耗時間0.747 h，GPU并行程序消耗時間0.140 h，因此加速比為5.3，這說明，將GPU并行加速引入類天宮飛行器無控飛行軌降再入解體過程氣動力系數一體化計算相當快速高效。

比較分析還可看出：

1）在單次求解過程中，CPU循環運算平均消耗75.5 ms，而優化后的GPU運算平均消耗5.8 ms，GPU并行計算加速比約為13?？梢姡?章引入GPU設備、發展的GPU并行優化策略是行之有效加速手段；

2）通過優化，單次求解過程中的數據傳輸耗時由原來的9.2 ms降為2.1 ms，減少了77%；數據歸約耗時由原來的2.0ms降為0.7 ms，減少了65%；GPU運算過程也由原來的4 ms降為3 ms；合計時間由優化前的GPU運行時間15.2 ms降為5.8 ms，減少了61%?？梢姡疚慕榻B的GPU優化手段成效顯著；

3）在試驗中共使用線程512×256=131 072個，而GPU計算核心僅調用640個，這是因為每個核心上被分配了多個線程，這增加了線程等待時間，因此，當調用核心數增加時，程序仍有進一步擴展的空間；

4）在單次求解過程中，GPU并行程序數據輸入輸出的平均耗時為2.6+6.6=9.2 ms，是單次求解過程平均耗時（15.2 ms）的近60%，即使在GPU并行優化后，輸入輸出的平均耗時為0.9+1.2=2.1ms，是單次串行求解過程平均耗時（5.8 ms）的近36.2%?？梢姡捎贕PU與CPU之間的傳輸速度過慢，輸入輸出占用了GPU調用的相當多時間，因此，如何發展高效快速的輸入輸出仍是限制GPU推廣應用需要進一步研究解決的瓶頸之一。

4.4 類天宮飛行器軌道衰降過程空氣動力特性計算分析

使用類天宮飛行器軌道衰降過程空氣動力特性一體化建模GPU并行算法程序，對天宮飛行器無控飛行軌道衰降過程340～120 km氣動特性詳細計算，圖7與圖8分別繪出太陽電池帆板面與天宮飛行器本體主軸夾角10°、迎角α=10°進行340～120 km不同高度340 km、260 km、180 km、140 km、120 km的氣動阻力、升力隨不同側滑角變化輪廓。可看出對于所計算的不同飛行繞流狀態，均在側滑角β=0°阻力系數最小、升力系數最大，說明天宮飛行器存在嚴重的板艙非規則構型，一旦無控飛行姿態出現側滑角飛行繞流現象，迎風面積、阻力系數都會迅速增大，升力系數快速減小，致失穩加速飛行繞流狀態，印證了天宮一號目標飛行器無控飛行軌降過程數值預報發現的，其自2017年9月初出現失穩后三軸旋轉加速直至最終再入大氣層解體的無控飛行繞流現象。

圖7 天宮飛行器軌道衰降過程氣動阻力系數GPU優化并行計算結果與CPU串行計算驗證Fig.7 Verification of GPU parallel program and CPU serial computation for aerodynam ic drag variation of Tiangong spacecraft during orbital declining at 340～120 km

圖9（a）、（b）分別繪出太陽電池帆板面與天宮飛行器本體主軸夾角10°、α=22°、β=0°時天宮一號目標飛行器軌降進入臨近再入段240～120 km不同高度240 km、220 km、200 km、180 km、160 km、140 km、120 km的氣動阻力、升力系數隨飛行高度變化趨勢。可看出對于所計算的不同飛行繞流狀態，隨著高度降低，阻力系數、升力系數呈現快速非線性下降趨勢，印證了天宮一號飛行航跡在此高度軌降過程中，將于20 d內到達120 km再入大氣層解體墜毀的無控飛行變化現象。

5 結論

1）通過分析近地軌道類天宮飛行器空氣動力特性當地化橋函數關聯算法及其對應程序具有多連續三層嵌套循環熱點代碼集中、循環間彼此獨立不存在數據依賴關系及循環內外數據傳輸需求少等特點，引入了CUDA架構的GPU并行加速設備，提出了基于多核處理單元的并行優化設計方案。

圖8 天宮飛行器軌道衰降過程氣動升力系數GPU優化并行計算結果與CPU串行計算驗證Fig.8 Verification of GPU parallel program and CPU serial computation for aerodynam ic lift variation of Tiangong spacecraft during orbital declining at 340～120 km

圖9 天宮飛行器軌道衰降過程氣動系數隨飛行高度GPU并行優化計算與CPU串行計算驗證Fig.9 Verification of GPU parallel program and CPU serial com putation for aerodynam ic drag and lift variation of Tiangong spacecraft during orbital declining at 340～120 km.

2）針對氣動融合軌道數值預報快速高效特點，發展了系統、算法、語句級別的3個層次并行優化方法。建立了基于GPU內存對齊訪問格式與數據傳輸函數方案、循環展開與合并、優化存儲器使用等CPU+GPU并行移植優化技術，快速提升了GPU數據傳輸速度和運算效率。

3）經本文GPU并行算法與CPU串行計算試驗驗證表明，類天宮飛行器空氣動力特性當地化工程算法程序經GPU并行移植優化后，運行程序消耗時間由CPU串行計算0.747 h，減小至GPU并行計算0.140 h，并行加速比為5.3，證實首次將GPU并行加速引入類天宮飛行器無控飛行軌降再入解體過程氣動力系數一體化計算相當快速高效。

4）GPU作為圖形處理工具，在架構上為了集成更多的計算核心而犧牲了存儲性能和單核計算能力，GPU雖然在并行加速能力上遠勝于CPU，但是在數據傳輸及單核處理能力方面較弱。計算程序只有在被改寫成數據依賴性好，傳輸需求少、計算量大的運算模式，才能發揮GPU設備優勢。本文提出的并行優化設計方案，可以為其他算法程序在大規模并行計算實現方面提供經驗和借鑒參考。

5）本文工作屬初步階段性，將GPU并行移植優化算法應用于其他跨流域空氣動力學數值模擬方法，特別是GPU與CPU之間傳輸速度慢，輸入輸出占用GPU調用較多時間，如何發展高效快速的輸入輸出是限制GPU推廣應用瓶頸，均有待進一步研究解決。