航天器著陸碰撞動力學CUDA并行計算技術

高樹義 谷良賢 張方 陳國平 李楊

（1 西北工業大學航天飛行動力學技術重點實驗室，西安 710072）（2 南京航空航天大學機械結構力學及控制國家重點實驗室，南京 210016）

1 引言

航天器著陸緩沖過程的瞬態分析問題，主要研究結構在碰撞過程的時間歷程響應，計算結構在振動時的動位移和動應力的大小，以及變化規律，是結構動力學的一個重要研究方向。物體間的接觸碰撞是復雜的力學問題，是引起物體失效和破壞的重要原因。接觸碰撞過程在力學上通常涉及三種非線性問題：材料非線性、幾何非線性和接觸非線性。接觸非線性來源于兩個方面：接觸界面隨時間發生變化，要在求解過程中確定；接觸條件的非線性。接觸條件又包括兩個方面：接觸體之間在接觸面的變形保持協調，即不可貫入性；切向接觸的摩擦條件。根據是否考慮接觸體的變形，可以將物體間的接觸分為剛性接觸和彈性接觸。在彈性體的接觸碰撞分析中，首先要將接觸體進行有限元離散；然后將約束條件引入離散系統的運動微分方程，再通過直接積分法或者振型疊加法對系統的響應進行求解。

由于接觸碰撞問題涉及非線性問題，并且需要求解時間歷程的瞬態響應，其計算過程復雜，計算量大，因此并行計算在這類問題的處理上具有較大的應用價值。文獻［1］開發了適應向量機的小球算法；文獻［2］設計了三維殼體結構的沖擊接觸并行算法；文獻［3］設計了接觸碰撞問題的并行算法；文獻［4］在多種并行機型上設計了接觸搜索的并行算法，該算法在大規模工程問題上得到了應用。統一計算架構（CUDA）是顯卡廠商NVIDIA 在2006年推出的基于圖形處理器（GPU）的并行運算平臺。其特點是能充分發揮圖形顯卡的潛力，通過顯卡對密集數據實施高效率的并行計算。CUDA 采用C 語言編程，這在很大程度上簡化了并行計算的編程，已經被推廣應用到某些領域，如分子動力學模擬、計算天體物理的N-Body算法模擬、數字天氣預報的加速處理、電磁模擬和地球物理數據處理等。

目前，CUDA 并行算法在航天器著陸碰撞仿真領域還處于研究起步階段。將并行計算應用到有限元計算上，要從現有的串行算法入手，分析各個環節的并行性，包括有限元模型的建立、計算分析的過程和計算結果的后處理，然后在并行計算平臺上設計優良的并行計算程序［5-9］。本文對接觸碰撞問題的顯式積分計算方法進行了討論，研究了航天器碰撞問題的并行算法，并通過算例驗證了該并行算法的高效性和準確度。

2 接觸碰撞問題的顯式積分過程

2.1 接觸碰撞問題的基本描述

如圖1所示，物體A 和B是兩個發生接觸碰撞的物體，虛線0VA和0VB是它們在接觸碰撞之前的位形，tVA和tVB是它們在t時刻發生接觸碰撞時的位形，tΓc是這兩個物體發生接觸碰撞時的接觸面，該接觸面在兩個物體中分別是tΓAc和tΓBc。通常，稱物體A 為接觸體（Contactor），物體B為目標體或靶體（Target），則tΓAc和tΓBc表示從接觸面和主接觸面。

圖1 物體A 和B在接觸碰撞時的位形Fig.1 Arrangement figure of object A and object B in contact

接觸物體在接觸面上要滿足不可貫穿性和法向接觸力為壓力的法向接觸條件，一般將這兩個法向接觸條件稱作接觸碰撞的運動學條件和動力學條件。

2.2 接觸問題計算方法

中心差分法是直接積分法的顯式求解方法，系統在每一時刻的運動狀態可以完全由上一時刻確定，因此可以避免每一增量步的迭代，這在非線性的瞬態分析中具有很重要的意義。由于中心差分法是條件穩定的算法，因此在利用其求解具體問題時，所選取的時間步長Δt必須不大于某個臨界步長Δtcr，否則，算法是不穩定的。

中心差分法的穩定性條件為

式中：ωn為離散系統的最高階固有圓頻率；Tn為離散系統的最小固有周期；n為階次。

在利用中心差分法求解接觸問題時，其迭代公式為

式中：tu和t＋Δtu分別為t和t＋Δt時刻的位移向量；和分別為t時刻的速度向量和加速度向量；和分別為t＋Δt時刻的速度向量和加速度向量。

