顏小婷，夏智勛，那旭東，沙本尚

(國防科技大學(xué) 空天科學(xué)學(xué)院, 湖南 長沙 410073)

固體火箭發(fā)動機(jī)具有成本低、結(jié)構(gòu)簡單、可靠性高等特點(diǎn)，在火箭和導(dǎo)彈中均有廣泛應(yīng)用。固體推進(jìn)劑是固體火箭發(fā)動機(jī)的主要能量來源，其燃燒特性可顯著影響固體火箭發(fā)動機(jī)的性能。固體推進(jìn)劑的燃燒特性會同時(shí)受到外部燃燒環(huán)境和內(nèi)部細(xì)觀結(jié)構(gòu)的影響，從細(xì)觀角度研究各種因素對固體推進(jìn)劑燃燒特性的影響，找出各種因素之間的制約關(guān)系，不僅可以掌握改善和提高推進(jìn)劑燃燒性能的方法，還可以對推進(jìn)劑配方的設(shè)計(jì)提供建議，大大降低推進(jìn)劑配方研制的成本和危險(xiǎn)。

當(dāng)對體積足夠大的復(fù)合固體推進(jìn)劑進(jìn)行空間平均分析時(shí)，推進(jìn)劑的化學(xué)性質(zhì)和物理性質(zhì)均趨于一致。然而，從細(xì)觀角度上來講，復(fù)合固體推進(jìn)劑中氧化劑顆粒和金屬顆粒在基體中隨機(jī)分布，會造成局部不均勻。推進(jìn)劑宏觀尺度上的燃燒特性和力學(xué)特性與推進(jìn)劑成分組成、粒徑分布以及細(xì)觀尺度上的顆粒排列均密切相關(guān)，如燃燒速率、金屬團(tuán)聚或彈性模量等均在一定程度上由細(xì)觀結(jié)構(gòu)形態(tài)控制[1]。要分析細(xì)觀結(jié)構(gòu)對推進(jìn)劑性質(zhì)的影響，首先要對推進(jìn)劑的細(xì)觀結(jié)構(gòu)有一個(gè)直觀詳細(xì)的認(rèn)識。目前可采用的試驗(yàn)手段包括光學(xué)分析方法[2]、微型計(jì)算機(jī)電子掃描(microComputer Tomography, micro-CT)[3-5]和X射線斷層攝影術(shù)(X-ray Computed Tomography, XCT)[6-7]等。由于試驗(yàn)方法耗時(shí)長、成本較高，且得到的推進(jìn)劑細(xì)觀結(jié)構(gòu)不易直接用于后續(xù)計(jì)算分析，需要建立一個(gè)能夠描述給定成分組成的推進(jìn)劑異質(zhì)細(xì)觀結(jié)構(gòu)顆粒填充模型，以此來預(yù)測推進(jìn)劑特性。

利用數(shù)值算法生成顆粒填充模型的相關(guān)研究始于20世紀(jì)50年代，其所仿真的對象不僅限于固體推進(jìn)劑，還包括瀝青、混凝土、土壤等任意顆粒填充系統(tǒng)[8]，目前已成功實(shí)現(xiàn)三維填充模型的顆粒形狀包括球體、橢球體、圓柱體以及各種多面體等[9]。表1對目前已有的部分顆粒填充模型進(jìn)行了分類總結(jié)。復(fù)合固體推進(jìn)劑填充體積分?jǐn)?shù)較大且顆粒多采用多級配粒徑分布，因此采用分子動力學(xué)方法更為普遍。Lubachevsky等首次提出了事件驅(qū)動的并行顆粒隨機(jī)裝填算法，也被稱為LS算法[12]。之后又對算法進(jìn)行了改進(jìn)，使其可以用于二維圓盤顆粒和三維球狀顆粒的隨機(jī)裝填[13]。由于LS適應(yīng)性比較強(qiáng)，許多學(xué)者以其為基礎(chǔ)發(fā)展了自己的顆粒填充模型[14-22]。本文所研究的對象為復(fù)合固體推進(jìn)劑，因此采用的是事件驅(qū)動的分子動力學(xué)方法。本文工作是之后固體推進(jìn)劑細(xì)觀燃燒模型和鋁團(tuán)聚模型的基礎(chǔ)。與顆粒粒徑和顆粒質(zhì)量分?jǐn)?shù)相比，顆粒形狀對推進(jìn)劑燃燒性能的影響相對較小，因此模型中所有顆粒簡化為球體。

表1 顆粒填充模型分類

