一种基于Riemann-SPH的极地环境水下冲击破冰数值模拟方法

本发明涉及爆炸流固耦合数值模拟，具体涉及一种基于riemann-sph的极地环境水下冲击破冰数值模拟方法。

背景技术：

1、水下爆炸破冰过程涉及爆炸冲击波传播、爆轰气泡膨胀与收缩、水体运动、冰层变形破坏以及碎冰运动等多种物理过程，属于典型的强冲击流固耦合问题。在爆炸载荷作用下，冰层不仅受到瞬时冲击压力作用，还受到气泡运动、水体流动以及自由面变化产生的持续载荷影响，因此需要对爆炸载荷传播过程、冰层损伤演化过程以及碎冰运动过程进行连续数值模拟。

2、现有水下爆炸破冰数值模拟方法主要采用有限元方法、任意拉格朗日—欧拉方法（ale）以及传统流固耦合计算方法，通过建立流体区域和冰层区域的离散模型实现爆炸载荷与冰层响应的计算。然而，在强冲击、大变形以及冰层裂纹扩展过程中，上述方法依赖计算网格描述材料运动，容易产生网格畸变、界面穿透以及流固耦合边界处理复杂的问题。同时，在冰层失效后通常采用删除失效单元的方式处理裂纹和碎冰区域，导致失效区域的质量、动量信息丢失，使碎冰无法继续参与后续水体运动和气泡载荷作用过程。

3、无网格粒子法由于不依赖固定网格，适用于大变形和自由界面问题，因此逐渐应用于爆炸冲击及流固耦合数值模拟领域。其中，光滑粒子流体动力学方法（sph）能够通过粒子间相互作用描述材料运动过程。但是，在水—冰等不同介质界面处理中，现有sph方法通常采用压力平均或者人工黏性方式传递界面作用力，未充分考虑冰层内部偏应力对界面法向牵引的影响，使得固体侧载荷难以与流体侧压力在统一通量框架下进行传递。

4、另一方面，冰层属于具有明显非线性变形特征的脆性材料，在爆炸冲击作用下同时存在压缩塑性变形、拉伸开裂以及损伤扩展过程。现有冰层材料模型通常难以在无网格计算框架下同时描述轴对称条件下的周向变形、偏应力演化以及损伤失效后的碎冰状态转换。

5、此外，水下爆炸过程包含爆轰阶段、气泡膨胀阶段以及后续载荷传播阶段，计算过程中存在材料状态快速变化和强非线性响应问题，传统显式时间推进方法容易受到状态突变影响，导致计算稳定性不足。

技术实现思路

1、本发明旨在解决上述技术问题而提供一种基于riemann-sph的极地环境水下冲击破冰数值模拟方法，用于在统一计算模型中描述水、炸药、爆轰产物、空气以及冰层之间的相互作用，实现冰—流体界面的准确载荷传递，并对冰层损伤后的碎冰状态进行连续模拟。

2、本发明的目的通过以下技术方案来实现：

3、一种基于riemann-sph的极地环境水下冲击破冰数值模拟方法，包括：

4、s1、建立轴对称多介质粒子模型，将水、炸药、爆轰产物、空气以及冰层离散为无网格粒子；

5、s2、根据轴对称多介质粒子模型确定相邻粒子之间的粒子对，并针对粒子对建立局部法向黎曼问题；

6、s3、当粒子对包含完整冰粒子时，根据完整冰粒子的总法向牵引构造冰侧等效界面压力，并结合流体侧压力以及两侧声阻抗求解黎曼星区状态，得到粒子通量；

7、s4、基于冰层体积响应与偏应力响应分离的材料模型更新冰层粒子的应力状态，并根据冰层粒子的塑性状态和拉伸状态确定损伤状态；

8、s5、当冰层粒子达到失效状态时，保留冰层粒子的质量、位置、密度以及速度，将冰层粒子的偏应力状态转换为体积响应状态；

9、s6、采用预测—校正时间推进方式更新轴对称多介质粒子模型，并输出冰层粒子的损伤状态、等效塑性应变、等效应力、压力、速度以及位置数据。

10、进一步的，步骤s1包括：

11、建立由径向坐标和轴向坐标构成的轴对称坐标系，将无网格粒子设置为绕对称轴旋转形成的圆环体，并根据圆环体的体积确定无网格粒子的质量；

12、无网格粒子的初始质量满足：

13、

