[发明专利]基于移动粒子半隐式法的轴对称流场二维模拟方法

权利要求书

1.基于移动粒子半隐式法的轴对称流场二维模拟方法，其特征在于：步骤如下：

步骤1、对于流场具有旋转轴对称特征的计算区域，选择一个旋转面作为研究对象；在旋转面内根据流体类别采用不同粒子类型进行离散，并对粒子设置流体的物性参数，包括密度、粘性、表面张力；然后在旋转面粒子的作用半径内设置关于旋转轴的虚拟旋转粒子，虚拟旋转粒子所代表的流体类型和物性参数与对应的旋转面粒子相一致；公式(1)-(9)为将旋转面置于X-Z平面，旋转轴为Z轴时，虚拟旋转粒子的位置、速度和物性参数的设置方式：

ρ_jk＝ρ_j                            公式(8)

式中

x_j——旋转面粒子j的X轴方向坐标，m；

y_j——旋转面粒子j的Y轴方向坐标，m；

z_j——旋转面粒子j的Z轴方向坐标，m；

——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向坐标，m；

——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向坐标，m；

——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向坐标，m；

——旋转面粒子j所对应的虚拟旋转粒子k与旋转面的夹角，rad；

u_x,j——旋转面粒子j的X轴方向速度，m/s；

v_j——旋转面粒子j的Y轴方向速度，m/s；

w_j——旋转面粒子j的Z轴方向速度，m/s；

——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向速度，m/s；

——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向速度，m/s；

——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向速度，m/s；

Type_j——旋转面粒子j的类型；

ρ_j——旋转面粒子j的密度，kg/m³；

μ_j——旋转面粒子j的动力粘度，N·s/m²；

——旋转面粒子j所对应的虚拟旋转粒子k的类型；

——旋转面粒子j所对应的虚拟旋转粒子k的密度，kg/m³；

——旋转面粒子j所对应的虚拟旋转粒子k的粘度，N·s/m²；旋转轴粒子在计算中将其作为自由滑移边界；

步骤2、选取无奇点核函数，如公式(10)所示，然后计算旋转面初始粒子数密度，如公式(12)所示，并进行规则化处理，即根据旋转面粒子与旋转轴的距离，将虚拟旋转粒子扇形布置时计算出的旋转面粒子的粒子数密度转化为粒子等间距均匀布置时的粒子数密度，如公式(13)所示：

式中

r_e——粒子作用半径，m；

r——粒子间距，m；

w——核函数；

n₀——粒子等间距均匀布置时的粒子数密度；

——初始旋转面粒子i的粒子数密度；

——初始与旋转面粒子i在X方向相距l₀的旋转面粒子m的粒子数密度；

——第k步旋转面粒子i的粒子数密度；

——第k步规则化处理后的旋转面粒子i的粒子数密度；

——第k步旋转面粒子i的X方向坐标值，m；

——初始旋转面粒子i的X方向坐标值，m；

w(r_i-r_j)——j粒子对i粒子的核函数值，表达形式如公式(10)；

r_i⁰——初始i粒子位置矢量，m；

——初始旋转面粒子i邻近的旋转面粒子j的位置矢量，m；

——初始旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的位置矢量，m；

m_j——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子总数；

——初始邻近旋转面粒子j对旋转面粒子i的核函数值，表达形式如公式(10)；

——初始邻近旋转面粒子j所对应的虚拟旋转粒子k对旋转面粒子i的核函数值，表达形式如公式(10)；

步骤3、显式计算粘性项、表面张力项和重力项，计算旋转面粒子临时速度和位置；

采用连续表面张力模型来计算表面张力，其中颜色函数按公式(18)进行离散，且虚拟旋转粒子的颜色函数与所对应的旋转面粒子的颜色函数相等，如公式(19)所示，相界面处旋转面粒子法向量采用公式(21)进行离散，其中虚拟旋转粒子的法向量通过公式(22)-(24)计算：

C_jk＝C_j公式(19)

式中

r_i——旋转面粒子i的位置矢量，m；

r_j——旋转面粒子i邻近的旋转面粒子j的位置矢量，m；

d——维度；

C_i——旋转面粒子i的颜色函数；

