细观尺度下推进剂/衬层/绝热层粘接界面参数反演研究

肖云东,王玉峰,李何龙,赖帅光,伍 鹏,韩永恒

(1.海军航空大学, 山东 烟台 264001; 2. 91526部队, 广东 湛江 524000;3. 91468部队, 海南 陵水 572400;4.海军装备部驻北京地区军事代表局, 北京 100071)

0 引言

粘接界面的力学性能是影响固体火箭发动机寿命的关键因素之一[1]。随着对粘接界面损伤研究的逐渐深入,众多学者从细观角度出发,运用细观均匀化理论和有限元法定量分析界面的力学性能,这有利于掌握界面破坏的本质。

有限元法中通常采用内聚力模型描述界面的力学性能,准确的内聚力模型是模拟界面失效过程的前提。关于如何确定内聚力模型的问题中,大多预先定义模型种类,再通过经验法、反演识别(有限元模型修正技术)等方法确定模型参数。针对细观尺度下的界面脱粘问题,由于推进剂内部颗粒夹杂,粘接界面处的组成复杂,采用经验法则很难保证内聚力模型的准确性。

有学者采用反演识别方法,即基于仿真与实测结果构造目标函数,通过最小化目标函数,迭代识别细观尺度下的内聚力模型参数。封涛[2]与职世君[3]建立了复合固体推进剂二维细观填充模型,均以仿真与实测的应力-应变数据信息构造目标函数,分别通过Hooke-Jeeves优化算法与分步迭代法对颗粒/基体界面之间的内聚力模型参数进行了反演识别。但其目标函数信息中仅能反映出试件的力学行为,无法反映出试件的局部信息。

通过数字图像相关技术(digital image correlation,DIC)可获取试件全场二维/三维的变形信息,可实现位移和应变的量化分析[4]。李高春[5]使用SEM对复合固体推进剂拉伸变形破坏过程进行了观察,并通过DIC技术获取了推进剂细观表面位移场。基于DIC技术的反演识别方法中包含了大量的材料响应数据,在统计学角度上,该反演识别方法的准确性更高。Fedele[6]采用DIC技术与有限元法相结合的反演识别方法对眩光层压板的脱胶实验结果拟合,确定了指数型内聚力模型参数。杨思满[7]通过DIC技术获取试件表面的位移场,并与仿真位移场构造目标函数,采用遗传算法和序列二次规划法相结合的策略对材料参数进行反求。

为提高细观尺度下粘接界面内聚力模型的准确性,通过SEM原位拉伸实验,记录粘接试件在拉伸过程中的变形与破坏,根据界面附近颗粒分布情况建立了简化后的粘接界面细观模型,采用DIC技术结合Hooke-Jeeves 优化算法对界面所采用的双线性内聚力模型相关参数开展反演识别研究。

1 基于DIC技术的反演识别方法

1.1 反演思路

通过相关文献[8],预先定义粘接界面所采用的模型为双线性内聚力模型,经过大量的试算工作确定初始模型的相关参数;采用有限元法和DIC技术分别获取ROI的位移场信息;基于仿真与实测位移场信息构建目标函数,并采用Hooke-Jeeves优化算法沿着目标函数减小的趋势逐步修正待求参数。反演流程如图1所示。

图1 反演流程框图

1.2 Hooke-Jeeves优化算法

Hooke-Jeeves优化算法又称模式搜索算法,通过目标函数来探索设计空间,可对待求的参数按照正确的“趋势”赋值,具有适定性强、收敛速度快等优点。算法基于所赋的初值点开始,由“探测移动”沿着不同轴不同方向进行探测性移动,从而确定目标函数值下降的方向,再通过所设定的“模式移动”,使目标函数值快速沿着“山谷”方向,逐渐逼近目标函数最小值。恰当的初始基点可大幅度提高反演精度,降低计算成本。通过“探测移动”与“模式移动”多次交替运行,直至最终符合收敛条件时,计算终止,从而实现寻优目的[9]。

2 细观拉伸实验

制作了厚度为3 mm的微型非标准试件,如图2所示。其两端为有机玻璃,中间部分是通过切取矩形粘接试件得到的推进剂/衬层/绝热层粘接件,使用302胶将3部分黏结,并自然固化2 d。

图2 微形非标准试件及其尺寸

通过日本岛津JSM-5410LV型试验系统(全数字电液伺服加载系统和高精度的扫描电子显微镜)进行细观拉伸实验,并实时记录试件在拉伸载荷作用下的力学行为与表面细观结构的变化情况。

拉伸过程由伺服控制器自动加载,加载速率设定为1.2 mm·min-1。SEM放大倍数为50,调节对比度和亮度后,可清晰成像,且图片特征明显。拉伸位移与推进剂/衬层/绝热层粘接件边长的比值定义为拉伸应变ε,衬层附近区域的破坏过程如图3所示。

图3 粘接界面破坏过程(×50)

