基于單胞體模型的復合固體推進劑松弛模量衰減過程數(shù)值模擬

何 濤

(海軍航空裝備計量監(jiān)修中心，上海200436)

0 引 言

隨著世界軍事科技的發(fā)展，導彈性能的不斷提高，對固體火箭發(fā)動機性能的要求越來越高。如何在保證固體火箭發(fā)動機裝藥結(jié)構(gòu)完整性的前提下，盡可能地提高發(fā)動機的裝藥量，一直是固體火箭發(fā)動機設(shè)計過程中需要解決的關(guān)鍵問題。要解決該問題，首先必須要展開固體推進劑力學性能的研究。復合固體推進劑是一種多相(基體、增強相、界面相等)復合材料，其力學性能受增強相的體積分數(shù)及其組分材料性質(zhì)的影響比較大。因此，基于細觀力學理論研究復合固體推進劑非線性力學性能，不僅可以直觀地反映出組分材料的影響，還可以對其宏觀非線性本構(gòu)關(guān)系的研究提供理論支撐。

復合固體推進劑是一種高填充比的復合材料。近些年來，隨著計算機性能的大幅度提高，計算細觀力學得到了迅速發(fā)展。Matous 等［1-2］在復合固體推進劑顆粒與基體之間的界面層設(shè)置了粘結(jié)單元(Cohesive element)模擬了固體顆粒和基體之間損傷的產(chǎn)生及發(fā)展。Tan 等［3-4］通過數(shù)字圖象等相關(guān)技術(shù)獲得了高能炸藥PBX9501 緊湊拉伸試樣裂尖周圍的應力場及位移場，利用擴展的Mori-Tanaka 方法對試驗結(jié)果做了均勻化處理，并結(jié)合試驗研究得到的顆粒與基體之間的非線性粘結(jié)模型，對含不同尺寸顆粒的細觀模型進行了數(shù)值計算。國內(nèi)一些學者研究了復合固體推進劑的顆粒夾雜模型的建模方法［5-7］。在此基礎(chǔ)上，有的學者采用有限元法對復合固體推進劑進行直接數(shù)值模擬，對推進劑內(nèi)部界面脫粘過程進行了有限元分析［8］。有的學者根據(jù)數(shù)值仿真結(jié)果，結(jié)合細觀力學方法，如:Mori-Tanaka 方法或改進的Mori-Tanaka 方法，研究了固體推進劑的模量、界面脫粘對固體推進劑力學性能的影響等［9-11］，或采用多步法，通過將基體與部分顆粒均質(zhì)化為一種混合物，計算出較復雜的固體推進劑的有效模量［12］。

本文在之前研究［13-14］的基礎(chǔ)上，通過建立復合固體推進劑單胞體模型，根據(jù)復合固體推進劑松弛模量試驗值對基體材料的松弛模量進行參數(shù)反演，結(jié)合數(shù)值計算方法對復合固體推進劑松弛模量的衰減特性進行了研究。

1 物理模型和計算方法

1.1 顆粒夾雜模型及組分材料屬性

復合固體推進劑是一種典型的顆粒增強復合材料，由于顆粒粒徑的分布特征及顆粒的隨機分布，使得其微結(jié)構(gòu)特征非常復雜，如圖1 所示。

圖1 復合固體推進劑顆粒堆積模型及網(wǎng)格

從圖1 中可以看出，復合固體推進劑顆粒堆積模型可以更真實地反映其微結(jié)構(gòu)特征，但這種模型建模過程復雜，并且由于顆粒之間距離較近，使得網(wǎng)格密度較大，計算量較大，這在后續(xù)工作中將逐步展開研究。本文旨在研究復合固體推進劑松弛模量應變的衰減特性，因此首先從單胞體模型出發(fā)進行研究。根據(jù)復合固體推進劑的固體含量確定AP 顆粒的體積分數(shù)，本文所研究的復合固體推進劑AP 顆粒體積分數(shù)為65%，單胞體模型如圖2 所示。

圖2 復合固體推進劑單胞體模型

由之前研究［13］可知，復合固體推進劑的松弛特性主要來源于基體材料，在AP 顆粒的增強作用下，復合固體推進劑的松弛特性與基體材料基本一致，區(qū)別主要體現(xiàn)在瞬時模量E0的增大。因此，基于圖2 所示模型，根據(jù)復合固體推進劑的松弛模量對基體材料的松弛模量進行參數(shù)反演。首先確定基體瞬時模量E0的上、下限值和，并根據(jù)已確定的上、下限值進行有限元計算，結(jié)合細觀力學方法，得出基體瞬時模量E0的上、下限值所對應的復合固體推進劑瞬時模量和;然后對比復合固體推進劑的實際瞬時模量EP0，根據(jù)EP0在和之間的線性分布關(guān)系，在和之間確定下一步計算時新的基體瞬時模量取和的平均值為新的，再次進行有限元計算;反復進行迭代，最后確定基體的松弛模量。假設(shè)AP 顆粒為彈性體，取其彈性模量和泊松比分別為:E=32 450 MPa，ν=0.143 3［1］。