14、其中，表示第个无网格粒子的质量，表示第个无网格粒子的密度，表示第个无网格粒子的径向坐标，表示径向离散长度，表示轴向离散长度；

15、为无网格粒子设置位置、速度、质量、密度、压力、比内能、声速、光滑长度以及材料类型，为冰层粒子设置偏应力、周向偏应力、等效塑性应变以及损伤变量，并将损伤变量的初始值设置为0；

16、根据无网格粒子的质量和密度确定等效环体积，等效环体积满足：

17、

18、其中，表示第个无网格粒子的等效环体积。

19、进一步的，步骤s2包括：

20、采用网格桶对无网格粒子进行邻域搜索，确定支持域内的相邻粒子，并根据相邻粒子构成粒子对；

21、对于材料类型相同的液体对应的粒子对以及由完整冰粒子构成的粒子对，采用带限制器的muscl线性重构得到粒子对中点状态；

22、对于材料类型不同的无网格粒子构成的粒子对、空气对应的粒子对以及爆轰产物对应的粒子对，采用粒子中心状态建立局部法向黎曼问题；

23、根据粒子对中两个无网格粒子的位置确定单位法向，单位法向满足：

24、

25、其中，表示粒子对中第个无网格粒子的位置，表示粒子对中第个无网格粒子的位置；

26、将局部法向黎曼问题两侧的速度投影至单位法向，得到左侧法向速度和右侧法向速度，左侧法向速度和右侧法向速度分别满足：

27、

28、

29、其中，表示局部法向黎曼问题的左侧速度，表示局部法向黎曼问题的右侧速度。

30、进一步的，步骤s3包括：

31、获取粒子对的材料类型和损伤状态，并根据材料类型和损伤状态确定粒子对是否包含完整冰粒子；

32、当粒子对包含完整冰粒子时，根据完整冰粒子的压力、偏应力张量以及单位法向构造冰侧等效界面压力，冰侧等效界面压力满足：

33、

34、其中，表示冰侧等效界面压力，表示完整冰粒子的压力，表示完整冰粒子的偏应力张量，表示单位法向的转置；

35、将冰侧等效界面压力和流体侧压力分别作为局部法向黎曼问题的左右重构压力和；

36、根据局部法向黎曼问题两侧的密度和声速确定两侧声阻抗，两侧声阻抗满足：

37、

38、其中，表示左侧声阻抗，表示右侧声阻抗，表示左侧密度，表示右侧密度，表示左侧声速，表示右侧声速；

39、根据左右重构压力、、左侧声阻抗、右侧声阻抗、左侧法向速度以及右侧法向速度求解星区压力和星区法向速度，其中：

40、

41、

42、根据星区压力和星区法向速度确定黎曼星区状态，并根据黎曼星区状态形成粒子对的密度通量、动量通量以及能量通量，得到粒子通量。

43、进一步的，步骤s4包括：

44、根据冰层粒子的密度确定冰层粒子的体积压力，体积压力满足：

45、

46、根据弹性模量和泊松比确定体积模量，体积模量满足：

47、

48、根据体积模量和参考密度确定冰层粒子用于两侧声阻抗计算的声速，声速满足：

49、

50、其中，表示冰层粒子的压力，表示参考压力，表示冰层粒子的密度，表示参考密度；

51、根据弹性模量和泊松比确定剪切模量，剪切模量满足：

52、。

53、将剪切模量用于更新冰层粒子的偏应力状态。

54、进一步的，步骤s4还包括：

55、根据冰层粒子的速度梯度确定径向应变率、轴向应变率、剪切应变率以及周向应变率，其中周向应变率满足：

56、

57、在冰层粒子接近对称轴的位置，根据径向速度梯度确定周向应变率，周向应变率满足：

58、

59、采用含jaumann客观率的偏应力更新方式更新冰层粒子的偏应力状态，并满足零迹约束：

60、

61、根据径向偏应力、轴向偏应力、周向偏应力以及剪切偏应力确定第二偏应力不变量，第二偏应力不变量满足：

62、

63、根据第二偏应力不变量确定等效应力，等效应力满足：

64、

65、根据初始屈服强度、塑性模量以及等效塑性应变确定屈服强度，屈服强度满足：

66、

67、当等效应力大于屈服强度时，对偏应力状态进行径向回归，并累计等效塑性应变，其中塑性增量满足：

68、

69、进一步的，步骤s4还包括：

70、设置失效等效塑性应变和拉伸截断强度；

71、根据等效塑性应变与失效等效塑性应变之间的关系形成塑性失效判据，塑性失效判据满足：

