邸 克，楊月誠，張昆鵬

(第二炮兵工程大學(xué)601室，西安 710025)

0 引言

固體火箭發(fā)動機推進劑/襯層界面裂紋問題是影響固體發(fā)動機工作可靠性的關(guān)鍵因素。因此，該類問題一直受到學(xué)者們的廣泛關(guān)注。人們主要采用實驗和以有限元為主的數(shù)值方法對其研究。在實驗方面，陽建紅等［1］采用X射線觀測了推進劑的脫粘擴展試件燃燒過程。在有限元法方面，王至存等［2］對推進劑/襯層/絕熱層矩形粘結(jié)試件進行了三維粘彈性應(yīng)力分析，蒙上陽等［3］計算了推進劑/襯層中不同深度脫粘裂紋的應(yīng)力強度因子。但該類界面粘結(jié)區(qū)域的細觀結(jié)構(gòu)或模量對裂紋的影響研究尚未見報道，其根本原因是缺乏一個合適的推進劑/襯層界面力學(xué)模型。

為有效地研究界面力學(xué)問題，人們提出了多種界面模型，如理想界面模型［4］、粘著層界面模型［5］和接觸界面模型［6］等。而對真實情況下的NEPE推進劑/襯層界面，吳豐軍等［7］應(yīng)用納米壓痕儀證實了它是模量和硬度呈現(xiàn)梯度變化的結(jié)構(gòu)，這為功能梯度界面層模型在推進劑/襯層界面裂紋中的應(yīng)用提供了實驗基礎(chǔ)。功能梯度界面層模型［8］是將雙材料界面粘結(jié)區(qū)域視為功能梯度材料(FGM)，并充分考慮了細觀層面的界面粘結(jié)結(jié)構(gòu)對整個模型的影響。早期，該模型一般采用某一種函數(shù)［9］表示界面層材料屬性的變化。而在最近幾年，人們將界面層劃分成兩層［10］或多層［11］。其中，每一層內(nèi)分別采用了不同的函數(shù)描述材料屬性分布。

在多數(shù)情況下，人們對于功能梯度界面層模型的研究采用理論推導(dǎo)方法，通過求解奇異積分方程獲得裂紋的應(yīng)力強度因子。但這種理論研究難以獲取模型應(yīng)力應(yīng)變的分布規(guī)律。另外，在理論研究時，為了能推導(dǎo)出控制方程的解析解，模型中許多參數(shù)的設(shè)定均受到一定限制，這使得該模型在工程應(yīng)用方面還有一些欠缺。為此，本文將功能梯度界面層模型和有限元方法結(jié)合，為固體推進劑/襯層界面建立了一個新的界面層模型，應(yīng)用軟件ABAQUS及其二次開發(fā)程序，獲取了推進劑/襯層界面裂紋尖端附近的應(yīng)力應(yīng)變場，以及Ⅰ和Ⅱ型粘彈性應(yīng)力強度因子隨加載時間的變化規(guī)律。

1 界面層模型和參數(shù)

文獻［7］中說明，襯層和推進劑粘結(jié)在一起時，由于材料的相互擴散和滲透，粘接處形成了模量快速上升區(qū)和模量快速下降區(qū)，這兩個模量變化區(qū)域組成界面層區(qū)。因此，本文中界面層的厚度即可看成上述兩個模量變化區(qū)域的總厚度。對于真實的推進劑/襯層粘接試件，界面層相對于襯層和推進劑厚度很小，而界面層模量分布又較復(fù)雜。為了突出研究重點以及節(jié)約計算量，本文根據(jù)文獻［7］中的數(shù)據(jù)，作適當(dāng)簡化后，建立了推進劑/襯層界面的二維平面應(yīng)力模型。其中，各層均為粘彈性材料(示意圖和有關(guān)參數(shù)的取值分別如圖1和表1所示)。

在界面層左側(cè)的中線處設(shè)置一長度為c的裂紋，界面層兩側(cè)的襯層和推進劑初始彈性模量分別是E1和E3，且它們不隨厚度的變化而變化。模型將界面層劃分成上中下3個區(qū)域，設(shè)界面層中間區(qū)域材料的初始模量為Er(和厚度無關(guān))，通過變化Er體現(xiàn)界面層整體模量的變化。此外，界面層內(nèi)上下兩個區(qū)域中，材料初始彈性模量隨厚度線性變化。計算時，選取界面層上中下3個區(qū)域厚度比為3∶2∶3，并假設(shè)模型各層材料松弛模量均寫成初始模量和時間因子相乘的形式［12］:

為了保證各層材料界面處松弛模量的連續(xù)性，假設(shè)各層材料時間因子f(t)相同，這里將f(t)寫成Prony級數(shù)形式:

