肖云東,王玉峰,李何龍,賴帥光,伍 鵬,韓永恒

(1.海軍航空大學(xué), 山東 煙臺 264001; 2. 91526部隊, 廣東 湛江 524000;3. 91468部隊, 海南 陵水 572400;4.海軍裝備部駐北京地區(qū)軍事代表局, 北京 100071)

0 引言

粘接界面的力學(xué)性能是影響固體火箭發(fā)動機壽命的關(guān)鍵因素之一[1]。隨著對粘接界面損傷研究的逐漸深入,眾多學(xué)者從細觀角度出發(fā),運用細觀均勻化理論和有限元法定量分析界面的力學(xué)性能,這有利于掌握界面破壞的本質(zhì)。

有限元法中通常采用內(nèi)聚力模型描述界面的力學(xué)性能,準(zhǔn)確的內(nèi)聚力模型是模擬界面失效過程的前提。關(guān)于如何確定內(nèi)聚力模型的問題中,大多預(yù)先定義模型種類,再通過經(jīng)驗法、反演識別(有限元模型修正技術(shù))等方法確定模型參數(shù)。針對細觀尺度下的界面脫粘問題,由于推進劑內(nèi)部顆粒夾雜,粘接界面處的組成復(fù)雜,采用經(jīng)驗法則很難保證內(nèi)聚力模型的準(zhǔn)確性。

有學(xué)者采用反演識別方法,即基于仿真與實測結(jié)果構(gòu)造目標(biāo)函數(shù),通過最小化目標(biāo)函數(shù),迭代識別細觀尺度下的內(nèi)聚力模型參數(shù)。封濤[2]與職世君[3]建立了復(fù)合固體推進劑二維細觀填充模型,均以仿真與實測的應(yīng)力-應(yīng)變數(shù)據(jù)信息構(gòu)造目標(biāo)函數(shù),分別通過Hooke-Jeeves優(yōu)化算法與分步迭代法對顆粒/基體界面之間的內(nèi)聚力模型參數(shù)進行了反演識別。但其目標(biāo)函數(shù)信息中僅能反映出試件的力學(xué)行為,無法反映出試件的局部信息。

通過數(shù)字圖像相關(guān)技術(shù)(digital image correlation,DIC)可獲取試件全場二維/三維的變形信息,可實現(xiàn)位移和應(yīng)變的量化分析[4]。李高春[5]使用SEM對復(fù)合固體推進劑拉伸變形破壞過程進行了觀察,并通過DIC技術(shù)獲取了推進劑細觀表面位移場?；贒IC技術(shù)的反演識別方法中包含了大量的材料響應(yīng)數(shù)據(jù),在統(tǒng)計學(xué)角度上,該反演識別方法的準(zhǔn)確性更高。Fedele[6]采用DIC技術(shù)與有限元法相結(jié)合的反演識別方法對眩光層壓板的脫膠實驗結(jié)果擬合,確定了指數(shù)型內(nèi)聚力模型參數(shù)。楊思滿[7]通過DIC技術(shù)獲取試件表面的位移場,并與仿真位移場構(gòu)造目標(biāo)函數(shù),采用遺傳算法和序列二次規(guī)劃法相結(jié)合的策略對材料參數(shù)進行反求。

為提高細觀尺度下粘接界面內(nèi)聚力模型的準(zhǔn)確性,通過SEM原位拉伸實驗,記錄粘接試件在拉伸過程中的變形與破壞,根據(jù)界面附近顆粒分布情況建立了簡化后的粘接界面細觀模型,采用DIC技術(shù)結(jié)合Hooke-Jeeves 優(yōu)化算法對界面所采用的雙線性內(nèi)聚力模型相關(guān)參數(shù)開展反演識別研究。

1 基于DIC技術(shù)的反演識別方法

1.1 反演思路

