基于DIC 技術(shù)的固體火箭發(fā)動機裝藥粘接界面參數(shù)反演研究

肖云東，王玉峰，李高春，孔令澤，伍 鵬，賴帥光

（1. 海軍航空大學(xué)，山東 煙臺 264001；2. 91468 部隊，海南 陵水 572400；3. 91526 部隊，廣東 湛江 524000）

0 引言

固體火箭發(fā)動機作為各類飛行器的動力裝置，具有貯存時間長、使用方便、結(jié)構(gòu)簡單等優(yōu)點，在航空航天等領(lǐng)域被廣泛使用［1］。但固體火箭發(fā)動機的推進劑/襯層/絕熱層界面易脫粘［2］，所產(chǎn)生的故障約占發(fā)動機總故障的1/3［3］，這是影響固體火箭發(fā)動機能否正常工作的關(guān)鍵問題，也是對其開展有限元仿真的難點。

大多學(xué)者采用內(nèi)聚力模型對粘接界面的脫粘行為［4-6］進行有限元仿真，準確的內(nèi)聚力模型是獲得有效仿真結(jié)果的重要前提。當(dāng)前，關(guān)于如何確定內(nèi)聚力模型的相關(guān)研究較少，大多預(yù)先定義模型種類，模型參數(shù)再通過經(jīng)驗法、實驗法、反演識別（有限元模型修正技術(shù)）等方法確定。常見的實驗法有：J 積分法［7］、“含能力”的平衡法［8］、直接拉伸法［9］。前2 種方法僅適用于線彈性材料；后一種方法則需試樣在各處損傷一致，這在實驗中幾乎不可能實現(xiàn)。

學(xué)者們基于仿真與實測數(shù)據(jù)信息，采用反演識別方法獲取內(nèi)聚力模型相關(guān)參數(shù)的做法較為常見。例如：伍鵬等［4］基于仿真與實驗獲取的應(yīng)力-應(yīng)變數(shù)據(jù)信息構(gòu)建目標函數(shù)，再通過分步反演與Hooke-Jeeves 優(yōu)化算法［10］相結(jié)合的方法，求解了內(nèi)聚力模型相關(guān)參數(shù)。饒玉文等［11］基于仿真與實測的載荷-位移曲線信息構(gòu)建目標函數(shù)，通過人工蜂群優(yōu)化算法，確定了雙線性內(nèi)聚力模型［12］相關(guān)參數(shù)。Chen X 等［13］采用改進的Levenberg-Marquardt 算法［14］，確定了2024-T3 鋁合金斷裂損傷界面內(nèi)聚力模型的參數(shù)，使荷載-裂紋擴展關(guān)鍵點的模擬預(yù)測值與實測值的差異最小化。學(xué)者大多通過載荷-位移相關(guān)數(shù)據(jù)開展反演研究，僅能反映出試件的宏觀力學(xué)行為，無法反映局部信息，而且該方法的約束性較高，需基于應(yīng)變場均勻性假設(shè)。通過感興趣區(qū)域（Region of Interest，ROI）位移場的相關(guān)數(shù)據(jù)信息開展反演研究可對上述問題進行有效改善。同時，基于數(shù)字圖像相關(guān)技術(shù)（Digital Image Correlation，DIC）技術(shù)獲取的位移場/應(yīng)變場中包含了大量的材料響應(yīng)數(shù)據(jù)，在統(tǒng)計意義上使反演獲得的參數(shù)準確性高于常見的載荷-位移構(gòu)造方法。DIC 技術(shù)具有全場非接觸性、適用范圍廣等優(yōu)點，已成為一種有效的表面變形非接觸測量方法，可實現(xiàn)變形區(qū)域位移的量化分析［15］，已廣泛應(yīng)用于力學(xué)實驗領(lǐng)域。李高春等［16］采用DIC 技術(shù)對掃描電子顯微鏡獲取的復(fù)合固體推進劑拉伸變形圖片序列開展計算，獲得了推進劑細觀表面位移場。Gonzalez J 等［17］采用DIC 技術(shù)分析了復(fù)合固體推進劑裂紋尖端附近的應(yīng)變場特點。

