基于虚拟流体方法的深海耐压结构内爆计算方法及系统

本发明涉及流固耦合，涉及高压水下环境固体结构的模拟和动态响应的计算，具体地，涉及一种基于虚拟流体方法的深海耐压结构内爆计算方法及系统。

背景技术：

1、深海耐压结构液-固-气三相流固耦合的水下内爆，是流固耦合问题中的一种，常见于水下压力容器的测试、深海潜水器和潜艇的研发。对深海耐压结构液-固-气三相流固耦合的水下内爆问题进行数值模拟，可以预测深海耐压结构的动态响应和变形以有材料的各种，从而进一步为深海载人设备的设计和研发提供参考。

2、用于数值模拟深海耐压结构液-固-气三相流固耦合的现有技术有存在以下缺陷：流固边界处需要加密网格，对于流固界面的物理量，仅考虑使用当前时刻的值，不考虑使用下一时刻的值。

技术实现思路

1、针对现有技术中的缺陷，本发明提供一种基于虚拟流体方法的深海耐压结构内爆计算方法及系统。

2、根据本发明提供的一种基于虚拟流体方法的深海耐压结构内爆计算方法及系统，所述方案如下：

3、第一方面，提供了一种基于虚拟流体方法的深海耐压结构内爆计算方法，所述方法包括：

4、步骤s1：设定流场初始条件，得到深海耐压结构的计算环境；

5、步骤s2：设定固体初始条件，得到深海耐压结构的几何模型；

6、步骤s3：用步骤s1和s2中的流场初始条件以及固体初始条件，或者用上一时间步的已知结果，结合虚拟流体黎曼解算器，迭代求解流固边界，模拟出深海环境水和结构体的相互作用；

7、步骤s4：用虚拟流体方法对流场虚拟网格赋值，使整个流场的物理量连续；

8、步骤s5：用步骤s1中的流场初始条件，或者用上一时间步的已知结果，结合tvd格式，对流场求解，刻画波在流场中的传播过程；

9、步骤s6：用步骤s2中固体初始条件，或者使用上一时间步的已知结果，结合有限元法，对固体网格结点求解，得到结构体的动态响应和变形过程。

10、优选地，所述步骤s1包括：

11、深海耐压结构，所处的流场初始条件为：

12、pw＝ρwgh+p0                         1)

13、其中，ρw为海水密度，g为重力加速度，h；深海耐压结构所处的水深，p0为海平面大气压。

14、优选地，所述步骤s2包括：

15、设定固体初始条件：密度ρs，杨氏模量ys,泊松比γs，屈服应力sy0,剪切模量gs；

16、给出得到深海耐压结构的几何模型，耐压结构内部的空气物理状态参数满足：

17、pa＝(γa-1)ρaea                         2)

18、其中，pa为气体压强，γa为气体绝热指数，ρa为气体密度，ea为气体内能。

19、优选地，所述步骤s3包括：

20、步骤s3.1：在流固界面处，固体获得流体的压强，ps＝pl，ps为固体压强，pl为流体压强；代入固体求解器式3)～式16)求解新的固体物理状态；

21、采用八节点六面体有限单元，插值函数选择三维等参插值函数：

22、

23、其中，a，b和c分别为待定系数，不同的物理量对应不同的系数；对于任一物理量q(x,y,z)，采用八个单元节点的插值形式：

24、

25、其中，qi为第i个节点物理量；加速度ak、位移xk及其导数的变分表示为：

26、

27、

28、

29、其中，δ为变分符号，下标i表示爱因斯坦求和约定，下标k和上标k均表示第k个单元，m表示沿x轴、y轴或z轴的分量。加速度ak、位移xk及其导数添加下标i后,表示第k个单元在第i个节点的加速度；表示第k个单元在第i个节点的位移；表示第k个单元在m轴上速度；将整个区域划分为ne个单元，

30、

31、其中，f为待求物理量，vn为单元体积，sn为单元表面积，[a]是由插值函数组成的矩阵，其形式为，

32、

33、引入偏导数矩阵[b]，

34、

35、和应力向量[t]，

36、[t]＝(σ11,σ22,σ33,σ23,σ31,σ12)t11)

37、其中，σ11,σ22,σ33,σ23,σ31,σ12分别表示不同方向的应力，下标1对应x轴，下标2对应y轴，下标3对应z轴；对于时间的推进，采用显示时间差分，

38、vn+1＝vn+an(δtn+δtn+1)/2                  12)

39、(δr)n+1＝(δr)n+vn+1·δtn+1                  13)

40、其中，δtn和δtn+1分别为n和n+1时刻的时间步长，δr和v分别为位移和速度。

41、用沙漏矢量抑制沙漏的产生，

42、

43、

44、

45、其中，下标k表示网格节点，i表示x轴、y轴或z轴的方向，α表示沙漏的各个方向，a＝qρl2c/4，q为常数，c是固体声速；l是单元特征长度；γαk是沙漏基矢量，γαk为人工粘性系数，为人工粘性项，nk,i为形函数，giα为沙漏质量，表示速度。