通過相關(guān)文獻[8],預(yù)先定義粘接界面所采用的模型為雙線性內(nèi)聚力模型,經(jīng)過大量的試算工作確定初始模型的相關(guān)參數(shù);采用有限元法和DIC技術(shù)分別獲取ROI的位移場信息;基于仿真與實測位移場信息構(gòu)建目標(biāo)函數(shù),并采用Hooke-Jeeves優(yōu)化算法沿著目標(biāo)函數(shù)減小的趨勢逐步修正待求參數(shù)。反演流程如圖1所示。

圖1 反演流程框圖

1.2 Hooke-Jeeves優(yōu)化算法

Hooke-Jeeves優(yōu)化算法又稱模式搜索算法,通過目標(biāo)函數(shù)來探索設(shè)計空間,可對待求的參數(shù)按照正確的“趨勢”賦值,具有適定性強、收斂速度快等優(yōu)點。算法基于所賦的初值點開始,由“探測移動”沿著不同軸不同方向進行探測性移動,從而確定目標(biāo)函數(shù)值下降的方向,再通過所設(shè)定的“模式移動”,使目標(biāo)函數(shù)值快速沿著“山谷”方向,逐漸逼近目標(biāo)函數(shù)最小值。恰當(dāng)?shù)某跏蓟c可大幅度提高反演精度,降低計算成本。通過“探測移動”與“模式移動”多次交替運行,直至最終符合收斂條件時,計算終止,從而實現(xiàn)尋優(yōu)目的[9]。

2 細觀拉伸實驗

制作了厚度為3 mm的微型非標(biāo)準(zhǔn)試件,如圖2所示。其兩端為有機玻璃,中間部分是通過切取矩形粘接試件得到的推進劑/襯層/絕熱層粘接件,使用302膠將3部分黏結(jié),并自然固化2 d。

圖2 微形非標(biāo)準(zhǔn)試件及其尺寸

通過日本島津JSM-5410LV型試驗系統(tǒng)(全數(shù)字電液伺服加載系統(tǒng)和高精度的掃描電子顯微鏡)進行細觀拉伸實驗,并實時記錄試件在拉伸載荷作用下的力學(xué)行為與表面細觀結(jié)構(gòu)的變化情況。

拉伸過程由伺服控制器自動加載,加載速率設(shè)定為1.2 mm·min-1。SEM放大倍數(shù)為50,調(diào)節(jié)對比度和亮度后,可清晰成像,且圖片特征明顯。拉伸位移與推進劑/襯層/絕熱層粘接件邊長的比值定義為拉伸應(yīng)變ε,襯層附近區(qū)域的破壞過程如圖3所示。

圖3 粘接界面破壞過程(×50)

其中圖3(b)、圖3(c)是ε為8%、18%時,通過開源數(shù)字圖像相關(guān)分析軟件Ncorr獲取圖3(a)虛線內(nèi)沿著x方向的實測位移(Uexpx)云圖。圖3(a)中,2條白色虛線分別位于絕熱層/襯層界面與襯層/推進劑界面處,d0為初始襯層厚度。Δdi表示拉伸過程中的界面張開位移:

Δdi=di-d0

(1)

式中,i=1、2、3分別表示ε為8%、18%、30%的受載狀態(tài)。拉伸過程中的襯層厚度與界面張開位移情況,見表1所示。

表1 襯層厚度與界面張開位移Table 1 Interlayer thickness and opening displacement μm

3 有限元建模

3.1 細觀模型建立

推進劑/襯層/絕熱層粘接試件的代表性體積單元(representative volume element,RVE)中包括HTPB推進劑、HTPB/IPDI襯層、EPDM絕熱層,HTPB推進劑中又包含了大量的固體燃料等顆粒。RVE在宏觀上尺寸需足夠小,可看成一個物質(zhì)點,滿足連續(xù)介質(zhì)力學(xué)的基本假定;在細觀上尺寸需足夠大,包含足夠多的細觀元素和細觀結(jié)構(gòu)信息,可代表局部連續(xù)介質(zhì)的統(tǒng)計平均性質(zhì)。模型中的各部分尺寸,如表2所示。

