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

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

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

0 引言

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

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

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

1 模擬方法

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

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

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

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

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

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

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

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

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

2 計算方法驗證

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

表1 仿真計算參數(shù)設(shè)置Table 1 Parameter settings for simulation calculation

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

圖3 帶有保護涂層的Kapton 掏蝕深度本文計算結(jié)果與文獻[15]試驗結(jié)果對比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 仿真結(jié)果

針對表面形貌變化仿真，設(shè)置初始為光滑表面、邊長為100 μm 的正方形為計算區(qū)域。仿真的AO 累積通量為1.0×1021atom·cm-2，通量的仿真步長為2.0×1019atom·cm-2，共仿真50 步，得到的典型形貌變化如圖4 所示。為對比方便，將多個形貌疊加在了一起，圖中尖刺狀黑色線條表示材料本體。對 應(yīng) 累 積 通 量 分 別 為：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 作用導(dǎo)致了材料不斷被剝蝕，厚度逐漸減小；同時，形貌也有明顯變化，主要體現(xiàn)在峰和谷形貌的變化。對比對應(yīng)位置的峰的變化，發(fā)現(xiàn)其逐漸變高并且變尖，而后逐漸變矮直至消失；谷的變化則是逐漸變深和變寬，而后又逐漸變淺和變細。仿真獲得的形貌與圖5 中Shimamura等在地面試驗中獲得的典型形貌[11]極為相似。

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

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

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

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

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

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

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

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

表2 PI 薄膜形貌參數(shù)與AO 累積通量關(guān)系擬合結(jié)果Table 2 Fitting results of PI film morphology parameters with AO fluence

根據(jù)擬合結(jié)果可以看出，形貌參數(shù)以AO 累積通量的0.253 次冪規(guī)律增大。

3.2 討論

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

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

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

4 結(jié)束語

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

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

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