表1 模型中各參數(shù)的取值Table 1 Values of the parameters in the model

2 有限元模型和計算工況

對模型進行有限元網(wǎng)格劃分時，各處均采用八節(jié)點四邊形平面應(yīng)力單元(CPS8)，并在裂紋尖端附近加密網(wǎng)格，其整體和局部網(wǎng)格分別如圖2和圖3所示，模型共劃分為9 854個單元。

裂紋尖端處網(wǎng)格劃分成環(huán)狀結(jié)構(gòu)，以便于采用圍線積分技術(shù)計算應(yīng)力強度因子。計算時，各層材料的粘彈性本構(gòu)方程均可寫成以下積分形式:

其中，λ為拉梅常數(shù)，Gv為剪切松弛模量，它和彈性松弛模量Ev的關(guān)系為為泊松比)。對于該模型的平面應(yīng)力問題，部分應(yīng)力應(yīng)變分量均為0，應(yīng)力 σij中只須考慮分量 σxx、σxy和 σyy;應(yīng)變 εij中只需考慮分量 εxx、εyy、εzz以及 γxy。

采用子程序USDFLD定義空間場變量，以實現(xiàn)材料初始模量隨界面層厚度的梯度變化。本文根據(jù)功能梯度材料屬性變化，將y坐標定義成場變量，并設(shè)置材料的初始彈性模量是y坐標的函數(shù)。這時，界面層中各部分初始彈性模量分別表示為

以上3式中，Ea、Eb和Ec分別表示界面層內(nèi)上中下3部分。將它們寫入子程序USDFLD，即可有效地表征功能梯度界面層的模量變化。

考慮法向和剪切兩類加載條件，并在裂紋上下兩端面施加位移載荷。兩類加載條件下，均約束圖1中B和D兩頂點全部自由度，而且位移載荷均在1 s內(nèi)由0線性增加到額定值0.001 mm，隨后保持99 s。

3 計算結(jié)果與分析

根據(jù)文獻［7］的實驗結(jié)果，當(dāng)推進劑/襯層界面貯存時間延長時，界面層彈性模量明顯減小。計算時，令Er取3、18、40 MPa 3個特殊值，以體現(xiàn)不同界面層模量對結(jié)果的影響規(guī)律。這3個值也可分別看成是界面貯存后期、中期和前期的界面層模量。計算結(jié)果中，默認采用Von Mises等效應(yīng)力(MPa)和最大主應(yīng)變分別描述最終時刻的應(yīng)力和應(yīng)變場。

3.1 不同精細度模型結(jié)果比較

為了證明計算的有效性，本文將不同精細度模型的計算結(jié)果進行比較。以圖2的網(wǎng)格結(jié)構(gòu)為基準，重新劃分了疏密兩個網(wǎng)格模型，其單元個數(shù)分別為7 162和13 484。在法向位移加載條件下(Er=18 MPa)，計算了裂紋尖端節(jié)點法向應(yīng)力隨時間的變化，結(jié)果如圖4所示。從圖4可看出，隨著模型精細度的提高，各時刻應(yīng)力逐漸趨向于某一定值。這說明圖2構(gòu)建的模型精細度符合計算要求，網(wǎng)格劃分到此程度已經(jīng)收斂。

3.2 應(yīng)力應(yīng)變場

圖5和圖6分別顯示的是法向和剪切位移加載下的應(yīng)力場?？煽闯?，當(dāng)界面層模量增大時，裂紋尖端附近應(yīng)力整體分布及最大值明顯增大。在界面貯存前期，由于界面層模量明顯高于兩側(cè)推進劑和襯層的模量，在界面層處產(chǎn)生了一個應(yīng)力相對較大的帶狀分布。在界面貯存后期，法向和剪切位移加載下應(yīng)力最大值分別為貯存前期的8.01%和7.84%。因此，可得結(jié)論:隨貯存時間延長，裂紋尖端附近應(yīng)力場明顯減弱。

法向和剪切位移加載下，不同界面層模量條件下的應(yīng)變場如圖7和圖8所示。這時，隨著界面層模量增大，裂紋尖端附近應(yīng)變最大值略有減小。在界面貯存前期，界面層處存在一個應(yīng)變相對較小的帶狀分布，到界面貯存后期，這種帶狀分布基本消失。以上說明，隨著貯存時間的延長，裂紋尖端附近的應(yīng)變最大值略有增大，而界面層處應(yīng)變也逐漸變大。

3.3 應(yīng)力強度因子

法向位移加載下Ⅰ型和剪切位移加載下Ⅱ型應(yīng)力強度因子(KⅠ和KⅡ)隨加載時間的變化規(guī)律分別如圖9和圖10所示。