1.2 均勻化方法

當復合固體推進劑單胞體模型受載荷作用時，其體積平均應力和平均應變可由式(1)～(2)計算:

當AP 顆粒體積分數(shù)為65%時，圖2 模型中顆粒之間距離很小，基體材料部分若采用六面體網(wǎng)格，在顆粒之間臨近部分，網(wǎng)格易產(chǎn)生畸變。因此，對于AP 顆?？刹捎昧骟w網(wǎng)格，基體材料部分采用四面體網(wǎng)格。根據(jù)有限元高斯積分原理，分別計算單個六面體單元和四面體單元的平均應力、應變。模型的體積平均應力和應變可以通過單個六面體網(wǎng)格和四面體網(wǎng)格的平均應力和平均應變來確定，即

2 計算結(jié)果及討論

復合固體推進劑松弛模量一般是根據(jù)GJB -770B 在定應變5%下測量的。根據(jù)復合固體推進劑定應變?yōu)?%時的松弛模量，結(jié)合均勻化方法對基體材料的松弛模量進行參數(shù)反演，獲得基體材料的松弛模量。定應變?yōu)?%時單胞體模型拉伸方向上的應力分布如圖3 所示。從圖中可以看出，應力較大的區(qū)域主要集中在顆粒間距小的區(qū)域，隨著時間的增大，應力有明顯的松弛效應。

圖3 定應變?yōu)?%時單胞體模型的應力分布

為研究復合固體推進劑松弛模量的應變相關(guān)性，假設(shè)基體的松弛模量不變，應變小于8%時忽略界面損傷的影響，分別計算當定應變?yōu)?.01%，1%，2%，3%，4%，5%，6%，7%，8%時復合固體推進劑的松弛模量，如圖4 所示。從圖中可以明顯看出，當定應變較小時所計算的松弛模量較大，隨著定應變值的增大，復合固體推進劑松弛模量逐漸減小。所得結(jié)論與?züpek［15］根據(jù)不同定應變(0.5%，1%，2%，5%，10%和15%)下測得的松弛模量試驗結(jié)果得到的結(jié)論一致。

分別對單胞體模型及復合固體推進劑均質(zhì)材料模型的拉伸過程進行數(shù)值計算。其中，復合固體推進劑均質(zhì)材料模型采用不同定應變下的松弛模量。拉伸時間分別取0.05 s 和100 s，最終的真實應變均為7.16%。計算結(jié)果如圖5 和圖6 所示。從圖中可以看出，采用不同定應變水平下的松弛模量的拉伸曲線與單胞體模型的拉伸曲線有明顯差異，只有當真實應變約小于3%時，所選取定應變值越低松弛模量的計算結(jié)果與單胞體模型的計算結(jié)果越接近。

圖4 不同定應變水平下的復合固體推進劑松弛模量曲線

圖5 0.05s 內(nèi)單胞體模型及推進劑拉伸曲線

圖6 100 s 內(nèi)單胞體模型及推進劑拉伸曲線

根據(jù)Boltzmann 疊加原理，多個起因的總效應等于各個起因的效應之和。對于一維模型，當有n個應變增量順次在ζi時刻分別作用于物體，則在ζn以后某時刻t 的總應力可表示為

對于三維模型，當時間增量Δt 足夠小時，Δσij，Δεij應滿足廣義胡克定律，即

因此，可根據(jù)單胞體模型的拉伸曲線，選取一定的時間增量Δt，求解出t 時刻時的松弛模量E(t)，對比初始時刻的松弛模量曲線，從而獲得瞬時模量E0的變化。為驗證該方法及確定合適的時間增量Δt，以復合固體推進劑定應變?yōu)?%時的松弛模量為輸入條件，拉伸時間為0.05 s，取Δt 分別為0.001 s，0.000 5 s 和0.000 1 s 計算復合固體推進劑的松弛模量，并與輸入的松弛模量進行對比，相對誤差如圖7 所示。從圖中可以看出，相對誤差隨時間的增大而增大，這是因為隨著計算時間的增加，在t ～t+Δt 時刻內(nèi)，時間t 之前累計載荷的松弛效應逐漸增大，而在計算的過程中忽略了Δt 時刻內(nèi)松弛效應的影響。雖然相對誤差隨時間的增大而增大，但數(shù)值仍然較小，說明了該計算方法可行。從不同時間增量Δt 的計算結(jié)果來看，當Δt 從0.001 s 縮小10 倍后，相對誤差產(chǎn)生了震蕩，并且減小Δt 時相對誤差降低的幅度不大。輸入條件不變，拉伸時間改為100 s，Δt 取0.001 s，計算的松弛模量與理論值對比如圖8 所示。當拉伸時間為100 s 時，相對誤差為0.83%。采用該方法可以根據(jù)復合固體推進劑的拉伸曲線計算松弛模量的變化。

