原子氧作用下聚酰亞胺薄膜剝蝕形貌仿真方法

焦子龍，喬世英，姜海富，姜利祥，劉宇明，徐焱林，李 濤

（1.可靠性與環境工程技術重點實驗室；2.北京衛星環境工程研究所：北京 100094；3.中國航天科技集團有限公司，北京 100048）

0 引言

聚酰亞胺（PI）是性能優異的航天器用聚合物材料 ，廣泛應用于熱控材料[1-2]、柔性太陽電池基板和電路系統絕緣材料[3]等，作為帆面材料在太陽帆[4-5]和離軌帆[6]等新型空間應用方面也逐步得到重視。然而，在低地球軌道（LEO）原子氧（AO）的剝蝕作用[7-9]下，PI 表面粗糙度增大[10]造成PI 薄膜力學性能退化，可能引起PI 薄膜撕裂[11]。粗糙度增大使Kapton/Al 薄膜太陽吸收比增加，影響航天器的熱平衡[12]。粗糙度增大還可以使衛星光譜產生紅化效應，進而影響地基觀測在軌人造目標時對目標反射光譜的解析[13]。此外，對薄膜粗糙度變化所引起的透過率變化分析還是某些在軌原位測量AO通量技術的關鍵[14]。

目前對PI 薄膜剝蝕形貌的變化規律研究仍較少。Banks 等最早對PI 薄膜剝蝕形貌進行了研究[15-19]，仿真分析了在軌和地面熱等離子體試驗條件下AO 剝蝕深度標準差的變化規律，仿真設置的AO 累積通量為4.0×1020atom·cm-2。結果發現，在軌條件下該標準差隨著AO 累積通量增加而增大，但并非線性變化，且地面熱等離子體試驗條件下產生的剝蝕形貌遠不如在軌條件下的形貌粗糙[16]。Allegri 等基于LDEF 及Mir 在軌暴露試驗獲得的Kapton 薄膜透過率數據，結合仿真數據（AO 累積通量為5.3×1020atom·cm-2）發現，隨著AO 通量增加，Kapton 薄膜粗糙度先是快速增大，而后逐漸穩定在250 nm 左右[20]。有些研究人員在地面試驗中也發現粗糙度隨暴露累積通量增加而增大，但未能進一步深入研究其變化規律[12]。

因此，針對PI 薄膜表面粗糙度變化規律問題，本文改進了用于AO 掏蝕效應分析的常規計算方法，提出基于粒子運動路徑的局部網格邊界相交判斷的計算方法，并將其用于AO 作用下PI 薄膜剝蝕形貌變化模擬，研究了不同AO 束流參數下的PI 薄膜表面形貌變化，對結果進行了分析討論。

1 模擬方法

Banks 等提出了基于蒙特卡羅法的AO 與PI薄膜相互作用計算方法，并首次將其用于帶有保護涂層的PI 薄膜掏蝕形貌計算[15-19]。該方法的基本原理為：以1 個模擬粒子代表大量AO。采用二維模型，計算域以矩形網格表示。在計算域邊界隨機布置模擬粒子，速度以麥克斯韋速度分布抽樣。然后觀察模擬粒子前進軌跡上與其相交的邊界和網格，通過隨機抽樣判斷模擬粒子是否發生反應：如發生反應，則標記相應網格為已剝蝕，該粒子模擬終止，開始下一個粒子的模擬。如不發生反應，表明模擬粒子發生鏡面反射或漫反射，其反應概率發生相應改變，繼續追蹤粒子的運動軌跡，重復上述相交等過程計算；若粒子反射一定次數后仍未能發生反應，為提高計算效率，令該粒子模擬終止。

粒子與網格邊界相交計算是計算量最大的部分。為了加快計算速度，高劭倫等[21]、李濤等[22-23]、Liu 等[24]、Yang 等[25]采用四叉樹、雙鏈表等技術對粒子與網格邊界相交計算方法進行改進。但是上述改進中尋找相交的網格邊界是在整個計算域進行的，計算量仍較大。

本文提出基于粒子運動路徑的局部網格邊界相交判斷的計算方法，如圖1 所示。圖1(a)為粒子在計算域中運動的示意：計算域為長方形區域，寬度方向為x方向（向右為+x方向），厚度方向為y方向（向下為+y方向）。在模擬邊界處隨機產生模擬粒子的位置和速度。由于網格為長方形或正方形，所以很容易計算得到其所在網格。然后，根據模擬粒子運動軌跡信息和網格頂點坐標計算得到相交邊界。若相交邊界所屬網格已被剝蝕，則模擬粒子保持原速度繼續運動。若所屬網格為PI 材料，判斷模擬粒子是否與其發生反應，如不發生反應，粒子將按照鏡面反射或漫反射規律從邊界反射，繼續運動。該計算方法將AO 的運動計算僅局限于當前網格，避免在整個計算域進行搜索判斷，與傳統計算方法相比有望顯著減小計算量，提高計算效率。

