0 引言
固体火箭发动机具有结构简单[1]、可靠性高[2]、贮存和运输方便[3]等显著优点，在运载火箭与导弹武器系统等方面应用十分广泛。满装填技术是未来固体动力系统发展的重要方向之一[4]。端面燃烧满装填固体发动机是指采用贴壁浇注工艺实现燃面沿轴向稳定退移的发动机，这种发动机通过高装填比技术最大限度利用燃烧室容积，显著提升发动机质量比与总冲。这类发动机的一个显著特点是：在发动机工作过程中燃面退移速度与包覆套烧蚀消融速度不匹配，随着发动机工作时间的增加，暴露于高温高压燃气中的包覆套/推进剂界面内外存在压差并产生额外应力，当应力超过界面粘接强度时，界面将发生损伤甚至剥离，导致端燃药柱外径区域燃面异常暴露、引起整体燃面突然扩展，这可能会导致燃烧室压强随之迅速升高直至发动机爆炸，这也是界面脱粘现象一直为研究人员所关注的原因[5]。该界面在服役载荷下的响应特性与失效过程极为复杂，除热-力载荷外，还受到内流场冲刷及压力振荡等多场耦合作用[6-7]。特别是在长时间工作的端燃发动机中，界面损伤是一个随时间累积的动态演化过程[8]。
目前，光学[9]、超声[10]、X射线[11]和脉冲热成像[12]等多种无损检测方法均可对推进剂、绝热层、包覆套以及界面处的缺陷进行检测。刘斌等[13]通过剪切散斑干涉技术开展了包覆药柱的粘接质量评估。李光保等[14]研发了一种用于包覆药柱缺陷检测的双轴运动激光散斑检测系统。邓亮亮等[15]提出了一种基于局部共振谱无监督深度学习算法的大型固体火箭发动机界面脱粘检测方法。结合激光剪切散斑干涉技术和有限元数值模拟方法，王硕等[16]对比分析了三种类型缺陷在真空负压加载及热加载条件下的变形机理。在数值模拟辅助检测方面，蒙上阳等[17]过三维有限元模型计算不同深度脱粘的应力强度因子，获得其随深度变化规律，从而判断界面裂纹稳定性；而在界面力学行为预测上，基于内聚力模型的渐进损伤分析已成为研究界面脱粘扩展的有效手段[18]。此外，针对界面失效过程的精细化实验观测也为数值模型的验证提供了重要支撑[19]。
固体火箭发动机燃烧室内界面失效研究方面，王松柏[20]通过对某端面燃烧固体火箭发动机的早期失败案例分析认为，失败原因是金属壳体与复合推进剂线膨胀系数不同导致界面脱粘。庞陶行等[21]提出整体人工脱粘法用于满装填发动机的新型防护。夏志全[22]分析了贴壁浇注满装填发动机尾端部人工脱粘层脱粘问题，提出了有效的解决方法。邢耀国等[23]根据脱粘面未扩展时内弹道仿真曲线，预估了某端燃发动机脱粘面可能发生扩展时的临界尺寸。姜利等[24]和曹熙炜等[25]开展了自由装填式发动机热匹配失效模式机理研究，确定了包覆套内外压差将导致界面发生损伤甚至破坏，直至发动机爆炸这一现象，并进行了试验验证。Liu等[26]研究了满装填发动机的力学环境适应性，并提出了相应的关键技术。近年来，基于扩展有限元法（XFEM）的裂纹扩展模拟方法也被应用于推进剂界面损伤分析中[27]。同时，考虑温度场与结构场双向耦合的数值方法能够更准确地再现发动机工作过程中的界面热力行为[28]。值得注意的是，随着计算能力的提升，考虑流动与结构相互作用的流固耦合数值方法正被逐步应用于固体发动机工作过程模拟中[29]，这为揭示界面在真实工况下的动态失效机理开辟了新的途径。
总体来看，现阶段相关研究多集中于发动机工作前界面粘接性能评价及脱粘检测等方面，而针对真实工况下端面燃烧满装填固体发动机包覆套与推进剂界面在高温高压环境中、考虑流固耦合效应的响应与损伤特性研究仍显不足。在发动机工作过程中，包覆套与推进剂粘接界面的受力状态、变形特征以及局部脱粘损伤的发展情况，直接关系到发动机点火成败和工作可靠性。同时，包覆套与绝热层间缝隙宽度的变化会改变包覆套与推进剂之间的相互作用关系，进而影响界面的载荷传递、受力分布及变形响应；包覆套弹性模量的差异则会进一步影响其在外部载荷作用下的响应特征及界面损伤行为。因此，为深入分析发动机工作过程中包覆套与推进剂界面的响应与损伤特性，本文基于 LS-DYNA 平台建立了考虑发动机内流场与结构件之间流固耦合作用的数值分析模型，对发动机工作过程中包覆套与推进剂界面的载荷传递、应力响应、变形特征及损伤情况进行了研究，从而为发动机结构设计、工艺设计及材料选取提供参考。
1 物理数学模型
1.1 界面损伤模型
固体火箭发动机由多个界面组成，界面的失效通常通过内聚力法则来描述。内聚力法则基于弹塑性断裂力学，假设界面满足牵引力的连续性关系，从而建立了裂纹面间的拉力与张开位移之间的联系。为了估计相邻层之间每个界面损伤开始的时间，在包覆套和推进剂之间引入了内聚力单元。该单元通过应力破坏准则来描述损伤的起始时刻，主要采用双线性内聚力模型来模拟开裂行为，其满足的牵引分离法则如式（1）所示。
σ（δ）=$\left\{\begin{array}{lll}\frac{{\sigma }_{max}}{{\delta }_{c}}\delta & 0\le \delta \le {\delta }_{c}& \\ {\sigma }_{max}\left(1-\frac{\delta -{\delta }_{c}}{{\delta }_{f}-{\delta }_{c}}\right)& {\delta }_{c}\delta \le {\delta }_{f}& \\ 0& \delta {\delta }_{f}& \end{array}\right.$
式中：σ（δ）为界面处的应力；σmax为界面的破坏应力；δc为临界位移，表示损伤开始的位移；δf为失效位移，表示界面完全失效时的位移。根据上述牵引-分离法则，当界面单元的应力达到σmax时，界面开始发生损伤。
1.2 发动机模型与网格划分
发动机工作过程中，燃气向绝热层与包覆套之间的间隙中填充。然而，由于流通截面限制和包覆套变形等因素，包覆套内外存在压差。为了将燃烧室内的燃气沿轴向导入药柱包覆套与壳体绝热层间的环形缝隙区域以减小包覆套内外压差，在包覆套与绝热层之间均匀布置4根导气条，其结构尺寸如图3所示。
表1 固体域材料参数表Table 1 Material parameters of solid domain |
| 材料 参数 | 包覆套 | 推进剂 | 其他 （刚体） |
|---|---|---|---|
| 密度（kg/m3） | 1260 | 1795 | — |
| 弹性模量（MPa） | 3.5 | 7 | — |
| 泊松比 | 0.4995 | 0.4995 | — |
表2 流体域材料参数表Table 2 Material parameters of fluid domain |
| 密度（kg/m3） | 1.225 |
|---|---|
| 定容比热J/（kg·K） | 713.5 |
| 定压比热J/（kg·K） | 1001.5 |
表3 缝隙宽度及包覆套弹性模量参数表Table 3 Gap width and elastic modulus of coating sleeve |
| 编号 | 缝隙宽度（mm） | 包覆套弹性模量（MPa） |
|---|---|---|
| H1 | 0.5 | 20 |
| H2 | 1 | 20 |
| H3 | 1.5 | 20 |
| H4 | 2 | 20 |
| H5 | 2.5 | 20 |
| P1 | 2 | 20 |
| P2 | 2 | 40 |
| P3 | 2 | 60 |
| P4 | 2 | 80 |
| P5 | 2 | 100 |
表4 界面内聚力单元参数表Table 4 Interface cohesive element parameters |
| 法向强度 （MPa） | 剪切强度 （MPa） | 极限位移 （mm） | 临界位移 （mm） |
|---|---|---|---|
| 1.47 | 1.52 | 3.098 | 0.1 |
固体域及流体域网格均采用Hypermesh软件进行划分，在壁面附近对流体域网格进行了局部加密以提高边界层区域流体参数的计算精度，流体域和固体域的网格划分如图4所示。
1.3 数值模拟方法与边界条件设置
LS-DYNA中的守恒元/解元方法（Conservation Element/Solution Element，CE/SE）把时间和空间统一起来并通过引进守恒元CE和解元SE来使格式在时间和空间上始终保持守恒。三维空间中守恒定律的积分形式为
$\frac{\partial }{\partial t}{\int }_{\Omega }^{}\rho \varphi dV+{\int }_{\partial \Omega }^{}\rho \varphi \overrightarrow{u}·\overrightarrow{n}dS={\int }_{\Omega }^{}{S}_{\varphi }dV$
其中，Ω是控制体积，∂Ω是边界，ρ是密度，ϕ是守恒变量，$\overrightarrow{u}$=（u，v，w）是速度矢量，$\overrightarrow{n}$是边界∂Ω的单位外法向向量，Sϕ是源项。
将三维时空域离散为守恒元CE和解元SE。对于一个三维矩形守恒元C${E}_{i，j，k}^{n}$，在空间上边长分别为Δx、Δy、Δz，时间步长为Δt。守恒元上的守恒方程离散形式为
$\begin{array}{l}\frac{{\varphi }_{i,j,k}^{n+1}-{\varphi }_{i,j,k}^{n}}{\Delta t}+\frac{{F}_{i+1/2,j,k}^{n}-{F}_{i-1/2,j,k}^{n}}{\Delta x}+\\ \frac{{G}_{i,j+1/2,k}^{n}-{G}_{i,j-1/2,k}^{n}}{\Delta y}+\frac{{H}_{i,j,k+1/2}^{n}-{H}_{i,j,k-1/2}^{n}}{\Delta z}={S}_{i,j,k}^{n}\end{array}$
其中，${\varphi }_{i，j，k}^{n}$是守恒变量ϕ在（i，j，k）位置，n时刻的值，F、G、H分别是x、y、z方向上的数值通量，${S}_{i，j，k}^{n}$是源项在（i，j，k）位置，n时刻的值。
x、y、z方向上数值通量计算公式分别为
$\begin{array}{l}{F}_{i+1/2,j,k}^{n}=\frac{1}{2}[{\left(\rho \varphi u\right)}_{i+1,j,k}^{n}+{\left(\rho \varphi u\right)}_{i,j,k}^{n}]-\frac{1}{2}\alpha ({\varphi }_{i+1,j,k}^{n}-{\varphi }_{i,j,k}^{n})\\ {G}_{i,j+1/2,k}^{n}=\frac{1}{2}[{\left(\rho \varphi v\right)}_{i,j+1,k}^{n}+{\left(\rho \varphi v\right)}_{i,j,k}^{n}]-\frac{1}{2}\beta ({\varphi }_{i,j+1,k}^{n}-{\varphi }_{i,j,k}^{n})\\ {H}_{i,j,k+1/2}^{n}=\frac{1}{2}[{\left(\rho \varphi w\right)}_{i,j,k+1}^{n}+{\left(\rho \varphi w\right)}_{i,j,k}^{n}]-\frac{1}{2}\gamma ({\varphi }_{i,j,k+1}^{n}-{\varphi }_{i,j,k}^{n})\end{array}$
其中，α、β、γ是流动特性相关的参数，用于保证算法的稳定性和准确性。
1.4 仿真模型验证
流体域内测点示意如图7所示，其中测点1与燃面后端齐平，各测点均匀分布，间隔距离100mm。通过将仿真获得的测点压强数据与试验结果进行对比来验证模型的有效性。测点1和包覆套与药柱界面第一个单元（见图4）位置接近，位于包覆套与药柱的初始界面位置，具有代表性，与试验结果进行对比能够对仿真模型的正确性进行验证。测点1处压强随时间变化曲线的计算值与试验值对比如图8所示，可以看出，计算结果与试验值变化趋势接近，最大相对误差为15.4%，这表明本文建立的模型能正确模拟点火建压过程中流体域和固体域之间的流固耦合过程。0.4s以前模拟值小于试验值的原因是由于建模时将包覆套简化为理想弹性体，忽略了其实际粘弹性在瞬间压力冲击下的高模量响应，导致模拟计算出的缝隙变形偏大、压力上升偏缓。0.4s以后模拟值大于试验值是因为数值模拟未考虑实际工作中包覆套烧蚀剥落导致的边界改变与容积增大效应。后续的研究中我们将进一步考虑燃面退移、包覆套烧蚀脱落及材料粘弹性等复杂因素，以提高模型的预测精度。
2 结果分析与讨论
2.1 不考虑流固耦合时发动机内流场分析
首先开展不考虑流固耦合时发动机内流场分析计算，其中固体域材料均视为刚体。
2.1.1 包覆套与绝热层间有缝隙
包覆套与绝热层之间有缝隙时，燃气可通过导气条自身的沟槽及包覆套与绝热层之间的缝隙向包覆套与绝热层之间的空腔内填充。包覆套与绝热层之间有缝隙时不同测点处的压强随时间变化过程如图9所示，可以看出，此时各测点处压强变化数据几乎完全一致，这说明包覆套内外基本不存在压差。
2.1.2 包覆套与绝热层间无缝隙
包覆套与绝热层之间无缝隙时，燃气仅能通过导气条自身的沟槽向包覆套与绝热层之间的空腔内填充。不同测点处的压强随时间变化过程如图10所示。可以看出，此时包覆套内外持续存在压差，压差约0.3MPa。这是由于导气条的沟槽横截面积较小，燃气向包覆套与绝热层之间的空腔填充时导气条沟槽内存在压强损失所致。
由于当前的数值模拟未考虑流固耦合效应，固体结构与流体之间的相互作用未能得到充分体现，因此在这种工况下，流固耦合作用对压强分布的影响被忽略。缺乏流固耦合分析使得模拟结果与实际情况及理论分析之间存在一定偏差。在实际工况中，固体结构的变形、热应力以及流体流动的相互作用可能会改变包覆套内外的压强分布，并进一步影响气流的填充过程。因此，未考虑流固耦合效应的数值模拟可能导致对压强分布和流场特性的分析不够精准。引入流固耦合效应后，将能更准确地反映实际工况下流场与结构之间的相互作用，从而提高数值模拟的精度与可靠性。
2.2 考虑流固耦合时发动机内流场及结构场分析
本节采用流固耦合的方法计算发动机内流场与结构件的相互作用，其中固体域材料为可变形体。计算过程中不考虑燃面的退移及包覆套脱落。
2.2.1 包覆套运动情况分析
不考虑燃面退移及包覆套脱落过程时不同测点处的压强随时间变化过程如图12所示。可以看出，测点1~2的压强值接近，测点3~7的压强值接近，根据前文中测点位置可知，测点1的压强数据可代表包覆套内的压强，测点3~7的压强数据可代表包覆套外的压强，由此可知，在0~0.8s范围内，包覆套内外压差先增大后减小；在0.4s时刻压差达到最大值0.604MPa。在压差作用下包覆套和推进剂界面处可能存在额外应力并导致界面的损伤和扩展。
2.2.2 不同缝隙宽度下包覆套受力情况分析
图13 不同缝隙宽度下界面单元的拉应力-时间曲线Fig.13 Tensile stress-time curves of interface elements under different gap widths |
图14 不同缝隙宽度下界面单元的剪应力-时间曲线Fig.14 Shear stress-time curves of interface elements under different gap widths |
由图13可见，随着点火瞬态过程的发展，各工况下界面单元法向应力均表现出由初始拉伸状态逐步向压缩状态转化的特征。这表明界面局部接触状态并非固定不变，而是随着内压加载、包覆套变形及药柱响应的演化，由早期的张开趋势逐步转变为后期的挤压接触状态。进一步比较不同缝隙宽度工况可以发现，随着缝隙宽度增大，界面单元所承受的拉应力整体呈减小趋势。其原因在于，缝隙增大后，包覆套与药柱之间的直接法向约束作用减弱，界面法向载荷传递能力下降，从而削弱了局部“剥离—拉开”效应，导致界面拉应力水平降低。换言之，较小缝隙条件下，界面更易形成直接挤靠与局部约束，包覆套变形更容易在界面起始处转化为法向拉伸作用，因此更易产生局部拉应力集中。
由图14可知，随着缝隙宽度增大，界面单元的剪应力水平逐渐升高。其力学本质在于：缝隙宽度增大后，包覆套局部具有更大的变形空间，界面法向接触约束被削弱，但切向相对位移更容易发展，致使载荷传递方式由法向挤压和拉伸为主，逐步向切向摩擦/剪切传递增强转变。因此，较大缝隙并非单纯改善界面受力状态，而是使界面由“法向拉伸主导”向“剪切作用增强”的受力模式转变，进而提高了界面发生滑移和脱粘的可能性。
由图15中不同缝隙宽度下最后时刻的应变云图可见，最大应变主要集中在包覆套与药柱界面邻近区域，说明该区域是流固耦合作用下的主要变形协调区。随着缝隙宽度增大，应变分布范围有所扩展，局部高应变区的形态也发生变化。这表明缝隙变化不仅影响应力大小，还改变了界面附近的变形传递路径。当缝隙较小时，包覆套与药柱之间相互约束较强，局部变形更容易集中在接触起始区域；而当缝隙增大时，包覆套在压力作用下更容易发生局部偏移和切向滑移，界面变形由局部集中逐渐转向更明显的剪切协调过程。
综合上述分析，缝隙宽度对界面受力模式具有显著影响。随着缝隙宽度增大，界面法向拉应力减小，而剪应力增大，其本质是界面载荷传递路径由法向约束主导逐步向切向滑移主导转变。因此，在结构设计中，缝隙宽度的选取不宜仅依据某一类应力的增减进行判断，而应综合考虑法向剥离风险与切向滑移风险之间的平衡，以实现界面受力状态与损伤倾向的协同优化。
2.2.3 不同弹性模量下包覆套受力情况分析
图16 不同弹性模量下界面单元的拉应力-时间曲线Fig.16 Tensile stress-time curves of interface elements under different elastic moduli |
图17 不同弹性模量下界面单元的剪应力-时间曲线Fig.17 Shear stress-time curves of interface elements under different elastic moduli |
由图16可见，在点火瞬态早期（约0.4s之前），随着包覆套弹性模量减小，界面单元所承受的拉应力整体呈增大趋势。其主要原因在于，较低弹性模量意味着包覆套刚度较小，在相同内压载荷作用下更容易发生径向和轴向变形，从而增大与药柱之间的变形不协调程度。由于药柱与包覆套在界面处存在粘接与接触约束，这种变形不匹配会在界面局部转化为更强的法向拉开作用，因此拉应力升高。随着时间推移，点火初期压差逐步衰减，结构响应趋于稳定，界面法向张开趋势减弱，因此拉应力逐步下降。当弹性模量增大至较高水平时，包覆套整体变形受到抑制，界面法向开裂位移减小，拉应力变化趋于平缓，说明提高包覆套模量有利于减弱界面剥离型受力响应。
由图17可知，不同弹性模量条件下界面单元的剪应力时程存在明显差异。总体上看，弹性模量较小时，包覆套更易通过自身变形释放部分载荷，界面切向约束作用相对较弱，因此剪应力响应变化相对有限；而随着弹性模量增大，包覆套刚度提高，结构变形自由度减小，外部载荷更容易通过界面直接传递，使界面切向约束作用增强，剪应力响应更加明显。这说明弹性模量变化不仅影响界面的法向拉伸水平，也会改变切向载荷在界面处的分配方式。
由图18可见，不同工况下最大应变始终集中于药柱与包覆套界面附近，说明该区域仍是结构响应和损伤敏感的主要区域。随着弹性模量增大，应变集中程度整体减弱，表明包覆套刚度提高后，局部变形被有效抑制，界面附近的协调变形能力增强。相反，在较低弹性模量条件下，包覆套更容易发生局部大变形，使界面处出现更大的相对位移和更明显的应变集中，这将提高界面发生脱粘和局部损伤的风险。因此，从变形协调角度来看，适当提高包覆套弹性模量有助于减弱界面局部高应变集中的不利影响。
综合上述分析，弹性模量通过改变包覆套的结构刚度和界面两侧的变形协调关系，显著影响界面受力模式。当弹性模量减小时，包覆套更易产生较大变形，界面法向拉开作用增强，因而拉应力水平升高；当弹性模量增大时，界面法向剥离趋势得到抑制，但切向载荷传递作用相对增强。由此可见，弹性模量的选取本质上影响的是界面法向开裂倾向与切向滑移倾向之间的分配关系。在工程设计中，应结合界面失效模式和结构服役需求，合理匹配包覆套刚度参数，以改善界面受力状态并降低损伤风险。
3 结论
基于LS-DYNA平台，本文采用流固耦合数值模拟方法，对满装填固体火箭发动机起始工作段包覆套与推进剂界面的损伤风险进行了研究，主要结论如下：
（1）数值仿真与试验获得的压强数据最大相对误差为15.4%，吻合良好，这表明本文建立的流固耦合模型能够较为准确地模拟发动机工作过程中的压强载荷分布。
（2）考虑流固耦合的计算结果表明，包覆套在燃气作用下的动态变形会显著改变流道结构，进而反作用于流场，使内外压差呈现复杂的动态演化规律，这一动态载荷导致界面产生额外应力、进而引发损伤甚至剥离，这与传统刚体假设存在本质差异。
（3）随着缝隙宽度增加，药柱最大剪应力增大而拉应力减小；随着包覆套弹性模量减小，拉应力增大且变化剧烈，剪应力减小且变化平缓。设计时应综合权衡这些参数，以优化界面力学性能，确保结构稳定性和可靠性。
本文的研究深化了对满装填固体发动机界面损伤特性的理解，可为同类发动机的结构优化、防护设计及可靠性评估提供理论依据。后续可进一步考虑燃面退移、包覆套烧蚀脱落及材料粘弹性等复杂因素，以构建更为精确的预测模型。