利用顆粒填充模型得到的推進(jìn)劑細(xì)觀結(jié)構(gòu)常用于后續(xù)對推進(jìn)劑的燃燒性能或力學(xué)性能相關(guān)研究。然而，由于結(jié)構(gòu)中顆粒數(shù)目大、計(jì)算量驚人，無法直接利用與實(shí)際試驗(yàn)所用尺寸相同的顆粒堆積體來進(jìn)行數(shù)值計(jì)算。因此，研究者通常利用顆粒堆積體的代表體元(Representative Volume Element, RVE)來進(jìn)行后續(xù)研究。此時(shí)，可以認(rèn)為作為研究對象的復(fù)合材料的微結(jié)構(gòu)為周期性結(jié)構(gòu)，即為代表體元的重復(fù)性排列，代表體元必須可以在統(tǒng)計(jì)意義上代表整個(gè)顆粒堆積體，即與顆粒堆積體具有相同的統(tǒng)計(jì)相關(guān)函數(shù)值。顯然，代表體元的尺寸越大，計(jì)算結(jié)果越準(zhǔn)確，然而所需的計(jì)算成本也越高。因此，研究者必須要在保證計(jì)算精度的前提下找到尺寸最小的代表體元。

材料微觀結(jié)構(gòu)的統(tǒng)計(jì)相關(guān)函數(shù)的量化對于其結(jié)構(gòu)、特性等相關(guān)研究非常重要，不同統(tǒng)計(jì)相關(guān)函數(shù)可以定義特定類別的無序材料的材料特性[23]，還可以用于計(jì)算異質(zhì)推進(jìn)劑的局部密度波動的固有長度，進(jìn)而研究其對燃燒不穩(wěn)定的影響[24]。Smith等對復(fù)合材料的統(tǒng)計(jì)特性進(jìn)行了長期研究。起初的研究對象為由基體和單一半徑顆粒組成的二維可穿透同心殼(Penetrable-Concentric-Shell，PCS)模型[25]，研究了顆粒不同滲透性下兩點(diǎn)函數(shù)與兩點(diǎn)距離和顆粒體積分?jǐn)?shù)之間的關(guān)系。之后對單分散球體進(jìn)行了拓展，針對“最近顆粒表面”建立了“最近表面分布函數(shù)”[26]。Hlushkou等[27]研究了三種不同粒徑分布的顆粒填充結(jié)構(gòu)的統(tǒng)計(jì)特性，計(jì)算得到多分散硬球隨機(jī)填充結(jié)構(gòu)的擴(kuò)散系數(shù)。Lochmann等[28]研究了四種不同顆粒粒徑分布的隨機(jī)填充結(jié)構(gòu)的統(tǒng)計(jì)特性，發(fā)現(xiàn)顆粒粒徑分布對于結(jié)構(gòu)的幾何特性和統(tǒng)計(jì)特性具有較大影響。Kumar等[29]使用Rocpack生成了固體推進(jìn)劑的微觀結(jié)構(gòu)，通過對一點(diǎn)和兩點(diǎn)概率函數(shù)的計(jì)算描述了其統(tǒng)計(jì)形態(tài)，并引入增強(qiáng)模擬退火對周期單元內(nèi)的顆粒位置進(jìn)行了優(yōu)化。Collins等[30]同樣利用Rocpack生成固體推進(jìn)劑微觀結(jié)構(gòu)并對其一點(diǎn)、兩點(diǎn)和三點(diǎn)概率函數(shù)進(jìn)行了研究，研究非線性過程時(shí)，更高階統(tǒng)計(jì)可以對材料給出更精確的描述。

本文采用分子動力學(xué)方法，以硝酸酯增塑聚醚(Nitrate Ester Plasticized Polyether, NEPE)固體推進(jìn)劑作為研究對象，生成了三維球狀顆粒的隨機(jī)填充結(jié)構(gòu)體。利用Monte-Carlo方法對顆粒填充結(jié)構(gòu)體的兩點(diǎn)相關(guān)函數(shù)進(jìn)行了計(jì)算，研究了細(xì)觀特性對顆粒填充結(jié)構(gòu)體的統(tǒng)計(jì)特性的影響。此外，基于得到的兩點(diǎn)相關(guān)函數(shù)計(jì)算了不同顆粒相的最小RVE尺寸，進(jìn)而得到了該NEPE固體推進(jìn)劑的最小RVE尺寸，為后續(xù)NEPE固體推進(jìn)劑細(xì)觀燃燒模型相關(guān)研究提供基礎(chǔ)。

1 復(fù)合固體推進(jìn)劑顆粒填充幾何模型

1.1 算法基本原理

NEPE推進(jìn)劑由基體、氧化劑顆粒和金屬顆粒構(gòu)成，其中氧化劑顆粒和金屬顆粒在基體中隨機(jī)分布。在利用分子動力學(xué)生成推進(jìn)劑顆粒填充模型過程中，首先將氧化劑顆粒和金屬顆粒的中心隨機(jī)分布在基體內(nèi)建立初始構(gòu)型，然后給定顆粒速度和粒徑增長速度，其中粒徑增長速度與顆粒粒徑成正比。采用分子動力學(xué)算法計(jì)算基體內(nèi)顆粒運(yùn)動、顆粒間碰撞以及顆粒變大過程。隨著顆粒粒徑不斷增大，顆粒體積分?jǐn)?shù)不斷上升，最后達(dá)到給定體積分?jǐn)?shù)。顆粒粒徑增大和顆粒運(yùn)動會導(dǎo)致顆粒間碰撞的發(fā)生，從而使得顆粒位置不斷調(diào)整，以避免出現(xiàn)顆粒重疊的情況。當(dāng)顆粒體積分?jǐn)?shù)到達(dá)給定體積分?jǐn)?shù)時(shí)，計(jì)算停止，生成給定配方的推進(jìn)劑顆粒填充模型。由此可見，分子動力學(xué)顆粒填充算法的核心問題是預(yù)測顆粒間的碰撞時(shí)間和碰撞后顆粒運(yùn)動速度的改變量。所有顆粒被簡化為球形顆粒進(jìn)行處理。