圖1 局部網格邊界相交判斷的計算方法Fig.1 Calculation method for determining local grid boundary intersections

判斷上述t值中的最小正值，記為tmin。若tmin=t1，則相交邊界為①，粒子運動至當前網格左側網格；若tmin=t2，則相交邊界為②，粒子運動至當前網格右側網格；若tmin=t3，則相交邊界為③，粒子運動至當前網格上方網格；若tmin=t4，則相交邊界為④，粒子運動至當前網格下方網格。

計算機模擬中會出現如圖2 所示邊界現象，即模擬區域左右邊界處剝蝕不完全。這是由于仿真計算時左右邊界處AO 只能來自一側，使得剝蝕出現不符合實際物理規律的結果。為此，須采用周期性邊界處理方法，即以x方向模擬區域尺寸Lx為周期，將x坐標向下取整數，得到新坐標x′，

圖2 仿真計算中與實際不符的邊界現象示意Fig.2 Schematic diagram of boundary effect in simulation not matching reality

而y方向坐標及其速度保持不變。

由于該方法為統計模擬方法，需要進行多次統計抽樣以減小誤差。一般而言，統計誤差與模擬次數N的1/2 次方成反比[24]，即，過度增大模擬次數并不能有效減小統計誤差。本文將模擬次數設置為50 次。

2 計算方法驗證

采用文獻[15]中帶有保護涂層的Kapton 薄膜AO 掏蝕效應試驗結果對本文提出的計算方法進行驗證。由于Kapton 薄膜的保護涂層存在裂紋等缺陷，當AO 穿過缺陷作用于底層材料時形成掏蝕。為模擬保護涂層，將計算區域頂部一部分邊界設置為不與AO 發生反應。AO 與Kapton 材料相互作用的參數取自文獻[15]，如表1 所示。計算區域設置為深50 μm、寬50 μm，缺陷寬度為2 μm，AO 累積通量為5.77×1021atom·cm-2。

表1 仿真計算參數設置Table 1 Parameter settings for simulation calculation

計算得到的掏蝕形貌如圖3 所示，圖3(a)為本文計算結果，圖3(b)為本文計算結果與文獻[15]中的試驗結果進行對比，圖中紅色曲線為文獻[15]結果，黑色曲線為本文仿真結果。本文得到的掏蝕深度為36.5 μm，文獻[15]的結果為38.0 μm。兩者相對誤差為3.94%，表明利用本文算法對AO 與Kapton 的相互作用進行仿真的結果與試驗結果接近。

圖3 帶有保護涂層的Kapton 掏蝕深度本文計算結果與文獻[15]試驗結果對比Fig.3 Comparison of undercutting depth of protected Kapton film of calculated result in this work with experimental result in Ref.[15]

3 表面剝蝕形貌仿真分析與討論

3.1 仿真結果

針對表面形貌變化仿真，設置初始為光滑表面、邊長為100 μm 的正方形為計算區域。仿真的AO 累積通量為1.0×1021atom·cm-2，通量的仿真步長為2.0×1019atom·cm-2，共仿真50 步，得到的典型形貌變化如圖4 所示。為對比方便，將多個形貌疊加在了一起，圖中尖刺狀黑色線條表示材料本體。對 應 累 積 通 量 分 別 為：2.0×1019atom·cm-2、2.0×1020atom·cm-2、4.0×1020atom·cm-2、6.0×1020atom·cm-2、8.0×1020atom·cm-2和 1.0×1021atom·cm-2。

從圖4 可以看出，AO 作用導致了材料不斷被剝蝕，厚度逐漸減小；同時，形貌也有明顯變化，主要體現在峰和谷形貌的變化。對比對應位置的峰的變化，發現其逐漸變高并且變尖，而后逐漸變矮直至消失；谷的變化則是逐漸變深和變寬，而后又逐漸變淺和變細。仿真獲得的形貌與圖5 中Shimamura等在地面試驗中獲得的典型形貌[11]極為相似。

圖5 AO 地面試驗中Kapton 薄膜截面典型形貌[11]Fig.5 Morphology of Kapton film under AO in ground experiment[11]