針對以載荷-位移曲線信息開展反演研究時，存在的局限性以及數(shù)據(jù)類型單一、數(shù)據(jù)量小的問題。本研究采用基于DIC 技術(shù)結(jié)合Hooke-Jeeves 優(yōu)化算法，以固體火箭發(fā)動機矩形粘接試件單向拉伸實驗中的實測數(shù)據(jù)（ROI 位移場和應(yīng)力-應(yīng)變曲線）為依據(jù)，建立表征粘接界面真實本構(gòu)關(guān)系的雙線性內(nèi)聚力模型，為粘接界面的力學(xué)行為分析提供理論基礎(chǔ)和參考依據(jù)。

1 基本理論

1.1 數(shù)字圖像相關(guān)技術(shù)

DIC 技術(shù)依據(jù)圖像中各像素點周圍區(qū)域的散斑分布情況不同的特點，按照定義的相關(guān)函數(shù)和搜索方法在目標圖像中確定與參考子集相關(guān)性最大的目標子集，從而實現(xiàn)物體表面位移的測量［18-19］。一階位移模式下，參考圖像像素點與目標圖像像素點關(guān)系示意圖，如圖1 所示。

圖1 一階位移模式示意圖Fig.1 Schematic diagram of first-order displacement pattern

參考圖像像素點與目標圖像像素點坐標對應(yīng)關(guān)系，如式（1）：

目標圖像中需確定與參考圖像ROI 相關(guān)性最大的目標區(qū)域，通過變形前后子集的相關(guān)性計算，作為匹配程度的標準，采用歸一化最小二乘相關(guān)函數(shù)［20］表征匹配程度，如式（2）：

式中，C為歸一化最小二乘相關(guān)函數(shù)，數(shù)值越趨近于0，子集匹配程度越高；S 表示子集區(qū)域；f（xi，yi）為參考圖像中像素點（xi，yi）的灰度值；g（x′i，y′i）為目標圖像中像素點（x′i，y′i）的灰度值；fm與gm分別為參考子集與目標子集的平均灰度值。子集的變形矢量P，如式（3）：

在運算求解過程中，目標子集與參考子集不會完全匹配，通過搜索算法尋找相關(guān)函數(shù)的極值，從而確定匹配程度最高的目標子集，對相關(guān)函數(shù)C（p）求關(guān)于ΔP的偏導(dǎo)數(shù)，并令其等于0，如式（4）：

式中，P0為經(jīng)整像素級別的搜索，獲得的初始估計值；??C(p0)為相關(guān)函數(shù)C（p）在初始估計值處的二次偏導(dǎo)，又稱為Hessian 矩陣，可通過牛頓迭代得到方程解。

1.2 內(nèi)聚力模型

粘接界面的本構(gòu)關(guān)系采用雙線性內(nèi)聚力模型［5］，其依次為上升段、下降段和完全脫粘段，上升段和下降段均為線性假設(shè)［12］，如圖2 所示。

圖2 中，σmax為最大粘接強度，MPa，δe、δf分別為損傷起始位移、最大失效位移，mm。

圖2 雙線性內(nèi)聚力模型Fig.2 Bilinear cohesive zone model

粘性界面單元簡化為僅考慮垂直界面的法向變形和沿界面方向的剪切變形，如圖3 所示。

圖3 界面單元變形示意圖Fig.3 Deformation schematic diagramof interface element

δⅠ、δⅡ分別為沿界面方向的法向變形和剪切變形，mm。

單元的總變形，如式（5）：

式中，δm為單元的總變形，mm。

以損傷因子d表征界面損傷情況，損傷因子與界面損傷起始位移的關(guān)系，如式（9）：

可見，雙線性內(nèi)聚力模型的形狀可由最大粘接強度、模量、失效斷裂能所決定。為方便后續(xù)開展仿真計算，將粘性界面單元簡化為各向同性，各物理量的法向、切向數(shù)值設(shè)置相同。

2 力學(xué)實驗

2.1 矩形粘接試件拉伸實驗

2.1.1 試件與實驗過程