從式（2）可以看出，下一增量步的增量解可以完全由上一步確定。

2.3 基于罰函數法的有限元方程計算機實施步驟

在全Lagrange格式下，接觸碰撞問題的有限元方程為

式中：M為質量矩陣；t＋ΔtQL和t＋ΔtQc分別為系統在t＋Δt時刻的等效節點力向量和接觸力向量；為系統在t＋Δt時刻的節點內力向量。

計算機實施流程簡圖見圖2，具體步驟如下。

（1）令t＝0，并選擇積分步長Δt、等效節點力QL、接觸力Qc和節點內力F的初值；

（2）形成集中質量矩陣M，并計算M-1；

（5）計算t＋ΔtQL；

（6）由節點內力公式fe＝kTδe計算節點內力的初始值，然后再計算節點應力和節點內力，其中kT為剛度矩陣，δe為節點位移；

（7）搜尋接觸點并判斷接觸點的接觸狀態；

（8）計算接觸力t＋ΔtQc；

（11）令t＝t＋Δt；

（12）若t＜T，返回第4步，否則計算結束。

圖2 實施流程Fig.2 Implementation flow chart

3 接觸碰撞問題的并行計算

多核CPU 和多核GPU 的出現，使“分而治之”的并行處理思想獲得硬件上的支持，隨著處理器核心的不斷擴展與成熟，并行系統已經開始成為主流。CUDA 是利用GPU 進行科學計算的全新并行計算體系架構，成為并行計算的一種實現模式，其主機、GPU 數據之間的傳遞關系如圖3所示［10］。

圖3 并行計算數據傳遞Fig.3 Data transfer of parallel computing

在接觸碰撞問題的處理上，除了一般的有限元分析步驟，還增加了接觸的搜索和接觸力的計算兩個過程。此外，結構的質量矩陣采用集中質量矩陣形式，這就使得運動微分方程得以解耦，簡化方程組的求解。其有限元計算的主要步驟如下。

（1）輸入相關參數并計算集中質量矩陣。

（2）時間歷程的循環：①節點內力的計算；②接觸搜索并計算接觸力；③根據運動微分方程求解每一迭代步的運動狀態；④根據需要輸出中間結果。

（3）結果的后處理。

瞬態分析采用的積分步長一般較小，完成瞬態分析過程需要較長的時間，因此第2步是并行計算的主要步驟［11］。在這個步驟中，并行計算要完成節點力、接觸力的搜索和計算，由于主機和設備之間數據的傳輸耗時取決于PCI-E 接口的帶寬，因此這一步要盡量減少主機和設備的數據傳輸。

在有限元計算中，接觸碰撞問題并行計算步驟如下。

1）質量矩陣和剛度矩陣的并行生成

采用CUDA 提供的原子函數atomicAdd（），其函數原型為：

unsigned long long int atomicAdd

（unsigned long long int＊address，unsigned long long int val）；

從函數原型中可以看到，atomicAdd（）只能進行整形數據的操作，由于計算能力支持64位數據的操作，因此可以先將質量矩陣的元素乘以一個大數（相乘之后的數值不能超過263-1），然后在完成質量矩陣的疊加之后除以這個大數。質量矩陣采用集中質量矩陣形式，按照單元數目為并行計算配置線程，一個線程控制一個單元，在每個單元內進行相關計算。在完成各個單元計算之后，要按照單元和節點對質量矩陣進行疊加。由于各個線程并行執行，在質量矩陣的疊加過程中會發生線程沖突，使矩陣在相應的疊加位置上只計入最后一個線程所計算的數值，這就要調用CUDA 的原子函數atomicAdd（）對疊加計算進行處理。其函數原型為：

int atomicAdd（int＊address，int val）；

unsigned int atomicAdd（unsigned int＊address，unsigned int val）；

unsigned long long int atomicAdd（unsigned long long int＊address，unsigned long long int val）；

從上面的函數原型可以看到，atomicAdd（）不支持雙精度浮點數。為了有效利用這一函數，可以將質量矩陣的元素同時乘以較大的數值，所乘的數值一方面要滿足計算的精度，另一方面要保證最后的計算數值不超過263，這一數值是由unsigned long long int所占的字節數決定的。在疊加完成之后，再將各個矩陣元素除以這個大數。由于剛度矩陣的某些元素值較大，在采用原子函數進行操作時會產生較大的誤差，因此剛度矩陣只能按照上述方式生成。