表2 模型尺寸

其中,基體尺寸滿足推進劑代表性體積單元最小尺寸680 μm×680 μm的要求[10]。Al顆粒粒徑一般為10～30 μm,AP顆粒粒徑總體呈雙峰分布規(guī)律,小顆粒平均粒徑約為20 μm,大顆粒平均粒徑在100～300 μm,各顆粒的粒徑相差大,若將其全部構(gòu)造,會極大程度上增加計算量,且不易收斂。因此,采用多步等效法將粒徑小于80 μm的顆粒等效為基體的一部分,其余顆粒簡化為圓形或橢圓形,通過二維圖像重構(gòu)的方式[11]建立粘接界面細觀模型,如圖4所示。

圖4 粘接界面細觀模型

粘接界面的脫粘行為與顆粒/基體的脫濕行為分別通過黏結(jié)單元、黏結(jié)接觸開展仿真計算。黏結(jié)單元較為常見,不過多介紹。黏結(jié)接觸是通過從屬面與主控面之間的接觸從而定義損傷模型,由從屬面上的點對主控面進行投影,從而確定接觸點對,若對應(yīng)的接觸點不在網(wǎng)格節(jié)點處,可通過鄰近節(jié)點的計算結(jié)果進行插值獲取,具有建模簡單、易收斂等優(yōu)點[12]。

3.2 本構(gòu)關(guān)系

雙線性內(nèi)聚力模型的形式簡單、且易于實現(xiàn),在界面脫粘問題中被廣泛采用[13]。粘接界面脫粘行為與顆粒/基體脫濕行為均采用雙線性內(nèi)聚力模型,如圖5所示。

圖5 雙線性內(nèi)聚力模型

n為垂直界面的法向方向,s為沿界面方向的切向方向,δ0為界面損傷起始位移,δf為界面失效位移。

當(dāng)δ≤δ0時,界面未發(fā)生損傷,界面力T隨界面張開位移的增大而線性增大,其關(guān)系如式(2):

(2)

式中,K為界面剛度。

界面損傷起始準(zhǔn)則采用最大名義應(yīng)力準(zhǔn)則,認為界面力達到界面的最大名義應(yīng)力時,界面開始發(fā)生損傷:

(3)

式中:σmax、τmax為法向與切向的最大名義應(yīng)力;〈〉為Macaulay括號,表示純壓應(yīng)力狀態(tài)不會引發(fā)損傷。

采用基于位移損傷演化準(zhǔn)則,當(dāng)δ0≤δ≤δf時,界面發(fā)生損傷,界面力隨著界面張開位移增加而線性減小:

(4)

式中,D表征界面損傷狀態(tài),對應(yīng)仿真軟件的輸出量為SDEG。當(dāng)D=0時表明界面未損傷,D介于0～1時表明界面處于損傷階段,D=1時表明界面完全失效開裂:

(5)

式中,δm表示混合模式下的界面張開位移:

(6)

可以看出,雙線性內(nèi)聚力模型的形狀可由界面剛度、最大名義應(yīng)力、界面失效位移所決定,為了方便后續(xù)仿真計算,模型簡化為各向同性,將各物理量法向與切向的數(shù)值設(shè)置相同?；谖墨I[2]中的部分界面模型參數(shù),又經(jīng)過大量的有限元試算,最終各界面所采用的模型參數(shù),如表3所示。

表3 各界面的模型參數(shù)

將絕熱層、基體、顆粒簡化為彈性材料,對絕熱層試件進行單向拉伸實驗,取其初始模量作為彈性模量,其泊松比以及基體和顆粒的力學(xué)性能參數(shù)參考文獻[14],如表4所示。

表4 粘接結(jié)構(gòu)各部分力學(xué)性能參數(shù)

3.3 網(wǎng)格劃分與邊界條件

