防雷/防熱涂層的雷擊熱解損傷的跨尺度模擬

董 琪，趙致遠，李 婷，賈玉璽，張如良

（1.山東科技大學，青島，266590；2.航天材料及工藝研究所，北京，100076；3.山東大學，濟南，250061）

0 引言

運載火箭、導彈等航天器飛行將面臨超高溫熱流環境，其鋁合金或碳纖維復合材料殼體不能直接裸露在外，需涂覆硅橡膠類燒蝕防熱涂層［1］，以阻止高溫破壞殼體并向內傳遞。而目前防熱涂層通常采用絕緣材料，在雷電環境中飛行會面臨雷擊破壞，致使航天器在后續飛行中失去防熱保護，導致任務失敗。因此，開展防熱涂層雷擊損傷研究十分必要。

雷電作用下防雷/防熱涂層產生的直接效應具有顯著的電-熱-化學-力多物理場耦合特征，涉及涂層材料的絕緣擊穿、電阻生熱、電弧燒蝕、樹脂熱解及陶瓷化以及力學破壞等。然而，雷擊試驗的成本高、危險性高、復雜瞬態過程難以表征，因此，使用計算機進行模擬計算成為目前雷擊損傷研究的重要手段。其中，以有限元模擬為代表，能夠考慮雷擊損傷復雜的電-熱-化學-力多物理場耦合特征，量化不同耦合場作用下不同類型的損傷情況［2］，可以有效闡明防雷/防熱涂層雷擊損傷的機理，為雷擊損傷預測提供理論基礎。Dong 等［3］建立了碳纖維增強復合材料（Carbon Fiber Reinforced Plastic，CFRP）雷擊損傷過程的電-熱-化學分析模型并加以模擬分析；Alexandros 等［4］建立了CFRP 層合板雷擊損傷過程的電-熱耦合模型，對比分析了層合板由螺栓連接和粘貼連接時的雷擊損傷情況；Zhu 等［5］則利用CFRP 復合材料雷電損傷有限元模型討論了其幾何形狀、纖維排布順序等因素對損傷程度的影響。

但以有限元為代表的多物理場耦合技術由于時間尺度的限制，在模擬涂層在雷擊熱效應下發生熱解這一過程時尚存在不足，采用的電-熱-化學耦合有限元模型計算熱解損傷時缺乏可靠的熱解反應動力學方程。這是由于熱解反應動力學方程的獲取通常來源于對材料熱重試驗曲線的熱分析，而熱重試驗條件遠達不到實際雷擊所造成的極端升溫速率和極端高溫情況，所得的動力學方程不能很好地與實際吻合。因此，為更合理地研究防雷/防熱涂層雷擊損傷機理，進而準確預測其雷擊燒蝕損傷，探明涂層材料在雷擊極端溫度條件下的熱解動力學機制及行為具有重要意義。

目前針對此類問題的研究較為少見，然而，針對這一問題，分子模擬手段能夠從原子-分子的微觀尺度模擬聚合物分子在極端升溫速率和極端高溫下的熱解行為。Zhang等［6］利用分子模擬方法對PVC和巴塞木進行了熱解模擬，并分析了模擬結果與雷擊損傷的關系；Paajanen等［7］模擬了高溫下纖維素的熱解并分析了其反應機理與動力學。

本文采用反應分子動力學模擬方法研究甲基苯基硅橡膠在雷電極端升溫速率下的熱解反應機制，進而通過雷擊過程電-熱-化學耦合有限元模擬，跨尺度計算甲基苯基硅橡膠基防雷/防熱涂層材料體系的雷擊熱解損傷，并將模擬結果與雷擊試驗損傷結果作對比，旨在構建能更準確反映涂層雷擊實際過程的跨尺度模擬方法，為涂層雷擊熱解損傷預測以及防熱層防雷材料開發提供理論支持。

1 基礎理論與模擬方法

1.1 涂層樹脂熱解反應分子動力學模擬

選取甲基苯基硅橡膠的交聯網絡狀結構作為研究體系，其分子的通用結構式如圖1所示。

苯基為硅橡膠提供卓越的耐熱、耐低溫及耐輻照性能，本文以30%苯基含量構建甲基苯基硅橡膠分子鏈模型作為交聯前的單體模型，如圖2所示。