為了減少數據在主機和設備端的傳輸耗時，質量矩陣和剛度矩陣應直接在設備端產生。由于系統只需要一次性生成質量矩陣和剛度矩陣，因此該過程的耗時在整個瞬態分析中占的比重較小，如果分析的時間較長，也可以直接在CPU 端生成剛度矩陣，然后再將數據傳輸到GPU 端。

2）節點內力的并行計算

在涉及幾何非線性的接觸碰撞問題中，節點內力fe計算公式［12］如下。

式中：幾何矩陣可分解為B0和BL兩部分，即

式中：B0是幾何矩陣中與單元節點位移無關的部分；BL是與單元節點位移相關的部分。

3）接觸搜索和接觸力的計算

接觸搜索的并行計算主要集中在全局搜索的過程中，而局部搜索主要是計算接觸間隙并確定接觸力的過程。在主從接觸面算法中，首先要指定主接觸面和從接觸面，然后按照節點數目配置線程。所有線程并行搜索與該線程對應的主擊節點的最近從節點，找出包含該從節點的距離主擊節點最近的靶單元面。形成接觸測試對后，計算主擊節點與靶單元面的距離，判斷接觸狀態并計算接觸力。

4）數據的傳輸

主機和設備數據的傳輸速率取決于PCI-E 的帶寬，PCI Express 2.0×16專門為顯卡設計，理想的傳輸速率為10Gbyte／s，而PCI Express 3.0 的設計帶寬為20Gbyte／s。但是，在有限元計算中，前后端的數據量很大。為了盡可能地減少數據傳輸耗時，一方面，將中間數據直接在顯卡中產生；另一方面，只輸出特定節點的運動狀態，并且在輸出選定節點運動狀態時間隔一定時間步長來傳輸。

4 算例分析

4.1 氣囊模型緩沖性能的仿真驗證

氣囊緩沖是航天器著陸緩沖的重要方式之一，如圖4所示。簡單氣囊模型的負載施加6m／s的初速度，振動響應點為負載的中心點。上部負載和下部地面均為剛性體，為四邊形的殼單元，其中負載總質量為50kg。氣囊為三角形膜單元，1／4的尺寸坐標可由式（6）計算，厚度為0.254 mm，彈性模量為6GPa，泊松比為0.33，密度為875kg／m3。

式中：線性坐標x和y的單位為m；θ為角度，0≤θ≤60°。

算例仿真的程序運行在同一平臺，在計算分析過程中沒有其他軟件的干擾。計算仿真過程的計算機硬件平臺如下。

由圖5～7的計算仿真結果可以看出，采用分析方法的串行計算程序和并行計算程序計算結果，與MSC.Dytran的計算結果一致，表明分析計算程序的正確性和可靠性。計算分析誤差首先來源于選取單元的類型和合理性，算例分析采用的是殼單元（Shell），而MSC.Dytran分析采用的是三角元（Tria3）等參元；另外，第2積分步長的選取也會產生計算分析結果的差異。

表1對比了本算例串行計算和并行計算的計算效率。結果表明：并行計算對串行計算的加速倍數，隨著自由度數的增加而不斷提高，提高的倍數會不斷趨近于一個峰值。

圖4 氣囊模型Fig.4 Airbag model

圖5 氣囊回收物重心處的位移響應曲線Fig.5 Displacement response curve of airbag recovery body center of gravity

圖6 氣囊回收物重心處的速度響應曲線Fig.6 Velocity response curve of airbag recovery body center of gravity

圖7 氣囊回收物重心處的加速度響應曲線Fig.7 Acceleration response curve of airbag recovery body center of gravity

表1 瞬態分析計算效率Table 1 Computation efficiency of transient analysis

4.2 典型圓柱筒接觸碰撞瞬態分析

月球軟著陸支架的主腿為一端開口的圓柱結構，如圖8所示。結構仍然采用三角形單元離散，其底部直徑為0.2 m，高度為1 m。材料的彈性模量為70GPa，泊松比為0.3，密度為2800kg／m3。假設該圓柱筒以10m／s的速度沿著垂直地面的方向與某一剛性墻相撞，計算結果如圖9～11所示。

由圖9～11可以看出：利用本文的串行算法和并行算法計算得到的薄壁筒上緣某節點的位移和速度響應曲線，基本與MSC.Dytran 的吻合，加速度響應曲線略有差別。產生差別的部分原因，與第4.1節一般瞬態分析問題一致；此外，還因為采用的接觸搜索算法和接觸力的計算方法不盡相同，加之MSC.Dytran采用變積分步長來進行數值積分。不過，利用本文算法計算的結果基本反映了接觸碰撞的規律，具備一定的可靠性。