由于顆粒在基體中隨機(jī)分布，顆粒的初始位置和初始速度可以由隨機(jī)數(shù)發(fā)生器隨機(jī)給出。在初始時(shí)刻指定顆粒所占基體體積分?jǐn)?shù)為1%，根據(jù)配方顆粒粒徑分布來給定顆粒初始粒徑。指定顆粒i的直徑為di，粒徑增長速度為ai，則t時(shí)刻該顆粒的粒徑為di(t)=ait,(t>di0/ai)，顆粒狀態(tài)變量為S(t)=(Ri(t),Vi(t))，其中Ri(t)和Vi(t)分別為顆粒i的顆粒中心位置和運(yùn)動速度。根據(jù)分子動力學(xué)算法基本原理可以看出，顆粒初始狀態(tài)對最終形成的推進(jìn)劑顆粒填充模型的細(xì)觀結(jié)構(gòu)會有所影響，初始狀態(tài)不同，最終形成的推進(jìn)劑細(xì)觀結(jié)構(gòu)也不同，但是不會影響后續(xù)從宏觀上對推進(jìn)劑屬性的分析和計(jì)算。

在從宏觀角度研究推進(jìn)劑特性時(shí)，可將整個(gè)復(fù)合材料視為代表性體積單元的周期性堆積。由于顆粒的運(yùn)動與初始位置都是隨機(jī)的，可能會出現(xiàn)顆粒越過RVE邊界的情況，此時(shí)需要給出邊界的顆粒施加周期性邊界條件。對于D維問題，認(rèn)為其周邊有3D-1個(gè)計(jì)算域與其相鄰，計(jì)算域邊界設(shè)置為周期性邊界，具有可穿透性，以此來確保填充模型體積分?jǐn)?shù)的一致性和顆粒個(gè)數(shù)的恒定性。

1.2 NEPE推進(jìn)劑顆粒填充幾何模型構(gòu)建

研究對象為NEPE推進(jìn)劑，其主要成分粒徑分布及質(zhì)量分?jǐn)?shù)如表2所示。其中黏合劑包括25%的聚乙二醇(PolyEthylene Glycol, PEG)和75%的塑化劑，塑化劑主要成分組成為50%的硝酸甘油(NitroGlycerin, NG)和50%的1,2,4-丁三醇三硝酸酯(ButaneTriol TriNitrate, BTTN)。

表2 NEPE推進(jìn)劑主要成分粒徑分布及質(zhì)量分?jǐn)?shù)

利用上文介紹的顆粒填充算法，根據(jù)所給配方和給定基體尺寸計(jì)算生成各顆粒的個(gè)數(shù)，計(jì)算得到NEPE推進(jìn)劑的顆粒填充結(jié)構(gòu)體。由于顆粒粒徑之間差距較大，為了提高計(jì)算效率，生成兩種不同邊長的推進(jìn)劑填充體用來分別計(jì)算不同粒徑顆粒相的兩點(diǎn)相關(guān)函數(shù)。實(shí)際生成的NEPE推進(jìn)劑結(jié)構(gòu)體邊長分別為400 μm和1 000 μm，各相的顆粒個(gè)數(shù)和體積分?jǐn)?shù)見表3。

表3 NEPE推進(jìn)劑顆粒填充結(jié)構(gòu)體顆粒組成及體積分?jǐn)?shù)Tab.3 Composites and volume fraction of NEPE propellant pack

由于填充體中顆粒個(gè)數(shù)須為整數(shù)，無法保證不同邊長填充體中各顆粒相體積分?jǐn)?shù)完全相同，在此研究中近似認(rèn)為生成的兩填充體中同種顆粒相的體積分?jǐn)?shù)相同。圖1給出了構(gòu)建的邊長為400 μm的顆粒填充結(jié)構(gòu)體的三維模型。其中，紅色顆粒分別代表AP130和AP13，藍(lán)色顆粒分別代表HMX86和HMX12，白色顆粒代表Al28。

圖1 NEPE推進(jìn)劑顆粒填充結(jié)構(gòu)體Fig.1 Pack of NEPE propellant

圖2 掃描電子顯微鏡下NEPE推進(jìn)劑表面形態(tài)Fig.2 Surface morphology of NEPE propellant observed by the scanning electron microscope (SEM)