按 照QJ2038.1A-2004 標 準［21］的 固 體 火 箭 發(fā) 動機矩形粘接試件為研究對象，該試件由鋼件、三元乙丙橡膠（EPDM）絕熱層、端羥基聚丁二烯（HTPB）/異佛爾酮二異氰酸酯（IPDI）襯層、HTPB 推進劑組成，在下側(cè)鋼件與絕熱層間的左右兩端，預(yù)制了長度為20 mm的人工脫粘層，以達到緩解兩端應(yīng)力集中的目的。

DIC 技術(shù)要求變形體表面須具有豐富的特征，這些特征會對測量的精度有著重要的影響，隨機散斑可提供足夠的特征信息［15］，為了獲得高對比度的隨機散斑，實驗前使用口徑0.5 mm 的水性筆在推進劑表面繪制網(wǎng)格，網(wǎng)格尺寸大致為3.0 mm×1.3 mm。

采用位移加載方式，拉伸角度為0°、速率為5 mm·min-1，使用冷光源持續(xù)提供照明，CCD 相機拍攝頻率為0.5 Hz，圖像的分辨率為4000 pixel×3000 pixel，記錄拉伸載荷與位移數(shù)據(jù)，直至拉伸實驗結(jié)束。拉伸載荷與試件上表面的面積（2000 mm2）比值定義為拉伸應(yīng)力σ，拉伸位移與試件高度（62 mm）比值定義為拉伸應(yīng)變ε。

2.1.2 實驗結(jié)果

粘接界面的面積相對較小，可獲取的信息少且隨機性強，因此采用推進劑本體區(qū)域的位移信息開展反演研究。同時，下側(cè)鋼件與絕熱層兩端，預(yù)制了長度為20 mm 的人工脫粘層，導(dǎo)致相對應(yīng)的推進劑位移信息無法有效反映粘接界面的本構(gòu)關(guān)系。因此以推進劑區(qū)域（20 mm≤x≤80 mm）作為ROI。

應(yīng)力-應(yīng)變曲線的加載段可大致分為3 個階段：加載初期，粘接試件與夾具間存在縫隙，曲線存在一定的非線性，ROI 位移場受縫隙的影響較大，不利于采用該階段的圖像開展ROI 位移場分析；隨后，曲線近似呈線性關(guān)系增長，界面發(fā)生了微小損傷，對試件的力學(xué)性能影響較??；達到一定載荷后，界面出現(xiàn)了不可逆損傷，曲線上升速率逐漸放緩，直至達到應(yīng)力峰值。使用開源數(shù)字圖像相關(guān)分析軟件Ncorr［19］對CCD 相機拍攝的拉伸變形圖片序列中的ROI 開展計算，選擇二、三階段對應(yīng)的ε為0.05、0.08 拉伸狀態(tài)開展ROI 位移場情況分析，如圖4 所示。

兩側(cè)人工脫粘層尖端處產(chǎn)生應(yīng)力集中現(xiàn)象，左側(cè)人工脫粘層尖端處承載能力較右側(cè)弱，ε為0.05 時，左側(cè)人工脫粘層尖端處出現(xiàn)微裂紋，ε為0.08 時，左側(cè)人工脫粘層尖端處的裂紋已大幅度向右擴展。ROI 沿x方向的變形程度小，圖4a、圖4b 中沿x方向的位移呈近似對稱分布，受左側(cè)人工脫粘層尖端處的微裂紋影響不明顯。0°拉伸條件下，試件主要在y方向上產(chǎn)生變形，ε為0.05 時，左側(cè)推進劑沿y方向的位移略高于右側(cè)，圖4c 中沿y方向的位移不呈對稱分布，ε為0.08時，左側(cè)推進劑沿y方向的位移明顯高于右側(cè)，圖4d中的紅色（沿y方向位移較?。﹨^(qū)域向右下側(cè)偏移明顯。

圖4 推進劑區(qū)域（20 mm≤x≤80 mm）的位移云圖Fig.4 Displacement cloud diagrams of propellant region（20 mm≤x≤80 mm）

2.2 HTPB 推進劑松弛實驗

2.2.1 試件與實驗過程

按 照QJ2487-1993 標 準［22］的HTPB 推 進 劑 試 件為研究對象。在室溫20 ℃條件下，對其開展松弛實驗，如圖5 所示。

