CN116992796A - 一种自适应低耗散的sph-hllc黎曼求解器耦合算法 - Google Patents

一种自适应低耗散的sph-hllc黎曼求解器耦合算法 Download PDF

Info

Publication number
CN116992796A
CN116992796A CN202311267049.9A CN202311267049A CN116992796A CN 116992796 A CN116992796 A CN 116992796A CN 202311267049 A CN202311267049 A CN 202311267049A CN 116992796 A CN116992796 A CN 116992796A
Authority
CN
China
Prior art keywords
equation
particle
hllc
dissipation
density
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
CN202311267049.9A
Other languages
English (en)
Other versions
CN116992796B (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.)
University of Science and Technology of China USTC
Original Assignee
University of Science and Technology of China USTC
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 University of Science and Technology of China USTC filed Critical University of Science and Technology of China USTC
Priority to CN202311267049.9A priority Critical patent/CN116992796B/zh
Publication of CN116992796A publication Critical patent/CN116992796A/zh
Application granted granted Critical
Publication of CN116992796B publication Critical patent/CN116992796B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • 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
    • 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)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Fluid Mechanics (AREA)
  • Computing Systems (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种自适应低耗散的SPH‑HLLC黎曼求解器耦合算法，属于耦合算法技术领域。计算过程首先是读入初始生成的粒子坐标数据，根据位置将它们划分区域，大区域进行并行计算，大区域内的小区域作为存储粒子各种信息的单位区域，计算过程使用蛙跳算法：应用控制方程的动量方程计算出初始加速度并且将速度和坐标推进半个时间步，然后进行时间迭代循环；其中，控制方程中加入了耗散修正项。本发明采用上述一种自适应低耗散的SPH‑HLLC黎曼求解器耦合算法，构造简单，耗散较低，能够准确捕捉激波；并且既保留了黎曼算法的自适应性，避免人工粘性对于不同算例需要人工调整的问题，又保留了SPH构造格式的简单性以保持较高的计算效率和计算准确性。

Description

一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法
技术领域
本发明涉及耦合算法技术领域，尤其是涉及一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法。
背景技术
光滑粒子流体动力学，英文全称为Smoothed Particle Hydrodynamics，简称SPH。该方法是一种拉格朗日无网格粒子方法，用于模拟各种流体的流动。SPH方法由Lucy和Gingold等在1977年提出，一开始是服务于天体物理学领域，用来求解相关问题。在此方法中，材料由一系列可自由运动的粒子组成，任意粒子间无需网格连接，这些粒子具有材料的质量，密度，速度，压力等物理性质，并遵守守恒控制方程的规律运动。该方法可以模拟任意变形问题，由于采用拉格朗日方法描述流体，无对流项离散，理论上无数值耗散，且无需在网格上追踪界面，适合复杂界面跟踪、大变形等应用。
2000年，Parshikov提出粒子接触算法，此算法通过黎曼求解器求出接触间断处的应力以及速度的黎曼解，来分别替换原粒子间应力与速度的平均值，以此来表征粒子之间的相互作用，为复杂冲击界面演化模拟提供了新的解决思路，近年来已有较多的应用开展。
但是，SPH黎曼求解器耦合时会引入过多耗散导致整体结果计算准确度不高。由于在传统标准SPH算法中对于稀疏波的捕捉不太需要引入额外的人工粘性，因此引入等效耗散的黎曼求解器在对于稀疏波的处理中引入了过多的耗散导致计算准确度下降，而对于激波的捕捉黎曼求解器也需要进行一定的限制以提高计算准确性。这个问题可以通过与本质加权无震荡方法(WENO)等网格高精度方法结合来部分解决，但是由于SPH方法的拉格朗日粒子是自由运动的，构造与网格类似的对应于每一个粒子对的模板非常耗时，对于较大规模的物理计算而言太过浪费资源。
发明内容
本发明的目的是提供一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，构造简单，耗散较低，能够准确捕捉激波；并且既保留了黎曼算法的自适应性，避免人工粘性对于不同算例需要人工调整的问题，又保留了SPH构造格式的简单性以保持较高的计算效率和计算准确性。
为实现上述目的，本发明采用了如下技术方案：
一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，包括以下步骤：
S1、启动程序，读入初始生成的粒子坐标数据；
S2、根据位置将步骤S1中的数据划分为大区域、小区域、邻区域三种区域，大区域进行并行计算，大区域内的小区域作为存储粒子各种信息的单位区域；
S3、计算过程使用蛙跳算法：首先应用控制方程的动量方程计算出初始加速度并且将速度和坐标推进半个时间步，然后进行时间迭代循环；
其中，控制方程中加入了耗散修正项，控制方程包括密度更新方程、连续性方程、动量方程、能量方程、状态方程。
优选的，步骤S3中时间迭代循环过程如下：
S31、划分大区域、小区域、邻区域三种区域；
S32、用控制方程中的连续性方程或密度更新方程更新粒子的密度，用能量方程更新粒子的内能，同时计算新光滑长度并储存，更新一个时间步；
S33、推进粒子位移半个时间步并且重新划分区域；
S34、计算状态方程得到压力值；
S35、通过得到的压力信息用动量方程得到加速度，将速度推进半个时间步；
S36、在满足输出条件后输出所有粒子的数据，记时间步+1；
S37、最后推进速度和位移半个时间步，并更新先前储存的光滑长度。
优选的，用标记界面，支持域内有异相粒子/>标记为1，同相粒子/>标记为0，同相粒子使用密度求和法更新密度，异相粒子使用连续性方程更新密度，
密度更新方程，通过支持域内粒子质量与其权重的乘积累加直接得到更新后的密度：
；
其中，表示粒子/>的密度，/>表示粒子/>的密度，/>表示粒子对/>对应的核函数。
连续性方程，通过速度的散度乘密度得到一个时间步密度变化量，表示密度的耗散修正项，加在连续性方程的末尾：
；
；
其中，下标表示对应黎曼解，上标/>表示该黎曼解是用HLLC黎曼求解器；/>为使用HLLC黎曼求解器解得的速度解，/>为粒子i的速度，/>为粒子i的空间求导，/>表示粒子j的体积，/>为自定义参量并设置为0.2，/>为粒子对/>的界面参量，/>为粒子对/>的声速，/>为粒子对/>的光滑长度，/>为粒子对/>密度的中间参量。
优选的，动量方程，通过压强来计算出粒子的加速度，/>表示加速度的耗散修正项，加在动量方程的末尾：
；
；
其中，表示使用HLLC黎曼求解器解得的压强解，/>为/>两粒子的Roe平均密度，为粒子对/>速度的中间参量。
优选的，能量方程，通过压强和速度的信息计算出粒子内能一个时间步的变化，/>表示内能的耗散修正项，加在能量方程的末尾：
；
；
其中，为粒子对/>内能的中间参量，/>为粒子j的速度，/>为自定义参量并设置为0.02。
优选的，状态方程，通过密度和内能来更新压强的大小：
。
其中，为压强，/>表示密度，/>表示粒子内能，/>表示气体常数随气体不同而变化。
优选的，在步骤S3的控制方程中使用了限制因子以限制HLLC黎曼求解器引入的耗散，重构后的黎曼解公式为：
；
将修正后得到的替代HLLC-SPH控制方程中的/>或/>进行计算；
其中，为使用重构后HLLC黎曼求解器解得的速度或压强，/>表示新算法赋予的平衡传统HLLC黎曼算法的权参数，/>为使用重构前HLLC黎曼求解器解得的速度或压强，/>为物理量速度或者压强的算术平均。
优选的，参数的计算公式为：
；
；
；
；
；
；
；
；
其中，、/>为中间变量参数，下标/>为压强的参数，下标/>为速度的参数，下标/>、/>为粒子/>或/>在局部坐标系下标记为/>或/>，下标/>或/>为针对压强黎曼解处于激波或稀疏波状态的处理，下标g为对于物理性质的处理。
优选的，表示有关密度内能和速度的中间参量，/>表示界面的参量，下标/>表示该参量是属于/>两个粒子组合的粒子对的变量：
；
；
；
。
优选的，步骤S3中，光滑长度通过光滑长度平滑过渡进行更新：
；
核函数内代入的光滑长度使用上一时间步的光滑长度信息；
其中，为光滑长度的无量纲参数，上标/>为计算时刻的时间步数，/>表示第/>个时间步下粒子/>的光滑长度，/>表示第/>个时间步下粒子/>的光滑长度，/>为计算的维度。
因此，本发明采用上述结构的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，实现的有益效果为：
1、首先是同相粒子使用密度求和法，异相粒子使用连续性方程更新密度。因为使用连续性方程更新密度总会在密度强烈变化时低估密度，而密度求和方法可以将密度的误差限制在一个有限的范围内不会造成误差累积导致计算失效。对于界面附近的异相粒子，使用密度求和方法更新密度会导致压强剧烈非物理震荡，因此同相粒子和异相粒子作更新密度的区分很好地兼顾了密度误差限制和界面压强震荡的弱点。同样，对于界面附近的光滑长度使用空间分布信息光滑函数插值得到，可以有效防止密度大变化时界面两侧或者接触简短处光滑长度相差过大的缺陷；
2、黎曼求解器引入了黎曼解的独特限制，重构了速度和压强的黎曼解。在参数中，用于防止接触间断处过大的压强震荡，此时如果界面被黎曼解识别为近似接触间断，则/>接近为1，使得重构黎曼解接近原先的HLLC黎曼解。/>的值取决于黎曼解的状态，如果显示为双稀疏波状态则/>的值较小，意味着引入了较小的耗散，趋近于无耗散的传统解，如果为双激波状态/>则接近为1，重构解接近原HLLC黎曼解；
3、对于接触间断除了上面提到的密度更新和光滑长度更新之外，还有弱耗散的引入以抹平界面震荡。它采用了常见的耗散引入形式，但是用了上面识别到的粒子同相异相状态来控制连续性方程、动量方程和能量方程的引入耗散量。如果粒子对异相，则不引入耗散以保证界面清脆间断；如果粒子对同相且相同，引入少量耗散稳定计算；如果粒子对同相但/>不同，引入相对上述更少的耗散以使得用密度求和法和连续性方程更新密度的区域之间平滑过渡。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1为本发明一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法的流程图；
图2为本发明一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法实施例2中的一个二维激波管算例图；（a）为直径-密度结果图，（b）为直径-压力结果图；
图3为本发明一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法实施例2中的另一个二维激波管算例图；（a）为直径-密度结果图，（b）为直径-压力结果图；
图4为本发明一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法实施例3中的空气六氟化硫的界面演化图；（a）为初始设置概念简图，右侧（b）为实际建模演示图；
图5为本发明一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法实施例4中的不同微秒时刻的计算结果图；（a）为传统SPH算法结果图，（b）为ROE黎曼求解器耦合SPH算法结果图，（c）为本发明的新算法的结果图。
具体实施方式
以下通过附图和实施例对本发明的技术方案作进一步说明。
除非另外定义，本发明使用的技术术语或者科学术语应当为本发明所属领域内具有一般技能的人士所理解的通常意义。本发明中使用的“第一”、“第二”以及类似的词语并不表示任何顺序、数量或者重要性，而只是用来区分不同的组成部分。“包括”或者“包含”等类似的词语意指出现该词前面的元件或者物件涵盖出现在该词后面列举的元件或者物件及其等同，而不排除其他元件或者物件。术语“设置”、“安装”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。“上”、“下”、“左”、“右”等仅用于表示相对位置关系，当被描述对象的绝对位置改变后，则该相对位置关系也可能相应地改变。
实施例1
如图1所示，本发明提供了一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，既保留了黎曼算法的自适应性避免人工粘性对于不同算例需要人工调整的问题，又保留SPH构造格式的简单性以保持较高的计算效率。
一、包括以下步骤：
S1、启动程序，读入初始生成的粒子坐标数据；
S2、根据位置将步骤S1中的数据划分为大区域、小区域、邻区域三种区域，大区域进行并行计算，大区域内的小区域作为存储粒子各种信息的单位区域；
S3、计算过程使用蛙跳算法：首先应用控制方程的动量方程计算出初始加速度并且将速度和坐标推进半个时间步，然后进行时间迭代循环：
；
；
循环过程如下：
S31、划分大区域、小区域、邻区域三种区域；
循环中划分区域时对于跑到另一个区域的粒子的处理方法是在本区域内删除这个粒子信息然后添加到跑到的那个区域，而对于还在同一个小区域的粒子就无需做处理。
S32、用控制方程中的连续性方程或密度更新方程更新粒子的密度，用能量方程更新粒子的内能，同相粒子使用密度求和法更新密度，异相粒子使用连续性方程更新密度，同时计算新光滑长度并储存，更新一个时间步，公式如下：
；
。
由于强可压缩现象中出现巨大的粒子支持域半径间断现象，对于光滑长度的更新使用如下更新方式来实现光滑长度平滑过渡：
；
核函数内代入的光滑长度使用上一时间步的光滑长度信息；
其中，为光滑长度的无量纲参数（一般根据不同工况采用1.2~1.4之间的值），上标/>为计算时刻的时间步数，/>表示第/>个时间步下粒子/>的光滑长度，/>表示第/>个时间步下粒子/>的光滑长度，/>为计算的维度；该公式通过使用空间分布信息光滑函数插值得到跨越界面较为连续的光滑长度，在计算控制方程内的密度和能量时同时计算。
S33、推进粒子位移半个时间步并且重新划分区域，公式如下：
；
S34、计算状态方程得到压力值，公式如下：
；
S35、通过得到的压力信息用动量方程得到加速度，将速度推进半个时间步，公式如下：
；
S36、在满足某输出条件（例如每规定时间步数或者每规定物理时间，视不同工况要求和数据分析而定）后输出所有粒子的数据，记时间步+1；
S37、最后推进速度和位移半个时间步，并更新先前储存的光滑长度：
；
二、步骤S3中的控制方程中加入了耗散修正项，包括：
（1）密度更新方程：
；
是迭代密度的一种方式，非界面的单相区域即使用密度求和法求解，通过支持域内粒子质量与其权重的乘积累加直接得到更新后的密度。
其中，表示粒子/>的密度，/>表示粒子/>的密度，/>表示粒子对/>应的核函数。
（2）连续性方程：
；
；
是迭代密度的另一种方式，在界面附近即的粒子密度使用连续性方程求解，通过速度的散度乘密度得到一个时间步密度变化量/>，/>表示密度的耗散修正项，加在连续性方程的末尾。
其中，下标表示对应黎曼解，上标/>表示该黎曼解是用HLLC黎曼求解器；/>为使用HLLC黎曼求解器解得的速度解，/>为粒子i的速度，/>为粒子i的空间求导，/>表示粒子j的体积，/>为自定义参量并设置为0.2，/>为粒子对/>的界面参量，/>为粒子对/>的声速，/>为粒子对/>的光滑长度，/>为粒子对/>密度的中间参量。
（3）动量方程：
；
；
通过压强来计算出粒子的加速度，即一个时间步内速度的变化量，/>表示加速度的耗散修正项，加在动量方程的末尾。
其中，表示使用HLLC黎曼求解器解得的压强解，/>为/>两粒子的Roe平均密度，为粒子对/>速度的中间参量。
（4）能量方程：
；
；
通过压强和速度的信息计算出粒子内能一个时间步的变化，/>表示内能的耗散修正项，加在能量方程的末尾。
其中，为粒子对/>内能的中间参量，/>为粒子j的速度，/>为自定义参量并设置为0.02。
（5）状态方程：
；
通过密度和内能来更新压强的大小，状态方程根据物质的不同可以有很多表达形式，并不一定局限于方程里面那一种形式。
其中，用标记界面，支持域内有异相粒子/>标记为1，同相粒子/>标记为0，/>为压强，/>表示密度，/>表示粒子内能，/>表示气体常数随气体不同而变化。
上述方程中有关密度内能和速度的中间参量以及/>的界面的参量由下式给出，下标/>表示该参量是属于/>两个粒子组合的粒子对的变量：
；
；
；
。
三、在步骤S3的控制方程中使用了限制因子以限制HLLC黎曼求解器引入的耗散，重构后的黎曼解公式为：
；
将修正后得到的替代HLLC-SPH控制方程中的/>或/>进行计算；
其中，为使用重构后HLLC黎曼求解器解得的速度或压强，/>表示新算法赋予的平衡传统HLLC黎曼算法的权参数，/>为使用重构前HLLC黎曼求解器解得的速度或压强，/>为物理量速度或者压强的算术平均。
参数的计算公式为：
；
；
；
；
；
；
；
；
其中，、/>为中间变量参数，下标/>为压强的参数，下标/>为速度的参数，下标/>、/>为粒子/>或/>在局部坐标系下标记为/>或/>，下标/>或/>为针对压强黎曼解处于激波或稀疏波状态的处理，下标g为对于物理性质的处理。
实施例2
本发明提供了一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，以两个二维激波管为算例，用于验证本发明有效性。
计算域是一个半径为1的圆面，在半径0.5以内的区域内置高压气体，外侧圆环放置低压气体，其初始条件分别为：
；
如图2、图3所示的两个激波管测试算例沿直径截面的计算结果，可以看到我们的新算法的结果相对于传统算法有极大的优越性，对于稀疏波和激波的捕捉更加精确，而对于界面附近的震荡也有极强的削弱效果。
实施例3
本发明还用了空气六氟化硫界面在圆柱形汇聚激波冲击下的界面演化问题来进一步验证方法优越性。
如图4所示，在半径为300mm的二维圆面区域里面，70mm外侧充满了高压气体，内侧是低压气体，而最里面有径长20mm的八边形六氟化硫气体。图4中的图（a）为初始设置概念简图，图4中的图（b）为实际建模演示图。
实施例4
如图5所示，该算例分别在80、160、240、320、400、480、560、640微秒时刻的计算结果，图5中图（a）、（b）、（c）第一行的右半部分为实验结果。图5中（a）为传统SPH算法结果图，图5中（b）为ROE黎曼求解器耦合SPH算法结果图，图5中（c）为本发明的新算法结果图。
可以看到在演化前期所有算法显示的界面演化都与实验差别不大，而到了后期400微妙之后的时刻里面，新算法明显比另外两个算法展现了更强的捕捉界面演化微小结构的能力。
因此，本发明采用上述一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，构造简单，耗散较低，能够准确捕捉激波；并且既保留了黎曼算法的自适应性，避免人工粘性对于不同算例需要人工调整的问题，又保留了SPH构造格式的简单性以保持较高的计算效率和计算准确性。
最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。

Claims (9)

1.一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于，包括以下步骤：
S1、启动程序，读入初始生成的粒子坐标数据；
S2、根据位置将步骤S1中的数据划分为大区域、小区域、邻区域三种区域，大区域进行并行计算，大区域内的小区域作为存储粒子各种信息的单位区域；
S3、计算过程使用蛙跳算法：首先应用控制方程的动量方程计算出初始加速度并且将速度和坐标推进半个时间步，然后进行时间迭代循环；
其中，控制方程中加入了耗散修正项，控制方程包括密度更新方程、连续性方程、动量方程、能量方程、状态方程。
2.根据权利要求1所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：步骤S3中时间迭代循环过程如下：
S31、划分大区域、小区域、邻区域三种区域；
S32、用控制方程中的连续性方程或密度更新方程更新粒子的密度，用能量方程更新粒子的内能，同时计算新光滑长度并储存，更新一个时间步；
S33、推进粒子位移半个时间步并且重新划分区域；
S34、计算状态方程得到压力值；
S35、通过得到的压力信息用动量方程得到加速度，将速度推进半个时间步；
S36、在满足输出条件后输出所有粒子的数据，记时间步+1；
S37、最后推进速度和位移半个时间步，并更新先前储存的光滑长度。
3.根据权利要求2所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：用标记界面，支持域内有异相粒子/>标记为1，同相粒子/>标记为0，同相粒子使用密度求和法更新密度，异相粒子使用连续性方程更新密度；
密度更新方程，通过支持域内粒子质量与其权重的乘积累加直接得到更新后的密度：
；
其中，表示粒子/>的密度，/>表示粒子/>的密度，/>表示粒子对/>对应的核函数；
连续性方程，通过速度的散度乘密度得到一个时间步密度变化量，/>表示密度的耗散修正项，加在连续性方程的末尾：
；
；
其中，下标表示对应黎曼解，上标/>表示该黎曼解是用HLLC黎曼求解器；/>为使用HLLC黎曼求解器解得的速度解，/>为粒子i的速度，/>为粒子i的空间求导，/>表示粒子j的体积，/>为自定义参量并设置为0.2，/>为粒子对/>的界面参量，/>为粒子对/>的声速，/>为粒子对/>的光滑长度，/>为粒子对/>密度的中间参量。
4.根据权利要求3所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：
动量方程，通过压强来计算出粒子的加速度，/>表示加速度的耗散修正项，加在动量方程的末尾：
；
；
其中，表示使用HLLC黎曼求解器解得的压强解，/>为/>两粒子的Roe平均密度，/>为粒子对/>速度的中间参量。
5.根据权利要求4所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：
能量方程，通过压强和速度的信息计算出粒子内能一个时间步的变化，/>表示内能的耗散修正项，加在能量方程的末尾：
；
；
其中，为粒子对/>内能的中间参量，/>为粒子j的速度，/>为自定义参量并设置为0.02。
6.根据权利要求5所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：
状态方程，通过密度和内能来更新压强的大小：
；
其中，为压强，/>表示密度，/>表示粒子内能，/>表示气体常数随气体不同而变化。
7.根据权利要求6所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：在步骤S3的控制方程中使用了限制因子以限制HLLC黎曼求解器引入的耗散，重构后的黎曼解公式为：
；
将修正后得到的替代HLLC-SPH控制方程中的/>或/>进行计算；
其中，为使用重构后HLLC黎曼求解器解得的速度或压强，/>表示新算法赋予的平衡传统HLLC黎曼算法的权参数，/>为使用重构前HLLC黎曼求解器解得的速度或压强，/>为物理量速度或者压强的算术平均。
8.根据权利要求7所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：参数的计算公式为：
；
；
；
；
；
；
；
；
其中，、/>为中间变量参数，下标/>为压强的参数，下标/>为速度的参数，下标/>、/>为粒子/>或/>在局部坐标系下标记为/>或/>，下标/>或/>为针对压强黎曼解处于激波或稀疏波状态的处理，下标/>为对于物理性质的处理。
9.根据权利要求8所述的一种自适应低耗散的SPH-HLLC黎曼求解器耦合算法，其特征在于：步骤S3中，光滑长度通过光滑长度平滑过渡进行更新：
；
核函数内代入的光滑长度使用上一时间步的光滑长度信息；
其中，为光滑长度的无量纲参数，上标/>为计算时刻的时间步数，/>表示第/>个时间步下粒子/>的光滑长度，/>表示第/>个时间步下粒子/>的光滑长度，/>为计算的维度。
CN202311267049.9A 2023-09-27 2023-09-27 一种自适应低耗散的sph-hllc黎曼求解器耦合方法 Active CN116992796B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202311267049.9A CN116992796B (zh) 2023-09-27 2023-09-27 一种自适应低耗散的sph-hllc黎曼求解器耦合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202311267049.9A CN116992796B (zh) 2023-09-27 2023-09-27 一种自适应低耗散的sph-hllc黎曼求解器耦合方法

Publications (2)

Publication Number Publication Date
CN116992796A true CN116992796A (zh) 2023-11-03
CN116992796B CN116992796B (zh) 2023-12-22

Family

ID=88534334

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202311267049.9A Active CN116992796B (zh) 2023-09-27 2023-09-27 一种自适应低耗散的sph-hllc黎曼求解器耦合方法

Country Status (1)

Country Link
CN (1) CN116992796B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117951973A (zh) * 2024-03-26 2024-04-30 中国科学技术大学 一种稳定的弹塑性实体光滑粒子动力学模拟方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012243288A (ja) * 2011-05-24 2012-12-10 Fujitsu Ltd シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
KR20120137827A (ko) * 2011-06-13 2012-12-24 한국과학기술원 Ｓｐｈ 유체를 위한 서브 입자 스케일 난류 시뮬레이션 방법, 시스템 및 이를 위한 기록 매체
US20130338977A1 (en) * 2011-02-15 2013-12-19 Fujitsu Limited Simulation device and simulation method
US20140365185A1 (en) * 2013-06-11 2014-12-11 Fujitsu Limited Numerical calculation method and apparatus
US20160004802A1 (en) * 2014-07-03 2016-01-07 Arizona Board Of Regents On Behalf Of Arizona State University Multiscale Modelling of Growth and Deposition Processes in Fluid Flow
CN109214082A (zh) * 2018-09-03 2019-01-15 北京理工大学 一种近场水下爆炸冲击波载荷的高精度数值模拟方法
CN111709197A (zh) * 2020-06-17 2020-09-25 福州大学 基于黎曼不变量的sph入流边界处理方法
CN112861445A (zh) * 2020-12-29 2021-05-28 中国船舶重工集团公司第七研究院 一种无网格数值模型的模拟方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130338977A1 (en) * 2011-02-15 2013-12-19 Fujitsu Limited Simulation device and simulation method
JP2012243288A (ja) * 2011-05-24 2012-12-10 Fujitsu Ltd シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
KR20120137827A (ko) * 2011-06-13 2012-12-24 한국과학기술원 Ｓｐｈ 유체를 위한 서브 입자 스케일 난류 시뮬레이션 방법, 시스템 및 이를 위한 기록 매체
US20140365185A1 (en) * 2013-06-11 2014-12-11 Fujitsu Limited Numerical calculation method and apparatus
US20160004802A1 (en) * 2014-07-03 2016-01-07 Arizona Board Of Regents On Behalf Of Arizona State University Multiscale Modelling of Growth and Deposition Processes in Fluid Flow
CN109214082A (zh) * 2018-09-03 2019-01-15 北京理工大学 一种近场水下爆炸冲击波载荷的高精度数值模拟方法
CN111709197A (zh) * 2020-06-17 2020-09-25 福州大学 基于黎曼不变量的sph入流边界处理方法
CN112861445A (zh) * 2020-12-29 2021-05-28 中国船舶重工集团公司第七研究院 一种无网格数值模型的模拟方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
QIUZU YANG ET AL.: "A multi-phase SPH model based on Riemann solvers for simulation of jet breakup", ENGINEERING ANALYSIS WITH BOUNDARY ELEMENTS 111 (2020), pages 134 - 147 *
徐建于等: "圆柱形汇聚激波诱导Richtmyer-Meshkov不稳定的SPH模拟", 力学学报, pages 998 - 1011 *
徐志宏;汤文辉;张若棋;: "基于黎曼解的粒子间接触算法在SPH中的应用", 高压物理学报, no. 01, pages 93 - 96 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117951973A (zh) * 2024-03-26 2024-04-30 中国科学技术大学 一种稳定的弹塑性实体光滑粒子动力学模拟方法
CN117951973B (zh) * 2024-03-26 2024-06-07 中国科学技术大学 一种稳定的弹塑性实体光滑粒子动力学模拟方法

Also Published As

Publication number Publication date
CN116992796B (zh) 2023-12-22

Similar Documents

Publication Publication Date Title
CN116992796B (zh) 一种自适应低耗散的sph-hllc黎曼求解器耦合方法
Macklin et al. Unified particle physics for real-time applications
Feng et al. Layered infill area generation from triply periodic minimal surfaces for additive manufacturing
JP5045853B2 (ja) 計算用データ生成装置、計算用データ生成方法及び計算用データ生成プログラム
Li et al. A convergent and accuracy preserving limiter for the FR/CPR method
Wang et al. A triangular grid generation and optimization framework for the design of free-form gridshells
JP2008152423A (ja) 粒子モデルを用いた変形挙動シミュレーション方法およびプログラム
Isola et al. Finite-volume solution of two-dimensional compressible flows over dynamic adaptive grids
Gerace et al. A model-integrated localized collocation meshless method for large scale three-dimensional heat transfer problems
Zhan et al. Meshfree method based on discrete gas-kinetic scheme to simulate incompressible/compressible flows
Tang et al. Honey, I Shrunk the Domain: Frequency‐aware Force Field Reduction for Efficient Fluids Optimization
Ciallella et al. Extrapolated DIscontinuity Tracking for complex 2D shock interactions
Esser et al. An extended finite element method applied to levitated droplet problems
Tang et al. A harmonic balance approach for designing compliant mechanical systems with nonlinear periodic motions
Zhao et al. An arbitrary Lagrangian-Eulerian RKDG method for multi-material flows on adaptive unstructured meshes
Huang et al. Physically-based smoke simulation for computer graphics: a survey
Wang et al. Axisymmetric Riemann–smoothed particle hydrodynamics modeling of high-pressure bubble dynamics with a simple shifting scheme
Lu et al. A Rigging‐Skinning Scheme to Control Fluid Simulation
CN111801672A (zh) 用于设计和制造稳定的网格结构的系统和方法
Qi et al. Acceleration strategies based on an improved bubble packing method
Miranda et al. A three-dimensional adaptive mesh generation approach using geometric modeling with multi-regions and parametric surfaces
Azevedo et al. Efficient smoke simulation on curvilinear grids
Rozenberg et al. Fluid structure interaction problems in turbomachinery using rbf interpolation and greedy algorithm
Altenhofen et al. Implicit Mesh Generation using Volumetric Subdivision.
Zeng et al. Self-adaptive triangular mesh generation framework for free-form single-layer reticulated shells based on virtual interaction forces

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