圖2給出了在掃描電子顯微鏡下觀察到的NEPE推進(jìn)劑的表面形態(tài)，可以看出，顆粒隨機(jī)分布在基體內(nèi)且整體分布均勻，小顆?？删鶆蚍植加诖箢w粒之間，顆粒所占體積分?jǐn)?shù)較高。這說明利用顆粒填充模型生成的顆粒填充體可以較好地體現(xiàn)出實(shí)際推進(jìn)劑的細(xì)觀結(jié)構(gòu)。為了更準(zhǔn)確直觀地驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確性和可行性，在表4中將利用顆粒填充模型計(jì)算得到的顆粒體積分?jǐn)?shù)與文獻(xiàn)[5]中通過微計(jì)算機(jī)斷層掃描技術(shù)重構(gòu)丁羥推進(jìn)劑細(xì)觀結(jié)構(gòu)得到的顆粒體積分?jǐn)?shù)進(jìn)行對比。對比結(jié)果顯示，本文建立的細(xì)觀填充模型與實(shí)際推進(jìn)劑的顆粒體積分?jǐn)?shù)誤差較小，主要誤差來自于細(xì)觀填充模型無法考慮推進(jìn)劑中可能存在的孔隙。

表4 丁羥推進(jìn)劑二維斷層圖像[5]與細(xì)觀填充體中各組分體積含量對比Tab.4 Comparison of the volume fraction of each component between the two-dimensional tomographic image of butyl hydroxide propellant [5] and the pack

2 復(fù)合固體推進(jìn)劑細(xì)觀結(jié)構(gòu)統(tǒng)計(jì)特性分析

2.1 兩點(diǎn)概率函數(shù)及其Monte-Carlo模擬

給定顆粒條件和顆粒堆積體的邊長時(shí)，顆粒隨機(jī)填充模型可以有多種實(shí)現(xiàn)方式，最終生成N種顆粒堆積構(gòu)型，所有這些不同的顆粒隨機(jī)堆積構(gòu)型即構(gòu)成樣本總體Ω，每個(gè)構(gòu)型均為該樣本總體的成員，且空間為ω。從微觀角度來看，樣本總體中每個(gè)成員的微結(jié)構(gòu)都與其他成員不同，但同時(shí)從宏觀角度考慮則是相同的。宏觀尺度是一種長度尺度，此時(shí)材料可被假設(shè)為具有平均量的連續(xù)結(jié)構(gòu)。因此，可將整體均值定義為集合中所有成員上這些數(shù)量的期待值。即當(dāng)樣本總體中每個(gè)成員的概率和概率密度分別為P(ω)和p(ω)時(shí)，在填充模型中的x點(diǎn)處的函數(shù)F(x,ω)的總體均值可定義為：

(1)

此時(shí)為計(jì)算F(x,ω)的總體均值，需要計(jì)算每個(gè)成員的F(x,ω)值來求取平均值，計(jì)算量較大，因此可通過引入各態(tài)歷經(jīng)性來簡化求值過程。

各態(tài)歷經(jīng)性的基本假設(shè)是樣本總體中任一成員的所有狀態(tài)都可應(yīng)用于樣本總體中的所有成員。由于從宏觀尺度上來講，顆粒填充模型生成的所有構(gòu)型具有相同的物理性質(zhì)和化學(xué)性質(zhì)，因此可以認(rèn)為樣本總體Ω符合各態(tài)歷經(jīng)性假設(shè)，樣本成員的體積為Vω。此時(shí)，可以通過研究單個(gè)樣本成員的體積平均值來得到樣本總體均值，即

(2)

由于在宏觀尺度下材料可被視為具有周期性結(jié)構(gòu)，因此可利用體積為V0的單個(gè)周期單元來計(jì)算樣本總體均值，即

(3)

為描述多相顆粒填充模型的統(tǒng)計(jì)特性，首先要定義特征函數(shù)Iγ(x,ω)。

(4)

即當(dāng)投點(diǎn)x落入γ相顆粒中時(shí)，特征函數(shù)取值為1，否則取值為0。

根據(jù)前文，其總體均值可由式(5)計(jì)算得到。

(5)

n點(diǎn)相關(guān)函數(shù)可以用來研究異質(zhì)顆粒填充模型的微觀結(jié)構(gòu)的統(tǒng)計(jì)特性，其定義式可表示為：

(6)

即n點(diǎn)相關(guān)函數(shù)等于n個(gè)點(diǎn)x1，x2，…，xn隨機(jī)投入空間時(shí)分別同時(shí)處于γ1，γ2，…，γn相的概率。通常來說，n點(diǎn)相關(guān)函數(shù)很難計(jì)算。研究對象僅為顆粒填充模型的幾何參數(shù)即代表體單元的長度，不涉及顆粒間作用力對生成的推進(jìn)劑的力學(xué)性質(zhì)的影響，因此只需研究常用的一點(diǎn)相關(guān)函數(shù)和兩點(diǎn)相關(guān)函數(shù)即可。