圖2 甲基苯基硅橡膠的單體分子模型Fig.2 Molecular model of methyl phenyl silicone rubber monomer

將上述16 個硅橡膠單體分子鏈裝入大小為30.4 ?×30.4 ?×30.4 ? 的具有用周期性邊界條件的立方盒子中，構建交聯體系，初始建模密度為1 g/cm3。根據氧化脫氫的硫化反應法，使用交聯腳本程序［8］代替催化劑的作用以簡化交聯過程，將硅橡膠鏈上甲基的C 原子標記為反應位點直接進行交聯，形成C—C 鍵的同時刪除甲基上的一個氫原子。模擬交聯過程采用COMPASS Ⅱ力場，目標交聯度為100%。為增加分子的移動性，設置溫度為500 K，最終得到交聯度為95.1%的硅橡膠交聯結構模型，如圖3所示。

圖3 甲基苯基硅橡膠交聯分子模型Fig.3 Molecular model of cross-linked methyl phenyl silicone rubber

基于上述硅橡膠分子模型，本文通過Lammps 2020 軟件，采用Adri van Duin 等［9］開發的ReaxFF 反應力場模擬硅橡膠的熱解過程。

在實際的雷擊過程中，復合材料雷擊作用點處升溫速率可達1011K/min以上，極短時間內可達3 000 K以上的極高溫度［10］。為使熱解反應模擬的條件盡可能接近雷擊的實際情況，同時兼顧計算成本和計算機運行效能，本文進行了升溫速率位于1012～1014K/s 區間范圍內的多次熱解反應分子模擬，且將模擬的終止溫度設置在7 000 K以上，以保證分子模型完全熱解，具體算例的溫度設置條件如表1所示。此外，模擬在NPT系綜下進行，設置步長為0.25 fs。

表1 不同極端升溫速率下硅橡膠熱解的分子模擬條件Tab.1 Molecular simulation conditions for the pyrolysis of silicone rubber under different extreme heating rates

根據數據處理需要，本文設置每100步輸出一次當前時刻模擬系統內產物信息，該信息包括當前時刻下系統中所有存在的產物的分子式及其相應數量。

1.2 樹脂熱解反應分子模擬的動力學方程

基于上述硅橡膠熱解反應分子模擬輸出的不同時刻產物信息，本文根據式（1）計算不同時刻下固相產物所占的質量比：

式中W為質量保留率；x為當前體系中存在的化合物種類總數量；N0為反應開始時體系的微粒數；M0為反應開始時體系的相對分子質量；Nti為t時刻下第i種剩余固相產物的微粒數；Mti為t時刻下第i種剩余固相產物的相對分子質量。

由此可得到硅橡膠熱解反應的虛擬熱解度曲線。通過不同升溫速率下的虛擬熱解度曲線，進一步計算動力學參數（活化能Ea、指前因子A、反應級數n）。本文采用的熱解動力學參數計算方法如下：

a）Starink法。

Starink［11］對比了熱分析動力學方法中經典的Kissinger 方程、Ozawa 方程和Boswell 方程的計算精度，并提出了較高計算精度的Starink方程，其具體形式如下：

式中β為升溫速率；T為熱力學溫度；Rg為普適氣體常量；G(α)為反應機理函數的積分式。

根據該方程進行擬合可計算出活化能Ea，且在確定機理函數的前提下，可進一步得到指前因子A。

b）Malek法。

硅橡膠的熱解過程是聚合物受熱解聚、無規則斷鏈的結構轉變，是分子量變小的反應，符合n級化學反應機理，Malek法［12］能夠在得到活化能Ea的基礎上計算反應級數n。對于n級化學反應機理，Malek法計算反應級數n的具體形式如下：

式中αp為熱解速率最大時的熱解度；up表示Ea/RgTp；π(up)為Luke近似式［13］中引入的有理函數。

將αp、Tp、Ea、Rg代入式（3）中迭代可求出反應級數n。

在得到活化能Ea、指前因子A、反應級數n后，代入式（4）中，即得到熱解反應分子模擬動力學方程。

1.3 防雷/防熱涂層雷擊熱解損傷的有限元模擬