KⅠ和KⅡ的變化規(guī)律相似，均隨界面層模量的增大而增大，這表明裂紋處材料的相對剛度越大，Ⅰ和Ⅱ型應(yīng)力強度因子越大。對于此變化規(guī)律，本文的數(shù)值計算和相應(yīng)的理論分析［13］是一致的。此外，在位移載荷剛好增加到額定值時(t=1s)，KⅠ和KⅡ均達到最大值，隨后緩慢下降，體現(xiàn)出粘彈性材料的應(yīng)力松弛特征。其中，界面層模量越大，這種下降趨勢越明顯。由上述分析可知，法向位移加載下的KⅠ和剪切位移加載下的KⅡ均隨著界面貯存時間的延長明顯減小。

綜上計算結(jié)果可看出，界面層模量隨界面貯存時間的延長而減小時，裂紋尖端附近的應(yīng)力應(yīng)變場以及裂紋應(yīng)力強度因子均產(chǎn)生一定的變化。因此，在求解推進劑/襯層界面裂紋問題時，應(yīng)充分考慮界面層模量(或界面貯存時間)對其影響，這也同時說明本文模型的建立是十分必要的。

4 結(jié)論

(1)在實驗的基礎(chǔ)上，建立了一個固體推進劑/襯層界面裂紋的界面層有限元模型。相對于其他界面力學(xué)模型，該模型可充分考慮到界面層模量的分布及變化趨勢，進而揭示推進劑/襯層界面貯存時間對斷裂參量的影響規(guī)律。相對于功能梯度界面層的理論模型，該模型可更方便地獲得應(yīng)力和應(yīng)變的分布規(guī)律。

(2)計算結(jié)果表明，隨著界面貯存時間的延長，裂紋尖端附近應(yīng)變的最大值略有增大，而整體應(yīng)力和應(yīng)力強度因子明顯減小。其中，在界面貯存后期，法向和剪切位移加載下應(yīng)力最大值分別減小到貯存前期的8.01%和 7.84%。

(3)本文工作對功能梯度界面層模型和數(shù)值方法的結(jié)合，以及該模型在工程方面的應(yīng)用有一定的借鑒意義。另外，本文也為精確分析推進劑/襯層界面裂紋問題提供了有效的參考。在此后的工作中，該模型可擴展至三維情況，用來分析固體發(fā)動機推進劑/襯層界面裂紋問題。

［1］陽建紅，侯根良，劉朝豐，等.燃燒條件下藥柱脫粘擴展的設(shè)計［J］.推進技術(shù)，2004，25(3):274－277.

［2］王至存，朱祖念，張善祁.固體發(fā)動機藥柱粘結(jié)試件的三維應(yīng)力分析［J］.宇航學(xué)報，1996，17(4):34－38.

［3］蒙上陽，唐國金，雷勇軍.固體發(fā)動機包覆層與推進劑界面脫粘裂紋穩(wěn)定性分析［J］.固體火箭技術(shù)，2004，27(1):46－49.

［4］Williams M L.The stresses around a fault or crack in dissimilar media［J］.Bulletin of the Seismological Society of America.1959，49(2):199－204.

［5］Atkinson C.On stress singularities and interfaces in linear elastic fracture mechanics［J］.International Journal of Fracture，1977，13(6):807－820.

［6］Comninou M.The interface crack［J］.Journal of Applied Mechanics，1977，44(4):631－636.

［7］吳豐軍，彭松，池旭輝.NEPE推進劑/襯層粘接界面層厚度表征方法研究［J］.固體火箭技術(shù)，2008，31(6):650－652.

［8］Delale F，Erdogan F.On the mechanical modeling of the interfacial region in bonded half－planes［J］.Journal of Applied Mechanics，1988，55(2):317－324.

［9］Ozturk M，Erdogan F.Antiplane shear crack problem in bonded materials with a graded interfacial zone［J］.International Journal of Engineering Science，1993，31(12):1641－1657.

［10］Li Y D，Lee K Y.Mechanical modeling and transient antiplane fracture analysis for the graded inter－diffusion regions in a bonded structure［J］.Meccanica，2009，44:389－408.

［11］Huang G Y，Wang Y S，Yu S W.Fracture analysis of a functionally graded interfacial zone under plane deformation［J］.International Journal of Solids and Structures，2004，41:731－743.

［12］Jin Z H.Some notes on the linear viscoelasticity of functionally graded materials［J］.Mathematics and Mechanics of Solids，2006，11(11):216－224.

［13］Cheng Z Q，Gao D Y，Zhong Z.Interface crack of two dissimilar bonded functionally graded strips with arbitrary distributed properties under plane deformations［J］.International Journal of Mechanical Sciences，2012，54:287－293.