一點(diǎn)相關(guān)函數(shù)為最簡單的多點(diǎn)概率統(tǒng)計(jì)，計(jì)算可得到生成推進(jìn)劑顆粒填充模型中x點(diǎn)位于γ相的概率，其定義如下：

(7)

引入各態(tài)歷經(jīng)性，則一點(diǎn)相關(guān)函數(shù)可表示為：

(8)

式中，φγ為顆粒填充模型中γ相所占體積分?jǐn)?shù)。

兩點(diǎn)相關(guān)函數(shù)計(jì)算得到的是推進(jìn)劑模型中兩點(diǎn)x1，x2分別同時(shí)位于γ1相和γ2相的概率，其定義為：

(9)

根據(jù)兩點(diǎn)相關(guān)函數(shù)的定義，在極端條件下兩點(diǎn)相關(guān)函數(shù)滿足以下取值式：

(10)

(11)

其中，δγ1,γ2為克羅內(nèi)克符號(Kronecker symbol)，即

(12)

x12趨于0時(shí)，兩點(diǎn)趨于一點(diǎn)。由于同一點(diǎn)不可能同時(shí)位于兩個(gè)不同相內(nèi)，因此γ1和γ2為不同相時(shí)，兩點(diǎn)相關(guān)函數(shù)取值為零；γ1和γ2為相同相時(shí)，兩點(diǎn)相關(guān)函數(shù)轉(zhuǎn)變?yōu)橐稽c(diǎn)概率函數(shù)問題。x12趨于無窮時(shí)，由于兩點(diǎn)之間距離足夠大，此時(shí)x1和x2具有統(tǒng)計(jì)獨(dú)立性，兩點(diǎn)相關(guān)函數(shù)取值為常數(shù)，與兩點(diǎn)坐標(biāo)及相對位置無關(guān)。

由于一點(diǎn)相關(guān)函數(shù)和兩點(diǎn)相關(guān)函數(shù)較難直接計(jì)算，通常采用Monte-Carlo積分方法對其進(jìn)行估算，即利用無窮和來近似積分。因此一點(diǎn)相關(guān)函數(shù)和兩點(diǎn)相關(guān)函數(shù)可分別由式(13)、式(14)計(jì)算。

(13)

(14)

由上文可知，兩點(diǎn)相關(guān)函數(shù)值會隨兩點(diǎn)之間距離增大最終收斂趨近于φγ1φγ2。本文目的在于得到最小周期單元體的尺寸，當(dāng)隨機(jī)選取兩點(diǎn)x1和x2的兩點(diǎn)相關(guān)函數(shù)值與φγ1φγ2足夠接近時(shí)，可以認(rèn)為尺寸不小于此時(shí)兩點(diǎn)距離的推進(jìn)劑顆粒填充模型與宏觀推進(jìn)劑在幾何統(tǒng)計(jì)角度上具有相同特性。由于算法生成的推進(jìn)劑顆粒填充模型具有各向同性，周期單元體的最小尺寸在各個(gè)方向具有相同值。因此，Monte-Carlo方法的實(shí)施主要通過在顆粒填充模型中多次置入距離為rd的兩點(diǎn)來進(jìn)行，通過計(jì)算不同rd下的兩點(diǎn)概率函數(shù)來得到周期單元體的最小尺寸。

為了兼顧計(jì)算效率和計(jì)算精確性，采用介于投針方法和采樣模板(Sampling Template, ST)方法之間的偽投針方法，即在顆粒填充模型中首先隨機(jī)選取點(diǎn)x1作為球心，然后隨機(jī)給出N2組天頂角θ和方位角φ，以rd作為半徑在球面上隨機(jī)選取N2個(gè)點(diǎn)作為x2，該過程重復(fù)進(jìn)行N1次。對x1的取值施加邊界限制，若x2位于顆粒填充體之外，則對其相應(yīng)坐標(biāo)施加周期性邊界條件以保證計(jì)算結(jié)果的完整性。

2.2 NEPE推進(jìn)劑細(xì)觀結(jié)構(gòu)統(tǒng)計(jì)特性及其RVE最小尺寸

根據(jù)上文可知，在極端條件下兩點(diǎn)相關(guān)函數(shù)取值僅與兩點(diǎn)所在兩相的體積分?jǐn)?shù)相關(guān)，為兩相體積分?jǐn)?shù)之積。因此可以將該值作為收斂值，隨著兩點(diǎn)距離增加，兩點(diǎn)相關(guān)函數(shù)的值會逐漸趨近于該收斂值。當(dāng)兩點(diǎn)相關(guān)函數(shù)值與收斂值的誤差始終小于5%時(shí)，認(rèn)為其達(dá)到收斂值，此后兩點(diǎn)之間距離繼續(xù)增大對于兩點(diǎn)相關(guān)函數(shù)的值不再有影響。達(dá)到收斂值時(shí)的兩點(diǎn)之間距離則為對應(yīng)顆粒相具有統(tǒng)計(jì)均勻性的微觀結(jié)構(gòu)的最小代表體元尺寸。所有顆粒相的最小代表體元尺寸的最大值即為該顆粒填充體的最小代表體元尺寸。

