張曉天, 賈光輝

(北京航空航天大學(xué) 宇航學(xué)院, 北京 100191)

0 引 言

航天器解體碎片生成研究是空間碎片環(huán)境建模的關(guān)鍵。目前的主要研究手段是對(duì)解體碎片數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，建立經(jīng)驗(yàn)公式。解體碎片數(shù)據(jù)來(lái)源于兩個(gè)方面：軌道碎片雷達(dá)觀測(cè)和地面解體實(shí)驗(yàn)。

美國(guó)太空觀測(cè)網(wǎng)(SSN，Space Surveillance Network)對(duì)地球軌道的物體進(jìn)行跟蹤和編目。以目前的觀測(cè)能力，能夠?qū)Φ偷厍蜍壍捞卣鏖L(zhǎng)度10cm以上的物體和高地球軌道特征長(zhǎng)度1m以上的物體進(jìn)行可靠跟蹤，難以獲得小尺寸碎片的數(shù)據(jù)。

在空間碎片群體中，小碎片的數(shù)目要遠(yuǎn)大于大碎片。為了研究小碎片的特性，NASA 開(kāi)展了SOCIT(Satellite Orbital Debris Characterization Impact Test)地面解體實(shí)驗(yàn)[1]。解體實(shí)驗(yàn)中使用了OSCAR 小衛(wèi)星，通過(guò)加速鋁彈丸對(duì)其碰撞造成解體。日本T. Hanada 研究小組與NASA 約翰遜研究中心合作，開(kāi)展了復(fù)合材料小衛(wèi)星地面解體實(shí)驗(yàn)研究[2]。柳森、黃潔、蘭勝威和李毅等開(kāi)展了小型衛(wèi)星模型的解體實(shí)驗(yàn)研究[3-5]，進(jìn)行了三種衛(wèi)星模型的撞擊解體實(shí)驗(yàn)，并提出了新的解體閾值。王海福、余慶波和金學(xué)科等開(kāi)展了圓柱殼結(jié)構(gòu)爆炸解體實(shí)驗(yàn)和數(shù)值仿真研究[6-7]。地面解體實(shí)驗(yàn)受到加速設(shè)備能力限制，難以構(gòu)造常規(guī)尺寸衛(wèi)星解體的碰撞條件，地面解體實(shí)驗(yàn)難度很大。

本文通過(guò)超高速碰撞數(shù)值模擬技術(shù)對(duì)航天器解體過(guò)程進(jìn)行研究，從而獲得解體碎片的特性。數(shù)值模擬技術(shù)是一種相對(duì)成本低的技術(shù)，一旦建立了有效的模擬算法，就可以大量開(kāi)展。由于目前廣泛使用的SPH方法在單個(gè)碎片的識(shí)別和置信度方面存在不足，本文通過(guò)一種SPH和有限元結(jié)合的方法來(lái)實(shí)現(xiàn)對(duì)碎片的識(shí)別和統(tǒng)計(jì)，并與衛(wèi)星模型地面解體實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。

1 衛(wèi)星解體仿真方法

目前的超高速撞擊數(shù)值模擬研究主要針對(duì)航天器防護(hù)結(jié)構(gòu)分析，廣泛使用的數(shù)值模擬方法是SPH(Smoothed Particle Hydrodynamics)方法[8]。但在航天器防護(hù)結(jié)構(gòu)仿真中，由于撞擊體較小，航天器一般不會(huì)發(fā)生解體；另外這種分析的關(guān)注點(diǎn)為防護(hù)結(jié)構(gòu)是否被碎片穿透，并不關(guān)注撞擊產(chǎn)生的二次碎片云中單個(gè)碎片的信息(如碎片數(shù)目、大小和速度矢量等)。而解體模型就是要根據(jù)解體條件給出解體碎片特征信息，解體碎片特性是解體動(dòng)力學(xué)數(shù)值模擬最關(guān)心的結(jié)果。因此，將現(xiàn)有數(shù)值模擬算法成果應(yīng)用于航天器解體動(dòng)力學(xué)數(shù)值模擬，首先需要對(duì)單個(gè)碎片的識(shí)別性和置信度等方面進(jìn)行分析。