C_j——旋转面粒子i邻近的旋转面粒子j的颜色函数；

——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的颜色函数；

——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的位置矢量，m；

n_i'——规则化处理后的旋转面粒子i的粒子数密度；

——邻近的旋转面粒子j所对应的虚拟旋转粒子k对旋转面粒子i的核函数值，表达形式如公式(10)；

式中

n_i——相界面处旋转面粒子i的法向量；

n_j——相界面处旋转面粒子i邻近的旋转面粒子j的法向量；

——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的法向量；

式中

n_i,x——相界面处旋转面粒子i的X方向分法向量；

n_i,y——相界面处旋转面粒子i的Y方向分法向量；

n_i,z——相界面处旋转面粒子i的Z方向分法向量；

——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的X方向分法向量；

——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Y方向分法向量；

——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Z方向分法向量；

对于粘性项则采用公式(25)进行计算：

式中

u_i——旋转面粒子i的速度矢量，m/s；

u_j——旋转面粒子i邻近的旋转面粒子j的速度矢量，m/s；

——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的速度矢量，m/s；

公式(25)中，

为了保证旋转面粒子仅在所建立的二维X-Z平面内运动，所以规定在Y轴方向旋转面粒子的位移为零，但在Y轴方向的速度不为零；

步骤4、根据旋转面粒子的临时速度和位置按照公式(1)-(6)计算与之对应的虚拟旋转粒子的临时速度和位置；

步骤5、按公式(13)规则化处理旋转面粒子移动后的粒子数密度，求解压力泊松方程来计算旋转面粒子压力；

其中压力泊松方程的源项采用临时速度散度结合粒子数密度的相对变化的形式，如公式(28)所示，对公式(28)左侧的拉普拉斯算子通过公式(29)进行离散，

式中

——旋转面粒子i规则化处理后的粒子数密度；

γ——截断误差补偿系数，0.01≤γ≤0.05；

式中

P_i^k+1——第k+1步旋转面粒子i的压力，Pa；

——第k+1步旋转面粒子i邻近的旋转面粒子j的压力，Pa；

——第k+1步旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a的压力，Pa；

——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a移动后的位置矢量，m；

结合公式(28)和(29)，可得公式(30)：

式中

a_i1——旋转面粒子i的第1个邻近粒子的压力系数；

a_i2——旋转面粒子i的第2个邻近粒子的压力系数；

a_ii——旋转面粒子i的压力系数；

a_iN+M——旋转面粒子i的第N+M个邻近粒子的压力系数；

P₁^k+1——第k+1步旋转面粒子i的第1个邻近粒子的压力，Pa；

——第k+1步旋转面粒子i的第2个邻近粒子的压力，Pa；

——第k+1步旋转面粒子i的第N+M个邻近粒子的压力，Pa；

b_i——旋转面粒子i的源项；

公式(30)中N为旋转面粒子总个数，M为虚拟旋转粒子总个数，其中：

公式(31)

式中

w(|r_j'-r_i|)——计算区域内任一粒子j'对旋转面粒子i的核函数值，表达形式如公式(10)；

将公式(30)变为矩阵表达式，如公式(33)所示：

由于虚拟旋转粒子的粒子数密度无法计算,导致公式(33)所述方程中未知变量个数多于方程个数，所以无法求解该压力泊松方程；因此，规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即P_jk＝P_j，则压力泊松方程系数(公式(31))变为：

则压力泊松方程矩阵表达式(公式(33))变为:

由公式(35)可知，压力泊松方程的系数矩阵为非对称矩阵，所以采用BICG求解器对压力泊松方程进行求解；

步骤6、通过步骤5更新流场压力后，求解压力梯度，然后修正旋转面粒子速度和位置；

其中压力梯度按公式(36)计算：

式中

——旋转面粒子j所对应的虚拟旋转粒子k的压力，Pa；

P_j——旋转面粒子j的压力，Pa；

因规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即上式变形为：

修正后的旋转面粒子在Y轴方向的位移被重新定义为零；

步骤7、根据旋转面粒子修正后的速度和位置按照公式(1)-(6)修正与之对应的虚拟旋转粒子的速度和位置。