圖5 HTPB 推進劑試件與夾具安裝裝置Fig.5 HTPB propellant specimen and fixture installation device

2.2.2 實驗結(jié)果

HTPB 推進劑的松弛模量可以由Prony 級數(shù)表示［5］，如式（11）：

式中，E∞為平衡模量，MPa。采用最小二乘法對實驗得到的松弛曲線進行擬合，擬合后的推進劑松弛曲線，如圖6 所示，得到7 級Prony 級數(shù)的各項參數(shù)，見表1。

圖6 擬合后的推進劑松弛模量曲線Fig.6 The relaxation modulus curve of propellant after fitting

表1 推進劑Prony 級數(shù)參數(shù)Table 1 Prony parameters of propellant

3 反演識別與效果分析

3.1 反演識別

在仿真過程中，粘接界面各處的本構(gòu)關(guān)系簡化一致，兩側(cè)人工脫粘層尖端處的破壞程度相同；而對矩形粘接試件的拉伸實驗中，其通常表現(xiàn)為單側(cè)起裂。在圖4 中發(fā)現(xiàn)，ε＞0.08 時，左側(cè)的人工脫粘層尖端處脫粘嚴重，位移云圖不再呈對稱分布，若仍以整個拉伸過程中的ROI 位移信息構(gòu)建目標函數(shù)，則會導(dǎo)致反演精度低。因此采用ROI 位移場信息與應(yīng)力-應(yīng)變曲線數(shù)據(jù)信息共同構(gòu)造目標函數(shù)的方法降低誤差：0≤ε≤0.105 的 實 測 數(shù) 據(jù) 信 息 采 用ε為0.05、0.08 的 位 移 場信息進行表征，0.105＜ε≤0.152 的實測數(shù)據(jù)信息采用應(yīng)力-應(yīng)變曲線（0.105≤ε≤0.152）數(shù)據(jù)信息進行表征。研究采用Hooke-Jeeves 優(yōu)化算法逐步修正待求參數(shù)，反演流程如圖7 所示。主要包括有限元模型的建立、目標函數(shù)的構(gòu)造、內(nèi)聚力模型參數(shù)迭代更新3 部分。

圖7 反演流程圖Fig.7 Inversion flow chart

（1）有限元模型的建立

按照矩形粘接試件的結(jié)構(gòu)和尺寸建立有限元模型，其中推進劑、絕熱層網(wǎng)格的單元類型為C3D8H，鋼件網(wǎng)格的單元類型為C3D8，襯層網(wǎng)格的單元類型為COH3D8，模型中共2556 個節(jié)點，如圖8 所示。

圖8 有限元模型Fig.8 Finite element model

絕熱層和鋼件在拉伸實驗中的變形小，視為線彈性材料，其力學(xué)性能參數(shù)見表2。

表2 力學(xué)性能參數(shù)［5］Table 2 Mechanical properties parameters

HTPB 推進劑為粘彈性材料，其松弛模量可以由Prony 級數(shù)表示，采用2.2 節(jié)的松弛實驗數(shù)據(jù)。

（2）目標函數(shù)的構(gòu)造

Hooke-Jeeves 優(yōu)化算法由目標函數(shù)探索設(shè)計空間，通過目標函數(shù)表征仿真結(jié)果與實測結(jié)果的相近程度，數(shù)值越小，則表明二者越接近。設(shè)置相同時間間隔對仿真、實測應(yīng)力-應(yīng)變曲線進行離散，進而構(gòu)造目標函數(shù)R（p），如式（12）：

式中，p為內(nèi)聚力模型相關(guān)參數(shù)組成的向量集合；U為位移，mm；x、y分別為沿界面平行、垂直的方向；a、b為放大/縮減因子；m、n分別為ROI 位移場與應(yīng)力-應(yīng)變曲線目標點的數(shù)量；i、j為目標點的編號；sim、exp 代表仿真情況與實測情況；k 表示拉伸狀態(tài)（ε=0.05、0.08）。

（3）內(nèi)聚力模型參數(shù)迭代更新