對(duì)于超高速碰撞問(wèn)題，SPH方法的仿真結(jié)果常表現(xiàn)為離散的粒子云，難以從中識(shí)別出哪些粒子構(gòu)成一個(gè)碎片，也難以分辨破碎材料和未破碎材料的界面。雖然可以通過(guò)搜索算法直接將結(jié)果中鄰近的粒子劃分為一個(gè)碎片，但是這種方法存在幾個(gè)方面的不足。第一，置信度問(wèn)題。在初始模型中，粒子被以一定間距δ填充入物體所占據(jù)的空間。此時(shí)粒子作為1個(gè)群體代表了1個(gè)物體，其可信性較高。但是單個(gè)粒子代表了1塊δ立方的物質(zhì)微元，其作為孤立個(gè)體可信性較低。因?yàn)樵撐镔|(zhì)微元有可能進(jìn)一步分裂破碎，而一個(gè)粒子不可能再被分裂。所以仿真結(jié)果中的孤立粒子，甚至是少量幾個(gè)粒子都不能作為置信的單個(gè)碎片存在。第二，流動(dòng)性問(wèn)題。超高速碰撞過(guò)程中，物質(zhì)流動(dòng)范圍很大，初始不相鄰的物質(zhì)微元完全有可能在撞擊過(guò)程中被“擠”到一起。因此僅僅根據(jù)某個(gè)時(shí)間點(diǎn)上粒子的空間坐標(biāo)是否相鄰來(lái)判斷這兩個(gè)物質(zhì)微元是否屬于同一個(gè)碎片是不夠可靠的。第三，鄰近閾值問(wèn)題。判斷粒子是否鄰近，一般是判斷兩者間距是否小于給定閾值。這個(gè)閾值的給定缺乏力學(xué)依據(jù)，難以制定統(tǒng)一的規(guī)則，這會(huì)影響不同問(wèn)題之間的一致性。

SPH方法在航天器解體碎片仿真問(wèn)題中存在上述不足；另一方面有限元方法材料界面清晰，碎片效果直觀，單個(gè)碎片置信度高，但是在大變形問(wèn)題中會(huì)發(fā)生網(wǎng)格畸變，使得計(jì)算無(wú)法進(jìn)行。將SPH與有限元方法相結(jié)合的方法，稱為有限元重構(gòu)方法，提高了對(duì)單個(gè)碎片的識(shí)別性和置信度。在有限元重構(gòu)方法中首先通過(guò)立方體填充，將任意形狀的幾何體離散為粒子-形心重合模型；進(jìn)而通過(guò)顯式積分迭代步中的SPH積分、有限元重構(gòu)與失效分析三個(gè)步驟，獲得包含F(xiàn)E單元與SPH粒子的混合仿真結(jié)果。

1.1粒子-形心重合模型

圖1是一個(gè)離散化的實(shí)例。A代表一個(gè)任意形狀的物體，B表示使用立方體FE單元對(duì)物體區(qū)域進(jìn)行填充，C表示在每個(gè)單元的形心處創(chuàng)建一個(gè)SPH粒子，相應(yīng)的粒子間距等于FE的單元邊長(zhǎng)。此時(shí)的模型稱為“粒子-形心重合模型”，包括FE網(wǎng)格與SPH粒子兩部分。

A B C D E

1.2顯式積分迭代步

對(duì)于粒子-形心重合模型中的粒子部分，按照SPH控制方程進(jìn)行顯示時(shí)間積分迭代，在每個(gè)時(shí)間步后能夠獲得所有粒子的位置、速度信息。

AEi={EH|Ni∈EH}

(1)

在立方體填充模型中，根據(jù)節(jié)點(diǎn)被多少個(gè)單元所共有——也就是AEi集合中單元的個(gè)數(shù)——共分為8種情況。如果1個(gè)節(jié)點(diǎn)被8個(gè)單元共有，則稱其為內(nèi)點(diǎn)，否則稱其為邊界點(diǎn)。有限元重構(gòu)分為3個(gè)步驟依次執(zhí)行：內(nèi)點(diǎn)重構(gòu)、第一類(lèi)邊界點(diǎn)重構(gòu)，第二類(lèi)邊界點(diǎn)重構(gòu)。

