2022 年 8 期目录

在固体火箭发动机(solid rocket motor, SRM)研制、生产、储存过程中, 会产生如装药界面脱粘、浇注气泡、拔模表面损伤、运输震动裂纹和存储老化裂纹等缺陷, 这些小缺陷对发动机性能的影响复杂:有时没有明显变化, 有时引起大小不同的推力异常, 导致装药结构破坏, 甚至引发爆炸事故。因此, 这些缺陷对结构的安全有重要影响。对激波与裂纹的流固耦合问题已有大量的研究^[1-7]:何国强^[1]针对缺陷对固体发动机性能影响的研究成果进行较细致的综述, C.T.Liu^[2-3]通过实验确定了裂纹增长与应变速率之间的关系, 职世君等^[4]采用数值模拟研究了对流燃烧对裂纹扩展的影响, 韩波等^[5]对激波扫过缺陷时裂纹壁面的变形情况进行了数值模拟和实验研究, 徐春光等^[6]和刘君等^[7]采用数值模拟研究了激波在不同尺度缺陷中的绕射规律。

强激波扫过缺陷时形成的流场极其复杂, 缺陷口存在激波的反射、绕射现象, 缺陷内存在马赫反射并重新生成平面激波的现象, 在缺陷内还存在激波诱导涡运动及激波与漩涡的相互干扰等现象^[6-7]。研究这种激波与漩涡占主导地位的高度非线性非定常流场, 探索其物理机理, 进一步了解激波在缺陷中的传播规律、激波冲击作用下含缺陷装药的动态响应以及应力强度因子的变化规律, 有重要的学术价值和应用价值。本文中, 对激波作用下固体火箭装药缺陷的流固耦合相互作用进行数值模拟。

流固耦合问题的研究方法主要有紧耦合和弱耦合。前者将流体域、固体域和耦合作用构造在同一控制方程中, 在单一时间步内同时求解所有变量; 后者在每一时间步分别采用计算流体动力学(computational fluid dynamics, CFD)和计算固体动力学(computational structural dynamics, CSD)对流体域和固体域依次求解, 通过搭建中间数据交换平台交换固体域和流体域的结果信息, 实现耦合求解。具体的弱耦合包括结构和界面共同求解的弱耦合、流体和界面共同求解的弱耦合和流体和固体都考虑界面的弱耦合。弱耦合方式立足于各自求解方法, 充分发挥各自求解的优势进行学科交叉, 有独特的应用优势^[8-9]。本文中, 采用流体和固体都考虑界面的弱耦合方法分析激波作用下含缺陷固体火箭装药的流固耦合问题。

1. 流体计算数值方法

在任意拉格朗日-欧拉(arbitrary Lagrangian-Eulerian, ALE)有限体积描述下, 二维可压缩非定常Euler方程的积分形式可表示为:

$\frac{\partial}{\partial t} \iint_{\Omega(t)} \boldsymbol{Q} \mathrm{~d} A+\oint_{\partial \Omega(t)} \boldsymbol{F}(\boldsymbol{Q}, \boldsymbol{u}) \cdot \boldsymbol{n} \mathrm{d} l=\boldsymbol{0}$

(1)

式中:Ω(t)为运动的流场区域, ∂Ω(t)为其边界; A为流场区域的面积, l为其边界长度; u为网格运动速度; n为控制体外边界的外法线单位矢量; Q和F分别为守恒变量和包含相对运动速度的对流项:

$\left\{\begin{array}{l} Q=[\rho, \rho u, \rho v, \rho e]^{\mathrm{T}} \\ \boldsymbol{F} \cdot \boldsymbol{n}=(\boldsymbol{U} \cdot \boldsymbol{n})[\rho, \rho u, \rho v, \rho e]^{\mathrm{T}}+p\left[0, n_{x}, n_{y}, u_{\mathrm{n}}\right]^{\mathrm{T}} \end{array}\right.$

(2)

式中:ρ为流体密度; u和v分别为x和y方向的流动速度; e为流体质量内能; p为流体压力; U为流体在ALE坐标系下的流体速度; u_n为网格运动的法向速度; n_x和n_y分别为n的2个分量。分别用u_x和u_y表示u的2个分量, 则:

$\boldsymbol{U}=\left[u-u_{x}, v-u_{y}\right], \quad u_{\mathrm{n}}=u_{x} n_{x}+u_{y} n_{y}$

(3)

采用格心格式的有限体积方法对非结构单元进行空间离散，由格林公式求得单元内的变量梯度，重建单元边界两侧的变量，形成近似的一维Riemann问题，相对通量的分裂采用Van Leer分裂。为了抑制流场中物理量间断处可能出现的数值振荡，采用Venkatakrishman通量限制器。网格单元编号为i的半离散形式为：

$\frac{\partial\left(\boldsymbol{Q}_{i} V_{i}\right)}{\partial t}=-\sum_{k=1}^{N_{\mathrm{f}}} \boldsymbol{F}_{k} \cdot \boldsymbol{n}_{k} S_{k}$

(4)

式中:Q_i为单元i的守恒变量在单元中心处的值, V_i为单元i的体积, F_k、n_k和S_k为单元i第k个表面边界的通量、外法向单位矢量和面积, N_f为控制体的表面个数。为保证时间精度, 流场方程采用四步Runge-Kutta方程进行求解。

2. 固体计算数值方法

在裂纹的数值模拟方面, 有限元法应用最广泛, T.Belytschko等^[10-11]最早将扩展有限元(extended finite element method, XFEM)应用于裂纹扩展问题。XFEM在裂纹扩展的数值模拟中能避免网格划分, 节省求解时间, 因此XFEM成为目前研究的热点。与常规有限元法相比, XFEM不需要重新划分网格, 通过富集自由度描述裂纹及其生长, 大大提高了计算效率。

XFEM中常规有限元的位移近似解为:

$u_{\mathrm{h}}(x)=\sum_{i} N_{i} s_{i}$

(5)

式中:N_i为第i个节点的形函数; s_i为第i个节点的位移。对于弹性问题含裂纹单元的近似位移插值函数可表示为:

$\begin{aligned} & u_{\mathrm{h}}(x)=\sum_{i \in N_{\Gamma}} N_{i}(x)\left[s_{i}+H(x) a_{i}\right]+\sum_{i \in N} N_{i}(x)\left[s_{i}+\sum_{\alpha=1}^{4} \Phi_{\alpha} b_{i}^{(\alpha)}\right]\end{aligned}$

(6)

$\begin{aligned} \left(\Phi_{1}, \Phi_{2}, \Phi_{3}, \Phi_{4}\right)=\sqrt{r}\left(\sin \frac{\theta}{2}, \cos \frac{\theta}{2}, \sin \frac{\theta}{2} \sin \theta, \cos \frac{\theta}{2} \sin \theta\right) \end{aligned}$

式中:a_i为裂纹穿过单元富集节点位移, b_i^(α)(α=1~4)为裂尖所在单元富集节点位移; N_Γ为被裂纹穿过但不包含裂尖的单元节点的集合; N_^为含裂尖的单元节点集合; r为裂尖区域节点到裂尖的距离; θ为裂尖区域节点与裂尖的连线与裂纹的夹角; H(x)为Heaviside函数, 是反映裂纹面位移不连续的附加函数, 用以反映位移不连续情况, 表达式为:

$H(x)=\operatorname{sgn}[\phi(x)]=\left\{\begin{array}{ll} 1 & \phi(x) \gt 0 \\ -1 & \phi(x) \lt 0 \end{array}\right.$

(7)

动力学支配方程为:

$\boldsymbol{M} \ddot{\boldsymbol{d}}(t)+\boldsymbol{C} \dot{\boldsymbol{d}}(t)+\boldsymbol{K} \boldsymbol{d}(t)=\boldsymbol{f}(t) \quad t \in(0, T)$

(8)

式中：M、C和K分别为结构整体质量阵、阻尼阵和刚度阵； $\boldsymbol{d}(t)、 \dot{\boldsymbol{d}}(t)$ 和 $\ddot{\boldsymbol{d}}(t)$ 分别为节点位移列矢量、节点速度矢量和节点加速度矢量，T为时域上限，f(t)为等效节点力矢量。对于常规节点i，位移d_i=[s_i]；对于裂尖所在单元节点i，d_i=[s_i, b_i⁽¹⁾, b_i⁽²⁾, b_i⁽³⁾, b_i⁽⁴⁾]；对于Heaviside函数加强的节点i，d_i=[s_i, a_i]。节点i的节点力矢量为：

$\begin{aligned} \left\{\begin{array}{l}\boldsymbol{f}_{i}=\left(f_{i}^{(s)}, f_{i}^{(a)}, f_{i}^{\left(b_{1}\right)}, f_{i}^{\left(b_{2}\right)}, f_{i}^{\left(b_{3}\right)}, f_{i}^{\left(b_{4}\right)}\right) \end{array}\right. \end{aligned}$

(9)

$\begin{aligned} \left\{\begin{array}{l} f_{i}^{(s)}=\int_{\Omega_{\mathrm{c}}} N_{i} \bar{b} \mathrm{~d} \Omega+\int_{\Gamma_{\mathrm{c}}} N_{i} \bar{t} \mathrm{~d} \Gamma \\ f_{i}^{(a)}=\int_{\Omega_{\mathrm{c}}} N_{i} H \bar{b} \mathrm{~d} \Omega+\int_{\Gamma_{\mathrm{c}}} N_{i} H_{t}^{-} \mathrm{d} \Gamma \\ f_{i}^{\left(b_{\alpha}\right)}=\int_{\Omega_{\mathrm{c}}} N_{i} \Phi_{\alpha} \bar{b} \mathrm{~d} \Omega+\int_{\Gamma_{\mathrm{c}}} N_{i} \Phi_{\alpha} \bar{t} \mathrm{~d} \Gamma \end{array}\right. \end{aligned}$

(10)

式中:b和t分别为体力和面力; Ω_c和Γ_c分别为单元的积分域及其边界; dΓ为边界的积分微元。

对于节点i和节点j, 单元刚度矩阵K_ij^(e)可统一表示为:

$\boldsymbol{K}_{i j}^{(e)}=\left[\begin{array}{lll} \boldsymbol{k}_{i j}^{(s s)} & \boldsymbol{k}_{i j}^{(s a)} & \boldsymbol{k}_{i j}^{(s b)} \\ \boldsymbol{k}_{i j}^{(a s)} & \boldsymbol{k}_{i j}^{(a a)} & \boldsymbol{k}_{i j}^{(a b)} \\ \boldsymbol{k}_{i j}^{(b s)} & \boldsymbol{k}_{i j}^{(b)} & \boldsymbol{k}_{i j}^{(b)} \end{array}\right]$

(11)

式中: $\boldsymbol{k}_{i j}^{(\mu \nu)}=\int_{\Omega}\left(\boldsymbol{B}_{i}^{(\mu)}\right)^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B}_{j}^{(\nu)} \mathrm{d} \Omega(\mu, \nu=s, a, b)$ ; D为弹性本构矩阵; B_i^(s)、B_i^(a)和B_i^(b)为形函数的变分矩阵, $\boldsymbol{B}_{i}^{(b)}=\left[\boldsymbol{B}_{i}^{\left(b_{1}\right)}, \boldsymbol{B}_{i}^{\left(b_{2}\right)}, \boldsymbol{B}_{i}^{\left(b_{3}\right)}, \boldsymbol{B}_{i}^{\left(b_{4}\right)}\right]$ , 具体地

$\boldsymbol{B}_{i}^{(s)}=\left[\begin{array}{cc} N_{i, x} & 0 \\ 0 & N_{i, y} \\ N_{i, y} & N_{i, x} \end{array}\right], \quad \boldsymbol{B}_{i}^{(a)}=\left[\begin{array}{cc} \left(N_{i} H\right)_{, x} & 0 \\ 0 & \left(N_{i} H\right)_{, y} \\ \left(N_{i} H\right)_{, y} & \left(N_{i} H\right)_{, x} \end{array}\right], \quad \boldsymbol{B}_{i}^{\left(b_{\alpha}\right)}=\left[\begin{array}{cc} \left(N_{i} \Phi_{\alpha}\right)_{, x} & 0 \\ 0 & \left(N_{i} \Phi_{a}\right), y \\ \left(N_{i} \Phi_{\alpha}\right)_{, y} & \left(N_{i} \Phi_{a}\right)_{, x} \end{array}\right]$

(12)

式中:下标x和y分别表示对x和y求偏导。

在动态裂纹的数值模拟过程中, 时域处理上通常利用差分法诸如Newmark、Wilson-θ等时域连续的方法对问题进行求解^[10-11]。但在高频强脉动冲击荷载作用下, 传统时程方法难以捕捉波阵面的间断, 求解过程往往带来虚假的数值振荡。针对传统时程算法求解过程的上述问题, X.Li等^[12]和W.Wu等^[13]发展了时域间断Galerkin方法(discontinuous Galerkin finite element method, DGFEM), 即在空间域有限元方法分析基础上, 在时域对位移进行3次Hermite插值, 对位移的导数进行时间的线性插值离散。已有的研究结果表明^[8-12-13], DGFEM具备自动引入数值耗散和滤去虚假的高阶模式和数值振荡效应的能力, 特别对在脉动和冲击荷载作用下结构中波传播过程的数值模拟, 与Newmark算法相比, DGFEM具有良好的捕捉波的间断并消除虚假数值振荡的能力。

考虑动力学方程空间域离散(见(8)式), DGFEM在时域离散时允许基本未知函数φ与其导数ψ在离散点处间断, 即在时刻t_n, 有:

$\Delta \varphi_{n}=\varphi_{n}^{+}-\varphi_{n}^{-}, \quad \Delta \psi_{n}=\psi_{n}^{+}-\psi_{n}^{-}$

(13)

式中:φ_n和ψ_n分别为离散的基本位置函数及其导数, φ_n⁺和φ_n^-分别为φ_n的右极限和左极限, ψ_n⁺和ψ_n^-分别为ψ_n的右极限和左极限, Δφ_n和Δψ_n为间断点处的跃变。DGFEM具体推导分析过程见文献[11-12]。本文中直接给出基本求解公式:

$\left[\begin{array}{cc}\boldsymbol{M}+\frac{\Delta t}{6} \boldsymbol{C}-\frac{(\Delta t)^{2}}{12} \boldsymbol{K} & -\frac{\Delta t}{6} \boldsymbol{C}-\frac{(\Delta t)^{2}}{12} \boldsymbol{K} \\ \frac{\Delta t}{2} \boldsymbol{C}+\frac{(\Delta t)^{2}}{3} \boldsymbol{K} & \boldsymbol{M}+\frac{\Delta t}{2} \boldsymbol{C}+\frac{(\Delta t)^{2}}{6} \boldsymbol{K}\end{array}\right]+\left[\begin{array}{c}\psi_{n} \\ \psi_{n+1}\end{array}\right]=\left[\begin{array}{c}\boldsymbol{Q}_{1}-\boldsymbol{Q}_{2}+\boldsymbol{M} \psi_{n}^{-} \\ \boldsymbol{Q}_{1}+\boldsymbol{Q}_{2}+\boldsymbol{M} \psi_{n}^{-}-\Delta t \boldsymbol{K} \varphi_{n}^{-}\end{array}\right]$

(14)

$\varphi_{n+1}=\varphi_{n}^{-}+\frac{\Delta t}{2}\left(\psi_{n}+\psi_{n+1}\right)$

(15)

式中:Δt为时间步长, φ_n=φ_n⁺=φ_n^-, Q_n为离散守恒变量, 满足:

$\boldsymbol{Q}_{1}=\frac{\Delta t}{3} \boldsymbol{Q}_{n}+\frac{\Delta t}{6} \boldsymbol{Q}_{n+1}, \quad \boldsymbol{Q}_{2}=\frac{\Delta t}{6} \boldsymbol{Q}_{n}+\frac{\Delta t}{3} \boldsymbol{Q}_{n+1}$

(16)

本文中，针对冲击波作用的间断、作用时间短的特点，采用由DGFEM与XFEM结合起来的DG_XFEM对含裂纹缺陷的固体火箭装药问题进行数值模拟。

3. 模型与算例分析

模拟算例网格如图 1所示。初始时刻, 马赫数为2.074 61的正激波入口位于裂纹前2cm处, 出口位于6m处, 如图 2所示。波前气体静止, 气体视为量热完全气体, 比热比γ为1.4。参考速度为1km/s。固体网格划分见图 3, 节点分布为20×20, 实心点为压力监测点所在节点, 星号为固支边界, 空圈为裂纹面所在单元的富集节点, 四方块为裂纹尖端所富集节点。静态强度因子采用平面单边裂纹受单向拉伸问题中的强度因子计算公式(K_I=Fσ(πl_a)^1/2)进行计算, 其中F=1.12-0.23(l_a/l_b)+10.6(l_a/l_b)²-21.7(l_a/l_b)³+30.4(l_a/l_b)⁴, l_a 为裂纹长度, l_b为固体板的宽度, 取σ=10⁶。

图 1 网络模型划分

Figure 1. Model and mesh

下载: 全尺寸图片幻灯片

图 2 初始时刻激波流线

Figure 2. Shock wave at initial time

下载: 全尺寸图片幻灯片

图 3 含裂纹的固体模型及网格划分

Figure 3. Mesh of solid with crack

下载: 全尺寸图片幻灯片

沿固体边界及裂纹面选择38个测压点, 其中测点P₁、P₂₀、P₃₈的坐标分别为(-5cm, 5cm)、(-5cm, -5cm)和(-5cm, -0.26cm)。3个测压点压力曲线如图 4所示。从图 4可以看出, P₂₀处于波后, 压力值不变; P₁在激波经过后, 压力有明显的波动; P₃₈点处于裂纹尖端附近, 激波到达裂纹尖端附近发生多次反射, 经过多次反射, 压力最终达到稳定。不同时刻的流场密度等值线如图 5所示。由图 4~5可以看出, 4 000步(2 869.1μs)后, 裂纹中流场压力趋于稳定。在前5 000步荷载样本中每隔100步取样, 对裂纹的动态响应进行计算, P₁处的位移d随荷载步n变化如图 6所示。分别采用位移反推法和J积分法对裂纹尖端应力强度因子K进行计算, 计算结果如图 7所示, 其中J积分的积分半径1为0.4cm, 积分半径2为1cm。

图 4 监控点压力曲线

Figure 4. Pressure of monitoring points

下载: 全尺寸图片幻灯片

图 5 激波扫裂纹的密度等位线

Figure 5. Density contour of shock wave

下载: 全尺寸图片幻灯片

图 6 P₁处y方向位移曲线

Figure 6. Displacement curve along y axis at P₁

下载: 全尺寸图片幻灯片

图 7 应力强度因子计算比较

Figure 7. Contrast of the stress intensity factors

下载: 全尺寸图片幻灯片

从图 6~7可以看出, 应力强度因子与裂纹张开位移的变化趋势基本一致。采用位移反推法与2种积分半径的J积分计算结果基本吻合, 验证了该方法的可靠性。

4. 结论

基于弱耦合算法实现了流固耦合分析程序框架, 并对冲击波作用下固体火箭装药缺陷动态响应进行数值模拟。从结果可以看出, 激波扫过裂纹的过程中, 在裂纹中经过几次反射, 使得裂纹内压力表现出振荡的现象, 但峰值不断升高最后趋于稳定, 体现了高度非定常非线性的特点; 从应力强度因子的计算看出, 装药缺陷裂纹的应力强度因子变化与冲击波作用下位移响应一致, 验证了该方法的正确性。