在試算時發(fā)現(xiàn)：最大粘接強度影響應(yīng)力-應(yīng)變曲線的峰值、模量影響曲線的上升段、失效斷裂能影響曲線的下降段。經(jīng)有限元試算工作，確定最大粘接強度、模量、失效斷裂能的初值，如表3 所示。

表3 相關(guān)參數(shù)的初值與最優(yōu)值Table 3 Initial and optimal values of relevant parameters

Hooke-Jeeves 優(yōu)化算法的初始步長為0.05，閾值為0.01，縮減因子為0.5，步長加速因子為1.1，目標函數(shù)閾值設(shè)為2.5。反演過程中，相關(guān)參數(shù)以及目標函數(shù)的變化情況，如圖9 所示。

圖9 相關(guān)參數(shù)與目標函數(shù)的變化情況Fig.9 Changes of relevant parameters and objective function

反演過程共開展19 輪“探測移動”，合計124 次計算。最終，增量步長縮減為6.25×10-3，小于設(shè)定的步長閾值。目標函數(shù)曲線在第85 次計算后基本保持平穩(wěn)，收斂情況較好，在第116 次賦值試算中的目標函數(shù)值最小，表明該組參數(shù)下的仿真結(jié)果與實測結(jié)果最為接近，即最優(yōu)值，如表3 所示。

3.2 反演效果分析

基于3.1 節(jié)的反演識別工作，獲取了表征粘接界面力學(xué)性能的最優(yōu)雙線性內(nèi)聚力模型。需通過與相應(yīng)的實驗數(shù)據(jù)進行對比，從而分析最優(yōu)雙線性內(nèi)聚力模型的準確性。通過單向拉伸過程中的應(yīng)力-應(yīng)變曲線信息和ROI 位移誤差云圖信息，對最優(yōu)雙線性內(nèi)聚力模型的準確性展開分析。

基于初始雙線性內(nèi)聚力模型和最優(yōu)雙線性內(nèi)聚力模型開展仿真計算，相對應(yīng)的仿真應(yīng)力-應(yīng)變曲線以及實測應(yīng)力-應(yīng)變曲線，如圖10a 所示。初始雙線性內(nèi)聚力模型和最優(yōu)雙線性內(nèi)聚力模型的應(yīng)力數(shù)據(jù)與實測應(yīng)力數(shù)據(jù)的比值變化圖，如圖10b 所示。

圖10 仿真和實測應(yīng)力-應(yīng)變相關(guān)曲線Fig. 10 Simulated and measured stress-strain correlation curves

圖10a中，初始應(yīng)力-應(yīng)變曲線的上升斜率與實驗曲線的差距較小，但在整個拉伸過程中，應(yīng)力數(shù)據(jù)普遍高于實測數(shù)據(jù)，可見，初始雙線性內(nèi)聚力模型并不能有效的表征粘接界面的真實力學(xué)性能。最優(yōu)雙線性內(nèi)聚力模型的應(yīng)力-應(yīng)變曲線與實測曲線具有相同的變化趨勢與規(guī)律，相似度較高。圖10b中，在拉伸階段（0≤ε＜0.11）時，初始與最優(yōu)的應(yīng)力數(shù)據(jù)比值均近似為1；在拉伸階段（0.11≤ε＜0.15）時，初始應(yīng)力數(shù)據(jù)比值大幅度提高，最大值達到14.6；在拉伸階段（0.11≤ε＜0.14）時，最優(yōu)的應(yīng)力數(shù)據(jù)比值近似為1，在拉伸階段（0.14≤ε＜0.15）時，最大比值僅為3.0。

以相對誤差r反映仿真、實測應(yīng)力-應(yīng)變曲線的相近程度，如式（13）：

式中，s為應(yīng)力-應(yīng)變曲線與橫坐標軸圍成的面積，MPa。

初始曲線的相對誤差為44.7%，最優(yōu)曲線的相對誤差降低為4.3%?？梢?，該反演識別方法可以對內(nèi)聚力模型進行有效修正，反演效果理想。

ε為0.05、0.08 拉伸狀態(tài)下的仿真位移云圖，如圖11 所示。

圖11 仿真位移云圖Fig.11 Simulation displacement cloud images

