基于粘彈性VCFEM的復(fù)合固體推進(jìn)劑等效松弛模量預(yù)示方法研究

申柳雷，李海陽，申志彬

(國防科技大學(xué)，長沙 410073)

0 引言

復(fù)合固體推進(jìn)劑等效松弛模量是進(jìn)行固體發(fā)動機(jī)藥柱結(jié)構(gòu)完整性分析重要的輸入?yún)?shù)之一，很多學(xué)者針對該問題進(jìn)行了長期研究。彭威等基于Eshelby-Mori-Tanaka等效模型，分別推導(dǎo)了含球形[1]和非球形[2]顆粒的推進(jìn)劑等效松弛模量，認(rèn)為推進(jìn)劑的等效模量可表示為增強(qiáng)系數(shù)與基體相松弛模量的乘積形式，并給出了增強(qiáng)系數(shù)與體積分?jǐn)?shù)、顆粒形狀及基體材料參數(shù)的關(guān)系；李高春等[3]結(jié)合Mori-Tanaka模型和有限元法，計(jì)算了含兩種粒徑顆粒的推進(jìn)劑等效模量；趙玖玲等[4]研究了鍵合劑和AP級配對推進(jìn)劑等效模量的影響規(guī)律；張建偉等[5]采用分子動力學(xué)和有限元法，分析了顆粒體積分?jǐn)?shù)和基體相材料特性與推進(jìn)劑的松弛模量的聯(lián)系。

復(fù)合固體推進(jìn)劑的力學(xué)性能主要由各組分的材料特性、夾雜相的含量與級配等細(xì)觀特征參數(shù)所決定。采用傳統(tǒng)的位移有限元法從細(xì)觀層面預(yù)示推進(jìn)劑力學(xué)性能時(shí)，需要劃分稠密的有限元網(wǎng)格，其帶來的巨大計(jì)算量難以滿足高效率的設(shè)計(jì)需要。美國Ohio州立大學(xué)的Ghosh等[6]提出了Voronoi單元有限元法(Voronoi Cell Finite Element Method,VCFEM)，該方法基于雜交應(yīng)力元的思想，利用Voronoi多邊形進(jìn)行有限單元離散，只需要用少量單元，就可模擬復(fù)合材料中夾雜相形狀和空間分布，已經(jīng)成為一種求解復(fù)合材料等效力學(xué)性能的重要數(shù)值手段。郭然等[7-8]推導(dǎo)了考慮界面脫層的Voronoi單元有限元列式，對顆粒增強(qiáng)復(fù)合材料的損傷過程進(jìn)行了模擬。鄭寧昆[9]將時(shí)域自適應(yīng)精細(xì)算法與VCFEM相結(jié)合，建立了含有夾雜相粘彈性單元的有限元遞推格式。申柳雷[10-12]分別采用彈性VCFEM和不含夾雜相粘彈性VCFEM研究了推進(jìn)劑細(xì)觀特征參數(shù)與等效力學(xué)性能之間的影響關(guān)系。

本文將粘彈性本構(gòu)模型在時(shí)域內(nèi)進(jìn)行離散，推導(dǎo)了含夾雜相粘彈性VCFEM有限元列式的遞推形式。結(jié)合二分法計(jì)算了位移荷載下各時(shí)刻推進(jìn)劑的松弛模量變化，并與ABAQUS對比驗(yàn)證了算法的正確性。運(yùn)用該預(yù)示模型分析了推進(jìn)劑粘彈性力學(xué)性能隨細(xì)觀特征參數(shù)的變化規(guī)律。

1 推進(jìn)劑力學(xué)性能粘彈性VCFEM預(yù)示模型

1.1 Prony級數(shù)粘彈性本構(gòu)模型的離散

如圖1所示，廣義Maxwell模型由多個(gè)Maxwell模型和一個(gè)彈簧原件并聯(lián)組成，假設(shè)廣義Maxwell模型兩端受到應(yīng)力荷載σ的作用，并定義第i個(gè)Maxwell單元的應(yīng)力和應(yīng)變分別為σi和εi，其彈簧原件的彈性系數(shù)為Ei，其粘壺元件的粘性系數(shù)為ηi。由于各Maxwell單元之間相互并聯(lián)，因此每個(gè)單元具有相同的應(yīng)變響應(yīng)，而總應(yīng)力σ等于各單元的應(yīng)力之和。

