一種考慮顆粒表面形貌的固體推進劑細觀模型構建方法①

袁 嵩,蔣秋梅,張鎮國*,翁 琳

(1.海軍裝備部,西安 710025;2.西安航天動力技術研究所,固體推進全國重點實驗室,西安 710025;3.上海工程技術大學 城市軌道交通學院,上海 201620)

0 引言

固體推進劑是固體發動機的動力之源,其能量特性直接關系到發動機乃至整彈的性能水平,并且其成型藥柱的結構完整性也直接關系到發動機乃至整彈的工作成敗[1]。當前,火箭總體對固體發動機的精細化設計要求越來越高,固體推進劑藥柱的細觀仿真分析已成為研究熱點之一。固體推進劑通常由金屬燃料Al顆粒、氧化劑AP顆粒、高分子粘合劑基體及防老劑、鍵合劑等組成,是典型的顆粒增強復合材料,構建能夠有效表征推進劑細觀結構特征的代表體積單元(Representative Volume Element,RVE)是進行細觀仿真分析的前提和基礎[2-3]。

目前,細觀結構模型的構建主要基于幾何拓撲方法和CT(Computed Tomography)掃描重構技術。幾何拓撲方法方面,LUBACHEVSKY等[4]最早研究了隨機顆粒堆積的幾何特性及最密堆積方法,對后來建立隨機細觀模型的研究起到奠基性作用。KOCHEVETS等[5-6]進一步發展了LUBACHEVSKY的方法,并最終建立了較為系統的分子動力學算法(MD:Molecular Dynamics)用于顆粒復合材料RVE模型生成。SEGURADO等[7]發展了隨機序列吸附(RSA:Random Sequential Adsorption)方法,并用于顆粒復合材料的細觀界面“脫濕”研究。目前,大多數學者主要采用MD方法或RSA方法建立推進劑的細觀模型,開展推進劑的界面“脫濕”、損傷演化等研究[8-11]。由于三維的RVE模型計算量大且收斂難度高,以往大多數研究均以二維模型為主,采用三維特別是考慮表面特征的非球形三維細觀模型的仿真研究剛剛起步[12-14]。此外,也有學者采用μCT對固體推進劑進行了掃描,重建得到了試樣的三維結構,但是受限于μCT掃描儀器精度,未能分辨出3.5 μm以下的顆粒[15-17]。綜合來看,采用μCT掃描重構技術可以獲取固體推進劑內部的真實結構,但是重構模型的精度受限于μCT掃描設備的精度及重構算法的準確性。目前,采用幾何拓撲方法開展固體推進劑的細觀研究仍是主流,但當前MD算法和RSA算法僅可以構建二維圓形或三維球形RVE模型,且在固體顆粒填充比較高時(一般≥70%),模型生成效率較低,尚缺乏高體積含量非球型RVE模型生成算法方面的研究。

本文基于Voronoi結構,通過Python對Abaqus二次開發,構建一種考慮顆粒復雜形貌的固體推進劑細觀模型構建方法,可快速生成圓形/球形、多邊形/凸多面體二/三維細觀模型,多邊形/凸多面體狀態的填充比理論上可以達到1,豐富推進劑細觀模型構建方法體系。

1 基于Voronoi結構的復雜形貌顆粒RVE模型生成算法

1.1 Voronoi結構介紹

維諾結構(Voronoi structure)又叫狄利克雷鑲嵌(Dirichlet tessellation),廣泛應用在幾何學、地理學、晶體學等學科[18]。典型的Voronoi結構如圖1所示。圖中,紅色點為種子點,藍色線條為Delaunay三角網,紅色線條為Voronoi結構。

圖1 Voronoi結構構造示意圖Fig.1 Schematic diagram of Voronoi structure

目前已有很多開源工具可以很方便地生成二維或三維Voronoi結構,由于Voronoi結構的每個晶胞都具有外凸多邊形的特點,在固體推進劑細觀結構建模中,可以將Voronoi結構離散之后,通過對單個晶胞進行離散、縮放等操作,建立具有復雜形貌的高填充比推進劑RVE模型。

1.2 多邊形/凸多面體顆粒RVE模型生成算法

基于Voronoi結構,實現高填充比的多邊形/凸多面體顆粒RVE模型的具體流程如圖2所示。

圖2 基于Voronoi結構構造多邊形/凸多面體 顆粒RVE模型流程圖Fig.2 Flow chart of construcing polygonal/convex polyhedral particle RVE model based on Voronoi structure

1.2.1 初始顆粒信息庫生成

目前推進劑細觀分析主要考慮AP、Al顆粒。AP顆粒一般有粗粒徑(平均粒徑約150 μm)和細粒徑(平均粒徑約20 μm)兩種規格;Al顆粒平均粒徑約20 μm,每個規格的顆粒實際粒徑均服從正態分布。對于AP顆粒,其分布可以表示為

(1)

對于Al顆粒,其分布為

(2)