在仿真過程中，由于將粘接界面各處的本構(gòu)關(guān)系簡化一致，兩側(cè)人工脫粘層尖端處的損傷情況相同，并且沿x、y方向的仿真位移云圖均呈對稱分布，與實測位移云圖的分布特征均基本一致。ε為0.05 時，兩側(cè)的人工脫粘層尖端已大幅度沿拉伸方向運動，粘接界面附近區(qū)域沿y方向大致進行了0.5 mm 位移，兩側(cè)的人工脫粘層有效的緩解應(yīng)力集中現(xiàn)象的發(fā)生。ε為0.08時，粘接界面附近區(qū)域沿y 方向大致發(fā)生了1.4 mm 位移，此時粘接界面損傷嚴重。

以絕對誤差r′反映仿真、實測ROI 位移場的相近程度，如式（14）。ROI 位移誤差云圖，如圖12 所示。

圖12 ROI 位移誤差云圖Fig.12 Cloud maps of ROI displacement error

ε為0.05 時，ROI 位 移 誤 差 云 圖 的 最 大 誤 差 為0.64 mm、平均誤差為0.38 mm，其誤差分布相對均勻。ε為0.08 時，ROI 位移誤差云圖的最大誤差為1.76 mm、平均誤差為0.45 mm，上半?yún)^(qū)域（7≤y≤31）的誤差相對較小，且誤差分布相對均勻，由左側(cè)的人工脫粘層處最先發(fā)生破壞，導(dǎo)致ROI 位移誤差云圖中下半?yún)^(qū)域（31＜y≤55）左側(cè)的誤差相對較大。

以上分析均表明采用反演方法獲取的雙線性內(nèi)聚力模型可以用于表征粘接界面的真實本構(gòu)關(guān)系。

經(jīng)分析，誤差主要來源于5 方面：1、推進劑實際本構(gòu)關(guān)系復(fù)雜，對于拉伸載荷較大時，推進劑所采用本構(gòu)參數(shù)不能完全表征其真實的力學(xué)性能，一定程度上影響了反演精度；2、對鋼件和絕熱層的力學(xué)性能予以簡化處理，通過引用相關(guān)文獻的方式確定二者的力學(xué)性能參數(shù)，會與實際的力學(xué)性能存在偏差，但較第一條誤差來源相比占次要因素；3、試件表面不平整，同時，CCD 相機在拍攝過程中，各圖片的亮度和對比度很難保證完全一致等因素，會影響圖片的灰度分布，進而影響位移場測量結(jié)果。4、仿真過程中，粘接界面作簡化處理，各處的本構(gòu)關(guān)系一致，而在矩形粘接試件的拉伸實驗中，其通常為單側(cè)起裂，粘接界面兩側(cè)的力學(xué)性能實際不完全相同。5、可通過適當(dāng)增加其他拉伸狀態(tài)下的位移場信息的方式，豐富拉伸階段（0≤ε≤0.105）中的實測數(shù)據(jù)信息，以達到更好的表征的效果。

4 結(jié)論

（1）基于DIC 技術(shù)與Hooke-Jeeves 優(yōu)化算法相結(jié)合的反演識別方法，可以準確地獲取固體火箭發(fā)動機粘接界面模型，拉伸速率為5 mm·min-1時，雙線性內(nèi)聚力模型的最大粘接強度、模量、失效斷裂能分別為0.55 MPa、0.57 MPa、2.26 kJ·m-2。

（2）仿真與實測應(yīng)力-應(yīng)變曲線的相對誤差由初始44.7%修正為4.3%，拉伸應(yīng)變?yōu)?.05、0.08 狀態(tài)下的ROI 最大位移誤差分別為0.64 mm、1.76 mm，平均位移誤差分別為0.38 mm、0.45 mm，均表明所采用的反演識別方法的精度較高，建立的雙線性內(nèi)聚力模型可以一定程度上表征粘接界面的真實本構(gòu)關(guān)系。

（3）本研究為粘接界面模型的精準獲取提供了一個新的思路，并對粘接界面的力學(xué)行為分析提供理論基礎(chǔ)和參考依據(jù)。在后續(xù)工作中，通過對五點誤差來源方式加以改進，會使反演精度有效提高。