對于第i個(gè)單個(gè)Maxwell單元，其應(yīng)力應(yīng)變關(guān)系滿足：

(1)

對該微分方程，求解得

σi(t)=Yi(t)εi(t)

(2)

其中，Yi(t)即為松弛模量，寫作

σi(t)=Eie-t/τiεi(t) (i=1,2,…,n)

(3)

其中，τi=ηi/Ei。

由于應(yīng)力值ε相等，所以模型的總應(yīng)變響應(yīng)為

(4)

即

σ(t)=Y(t)ε(t)

(5)

其中，Y(t)即為廣義Maxwell模型的松弛模量。

(6)

采用時(shí)域自適應(yīng)算法[9]，得到第i個(gè)Maxwell單元應(yīng)力應(yīng)變關(guān)系在時(shí)間步ΔTk內(nèi)的m階展開式：

(7)

式中s為時(shí)間步ΔTk內(nèi)一個(gè)無量綱的變量。

對于各向同性材料，將上式拓展到二維情形，得

(8)

其中

將式(8)中的m+1記作m，再根據(jù)s的階次進(jìn)行對應(yīng)，整理可得

(9)

廣義Maxwell模型應(yīng)力應(yīng)變關(guān)系在時(shí)間步ΔTk內(nèi)的m階展開式可寫作

(10)

整理可得

(11)

其中

(12)

(13)

1.2 粘彈性VCFEM有限元列式

對于含有彈性夾雜相顆粒的單個(gè)Voronoi單元(見圖2)，將各物理量在時(shí)域內(nèi)進(jìn)行展開，并根據(jù)s的階次進(jìn)行對應(yīng)，可得到時(shí)間步ΔTk內(nèi)的m階展開遞推格式：

(14)

對于含夾雜相的粘彈性Voronoi單元，其修正余能泛函在時(shí)間步ΔTk內(nèi)的m階的遞推格式為

(15)

應(yīng)力和位移m階的遞推格式寫作

(16)

(17)

圖2 基本Voronoi單元

將以上兩式代入式(15)，整理可得

(βm)TDm+QTqm

(18)

其中

(19)

(20)

Hβm-Gqm+Dm=0

(21)

βm=H-1(Gqm-Dm)

(22)

(23)

將βm向量代回到式(18)，有

(24)

其中

(25)

由余能最小原理和泛函的駐值條件，有

即

Kqm=-Rm+Q

(27)

將式(27)展開，寫作

(28)

其中

(29)

同樣，對每一個(gè)單元進(jìn)行自由度凝聚和剛體位移抑制，得到

(30)

其中

累加每一階的單元邊界節(jié)點(diǎn)位移，即可得到每一個(gè)時(shí)刻的單元邊界節(jié)點(diǎn)位移，進(jìn)一步可得到每一個(gè)時(shí)刻的內(nèi)部節(jié)點(diǎn)位移和各相的應(yīng)力場結(jié)果。

2 粘彈性VCFEM算法

2.1 等效松弛模量計(jì)算方法

每一步時(shí)間步長的根據(jù)上一步的迭代步數(shù)m決定，當(dāng)上一步迭代步數(shù)m小于步數(shù)閾值mlim時(shí)，則將時(shí)間步長擴(kuò)大2倍。為了得到定位移下推進(jìn)劑的等效松弛模量，引入二分法查找符合位移邊界條件的節(jié)點(diǎn)荷載。為此，在每一個(gè)時(shí)間步的內(nèi)循環(huán)基礎(chǔ)上增加了一個(gè)查找的循環(huán)，該循環(huán)通過比較施加的力荷載得到的位移結(jié)果與給定的節(jié)點(diǎn)位移之間的比值，決定下一次節(jié)點(diǎn)力荷載的大小。如果某節(jié)點(diǎn)位移結(jié)果大于給定節(jié)點(diǎn)位移，則取節(jié)點(diǎn)力荷載與給定下限的1/2；反之，如果某節(jié)點(diǎn)位移結(jié)果小于給定節(jié)點(diǎn)位移，則取節(jié)點(diǎn)力荷載與給定上限的1/2。如此反復(fù)迭代，最終可得到節(jié)點(diǎn)位移符合要求的對應(yīng)節(jié)點(diǎn)荷載，進(jìn)而得到應(yīng)力、應(yīng)變場。采用文獻(xiàn)[10]中的直接均勻化方法計(jì)算每一時(shí)刻的等效松弛模量。該計(jì)算方法的具體實(shí)施流程如圖3所示。