絕熱層與顆粒的楊氏模量較基體相比要大得多,載荷較小的情況下可近似認為二者不發(fā)生變形,對其采用四結(jié)點雙線性平面應(yīng)變四邊形單元CPE4;基體作為近似不可壓材料,受載時體積不發(fā)生變化,采用四節(jié)點四邊形線性積分的雜交單元CPE4H來模擬其在大變形下的力學(xué)響應(yīng);襯層采用四節(jié)點二維粘接單元COH2D4。各界面附近的網(wǎng)格需適當(dāng)細化以提高計算的準(zhǔn)確性和收斂性。

粘接試件在拉伸載荷的作用下,Y方向變形小,對邊界條件進行簡化,仿真過程中不考慮試件沿Y方向的變形。將A、C邊設(shè)置為沿Y方向固定,D邊設(shè)置為完全固定,B邊沿X正方向受均布拉伸載荷,速率為1.2 mm/min。網(wǎng)格劃分與邊界條件的設(shè)置情況,如圖6所示。

圖6 網(wǎng)格劃分與邊界條件

4 參數(shù)反演

4.1 反演過程與結(jié)果

將部分襯層與推進劑區(qū)域作為ROI,為圖3(a)藍色虛線內(nèi)的區(qū)域。文獻[15]中指出界面模型會對試件的力學(xué)行為產(chǎn)生影響,界面本構(gòu)關(guān)系也勢對影響推進劑表面的變形特征,同時,通過細觀拉伸實驗發(fā)現(xiàn):界面的本構(gòu)關(guān)系主要影響ROI沿x方向的位移場分布,對沿y方向的位移場分布的影響可忽略。

以ε為8%、18%時的位移場(Ux)信息共同構(gòu)造目標(biāo)函數(shù),構(gòu)造方法參考文獻[16]并改進:

(7)

式中:p為所有界面參數(shù)組成的向量集合;k為拉伸狀態(tài);i為目標(biāo)點的編號;m為目標(biāo)點的數(shù)量。

反演程序中的增量步長設(shè)為0.05、步長縮減因子設(shè)為0.5、步長加速因子設(shè)為1、增量步長閾值設(shè)為5×10-3,粘接界面的內(nèi)聚力模型參數(shù)初值,見表1所示。

迭代過程中,界面剛度、最大名義應(yīng)力、界面失效位移以及目標(biāo)函數(shù)值的變化情況,如圖7所示。

圖7 迭代過程曲線

迭代過程共計100次賦值計算,目標(biāo)函數(shù)曲線在第83次計算后,開始逐漸平穩(wěn),各曲線的收斂情況較好。在第16輪的“探測移動”中,增量步長為3.12×10-3小于設(shè)定的增量步閾值,完成此輪“探測移動”后,計算終止。第98次賦值試算中的目標(biāo)函數(shù)值最小,認為此界面模型下的ROI的仿真位移與實測位移最為接近,即最終的反演結(jié)果,如表5所示。

表5 反演結(jié)果值

4.2 反演效果分析

反演獲取的模型參數(shù)的準(zhǔn)確性需通過相應(yīng)的實驗數(shù)據(jù)進行驗證。基于反演結(jié)果開展仿真計算,將仿真獲取的ROI位移場與實測數(shù)據(jù)進行比較。ε為8%、18%時,ROI的仿真位移(Usimx)云圖,如圖8所示,相對應(yīng)的實測位移(Uexp x)云圖,見圖3(b)、圖3(c)。可見2種邊界條件下的位移(Ux)云圖分布特征均基本一致,相似度極高。

各節(jié)點處的位移誤差ηU:

(8)

ε為8%、18%時,ROI的位移(Ux)誤差云圖,如圖9所示。

圖9 位移(Ux)誤差云圖

ROI沿x方向的位移誤差統(tǒng)計情況,見表6所示。