3 兩點(diǎn)相關(guān)函數(shù)計(jì)算結(jié)果及影響因素分析

3.1 取樣次數(shù)和球面取點(diǎn)數(shù)

研究對象為六相體，具有21個(gè)不同的兩點(diǎn)相關(guān)函數(shù)，由于整體來說不同的兩點(diǎn)概率函數(shù)具有相同的分布趨勢，僅數(shù)值不同，因此這里只給出AP13-AP13的兩點(diǎn)概率函數(shù)作比較。

圖3給出了球面取點(diǎn)數(shù)N2均為100時(shí)不同N1下的兩點(diǎn)相關(guān)函數(shù)，宏觀來看曲線基本重合。圖4給出了5次試驗(yàn)下兩點(diǎn)相關(guān)函數(shù)的標(biāo)準(zhǔn)差，可以看出，標(biāo)準(zhǔn)差隨N1的增大逐漸降低，當(dāng)N1取50 000時(shí)，標(biāo)準(zhǔn)差基本比較穩(wěn)定，并始終處于較低水平。

圖3 N2為100時(shí)的兩點(diǎn)相關(guān)函數(shù)Fig.3 Two-point correlation function when N2 is 100

圖4 N2為100時(shí)兩點(diǎn)相關(guān)函數(shù)的標(biāo)準(zhǔn)差Fig.4 Standard deviation of two-point correlation function when N2 is 100

圖5和圖6給出了取樣次數(shù)為50 000時(shí)，不同N2下的兩點(diǎn)相關(guān)函數(shù)及標(biāo)準(zhǔn)差，可發(fā)現(xiàn)此時(shí)N2對取值準(zhǔn)確度的影響已基本可以忽略。N1取50 000時(shí)，無論N2如何取值，兩點(diǎn)相關(guān)函數(shù)達(dá)到穩(wěn)定值以后標(biāo)準(zhǔn)差始終低于0.000 5。因此，為同時(shí)保證計(jì)算精度和計(jì)算效率，在本文中取樣次數(shù)N1和球面取點(diǎn)數(shù)N2分別取50 000和25。

圖5 N1為50 000時(shí)的兩點(diǎn)相關(guān)函數(shù)Fig.5 Two-point correlation function when the N1 is 50 000

圖6 N1為50 000時(shí)的兩點(diǎn)相關(guān)函數(shù)的標(biāo)準(zhǔn)差Fig.6 Standard deviation of two-point correlation function when the N1 is 50 000

3.2 顆粒粒徑

顆粒體積分?jǐn)?shù)相同時(shí)，顆粒粒徑越大，個(gè)數(shù)越少。圖7和圖8分別給出了HMX12和HMX86的兩點(diǎn)相關(guān)函數(shù)和標(biāo)準(zhǔn)差隨兩點(diǎn)距離的變化情況，為提高計(jì)算效率，HMX12的兩點(diǎn)相關(guān)函數(shù)利用邊長為400 μm的填充體進(jìn)行計(jì)算。

圖7 HMX12和HMX86的兩點(diǎn)相關(guān)函數(shù)與兩點(diǎn)距離的關(guān)系Fig.7 Two-point correlation function as a function of rd of HMX12 and HMX86

圖8 HMX12和HMX86的兩點(diǎn)相關(guān)函數(shù)的標(biāo)準(zhǔn)偏差與兩點(diǎn)距離的關(guān)系Fig.8 Standard deviation of two-point correlation functions as a function of rd of HMX12 and HMX86

由圖7～8可以發(fā)現(xiàn)，顆粒體積分?jǐn)?shù)相同時(shí)，兩種顆粒的兩點(diǎn)相關(guān)函數(shù)收斂于相同值，粒徑小的顆粒的兩點(diǎn)相關(guān)函數(shù)達(dá)到收斂時(shí)所需兩點(diǎn)之間的距離更小。由于兩種顆粒個(gè)數(shù)不同，用變異系數(shù)(Coefficient of Variation, CV)來比較數(shù)據(jù)的離散程度。HMX12和HMX86的變異系數(shù)隨兩點(diǎn)距離的變化情況如圖9所示，標(biāo)準(zhǔn)差和變異系數(shù)隨兩點(diǎn)距離與顆粒粒徑比的變化情況分別如圖10、圖11所示。

圖9 HMX12和HMX86的變異系數(shù)與兩點(diǎn)距離的關(guān)系Fig.9 Coefficient of variation comparison of two-point correlation functions as a function of rd of HMX12 and HMX86

圖10 HMX12和HMX86的兩點(diǎn)相關(guān)函數(shù)的標(biāo)準(zhǔn)差與兩點(diǎn)距離與顆粒粒徑比的關(guān)系Fig.10 Standard deviation of two-point correlation functions as a function of rd/d of HMX12 and HMX86