圖3 粘彈性VCFEM等效松弛模量計(jì)算流程圖

2.2 粘彈性VCFEM預(yù)示模型正確性檢驗(yàn)

為了驗(yàn)證粘彈性VCFEM算法的正確性，本節(jié)設(shè)計(jì)一個(gè)含30個(gè)顆粒的連續(xù)級配RVE模型作為驗(yàn)證模型。其中，顆粒直徑為100～250 μm的等差數(shù)列，公差為10 μm，共有16級粒徑分布，各級級配均有2個(gè)顆粒。RVE模型的尺寸為1200 μm×1200 μm，夾雜相的面積含量為56.1%，該模型的構(gòu)造和網(wǎng)格劃分均通過第二章所提方法實(shí)現(xiàn)，使用15 μm作為控制間距，顆?？臻g分布見圖4。

在模型頂部施加120 μm的位移荷載，模型下邊界和左邊界布置對稱約束。基體相主曲線參考文獻(xiàn)[5]中的Prony級數(shù)粘彈性模型，其松弛模量為

E(t)=0.244+0.275e-t/1000+0.21e-t/100+

0.23e-t/10+0.41e-tMPa

(31)

基體相泊松比νM=0.495。夾雜相：EI=32.4 GPa，νI=0.14。

如圖5所示，VCFEM共劃分為30個(gè)單元，ABAQUS共劃分了30 879個(gè)四邊形單元，包括14 989個(gè)基體單元和15 890個(gè)夾雜單元。

采用上文所述二分法計(jì)算材料的松弛模量，并與ABAQUS計(jì)算得到的結(jié)果進(jìn)行對比(圖6)，可發(fā)現(xiàn)兩者得到的結(jié)果十分接近，相對誤差在3%以內(nèi)，達(dá)到了較高的精度要求?？紤]到實(shí)際數(shù)值實(shí)現(xiàn)方法上的差異，本文并未對兩者的CPU時(shí)間進(jìn)行比較，但從每次迭代的剛度矩陣規(guī)?？煽闯?，VCFEM具有明顯的計(jì)算效率上的優(yōu)勢。

圖4 多級配RVE模型

(a)VCFEM網(wǎng)格 (b)ABAQUS網(wǎng)格

圖6 VCFEM和ABAQUS松弛模量結(jié)果對比

3 細(xì)觀特征參數(shù)對推進(jìn)劑力學(xué)性能影響分析

為了對推進(jìn)劑粘彈性力學(xué)性能的影響進(jìn)行深化分析，本節(jié)參考工程推進(jìn)劑的常用配方，設(shè)計(jì)了三級配RVE模型作為分析模型(圖7(a))，顆粒級配的粒徑分別為80～120目(120～180 μm)、60～80目(180～250 μm)以及40～60目(250～425 μm)，顆粒之間的控制間距為5 μm。在模型頂部添加10 μm的位移均布荷載，約束模型左側(cè)x方向的位移和模型底部y方向的位移。各級配的面積分?jǐn)?shù)均為20%，故總面積分?jǐn)?shù)為60%，換算為質(zhì)量分?jǐn)?shù)等于76.31%。