首先，構建與雷擊試驗中試樣一致的涂層平板幾何模型，如圖4所示。模型由三層材料構成，表層為130 μm的防雷層、次表層為4 mm的防熱層、底層為3 mm 的鋁基板。為提高計算效率，構建的幾何模型為1/4的平板模型，平面尺寸為250 mm×250 mm。

圖4 雷擊損傷電-熱-化學耦合有限元分析1/4涂層平板幾何模型Fig.4 Quarter geometric model for electrical-thermal-chemical coupling analysis of lightning damage of the coating

a）有限元邊界條件與載荷。

在模型的兩個非對稱側面設置零電勢，并在模型上表面施加對流和輻射熱邊界；將雷電弧的電流密度和熱流密度定義為面載荷，施加在模型上表面。雷電流密度分布方程定義如下：

式中J為電流密度；r為徑向坐標；c為常數；Q為熱流密度；R(t)為雷電通道半徑，見式（7）。

式中Ip為峰值電流強度；t為時間。

電流載荷遵循人工模擬雷電試驗SAE APR5412［14］標準，該標準依據自然雷電現象中的電流變化定義了用于試驗研究的電流波形。為使模擬中加載的雷電流波形與試驗電流波形保持一致，提取文獻［1］中雷擊試驗所用雷電流波形曲線上測得的數據點進行擬合，見圖5。

圖5 雷電流波形Fig.5 Lightning current waveform

模型施加電流載荷分兩個階段，第1 階段施加500 μs 的雷電分量A，第2 階段不施加雷電流，散熱至60 s。其中，第1階段電流載荷為圖5b中擬合得到的紅色曲線，其具體函數形式如下：

b）控制方程。

涂層中雷電流傳導過程由如下方程控制：

式中V為控制單元面積；S為控制單元表面；rc為單位體積內體積電流源；m表示S的法向量；σ為電導率矩陣；φ為電勢。

雷電流在涂層中傳導過程中電能會以焦耳熱的形式部分轉換為熱能，相應的溫度場控制方程為

式中ρ為密度；cp為比熱；k為熱導率矩陣；θ為內熱源。

其中，內熱源又分為焦耳熱θe和熱解反應熱Qp，兩者分別由式（11）與式（12）表示。

式中η為能量轉換系數；Hs為聚合物熱解反應的焓。

式（12）的熱解反應熱可由1.2 節中熱解反應分子模擬動力學方程式（4）中定義的熱解速率代入計算得到。

上述涂層樹脂熱解反應在有限元熱-化學反應模塊的定義是通過用戶自定義子程序HETVAL 和USDFLD實現的。

本文利用子程序USDFLD計算并定義樹脂材料的熱解反應和熱解程度，從而控制雷擊過程中材料的電學、熱學參數隨溫度和熱分解反應程度的演變；另一方面利用子程序HETVAL定義了材料的升華吸熱、熱分解吸熱等熱通量。在仿真計算的每一個增量步開始時，Abaqus求解器會將當前的溫度、時間和狀態變量等數據傳遞給子程序USDFLD，隨后USDFLD根據這些數據和定義的算法等計算增量步結束時狀態變量的值，隨后，這些變量的值將被傳遞給Abaqus 求解器和HETVAL，由HETVAL計算當前增量步所產生的熱通量，并將計算結果返回給Abaqus 求解器。Abaqus求解器首先根據USDFLD返回的狀態變量確定當前的熱導率、電導率、比熱容等材料參數，而后再結合HETVAL返回的熱通量計算增量步結束時的溫度。如此不斷迭代，直至仿真結束。

c）材料性能參數。

防熱層、防雷層和鋁基板的材料屬性設置見表2至表4。其中，防雷層的電導率隨溫度變化的規律如圖6所示。

表2 防熱層材料性能參數Tab.2 Material performance parameters of the thermal insulation layer

圖6 防雷層電導率隨溫度變化的曲線Fig.6 Curve of electrical conductivity of the lightning protection layer as a function of temperature