(1)內(nèi)點(diǎn)重構(gòu)

用IN表示立方體填充網(wǎng)格中的內(nèi)點(diǎn)集合，它包含了所有被8個(gè)單元共享的節(jié)點(diǎn)。內(nèi)點(diǎn)的重構(gòu)坐標(biāo)定義為這8個(gè)單元對(duì)應(yīng)的SPH粒子當(dāng)前坐標(biāo)的算術(shù)平均：

(2)

其中，pH表示EH形心對(duì)應(yīng)的SPH粒子的坐標(biāo)向量，ri表示Ni的坐標(biāo)向量。

(2)第一類(lèi)邊界點(diǎn)重構(gòu)

定義一個(gè)與Ni關(guān)聯(lián)的單元標(biāo)識(shí)集合：

(3)

第一類(lèi)邊界點(diǎn)集合定義為：

(4)

第一類(lèi)邊界點(diǎn)的重構(gòu)坐標(biāo)按下式進(jìn)行計(jì)算：

(5)

(3)第二類(lèi)邊界點(diǎn)重構(gòu)

定義一個(gè)與Ni關(guān)聯(lián)的單元標(biāo)識(shí)集合：

(6)

第二類(lèi)邊界點(diǎn)集合定義為：

(7)

第二類(lèi)邊界點(diǎn)的坐標(biāo)按下式進(jìn)行計(jì)算：

(8)

如果一個(gè)邊界點(diǎn)所有的背點(diǎn)都是第二類(lèi)邊界點(diǎn)，將陷入一個(gè)死循環(huán)，這些點(diǎn)無(wú)法被更新。為避免這種極端情況，在幾何體離散的過(guò)程中需要對(duì)最薄的部位采取兩層以上的單元來(lái)填充。

在大變形問(wèn)題中，SPH粒子流動(dòng)范圍很大，重構(gòu)后的單元不一定能維持常規(guī)的形狀，很多單元會(huì)發(fā)生畸變，此時(shí)失效單元將被刪除，單元形心對(duì)應(yīng)的SPH粒子被保留；未失效單元被保留，對(duì)應(yīng)的SPH粒子被刪除。本文采用Johnson-Cook失效準(zhǔn)則。FE重構(gòu)分析和失效分析按照SPH計(jì)算的時(shí)間步逐步進(jìn)行，第n+1個(gè)時(shí)間步的分析繼承第n個(gè)時(shí)間步的結(jié)果。在結(jié)果中FE單元代表未失效材料，SPH粒子代表失效材料。

2 兩體碰撞解體算例

使用上述方法對(duì)實(shí)心球體撞擊實(shí)心立方體解體問(wèn)題進(jìn)行仿真，兩物體相對(duì)位置關(guān)系如圖2所示，從右上角局部放大圖中可以看出粒子和單元形心的重合關(guān)系。球體的直徑為1m，立方體邊長(zhǎng)1m。兩物體的材料密度均為0.2g/cm3，球體質(zhì)量為104.72kg，立方體質(zhì)量為200kg。圖中的撞擊速度v=7.6km/s，相對(duì)撞擊速度為15.2km/s。建模中有限元單元邊長(zhǎng)定為2cm，因此在仿真分析中可以得到邊長(zhǎng)大于2cm的置信碎片結(jié)果。

圖2 球體與立方體碰撞解體仿真模型