72、

73、根据拉伸截断强度确定拉伸截断压力，拉伸截断压力满足：

74、

75、根据压力与拉伸截断压力之间的关系形成拉伸失效判据，拉伸失效判据满足：

76、

77、当塑性失效判据成立时，将损伤变量设置为1；

78、当拉伸失效判据成立时，将损伤变量设置为1；

79、根据损伤变量确定损伤状态，损伤变量为1时确定冰层粒子达到失效状态。

80、进一步的，步骤s5包括：

81、当损伤变量表征冰层粒子达到失效状态时，将冰层粒子的径向偏应力、轴向偏应力、周向偏应力、剪切偏应力以及偏应力做功清零；

82、保留达到失效状态的冰层粒子的质量、位置、密度、压力以及速度；

83、根据冰层体积响应与偏应力响应分离的材料模型，保留达到失效状态的冰层粒子的体积响应，并根据体积响应将冰层粒子的偏应力状态转换为体积响应状态；

84、将达到失效状态的冰层粒子继续作为无网格粒子参与邻域搜索、局部法向黎曼问题求解以及粒子通量计算。

85、进一步的，步骤s6包括：

86、保存当前时刻的轴对称多介质粒子模型；

87、基于当前时刻的轴对称多介质粒子模型依次执行状态方程计算、爆轰转换、镜像边界处理、邻域搜索、针对粒子对的局部法向黎曼问题求解、冰层粒子的应力状态和损伤状态更新以及轴对称几何源项计算，并按照半个时间步得到预测状态；

88、爆轰转换包括：将炸药设置为未反应炸药状态，当起爆点至对应无网格粒子的距离不大于爆速与起爆后时间的乘积时，将对应无网格粒子的材料状态转换为爆轰产物；当爆轰产物的jwl状态方程指数项的绝对值低于预设阈值时，将爆轰产物转换为理想气体支路；

89、基于预测状态再次执行状态方程计算、爆轰转换、镜像边界处理、邻域搜索、针对粒子对的局部法向黎曼问题求解、冰层粒子的应力状态和损伤状态更新以及轴对称几何源项计算，得到预测状态变化率；

90、根据当前时刻保存的轴对称多介质粒子模型和一个时间步对应的预测状态变化率得到下一时刻的轴对称多介质粒子模型。

91、进一步的，步骤s6还包括：

92、根据声速—速度约束、加速度约束以及密度下降率约束确定时间步，时间步满足：

93、

94、其中，表示声速—速度约束系数，表示加速度约束系数，表示密度下降率约束系数，表示最小光滑长度，表示最大声速，表示最大速度，表示最大加速度，表示无网格粒子的密度，表示无网格粒子的密度变化率；

95、当无网格粒子的密度下降时，采用密度下降率约束参与确定时间步；

96、在完成一个时间步的预测—校正时间推进方式后，对水、空气以及爆轰产物对应的无网格粒子执行粒子位置修正，保持完整冰粒子的位置不执行粒子位置修正；

97、根据更新后的轴对称多介质粒子模型输出冰层粒子的损伤状态、等效塑性应变、等效应力、压力、速度以及位置数据。

98、与现有技术相比，本发明具有以下有益效果：

99、1、本发明通过建立包含水、炸药、爆轰产物、空气以及冰层的轴对称多介质粒子模型，并基于riemann-sph方法对不同介质之间的相互作用过程进行计算；针对冰层与流体之间的界面作用，根据冰层粒子的总法向牵引构造冰侧等效界面压力，并结合流体侧压力求解黎曼星区状态，使冰层载荷状态参与粒子通量计算。

100、2、本发明通过建立体积响应与偏应力响应分离的冰层材料模型，基于包含周向分量的j2弹塑性模型更新冰层粒子的应力状态，并根据塑性状态和拉伸状态确定冰层粒子的损伤状态；在冰层粒子达到失效状态后，保留冰层粒子的质量、位置、密度以及速度，并将冰层粒子的偏应力状态转换为体积响应状态，使失效冰粒子继续参与后续粒子计算过程。

101、3、本发明采用爆炸状态转换和预测—校正时间推进方式，对炸药爆轰过程、爆轰产物膨胀过程以及后续气泡运动过程进行连续计算，并根据更新后的粒子状态输出冰层损伤状态、等效塑性应变、等效应力、压力、速度以及位置数据，以获得破冰过程中的损伤区域和碎冰区域信息。