表2中防熱層材料性能參數取自參考文獻［15］。由于防雷層是在與防熱層相同的樹脂基體中摻雜了金屬銅粉制備得到，因此表3中給出的防雷層材料性能是基于表2中硅橡膠的性能參數，根據銅、硅橡膠的比例，按照混合法則計算得到。圖6中防雷層電導率隨溫度變化的數據為試驗測試得到。由于超過300 ℃后，防雷層材料的樹脂基體將開始發生熱分解反應，繼續升溫過程（高于300 ℃）中的電導率測試難以進行，而在實際仿真中，雷擊區域的溫度將達到上千攝氏度，電導率的實測數據顯然不滿足仿真要求。因此，考慮到防雷層中樹脂基體的熱分解會使其中導電金屬粒子接觸電阻降低，本文通過唯象的方法，在防雷材料溫度超過300 ℃后，將其電導率設置為熱解度的線性函數，假設防雷材料完全熱解時的電導率為純銅的電導率（5.714×107S/m）；而在達到銅的沸點（2 562 ℃）時，防雷材料不再導電，電導率退化至1×10-14S/m。

表3 防雷層材料性能參數Tab.3 Material performance parameters of the lightning protective layer

由于雷電流并未擊穿防熱層到達鋁基板，即鋁基板并不會參與導電，所以其溫度變化不明顯，表4中鋁基板的材料參數無需考慮溫度的影響。

表4 鋁基板材料性能參數Tab.4 Material performance parameters of the aluminum substrate

2 結果與分析

2.1 涂層樹脂熱解反應分子模擬的動力學參數計算

涂層樹脂熱解過程分子模擬結果可視化處理如圖7所示，可以發現，在溫度上升的過程中交聯甲基苯基硅橡膠分子不斷解聚，熱分解產生的小分子體現出氣體擴散現象，系統體積也隨之擴大。

圖7 1×1013 K/s升溫速率下硅橡膠熱解反應分子模擬可視化結果Fig.7 Visualization results of molecular simulation of silicone rubber pyrolysis reaction at a heating rate of 1×1013 K/s

基于硅橡膠熱解過程產物數據，根據式（1）計算表1中所列算例的剩余固相產物，進而得到硅橡膠在不同升溫速率下的熱解度隨溫度的變化曲線，部分升溫速率下所得的虛擬熱解度曲線如圖8所示。

圖8 在不同極端升溫速率下硅橡膠熱解反應分子模擬所得虛擬熱解度曲線Fig.8 Virtual pyrolysis degree curves obtained from molecular pyrolysis simulations of silicone rubber at different extreme heating rates

取圖8 中每條熱解度曲線上α為0.1、0.2、…、0.9時的數據，依據式（2）的Starink方程，分別計算ln（β/T1.92）和-1.000 8/RgT，繪制散點圖并擬合，如圖9 所示。圖9 中每組擬合直線的斜率為對應不同熱解度下的活化能Ea，其數值與擬合相關系數Rc如表5所示。

表5 Starink法所得硅橡膠熱解反應分子模擬對應不同熱解度下的活化能Tab.5 Activation energies at different pyrolysis degrees obtained from molecular pyrolysis simulation of silicone rubber by Starink method

圖9 Starink法所得硅橡膠熱解分子模擬對應不同熱解度的活化能擬合結果Fig.9 Fitting results of activation energies for different pyrolysis degrees obtained from molecular pyrolysis simulation of silicone rubber by Starink method

根據式（3）的Malek法，將表5計算得到的活化能Ea分別代入不同升溫速率下的αp，Tp，Ea，Rg，迭代計算得到反應級數n如表6 所示，取算數平均值，得到的反應級數n為1.041。

表6 Malek法所得硅橡膠熱解反應分子模擬的反應級數nTab.6 Reaction order n obtained from the molecular pyrolysis simulation of silicone rubber by Malek method

將活化能Ea、n級化學反應機理函數的積分式G(α)代入式（2）可得到指前因子A。其中，n級化學反應機理函數的積分式G(α)代入反應級數n后形式如下：

將該式與活化能Ea代入式（2），得到指前因子A如表7所示。

表7 Starink法所得硅橡膠熱解反應分子模擬對應不同熱解度的指前因子ATab.7 Pre-exponential factors A at different pyrolysisi degrees obtained from molecular pyrolysis simulations of silicone rubber by Starink method

將活化能Ea、指前因子A、反應級數n代入式（4），就得到了涂層熱解分子模擬的反應動力學方程。

2.2 有限元模擬結果及試驗驗證