圖3給出了350μs時(shí)刻碎片云的對(duì)比。左圖給出了SPH仿真結(jié)果，從圖中可以看出，SPH的結(jié)果表現(xiàn)為粒子云。由局部放大圖中難以分辨出哪些粒子構(gòu)成一個(gè)碎片，破碎材料和未破碎材料的界面難以清晰定義。中圖為有限元重構(gòu)的結(jié)果，結(jié)果由SPH粒子和有限元單元混合表達(dá)。在局部放大圖中，由單元表達(dá)的單個(gè)碎片清晰可見(jiàn)。右圖為有限元重構(gòu)結(jié)果中的單元部分，也就是置信碎片部分。SPH粒子代表了充分破碎的材料，作為整體表達(dá)了碎片云的形狀、演化與材料分布，但是作為單個(gè)碎片是不夠可信的。

圖3 球體與立方體碰撞解體數(shù)值模擬結(jié)果

固定單元是否失效的狀態(tài)不變，將節(jié)點(diǎn)和粒子坐標(biāo)反推回0時(shí)刻，就得到了反推圖。圖4所示為350μs的反推圖。其中由單元表示的部分在撞擊過(guò)程結(jié)束后會(huì)轉(zhuǎn)變?yōu)橹眯潘槠?，或者說(shuō)是較大的碎片；由粒子表示的部分在撞擊結(jié)束后會(huì)轉(zhuǎn)變?yōu)榉侵眯潘槠?，這些部分會(huì)充分破碎為小碎片云。由圖中可見(jiàn)，撞擊區(qū)域內(nèi)的材料大部分都充分破碎了，大碎片主要由遠(yuǎn)離撞擊區(qū)的材料形成。

圖4 球體與立方體碰撞結(jié)果反推圖

有限元重構(gòu)方法識(shí)別了SPH仿真結(jié)果中的碎片，確定了清晰的破碎與未破碎材料的界面，對(duì)SPH仿真結(jié)果中單個(gè)碎片的置信部分和非置信部分進(jìn)行了清晰的劃分。首先，在解體過(guò)程中如果單元失效，它將被轉(zhuǎn)換為粒子；而粒子繼續(xù)受變形破裂，就無(wú)法描述，因?yàn)槠洳豢稍俜?。因此在解體數(shù)值模擬結(jié)果中，以單元形式存在的物質(zhì)就構(gòu)成了置信碎片；而以粒子形式存在的失效單元?jiǎng)t屬于非置信碎片。對(duì)于單個(gè)碎片的分析應(yīng)該基于置信碎片或者說(shuō)基于單元部分進(jìn)行。第二，由于引入了單元，刻畫(huà)了單元之間的物質(zhì)連通性，初始相鄰的單元如果解體結(jié)束后仍然相鄰則屬于同一個(gè)碎片；如果初始不相鄰，而在解體結(jié)束后相鄰，則說(shuō)明它們的鄰近是由于材料的運(yùn)動(dòng)引起的，因此不構(gòu)成同一個(gè)碎片。這樣，哪些單元構(gòu)成一個(gè)碎片將由單元部分的拓?fù)浣Y(jié)構(gòu)決定。第三，單元轉(zhuǎn)換粒子判據(jù)是材料失效或相變，而不是人為指定幾何閾值，因此置信與否的分界線從根本上是由力學(xué)機(jī)理決定的。

3 碎片統(tǒng)計(jì)方法

在有限元重構(gòu)的仿真結(jié)果中，置信碎片由單元部分表達(dá)，碎片的總數(shù)以及碎片的單元構(gòu)成是由單元部分的拓?fù)浣Y(jié)構(gòu)決定的。碎片統(tǒng)計(jì)就是要設(shè)計(jì)一種自動(dòng)的算法，能夠統(tǒng)計(jì)出碎片的數(shù)目、碎片的單元構(gòu)成以及每個(gè)碎片的尺寸、質(zhì)量和速度矢量等信息。

把置信部分的單元看作頂點(diǎn)，如果兩個(gè)單元之間存在一個(gè)公共節(jié)點(diǎn)，則在兩單元之間定義一條邊，這樣頂點(diǎn)和邊的集合可以表示為：

(9)

