混合藥型內(nèi)孔燃燒藥柱結構完整性快速分析*

張 路,余 劍,鄧康清,余 瑞,楊 玲,龐愛民

(湖北航天化學技術研究所,湖北襄陽 441003)

0 引言

固體發(fā)動機的固化降溫過程主要依靠發(fā)動機殼體散熱,溫度變化過程中藥柱和殼體會發(fā)生變形,由于它們的熱膨脹系數(shù)不同,導致變形大小也不同,于是在它們內(nèi)部就會產(chǎn)生熱應力和熱應變。當變形過大或者熱應力、應變導致藥柱產(chǎn)生裂紋或脫粘時,會使固體發(fā)動機失效,因此對固化降溫過程中藥柱結構完整性的分析十分重要。而對于混合藥型內(nèi)孔燃燒發(fā)動機,其藥型比較復雜,既有圓管型內(nèi)孔也有星型內(nèi)孔,在進行藥柱結構完整性分析時,由于很難利用試驗去測量實際的結構響應,因此往往需要基于藥柱的本構模型借助于計算機仿真軟件來實現(xiàn)。

固體推進劑是一種非線性粘彈性材料,其力學行為表現(xiàn)出非線性的特性。針對該特性,根據(jù)載荷作用下粘彈性材料的變形大小和應變率的高低,國內(nèi)外已有很多學者[1-6]提出了能解釋粘彈性材料非線性行為的非線性本構模型。但在提出的眾多非線性本構模型中,每一種模型只是針對特定條件下給出,仍然缺乏能夠描述動態(tài)加載條件下固體推進劑大變形的非線性本構模型,且沒有找到一種考慮全部損傷的適合藥柱在復雜工況下力學分析的本構模型。在應用非線性粘彈性本構模型對推進劑藥柱進行完整性分析時,常常要基于有限元軟件二次開發(fā)出相應的子程序。如R.M.Hinterhoelzl[7]、許進升[8]、鄧斌[9-10]等人基于有限元軟件ABAQUS的二次開發(fā)技術,實現(xiàn)了應用非線性粘彈性本構模型對藥柱結構響應的有限元計算。雖然可以借助有限元仿真軟件將含損傷非線性本構模型應用于藥柱結構完整性分析中,但需對現(xiàn)有仿真軟件進行二次開發(fā),技術門檻高,參數(shù)多,成本高,仍需要進一步的研究。

由于線粘彈性模型結構形式簡單,并且已經(jīng)在商業(yè)有限元軟件中進行了開發(fā),具有計算方便,便于操作,在小變形情況下計算精度高的優(yōu)點,現(xiàn)多用于工程應用中。如Shiang-Woei Chyuan[11-12]基于線粘彈本構方程,分析了點火增壓下藥柱的瞬態(tài)應力、應變分布和泊松比對藥柱最大應力、應變的影響;國內(nèi)有學者[13-17]應用線粘彈性模型分析了固化降溫過程中推進劑參數(shù)(如模量、泊松比和熱膨脹系數(shù)等)對藥柱結構完整性的影響。但在實際應用中仍然存在一些缺點,如對固化降溫的不同降溫條件下藥柱的結構完整性進行分析時,要設置不同的計算條件,需要多次計算,使得工作效率低下。

文中針對某混合藥型內(nèi)孔燃燒發(fā)動機藥柱的結構完整性分析,對基于粘彈性分析和波爾茲曼疊加原理提出的快速評估方法進行二次開發(fā),將整個評估過程集成一個操作界面中,可得到固化降溫過程中藥柱危險位置等效應力、等效應變以及藥柱等效模量的實時變化情況,實現(xiàn)藥柱的結構完整性進行分析。最后與基于線粘彈性本構模型的有限元仿真計算結果進行對比,以驗證基于該評估方法開發(fā)的操作系統(tǒng)的有效性和高效性,實現(xiàn)了實際復雜發(fā)動機藥柱結構完整性的快速評估。

