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

（1 西北工業(yè)大學(xué)航天飛行動力學(xué)技術(shù)重點實驗室，西安 710072）（2 南京航空航天大學(xué)機械結(jié)構(gòu)力學(xué)及控制國家重點實驗室，南京 210016）

1 引言

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

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

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

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

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

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

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

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

2.2 接觸問題計算方法

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

中心差分法的穩(wěn)定性條件為

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

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

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

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

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

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

式中：M為質(zhì)量矩陣；t＋ΔtQL和t＋ΔtQc分別為系統(tǒng)在t＋Δt時刻的等效節(jié)點力向量和接觸力向量；為系統(tǒng)在t＋Δt時刻的節(jié)點內(nèi)力向量。

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

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

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

（5）計算t＋ΔtQL；

（6）由節(jié)點內(nèi)力公式fe＝kTδe計算節(jié)點內(nèi)力的初始值，然后再計算節(jié)點應(yīng)力和節(jié)點內(nèi)力，其中kT為剛度矩陣，δe為節(jié)點位移；

（7）搜尋接觸點并判斷接觸點的接觸狀態(tài)；

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

（11）令t＝t＋Δt；

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

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

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

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

圖3 并行計算數(shù)據(jù)傳遞Fig.3 Data transfer of parallel computing

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

（1）輸入相關(guān)參數(shù)并計算集中質(zhì)量矩陣。

（2）時間歷程的循環(huán)：①節(jié)點內(nèi)力的計算；②接觸搜索并計算接觸力；③根據(jù)運動微分方程求解每一迭代步的運動狀態(tài)；④根據(jù)需要輸出中間結(jié)果。

（3）結(jié)果的后處理。

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

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

1）質(zhì)量矩陣和剛度矩陣的并行生成

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

unsigned long long int atomicAdd

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

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

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）；

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

為了減少數(shù)據(jù)在主機和設(shè)備端的傳輸耗時，質(zhì)量矩陣和剛度矩陣應(yīng)直接在設(shè)備端產(chǎn)生。由于系統(tǒng)只需要一次性生成質(zhì)量矩陣和剛度矩陣，因此該過程的耗時在整個瞬態(tài)分析中占的比重較小，如果分析的時間較長，也可以直接在CPU 端生成剛度矩陣，然后再將數(shù)據(jù)傳輸?shù)紾PU 端。

2）節(jié)點內(nèi)力的并行計算

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

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

式中：B0是幾何矩陣中與單元節(jié)點位移無關(guān)的部分；BL是與單元節(jié)點位移相關(guān)的部分。

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

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

4）數(shù)據(jù)的傳輸

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

4 算例分析

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

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

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

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

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

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

圖4 氣囊模型Fig.4 Airbag model

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

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

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

表1 瞬態(tài)分析計算效率Table 1 Computation efficiency of transient analysis

4.2 典型圓柱筒接觸碰撞瞬態(tài)分析

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

（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，F(xiàn)ulton E.An investigation of parallel numerical integration methods for nonlinear dynamics［J］.Computers ＆Structures，1988，30（1）：403-409

［7］余天堂.一種求解結(jié)構(gòu)動力響應(yīng)的并行解法［J］.河海大學(xué)學(xué)報（自然科學(xué)版），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，F(xiàn)ulton 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］金先龍，李淵印.結(jié)構(gòu)動力學(xué)并行計算方法及應(yīng)用［M］.北京：國防工業(yè)出版社，2008 Jin Xianlong，Li Yuanyin.The parallel algorithm and its application for structure dynamics［M］.Beijing：National Defense Industry Press，2008（in Chinese）

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