其中图3(b)、图3(c)是ε为8%、18%时,通过开源数字图像相关分析软件Ncorr获取图3(a)虚线内沿着x方向的实测位移(Uexpx)云图。图3(a)中,2条白色虚线分别位于绝热层/衬层界面与衬层/推进剂界面处,d0为初始衬层厚度。Δdi表示拉伸过程中的界面张开位移:

Δdi=di-d0

(1)

式中,i=1、2、3分别表示ε为8%、18%、30%的受载状态。拉伸过程中的衬层厚度与界面张开位移情况,见表1所示。

表1 衬层厚度与界面张开位移Table 1 Interlayer thickness and opening displacement μm

3 有限元建模

3.1 细观模型建立

推进剂/衬层/绝热层粘接试件的代表性体积单元(representative volume element,RVE)中包括HTPB推进剂、HTPB/IPDI衬层、EPDM绝热层,HTPB推进剂中又包含了大量的固体燃料等颗粒。RVE在宏观上尺寸需足够小,可看成一个物质点,满足连续介质力学的基本假定;在细观上尺寸需足够大,包含足够多的细观元素和细观结构信息,可代表局部连续介质的统计平均性质。模型中的各部分尺寸,如表2所示。

表2 模型尺寸

其中,基体尺寸满足推进剂代表性体积单元最小尺寸680 μm×680 μm的要求[10]。Al颗粒粒径一般为10～30 μm,AP颗粒粒径总体呈双峰分布规律,小颗粒平均粒径约为20 μm,大颗粒平均粒径在100～300 μm,各颗粒的粒径相差大,若将其全部构造,会极大程度上增加计算量,且不易收敛。因此,采用多步等效法将粒径小于80 μm的颗粒等效为基体的一部分,其余颗粒简化为圆形或椭圆形,通过二维图像重构的方式[11]建立粘接界面细观模型,如图4所示。

图4 粘接界面细观模型

粘接界面的脱粘行为与颗粒/基体的脱湿行为分别通过黏结单元、黏结接触开展仿真计算。黏结单元较为常见,不过多介绍。黏结接触是通过从属面与主控面之间的接触从而定义损伤模型,由从属面上的点对主控面进行投影,从而确定接触点对,若对应的接触点不在网格节点处,可通过邻近节点的计算结果进行插值获取,具有建模简单、易收敛等优点[12]。

3.2 本构关系

双线性内聚力模型的形式简单、且易于实现,在界面脱粘问题中被广泛采用[13]。粘接界面脱粘行为与颗粒/基体脱湿行为均采用双线性内聚力模型,如图5所示。

图5 双线性内聚力模型

n为垂直界面的法向方向,s为沿界面方向的切向方向,δ0为界面损伤起始位移,δf为界面失效位移。

当δ≤δ0时,界面未发生损伤,界面力T随界面张开位移的增大而线性增大,其关系如式(2):

(2)

式中,K为界面刚度。

界面损伤起始准则采用最大名义应力准则,认为界面力达到界面的最大名义应力时,界面开始发生损伤:

(3)

式中:σmax、τmax为法向与切向的最大名义应力;〈〉为Macaulay括号,表示纯压应力状态不会引发损伤。

采用基于位移损伤演化准则,当δ0≤δ≤δf时,界面发生损伤,界面力随着界面张开位移增加而线性减小:

(4)

式中,D表征界面损伤状态,对应仿真软件的输出量为SDEG。当D=0时表明界面未损伤,D介于0～1时表明界面处于损伤阶段,D=1时表明界面完全失效开裂:

(5)

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

(6)

可以看出,双线性内聚力模型的形状可由界面刚度、最大名义应力、界面失效位移所决定,为了方便后续仿真计算,模型简化为各向同性,将各物理量法向与切向的数值设置相同。基于文献[2]中的部分界面模型参数,又经过大量的有限元试算,最终各界面所采用的模型参数,如表3所示。

表3 各界面的模型参数

将绝热层、基体、颗粒简化为弹性材料,对绝热层试件进行单向拉伸实验,取其初始模量作为弹性模量,其泊松比以及基体和颗粒的力学性能参数参考文献[14],如表4所示。

表4 粘接结构各部分力学性能参数

3.3 网格划分与边界条件

绝热层与颗粒的杨氏模量较基体相比要大得多,载荷较小的情况下可近似认为二者不发生变形,对其采用四结点双线性平面应变四边形单元CPE4;基体作为近似不可压材料,受载时体积不发生变化,采用四节点四边形线性积分的杂交单元CPE4H来模拟其在大变形下的力学响应;衬层采用四节点二维粘接单元COH2D4。各界面附近的网格需适当细化以提高计算的准确性和收敛性。

粘接试件在拉伸载荷的作用下,Y方向变形小,对边界条件进行简化,仿真过程中不考虑试件沿Y方向的变形。将A、C边设置为沿Y方向固定,D边设置为完全固定,B边沿X正方向受均布拉伸载荷,速率为1.2 mm/min。网格划分与边界条件的设置情况,如图6所示。

图6 网格划分与边界条件

4 参数反演

4.1 反演过程与结果