則根據推進劑配方規定的具體體積分數,按照式(1)和式(2)可以分別生成每個組分的顆粒數量,每種組分的顆粒數量、組分分布共同組成初始顆粒信息庫。

1.2.2 Voronoi結構生成

目前已有很多開源工具(Neper、Python Scipy等)可以方便地生成服從一定分布的Voronoi結構,本文將初始顆粒信息庫及RVE模型邊界尺寸作為輸入,采用成熟開源軟件Neper生成周期性Voronoi結構,其總體積等于RVE的體積,如圖3所示。Voronoi結構的節點Vi的坐標 (xi,yi,zi)、面Fj幾何結構(組成面的節點,Vj1,Vj2,Vj3,…)、晶胞Ck的幾何結構(組成晶胞的面,Fk1,Fk2,Fk3,…),可以導出作為下一步操作的輸入數據,其中i、j、k分別為節點、面和胞元編號。

圖3 采用Neper生成的典型Voronoi結構Fig.3 Typical Voronoi structure generated by Neper

1.2.3 目標RVE生成

采用Python語言編寫程序調用Abaqus軟件中的幾何模塊功能,根據Neper的幾何輸出文件(.ply文件),按照點-線-面-體之間的拓撲關系,在Abaqus中生成顆粒的幾何模型,基本過程如圖4所示。

圖4 Abaqus中顆粒構造過程Fig.4 Process of constructing particles in Abaqus

采用Python語言編寫導入ply文件的接口程序,主要步驟如下:

(1)讀取ply文件,根據面Fj(j為面編號)包含的點組{Vj1,Vj2,Vj3,…}(其中,j1,j2,j3,…為組成面Fj的點的順序編號)的坐標,在Abaqus中建立幾何點組;用順序相鄰兩點以及首尾兩點連成線,用封閉的線組生成面部件Fj。

(2)根據單個胞元的面的索引,挑出組成胞元的面部件,在Abaqus的“Assembly”模塊里合并生成單個胞元。

(3)重復步驟(2)生成全部胞元。

使用Abaqus提供的內置函數(Abaqus Scripting Reference Guide,7.1.1),可以直接獲取每個晶胞的體積,然后計算各個組分對應的晶胞總體積分數Vf(i),然后根據各組分目標體積分數VfD(i),按照式(3)計算各組分的平均縮放系數Zf(i):

(3)

按照平均縮放系數對各晶胞進行縮放,各組分的實際體積分數已經達到目標體積分數值,但為了增加顆粒間縮放間距的隨機性,對每個組分的每個具體晶胞,采用如式(4)所示的隨機縮放系數進行縮放操作。

(4)

式中Zf(i,j)為第i個組分的第j個具體晶胞;random(0.8,1.2)為0.8～1.2之間的隨機數,隨機數產生的范圍可以根據實際情況進行微調,但需確保每個Zf(i,j)具體值小于1。

采用式(4)進行隨機縮放操作后,再次計算更新各組分對應的晶胞總體積分數Vf(i),然后按照式(3)執行一次縮放操作對體積含量進行修正,即可得到RVE模型中的目標顆粒。

