CN116384288B - 一种可压缩流动高分辨率数值模拟方法、介质及设备 - Google Patents

一种可压缩流动高分辨率数值模拟方法、介质及设备 Download PDF

Info

Publication number
CN116384288B
CN116384288B CN202310654382.9A CN202310654382A CN116384288B CN 116384288 B CN116384288 B CN 116384288B CN 202310654382 A CN202310654382 A CN 202310654382A CN 116384288 B CN116384288 B CN 116384288B
Authority
CN
China
Prior art keywords
flux
point
grid
variable
time
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.)
Active
Application number
CN202310654382.9A
Other languages
English (en)
Other versions
CN116384288A (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.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
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 Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202310654382.9A priority Critical patent/CN116384288B/zh
Publication of CN116384288A publication Critical patent/CN116384288A/zh
Application granted granted Critical
Publication of CN116384288B publication Critical patent/CN116384288B/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/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • 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

Abstract

本申请涉及计算流体力学中数值计算方法领域，公开了一种可压缩流动高分辨率数值模拟方法、介质及设备，包括：读取初始流场数据，得到当前时刻计算空间中各网格点的变量和通量；根据得到的变量，采用非线性插值计算半点处的左值和右值，并根据左值和右值，计算半点处的低阶通量；根据得到的通量，计算半点处的修正通量；根据低阶通量和修正通量，获取计算空间中不同方向各节点处的一阶导数，以完成空间离散；采用龙格‑库塔法将当前时刻的变量推进到下一时刻的变量，以完成时空离散；将时间推进至设定时刻结束计算，得到设定时刻的流场数据。这样可以得到合理的可压缩流动数值模拟结果，稳定捕捉相应的可压缩流动，显著提高对流动结构的分辨能力。

Description

一种可压缩流动高分辨率数值模拟方法、介质及设备
技术领域
本发明涉及计算流体力学中数值计算方法领域，特别是涉及一种可压缩流动高分辨率数值模拟方法、介质及设备。
背景技术
可压缩流动在高速列车、空天飞行器等绕流中广泛存在。欧拉方程是可压缩流动在无粘假设条件下的控制方程。有限差分法能以较小计算代价实现多维问题的高阶精度计算，是可压缩流动数值模拟的主要方法之一。
针对复杂几何外形绕流问题的可压缩流动数值模拟时，采用有限差分方法模拟主要面临两个问题：一是复杂几何外形的网格会存在局部不光滑，在网格不光滑处坐标变换时所引入的几何诱导误差，会带来分辨能力的降低，严重时甚至会导致计算得到的流动结构失真；二是工程计算中无粘通量离散通常采用总模板为四点的总变差减小类格式或加权本质无振荡格式，该类数值格式虽然具有良好的数值鲁棒性，但存在数值耗散误差过大，这将会带来流动结构分辨能力的降低。
因此，如何提高对流动结构的分辨能力，是本领域技术人员亟待解决的技术问题。
发明内容
有鉴于此，本发明的目的在于提供一种可压缩流动高分辨率数值模拟方法、介质及设备，可以显著提高对流动结构的分辨能力。其具体方案如下：
一种可压缩流动高分辨率数值模拟方法，包括：
读取初始流场数据，得到当前时刻计算空间中各网格点的变量和通量；
根据得到的所述各网格点的变量，采用非线性插值计算半点处的左值和右值，并根据所述左值和右值，计算半点处的低阶通量；
根据得到的所述各网格点的通量，计算半点处的修正通量；
根据计算出的所述低阶通量和所述修正通量，获取所述计算空间中不同方向各节点处的一阶导数，以完成空间离散；
采用龙格-库塔法将当前时刻的变量推进到下一时刻的变量，以完成时空离散；
将时间推进至设定时刻结束计算，得到所述设定时刻的流场数据。
优选地，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，读取初始流场数据，得到当前时刻计算空间中各网格点的变量和通量，包括：
读取t 0时刻流场数据，采用四阶中心差分离散方法计算所述t 0时刻流场数据对应的网格度量和雅克比，并将t n 时刻物理空间中各网格点的变量/>和通量/>进行转换，得到t n 时刻计算空间/>中各网格点的变量/>和通量/>。
优选地，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，采用四阶中心差分离散方法计算所述t 0时刻流场数据对应的网格度量和雅克比，并将t n 时刻物理空间中各网格点的变量/>和通量/>进行转换，得到t n 时刻计算空间中各网格点的变量/>和通量/>，包括：
在无粘假设条件下，将可压缩流动的控制方程设定为欧拉方程，在笛卡尔坐标系下，无量纲形式三维守恒型的欧拉方程为：
；
其中，为守恒变量，/>分别为直角坐标系x,y,z方向的无粘通量；
在固定网格下，进行坐标变换：
；
得到计算坐标下的欧拉方程为：
；
其中，
；
其中雅克比和网格度量公式采用对称守恒计算形式：
；
；
；
；
其中，、/>、/>、/>、/>、/>、/>、/>、/>分别为一阶偏导数；
在所述计算空间进行等距离散，若以f表示为某一特定物理变量，以方向为例，上述公式中的一阶偏导数离散方式，采用下述公式离散：
；
其中，下标i为网格节点指标、i+1/2为离散后半点上指标、h为网格间距；
半点处i+1/2的物理量f的值由四阶中心差分插值得到：
；
根据计算得到的网格度量和雅克比，得到计算空间中各网格点上的守恒变量和通量/>。
优选地，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，根据得到的所述各网格点的变量，采用非线性插值计算半点处的左值和右值，包括：
考虑方向，采用四阶中心型非线性插值的方法将节点处物理量/>插值到半点处；半点处的左值记为/>，半点处的右值记为/>：
；
；
；
其中非线性权值求解如下：
；
；
；
其中全局光滑度量为：
；
当地光滑度量为：
；
。
优选地，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，采用下述公式计算半点处的修正通量：
；
其中，为半点处的修正通量，/>为/>方向网格间距，/>为/>在半点i+1/2处/>方向的二阶偏导数。
优选地，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，根据计算出的所述低阶通量和所述修正通量，获取方向各节点处的一阶导数，包括：
将计算出的所述低阶通量和所述修正通量/>求和，得到半点处的高阶通量/>：
；
采用下述公式计算方向各节点处一阶导数/>：
。
优选地，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，采用下述公式计算方向和/>方向每个节点处一阶导数/>和/>：
；
；
欧拉方程的半离散形式为：
。
优选地，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，采用龙格-库塔法将当前时刻的变量推进到下一时刻的变量，包括：
采用三阶龙格-库塔法进行时间推进：
；
其中，上标n为第n时刻步的值，上标n+1为第n+1时刻步的值，RHS为右端项格式空间离散的值。
本发明实施例还提供了一种计算机可读存储介质，用于存储计算机程序，其中，所述计算机程序被处理器执行时实现如本发明实施例提供的上述可压缩流动高分辨率数值模拟方法。
本发明实施例还提供了一种可压缩流动高分辨率数值模拟设备，包括处理器和存储器，其中，所述处理器执行所述存储器中存储的计算机程序时实现如本发明实施例提供的上述可压缩流动高分辨率数值模拟方法。
从上述技术方案可以看出，本发明所提供的一种可压缩流动高分辨率数值模拟方法，包括：读取初始流场数据，得到当前时刻计算空间中各网格点的变量和通量；根据得到的各网格点的变量，采用非线性插值计算半点处的左值和右值，并根据左值和右值，计算半点处的低阶通量；根据得到的各网格点的通量，计算半点处的修正通量；根据计算出的低阶通量和修正通量，获取计算空间中不同方向各节点处的一阶导数，以完成空间离散；采用龙格-库塔法将当前时刻的变量推进到下一时刻的变量，以完成时空离散；将时间推进至设定时刻结束计算，得到设定时刻的流场数据。
本发明提供的上述可压缩流动高分辨率数值模拟方法，通过执行上述步骤，利用欧拉方程作为模型方程进行逐步数值离散，可以得到合理的可压缩流动的数值模拟结果，稳定地捕捉相应的可压缩流动，能够解决现有技术中畸变网格适应性差且数值耗散较大的缺点，显著提高对流动结构的分辨能力。
此外，本发明还针对可压缩流动高分辨率数值模拟方法提供了相应的计算机可读存储介质及设备，进一步使得上述方法更具有实用性，该计算机可读存储介质及设备具有相应的优点。
附图说明
为了更清楚地说明本发明实施例或相关技术中的技术方案，下面将对实施例或相关技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1为本发明实施例提供的可压缩流动高分辨率数值模拟方法的流程图；
图2为本发明实施例提供的可压缩流动高分辨率数值模拟方法的具体流程图；
图3a为采用经典三阶WENO的方法在规则网格下等熵涡传播问题的密度分布图；
图3b为采用本发明实施例提供的方法在规则网格下等熵涡传播问题的密度分布图；
图4a为采用经典三阶WENO的方法在随机网格下等熵涡传播问题的密度分布图；
图4b为采用本发明实施例提供的方法在随机网格下等熵涡传播问题的密度分布图；
图5为本发明实施例提供的双马赫反射问题的计算网格和密度分布图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明提供一种可压缩流动高分辨率数值模拟方法，如图1所示，包括以下步骤：
S101、读取初始流场数据，得到当前时刻计算空间中各网格点的变量和通量；
在具体实施时，步骤S101可以包括：读取t 0时刻流场数据作为初始流场数据，采用四阶中心差分离散方法计算t 0时刻流场数据对应的网格度量和雅克比，并将t n 时刻物理空间中各网格点的变量/>和通量/>进行转换，得到t n 时刻计算空间/>中各网格点的变量/>和通量/>。
S102、根据得到的各网格点的变量，采用非线性插值计算半点处的左值和右值，并根据左值和右值，计算半点处的低阶通量；
具体地，可以先考虑方向，根据步骤S101得到的当前时刻的各个网格点物理量，采用非线性插值计算得到半点处的左值/>和右值/>；根据得到的半点处左值/>和右值/>，计算半点处的低阶黎曼通量。
S103、根据得到的各网格点的通量，计算半点处的修正通量；
具体地，根据步骤S101得到的各网格点的，计算半点处的修正通量。
S104、根据计算出的低阶通量和修正通量，获取计算空间中不同方向各节点处的一阶导数，以完成空间离散；
具体地，根据步骤S102得到的低阶通量和步骤S103得到的修正通量计算半点处的高阶精度通量，并计算/>方向各节点处的一阶导数/>。
针对方向和/>方向的一阶导数/>和/>，重复/>方向步骤S102至步骤S104的操作，完成方程空间离散。
S105、采用龙格-库塔法将当前时刻的变量推进到下一时刻的变量，以完成时空离散；
具体地，采用龙格-库塔法将第n时刻推进到第n+1时刻，完成方程的时空离散。
S106、将时间推进至设定时刻结束计算，得到设定时刻的流场数据。
需要说明的是，本发明实施例提供的上述可压缩流动高分辨率数值模拟方法是面向可压缩流动，可压缩流动在高速列车、空天飞行器等绕流中广泛存在，这些都可以是应用场景。
在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，通过执行上述步骤，利用欧拉方程作为模型方程进行逐步数值离散，可以得到合理的可压缩流动的数值模拟结果，稳定地捕捉相应的可压缩流动，能够解决现有技术中畸变网格适应性差且数值耗散较大的缺点，显著提高对流动结构的分辨能力。
在具体实施时，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，在执行步骤S101时，具体可以包括以下步骤：
可压缩流动的控制方程无粘假设条件下可考虑为欧拉方程。在笛卡尔坐标系下，无量纲形式三维守恒型如下：
； （1）
其中，为守恒变量，/>分别为直角坐标系x,y,z方向的无粘通量；
在固定网格下，进行坐标变换：
； （2）
得到计算坐标下的欧拉方程为：
； （3）
其中，
； （4）
其中雅克比和网格度量公式采用对称守恒计算形式，具体如下：
；
；
；
；
其中，、/>、/>、/>、/>、/>、/>、/>、/>分别为一阶偏导数；
需要注意的是，在计算空间进行等距离散，若以f表示为某一特定物理变量，以方向为例，上述公式（5）中的一阶偏导数离散方式，采用下述公式（6）进行离散：
；
其中，下标i为网格节点指标、i+1/2为离散后半点上指标、h为网格间距；
半点处i+1/2的物理量f的值由四阶中心差分插值得到：
； （7）
将计算得到的网格度量和雅克比代入上述公式（4）中，得到计算空间中各网格点上的守恒变量和通量/>。由于有限差分法可分维求解，之后为了方便叙述，会省略相同的下标，如在/>方向进行求解，/>和/>方向为某一常数，则/>简记为/>。
在具体实施时，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，在执行步骤S102时，如图2所示，具体可以包括以下步骤：
首先考虑方向，采用四阶中心型非线性插值（即四点模板下半点处的中心型非线性插值）的方法将节点处物理量/>插值到半点处；半点处的左值记为/>，半点处的右值记为/>，具体如下：
；
；
； （8）
其中非线性权值求解如下：
；
；/>
； （9）
其中全局光滑度量为：
； （10）
当地光滑度量为：
；
。 （11）
进一步地，在具体实施时，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，在执行步骤S102时，如图2所示，具体可以还包括以下步骤：
将前述步骤得到的半点处左值和右值/>，计算半点处的低阶通量，可以选取的黎曼通量由很多种，这里以Lax-Friedrich分裂为例，具体如下：
；
（12）
其中，是通量雅克比矩阵/>的全局谱半径。
在具体实施时，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，在执行步骤S103时，如图2所示，具体可以包括以下步骤：
根据当前时刻物理量，计算半点处的修正通量/>。具体可以根据采用下述公式计算半点处的修正通量：/>
；（13）
其中，为半点处的修正通量，/>为/>方向网格间距，/>为/>在半点i+1/2处/>方向的二阶偏导数。
在具体实施时，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，在执行步骤S104时，如图2所示，具体可以包括以下步骤：
将步骤S102计算出的低阶通量和步骤S103计算出的修正通量/>求和，得到半点处的高阶通量/>，具体如下：
；
（14）
采用下述公式计算方向各节点处一阶导数/>：
。（15）
另外，针对方向和/>方向的一阶导数/>和/>，重复/>方向步骤S102至步骤S104的操作，具体公式如下：
；
；
（16）
将上述公式（15）和（16）代入上述公式（3），得到半离散形式：
。 （17）
至此，欧拉方程的空间离散完成。
需要说明的是，为了提高四点模板上经典WENO-JS格式的性能，并进一步降低网格度量计算所带来的数值误差，本发明在四点模板下设计了四阶中心型非线性插值方法，见公式（8）-（11）；另外，针对网格度量、雅克比及通量中的一阶导数构造了四阶中心差分离散方法，见公式（6）、（15）、（16）。
在具体实施时，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，在执行步骤S105时，具体可以包括以下步骤：
由于此时方程（3）变为常微分方程组，采用三阶龙格-库塔法进行时间推进，如下：
； （18）
其中，上标n为第n时刻步的值，上标n+1为第n+1时刻步的值，RHS为右端项格式空间离散的值。至此，欧拉方程的时空离散完成。
在具体实施时，在本发明实施例提供的上述可压缩流动高分辨率数值模拟方法中，步骤S106将时间推进至设定时刻结束计算，得到设定时刻的流场数据，如图2所示，具体可以包括：对于空间离散和时空离散后的欧拉方程，采用R-K时间推进法得到t n-1时刻的流场数据；判断t n-1-t N 是否小于，若是，则表示t n-1时刻的流场数据为最终时刻t N 的流场数据；若否，则继续进行时间推进，计算流场数据。
下面以第一测试算例和第二测试算例为例，对本发明实施例提供的上述可压缩流动高分辨率数值模拟方法进行测试。其中第一测试算例是等熵涡传播流动，采用规则网格和随机网格进行计算，用来测试本发明离散方法的分辨能力；第二测试算例是Ms=10双马赫反射问题，采用随机网格进行计算，用来测试计算方法对强激波的捕捉能力。
从计算结果可以看出：本发明实施例提供的上述可压缩流动高分辨率数值模拟方法能得到较好的计算结果。具体表现在：如图3a和3b所示，从等熵涡传播流动看，在规则网格下本发明方法分辨率有明显提高；另外，如图4a和4b所示，在有畸变的随机网格下经典WENO-JS格式容易发散，而本发明方法有更好的涡保持特性。如图5所示，从双马赫反射算例可以看到本发明方法能够稳定的捕捉包含强激波的流动。
相应地，本发明还公开了一种计算机可读存储介质，用于存储计算机程序；计算机程序被处理器执行时实现前述公开的可压缩流动高分辨率数值模拟方法。关于上述方法更加具体的过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
进一步地，本发明实施例还公开了一种可压缩流动高分辨率数值模拟设备，包括处理器和存储器；其中，处理器执行存储器中存储的计算机程序时实现前述实施例公开的可压缩流动高分辨率数值模拟方法。关于上述方法更加具体的过程可以参考前述实施例中公开的相应内容，在此不再进行赘述。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其它实施例的不同之处，各个实施例之间相同或相似部分互相参见即可。对于实施例公开的存储介质、设备而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
专业人员还可以进一步意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本申请的范围。
结合本文中所公开的实施例描述的方法或算法的步骤可以直接用硬件、处理器执行的软件模块，或者二者的结合来实施。软件模块可以置于随机存储器（RAM）、内存、只读存储器（ROM）、电可编程ROM、电可擦除可编程ROM、寄存器、硬盘、可移动磁盘、CD-ROM、或技术领域内所公知的任意其它形式的存储介质中。
最后，还需要说明的是，在本文中，诸如第一和第二等之类的关系术语仅仅用来将一个实体或者操作与另一个实体或操作区分开来，而不一定要求或者暗示这些实体或操作之间存在任何这种实际的关系或者顺序。而且，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、物品或者设备中还存在另外的相同要素。
以上对本发明所提供的可压缩流动高分辨率数值模拟方法、介质及设备进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (5)

1.一种可压缩流动高分辨率数值模拟方法，其特征在于，包括：
读取初始流场数据，得到当前时刻计算空间中各网格点的变量和通量，包括：读取t 0时刻流场数据，采用四阶中心差分离散方法计算所述t 0时刻流场数据对应的网格度量和雅克比，并将t n 时刻物理空间中各网格点的变量/>和通量/>进行转换，得到t n 时刻计算空间/>中各网格点的变量/>和通量/>；
在无粘假设条件下，可压缩流动的控制方程为欧拉方程，在笛卡尔坐标系下，无量纲形式三维守恒型的欧拉方程为：
；
其中，为守恒变量，/>分别为直角坐标系下x, y, z方向的无粘通量；
在固定网格下，进行坐标变换：
；
得到计算坐标下的欧拉方程为：
；
其中，
；
其中雅克比和网格度量公式采用对称守恒计算形式：
；
；
；
；
其中，、/>、/>、/>、/>、/>、/>、/>、/>分别为一阶偏导数；
在所述计算空间进行等距离散，若以f表示为某一特定物理变量，以方向为例，上述公式中的一阶偏导数离散方式，采用下述公式离散：
；
其中，下标i为网格节点指标、i+1/2为离散后半点上指标、h为网格间距；
半点处i+1/2的物理量f的值由四阶中心差分插值得到：
；
根据计算得到的网格度量和雅克比，得到计算空间中各网格点上的守恒变量和通量/>；
根据得到的所述各网格点的变量，采用非线性插值计算半点处的左值和右值，并根据所述左值和右值，计算半点处的低阶通量，包括：
考虑方向，采用四阶中心型非线性插值的方法将节点处物理量/>插值到半点处；半点处的左值记为/>，半点处的右值记为/>：
；
；
；
其中非线性权值求解如下：
；
；
；
其中全局光滑度量为：
；
当地光滑度量为：
；
；
根据得到的所述各网格点的通量，采用下述公式计算半点处的修正通量：
；
其中，为半点处的修正通量，/>为/>方向网格间距，/>为/>在半点i+1/2处/>方向的二阶偏导数；
根据计算出的所述低阶通量和所述修正通量，获取所述计算空间中不同方向各节点处的一阶导数，以完成空间离散；其中，获取方向各节点处的一阶导数，包括：
将计算出的所述低阶通量和所述修正通量/>求和，得到半点处的高阶通量：
；
采用下述公式计算方向各节点处一阶导数/>：
；
采用龙格-库塔法将当前时刻的变量推进到下一时刻的变量，以完成时空离散；
将时间推进至设定时刻结束计算，得到所述设定时刻的流场数据。
2.根据权利要求1所述的可压缩流动高分辨率数值模拟方法，其特征在于，采用下述公式计算方向和/>方向每个节点处一阶导数/>和/>：
；
；
欧拉方程的半离散形式为：
。
3.根据权利要求2所述的可压缩流动高分辨率数值模拟方法，其特征在于，采用龙格-库塔法将当前时刻的变量推进到下一时刻的变量，包括：
采用三阶龙格-库塔法进行时间推进：
；
其中，上标n为第n时刻步的值，上标n+1为第n+1时刻步的值，RHS为右端项格式空间离散的值。
4.一种计算机可读存储介质，其特征在于，用于存储计算机程序，其中，所述计算机程序被处理器执行时实现如权利要求1至3任一项所述的可压缩流动高分辨率数值模拟方法。
5.一种可压缩流动高分辨率数值模拟设备，其特征在于，包括处理器和存储器，其中，所述处理器执行所述存储器中存储的计算机程序时实现如权利要求1至3任一项所述的可压缩流动高分辨率数值模拟方法。
CN202310654382.9A 2023-06-05 2023-06-05 一种可压缩流动高分辨率数值模拟方法、介质及设备 Active CN116384288B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310654382.9A CN116384288B (zh) 2023-06-05 2023-06-05 一种可压缩流动高分辨率数值模拟方法、介质及设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310654382.9A CN116384288B (zh) 2023-06-05 2023-06-05 一种可压缩流动高分辨率数值模拟方法、介质及设备

Publications (2)

Publication Number Publication Date
CN116384288A CN116384288A (zh) 2023-07-04
CN116384288B true CN116384288B (zh) 2023-08-25

Family

ID=86971634

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310654382.9A Active CN116384288B (zh) 2023-06-05 2023-06-05 一种可压缩流动高分辨率数值模拟方法、介质及设备

Country Status (1)

Country Link
CN (1) CN116384288B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103914602A (zh) * 2012-12-30 2014-07-09 西安远景动力模拟技术有限公司 可压缩旋流场的数值模拟方法
CN112214869A (zh) * 2020-09-03 2021-01-12 空气动力学国家重点实验室 一种求解欧拉方程的改进型高阶非线性空间离散方法
CN113486453A (zh) * 2021-09-07 2021-10-08 中国空气动力研究与发展中心低速空气动力研究所 一种可压缩流体的数值模拟方法
CN114638173A (zh) * 2022-01-25 2022-06-17 中国空气动力研究与发展中心计算空气动力研究所 一种高阶非线性激波捕捉空间离散方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102203782B (zh) * 2010-09-09 2013-03-13 天津空中代码工程应用软件开发有限公司 求解拉格朗日形式的欧拉方程的数值方法
US20160306907A1 (en) * 2012-09-18 2016-10-20 Ming Lu Numerical method for solving the two-dimensional riemann problem to simulate inviscid subsonic flows

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103914602A (zh) * 2012-12-30 2014-07-09 西安远景动力模拟技术有限公司 可压缩旋流场的数值模拟方法
CN112214869A (zh) * 2020-09-03 2021-01-12 空气动力学国家重点实验室 一种求解欧拉方程的改进型高阶非线性空间离散方法
CN113486453A (zh) * 2021-09-07 2021-10-08 中国空气动力研究与发展中心低速空气动力研究所 一种可压缩流体的数值模拟方法
CN114638173A (zh) * 2022-01-25 2022-06-17 中国空气动力研究与发展中心计算空气动力研究所 一种高阶非线性激波捕捉空间离散方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
可压缩多相流问题的高精度高分辨率数值模拟方法研究;田保林 等;《第九届全国流体力学学术会议论文摘要集》;278 *

Also Published As

Publication number Publication date
CN116384288A (zh) 2023-07-04

Similar Documents

Publication Publication Date Title
Zhuang et al. Fracture modeling using meshless methods and level sets in 3D: framework and modeling
Chakraverty et al. Advanced numerical and semi-analytical methods for differential equations
US20200043231A1 (en) Method for constructing three-dimensional solid model with geometric error and computer-readable storage medium
CN112214869B (zh) 一种求解欧拉方程的改进型高阶非线性空间离散方法
Zhou et al. Shape variable radial basis function and its application in dual reciprocity boundary face method
Li The meshless Galerkin boundary node method for Stokes problems in three dimensions
CN109858158B (zh) 一种计算流体力学模拟的参数配置方法及系统
CN107577857B (zh) 一种基于热辐射边界条件的三维有限元模拟方法
Michoski et al. Foundations of the blended isogeometric discontinuous Galerkin (BIDG) method
Qian et al. The characteristic variational multiscale method for convection-dominated convection–diffusion–reaction problems
Cochez‐Dhondt et al. Equilibrated error estimators for discontinuous Galerkin methods
CN113792445B (zh) 一种基于积分方程法的三维大地电磁数值模拟方法
CN115932960A (zh) 全矩张量震源机制反演方法、系统、设备及存储介质
CN116384288B (zh) 一种可压缩流动高分辨率数值模拟方法、介质及设备
CN114139109A (zh) 一种目标跟踪方法、系统、设备、介质及数据处理终端
CN107886573B (zh) 一种复杂地质条件下边坡三维有限元网格生成方法
Juhl A note on the convergence of the direct collocation boundary element method
Benkhaldoun et al. A family of finite volume Eulerian–Lagrangian methods for two-dimensional conservation laws
CN115439330A (zh) 高精度局部三维坐标转换方法及其装置
CN111158059A (zh) 基于三次b样条函数的重力反演方法
Fernàndez Méndez Mesh-free methods and finite elements: friend or foe?
JP5039003B2 (ja) 高速磁場解析方法、高速磁場解析プログラム、および記録媒体
Hu et al. An efficient high‐precision recursive dynamic algorithm for closed‐loop multibody systems
CN113947628A (zh) 一种基于不同比尺河道地形图计算冲淤量的方法
OKHOVATI et al. Numerical coupling of two scalar conservation laws by a RKDG method

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