1. 引言

随着人类科技水平的提高，太空探索逐渐成为未来各国新的发展方向，目前航天多以化学能、太阳能为主要能源。但化学能受质量和体积的限制、太阳能受光照影响，均无法满足太空探索的需要。空间堆具有寿命长、不依赖太阳光照、能量密度高等特点。因此，核能在航天领域有着广阔的应用前景[1]。

然而受到尺寸及重量的限制，空间堆无法做到与地面核设施同等的冗余性和多样性。同时空间堆一旦发生事故将影响更为恶劣。截至目前，已经开展的空间核动力源应用中共发生了9次核事故。其中最为严重的是1978年苏联的“宇宙-954号”(SPACE-954)坠落在加拿大，引起了世界各国对空间核动力应用安全的高度关注[2]。海洋占据了地球表面约71%的面积，当反应堆于发射后期或在轨期阶段坠落后有很大可能落入海洋。因此，对空间堆撞击水体问题的研究具有重要意义。

目前国内外对空间堆撞击水体研究开展较少，由于空间堆跌落实验不现实，研究主要从理论研究及仿真模拟两个方面开展。理论研究集中在计算空间堆在各工况下的k_eff，安伟健等[3]针对Kilopower空间堆研究了其在堆芯被水体侵入等工况下的k_eff，发现在水体侵入时空间堆可能达到瞬发超临界状态。肖启冬等[4]针对TOPAZ-II空间堆(Thermionic Operational Reactor in Power Assembly-II)进行了特殊临界安全分析，基于实验数据给出了TOPAZ-II空间堆掉落后的最危险事故情景。仿真模拟方面，国外Kim等[5]建立了空间反应堆高速撞击地面的流体力学模型，对反应堆以230 m/s和150 m/s的速度撞击干砂的工况进行了计算并得到了相应的破坏结果。

理论分析只能分析简单假设下的空间堆反应性，无法包络空间堆撞击水体后的各种破坏结果，仿真模拟则暂无对空间堆撞击水体的相关研究，而SPH法已广泛应用于解决入水问题，丁宁等[6]基于SPH法对返回舱的入水过程进行模拟。结果表明，SPH法能够较好地模拟返回舱入水过程。高英杰等[7]采用SPH法进行100 m/s初速度条件下回转体垂直入水数值计算并对方法的有效性进行验证。证实SPH法可以有效解决入水冲击流固耦合数值仿真中的计算适应性问题。综合上述分析，目前已有研究未涉及空间堆撞击水体问题，因此本文通过SPH法与欧拉法，以TOPAZ-II型空间堆结构为参考研究空间堆高速撞击水体问题。

2. 计算模型

2.1. 几何模型及网格划分

1975年，前苏联开发基于热离子系统的TOPAZ-II型空间核反应堆电源[8]，该空间堆以UO₂为燃料，NaK液态金属为冷却剂。TOPAZ-II堆的堆芯由热离子燃料元件、氢化锆慢化剂、端部铍反射层、侧铍反射层、12个控制鼓以及堆容器等组成，本文反应堆以TOPAZ-II堆为基础进行适当简化，其横截面和剖视图如图1所示。

Figure 1. Space reactor sectional view and cross section

图1. 空间堆横截面和剖视图

空间堆由多个部件组成，其主体结构材料为2A12铝合金，反射层、控制鼓筒体材料为不锈钢，填充为铍材，安全棒材料为B₄C，燃料元件内层材料为UO₂，外层材料为不锈钢。

堆本体有限元模型如图2所示。堆芯筒体、安全棒套管及反射层、控制鼓外壳采用壳单元，单元类型为S4R。其余部件均为实体单元，单元类型为C3D8R。反应堆有限元模型单元总数为475,636，结点总数为698,944。

Figure 2. Finite element model of space reactor

图2. 空间堆有限元模型

燃料元件简化为内外两层及上下封盖，外层材料为不锈钢；内层材料为UO₂。反应堆筒体与栅板及上下封盖均绑定连接。燃料元件外层与上下栅板绑定连接。连接螺栓分别将反射层、控制鼓与上下栅板连接，螺栓不施加预紧力。

2.2. SPH方法基本理论