圖11 HMX12和HMX86的變異系數(shù)與兩點(diǎn)距離與顆粒粒徑比的關(guān)系Fig.11 Cofficient of variation comparison of two-point correlation functions as a function of rd/d of HMX12 and HMX86

N1和N2相同時(shí)，從圖9可知HMX12的變異系數(shù)似乎低于HMX86。在圖10和圖11，以兩點(diǎn)距離與顆粒粒徑的比值作為橫坐標(biāo)，可發(fā)現(xiàn)在趨于收斂之后兩種顆粒的標(biāo)準(zhǔn)差和變異系數(shù)均基本處于相同水平。因此可認(rèn)為，兩種顆粒體積分?jǐn)?shù)相同時(shí)，其顆粒粒徑對于N1和N2的選取沒有什么影響。

圖12給出了HMX12和HMX86的兩點(diǎn)相關(guān)函數(shù)相對于收斂值的誤差，可看出HMX12的兩點(diǎn)相關(guān)函數(shù)更快達(dá)到收斂。為更直觀地觀察達(dá)到收斂值時(shí)兩點(diǎn)之間的距離，選取局部圖進(jìn)行放大，如圖13所示。通過觀察和計(jì)算可得到在此推進(jìn)劑填充模型中HMX12的代表體元尺寸為60 μm，HMX86的代表體元尺寸為300 μm。即體積分?jǐn)?shù)相同時(shí)，粒徑小的顆粒所需的RVE尺寸更小。

圖12 HMX12和HMX86的兩點(diǎn)相關(guān)函數(shù)關(guān)于收斂值的誤差Fig.12 Error between two-point correlation function and convergence value of HMX12 and HMX86

圖13 HMX12和HMX86兩點(diǎn)相關(guān)函數(shù)關(guān)于收斂值誤差的局部放大圖Fig.13 Error comparison between two-point correlation functions of HMX12 and HMX86 and the convergence value within a certain range

3.3 顆粒體積分?jǐn)?shù)

圖14比較了AP13和HMX12的兩點(diǎn)相關(guān)函數(shù)，可以發(fā)現(xiàn)在粒徑相近時(shí)，兩種顆粒收斂的速度近似相同，需要分別計(jì)算其最小RVE尺寸才能發(fā)現(xiàn)體積分?jǐn)?shù)對其的影響。圖15給出了AP13和HMX12的變異系數(shù)隨兩點(diǎn)距離的變化情況。從圖15可發(fā)現(xiàn)，AP13的變異系數(shù)基本上均大于HMX12，說明當(dāng)顆粒粒徑相同時(shí)，其在顆粒填充體內(nèi)所占體積分?jǐn)?shù)越大，利用Monte-Carlo方法計(jì)算其兩點(diǎn)相關(guān)函數(shù)時(shí)所需取點(diǎn)越少。

圖14 AP13和HMX12的兩點(diǎn)相關(guān)函數(shù)Fig.14 Two-point correlation functions of AP13 and HMX12

圖15 AP13和HMX12的變異系數(shù)Fig.15 Coefficient of variation comparison of two-point correlation functions of AP13 and HMX12

圖16給出了AP13和HMX12的兩點(diǎn)相關(guān)函數(shù)值對于收斂值的誤差，根據(jù)2.2節(jié)給出的RVE最小尺寸判斷依據(jù)可看出，AP13和HMX12均可選取60 μm作為其代表體元尺寸，由于AP13和HMX12的粒徑較小，這種結(jié)果可能是由取點(diǎn)不夠密集造成。因此，在兩點(diǎn)距離50～70 μm之間重新密集取點(diǎn)計(jì)算兩點(diǎn)相關(guān)函數(shù)及其相對于收斂值的誤差，進(jìn)而計(jì)算AP13和HMX12的代表體元尺寸。圖17為圖16的局部放大圖。根據(jù)RVE最小尺寸判斷依據(jù)，由圖17可得，HMX12在兩點(diǎn)距離為54 μm時(shí)達(dá)到收斂，AP13在兩點(diǎn)距離為60 μm時(shí)達(dá)到收斂，且AP13相對于收斂值的誤差波動更大。說明顆粒粒徑相近時(shí)，體積分?jǐn)?shù)越大，其代表體元長度越小。

圖16 AP13和HMX12的兩點(diǎn)相關(guān)函數(shù)關(guān)于收斂值的誤差Fig.16 Error comparison between two-point correlation functions of AP13 and HMX12 and the convergence value

圖17 AP13和HMX12的兩點(diǎn)相關(guān)函數(shù)關(guān)于收斂值誤差的局部放大圖Fig.17 Error comparison between two-point correlation functions of AP13 and HMX12 and the convergence value within a certain range

3.4 不同顆粒粒徑分布