46、步骤s3.2：将新的固体加速度as代入虚拟流体黎曼解算器式14)～16)，求解新的流体压强pl；

47、沿着右特征线

48、

49、的相容关系为：

50、

51、将界面两侧的ρc视为常数，

52、

53、其中，ui、pi、ρil分别为界面流体速度、界面流体压力、界面流体密度，ul为流体相对于流固界面的速度；cil表示界面流体声速；pl表示界面附近流体压强；ρl表示表示界面附近流体密度；cl表示表示界面附近流体声速。

54、步骤s3.3：回到步骤s3.1，直到误差小于设定值。

55、优选地，所述步骤s4包括：

56、定义为符号距离函数，

57、

58、其中，ω+表示虚拟区域，ω-表示真实区域，表示物质界面，dis表示距离，表示位置矢量。

59、将步骤s3得到的流固界面压强和流固界面速度，用下述方程外推：

60、

61、其中，q为需要外推的函数；为流固界面法向单位向量，从流体指向固体；t为人工时间；把i从ω-外推到ω+时，方程18)取+；把i从ω+外推到ω-时，方程18)取-。

62、优选地，所述步骤s5包括：

63、使用二阶tvd格式，设待求网格节点为(i0,j,k)，流场三个方向的通量f，g，h都需要离散，

64、f＝(ρu,ρu2+p,ρuv,ρuw,u(e+p))t             22)

65、g＝(ρv,ρvu,ρv2+p,ρvw,v(e+p))t             23)

66、h＝(ρu,ρwu,ρwv,ρw2+p,w(e+p))t             24)

67、以f为例，

68、

69、其中，和表示半节点的值；ρu、ρσ、ρu、ρuv、ρuw、ρvw均为与速度相关的动量，p为压强，e为单位体积流体的能量。

70、

71、这里，r为右特征矩阵；看作一个限制器，求解方程组如下：

72、

73、

74、

75、

76、

77、

78、

79、

80、这里，ε为一小量，λ为特征值，l为左特征矩阵，αm为特征向量，z为某一特征值，ψ，q和σ均为限制函数，δ为平均速度；

81、求解和时，需要用到半节点的值；采用roe平均求解下列物理量：

82、

83、

84、

85、

86、

87、其中，为平均值，

88、

89、和分别为半节点网格的密度、三个方向的速度和焓，接下来可间接求得能量和声速和

90、同理，对于任一网格节点(i,j,k)，依次求得

91、对于时间步长δt的选择，取

92、

93、cfl为收敛系数，cfl≤1。

94、优选地，所述步骤s6中使用公式3)至16)求解结构体动态响应曲线和结构体的变形。

95、第二方面，提供了一种基于虚拟流体方法的深海耐压结构内爆计算系统，所述系统包括：

96、模块m1：设定流场初始条件，得到深海耐压结构的计算环境；

97、模块m2：设定固体初始条件，得到深海耐压结构的几何模型；

98、模块m3：用模块m1和m2中的流场初始条件以及固体初始条件，或者用上一时间步的已知结果，结合虚拟流体黎曼解算器，迭代求解流固边界，模拟出深海环境水和结构体的相互作用；

99、模块m4：用虚拟流体方法对流场虚拟网格赋值，使整个流场的物理量连续；

100、模块m5：用模块m1中的流场初始条件，或者用上一时间步的已知结果，结合tvd格式，对流场求解，刻画波在流场中的传播过程；

101、模块m6：用模块m2中固体初始条件，或者使用上一时间步的已知结果，结合有限元法，对固体网格结点求解，得到结构体的动态响应和变形过程。

102、优选地，所述模块m1包括：

103、深海耐压结构，所处的流场初始条件为：

104、pw＝ρwgh+p0                         1)

105、其中，ρw为海水密度，g为重力加速度，h；深海耐压结构所处的水深，p0为海平面大气压；

106、所述模块m2包括：

107、设定固体初始条件：密度ρs，杨氏模量ys,泊松比γs，屈服应力sy0,剪切模量gs；

108、给出得到深海耐压结构的几何模型，耐压结构内部的空气物理状态参数满足：

109、pa＝(γa-1)ρaea                         2)

110、其中，pa为气体压强，γa为气体绝热指数，ρa为气体密度，ea为气体内能；

111、所述模块m3包括：

112、模块m3.1：在流固界面处，固体获得流体的压强，ps＝pl，ps为固体压强，pl为流体压强；代入固体求解器式3)～式16)求解新的固体物理状态；

113、采用八节点六面体有限单元，插值函数选择三维等参插值函数：

114、

115、其中，a，b和c分别为待定系数，不同的物理量对应不同的系数；对于任一物理量q(x,y,z)，采用八个单元节点的插值形式：

116、

117、其中，qi为第i个节点物理量；加速度ak、位移xk及其导数的变分表示为：