将部分衬层与推进剂区域作为ROI,为图3(a)蓝色虚线内的区域。文献[15]中指出界面模型会对试件的力学行为产生影响,界面本构关系也势对影响推进剂表面的变形特征,同时,通过细观拉伸实验发现:界面的本构关系主要影响ROI沿x方向的位移场分布,对沿y方向的位移场分布的影响可忽略。

以ε为8%、18%时的位移场(Ux)信息共同构造目标函数,构造方法参考文献[16]并改进:

(7)

式中:p为所有界面参数组成的向量集合;k为拉伸状态;i为目标点的编号;m为目标点的数量。

反演程序中的增量步长设为0.05、步长缩减因子设为0.5、步长加速因子设为1、增量步长阈值设为5×10-3,粘接界面的内聚力模型参数初值,见表1所示。

迭代过程中,界面刚度、最大名义应力、界面失效位移以及目标函数值的变化情况,如图7所示。

图7 迭代过程曲线

迭代过程共计100次赋值计算,目标函数曲线在第83次计算后,开始逐渐平稳,各曲线的收敛情况较好。在第16轮的“探测移动”中,增量步长为3.12×10-3小于设定的增量步阈值,完成此轮“探测移动”后,计算终止。第98次赋值试算中的目标函数值最小,认为此界面模型下的ROI的仿真位移与实测位移最为接近,即最终的反演结果,如表5所示。

表5 反演结果值

4.2 反演效果分析

反演获取的模型参数的准确性需通过相应的实验数据进行验证。基于反演结果开展仿真计算,将仿真获取的ROI位移场与实测数据进行比较。ε为8%、18%时,ROI的仿真位移(Usimx)云图,如图8所示,相对应的实测位移(Uexp x)云图,见图3(b)、图3(c)。可见2种边界条件下的位移(Ux)云图分布特征均基本一致,相似度极高。

各节点处的位移误差ηU:

(8)

ε为8%、18%时,ROI的位移(Ux)误差云图,如图9所示。

图9 位移(Ux)误差云图

ROI沿x方向的位移误差统计情况,见表6所示。

主要误差因素有4个方面:① 界面附近存在大量微裂纹与微孔洞等缺陷,部分颗粒存在破碎与包裹不完全等缺陷,导致ROI的下区域平均位移(Uexpx)略高于上区域,而细观模型中未对其考虑;② 各图片亮度和对比度在获取过程中很难保证完全一致,会影响图片的灰度分布,进而影响Ncorr计算结果;③ (1～7)颗粒附近区域的变形梯度大,Ncorr软件中出现匹配相关性差、测量精度低的问题;④ 仿真过程中,将各材料简化为弹性,部分材料的力学性能参数直接引用相关文献,与实际力学性能存在偏差。

表6 ROI位移误差情况(Ux)

5 界面损伤分析

采用界面单元的SDEG对界面的损伤程度进行表征,使用Python脚本语言对仿真软件进行2次开发,提取内聚力界面单元高斯积分点的SDEG数据,绘制SDEG变化曲线,如图10所示。

图10 损伤因子变化曲线

起初在外载荷的作用下,界面力随着界面张开位移的增加呈线性增长,SDEG均为0,未进入损伤阶段;界面力达到最大值后,界面处损伤开始萌生并扩展;随着界面张开位移的逐渐增加,裂纹逐步形成扩展,界面承受载荷的能力降低;当界面力减小至0时,界面损伤达到高峰,SDEG为1,此时界面完全损伤。

结合粘接界面破坏过程(见图3),以及获取的法向内聚力模型可知:当ε=8%时,由于推进剂的杨氏模量较小,其在载荷作用下的变形明显,此时界面张开位移为75 μm,SDEG为0,界面力为0.48 MPa;当界面张开位移为112.5 μm时,界面力达到最大值0.72 MPa;当ε=18%时,(1～7)颗粒脱湿现象明显,且6号、7号颗粒脱湿情况额外严重,脱湿处的承载能力降低,界面附近区域产生微裂纹,导致(Uexpx)云图中产生像素点不匹配的无位移输出区域,此时界面张开位移为150 μm,大于界面损伤起始位移,界面进入损伤阶段,表现出线性软化行为,SDEG为0.26,界面力为0.70 MPa;当ε=30%时,界面与颗粒/基体附近的微裂纹汇合产生宏观裂纹,且宏观裂纹的扩展方向几乎与拉伸方向垂直,此时界面张开位移为375 μm,SDEG为0.75,界面力为0.61 MPa。

6 结论

1) 根据双线性内聚力模型中的参数特性,采用基于DIC技术与Hooke-Jeeves优化算法相结合的反演识别方法,实现了对细观尺度下粘接面损伤模型(双线性内聚力模型)中的界面刚度、最大名义应力、界面失效位移的求解。当加载速率为1.2 mm·min-1时,模型参数分别为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%的区域占比分别为78.6%、85.7%,均表明仿真结果与实验吻合较好,所建立的界面模型可以在一定程度上表征界面真实的本构关系。

3) 获取了界面单元的损伤因子变化曲线,与粘接试件拉伸过程的损伤情况有较好地吻合,验证了所建立的界面模型可正确地反映出粘接界面在拉伸过程中脱粘规律。