1 藥柱完整性快速評估的理論分析

1.1 降溫過程分段

固化降溫過程中藥柱的溫度變化方程為T(t),隨時間變化曲線如圖1所示?，F(xiàn)將溫度作用時間t劃分為n份,每一份時間內(nèi)的溫度都為一個恒定值,其中第i份的作用時間為τi-τi-1,溫度為Ti。這樣,溫度變化方程就可以用一個分段階躍函數(shù)表示,如式(1)所示,曲線也如圖1所示。當時間步長τi-τi-1越小,式(1)溫度階躍函數(shù)就越接近真實的溫度變化曲線,此時由溫度所引起的誤差也就越小;在誤差允許范圍內(nèi),認為兩種曲線重合,該溫度階躍函數(shù)可代替真實的溫度變化函數(shù)。

(1)

溫度變化會導致推進劑藥柱發(fā)生應變,藥柱應變僅與溫差和推進劑的物理特性有關,而與推進劑的彈性模量無關[18]。當推進劑的物理特性不變時,推進劑的應變與溫差成線性關系。令每一個時間步內(nèi)溫度Ti保持不變,該溫度產(chǎn)生的應變?yōu)棣舏,應變隨時間分布如圖2所示。由圖可知,當溫度變化為ΔTi時,產(chǎn)生的相應應變?yōu)棣う舏。當溫度變化時間為t時,藥柱產(chǎn)生的總應變?yōu)椋?/p>

ε(t)=Δε1+Δε2+Δε3+…+Δεn

(2)

另外,由于固體推進劑藥柱是一種粘彈性材料,它的模量與溫度和時間相關。通常情況下,在對推進劑進行完整性分析時,首先要得到推進劑的松弛模量主曲線。推進劑在參考溫度Ts下的松弛模量可表示為：

(3)

參考溫度Ts下的松弛模量如圖3所示。

得到推進劑在參考溫度下的松弛模量主曲線后,可根據(jù)溫度-時間轉換因子αT得到推進劑在各個溫度不同時間的松弛模量,其WLF方程為：

(4)

藥柱的變溫條件近似等效于溫度階躍變化,當溫度階躍變化時,藥柱的模量隨溫度變化曲線如圖4所示。

另外,變溫過程中所用的時間t在參考溫度Ts下的等效時間為ξ,可表示為：

(5)

將等效時間代入式(3),可得到對應的松弛模量為：

(6)

變溫過程中等效時間ξ與變溫時間t的關系如圖5所示,其中ξ=ξn。

1.2 降溫過程藥柱等效模量的理論推導

固化降溫過程中,藥柱整個降溫過程所用時間為t時,將時間t劃分為τ1、τ2、τ3、…、τi-1、τi、…、τn,時間步長ti-ti-1內(nèi)的溫度變化ΔTi所產(chǎn)生的應變Δεi對固體推進劑藥柱的作用時間為t-ti-1,則該段作用時間在參考溫度Ts下的等效時間為ξ-ξi-1。將藥柱在任一時刻總應力、應變的比值稱為等效模量,現(xiàn)對時間t內(nèi)t′(0≤t′≤t)時刻藥柱的等效模量進行分析:

當t′=τ1時,如圖1、圖2所示,藥柱溫度為T1,溫度變化為ΔT1,產(chǎn)生的應變?yōu)棣う?,該Δε1對t′時刻推進劑藥柱的作用時間等效至參考溫度下的等效時間為ξ(T1,τ1),則：

(7)

設應變Δε1在τ1時刻產(chǎn)生的應力為σ1,則

σ1=E(ξ(T1,τ1))·Δε1=E(ξ1)·Δε1

(8)

此時,藥柱的等效模量E(t′)為：

E(t′)=E(τ1)=σ1/Δε1=E(ξ1)

(9)

當t′=τ2時,藥柱溫度為T2,溫度變化為ΔT2,產(chǎn)生的應變?yōu)棣う?,該Δε2對t′時刻推進劑藥柱的作用時間等效至參考溫度下的等效時間為ξ(T2,τ2-τ1),則

(10)

設τ2時刻藥柱的總應力為σ2,則σ2為應變Δε1和Δε2所產(chǎn)生的應力之和,即

σ2=E(ξ(T1,τ2))·Δε1+E(ξ(T2,τ2-τ1))·

Δε2=E(ξ2)·Δε1+E(ξ2-ξ1)·Δε2

(11)

此時,藥柱的等效模量E(t′)為：

E(t′)=E(τ2)=σ2/(Δε1+Δε2)=

(12)

當t′=τn時,藥柱溫度為Tn,溫度變化為ΔTn,產(chǎn)生的應變?yōu)棣う舗,該Δεn對t′時刻推進劑藥柱的作用時間等效至參考溫度下的等效時間為ξ(Tn,τn-τn-1),則

(13)

設τn時刻藥柱的總應力為σn,則σn為應變Δε1、Δε2、Δε3、…、Δεn所產(chǎn)生的應力之和,即

σn=E(ξ(T1,τn))·Δε1+E(ξ(T2,τn-τ1))·Δε2+

E(ξ(T3,τn-τ2))·Δε3+…+E(ξ(Tn,τn-1))·

Δεn=E(ξn)·Δε1+E(ξn-ξ1)·Δε2+E(ξn-ξ2)·

Δε3+…+E(ξn-ξn-1)·Δεn

(14)

此時,藥柱的等效模量E(t′)為：

E(t′)=E(τn)=σn/(Δε1+Δε2+Δε3+…+Δεn)=

(15)

綜上所述,當把變溫時間t分成n個時間步時,每個時間步區(qū)間的時刻點的等效模量都可以求出。時間t內(nèi)任意時刻τn藥柱的等效模量表達式為：

(ξ0=0)

(16)

在計算中為使以上公式處理方便,通常令：

ΔT1=ΔT2=ΔT3=…=ΔTn

(17)

在推進劑的物理性質不變的情況下,可得到[18]：

Δε1=Δε2=Δε3=…=Δεn

(18)

此時,式(19)可化簡為：

(19)

綜上所述,式(19)是變溫條件下推進劑藥柱在不同時刻的模量。

1.3 溫度應變方程

由于在溫度載荷作用下,藥柱的應變僅與溫差和推進劑的物理特性有關,對于已知的推進劑藥柱,其物理特性不變,藥柱的應變與溫差成線性關系,可表示為：

ε(t)=St1·T+St2

(20)

式中,St1、St2為溫度應變系數(shù),它們與發(fā)動機材料的幾何參數(shù)和物理特性有關,與藥柱的模量無關,對于給定的發(fā)動機,St1、St2為常數(shù)(注:不同應變形式的溫度應變系數(shù)不同)。

文中將以上理論分析過程通過計算機語言編程,二次開發(fā)集成了一個藥柱的結構完整性進行快速評估操作界面。

2 藥柱結構完整性快速評估界面開發(fā)

基于C#對以上理論推導公式進行編程,二次開發(fā)建立了一種藥柱結構完整性分析用戶操作界面,操作界面包括固化降溫過程中藥柱等效時間和等效模量的計算、等效應力、應變的求解。操作界面中固體推進劑藥柱結構完整性分析流程如圖6所示。該流程主要包括藥柱材料參數(shù)設置、降溫曲線選擇、溫度應變方程的輸入、求解時間的設定、提交計算以及結果輸出和后處理等步驟。

該界面通過對文中提出的理論分析過程進行集成,可以實現(xiàn)通過輸入相應的降溫方程更加直觀的觀察藥柱等效模量的變化情況,可用于真實發(fā)動機藥柱結構完整性分析。

3 某發(fā)動機藥柱結構完整性分析

3.1 計算模型和材料參數(shù)

XX型發(fā)動機為大長徑比雙燃速發(fā)動機,其部件包括殼體、整體絕熱層、前后人工脫粘層、中燃速藥柱、低燃速藥柱。由于結構的對稱性,對發(fā)動機的1/16進行三維有限元建模,在進行有限元模型建立時,首先對該模型進行簡化,絕熱層的力學性能與人工脫粘層的力學性能相同、中燃速藥柱和低燃速藥柱所用推進劑相同,因此在做前處理時,可將絕熱層和人工脫粘層、中燃速藥柱和低燃速藥柱分別合并成一個整體。網(wǎng)格劃分時,由于該發(fā)動機結構比較復雜,在劃分網(wǎng)格時,為了提高網(wǎng)格質量,減少結構不規(guī)整帶來的困難,將由前向后依次分為3個部分,分別為發(fā)動機頭部、中間段、尾部。XX型發(fā)動機模型共劃分259098個節(jié)點,56761個單元,如圖7所示。在分開的部位和各接觸面分別設置bonded粘接約束。

圓筒發(fā)動機殼體、絕熱層和藥柱的基本材料參數(shù)如表1。

表1 發(fā)動機材料基本性能參數(shù)

固體推進劑藥柱固化溫度為52 ℃,零應力溫度為60 ℃,推進劑在0 ℃下的松弛模量的Prony級數(shù)為：

(21)

其中,推進劑的初始模量E0=21.703 MPa,平衡模量E∞=2.026 33 MPa。

推進劑的相對模量和松弛時間如表2所示。

表2 推進劑相對模量和松弛時間

其中,該推進劑的WLF(Williams-Landel-Ferry)方程的參數(shù)為:參考溫度Ts=273.15 K,C1=13.97,C2=253.7。

3.2 載荷工況及邊界條件

降溫載荷工況:固體發(fā)動機溫度從60 ℃降至-40 ℃,認為降溫過程中藥柱、溫度場在任意時刻都是均勻的,溫度降溫速率為10 ℃/h,降溫所用時間為10 h,降溫方程為T(t)=60-1/360·t。

位移邊界條件:固化降溫過程中,在發(fā)動機殼體兩個端面的一條邊上對發(fā)動機施加位移約束,對稱面滿足對稱邊界條件,不考慮應力過大產(chǎn)生的脫粘。

3.3 藥柱溫度應變方程擬合

由于溫度載荷下,藥柱的應變僅與溫度有關,可表示為式(20)的溫度應變方程,現(xiàn)首先通過有限元軟件仿真計算對結果進行擬合得到藥柱的溫度應變系數(shù)。分析恒定溫度-40 ℃、0 ℃、20 ℃下藥柱的結構響應,基于線彈性模型應用有限元軟件對發(fā)動機進行穩(wěn)態(tài)熱載荷仿真分析,由于溫度載荷作用下藥柱的應變僅與溫度變化有關,分析時取藥柱的彈性模量為2 MPa,分別設置溫度載荷為-40 ℃、0 ℃、20 ℃,得到藥柱內(nèi)部等效應變的分布如圖8、圖9和圖10所示。

由圖8、圖9、圖10可知,溫度載荷作用下藥柱的最大等效應變始終發(fā)生在藥柱圓筒段和星孔段過渡段處,選擇藥柱危險位置A點為研究對象,-40 ℃、0 ℃、20 ℃下A點等效應變及計算時間如表3所示。由表3可知,在三種恒定溫度下A點的等效應變分別為0.197 9、0.122 21、0.082 996,結合初始溫度下藥柱的等效應變?yōu)榱?擬合出該藥柱的溫度應變系數(shù)St1為-0.001 95,St2為0.119 93,因此固化降溫載荷下的藥柱溫度應變方程為：

εeq(t)=-0.001 95·T+0.119 93

(22)

另外,由表3可知,三種恒定溫度下藥柱結構響應分析所用總時間為1 h 54 min。

表3 不同恒定溫度下藥柱內(nèi)表面A點的等效應變

3.4 藥柱結構完整性快速評估

將推進劑材料參數(shù)、降溫曲線T(t)=60-1/360·t、溫度應變系數(shù)St1和St2、計算時間輸入藥柱完整性操作界面中,得到固化降溫過程中等效時間、等效模量、等效應力、等效應變隨降溫時間t分布圖的計算界面如圖11所示。

由圖11可知,藥柱的等效模量在降溫初期短時間內(nèi)會下降至平衡模量,然后隨著降溫時間逐漸增加。固化降溫過程中藥柱的等效應力、應變均隨降溫時間增加而增加,當溫度降到最低時得到藥柱A點的等效應力、等效應變最大,分別為0.645 13 MPa、0.197 93。

3.5 驗證

為了驗證基于該快速評估界面分析的準確性和快速性,基于ANSYS有限元軟件對固化降溫過程中藥柱的力學響應進行粘彈性分析,將得到的結果與評估結果進行對比。計算時選擇軟件內(nèi)自帶的粘彈性本構模型,對發(fā)動機進行熱機耦合瞬態(tài)分析,開啟求解器中的大變形選項,對溫度從零應力溫度降至-40 ℃過程中藥柱內(nèi)部等效應力、等效應變的分布進行研究。圖12為降溫到-40 ℃時藥柱內(nèi)部等效應力、等效應變場。

由圖12可知,基于粘彈性模型有限元仿真得到的藥柱最大等效應力、等效應變分別為0.639 81 MPa,0.199 64。固化降溫過程中藥柱內(nèi)A點等效應力、應變隨時間分布如圖13所示,圖14為A點的等效模量分布。然后將仿真分析得到的結果與快速評估結果進行對比分析,將基于操作界面計算得到的藥柱A點等效應力、應變和等效模量也分別繪制在圖13、圖14中(為觀察方便,只繪制松弛0.000 1 s后的等效模量)。

由圖13、圖14可知,基于該界面評估和基于有限元仿真得到的等效應力、應變曲線基本重合,兩種計算方法得到的等效模量也基本一致。固化降溫條件下藥柱的等效應變和等效應力的最大相對誤差分別發(fā)生在14 040 s和10 800 s時,大小分別為5.497%和5.877%。從而說明基于藥柱結構完整性快速評估界面的計算可以準確預測固化降溫過程中藥柱的等效模量變化及其等效應力、應變的大小,計算精度較好。

為驗證評估方法的快速性,將兩種計算方法所用的時間統(tǒng)計在表4。

表4 不同計算方法所用的時間

由表4可知,應用文中開發(fā)的快速評估界面對混合藥型內(nèi)孔燃燒藥柱結構完整性分析所用時間為1.9 h,相比于常用的基于線粘彈性模型的有限元仿真計算所用的51 h,基于快速評估界面可實現(xiàn)混合藥型內(nèi)孔燃燒藥柱結構完整性的快速分析。

4 總結

1)基于C#語言二次開發(fā)了固化降溫過程中藥柱結構完整性評估界面,該界面包括藥柱等效模量、危險位置等效應力和等效應變的變化情況。

2)應用該評估界面對某型混合藥型內(nèi)孔燃燒發(fā)動機藥柱的結構完整性進行評估,評估結果與基于粘彈性模型的有限元仿真結果吻合很好,等效應力和等效應變的誤差都在6%以內(nèi),但基于評估界面的計算方法所需時間則大大減少。

3)固化降溫過程中,利用該評估界面可計算出固體發(fā)動機藥柱危險位置的等效應力和等效應變的實時變化情況,常用的基于線粘彈性本構模型的有限元仿真方法則得到藥柱內(nèi)部所有點的應力、應變分布,而藥柱的結構完整性破壞往往發(fā)生在藥柱的危險位置,因此,文中提出的評估界面可實現(xiàn)實際復雜藥柱在固化降溫下結構完整性的快速分析。