夾雜相為AP顆粒，其模量EI=32.4 GPa，泊松比νI=0.14。將HTPB作為基體相材料，泊松比取νM=0.495，為了確定HTPB的松弛模量，根據(jù)《GJB 770B—2005火藥試驗(yàn)方法》等標(biāo)準(zhǔn)[13]，分別進(jìn)行了-50、-30、23、50、70 ℃下的松弛試驗(yàn)，根據(jù)時(shí)溫等效原理對試驗(yàn)結(jié)果進(jìn)行處理，得到如圖7(b)所示的23 ℃下基體的松弛模量主曲線，通過擬合得到Prony級數(shù)形式的松弛模量：

(32)

式中α=233.838 2，而各項(xiàng)系數(shù)如表1所示。

(a)連續(xù)三級配RVE模型

(b)基體相松弛模量23 ℃主曲線

表1 基體相松弛模量Prony級數(shù)系數(shù)

(33)

文獻(xiàn)[1-2,5]均指出，在單軸應(yīng)力作用下，顆粒的增強(qiáng)作用主要體現(xiàn)推進(jìn)劑松弛模量的初始模量變化上，而反映模量隨時(shí)域變化的α值變化較小。因此，本文采用等效初始模量和基體相初始模量的比值(即增強(qiáng)系數(shù)γeff)來表達(dá)顆粒的增強(qiáng)效應(yīng)，相比于等效初始模量，γeff可更直觀地體現(xiàn)出顆粒的增強(qiáng)效果。此時(shí)，推進(jìn)劑的松弛模量可近似表示為

(34)

由于復(fù)合固體推進(jìn)劑是由人工隨機(jī)合成的，AP顆粒在推進(jìn)劑中處于隨機(jī)分布狀態(tài)，而一個(gè)幾何仿真模型只能反映一種分布狀態(tài)，為了從統(tǒng)計(jì)學(xué)角度考慮幾何仿真模型不確定性對結(jié)果的影響，運(yùn)用順序投放算法系統(tǒng)地構(gòu)造了100個(gè)具有相同幾何參數(shù)的RVE模型進(jìn)行分析。

本節(jié)分析結(jié)果都將以含有誤差條的曲線形式呈現(xiàn)，每一個(gè)點(diǎn)的縱坐標(biāo)和誤差條分別為所有樣本結(jié)果的均值和標(biāo)準(zhǔn)偏差，這樣的結(jié)果呈現(xiàn)方式可更好地表征結(jié)果的離散程度和不確定性。

3.1 基體相材料參數(shù)的影響分析

工程實(shí)踐中，在不影響燃燒性能和能量特性的情況下，可通過添加不同劑量和不同種類的小組分添加劑，改變基體相的初始模量，這為推進(jìn)劑的設(shè)計(jì)提供了調(diào)節(jié)參數(shù)；另外，作為粘彈性材料，由于溫度效應(yīng)，基體相材料在不同溫度下面表現(xiàn)出不同的松弛模量。因此，有必要考察基體模量的改變對推進(jìn)劑等效松弛模量的影響。本節(jié)固定其他細(xì)觀參數(shù)，通過調(diào)整基體相初始模量EM0和泊松比νM，以分析基體相材料參數(shù)的影響，得到圖8所示的曲線。

為了表達(dá)清晰，圖8(a)選擇3條松弛模量曲線進(jìn)行分析。可看出，隨著EM0的減小，推進(jìn)劑的松弛模量隨著時(shí)間依然保持S形的趨勢，即由等效初始模量逐漸向等效平衡模量過渡，當(dāng)EM0取值較小時(shí)(EM0=0.5 MPa)，初始模量和平衡模量非常接近，在不同EM0取值的曲線中，該曲線的變化趨勢不是很明顯。由圖8(b)可知，隨著EM0增大，推進(jìn)劑的等效初始模量也呈線性增長，說明兩者呈線性正相關(guān)關(guān)系。而不同基體相泊松比νM(分別取0.01、0.3和0.495)取值對推進(jìn)劑等效模量的影響不大。

(a)不同基體相模量的松弛模量曲線

(b)等效初始模量

3.2 夾雜相材料參數(shù)的影響分析