SPH方法基于粒子，求解问题的本质就是选择合适的光滑核函数(插值函数)进行插值的过程．对于任意未知函数，用下式来近似x处的函数值：

(1)

式中： $f\left(x\right)$ 为三维坐标向量 $x$ 的函数； $\Omega$ 为包含 $x$ 的积分域； $\delta \left(x-x^{\prime }\right)$ 为狄拉克 $\delta$ 函数，具有如下性质：

$\delta \left(x-x^{\prime }\right)=\{\begin{array}{c}\infty ,x=x^{\prime }\\ 0,x\ne x^{\prime }\end{array}$ (2)

通过将 $\delta \left(x-x^{\prime }\right)$ 替换成核函数 $W\left(x-x^{\prime },h\right)$ 可得到场函数的核近似表达式

(3)

式中： $x$ 为计算点处的位置向量， $x^{\prime }$ 为问题域内其他点的位置向量， $\text{d}{x}^{\prime }$ 为 $x^{\prime }$ 点处微元体的体积，h为光滑长度， $\Omega$ 为计算点处核函数的支持域。粒子近似法是粒子叠加求和相对应的离散化过程，对粒子的求解与光滑函数有关，通过其支持域内的粒子表达式来进行求解，在粒子处的函数粒子近似式可用下式表达：

(4)

式中：N为搜索域内的粒子总数， $m_j$ 和 ${\rho }_j$ 分别表示对粒子i作用的粒子j的质量和密度， $W_{ij}$ 表示粒子j对粒子i产生影响的光滑函数。对式4求导得：

(5)

用SPH方法离散后得到的控制方程表示如下：

(6)

2.3. 工况设计及假设

空间堆撞击水体时的速度、角度受多种因素影响，包括事故发生时的发射阶段、初速度及再入大气层过程中的烧蚀作用等，美国对同位素电池(RHU)开展过相关跌落实验，4种试件的平均碰撞速度分别为：93.9、95.1、123.4和81.7 m/s [9]，而发射前空间堆自火箭端跌落时可视为自由落地，计算得初速度在25~40 m/s区间内，同时当撞击速度高于200 m/s时根据以往火箭发动机及一级、二级火箭掉落经验，认为堆体完全失效，即不存在临界风险。

因此，根据以上对空间堆发射各阶段可能发生事故情况的分析，将撞击地面初速度设定为25 m/s (可能发生撞击的最低速度)、200 m/s (对应极限工况)。反应堆撞击地面角度设定为0˚、45˚、90˚，分别对应反应堆正面撞击、斜45˚撞击以及侧面撞击。

2.4. 本构模型及失效模型

材料本构模型及失效准则对空间堆结构破坏有着很大影响，本文中使用Johnson-Cook本构模型、失效模型及JH-2模型对材料进行定义。Johnson-Cook本构模型于上个世纪八十年代提出，主要用于描述材料在高应变率、大应变以及不同环境温度下的力学性能。该模型通过大量实验提出，广泛应用于鸟撞击实验、汽车碰撞、霍普金森杆等冲击领域。Johnson-Cook模型的特点是形式简单，同时考虑了应变硬化、应变率强化以及温度软化效应的影响，Johnson-Cook本构模型可表示为：

$\bar{\sigma }=\left[A+B{\left({\bar{\epsilon }}^{pl}\right)}^n\right]\left[1+C\ln {\stackrel{\dot{}}{\epsilon }}^{\ast }\right]\left[1-{\left(\frac{T-T_r}{T_{\text{melt}}-T_r}\right)}^m\right]$ (7)

其中，A为参考应变率和参考温度下的初始屈服应力；B，n为材料应变硬化模量和硬化指数；C为材料应变率强化参数； ${\bar{\epsilon }}^{pl}$ 为等效塑性应变； ${\stackrel{\dot{}}{\epsilon }}^{*}$ 为无量纲化等效塑性应变率； $T_r$ 为室温； $T_{\text{melt}}$ 为融化温度；M为温度软化指数，本文中不考虑该效应。

Johnson-Cook失效准则考虑到应力三轴度、温度以及应变率效应：

${\epsilon }_{\text{f}}=\left[D_1+D_2\text{exp}\left(D_3{\sigma }^{\ast }\right)\right]\left(1+D_4\text{ln}{\stackrel{\dot{}}{\epsilon }}^{\ast }\right)\left(1+D_5{T}^{\ast }\right)$ (8)

其中， $D_1~D_5$ 为材料参数； ${\sigma }^{\text{*}}$ 为应力三轴度， ${\stackrel{\dot{}}{\epsilon }}^{\text{*}}$为等效塑性应变率； $T^{\text{*}}$ 为无量纲温度参数。J-C失效准则假设，损伤变量初始为0，材料在变量积累达到1时失效。损伤变量定义为 $D={\displaystyle \sum \left(\Delta {\epsilon }_{eq}/{\epsilon }_f\right)}$ ， $\text{Δ}{\epsilon }_{\text{eq}}$ 为一个时间步的等效塑性应变增量。

本文中堆容器主体结构材料2A12铝合金[10]、不锈钢外壳[11]及螺栓材料[12]使用Johnson-Cook本构模型、失效模型，Be本构方程为Johnson-Cook模型，失效模型为最大主应力失效模型，失效应力为460 Mpa [13]，以上材料输入参数均来自相关文献中对对应材料的实验及参数拟合[10]-[13]，数值列于表1。

Table 1. Materials and functions of each component

表1. 各组成部件材料输入参数

B₄C安全棒材料陶瓷具有较高的热中子吸收能力，常用于反应堆的控制棒材料，其本构及失效模型采用JH-2本构模型。

JH-2理论模型是一种专门用于表征脆性材料断裂行为的本构模型。该模型通过强度模型、损伤演化模型以及静水压力–体积应变状态方程三个核心组成部分，系统描述了脆性材料在受力过程中的弹性变形阶段、塑性软化阶段以及最终的脆性断裂阶段。

强度模型表示如下：

$\{\begin{array}{c}{\sigma }^{\ast }={\sigma }_i^{\ast }-D\left({\sigma }_i^{\ast }-{\sigma }_f^{\ast }\right)\\ {\sigma }_i^{\ast }=A{\left(P^{\ast }+T^{\ast }\right)}^N\left(1+C\ln \frac{\stackrel{\dot{}}{\epsilon }}{{\stackrel{\dot{}}{\epsilon }}_0}\right)\\ {\sigma }_f^{\ast }=B{\left(P^{\ast }\right)}^M\left(1+C\ln \frac{\stackrel{\dot{}}{\epsilon }}{{\stackrel{\dot{}}{\epsilon }}_0}\right)\end{array}$ (9)

其中： ${\sigma }^{\ast }$ 为材料当前标准化等效应力； ${\sigma }_i^{\ast }$ 为材料初始状态下标准化完整等效应力； ${\sigma }_f^{\ast }$ 为材料破碎后标准化断裂等效应力；D为损伤因子(0 ≤ D ≤ 1)；A，N，B，M，C为材料常数，可通过材料断裂等实验获得，这里取值分别为0.927，0.67，0.7，0.85，0.005。 $P^{*}$ 为标准化静水压力， $P^{*}=P/P_{\text{HEL}}$ ，P为材料实际静水压力， $P_{\text{HEL}}$ 为HEL处压力分量，取值为8.71 GPa； $T^{*}$ 为最大标准化净拉伸强度，取值为0.2； $\stackrel{\dot{}}{\epsilon }$ 为实际应变率， ${\stackrel{\dot{}}{\epsilon }}_0$ 为参考应变率，取1 s⁻¹ [14]。

B₄C模型其他参数如下，密度为2510 kg/m³，剪切模量为197 GPa，最大拉伸强度T为0.26 GPa，Hugoniot弹性极限HEL为8.71 GPa [14]。

水体使用us-up状态方程定义，参数如下，密度为1000 kg/m³，c0 = 1450，s、Gamma0均为0，动力粘度为0.001。

3. 计算结果分析

欧拉法利用abaqus软件进行计算，SPH法利用ansys/dyna软件进行计算。二者区别在于水体的网格，欧拉法需要建立水体域与空气域，计算网格固定在空间中，不随物体运动，而材料相对于网格运动。SPH法中，水体网格被转化为粒子从而避免了严重的单元变形问题。两种方法使用的空间堆有限元模型相同。

3.1. 欧拉法撞击水体过程

欧拉法水体模型高6 m，长宽为5 m × 5 m，水体及空气域各高3 m，有限元模型如图3所示。空间堆初始条件下位于空气域水平方向的正中位置，距液气分界线5 mm。空气域与水体域网格尺寸中心加密以提高计算的精确度。水体边界无约束。

Figure 3. Eulerian mesh for water domain

图3. 欧拉法水域网格

图4所示为欧拉法空间堆200 m/s，0˚撞击水体过程，安全棒首先接触水体，由于安全棒为脆性材料，接触后即破碎，在0.4 ms内下裙板及螺栓在撞击水体后失效，同时填充铍材的外层不锈钢壳也失去包裹作用；2.4 ms内空间堆持续下沉，反射层及安全鼓失去下部固定约束向外运动；4 ms时空间堆保持向下运动趋势但破坏效应基本结束，内部燃料元件散落，无临界风险。

Figure 4. Stress contours of space reactor impacting water at 200 m/s and 0˚ using the eulerian method

图4. 空间堆200 m/s，0˚欧拉法撞击水体应力云图

图5所示为45˚下空间堆以200 m/s撞击水体过程。0.8 ms内安全棒下端首先接触水体并从中间折断，下裙板及螺栓接触水面随即发生形变并破坏，同时空间堆在水体接触力作用下发生转动，空间堆与水面夹角减小，上栅板接触水面并轻微破坏；2.4 ms内空间堆下栅板与筒体部分破坏，燃料元件接触水体并变形；4 ms空间堆完全进入水体，动能大部分被水体吸收，空间堆破坏效应基本完成，近地侧反射层与控制鼓完全失去螺栓约束开始与空间堆分离。

Figure 5. Stress contours of space reactor impacting water at 200 m/s and 45˚ using the eulerian method

图5. 空间堆200 m/s，45˚欧拉法撞击水体应力云图

图6为90˚下空间堆撞击水体过程，0.4 ms内空间堆侧面先撞击水体，上下裙板撞击水面后弯曲变形，与筒体连接处失效；2.4 ms内空间堆完全进入水体，部分螺栓破坏失效，但仍随着空间堆向下运动；4 ms内三根安全棒中近地端安全棒从中部折断，其余两根部分破坏。

Figure 6. Stress contours of space reactor impacting water at 200 m/s and 90˚ using the eulerian method

图6. 空间堆200 m/s，90˚欧拉法撞击水体应力云图

50 m/s下三个角度对应的空间堆损伤情况都很轻微，各角度对应空间堆损伤情况如图7所示。空间堆筒体、反射层、燃料元件都未出现失效单元，损伤集中在安全棒、螺栓及下裙板。0˚时下裙板变形，裙板边缘破坏，安全棒最外端撞击水体后破坏但未折断；45˚时裙板撞击水体侧破坏，安全棒下端断裂，下层近水侧螺栓失效但反射层受上端螺栓约束未脱离空间堆；90˚时由于空间堆侧面首先与水体接触且侧面的反射层强度较高，故仅有个别螺栓及反射层不锈钢外壳出现破坏，其余部件完好。

Figure 7. Damage contours of space reactor impacting water at 50 m/s with different angles

图7. 50 m/s下不同角度空间堆撞击水体损伤云图

3.2. SPH法撞击水体过程

SPH法将整个水体转化为粒子，为了提高计算效率本文将水体分为内层及外层两个区域，内层尺寸为1 m × 1 m × 3 m，内层粒子位于空间堆撞击方向上，尺寸加密为20 mm，外层粒子尺寸设置为65 mm，水体边界无约束。图8为SPH法水体有限元模型。

Figure 8. Finite element model of water impact using SPH method

图8. SPH法撞击水体有限元模型

图9为空间堆以200 m/s各角度撞击水体过程。垂直撞击时空间堆运动过程与欧拉法相似，但结构破坏情况有很大不同，在2.4 ms后空间堆完全进入水体，欧拉法模拟中空间堆上封盖未受到水体冲击保持完好，而SPH法中空间堆上封盖受到水体粒子冲击破坏。45˚撞击水体模拟中两种方法过程及破坏情况相似，4 ms后空间堆下裙板及筒体失效，水体进入空间堆。水平撞击水体过程，当空间堆与液面接触后，空间堆周围附近液面有所上升。随着空间堆进入水体继续向下运动，空泡开始形成并扩张，空泡长度逐渐增长。4 ms内上封盖及下裙板破坏，残余部分飞出，燃料元件变形，反射层、控制鼓未脱离。

Figure 9. Stress contours of space reactor impacting water at 200 m/s with different angles using the SPH method

图9. 空间堆SPH法200 m/s各角度撞击水体应力云图

图10为空间堆以50 m/s撞击水体结果。图10中可以看出空间堆在此速度下撞击水体时破坏非常轻微且45˚倾角撞击时的整体应力相对更高。与欧拉法对比，0˚时裙板均出现了较大变形，45˚时安全棒下端均发生了断裂，90˚时空间堆破坏均较为轻微。因此在此速度下两种方法的模拟结果相近。

Figure 10. Stress contours of space reactor impacting water at 50 m/s and 1 ms with different angles

图10. 空间堆以50 m/s各角度下撞击水体1 ms时应力云图

3.3. 两种方法对比

从本节计算结果中可以看出，低速撞击时欧拉法与SPH法计算结果区别不明显而高速撞击时两者对应结果有着较大差异，需要确认哪种方法更适用于解决本问题。本文选用空间堆模型外形为圆柱形，因此参考有关圆柱体入水文献中的相关研究，本文两种方法的计算结果与文献进行比较如下。

Figure 11. Comparison of the process of a cylinder entering water [15] [16]

图11. 圆柱体入水过程对比[15] [16]

如图11为空间堆撞击水体过程对比，第一行为SPH法，第二行为文献中的撞击过程，第三行为欧拉法撞击过程。图11(a)为圆柱体低速水平撞击水体过程，图11(b)为圆柱体100 m/s垂直入水过程，观察液面空泡轮廓及溅起水体形状可以发现，SPH法仿真结果与实验结果相近，欧拉法受限于网格无法准确模拟空间堆入水过程。

表2中列出了0˚下两种仿真方法对应沙漏能及计算时间，通过对比可以发现SPH法具有计算效率高、计算过程产生沙漏能低的优点。综合以上比较，认为选用SPH法模拟水体更为适合。

Table 2. Comparison of hourglass energy and computation time corresponding to two simulation methods at 0˚

表2. 0˚下两种仿真方法对应沙漏能及计算时间对比

Figure 12. The energy variation curve of space reactor impacting water body horizontally at 200 m/s

图12. 空间堆以200 m/s水平撞击水体能量变化曲线

图12为200 m/s水平撞击时的能量变化曲线，可以看出沙漏能占动能比例远低于5%，即沙漏能对仿真结果影响不大且空间堆撞击水体4 ms时侵蚀能及内能基本稳定，即破坏效应基本完成。

4. 结论

本章在abaqus软件及ansys/dyna软件中使用欧拉法及SPH法两种方法分别在0˚、45˚、90˚三个角度及50 m/s、200 m/s两个速度下对空间堆撞击水体过程进行了模拟仿真。对比仿真结果并参考相关实验及文献中圆柱体入水过程认为SPH法相对欧拉法更适用于解决空间堆撞击水体问题，因此选用SPH法对运动参数在空间堆撞击水体过程中的影响进行了分析，得到以下结论：

欧拉法及SPH法对于接触面(空间堆底部)的模拟都能做到较为准确，但欧拉法对于不规则、自由、复杂边界问题受限于水体网格无法准确模拟。

SPH法相对于欧拉法有着更高的计算效率与更低的沙漏能，计算结果更为贴合实际，因此SPH法更适于解决空间堆撞击水体问题。

随着撞击速度的提高，空间堆破坏程度逐渐加深。低速撞击水体时反应堆基本保持完好；高速撞击水体时，45˚撞击时结构破坏最严重，0˚时次之，90˚时结构破坏最轻。

参考文献