主要誤差因素有4個方面:① 界面附近存在大量微裂紋與微孔洞等缺陷,部分顆粒存在破碎與包裹不完全等缺陷,導(dǎo)致ROI的下區(qū)域平均位移(Uexpx)略高于上區(qū)域,而細觀模型中未對其考慮;② 各圖片亮度和對比度在獲取過程中很難保證完全一致,會影響圖片的灰度分布,進而影響Ncorr計算結(jié)果;③ (1～7)顆粒附近區(qū)域的變形梯度大,Ncorr軟件中出現(xiàn)匹配相關(guān)性差、測量精度低的問題;④ 仿真過程中,將各材料簡化為彈性,部分材料的力學(xué)性能參數(shù)直接引用相關(guān)文獻,與實際力學(xué)性能存在偏差。

表6 ROI位移誤差情況(Ux)

5 界面損傷分析

采用界面單元的SDEG對界面的損傷程度進行表征,使用Python腳本語言對仿真軟件進行2次開發(fā),提取內(nèi)聚力界面單元高斯積分點的SDEG數(shù)據(jù),繪制SDEG變化曲線,如圖10所示。

圖10 損傷因子變化曲線

起初在外載荷的作用下,界面力隨著界面張開位移的增加呈線性增長,SDEG均為0,未進入損傷階段;界面力達到最大值后,界面處損傷開始萌生并擴展;隨著界面張開位移的逐漸增加,裂紋逐步形成擴展,界面承受載荷的能力降低;當(dāng)界面力減小至0時,界面損傷達到高峰,SDEG為1,此時界面完全損傷。

結(jié)合粘接界面破壞過程(見圖3),以及獲取的法向內(nèi)聚力模型可知:當(dāng)ε=8%時,由于推進劑的楊氏模量較小,其在載荷作用下的變形明顯,此時界面張開位移為75 μm,SDEG為0,界面力為0.48 MPa;當(dāng)界面張開位移為112.5 μm時,界面力達到最大值0.72 MPa;當(dāng)ε=18%時,(1～7)顆粒脫濕現(xiàn)象明顯,且6號、7號顆粒脫濕情況額外嚴重,脫濕處的承載能力降低,界面附近區(qū)域產(chǎn)生微裂紋,導(dǎo)致(Uexpx)云圖中產(chǎn)生像素點不匹配的無位移輸出區(qū)域,此時界面張開位移為150 μm,大于界面損傷起始位移,界面進入損傷階段,表現(xiàn)出線性軟化行為,SDEG為0.26,界面力為0.70 MPa;當(dāng)ε=30%時,界面與顆粒/基體附近的微裂紋匯合產(chǎn)生宏觀裂紋,且宏觀裂紋的擴展方向幾乎與拉伸方向垂直,此時界面張開位移為375 μm,SDEG為0.75,界面力為0.61 MPa。

6 結(jié)論

1) 根據(jù)雙線性內(nèi)聚力模型中的參數(shù)特性,采用基于DIC技術(shù)與Hooke-Jeeves優(yōu)化算法相結(jié)合的反演識別方法,實現(xiàn)了對細觀尺度下粘接面損傷模型(雙線性內(nèi)聚力模型)中的界面剛度、最大名義應(yīng)力、界面失效位移的求解。當(dāng)加載速率為1.2 mm·min-1時,模型參數(shù)分別為6.40 MPa·mm-1、0.72 MPa、1.80 mm。

2) 基于反演識別方法獲取的界面模型開展仿真計算。ε為8%、18%時,仿真與實測ROI位移(Ux)的平均誤差分別為7.1%、4.9%、最大誤差分別為16.5%、12.8%、誤差小于10%的區(qū)域占比分別為78.6%、85.7%,均表明仿真結(jié)果與實驗吻合較好,所建立的界面模型可以在一定程度上表征界面真實的本構(gòu)關(guān)系。

3) 獲取了界面單元的損傷因子變化曲線,與粘接試件拉伸過程的損傷情況有較好地吻合,驗證了所建立的界面模型可正確地反映出粘接界面在拉伸過程中脫粘規(guī)律。