將具有周期性的顆粒進行陣列,并用RVE模型的邊界在顆粒中進行隨機截取,可以得到具有周期性的RVE模型。但是模型邊界處的顆粒會被切開,可能存在很細小的顆粒結構,影響高質量網格的劃分,并可能最終影響模型的收斂性。為了盡可能減少RVE模型邊界處顆粒的細碎程度,編寫自動搜索程序,以模型邊界尺寸的1/n(建議10

圖5 RVE模型邊界顆粒精細化處理示意圖Fig.5 Schematic diagram of fine treatment of boundary particles in RVE model

1.3 多邊形/凸多面體顆粒RVE模型生成算法

基于Voronoi結構也可以生成圓形/球形顆粒的RVE模型,按照圖2所示流程和式(4)生成多邊形/凸多面體顆粒RVE模型后,提取每個晶胞顆粒的質心坐標及其面積(或體積),然后在該質心位置建立具有相同面積的圓或相同體積的球,如圖6所示。Neper提供“sphericity”參數,可以控制多邊形/多面體的球度。當球度較高時,通過多邊形/多面體生成同分布圓形/球形顆粒一般不會出現干涉。但是體積分數較高時,為防止干涉產生,在依次建立顆粒的時候,對新加入顆粒j和已有的顆粒(1～i)進行位置判定,新加入顆粒和前面任意一個顆粒之間的質心間離大于兩圓(或球)半徑之和,如式(5)所示:

(5)

圖6 圓形/球形顆粒RVE模型生成過程示意圖Fig.6 Schematic diagram of generating RVE model of circular/spherical particles

如果式(5)不滿足,則減小新加入的顆粒j的半徑,使任意兩個顆粒的質心間距大于其半徑之和。當所有顆粒均完成干涉檢查及處理后,按照式(3)對總體積分數進行修正,以滿足目標體積分數要求。

2 RVE模型生成及數值仿真實例

2.1 RVE模型生成實例

按照表1所示HTPB推進劑簡化配方,該配方中AP的體積分數為65.5%,其中粗AP顆粒占32.6%,細AP顆粒占32.9%,Al顆粒占10.8%,顆粒總體積分數約為76.3%。

表1 HTPB推進劑典型配方Table 1 Typical formulation of HTPB propellant

參考文獻[19]中的二級配粒徑范圍,本文按照表2所示的粒徑分布參數生成細觀模型。RVE模型的邊長L選為1000 μm。

表2 粒徑分布參數Table 2 Particle size distribution parameter μm

按照本文提出的方法在Abaqus中生成的二維圓形和多面體細觀模型如圖7所示。

(a)Circular particles (b)Polygon particles圖7 二維圓形顆粒、多邊形顆粒RVE模型(Vf =76.3%)Fig.7 2D RVE model composed of circular particles and polygon particles (Vf =76.3%)

由于生成考慮細小顆粒的三維RVE模型需要較高的計算資源,大多數學者在做具體細觀仿真時,首先將小顆粒的Al顆粒、AP顆粒與推進劑基體按照均質化理論進行處理,以均質化之后的基體作為復合基體,顆粒則只考慮大粒徑的AP。按照此方法生成只包含32.6%的二維圓形顆粒RVE模型和二維多邊形顆粒RVE模型如圖8所示,三維圓形顆粒RVE模型和三維多邊形顆粒RVE模型如圖9所示。

(a)Circular particles (b)Polygon particles圖8 二維圓形顆粒、多邊形顆粒RVE模型(Vf=32.6%)Fig.8 2D RVE model composed of circular particles and polygon particles (Vf=32.6%)

(a)Circular particles (b)Polygon particles圖9 三維圓形顆粒、多邊形顆粒RVE模型(Vf=32.6%)Fig.9 3D RVE model composed of circular particles and polygon particles (Vf=32.6%)

2.2 RVE模型數值仿真實例

考慮到實際的計算資源,本文采用2.1節生成的僅包含大粒徑AP的三維細觀模型進行仿真,驗證顆粒形貌對推進劑細觀仿真的影響。

在RVE模型施加準周期性邊界條件:

其中,L=1000 μm,為RVE模型邊長。在y=L面上施加u=500 μm位移載荷,并使用Abaqus中的多點約束方程,使得x=L面上的節點在x方向位移一致,從而x=L面和y=L面在加載過程中保持為平面。

基體采用Neo-Hookean超彈性本構,基體/顆粒界面采用雙線性內聚力模型模擬界面分離,參考文獻[20]中的復合材基體、AP顆粒及界面的材料參數選取,如表3和表4所示。

表3 材料參數選取Table 3 Material parameter selection

表4 界面參數選取Table 4 Interfacial parameter selection

Abaqus仿真結果如圖10所示。可以看到,圓形顆粒RVE模型和多邊形RVE模型均可以模擬出推進劑在受載作用下的界面失效特征,但當AP顆粒表面為多邊形形貌時,可以看到其顆粒附近應力集中現象更為明顯,圓形顆粒附近的最大Mises應力為5.947 MPa,多邊形顆粒附近的Mises應力最大達到7.698 MPa。

(a)Circular particles (b)Polygon particles圖10 ε=30%時Von Mises應力云圖Fig.10 Von Mises stress contours at 30% strain

為了分析兩種顆粒形貌狀態下的界面“脫濕”損傷發展情況,本文提出針對推進劑二維細觀模擬的折算空穴變化率公式,表征界面“脫濕”后形成的空穴,公式如下:

(6)

式中fvoid為模型中cohesive單元的折算空穴變化率;S0為模型初始面積;SDEGi為第i個cohesive單元的折算SDEG損傷值;Si為第i個cohesive單元的面積。

兩種顆粒折算空穴變化率對比如圖11所示。可以看到,在10%應變之前,兩種顆粒的空穴變化率曲線基本重合,表明初始小應變狀態下,界面大部分未失效,空穴變化率維持在3%以下,當應變進一步增大,兩種顆粒形貌的空穴變化率曲線逐漸分離,且多邊形顆粒的空穴變化率略大于圓形顆粒的空穴變化率,說明多邊形顆粒的損傷發展要快于圓形顆粒。

圖11 折算空穴變化率對比Fig.11 Comparison of converted hole area ratio

3 結論

(1) 提出了一種基于Voronoi結構的固體推進劑細觀模型構建方法,可以實現推進劑圓形/球形顆粒RVE模型及考慮顆粒表面復雜形貌特征的多邊形/多面體RVE模型快速生成,多邊形/凸多面體狀態的填充比理論上可以達到1,豐富了推進劑細觀模型構建方法體系。

(2) 數值仿真實例表明本文所建立細觀模型構建方法的有效性,仿真結果顯示當應變大于10%以后,多邊形顆粒的界面損傷發展逐步快于圓形顆粒,建議進一步開展考慮復雜形貌的固體推進劑在寬溫變速率條件下的細觀響應特性,以支撐固體推進劑多尺度研究不斷向精細化發展。