[发明专利]一种基于连续介质力学的流固耦合高精度数值模拟方法

权利要求书

1.一种基于连续介质力学的流固耦合高精度数值模拟方法，其特征在于：包括如下步骤，

步骤1：确定流体、固体计算区域，建立直角坐标系，在直角坐标系将该区域划分为m*n个网格，其中，m表示x方向网格数量，n表示y方向网格数量；

步骤2：定义计算区域内初始状态下的用于区分界面位置的Level-Set函数、各区域物理量、建立基于连续介质力学的流体、固体控制方程组、材料状态方程函数，根据用于区分界面位置的Level-Set函数划分的各区域，并为上述函数及方程组各未知物理量赋初始值；

步骤3：对步骤2划分的各区域设置边界条件；

步骤4：选取计算参数CFL计算时间步长；

步骤5：基于步骤2各物理量的初始值赋值对计算区域初始化；

步骤6：根据Level-Set函数分别求解各区域多物质界面法向量；

步骤7：求解流体计算区域界面处的黎曼问题并确定流固耦合界面附近网格的物理状态；

步骤8：在计算固体时，将流体计算区域设置为虚拟固体网格，按照步骤7的方法求解固体计算区域界面处的黎曼问题并确定流固耦合界面附近网格的物理状态；

步骤9：采用高精度的WENO格式对单一介质计算区域进行空间离散；

步骤10：采用TVD Runge-Kutta格式对单一介质计算区域进行时间离散；

步骤11：重复步骤9-步骤10，对各物质区域分别进行空间、时间离散，得到tⁿ⁺¹时刻各物质区域的物理状态；

步骤12：整合步骤11得到的不同区域流场物理状态，得到全流场物理状态，对Level-Set运动方程的空间导数项进行离散并进行求解，得到tⁿ⁺¹时刻各物质区域Level-Set函数；

步骤13：基于当前计算的Level-Set全场值，对步骤9-步骤12得到的固体、流体流场进行选取，整合全流场物理状态，得到n+1时刻全场的物理状态；

步骤14:判断当前计算时间tⁿ⁺¹与设定总时间t^total的关系：若tⁿ⁺¹＞t^total，则结束数值模拟，输出计算区域内网格单元的压力值；若tⁿ⁺¹＜t^total，则返回步骤3，按照步骤13得到的物理状态设置边界条件，继续进行流固耦合的高精度数值模拟；

其中，步骤2实现方法为，

步骤2.1：定义计算区域内初始状态下的用于区分界面位置的level-Set函数，根据实际工况，通过LS方法将待计算的流体、固体计算区域划分为多区域；

定义计算区域内初始状态下的用于区分界面位置的Level-Set函数，基于Level-Set函数求解多物质界面的方法为LS方法；LS方法采用符号距离函数实现在任何时刻对界面演化过程的追踪，界面始终处于符号距离函数为零的等值面位置：

在计算初始时刻，采用公式(2)定义符号距离函数的初值：

其中，x、y、z分别为网格节点的物理空间坐标，d表示网格节点(x,y,z)到界面的距离，Ω₁、Ω₂分别为两种物质界面两侧的区域；

根据实际工况，通过LS方法将待计算的流体、固体计算区域划分为多区域；

步骤2.2：定义步骤2.1划分的各区域初始状态物理量；

所述初始状态物理量包括初始状态下流体的密度ρ、速度u、熵S、变形梯度张量F，初始状态固流体的密度、速度、熵、变形梯度张量；其中，变形梯度张量用于衡量物质存在初始拉伸、压缩、扭转的初始变形状态；

步骤2.3：构建多区域连续介质力学控制方程组；

由于流体、固体在自然界中均表示为连续介质，所以，由连续介质力学表示的欧拉方程组为：

其中，E为单位质量材料的总能量，ε为比内能，×为旋度算符、为梯度算符；

步骤2.4：选取材料的状态方程函数，进而确定材料在变形、运动下的物理状态；

不考虑化学反应的情况下，由于材料的内能由材料的变形和熵唯一决定，步骤2.4选取材料内能表示的状态方程函数如公式(4)所示，并确定内能形式状态方程函数的相关参数；

ε＝ε(g,S)＝ε(F,S)＝ε(C,S)＝ε(G,S) (4)

所述相关参数包括逆变形梯度张量g、变形梯度张量F，左柯西-格林应变张量C，右柯西-格林应变张量G；

步骤2.4所述内能形式状态方程函数优选用左柯西-格林应变张量来统一描述弹性固体以及流体的内能形式状态方程函数，如公式(5)所示：

其中，γ为多方气体指数，κ为熵的函数，a、p_∞为材料常数，χ为固体材料剪切弹性模量，对于流体材料χ＝0。