根據上述分子模擬計算得到的反應動力學方程編寫用戶自定義子程序HETVAL 和USDFLD，在Abaqus軟件中模擬防雷/防熱涂層的雷擊熱解損傷。

圖10 為雷電流分量A 結束時刻和雷擊結束后自然冷卻60 s時的溫度分布云圖。結果表明，雷擊結束時，雷擊中心處溫度接近2 000 ℃，而散熱60 s 后導電層與防熱層溫度均降至30 ℃左右，遠低于硅橡膠的熱解溫度，因此認為雷擊結束后60 s時防護涂層的熱解損傷過程已經結束。

圖10 標準雷電波形A作用下1/4涂層平板模型溫度分布Fig.10 Temperature distribution of the quarter flat model under the standard lightning waveform A

圖11為標準雷電波形A作用下涂層結構在不同時刻的熱解云圖。結果顯示，標準雷電波形A作用下涂層的熱解損傷形貌呈圓形。表面導電防護層熱解損傷面積在雷擊作用開始的前100 μs 快速擴大，100～500 μs變化速率減慢，而在雷擊作用結束后的散熱過程中幾乎不再擴大。在厚度方向上，熱解損傷在雷擊作用階段及其結束后的前1 ms內主要集中在防雷層，此后才開始向防熱層擴展。熱解損傷深度在雷擊結束后約5 ms時達到最大值。

圖11 標準雷電波形A作用下1/4涂層平板模型在不同時刻下的熱解度云圖Fig.11 Pyrolysis degree contours of quarter plat model at different times under the standard lightning waveform A

這是由于防雷層的熱解損傷主要由雷電流在材料中傳導產生的焦耳熱導致，雷擊過程的前100 μs雷電流較大，在材料中產生了大量的焦耳熱，導致熱解區域快速擴張。在100 μs 后，隨著雷電流強度的減弱，熱解損傷的擴展速率降低。而防熱層不導電，其升溫熱源為防熱層與防雷層之間的熱傳遞，傳遞效率受兩者間溫差及材料導熱效率的影響。而由于防熱層具有較好的隔熱性能，因此需要其與防雷層溫差較大時溫度才會顯著升高并導致材料的熱解損傷，具有一定的滯后性。

為了驗證仿真模型與結果的合理性，對比雷擊熱解損傷有限元模擬結果與硅橡膠防護涂層人工雷擊試驗的損傷形貌，如圖12 所示。由圖12 可以發現，試驗結果與雷擊熱解損傷模擬結果均體現出圓形損傷形貌，但試驗得到的損傷形貌由靠近中心的完全熱解區和外側的放射狀未完全熱解區組成?？紤]到試驗中試樣的接地狀態，雷電流的傳導難以形成未完全熱解區的放射狀損傷形貌，因此，推測放射狀的損傷由雷電電弧掃略并燒蝕形成。而本文所建立的有限元模型僅考慮了雷電流在模型中傳導產生的焦耳熱及其導致的燒蝕損傷，暫未考慮電弧掃略注入的電弧熱等，因此，僅取試驗結果中的完全熱解區與仿真結果進行對比。試驗結果中完全熱解區的直徑約為220 mm，仿真結果得到的熱解區域直徑約為195 mm，兩者相對誤差約為11.4%。以上結果表明，仿真結果與試驗結果較為吻合，所建立的有限元仿真模型較為可靠。

圖12 防雷/防熱涂層雷擊損傷模擬與試驗結果對比Fig.12 Comparison of lightning damage simulation and experimental results for lightning protection/thermal insulation coatings

3 結束語

針對甲基苯基硅橡膠基防雷/防熱涂層，采用反應分子動力學模擬與有限元模擬兩種手段，建立了在雷擊極端溫度環境下涂層熱解損傷的跨尺度模擬方法。模擬所得的涂層完全熱解形貌及面積與人工模擬雷擊試驗結果吻合較好。但由于有限元模型未考慮電弧掃略造成的燒蝕作用，仿真結果未能體現試驗結果中涂層未完全熱解區的放射狀形貌。綜上，文中構建的跨尺度模擬方法能較可靠地預測硅橡膠基防雷/防熱涂層的雷擊熱解損傷，具有較好的工程應用價值。