圖7 拉伸時間為0.05 s 時的松弛模量相對誤差

圖8 拉伸時間為100 s 時的計算值與理論值對比

選取時間增量Δt 為0.001 s，最終真實應變?yōu)?.16%，拉伸時間分別取0.05 s 和100 s 時，復合固體推進劑單胞體模型松弛模量的衰減如圖9 和圖10 所示。從圖中可以看出，在復合固體推進劑單胞體的拉伸過程中，模量有明顯的衰減現(xiàn)象，衰減程度與載荷歷程相關(guān)。

圖9 拉伸時間為0.05 s 時的松弛模量衰減

圖10 拉伸時間為100 s 時的松弛模量衰減

3 總 結(jié)

(1)由于復合固體推進劑材料的顆粒夾雜效應，通過仿真計算可知其松弛模量具有明顯的應變相關(guān)性，即在不同的定應變水平下測得的試驗結(jié)果會存在差異，隨著定應變水平的增大，測得的松弛模量逐漸減小。計算結(jié)果所得結(jié)論與?züpek［15］試驗結(jié)果一致。

(2)根據(jù)單胞體模型不同拉伸速率下，不考慮界面損傷時的拉伸曲線，計算了拉伸過程中復合固體推進劑松弛模量的變化。計算結(jié)果表明，在拉伸過程中，復合固體推進劑的松弛模量有明顯的衰減現(xiàn)象，并且衰減程度與載荷歷程相關(guān)。

［1］Matous K，Inglis H M，Gu X F，et al. Multiscale Damage Modeling of Solid Propellants:Theory and Computational Framework［C］//41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference ＆ Exhibit，AIAA 2005 - 4347，2005:1 -14.

［2］Matous K，Geubelle P H. Multiscale Modeling of Particle Debonding in Reinforced Elastomers Subjected to Finite Deformation［J］. International Journal for Numerical Methods in Engineering，2006，65(2):190 -223.

［3］Tan H，Liu C，Huang Y，et al. The Cohesive Law for the Particle/Matrix Interfaces in High Explosives［J］. Journal of the Mechanics and Physics of Solids，2005，53(8):1892 -1917.

［4］Tan H，Huang Y，Liu C，et al. The Uniaxial Tension of Particulate Composite Materials with Nonlinear Interface Debonding［J］. International Journal of Solids and Structures. 2007，44(6):1809 -1822.

［5］史佩，李高春，李昊. 復合固體推進劑細觀力學模型研究［J］. 計算機仿真，2007，24(5):21 -24.

［6］趙玖玲，田先斌. 顆粒復合材料代表性體元并行建模算法研究［J］. 計算機仿真，2010，27(1):46 -49.

［7］劉著卿，李高春，邢耀國. 復合固體推進劑細觀損傷掃描電鏡實驗及數(shù)值模擬［J］. 推進技術(shù)，2011，32(3):412 -416.

［8］曲凱，張旭東，李高春. 基于內(nèi)聚力界面脫粘的復合固體推進劑力學性能研究［J］. 火炸藥學報，2008，31(6):77 -81.

［9］李高春，邢耀國，王玉峰. 基于細觀力學的復合固體推進劑模量預估方法［J］. 推進技術(shù)，2007，28(4):441 -444.

［10］李高春，邢耀國，戢治洪，等. 復合固體推進劑細觀界面脫粘有限元分析［J］. 復合材料學報，2011，28(3):229 -235.

［11］劉承武，陽建紅，陳飛. 改進的Mori-Tanaka 法在復合推進劑非線界面脫粘中的應用［J］. 固體火箭技術(shù)，2011，34(1):67 -70.

［12］馬昌兵，強洪夫，武文明，等. 顆粒增強復合材料有效彈性模量預測的多步法［J］. 固體力學學報，2010，31(S1):12 -16.

［13］Zhi Shijun，Sun Bing，Zhang Jianwei. Multiscale Modeling of Heterogeneous Propellants from Particle Packing to Grain Failure Using a Surface-Based Cohesive Approach［J］. Acta Mechanica Sinica，2012，28(3):746 -759.

［14］職世君，孫冰，張建偉. 基于表面粘結(jié)損傷的復合固體推進劑細觀損傷數(shù)值模擬［J］. 推進技術(shù)，2013，34(2):273 -279.

［15］Sebnem ?züpek. Constitutive Equations for Solid Propellants［D］. University of Texas at Austin，1997.