統計得到的剝蝕深度隨AO 累積通量的變化如圖6 所示。從圖中可以看出，雖然理論上航天器在軌運行速度（簡稱“軌道速度”）越快，AO 具有的相對平均動能越高，剝蝕深度越大；但是由于不同軌道速度下AO 相對能量變化不大，且根據表1 中反應概率和能量的關系，可知不同軌道速度下剝蝕深度相差不大。若考慮模擬誤差的影響，可忽略不計。

統計得到了PI 薄膜表面輪廓距離平均剝蝕深度的算術平均偏差和標準差隨AO 累積通量的變化，如圖7 所示。

圖7 PI 薄膜表面輪廓距離平均剝蝕深度隨AO 累積通量的變化Fig.7 Variation of surface profile relative to the average erosion depth of PI film with AO fluence

圖中還給出了軌道速度對仿真結果的影響，模擬設定軌道速度為7.8 km·s-1、7.7 km·s-1、7.6 km·s-1、7.5 km·s-1和7.4 km·s-1。從圖中可以看出：不同軌道速度下，表面輪廓距離平均剝蝕深度變化規律基本相同；AO 累積通量較小時，表面輪廓距離平均剝蝕深度的算術平均偏差和標準差迅速增大。隨著AO 累積通量增大，兩者變化逐漸平緩。軌道速度越高，兩者數值越大，但相差不大。這主要是因為AO 能量相差不大，因而與PI 反應特性基本相同。算術平均偏差和標準差變化趨勢相同。比較圖7(a)與圖7(b)的數值，發現兩者對應AO 累積通量下的比例固定，約為1.25。這說明PI 薄膜的表面形貌特征基本未變，只是表面峰/谷數值變大。

將上述計算結果與文獻[16]中模擬軌道速度7.8 km·s-1的算術平均偏差結果進行了比較，結果如圖8 所示。可以看出，AO 累積通量大于等于8.0×1019atom·cm-2時符合較好。

圖8 本文算術平均偏差與文獻[16]對比Fig.8 Comparison between the arithmetic mean deviation in this paper and Ref.[16]

針對軌道速度為7.8 km·s-1的PI 薄膜形貌參數計算結果進行其與AO 累積通量關系擬合，結果如表2 所示。

表2 PI 薄膜形貌參數與AO 累積通量關系擬合結果Table 2 Fitting results of PI film morphology parameters with AO fluence

根據擬合結果可以看出，形貌參數以AO 累積通量的0.253 次冪規律增大。

3.2 討論

目前對AO 剝蝕形貌微觀機理和形貌參數長期變化趨勢的研究仍存在不足。本文結果是基于文獻中AO 與PI 材料作用參數（如反應概率及其隨入射角度和能量的變化、熱同化概率等）計算獲得的，但上述參數是通過掏蝕形貌計算擬合在軌實驗數據所得，其合理性尚需進一步分析。

本研究中采用矩形網格表示材料單元的合理性仍有待進一步研究，也有必要探討采用其他幾何結構，如三角網格等進行剝蝕形貌計算分析——三角形網格可能更適合表示起伏的峰谷形貌，但編程實施難度會顯著增加。另外，有必要進行三維形貌的計算分析和試驗表征，因為其可對二維形貌計算中采用的作用參數進行驗證。

目前系統的剝蝕形貌定量研究仍較少，無法對仿真結果進行深入對比驗證。因此表面形貌的測量和表征試驗對于AO 剝蝕形貌研究極為重要。

4 結束語

為進一步研究PI 薄膜在AO 束流作用下的表面剝蝕形貌，參考現有蒙特卡羅掏蝕形貌仿真方法，提出了一種基于局部網格邊界相交判斷的方法，并應用周期性邊界處理方法獲得符合實際物理規律的表面形貌。采用本文方法仿真了已有的掏蝕形貌算例，計算結果與文獻[15]中的試驗結果相差小于4%，證明本文計算方法能夠正確描述AO 與PI 材料的相互作用。

采用剝蝕形貌與平均剝蝕深度的算術平均偏差和標準差作為材料形貌特征的描述參數。結果發現，算術平均偏差和標準差均隨著AO 累積通量增大而以0.253 次冪律增大；標準差和算術平均偏差的比值隨AO 累積通量增大基本保持不變，說明表面形貌特征在AO 作用下基本不變，只是表面峰值和谷值增大。此外，不同軌道速度條件下表面剝蝕形貌差異不大。

本文的仿真分析方法和計算結果將有助于更好地理解AO 作用下航天器用聚合物表面剝蝕形貌形成機理，為PI 等航天用聚合物薄膜材料更好的空間應用提供保障。