在比較并行計算的加速倍數時，以10個迭代步的計算時間為比較標準，計算時間和加速倍數見表2。由于接觸碰撞涉及到非線性計算，每一個迭代步的節點內力都要在更新位形后進行計算，這與線性計算中提前生成剛度矩陣，然后利用剛度矩陣計算節點內力的方式不同。從表2可以看出：接觸碰撞并行計算的加速倍數峰值約為5.70；在2000個自由度以內，加速倍數基本線性增長；當自由度在2000～5000范圍時，加速倍數基本進入峰值階段。由于主體的分析計算基本上為一個線性疊加的過程，沒有涉及較為復雜的計算過程，因此，隨著自由度的繼續增加，并行計算的加速倍數基本維持在峰值附近。

圖8 典型一端開口的圓柱筒Fig.8 Typical cylinder with open end

圖9 薄壁筒上緣節點的位移響應曲線Fig.9 Displacement response curve of edge node on thin wall cylinder

圖10 薄壁筒上緣節點的速度響應曲線Fig.10 Velocity response curve of edge node on thin wall cylinder

圖11 薄壁筒上緣節點的加速度響應曲線Fig.11 Acceleration response curve of edge node on thin wall cylinder

表2 接觸碰撞的計算效率Table 2 Computation efficiency of contact impact

5 結論

（1）航天器著陸碰撞過程極為復雜，并行算法可有效提高碰撞仿真計算效率，實現航天器著陸過程中動力學計算的快速和高精度，并解決航天器著陸碰撞瞬態分析的長耗時、費資源等問題。

（2）本文將CUDA 并行計算技術應用到線性瞬態分析問題的顯式求解上，研究了質量矩陣和剛度矩陣的并行生成方法，并將基于中心差分法的求解過程并行化實施。算例表明，本文算法具有較高的可靠性，并行計算效率隨著自由度數的增加而不斷提高。

（References）

［1］Belytschko T，Neal M O.Contact-impact by the pinball algorithm with Penalty and Langrangian methods［J］.International Journal for Numerical Methods in Engineering，1991，31（3）：547-572

［2］Malone J G，Johnson N L.A parallel finite-element contact-impact algorithm for nonlinear explicit transient analysis，part 1：the search algorithm and contact mechanics［J］.International Journal for Numerical Methods in Engineering，1994，37（4）：559-590

［3］Zhong Z H，Nilsson L.Contact-impact algorithm on parallel computers［J］.Nuclear Engineering and Design，1994，150（2／3）：253-263

［4］Attaway S W，Hendrickson B A，Plimpton S J，et al.Parallel contact detection algorithm for transient solid dynamics simulations using PRONTO3D［J］.Computational Mechanies，1998，22（2）：143-159

［5］Noor A K，Lambiotte J J.Finite element dynamics analysis on CDC Star 100computer［J］.Computers ＆Structures，1979，10（1）：7-19

［6］Ou Ronfu，Fulton E.An investigation of parallel numerical integration methods for nonlinear dynamics［J］.Computers ＆Structures，1988，30（1）：403-409

［7］余天堂.一種求解結構動力響應的并行解法［J］.河海大學學報（自然科學版），1999，27（3）：75-78 Yu Tiantang.A parallel solution of dynamic response of structure［J］.Journal of Hohai University （Natural Sciences），1999，27（3）：75-78（in Chinese）

［8］Hajjar J F，Abel J F.Parallel processing of central difference transient analysis for three-dimensional nonlinear framed structures ［J］. Communications in Applied Numerical Methods，1989，5（1）：39-46

［9］Chiang K N，Fulton R E.Structural dynamics methods for concurrent processing computers［J］.Computers ＆Structures，1990，36（6）：1031-1037

［10］NVIDIA Corporation.NVIDIA＿CUDA＿Programming Guide＿2.3 ［Z／OL］.（2009-07-01）.［2013-09-02］.http：／／www.nvidia.cn／object／cuda＿develop＿cn.html

［11］金先龍，李淵印.結構動力學并行計算方法及應用［M］.北京：國防工業出版社，2008 Jin Xianlong，Li Yuanyin.The parallel algorithm and its application for structure dynamics［M］.Beijing：National Defense Industry Press，2008（in Chinese）

［12］溫文源.機械結構計算的有限元法［M］.南京：東南大學出版社，1989 Wen Wenyuan.Finite element method for mechanical structure［M］.Nanjing：Southeast University Press，1989（in Chinese）