〈Vtx,Edg〉構(gòu)成一個(gè)無(wú)向圖。如果該圖中存在一個(gè)首尾相接的邊序列，能夠把圖中兩個(gè)頂點(diǎn)連接起來(lái)，則這兩個(gè)頂點(diǎn)是連通的。如果〈Vtx,Edg〉的一個(gè)子集中任意兩個(gè)頂點(diǎn)是連通的，則稱這個(gè)子集為該無(wú)向圖的一個(gè)連通分支。每一個(gè)連通分支就對(duì)應(yīng)了碎片云中的一個(gè)碎片。因此，碎片的統(tǒng)計(jì)與特征量計(jì)算就是對(duì)連通分支的分析。

本文使用了基于廣度優(yōu)先的連通分支遍歷算法，時(shí)間復(fù)雜度為o(n2)。對(duì)圖4中的碎片進(jìn)行分析，圖5給出了連通分支搜尋的結(jié)果，其中每一種顏色代表了一個(gè)連通分支，也就是一個(gè)碎片。該算例中共產(chǎn)生了大于2cm的碎片562個(gè)，其中178來(lái)自球體，384個(gè)來(lái)自立方體。

在獲得了碎片的單元構(gòu)成之后，可以對(duì)每個(gè)碎片的特征長(zhǎng)度進(jìn)行計(jì)算：

(10)

圖5 碎片統(tǒng)計(jì)結(jié)果

4 衛(wèi)星模型解體碎片數(shù)值模擬

利用上述方法對(duì)文獻(xiàn)[3]開(kāi)展的地面衛(wèi)星模型撞擊解體實(shí)驗(yàn)進(jìn)行數(shù)值模擬。

圖6給出了衛(wèi)星解體仿真模型，圖7為Case B的仿真結(jié)果。從圖7中可以看出，雖然彈丸尺寸較小，但是造成了很大的破壞。由于衛(wèi)星模型在撞擊法平面內(nèi)的兩個(gè)軸向并不對(duì)稱，所以結(jié)構(gòu)中的開(kāi)孔也不對(duì)稱，而是呈橢圓形。仿真結(jié)果中，衛(wèi)星結(jié)構(gòu)并沒(méi)有像實(shí)驗(yàn)中一樣板結(jié)構(gòu)之間互相脫離。這是因?yàn)樵趯?shí)驗(yàn)結(jié)構(gòu)中，相鄰板結(jié)構(gòu)采用了L型鈑金配合螺釘連接，連接強(qiáng)度較弱；在仿真中，連接形式進(jìn)行了近似處理，直接采用了相鄰邊焊接的方式。

圖6 衛(wèi)星模型撞擊解體仿真建模

圖8給出了Case A和Case B的累計(jì)碎片數(shù)目對(duì)碎片尺寸的分布曲線對(duì)比?？傮w來(lái)說(shuō)仿真結(jié)果與實(shí)驗(yàn)數(shù)據(jù)比較一致。仿真曲線橫坐標(biāo)值在0.001以下的缺失是由模型分辨率引起的，為了獲得更小碎片的特性，需要進(jìn)一步提高模型分辨率。另外，在小尺寸碎片區(qū)域，NASA模型的計(jì)算值明顯高于實(shí)驗(yàn)值，這可能是由于NASA模型的建模數(shù)據(jù)源與文獻(xiàn)[3]實(shí)驗(yàn)的差異造成的。NASA模型小碎片特性部分的建模數(shù)據(jù)主要來(lái)源于SOCITE4解體實(shí)驗(yàn)，該實(shí)驗(yàn)中用到的OSCAR衛(wèi)星的內(nèi)部結(jié)構(gòu)與材料和文獻(xiàn)[3]實(shí)驗(yàn)中用的衛(wèi)星模型相差較大。這也體現(xiàn)了衛(wèi)星的結(jié)構(gòu)和材料對(duì)解體碎片特性的影響，表明了基于解體動(dòng)力學(xué)過(guò)程分析、計(jì)算碎片特性的必要性。

圖7 Case B數(shù)值模擬結(jié)果

(a) Case A

(b) Case B

5 結(jié) 論

本文基于超高速碰撞數(shù)值模擬，建立了一種用于模擬衛(wèi)星解體碎片生成的有限元重構(gòu)方法，對(duì)衛(wèi)星模型解體問(wèn)題進(jìn)行了仿真研究，得到如下結(jié)論：