118、

119、

120、

121、其中，δ为变分符号，下标i表示爱因斯坦求和约定，下标k和上标k均表示第k个单元，m表示沿x轴、y轴或z轴的分量。加速度ak、位移xk及其导数添加下标i后,表示第k个单元在第i个节点的加速度；表示第k个单元在第i个节点的位移；表示第k个单元在m轴上速度；将整个区域划分为ne个单元，

122、

123、其中，f为待求物理量，vn为单元体积，sn为单元表面积，[a]是由插值函数组成的矩阵，其形式为，

124、

125、引入偏导数矩阵[b]，

126、

127、和应力向量[t]，

128、[t]＝(σ11,σ22,σ33,σ23,σ31,σ12)t            11)

129、其中，σ11,σ22,σ33,σ23,σ31,σ12分别表示不同方向的应力，下标1对应x轴，下标2对应y轴，下标3对应z轴；对于时间的推进，采用显示时间差分，

130、vn+1＝vn+an(δtn+δtn+1)/2                  12)

131、(δr)n+1＝(δr)n+vn+1·δtn+1                  13)

132、其中，δtn和δtn+1分别为n和n+1时刻的时间步长，δr和v分别为位移和速度。

133、用沙漏矢量抑制沙漏的产生，

134、

135、

136、

137、其中，下标k表示网格节点，i表示x轴、y轴或z轴的方向，α表示沙漏的各个方向，a＝qρl2c/4，q为常数，c是固体声速；l是单元特征长度；γαk是沙漏基矢量，γαk为人工粘性系数，为人工粘性项，nk,i为形函数，giα为沙漏质量，表示速度。

138、模块m3.2：将新的固体加速度as代入虚拟流体黎曼解算器式14)～16)，求解新的流体压强pl；

139、沿着右特征线

140、的相容关系为：

141、

142、将界面两侧的ρc视为常数，

143、

144、其中，ui、pi、ρil分别为界面流体速度、界面流体压力、界面流体密度，ul为流体相对于流固界面的速度；cil表示界面流体声速；pl表示界面附近流体压强；ρl表示表示界面附近流体密度；cl表示表示界面附近流体声速。

145、模块m3.3：回到模块m3.1，直到误差小于设定值；

146、所述模块m4包括：

147、定义为符号距离函数，

148、

149、其中，ω+表示虚拟区域，ω-表示真实区域，表示物质界面，dis表示距离，表示位置矢量。

150、将模块m3得到的流固界面压强和流固界面速度，用下述方程外推：

151、

152、其中，q为需要外推的函数；为流固界面法向单位向量，从流体指向固体；t为人工时间；把i从ω-外推到ω+时，方程18)取+；把i从ω+外推到ω-时，方程18)取-。

153、优选地，所述模块m5包括：

154、使用二阶tvd格式，设待求网格节点为(i0,j,k)，流场三个方向的通量f，g，h都需要离散，

155、f＝(ρu,ρu2+p,ρuv,ρuw,u(e+p))t             22)

156、g＝(ρv,ρvu,ρv2+p,ρvu,v(e+p))t             23)

157、h＝(ρu,ρwu,ρwv,ρw2+p,w(e+p))t             24)以f为例，

158、

159、其中，和表示半节点的值；ρu、ρv、ρu、ρuv、ρuw、ρvw均为与速度相关的动量，p为压强，e为单位体积流体的能量。

160、

161、这里，r为右特征矩阵；看作一个限制器，求解方程组如下：

162、

163、

164、

165、

166、

167、

168、

169、

170、这里，ε为一小量，λ为特征值，l为左特征矩阵，αm为特征向量，z为某一特征值，ψ，q和σ均为限制函数，δ为平均速度；

171、求解和时，需要用到半节点的值；采用roe平均求解下列物理量：

172、

173、

174、

175、

176、

177、其中，为平均值，

178、

179、和分别为半节点网格的密度、三个方向的速度和焓，接下来可间接求得能量和声速和

180、同理，对于任一网格节点(i,j,k)，依次求得

181、对于时间步长δt的选择，取

182、

183、cfl为收敛系数，cfl≤1；

184、所述模块m6中使用公式3)至16)求解结构体动态响应曲线和结构体的变形。

185、与现有技术相比，本发明具有如下的有益效果：

186、1、本发明基于虚拟流体方法的思想，对流固边界进迭代处理，能很好地计算深海耐压结构的动态响应；

187、2、本发明对固体结构使用有限元离散，能很好地模拟深海耐压结构液-固-气三相耦合的水下内爆的结构变形；

188、3、本发明对流场使用tvd差分格式，能很好地刻画流场中波的传播。

189、本发明的其他有益效果，将在具体实施方式中通过具体技术特征和技术方案的介绍来阐述，本领域技术人员通过这些技术特征和技术方案的介绍，应能理解所述技术特征和技术方案带来的有益技术效果。