本節(jié)計(jì)算不同夾雜相材料參數(shù)取值時(shí)的推進(jìn)劑的松弛模量，得到圖9所示的曲線。觀察圖9(a)可知，即使EI取值很小(10 MPa)，也可觀察到推進(jìn)劑的松弛過程，而隨著EI增大，等效松弛模量基本相等；結(jié)合圖9(b)可看出，當(dāng)夾雜相模量大于2000 MPa時(shí)，等效初始模量的結(jié)果基本相等，驗(yàn)證了夾雜相模量變化可以忽略不計(jì)的結(jié)論。同時(shí)，由圖9(b)還可發(fā)現(xiàn)，夾雜相泊松比νI的變化對推進(jìn)劑的等效初始模量也沒有明顯的影響。在實(shí)際工程中，夾雜相模量EI遠(yuǎn)大于基體相初始模量EM0。所以，EI的變化不會對推進(jìn)劑力學(xué)性能產(chǎn)生明顯的影響。

3.3 細(xì)觀結(jié)構(gòu)參數(shù)的影響分析

夾雜相含量是推進(jìn)劑主要細(xì)觀特征參數(shù)，通過同步改變模型中各級配顆粒的數(shù)量，得到不同質(zhì)量分?jǐn)?shù)下的RVE模型。由圖10(a)可知，隨著質(zhì)量分?jǐn)?shù)的增加，推進(jìn)劑的增強(qiáng)系數(shù)γeff呈現(xiàn)明顯的增大趨勢，且增長速度明顯加快。由于隨著質(zhì)量分?jǐn)?shù)的增加，推進(jìn)劑等效模量加速增長，為了作圖清晰，圖10(a)的誤差條的取值等于標(biāo)準(zhǔn)差與等效模量的比值?？煽闯觯S著質(zhì)量分?jǐn)?shù)增大，增強(qiáng)系數(shù)的結(jié)果越離散。

(a)等效松弛模量隨時(shí)間變化曲線

(b)等效初始模量

推進(jìn)劑不僅是結(jié)構(gòu)材料，更是一種燃燒劑，因而在設(shè)計(jì)中，顆粒含量常常根據(jù)能量特性的需要提前確定。因此，很多文獻(xiàn)將研究重點(diǎn)轉(zhuǎn)移到各級配的質(zhì)量分布上。保持ψ1+ψ2=50.87%和ψ3=25.44%不變，增加第一級級配的質(zhì)量分?jǐn)?shù)ψ1，并相應(yīng)減少第二級級配的質(zhì)量分?jǐn)?shù)ψ2，實(shí)現(xiàn)夾雜相級配質(zhì)量分布的改變。如圖10(b)所示，當(dāng)ψ1取不同值時(shí)，增強(qiáng)系數(shù)γeff的變化幅度較小(<10%)。

保持顆粒質(zhì)量分?jǐn)?shù)76.31%不變時(shí)，改變第一級級配與第二級級配的粒徑比例R1/R2，得到不同粒徑比的RVE模型。通過計(jì)算，得到圖10(c)所示的變化曲線。觀察可知，粒徑比變化對增強(qiáng)系數(shù) 的影響并不是很明顯。

(a)質(zhì)量分?jǐn)?shù)

(b)夾雜相粒徑質(zhì)量分?jǐn)?shù)

(c)夾雜相粒徑比

4 結(jié)論

(1)隨著基體相初始模量EM0的增大，等效初始模量Eeff呈線性增大；

(2)在給定的夾雜相模量EI變化范圍內(nèi)，EI和泊松比νI的變化對Eeff的影響可以忽略；

(3)隨著夾雜相質(zhì)量分?jǐn)?shù)ΨI的增加，模量增強(qiáng)系數(shù) 呈現(xiàn)加速遞增的趨勢，在保持總質(zhì)量不變的情況下，改變級配的粒徑比和質(zhì)量分布對γeff的影響較小。

通過本文預(yù)示工作，有助于篩選出可供優(yōu)化的細(xì)觀設(shè)計(jì)參數(shù)，為推進(jìn)劑細(xì)觀結(jié)構(gòu)優(yōu)化設(shè)計(jì)提供指導(dǎo)。