(1) 有限元重構(gòu)方法能夠?qū)πl(wèi)星解體生成的碎片進(jìn)行模擬，結(jié)合無(wú)向圖連通域分析方法，能夠獲得碎片分布以及單個(gè)碎片的數(shù)據(jù)；

(2) 對(duì)地面衛(wèi)星模型解體實(shí)驗(yàn)進(jìn)行了數(shù)值模擬分析。對(duì)于Case A和Case B兩種工況，模擬獲得的累計(jì)碎片數(shù)對(duì)碎片尺寸的分布曲線與實(shí)驗(yàn)曲線基本一致。表明了本文方法的有效性。另外，通過(guò)和NASA模型計(jì)算結(jié)果對(duì)比，表明解體衛(wèi)星的結(jié)構(gòu)形式和材料對(duì)解體碎片特性有較大影響，通過(guò)對(duì)解體動(dòng)力學(xué)過(guò)程的研究，可以更加有效的計(jì)算解體碎片特性。

參考文獻(xiàn):

[1]Krisko P H, Horstman M, Fudge M L. SOCIT4 collisional-breakup test data analysis: With shape and materials characterization[J]. Advances in Space Research, 2008, 41(7): 1138-1146.

[2]Hanada T, Liou J C, Nakajima T, et al. Outcome of recent satellite impact experiments[J]. Advances in Space Research, 2009, 44(5): 558-567.

[3]Lan S, Li Y, Huang J, et al. Fragment characteristic of simulated spacecraft under hypervelociy impact[C]//Proceeding of 62nd International Astronautical Congress. Cape Town, SA, 2011.

[4]蘭勝威, 李毅, 黃潔, 等. 模擬航天器超高速撞擊解體碎片特性研究[C]//第六屆全國(guó)空間碎片學(xué)術(shù)交流會(huì)論文集, 四川成都, 2011.

Lan Shengwei, Li Yi, Huang Jie, et al. Investigation on the characteristics of fragments generated by the simulated spacecraft hypervelocity impact test[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[5]李毅, 黃潔, 馬兆俠, 等. 一種新的衛(wèi)星超高速撞擊解體閾值模型研究[C]//第六屆全國(guó)空間碎片學(xué)術(shù)交流會(huì)論文集. 四川成都, 2011.

Li Yi, Huang Jie, Ma Zhaoxia, et al. Investigation on new

breakup threshold model of satellite breakup under hypervelocity impact[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[6]余慶波, 金學(xué)科, 王勤智, 等. 結(jié)構(gòu)強(qiáng)度對(duì)衛(wèi)星爆炸解體影響研究[C]//第六屆全國(guó)空間碎片學(xué)術(shù)交流會(huì)論文集. 四川成都, 2011.

Yu Qingbo, Jin Xueke, Wang Qinzhi, et al. The influence of structure strength on satellite breakup caused by explosion[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[7]金學(xué)科, 余慶波, 王勤智, 等. 燃料量對(duì)衛(wèi)星爆炸解體影響研究[C]//第六屆全國(guó)空間碎片學(xué)術(shù)交流會(huì)論文集. 四川成都, 2011.

Jin Xueke, Yu Qingbo, Wang Qinzhi, et al. The influence of the amount of fuel on satellite breakup caused by explosion[C]//Proceeding of 6thNational conference of space debris. Chengdu, Sichuan, 2011.

[8]Johnson G R, Stryk R A, Beissel S R. SPH for high velocity impact computations[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 57(4): 573-580.

作者簡(jiǎn)介:

張曉天(1984-)，男，山東濟(jì)南人，博士，講師。研究方向：高速碰撞動(dòng)力學(xué)，航天器防護(hù)結(jié)構(gòu)設(shè)計(jì)，空間碎片與空間環(huán)境。通信地址：北京航空航天大學(xué)宇航學(xué)院航天飛行器技術(shù)系(100191)。E-mail: zhangxiaotian@buaa.edu.cn