圖18給出了兩種顆粒相的三個(gè)兩點(diǎn)相關(guān)函數(shù)的分布，圖19比較了這三個(gè)兩點(diǎn)相關(guān)函數(shù)的變異系數(shù)?？梢园l(fā)現(xiàn)，AP130-AP130的變異系數(shù)最大且波動較大，HMX86-HMX86的變異系數(shù)始終處于較低水平且相對平緩，AP130-HMX86的變異系數(shù)基本處于二者之間。

圖18 AP130-AP130，HMX86-HMX86和AP130-HMX86的兩點(diǎn)相關(guān)函數(shù)Fig.18 Two-point correlation functions of AP130-AP130, HMX86-HMX86 and AP130-HMX86

圖19 AP130-AP130，HMX86-HMX86和AP130-HMX86的變異系數(shù)Fig.19 Coefficient of variation comparison of two-point correlation functions of AP130-AP130, AP130-HMX86 and AP130-HMX86

圖20給出了這兩種顆粒相的三個(gè)兩點(diǎn)相關(guān)函數(shù)關(guān)于收斂值的誤差。可發(fā)現(xiàn)，在給定兩點(diǎn)長度內(nèi)，兩點(diǎn)相關(guān)函數(shù)已達(dá)到較好收斂效果。為了更清晰地比較不同兩點(diǎn)相關(guān)函數(shù)收斂速度及其最小代表體元尺寸，圖21通過對數(shù)據(jù)范圍進(jìn)行篩選對圖20進(jìn)行了放大。根據(jù)RVE最小尺寸判斷依據(jù)可得出AP130-AP130的最小代表體元尺寸為750 μm，HMX86-HMX86的最小代表體元尺寸為130 μm，AP130-HMX86的最小代表體元尺寸為120 μm。由此可推，對于多級配模型，不同顆粒相的兩點(diǎn)相關(guān)函數(shù)比單顆粒相的兩點(diǎn)相關(guān)函數(shù)收斂更快，其最小代表體元尺寸比兩種單顆粒相的最小代表體元尺寸小。

圖20 AP130-AP130，HMX86-HMX86和AP130-HMX86的兩點(diǎn)相關(guān)函數(shù)關(guān)于收斂值的誤差Fig.20 Error between two-point correlation functions of AP130-AP130, HMX86-HMX86 and AP130-HMX86 and the convergence value

圖21 AP130-AP130，HMX86-HMX86和AP130-HMX86的兩點(diǎn)相關(guān)函數(shù)關(guān)于收斂值誤差的局部放大圖Fig.21 Error comparison between two-point correlation functions of AP130-AP130, HMX86-HMX86 and AP130-HMX86 and the convergence value within a certain rangee

4 結(jié)論

本文利用基于分子動力學(xué)的復(fù)合固體推進(jìn)劑顆粒填充模型，以NEPE固體推進(jìn)劑為研究對象，生成了指定成分的NEPE推進(jìn)劑的顆粒填充體。為了得到顆粒填充體的最小代表體元尺寸，為之后對于推進(jìn)劑燃燒特性的研究提供模型基礎(chǔ)，基于Monte-Carlo理論，采用偽投針方法對顆粒填充體的統(tǒng)計(jì)特性進(jìn)行了研究。研究發(fā)現(xiàn)：

1)兩種顆粒體積分?jǐn)?shù)相同時(shí)，其顆粒粒徑對于投點(diǎn)次數(shù)的選取幾乎沒有影響，但是對最小代表體元尺寸有較大影響。體積分?jǐn)?shù)相同時(shí)，最小代表體元尺寸與顆粒相直徑具有直接關(guān)系，直徑越大，實(shí)現(xiàn)統(tǒng)計(jì)均勻性所需的兩點(diǎn)距離越大。

2)當(dāng)顆粒粒徑固定時(shí)，其在顆粒填充體內(nèi)所占體積分?jǐn)?shù)越大，利用Monte-Carlo方法計(jì)算其兩點(diǎn)相關(guān)函數(shù)時(shí)所需取點(diǎn)的次數(shù)越少。當(dāng)顆粒相直徑相近時(shí)，顆粒在填充體中體積分?jǐn)?shù)越大，實(shí)現(xiàn)統(tǒng)計(jì)均勻性所需的兩點(diǎn)距離越小，即其最小代表體元尺寸也越小。

3)對于多級配粒徑分布模型，分別位于不同顆粒相的兩點(diǎn)相關(guān)函數(shù)比對應(yīng)的兩種單顆粒相的兩點(diǎn)相關(guān)函數(shù)收斂更快，即多級配顆粒的最小代表體元尺寸要小于對應(yīng)單顆粒的最小代表體元尺寸。為保證計(jì)算精度，顆粒填充體的最小代表體元尺寸應(yīng)不小于填充體內(nèi)所有顆粒相的最小代表體元尺寸中的最大值。因此，在實(shí)際應(yīng)用中，可考慮只計(jì)算單顆粒相的最小代表體元尺寸，其中的最大值即可用作該顆粒填充體的最小代表體元尺寸。