CN113807034A - 基于移动粒子半隐式法的轴对称流场二维模拟方法 - Google Patents

基于移动粒子半隐式法的轴对称流场二维模拟方法 Download PDF

Info

Publication number
CN113807034A
CN113807034A CN202111004067.9A CN202111004067A CN113807034A CN 113807034 A CN113807034 A CN 113807034A CN 202111004067 A CN202111004067 A CN 202111004067A CN 113807034 A CN113807034 A CN 113807034A
Authority
CN
China
Prior art keywords
particle
rotating
rotating surface
virtual
particles
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202111004067.9A
Other languages
English (en)
Other versions
CN113807034B (zh
Inventor
李�根
高金辰
王进仕
严俊杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN202111004067.9A priority Critical patent/CN113807034B/zh
Publication of CN113807034A publication Critical patent/CN113807034A/zh
Application granted granted Critical
Publication of CN113807034B publication Critical patent/CN113807034B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Fluid Mechanics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

基于移动粒子半隐式法的轴对称流场二维模拟方法，步骤如下：1、选择旋转面进行粒子初始布置，且在旋转面两侧布置虚拟旋转粒子，并设定参数；2、选取核函数，计算旋转面初始粒子数密度，并进行规则化处理；3、显式计算粘性项、表面张力项和重力项，预估旋转面粒子临时速度和位置；4、计算虚拟旋转粒子临时速度和位置；5、规则化处理移动后旋转面粒子的粒子数密度，并求解泊松方程来计算旋转面粒子压力；6、计算压力梯度，并修正旋转面粒子速度和位置；7、计算虚拟旋转粒子的速度和位置。本发明可将流场具有旋转轴对称特征的计算区域简化成二维旋转平面进行求解，且在保证准确性的前提下，极大地提高了移动粒子半隐式方法的计算效率。

Description

基于移动粒子半隐式法的轴对称流场二维模拟方法
技术领域
本发明涉及计算流体动力学技术领域，具体涉及一种能极大提高移动粒子半隐式法计算效率的计算方法。
背景技术
在计算流体动力学领域应用最广泛的是基于网格法的数值模拟方法，并且已经成为计算流体力学领域甚至固体力学领域中的主导方法。但传统的网格法在工程中的形变、相变、流固耦合及在复杂多相界面追踪和捕捉问题上存在许多困难和挑战。而粒子法在计算过程中，不需要任何网格参与，避免了所有基于网格法的局限性，并且计算节点之间没有固定的拓扑结构，这就使得无网格法在处理大变形问题时，展现了较网格法不可比拟的优势。这一点填补了网格法在数值模拟领域的巨大空白，使得高速撞击、动态断裂、塑性流动等涉及大变形情况下的数值模拟成为可能。
粒子法计算常采用拉格朗日坐标系和时间推进的描述方式，通过追踪目标粒子的运动轨迹及粒子携带的物理信息来掌握整个物理系统的动态变化信息，所以当参与计算的粒子数增多时，粒子法的计算量将会增大，导致计算效率较低，且随着计算区域的扩大和离散粒子数目的增加，这一情况将更加恶化。对于大规模复杂物理对象，开展全范围高分辨率的三维模拟所需的计算成本极高。目前传统基于网格划分的数值模拟方法开发了较为成熟的二维旋转轴对称模拟方法来近似三维模拟，但对于移动粒子半隐式方法中二维旋转轴对称计算方法还尚未见报道。
发明内容
为了克服上述现有技术存在的问题，本发明的目的在于提供一种基于移动粒子半隐式法的轴对称流场二维模拟方法，本发明针对具有旋转轴对称特征的流场，将其三维流场结构简化为二维旋转平面，利用本发明提出的方法进行计算，可以在保证一定精度的前提下，极大的缩减计算时间，从而提高计算效率。
为了达到上述目的，本发明采用的技术方案如下：
基于移动粒子半隐式法的轴对称流场二维模拟方法，该方法包括以下步骤：
步骤1、对于流场具有旋转轴对称特征的计算区域，选择一个旋转面作为研究对象；在旋转面内根据流体类别采用不同粒子类型进行离散，并对粒子设置流体的物性参数，包括密度、粘性、表面张力；然后在旋转面粒子的作用半径内设置关于旋转轴的虚拟旋转粒子，虚拟旋转粒子所代表的流体类型和物性参数与对应的旋转面粒子相一致；公式(1)-(9)为将旋转面置于X-Z平面，旋转轴为Z轴时，虚拟旋转粒子的位置、速度和物性参数的设置方式：
Figure BDA0003236532650000021
Figure BDA0003236532650000022
Figure BDA0003236532650000023
Figure BDA0003236532650000024
Figure BDA0003236532650000025
Figure BDA0003236532650000026
Figure BDA0003236532650000027
ρjk＝ρj 公式(8)
Figure BDA0003236532650000031
式中
xj——旋转面粒子j的X轴方向坐标，m；
yj——旋转面粒子j的Y轴方向坐标，m；
zj——旋转面粒子j的Z轴方向坐标，m；
Figure BDA0003236532650000032
——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向坐标，m；
Figure BDA0003236532650000033
——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向坐标，m；
Figure BDA0003236532650000034
——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向坐标，m；
Figure BDA0003236532650000035
——旋转面粒子j所对应的虚拟旋转粒子k与旋转面的夹角，rad；
uj——旋转面粒子j的X轴方向速度，m/s；
vj——旋转面粒子j的Y轴方向速度，m/s；
wj——旋转面粒子j的Z轴方向速度，m/s；
Figure BDA0003236532650000036
——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向速度，m/s；
Figure BDA0003236532650000037
——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向速度，m/s；
Figure BDA0003236532650000038
——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向速度，m/s；
Typej——旋转面粒子j的类型；
ρj——旋转面粒子j的密度，kg/m3；
μj——旋转面粒子j的动力粘度，N·s/m2；
Figure BDA0003236532650000039
——旋转面粒子j所对应的虚拟旋转粒子k的类型；
Figure BDA00032365326500000310
——旋转面粒子j所对应的虚拟旋转粒子k的密度，kg/m3；
Figure BDA00032365326500000311
——旋转面粒子j所对应的虚拟旋转粒子k的粘度，N·s/m2；
旋转轴粒子在计算中将其作为自由滑移边界；
步骤2、选取无奇点核函数，然后计算旋转面初始粒子数密度，如公式(12)所示，并进行规则化处理，即根据旋转面粒子与旋转轴的距离，将虚拟旋转粒子扇形布置时计算出的旋转面粒子的粒子数密度转化为粒子等间距均匀布置时的粒子数密度，如公式(13)所示：
Figure BDA0003236532650000041
Figure BDA0003236532650000042
式中
n0——粒子等间距均匀布置时的粒子数密度；
Figure BDA0003236532650000043
——初始旋转面粒子i的粒子数密度；
Figure BDA0003236532650000044
——初始与旋转面粒子i在X方向相距l0的旋转面粒子m的粒子数密度；
Figure BDA0003236532650000045
——第k步旋转面粒子i的粒子数密度；
Figure BDA0003236532650000046
——第k步规则化处理后的旋转面粒子i的粒子数密度；
Figure BDA0003236532650000047
——第k步旋转面粒子i的X方向坐标值，m；
Figure BDA0003236532650000048
——初始旋转面粒子i的X方向坐标值，m；
w(|ri-rj|)——j粒子对i粒子的核函数值，表达形式如公式(10)；
ri 0——初始i粒子位置矢量，m；
Figure BDA0003236532650000049
——初始旋转面粒子i邻近的旋转面粒子j的位置矢量，m；
Figure BDA00032365326500000410
——初始旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的位置矢量，m；
mj——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子总数；
Figure BDA00032365326500000411
——初始邻近旋转面粒子j对旋转面粒子i的核函数值，表达形式如公式(10)；
Figure BDA00032365326500000412
——初始邻近旋转面粒子j所对应的虚拟旋转粒子k对旋转面粒子i的核函数值，表达形式如公式(10)；
步骤3、显式计算粘性项、表面张力项和重力项，计算旋转面粒子临时速度和位置；
采用连续表面张力模型来计算表面张力，其中颜色函数按公式(18)进行离散，且虚拟旋转粒子的颜色函数与所对应的旋转面粒子的颜色函数相等，如公式(19)所示，相界面处旋转面粒子法向量采用公式(21)进行离散，其中虚拟旋转粒子的法向量通过公式(22)-(24)计算：
Figure BDA0003236532650000051
Figure BDA0003236532650000052
式中
ri——旋转面粒子i的位置矢量，m；
rj——旋转面粒子i邻近的旋转面粒子j的位置矢量，m；
d——维度；
Ci——旋转面粒子i的颜色函数；
Cj——旋转面粒子i邻近的旋转面粒子j的颜色函数；
Figure BDA0003236532650000053
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的颜色函数；
Figure BDA0003236532650000054
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的位置矢量，m；
n'i——规则化处理后的旋转面粒子i的粒子数密度；
Figure BDA0003236532650000055
——邻近的旋转面粒子j所对应的虚拟旋转粒子k对旋转面粒子i的核函数值，表达形式如公式(10)；
Figure BDA0003236532650000056
式中
ni——相界面处旋转面粒子i的法向量；
nj——相界面处旋转面粒子i邻近的旋转面粒子j的法向量；
Figure BDA0003236532650000061
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的法向量；
Figure BDA0003236532650000062
Figure BDA0003236532650000063
Figure BDA0003236532650000064
式中
ni,x——相界面处旋转面粒子i的X方向分法向量；
ni,y——相界面处旋转面粒子i的Y方向分法向量；
ni,z——相界面处旋转面粒子i的Z方向分法向量；
Figure BDA0003236532650000065
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的X方向分法向量；
Figure BDA0003236532650000066
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Y方向分法向量；
Figure BDA0003236532650000067
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Z方向分法向量；
对于粘性项则采用公式(25)进行计算：
Figure BDA0003236532650000068
式中
ui——旋转面粒子i的速度矢量，m/s；
uj——旋转面粒子i邻近的旋转面粒子j的速度矢量，m/s；
Figure BDA0003236532650000069
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的速度矢量，m/s；
公式(25)中，
Figure BDA0003236532650000071
为了保证旋转面粒子仅在在所建立的二维X-Z平面内运动，所以规定在Y轴方向旋转面粒子的位移为零，但在Y轴方向的速度不为零；
步骤4、根据旋转面粒子的临时速度和位置按照公式(1)-(6)计算与之对应的虚拟旋转粒子的临时速度和位置；
步骤5、按公式(13)规则化处理旋转面粒子移动后的粒子数密度，求解压力泊松方程来计算旋转面粒子压力；
其中压力泊松方程的源项采用临时速度散度结合粒子数密度的相对变化的形式，如公式(28)所示，对公式(28)左侧的拉普拉斯算子通过公式(29)进行离散，
Figure BDA0003236532650000072
式中
Figure BDA0003236532650000073
——旋转面粒子i规则化处理后的粒子数密度；
γ——截断误差补偿系数，0.01≤γ≤0.05；
Figure BDA0003236532650000074
式中
Pi k+1——第k+1步旋转面粒子i的压力，Pa；
Figure BDA0003236532650000075
——第k+1步旋转面粒子i邻近的旋转面粒子j的压力，Pa；
Figure BDA0003236532650000081
——第k+1步旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a的压力，Pa；
Figure BDA0003236532650000082
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a移动后的位置矢量，m；
结合公式(28)和(29)，可得公式(30)：
Figure BDA0003236532650000083
式中
ai1——旋转面粒子i的第1个邻近粒子的压力系数；
ai2——旋转面粒子i的第2个邻近粒子的压力系数；
aii——旋转面粒子i的压力系数；
aiN+M——旋转面粒子i的第N+M个邻近粒子的压力系数；
P1 k+1——第k+1步旋转面粒子i的第1个邻近粒子的压力，Pa；
Figure BDA0003236532650000084
——第k+1步旋转面粒子i的第2个邻近粒子的压力，Pa；
Figure BDA0003236532650000085
——第k+1步旋转面粒子i的第N+M个邻近粒子的压力，Pa；
bi——旋转面粒子i的源项；
公式(30)中N为旋转面粒子总个数，M为虚拟旋转粒子总个数，其中：
Figure BDA0003236532650000086
Figure BDA0003236532650000087
式中
w(|rj'-ri|)——计算区域内任一粒子j'对旋转面粒子i的核函数值，表达形式如公式(10)；
将公式(30)变为矩阵表达式，如公式(33)所示：
Figure BDA0003236532650000091
由于虚拟旋转粒子的粒子数密度无法计算,导致公式(33)所述方程中未知变量个数多于方程个数，所以无法求解该压力泊松方程；因此，规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即
Figure BDA0003236532650000092
则压力泊松方程系数(公式(31))变为：
Figure BDA0003236532650000093
则压力泊松方程矩阵表达式(公式(33))变为:
Figure BDA0003236532650000094
由公式(35)可知，压力泊松方程的系数矩阵为非对称矩阵，所以采用BICG求解器对压力泊松方程进行求解；
步骤6、通过步骤5更新流场压力后，求解压力梯度，然后修正旋转面粒子速度和位置；
其中压力梯度按公式(36)计算：
Figure BDA0003236532650000101
式中
Figure BDA0003236532650000102
——旋转面粒子j所对应的虚拟旋转粒子k的压力，Pa；
Pj——旋转面粒子j的压力，Pa；
因规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即
Figure BDA0003236532650000103
上式变形为：
Figure BDA0003236532650000104
修正后的旋转面粒子在Y轴方向的位移被重新定义为零；
步骤7、根据旋转面粒子修正后的速度和位置按照公式(1)-(6)修正与之对应的虚拟旋转粒子的速度和位置。
和现有技术相比，本发明提出的计算方法具备如下优点：
利用移动粒子半隐式方法模拟流场运动时，随着流场范围的扩大及分辨率的提高，参与计算的粒子数会明显增多，导致计算量会显著增加。而对于流场特性在任意轴向截面上均匀分布的旋转体流体，由于其具有明显的旋转轴对称特征，可以利用本发明提出的模拟方法将三维流场简化成二维旋转平面进行计算，这种方法能够大幅缩减粒子数量，减少计算成本，提高计算效率，并且具有较高准确性。
附图说明
图1是本发明提出的基于移动粒子半隐式法的轴对称流场二维模拟方法的流程图。
图2是旋转面及其虚拟旋转粒子布置示意图。
图3是圆柱型毛细射流旋转面粒子初始布置图。
图4a、图4b和图4c分别是τ＝0、τ＝12、τ＝15在密度比ρ2/ρ1＝0.8，粘度比μ2/μ1＝0.01，雷诺数Re＝10,扰动波数k＝0.7工况下，射流随时间断裂的模拟结果；τ为射流断裂的无量纲时间。
图5a和图5b是不同扰动波数下主、卫星液滴直径模拟值与实验值比较结果。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步详细说明。
本发明提出的基于移动粒子半隐式法的轴对称流场二维模拟方法的计算流程如图1所示，计算步骤如下：
步骤1、对于流场具有旋转轴对称特征的计算区域，选择一个旋转面作为研究对象；在旋转面内根据流体类别采用不同粒子类型进行离散，并对粒子设置流体的物性参数，含密度、粘性、表面张力；然后在旋转面粒子的作用半径内设置关于旋转轴的虚拟旋转粒子，虚拟旋转粒子所代表的流体类型和物性参数与对应的旋转面粒子相一致；公式(1)-(9)为将旋转面置于X-Z平面，旋转轴为Z轴时，虚拟旋转粒子的位置、速度和物性参数的设置方式：
Figure BDA0003236532650000111
Figure BDA0003236532650000112
Figure BDA0003236532650000113
Figure BDA0003236532650000121
Figure BDA0003236532650000122
Figure BDA0003236532650000123
Figure BDA0003236532650000124
ρjk＝ρj 公式(8)
Figure BDA0003236532650000125
式中
xj——旋转面粒子j的X轴方向坐标，m；
yj——旋转面粒子j的Y轴方向坐标，m；
zj——旋转面粒子j的Z轴方向坐标，m；
Figure BDA0003236532650000126
——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向坐标，m；
Figure BDA0003236532650000127
——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向坐标，m；
Figure BDA0003236532650000128
——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向坐标，m；
Figure BDA0003236532650000129
——旋转面粒子j所对应的虚拟旋转粒子k与旋转面的夹角，rad；
uj——旋转面粒子j的X轴方向速度，m/s；
vj——旋转面粒子j的Y轴方向速度，m/s；
wj——旋转面粒子j的Z轴方向速度，m/s；
Figure BDA00032365326500001210
——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向速度，m/s；
Figure BDA00032365326500001211
——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向速度，m/s；
Figure BDA00032365326500001212
——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向速度，m/s；
Typej——旋转面粒子j的类型；
ρj——旋转面粒子j的密度，kg/m3；
μj——旋转面粒子j的动力粘度，N·s/m2；
Figure BDA0003236532650000131
——旋转面粒子j所对应的虚拟旋转粒子k的类型；
Figure BDA0003236532650000132
——旋转面粒子j所对应虚的拟旋转粒子k的密度，kg/m3；
Figure BDA0003236532650000133
——旋转面粒子j所对应的虚拟旋转粒子k的粘度，N·s/m2；
对于图2中的旋转轴粒子在计算中将其作为自由滑移边界；
步骤2、在移动粒子半隐式方法中，采用核函数表征粒子间相互作用的强弱程度。本发明采用文献[1](KOSHIZUKA S,NOBE A,OKA Y.Numerical analysis of breakingwaves using the moving particle semi-implicit method[J].international Journalfor Numerical Methods in Fluids,1998,26(7):751-769)中无奇点核函数，如公式(10)所示：
Figure BDA0003236532650000134
式中
re——粒子作用半径，m；
r——粒子间距，m；
w——核函数；
且在移动粒子半隐式法中，是通过维持粒子数密度为常数来保证流体的不可压缩性。而在计算过程中，由于本发明中虚拟旋转粒子的扇形排布方式导致在压力作用下粒子数密度很难自主靠近初始粒子数密度n0，所以需根据旋转面粒子与旋转轴的距离，将虚拟旋转粒子扇形布置时计算出的旋转面粒子的粒子数密度转化为粒子等间距均匀布置时的粒子数密度，即规则化处理。首先计算初始时刻，粒子规则分布情况下的粒子数密度，即初始粒子数密度n0，如公式(11)所示；再计算初始时刻在旋转面两侧布置旋转虚拟粒子情况下的旋转面粒子的粒子数密度
Figure BDA0003236532650000135
如公式(12)所示；最后进行规则化处理，如公式(13)所示：
Figure BDA0003236532650000141
Figure BDA0003236532650000142
Figure BDA0003236532650000143
式中
ri——等间距均匀布置时初始i粒子位置矢量，m；
rj——等间距均匀布置时初始i粒子邻近的粒子j的位置矢量，m；
n0——粒子等间距均匀布置时的粒子数密度；
Figure BDA0003236532650000144
——初始旋转面粒子i的粒子数密度；
Figure BDA0003236532650000145
——初始与旋转面粒子i在X方向相距l0的旋转面粒子m的粒子数密度；
Figure BDA0003236532650000146
——第k步旋转面粒子i的粒子数密度；
Figure BDA0003236532650000147
——第k步规则化处理后的旋转面粒子i的粒子数密度；
Figure BDA0003236532650000148
——第k步旋转面粒子i的X方向坐标值，m；
Figure BDA0003236532650000149
——初始旋转面粒子i的X方向坐标值，m；
w(|ri-rj|)——j粒子对i粒子的核函数值，表达形式如公式(10)；
ri 0——初始i粒子位置矢量，m；
Figure BDA00032365326500001410
——初始旋转面粒子i邻近的旋转面粒子j的位置矢量，m；
Figure BDA00032365326500001411
——初始旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的位置矢量，m；
mj——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子总数；
Figure BDA00032365326500001412
——初始邻近旋转面粒子j对旋转面目标粒子i的核函数值，表达形式如公式(10)；
Figure BDA00032365326500001413
——初始邻近旋转面粒子j所对应的虚拟旋转粒子k对旋转面目标粒子i的核函数值，表达形式如公式(10)；
步骤3、估算临时速度场，首先通过显式求解动量方程公式(14)中的重力项、表面张力项及粘性项，对旋转面粒子的临时速度进行估算，再估算其临时位置：
Figure BDA0003236532650000151
式中
u——粒子速度矢量，m/s；
P——粒子压力，Pa；
F——粒子表面张力矢量，N/kg；
g——粒子重力加速度矢量，m/s2；
μ——粒子动力粘度，N·s/m2；
ρ——粒子密度，kg/m3；
t——时间，s；
D——随体导数；
显式求解动量方程公式(14)中的重力项，如公式(15)所示：
G＝ρg
公式(15)
式中
G——重力矢量，N；
显式求解动量方程公式(14)中的表面张力项，本发明采用文献[2](F.Yeganehdoust,M.Yaghoubi,H.Emdad,et al.Numerical study of multiphasedroplet dynamics and contact angles by smoothed particlehydrodynamics.Applied Mathematical Modelling,2016,40:8493-8512)中基于相界面法向量梯度的连续表面张力模型，连续表面张力模型将相界面上的表面力转化成体积力，则表面张力项可以写成如下形式，表达式如公式(16)所示：
Figure BDA0003236532650000161
式中
Fs——表面张力，N；
σ——表面张力系数；
κi——相界面处旋转面粒子i的曲率；
Figure BDA0003236532650000162
——旋转面粒子i的颜色函数梯度；
首先采用一个颜色函数C来标记不同相内的粒子，如公式(17)所示：
Figure BDA0003236532650000163
式中
Ci——i粒子的颜色函数
颜色函数通过梯度模型进行离散，如公式(18)所示，且虚拟旋转粒子的颜色函数与所对应的旋转面粒子的颜色函数相等，如公式(19)所示：
Figure BDA0003236532650000164
Figure BDA0003236532650000165
式中
ri——目标旋转面粒子i的位置矢量，m；
rj——目标旋转面粒子i邻近的旋转面粒子j的位置矢量，m；
d——维度；
Ci——目标旋转面粒子i的颜色函数；
Cj——目标旋转面粒子i邻近的旋转面粒子j的颜色函数；
Figure BDA0003236532650000166
——目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的颜色函数；
Figure BDA0003236532650000167
——目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k位置矢量，m；
n'i——规则化处理后的目标旋转面粒子i的粒子数密度；
Figure BDA0003236532650000171
——邻近的旋转面粒子j所对应的虚拟旋转粒子k对目标旋转面粒子i的核函数值，表达形式如公式(10)；
通过散度模型离散相界面处的粒子法向量来计算相界面处旋转粒子的曲率，如公式(20)所示：
Figure BDA0003236532650000172
Figure BDA0003236532650000173
式中
ni——相界面处旋转面粒子i的法向量；
nj——相界面处旋转面粒子i邻近的旋转面粒子j的法向量；
Figure BDA0003236532650000174
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的法向量；
虚拟旋转粒子的法向量通过公式(22)-(24)计算：
Figure BDA0003236532650000175
Figure BDA0003236532650000176
Figure BDA0003236532650000177
式中
ni,x——相界面处旋转面粒子i的X方向分法向量；
ni,y——相界面处旋转面粒子i的Y方向分法向量；
ni,z——相界面处旋转面粒子i的Z方向分法向量；
Figure BDA0003236532650000178
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的X方向分法向量；
Figure BDA0003236532650000181
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Y方向分法向量；
Figure BDA0003236532650000182
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Z方向分法向量；
显式求解动量方程公式(14)中的粘性项，通过上一时层计算出的旋转面粒子的临时速度及其虚拟旋转粒子的临时速度来显式求解动量方程公式(14)中的粘性项，采用拉普拉斯算子模型对公式(14)中的拉普拉斯算子进行离散，如公式(25)所示：
Figure BDA0003236532650000183
式中
ui——目标旋转面粒子i的速度矢量，m/s；
uj——目标旋转面粒子i邻近的旋转面粒子j的速度矢量，m/s；
Figure BDA0003236532650000184
——目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的速度矢量，m/s；
上式中
Figure BDA0003236532650000185
通过公式(15)、公式(16)和公式(25)，求解动量方程公式(14)中的重力项、表面张力项和粘性项，再通过公式(26)对旋转面粒子的临时速度进行估算：
Figure BDA0003236532650000186
式中
Figure BDA0003236532650000191
——旋转面粒子i当前的速度矢量，m/s；
Figure BDA0003236532650000192
——旋转面粒子i估算速度矢量，m/s；
通过计算出的旋转面粒子的临时速度来估算旋转面粒子的临时位置，如公式(27)所示，为了保证旋转面粒子仅在所建立的二维X-Z平面内运动，所以移动后的旋转面粒子在Y轴方向的位移被定义为零，但在Y轴方向的速度不为零:
Figure BDA0003236532650000193
式中
ri *——移动后旋转面粒子i的估算位置矢量，m；
ri k——当前时刻旋转面粒子i的位置矢量，m。
步骤4、预估虚拟旋转粒子的临时速度和位置：根据旋转面粒子的临时速度和位置按照公式(1)-(6)计算与之对应的虚拟旋转粒子的临时速度和位置；
步骤5、由于上述显示计算部分未考虑压力项，故而粒子在临时移动后，规范化处理后的粒子数密度
Figure BDA0003236532650000194
会偏移初始粒子数密度n0，为了保证流体的不可压缩性，即保证粒子数密度为常数，需要用质量守恒方程对其进行隐式修正。由不可压缩模型推导得流体速度场及压力场的耦合关系式——压力泊松方程，本发明中压力泊松方程的源项采用文献[3](Khayyer A,Gotoh H.Modified moving particle semi-implicit methods for theprediction of 2D wave impact pressure[J].Coastal engineering,2009,56(4):419-440)中临时速度的散度结合粒子数密度的相对变化的形式，如公式(28)所示，然后隐式求解压力泊松方程：
Figure BDA0003236532650000195
式中
Figure BDA0003236532650000196
——旋转面粒子i规则化处理后的粒子数密度；
γ——截断误差补偿系数，0.01≤γ≤0.05；
对公式(28)等式左侧的拉普拉斯算子通过拉普拉斯模型进行离散，如公式(29)所示：
Figure BDA0003236532650000201
式中
Pi k+1——第k+1步旋转面粒子i的压力，Pa；
Figure BDA0003236532650000202
——第k+1步旋转面粒子i邻近的旋转面粒子j的压力，Pa；
Figure BDA0003236532650000203
——第k+1步旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a的压力，Pa；
Figure BDA0003236532650000204
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a移动后的位置矢量，m；
结合公式(28)和(29)，可得公式(30)：
Figure BDA0003236532650000205
式中
ai1——旋转面粒子i的第1个邻近粒子的压力系数；
ai2——旋转面粒子i的第2个邻近粒子的压力系数；
aii——旋转面粒子i的压力系数；
aiN+M——旋转面粒子i的第N+M个邻近粒子的压力系数；
P1 k+1——第k+1步旋转面粒子i的第1个邻近粒子的压力，Pa；
Figure BDA0003236532650000206
——第k+1步旋转面粒子i的第2个邻近粒子的压力，Pa；
Figure BDA0003236532650000207
——第k+1步旋转面粒子i的第N+M个邻近粒子的压力，Pa；
bi——旋转面粒子i的源项；
公式(30)中N为旋转面粒子总个数，M为虚拟旋转粒子总个数,其中：
Figure BDA0003236532650000211
Figure BDA0003236532650000212
式中
w(|rj'-ri|)——计算区域内任一粒子j'对旋转面粒子i的核函数值，表达形式如公式(10)；
将公式(30)变为矩阵表达式，如公式(33)所示：
Figure BDA0003236532650000213
由于虚拟旋转粒子的粒子数密度无法计算,导致公式(33)所述方程中未知变量个数多于方程个数，所以无法求解该压力泊松方程。因此，本发明规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即
Figure BDA0003236532650000215
则压力泊松方程系数(公式(31))可变为：
Figure BDA0003236532650000214
则压力泊松方程矩阵表达式(公式(33))可变为:
Figure BDA0003236532650000221
由公式(35)可知，压力泊松方程的系数矩阵为非对称矩阵，所以本发明采用BICG求解器对压力泊松方程进行求解；
步骤6、通过步骤5更新压力场后，求解动量方程中压力梯度项，然后修正旋转面粒子的速度和位置，压力梯度项按公式(36)进行离散：
Figure BDA0003236532650000222
式中
Figure BDA0003236532650000223
——旋转面粒子j所对应的虚拟旋转粒子k的压力，Pa；
Pj——旋转面粒子j的压力，Pa；
因本发明规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即
Figure BDA0003236532650000224
上式可变形为：
Figure BDA0003236532650000225
利用求解出的压力梯度项对旋转面粒子的速度和位置进行修正，如公式(38)和公式(39)所示：
Figure BDA0003236532650000226
Figure BDA0003236532650000227
式中
Figure BDA0003236532650000231
——第k+1步旋转面粒子i的速度矢量，m/s；
Figure BDA0003236532650000232
——第k+1步旋转面粒子i的位置矢量，m；
修正后的旋转面粒子在Y轴方向的位移被重新定义为零；
步骤7、修正虚拟旋转粒子的速度和位置：根据旋转面粒子修正后的速度和位置按照公式(1)-(6)修正与之对应的虚拟旋转粒子的速度和位置。
为了验证本发明提出的基于移动粒子半隐式法的轴对称流场二维模拟方法的准确性及稳定性，本发明采用文献[4](Absar M.Lakdawala and Rochish Thaokar.DGLSMbased study of temporal instability and formation of satellite drop in acapillary jet breakup[J].Chemical Engineering Science,2015,130:239-253)中的数据来模拟毛细管射流断裂。
图3展示了圆柱型毛细管射流在扰动波数k＝0.7的旋转面初始粒子布置，初始均匀布置的粒子间距为0.00025m,粒子总数为15428。液体射流不受重力影响，液-液界面的表面张力系数恒定。在初始时刻射流处于静止状态，且对液-液界面施加轴对称余弦扰动，扰动波函数为Rint＝R1+A0cos(2πZ/λ)，R1＝1，扰动波幅值A0＝0.01R1，圆柱高λ＝2π/k，圆柱半径R＝λ/4，与参考文献[4]保持一致。
图4a、图4b和图4c分别展示了在密度比ρ2/ρ1＝0.8，粘度比μ2/μ1＝0.01，雷诺数Re＝10,扰动波数k＝0.7的工况下(参考文献[4]保持一致)，射流在无量纲时间τ＝0、12、15时的模拟结果，图中黑色射流线型为文献1中的结果，可以看出射流在断裂过程中，射流的液-液界面形状与参考文献中的结果吻合较好。
图5a和图5b分别定量比较了利用本发明提出的方法模拟得到的不同扰动波数下(k＝0.5、0.7、0.9)的主、卫星液滴直径模拟值与文献[4]中的模拟值，文献[5](Rutland,Eand Jameson,J.Theoretical prediction of the sizes of droplets formed in thebreakup of capillary jets.Chemical Engineering Science,1970.25:1689-1698)及其文献[6](Lafrance，P.Nonlinear breakup of a laminar liquid jet.The Physics ofFluids,1974.18:428-432)中的实验值，可以看出，利用本发明提出的方法模拟得到的主、卫星液滴直径的模拟值与实验值吻合较好。

Claims (1)

1.基于移动粒子半隐式法的轴对称流场二维模拟方法，其特征在于：步骤如下：
步骤1、对于流场具有旋转轴对称特征的计算区域，选择一个旋转面作为研究对象；在旋转面内根据流体类别采用不同粒子类型进行离散，并对粒子设置流体的物性参数，包括密度、粘性、表面张力；然后在旋转面粒子的作用半径内设置关于旋转轴的虚拟旋转粒子，虚拟旋转粒子所代表的流体类型和物性参数与对应的旋转面粒子相一致；公式(1)-(9)为将旋转面置于X-Z平面，旋转轴为Z轴时，虚拟旋转粒子的位置、速度和物性参数的设置方式：
Figure FDA0003236532640000011
Figure FDA0003236532640000012
zjk＝zj 公式(3)
Figure FDA0003236532640000013
Figure FDA0003236532640000014
Figure FDA0003236532640000015
Figure FDA0003236532640000016
ρjk＝ρj 公式(8)
Figure FDA0003236532640000017
式中
xj——旋转面粒子j的X轴方向坐标，m；
yj——旋转面粒子j的Y轴方向坐标，m；
zj——旋转面粒子j的Z轴方向坐标，m；
Figure FDA0003236532640000021
——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向坐标，m；
Figure FDA0003236532640000022
——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向坐标，m；
Figure FDA0003236532640000023
——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向坐标，m；
Figure FDA0003236532640000024
——旋转面粒子j所对应的虚拟旋转粒子k与旋转面的夹角，rad；
uj——旋转面粒子j的X轴方向速度，m/s；
vj——旋转面粒子j的Y轴方向速度，m/s；
wj——旋转面粒子j的Z轴方向速度，m/s；
Figure FDA0003236532640000025
——旋转面粒子j所对应的虚拟旋转粒子k的X轴方向速度，m/s；
Figure FDA0003236532640000026
——旋转面粒子j所对应的虚拟旋转粒子k的Y轴方向速度，m/s；
Figure FDA0003236532640000027
——旋转面粒子j所对应的虚拟旋转粒子k的Z轴方向速度，m/s；
Typej——旋转面粒子j的类型；
ρj——旋转面粒子j的密度，kg/m3；
μj——旋转面粒子j的动力粘度，N·s/m2；
Figure FDA0003236532640000028
——旋转面粒子j所对应的虚拟旋转粒子k的类型；
Figure FDA0003236532640000029
——旋转面粒子j所对应的虚拟旋转粒子k的密度，kg/m3；
Figure FDA00032365326400000210
——旋转面粒子j所对应的虚拟旋转粒子k的粘度，N·s/m2；旋转轴粒子在计算中将其作为自由滑移边界；
步骤2、选取无奇点核函数，然后计算旋转面初始粒子数密度，如公式(12)所示，并进行规则化处理，即根据旋转面粒子与旋转轴的距离，将虚拟旋转粒子扇形布置时计算出的旋转面粒子的粒子数密度转化为粒子等间距均匀布置时的粒子数密度，如公式(13)所示：
Figure FDA00032365326400000211
Figure FDA0003236532640000031
式中
n0——粒子等间距均匀布置时的粒子数密度；
Figure FDA0003236532640000032
——初始旋转面粒子i的粒子数密度；
Figure FDA0003236532640000033
——初始与旋转面粒子i在X方向相距l0的旋转面粒子m的粒子数密度；
Figure FDA0003236532640000034
——第k步旋转面粒子i的粒子数密度；
Figure FDA0003236532640000035
——第k步规则化处理后的旋转面粒子i的粒子数密度；
Figure FDA0003236532640000036
——第k步旋转面粒子i的X方向坐标值，m；
Figure FDA0003236532640000037
——初始旋转面粒子i的X方向坐标值，m；
w(|ri-rj|)——j粒子对i粒子的核函数值，表达形式如公式(10)；
ri 0——初始i粒子位置矢量，m；
Figure FDA0003236532640000038
——初始旋转面粒子i邻近的旋转面粒子j的位置矢量，m；
Figure FDA0003236532640000039
——初始旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的位置矢量，m；
mj——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子总数；
Figure FDA00032365326400000310
——初始邻近旋转面粒子j对旋转面粒子i的核函数值，表达形式如公式(10)；
Figure FDA00032365326400000311
——初始邻近旋转面粒子j所对应的虚拟旋转粒子k对旋转面粒子i的核函数值，表达形式如公式(10)；
步骤3、显式计算粘性项、表面张力项和重力项，计算旋转面粒子临时速度和位置；
采用连续表面张力模型来计算表面张力，其中颜色函数按公式(18)进行离散，且虚拟旋转粒子的颜色函数与所对应的旋转面粒子的颜色函数相等，如公式(19)所示，相界面处旋转面粒子法向量采用公式(21)进行离散，其中虚拟旋转粒子的法向量通过公式(22)-(24)计算：
Figure FDA0003236532640000041
Figure FDA0003236532640000042
式中
ri——旋转面粒子i的位置矢量，m；
rj——旋转面粒子i邻近的旋转面粒子j的位置矢量，m；
d——维度；
Ci——旋转面粒子i的颜色函数；
Cj——旋转面粒子i邻近的旋转面粒子j的颜色函数；
Figure FDA0003236532640000043
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的颜色函数；
Figure FDA0003236532640000044
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的位置矢量，m；
n′i——规则化处理后的旋转面粒子i的粒子数密度；
Figure FDA0003236532640000045
——邻近的旋转面粒子j所对应的虚拟旋转粒子k对旋转面粒子i的核函数值，表达形式如公式(10)；
Figure FDA0003236532640000046
式中
ni——相界面处旋转面粒子i的法向量；
nj——相界面处旋转面粒子i邻近的旋转面粒子j的法向量；
Figure FDA0003236532640000047
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的法向量；
Figure FDA0003236532640000051
Figure FDA0003236532640000052
Figure FDA0003236532640000053
式中
ni,x——相界面处旋转面粒子i的X方向分法向量；
ni,y——相界面处旋转面粒子i的Y方向分法向量；
ni,z——相界面处旋转面粒子i的Z方向分法向量；
Figure FDA0003236532640000054
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的X方向分法向量；
Figure FDA0003236532640000055
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Y方向分法向量；
Figure FDA0003236532640000056
——相界面处目标旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的Z方向分法向量；
对于粘性项则采用公式(25)进行计算：
Figure FDA0003236532640000057
式中
ui——旋转面粒子i的速度矢量，m/s；
uj——旋转面粒子i邻近的旋转面粒子j的速度矢量，m/s；
Figure FDA0003236532640000058
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子k的速度矢量，m/s；
公式(25)中，
Figure FDA0003236532640000061
为了保证旋转面粒子仅在在所建立的二维X-Z平面内运动，所以规定在Y轴方向旋转面粒子的位移为零，但在Y轴方向的速度不为零；
步骤4、根据旋转面粒子的临时速度和位置按照公式(1)-(6)计算与之对应的虚拟旋转粒子的临时速度和位置；
步骤5、按公式(13)规则化处理旋转面粒子移动后的粒子数密度，求解压力泊松方程来计算旋转面粒子压力；
其中压力泊松方程的源项采用临时速度散度结合粒子数密度的相对变化的形式，如公式(28)所示，对公式(28)左侧的拉普拉斯算子通过公式(29)进行离散，
Figure FDA0003236532640000062
式中
Figure FDA0003236532640000063
——旋转面粒子i规则化处理后的粒子数密度；
γ——截断误差补偿系数，0.01≤γ≤0.05；
Figure FDA0003236532640000064
式中
Pi k+1——第k+1步旋转面粒子i的压力，Pa；
Figure FDA0003236532640000065
——第k+1步旋转面粒子i邻近的旋转面粒子j的压力，Pa；
Figure FDA0003236532640000066
——第k+1步旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a的压力，Pa；
Figure FDA0003236532640000071
——旋转面粒子i邻近的旋转面粒子j所对应的虚拟旋转粒子a移动后的位置矢量，m；
结合公式(28)和(29)，可得公式(30)：
Figure FDA0003236532640000072
式中
ai1——旋转面粒子i的第1个邻近粒子的压力系数；
ai2——旋转面粒子i的第2个邻近粒子的压力系数；
aii——旋转面粒子i的压力系数；
aiN+M——旋转面粒子i的第N+M个邻近粒子的压力系数；
P1 k+1——第k+1步旋转面粒子i的第1个邻近粒子的压力，Pa；
Figure FDA0003236532640000073
——第k+1步旋转面粒子i的第2个邻近粒子的压力，Pa；
Figure FDA0003236532640000074
——第k+1步旋转面粒子i的第N+M个邻近粒子的压力，Pa；
bi——旋转面粒子i的源项；
公式(30)中N为旋转面粒子总个数，M为虚拟旋转粒子总个数，其中：
Figure FDA0003236532640000075
Figure FDA0003236532640000076
式中
w(|rj'-ri|)——计算区域内任一粒子j'对旋转面粒子i的核函数值，表达形式如公式(10)；
将公式(30)变为矩阵表达式，如公式(33)所示：
Figure FDA0003236532640000081
由于虚拟旋转粒子的粒子数密度无法计算,导致公式(33)所述方程中未知变量个数多于方程个数，所以无法求解该压力泊松方程；因此，规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即
Figure FDA0003236532640000082
则压力泊松方程系数(公式(31))变为：
Figure FDA0003236532640000083
则压力泊松方程矩阵表达式(公式(33))变为:
Figure FDA0003236532640000084
由公式(35)可知，压力泊松方程的系数矩阵为非对称矩阵，所以采用BICG求解器对压力泊松方程进行求解；
步骤6、通过步骤5更新流场压力后，求解压力梯度，然后修正旋转面粒子速度和位置；
其中压力梯度按公式(36)计算：
Figure FDA0003236532640000091
式中
Figure FDA0003236532640000092
——旋转面粒子j所对应的虚拟旋转粒子k的压力，Pa；
Pj——旋转面粒子j的压力，Pa；
因规定虚拟旋转粒子的压力值等于所对应旋转面粒子的压力值，即
Figure FDA0003236532640000093
上式变形为：
Figure FDA0003236532640000094
修正后的旋转面粒子在Y轴方向的位移被重新定义为零；
步骤7、根据旋转面粒子修正后的速度和位置按照公式(1)-(6)修正与之对应的虚拟旋转粒子的速度和位置。
CN202111004067.9A 2021-08-30 2021-08-30 基于移动粒子半隐式法的轴对称流场二维模拟方法 Active CN113807034B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111004067.9A CN113807034B (zh) 2021-08-30 2021-08-30 基于移动粒子半隐式法的轴对称流场二维模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111004067.9A CN113807034B (zh) 2021-08-30 2021-08-30 基于移动粒子半隐式法的轴对称流场二维模拟方法

Publications (2)

Publication Number Publication Date
CN113807034A true CN113807034A (zh) 2021-12-17
CN113807034B CN113807034B (zh) 2023-05-16

Family

ID=78942274

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111004067.9A Active CN113807034B (zh) 2021-08-30 2021-08-30 基于移动粒子半隐式法的轴对称流场二维模拟方法

Country Status (1)

Country Link
CN (1) CN113807034B (zh)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100049489A1 (en) * 2008-08-22 2010-02-25 Kabushiki Kaisha Toshiba Flow simulation method, flow simulation system, and computer program product
CN102867094A (zh) * 2012-09-19 2013-01-09 西安交通大学 一种移动粒子半隐式算法中自由表面流动模型的构建方法
CN107633123A (zh) * 2017-09-13 2018-01-26 浙江工业大学 一种用于光滑粒子流体动力学模拟出血及处理加速的方法
CN109490955A (zh) * 2018-11-14 2019-03-19 深圳市勘察研究院有限公司 一种基于规则网格的声波波动方程正演模拟方法及装置
WO2020192126A1 (zh) * 2019-03-22 2020-10-01 大连理工大学 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112507600A (zh) * 2020-11-24 2021-03-16 西安交通大学 一种移动粒子半隐式法的对称边界条件的构建方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100049489A1 (en) * 2008-08-22 2010-02-25 Kabushiki Kaisha Toshiba Flow simulation method, flow simulation system, and computer program product
CN102867094A (zh) * 2012-09-19 2013-01-09 西安交通大学 一种移动粒子半隐式算法中自由表面流动模型的构建方法
CN107633123A (zh) * 2017-09-13 2018-01-26 浙江工业大学 一种用于光滑粒子流体动力学模拟出血及处理加速的方法
CN109490955A (zh) * 2018-11-14 2019-03-19 深圳市勘察研究院有限公司 一种基于规则网格的声波波动方程正演模拟方法及装置
WO2020192126A1 (zh) * 2019-03-22 2020-10-01 大连理工大学 一种基于改进的移动粒子半隐式法和模态叠加方法求解强非线性时域水弹性问题设计方法
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112507600A (zh) * 2020-11-24 2021-03-16 西安交通大学 一种移动粒子半隐式法的对称边界条件的构建方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
JANDAGHIAN, M.1等: "An enhanced weakly-compressible MPS method for free-surface flows", 《COMPUTER METHODS IN APPLIED MECHANICS & ENGINEERING》 *
LI, GEN等: "2D MPS analysis of hydrodynamic fine fragmentation of melt drop with initial steam film during fuel–coolant interaction", 《ANNALS OF NUCLEAR ENERGY》 *
傅宁琪等: "基于改进SPH方法的手术烟雾模拟研究", 《通化师范学院学报》 *
周子棋 等: "基于移动粒子半隐式方法波传播模型的声传播数值求解方法研究", 《西安交通大学学报》 *
宗潇 等: "矩形通道内蒸汽射流凝结换热面积的实验研究", 《西安交通大学学报》 *
段广涛等: "一种适用于气泡运动模拟的半隐式粒子法", 《工程热物理学报》 *

Also Published As

Publication number Publication date
CN113807034B (zh) 2023-05-16

Similar Documents

Publication Publication Date Title
Wang Adaptive high-order methods in computational fluid dynamics
Liu Integrated modeling of insect flight: from morphology, kinematics to aerodynamics
Yuki et al. Efficient immersed boundary method for strong interaction problem of arbitrary shape object with the self-induced flow
Trunk et al. Towards the simulation of arbitrarily shaped 3D particles using a homogenised lattice Boltzmann method
Magagnato KAPPA-Karlsruhe parallel program for aerodynamics
Zhao et al. A viscous vortex particle model for rotor wake and interference analysis
Sato et al. A new contact line treatment for a conservative level set method
Ding et al. An efficient dynamic mesh generation method for complex multi-block structured grid
CN111259325A (zh) 一种基于局部曲率自适应修正的改进水平集方法
Zhan et al. Numerical study on the six-DOF anchoring process of gravity anchor using a new mesh update strategy
Cheng et al. A novel finite point method for flow simulation
CN113807034A (zh) 基于移动粒子半隐式法的轴对称流场二维模拟方法
JP2009193110A (ja) グリッドフリー手法を用いた固気二相流シミュレーションプログラム及びそれを記憶した記憶媒体並びに固気二相流シミュレーション装置
Gouravaraju et al. Estimating the Hausdorff–Besicovitch dimension of boundary of basin of attraction in helicopter trim
Anandhanarayanan et al. Development and applications of a gridfree kinetic upwind solver to multi-body configurations
CN114692446A (zh) 一种矢量转动球坐标参数化的非线性壳有限元建模方法
CN114021499B (zh) 基于fvm-tlbfs方法的飞行器热防护结构热传导计算方法
Lofthouse et al. Static and dynamic simulations of a generic UCAV geometry using the kestrel flow solver
Haeri et al. Granular flow modeling of robot-terrain interactions in reduced gravity
Zhu et al. Modeling and aerodynamic characteristics analysis of morphing aircraft
Huang et al. A new grid deformation technology with high quality and robustness based on quaternion
CN105677983B (zh) 基于软硬件实时交互优化的计算方法
Dingman et al. Particle tracking in three‐dimensional Stokes flow
CN109741428A (zh) 一种适用于二维流体模拟的三阶高精度对流插值算法
Lacoursière et al. Fast and stable simulation of granular matter and machines

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant