CN114638173B - 一种高阶非线性激波捕捉空间离散方法 - Google Patents

一种高阶非线性激波捕捉空间离散方法 Download PDF

Info

Publication number
CN114638173B
CN114638173B CN202210086558.0A CN202210086558A CN114638173B CN 114638173 B CN114638173 B CN 114638173B CN 202210086558 A CN202210086558 A CN 202210086558A CN 114638173 B CN114638173 B CN 114638173B
Authority
CN
China
Prior art keywords
flux
format
grid
discontinuity
grid point
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
CN202210086558.0A
Other languages
English (en)
Other versions
CN114638173A (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 CN202210086558.0A priority Critical patent/CN114638173B/zh
Publication of CN114638173A publication Critical patent/CN114638173A/zh
Application granted granted Critical
Publication of CN114638173B publication Critical patent/CN114638173B/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
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4007—Scaling of whole images or parts thereof, e.g. expanding or contracting based on interpolation, e.g. bilinear interpolation
    • 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)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computing Systems (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Complex Calculations (AREA)

Abstract

本发明适用于流体力学计算领域，提供了一种高阶非线性激波捕捉空间离散方法，包括如下步骤：读取流场数据，得到各网格点上的物理通量和特征通量；根据所述特征通量判断各网格点的间断性，所述间断性包括光滑和间断；当所述网格点的间断性为光滑时，通过k次样条插值格式得到半网格点上通量；当所述网格点的间断性为间断时，通过单调保持格式得到半网格点上通量；根据所述半网格点上通量，并结合各网格点上的物理通量，得到整体的空间离散格式；将时间推进至预设结束时间，得到预设结束时间的流场数据。本发明提供的一种高阶非线性激波捕捉空间离散方法具有低耗散、高精度、高分辨率和间断的精确捕捉等优势。

Description

一种高阶非线性激波捕捉空间离散方法
技术领域
本发明涉及流体力学计算领域，尤其是涉及一种高阶非线性激波捕捉空间离散方法。
背景技术
随着航空技术飞速发展，先进飞行器的气动噪声问题越来越受到广泛重视。航空界把飞行器的声舒适性和对环境噪声影响作为研究的出发点，尽量降低飞行器的气动噪声，提高商用飞机的竞争能力。求解包含激波的气动声学问题是一个重要的研究方向，数值模拟这类问题对计算格式提出了更高的要求：求解声场要求数值格式具有高阶精度、高分辨率和低耗散的特性，而模拟激波运动还要求数值格式必须能够捕捉激波。目前，还不存在比较理想的模拟这类激波噪声问题的数值格式，但激波噪声却是航天飞行器噪声的一个重要来源，具有重要的工程应用背景和理论研究价值。因此非常有必要研究适合于模拟激波噪声问题的数值格式，从而为我们开展激波噪声问题的研究提供帮助，并用于激波噪声产生、传播过程的数值模拟，促进对物理本质的认识，从而对激波噪声的降低与控制提供新的手段。
其中，控制方程中对流项的数值离散会直接影响无粘流区域的激波的分辨，并间接影响激波噪声的预测。由于中心紧致差分格式没有数值耗散并且色散误差小，所以对于计算多尺度问题有着天然的优势。但是直接把高阶中心紧致格式应用到非线性双曲守恒律也不现实，因为差分离散格式不能分辨的高频波将会累积误差，最后导致数值不稳定。
传统的激波捕捉格式WENO格式虽然能够达到较好的分辨率，但其计算模板上节点数据较多，为边界处理带来了一定程度的困难，且增加了计算工作量，耗散较高，在光滑解极值点并没有达到目标精度，且不能够很好的计算剪切层等精细结构，并且WENO格式对噪声的识别能力不够准确。
直接采用线性中心紧致差分格式应用到非线性双曲守恒律，则差分离散格式不能分辨的高频波将会累积误差，最后导致数值不稳定。为了克服中心紧致格式的数值不稳定缺陷，加入人工耗散来消去非物理的波动，这样的经验性很强，不好控制。
发明内容
本发明的目的是解决当流场存在激波间断时出现非物理振荡而导致的数值不稳定性的技术问题，提供了一种高阶非线性激波捕捉空间离散方法，具有低耗散、高精度、高分辨率和间断的精确捕捉等优势。
本发明提供了一种高阶非线性激波捕捉空间离散方法，包括如下步骤：
步骤S10：读取流场数据，得到各网格点上的物理通量和特征通量；
步骤S20：根据所述特征通量判断各网格点的间断性，所述间断性包括光滑和间断；
步骤S30：当所述网格点的间断性为光滑时，通过k次样条插值格式得到半网格点上通量，其中，k为样条插值格式的最高次数项，k为正整数且k≥3；当所述网格点的间断性为间断时，通过单调保持格式得到半网格点上通量；
步骤S40：根据所述半网格点上通量，并结合各网格点上的物理通量，得到整体的空间离散格式；
步骤S50：将时间推进至预设结束时间，得到预设结束时间的流场数据。
进一步的，步骤S20中，判断各网格点的间断性的方法为：同时满足
Figure RE-GDA0003640969000000021
且/>
Figure RE-GDA0003640969000000022
其中，j为网格点的编号，j为正整数，/>
Figure RE-GDA0003640969000000023
为半网格点上通量，uj为第j个网格点的特征通量，uMP为初始单调保持通量，/>
Figure RE-GDA0003640969000000031
为半网格点的光滑因子。
进一步的，所述初始单调保持通量uMP的计算公式为uMP＝uj+midmod[uj+1-uj，4(uj-uj-1)]，其中，uj+1为第j+1个网格点的特征通量，uj-1为第j-1个网格点的特征通量。
进一步的，所述半网格点的光滑因子
Figure RE-GDA0003640969000000032
的计算公式为
Figure RE-GDA0003640969000000033
其中，Mj+1为第j+1个网格点上k次样条插值函数的k-1次导数，Mj为第j个网格点上k次样条插值函数的 k-1次导数，k次样条插值函数为S(j)，a≤j≤b，b-a小于预设区间，h为相邻网格点之间的距离。
进一步的，步骤S30中，当所述网格点的间断性为光滑时，通过三次样条插值格式得到半网格点上通量。
进一步的，半网格点上通量
Figure RE-GDA0003640969000000034
的计算公式为/>
Figure RE-GDA0003640969000000035
其中，/>
Figure RE-GDA0003640969000000036
为三次样条插值。
进一步的，步骤S30中，当所述网格点的间断性为间断时，半网格点上通量
Figure RE-GDA0003640969000000037
的计算公式为/>
Figure RE-GDA0003640969000000038
其中，/>
Figure RE-GDA0003640969000000039
为第 j个半网格点特征通量的最小值，/>
Figure RE-GDA00036409690000000310
为第j个半网格点的特征通量的最大值。
进一步的，步骤S10中所述特征通量的计算方法为：根据控制方程计算各网格点上的初始通量，通过对所述初始通量分裂得到正负通量，对所述正负通量进行特征投影，得到特征通量。
进一步的，步骤S40中所述空间导数通过线性紧致差分格式得到，所述紧致差分格式为
Figure RE-GDA0003640969000000041
其中，a为第一系数，b为第二系数，α为第三系数，f′j-1为第j-1个网格点上物理通量的导数的差分近似式，f′j为第j个网格点上物理通量的导数的差分近似式，f′j+1为第j+1个网格点上物理通量的导数的差分近似式，/>
Figure RE-GDA0003640969000000042
为第j个网格点的右半网格点上通量/>
Figure RE-GDA0003640969000000043
的反投影，/>
Figure RE-GDA0003640969000000044
为第j个网格点的左半网格点上通量
Figure RE-GDA0003640969000000045
的反投影，fj+1为第j+1 个网格点上物理通量，fj-1为第j-1个网格点上物理通量；
当j为网格内点的编号时，第三系数
Figure RE-GDA0003640969000000046
此时紧致差分格式为四阶显式格式；当j为网格边界点的编号时，第三系数α＝0，此时紧致差分格式为四阶紧致格式；
所述第一系数a的计算公式为
Figure RE-GDA0003640969000000047
所述第二系数b的计算公式为
Figure RE-GDA0003640969000000048
进一步的，步骤S50中采用三阶龙格库塔方法对时间项进行离散，将时间推进至预设结束时间。
综上所述，本发明至少能够实现如下技术效果：
1.对于光滑流场，多次样条插值耦合中心差分格式可以很精确地计算流场。本发明在当网格点的间断性为光滑时，通过k次样条插值格式得到半网格点上通量，实现了流场单元界面处通量的的快速求解，同时提高了光滑流场的计算精度；
2.当流场存在间断时，非物理振荡可能导致数值不稳定性，为了能够计算包含激波等的复杂多尺度问题，将多次样条插值后的半网格点通量采用单调保持格式来修正，获得对激波等间断问题的精确捕捉，实现对超声速问题高精度、高分辨率的计算；
3.本发明还提出了一种新的网格点间断性的判断方法，能够直接由空间格式自身中计算得到，具有节省计算步骤，计算效率高的特点；
4.本发明通过单调保持非线性样条插值的办法来求解半网格点处的函数值，从而得到了高阶非线性紧致格式，所述高阶非线性紧致格式保持了原有中心紧致差分格式的高阶精度、高分辨率和低耗散特性，同时由于其单调保持非线性插值的机制，使得高阶非线性紧致格式能够捕捉强激波，一系列数值算例表明所述样条插值非线性紧致格式是高阶精度、高分辨率和低耗散数值计算方法，具有广泛的科研以及应用前景。
附图说明
为了更清楚地说明本发明实施例的技术方案，下面将对本发明实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面所描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明中高阶非线性激波捕捉空间离散方法的示意图；
图2是Exact、MPSI4、CNCS4和WENO5的分辨率曲线；
图3是本发明中MPSI4格式计算Lax问题；
图4是本发明中MPSI4格式计算激波双马赫反射；
图5是本发明中MPSI4格式计算前台阶问题；
图6是本发明中MPSI4格式计算气动噪声传播问题；
图7是本发明中MPSI4格式计算气动噪声传播的噪声脉动分布。
具体实施方式
以下的说明提供了许多不同的实施例、或是例子，用来实施本发明的不同特征。以下特定例子所描述的元件和排列方式，仅用来精简的表达本发明，其仅作为例子，而并非用以限制本发明。
下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。同时，在本发明的描述中，术语“第一”、“第二”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
实施例一：
如图1所示，本发明实施例一提供了一种高阶非线性激波捕捉空间离散方法，包括如下步骤：
步骤S10：读取流场数据，得到各网格点上的物理通量和特征通量；
步骤S20：根据所述特征通量判断各网格点的间断性，所述间断性包括光滑和间断；
步骤S30：当所述网格点的间断性为光滑时，通过k次样条插值格式得到半网格点上通量，其中，k为样条插值格式的最高次数项，k为正整数且k≥3；当所述网格点的间断性为间断时，通过单调保持格式得到半网格点上通量；
步骤S40：根据所述半网格点上通量，并结合各网格点上的物理通量，得到整体的空间离散格式；
步骤S50：将时间推进至预设结束时间，得到预设结束时间的流场数据。
对于光滑流场，多次样条插值耦合中心差分格式可以很精确地计算流场。所述k次样条插值格式的最高次数项k仅需要满足k为正整数且k≥3即可，具体的最高次数项k的数值越大对间断的捕捉精度越高，得到的分辨率也会越高，如五次样条插值格式的计算精度>四次样条插值格式的计算精度>三次样条差值格式的计算精度。
本实施例在当网格点的间断性为光滑时，通过k次样条插值格式得到半网格点上通量，实现了流场单元界面处通量的的快速求解，同时提高了光滑流场的计算精度。
当流场存在间断时，非物理振荡可能导致数值不稳定性，为了能够计算包含激波等的复杂多尺度问题，将多次样条插值后的半网格点通量采用单调保持格式来修正，获得对激波等间断问题的精确捕捉，实现对超声速问题高精度、高分辨率的计算。
进一步的，步骤S20中，判断各网格点的间断性的方法为：同时满足
Figure RE-GDA0003640969000000071
且/>
Figure RE-GDA0003640969000000072
其中，j为网格点的编号，j为正整数，/>
Figure RE-GDA0003640969000000073
为半网格点上通量，uj为第j个网格点的特征通量，uMP为初始单调保持通量，/>
Figure RE-GDA0003640969000000074
为半网格点的光滑因子。
所述网格点间断性的判断方法能够直接由空间格式自身中计算得到，具有节省计算步骤，计算效率高的特点。
进一步的，所述初始单调保持通量uMP的计算公式为 uMP＝uj+mid mod[uj+1-uj,4(uj-uj-1)]，其中，uj+1为第j+1个网格点的特征通量，uj-1为第j-1个网格点的特征通量。
进一步的，所述半网格点的光滑因子
Figure RE-GDA0003640969000000075
的计算公式为
Figure RE-GDA0003640969000000076
其中，Mj+1为第j+1个网格点上k次样条插值函数的k-1次导数，Mj为第j个网格点上k次样条插值函数的 k-1次导数，k次样条插值函数为S(j)，a≤j≤b，b-a小于预设区间，h为相邻网格点之间的距离。
进一步的，步骤S30中，当所述网格点的间断性为光滑时，通过三次样条插值格式得到半网格点上通量。
对于k次样条插值格式的最高次数项k，虽然最高次数项k的数值越大对间断的捕捉精度越高，得到的分辨率也会越高，但最高次数项k的数值越大计算过程也将更加复杂，运算效率也会降低。因而，k次样条插值格式的最高次数项k优选为3，此时，既能够实现高精度、高分辨率，也可确保计算过程的简单化和运算的高效率。
进一步的，半网格点上通量
Figure RE-GDA0003640969000000081
的计算公式为/>
Figure RE-GDA0003640969000000082
其中，/>
Figure RE-GDA0003640969000000083
为三次样条插值。
以半网格点j+1/2为例，本实施例中以三次样条插值为例，构造半网格点上通量的高阶混合计算方法计算步骤如下：
将三次样条插值定义为函数S(x)∈C2[a,b]，在每个小区间上是三次多项式，其中a＝x0＜x1＜…＜xn-1＜xn＝b是给定的网格节点。若函数S(x)在节点上满足以下公式，则称S(x)是三次样条函数：
S0(x0)＝f(x0),Sn-1(xn)＝f(xn)
Sj-1(xj)＝f(xj)＝Sj(xj),j＝1,…,n-1
S′j-1(xj)＝S′j(xj)
S″j-1(xj)＝S″j(xj)
上述方程共有4n-2个，为求解所有Sj(xj)的系数，还需两个方程，通常根据物理边界条件来得到，具体如下：
假设S″(xj)＝Mj，则得到：
Figure RE-GDA0003640969000000084
其中hj＝xj+1-xj。对上式积分两次，则可以得到如下两个函数：
Figure RE-GDA0003640969000000091
利用一阶导数的接触面条件S′j-1(xj)＝S′j(xj)可得：
Figure RE-GDA0003640969000000092
由上式可得三对角方程组：μjMj-1+2Mj+λjMj+1＝dj,j＝1,…,n-1
其中
Figure RE-GDA0003640969000000093
上述方程组联合边界条件可以求解出S″(xj)＝Mj，从而可以得到：
Figure RE-GDA0003640969000000094
当计算网格是均匀网格时，可以简化为：
Figure RE-GDA0003640969000000095
由基础数值分析可知上述逼近为4阶精度误差，即
Figure RE-GDA0003640969000000096
所述三次样条插值
Figure RE-GDA0003640969000000097
即采用Hermit插值多项式显式表达的紧致格式。在均匀网格上三次样条插值耦合四阶中心紧致差分格式可得到四阶差分格式。对于光滑流场，三次样条插值耦合中心差分格式/>
Figure RE-GDA0003640969000000098
能够实现流场的精确计算。
采用三次样条插值格式得到三次样条插值耦合中心差分格式
Figure RE-GDA0003640969000000101
之后，当流场处于光滑区域时，半网格点上通量即为三次样条插值耦合中心差分格式，即
Figure RE-GDA0003640969000000102
进一步的，步骤S30中，当所述网格点的间断性为间断时，半网格点上通量
Figure RE-GDA0003640969000000103
的计算公式为/>
Figure RE-GDA0003640969000000104
其中，/>
Figure RE-GDA0003640969000000105
为第j个半网格点特征通量的最小值，/>
Figure RE-GDA0003640969000000106
为第j个半网格点的特征通量的最大值。
当流场处于间断区域时，首先取线性插值为
Figure RE-GDA0003640969000000107
然后根据单调保持(Monotonicity–Preserving)限制器中的单调保持格式来进行非线性修正，具体步骤如下：
dj＝uj-1-2uj+uj+1；
Figure RE-GDA0003640969000000108
Figure RE-GDA0003640969000000109
Figure RE-GDA00036409690000001010
Figure RE-GDA00036409690000001011
Figure RE-GDA00036409690000001012
Figure RE-GDA0003640969000000111
Figure RE-GDA0003640969000000112
其中，为计算方便，令最小模函数
Figure RE-GDA0003640969000000113
中的
Figure RE-GDA0003640969000000114
则最小模函数的计算过程为：
Figure RE-GDA0003640969000000115
min mod(x1,…,xk)＝s·min(|x1|,|x2|,…,xk∣)
Figure RE-GDA0003640969000000116
Figure RE-GDA0003640969000000117
median(x,y,z)＝x+min mod(y-x,z-x)
根据上述步骤，可以得到半网格点上通量
Figure RE-GDA0003640969000000118
由于/>
Figure RE-GDA0003640969000000119
和/>
Figure RE-GDA00036409690000001110
是关于/>
Figure RE-GDA00036409690000001111
对称的，同样的计算步骤可以对称地计算出/>
Figure RE-GDA00036409690000001112
计算得到
Figure RE-GDA00036409690000001113
和/>
Figure RE-GDA00036409690000001114
之后再求和，然后从特征空间返回到物理空间，即可得到半网格点的特征通量/>
Figure RE-GDA00036409690000001115
进一步的，步骤S10中所述特征通量的计算方法为：根据控制方程计算各网格点上的初始通量，通过对所述初始通量分裂得到正负通量，对所述正负通量进行特征投影，得到特征通量。
控制方程为在不考虑质量力和源项的情况下，无量纲积分形式的三维Navier-Stokes控制方程：
Figure RE-GDA0003640969000000121
其中，U为守恒变量，Fe，Ge，Ge为无粘通量，
Figure RE-GDA00036409690000001212
为粘性通量。无粘通量的具体表达式如下：
Figure RE-GDA0003640969000000122
粘性通量的表达式为：
Figure RE-GDA0003640969000000123
Figure RE-GDA0003640969000000124
其中：
Figure RE-GDA0003640969000000125
Figure RE-GDA0003640969000000126
在牛顿流体中，粘性应力张量各个分量分别为：
Figure RE-GDA0003640969000000127
/>
Figure RE-GDA0003640969000000128
Figure RE-GDA0003640969000000129
其中，μ为动力粘性系数，对于空气，μ可由Sutherland公式计算得出。
Figure RE-GDA00036409690000001210
为热传导系数，其中cp为定压比热，Pr为普朗特数。
特征矩阵
Figure RE-GDA00036409690000001211
的具体形式如下：
Figure RE-GDA0003640969000000131
其中，
Figure RE-GDA0003640969000000132
Λk为特征矩阵的特征值矩阵，特征值代表了各类波的特征速度，其具体形式如下：
Figure RE-GDA0003640969000000133
Lk和Rk分别为左、右特征矩阵，具体形式为：
Figure RE-GDA0003640969000000134
/>
Figure RE-GDA0003640969000000135
其中，
Figure RE-GDA0003640969000000141
a为当地声速。
在数值计算中，根据特征传播方向，通常将无粘通量分为正、负两个部分，上标“+”表示这部分的无粘矢通量的特征速度指向计算坐标轴的正方向，上标“-”则表示该部分的特征速度指向计算坐标轴的负方向传播。
根据Local Lax-Friedrichs分裂方法，得到正负通量
Figure RE-GDA0003640969000000142
计算公式如下：
Figure RE-GDA0003640969000000143
其中，
Figure RE-GDA0003640969000000144
采用左特征矩阵
Figure RE-GDA0003640969000000145
对前面得到的/>
Figure RE-GDA0003640969000000146
进行投影，得到第j个半网格点的特征通量/>
Figure RE-GDA0003640969000000147
进一步的，步骤S40中所述空间导数通过线性紧致差分格式得到，所述紧致差分格式为
Figure RE-GDA0003640969000000148
其中，a为第一系数，b为第二系数，α为第三系数，f′j-1为第j-1个网格点上物理通量的导数的差分近似式，f′j为第j个网格点上物理通量的导数的差分近似式，f′j+1为第j+1 个网格点上物理通量的导数的差分近似式，/>
Figure RE-GDA0003640969000000149
为第j个网格点的右半网格点上通量/>
Figure RE-GDA00036409690000001410
的反投影，/>
Figure RE-GDA00036409690000001411
为第j个网格点的左半网格点上通量/>
Figure RE-GDA00036409690000001412
的反投影，fj+1为第j+1个网格点上物理通量，fj-1为第j-1个网格点上物理通量；/>
当j为网格内点的编号时，第三系数
Figure RE-GDA00036409690000001413
此时紧致差分格式为四阶显式格式；当j为网格边界点的编号时，第三系数α＝0，此时紧致差分格式为四阶紧致格式；
所述第一系数a的计算公式为
Figure RE-GDA0003640969000000151
所述第二系数b的计算公式为/>
Figure RE-GDA0003640969000000152
中心紧致格式中的网格中心处的值采用三次样条插值可以得到四阶紧致格式，进而得到整个的空间离散格式，将本实施例中的空间离散格式记为四阶 MPSI(Monotonicitypreserving Spline Interpolation)格式，简称MPSI4。
当k次样条插值格式的最高次数项k大于等于3时，则得到k+1阶紧致格式，此时空间离散格式记为k阶MPSI(Monotonicity preserving Spline Interpolation)格式，简称MPSIk。
进一步的，步骤S50中采用三阶龙格库塔方法对时间项进行离散，将时间推进至预设结束时间。
进行空间离散之后，得到关于时间导数的半离散形式为
Figure RE-GDA0003640969000000153
采用显式Runge-Kutta时间推进方法来离散时间导数项。
优化的三阶TVD Runge-Kutta时间离散格式为： U(1)＝Un+ΔtL(Un)
Figure RE-GDA0003640969000000154
Figure RE-GDA0003640969000000155
其中，n为第n时刻步的值。至此，完成了对欧拉或者纳维斯托克斯方程的高阶时空离散，从而将时间从开始时间tini一直推进至预设结束时间tend。
实施例二：
本发明实施例二采用标准算例来验证四阶MPSI格式的精度，令ut+ux＝0，其中-1≤x≤1。设定初始条件为u(x,t＝0)＝sin(πx)，精确解为u(x,t)＝sin(π(x-t))，给定周期性的边界条件，计算t＝1。
求出数值解和精确解在不同网格密度下的L1误差，L1数值精度阶数，L∞误差，L∞数值精度阶数。时间离散采用三阶TVD Runge-Kutta格式，数值精度的验证如表1所示，其中，N为网格数，从表1中可以看到MPSI4格式确实达到了理论的收敛精度。
表1四阶MPSI格式和三阶TVD Runge-Kutta格式的数值精度验证
N L1误差 L1数值精度阶数 L∞误差 L∞数值精度阶数
20 0.109E-03 -- 0.172E-03 --
40 0.678E-05 4.01 0.107E-04 4.01
60 0.134E-05 4.00 0.210E-05 4.00
80 0.423E-06 4.00 0.665E-06 4.00
100 0.173E-06 4.00 0.272E-06 4.00
同时，采用非线性格式的频谱分析方法来研究MPSI4格式的分辨率。如图 2所示为Exact、MPSI4、CNCS4和WENO5的分辨率曲线，表明MPSI4格式的分辨率远优于传统的紧致差分格式和激波捕捉WENO格式。
(1)Lax激波管问题
如图3所示为采用MPSI4格式计算Lax问题(N＝100)，该问题的初始条件为：
Figure RE-GDA0003640969000000161
计算区域取为-5≤x≤5，计算网格采用100个点，初始间断位于x＝0处，最终计算时刻取为t＝1.3。
采用MPSI4格式计算一维激波管Lax问题，计算结果在激波、膨胀波及接触间断附近都能达到基本无波动。
(2)激波双马赫反射问题
如图4所示为采用MPSI4格式计算激波双马赫反射(N＝960×240)，双马赫反射问题描述的是Mach数为10的强运动斜激波以与x轴方向呈60度角的方向入射，入射点在(1/6,0)，计算区域取[0,4]×[0,1]。激波波前静止空气的密度为1.4，压力为1，波后按激波关系给定初始条件。下边界在[1/6,4]的区域给定壁面边界条件，其它边界区域按照激波运动所在的位置，分别给定波前或波后的值。初始物理参数为：
Figure RE-GDA0003640969000000171
本发明实施例中采用960×240的均匀网格计算到无量纲时间t＝0.2的时刻。计算结果在激波附近都能达到基本无波动，对剪切层的分辨率也非常高。
(3)前台阶问题
如图5所示为采用MPSI4格式计算前台阶问题(N＝960×320)。前台阶问题的初始条件为Mach数为3的强运动激波沿着带有台阶的风洞向右传播。台阶高度为0.2，前缘位置在0.6，计算域为[0,3]×[0,1]。初始物理参数为：
(ρ,u,v,p)＝(1.4,3.0,0.0,1.0)
左、右边界为开边界，其余边界取壁面反射条件。采用960×320的均匀网格计算截止到无量纲时间t＝4的时刻。
计算结果在激波附近都能达到基本无波动，对剪切层的分辨率也非常高。
(4)二维气动噪声传播问题
计算区域取为x∈[-100,100]，y∈[-100,100]，x方向的马赫数Mx＝0.5，y方向的马赫数My＝0，初始时刻的物理量为：
Figure RE-GDA0003640969000000172
Figure RE-GDA0003640969000000181
/>
Figure RE-GDA0003640969000000182
Figure RE-GDA0003640969000000183
采用MPSI4格式计算二维气动噪声传播问题，如图6所示为采用MPSI4 格式计算气动噪声传播问题；如图7所示为采用MPSI4格式计算气动噪声传播的噪声脉动分布。从计算结果与精确解的比较可以看出：MPSI4格式耗散性非常小，长时间气动噪声模拟能力很强。
本发明提供的一种高阶非线性激波捕捉空间离散方法，通过单调保持非线性样条插值的办法来求解半网格点处的函数值，从而得到了高阶非线性紧致格式MPSIk，所述高阶非线性紧致格式MPSIk保持了原有中心紧致差分格式的高阶精度、高分辨率和低耗散特性，同时由于其单调保持非线性插值的机制，使得高阶非线性紧致格式能够捕捉强激波，一系列数值算例表明所述样条插值非线性紧致格式是高阶精度、高分辨率和低耗散数值计算方法，具有广泛的科研以及应用前景。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种高阶非线性激波捕捉空间离散方法，其特征在于，包括如下步骤：
步骤S10：读取流场数据，得到各网格点上的物理通量和特征通量；
步骤S20：根据所述特征通量判断各网格点的间断性，所述间断性包括光滑和间断；
步骤S30：当所述网格点的间断性为光滑时，通过k次样条插值格式得到半网格点上通量，其中，k为样条插值格式的最高次数项，k为正整数且k≥3；当所述网格点的间断性为间断时，通过单调保持格式得到半网格点上通量；
步骤S40：根据所述半网格点上通量，并结合各网格点上的物理通量，得到整体的空间离散格式；
步骤S50：将时间推进至预设结束时间，得到预设结束时间的流场数据。
2.如权利要求1所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，步骤S20中，判断各网格点的间断性的方法为：同时满足
Figure FDA0003640968990000011
且/>
Figure FDA0003640968990000012
其中，j为网格点的编号，j为正整数，/>
Figure FDA0003640968990000013
为半网格点上通量，uj为第j个网格点的特征通量，uMP为初始单调保持通量，/>
Figure FDA0003640968990000014
为半网格点的光滑因子。
3.如权利要求2所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，所述初始单调保持通量uMP的计算公式为uMP＝uj+midmod[uj+1-uj，4(uj-uj-1)]，其中，uj+1为第j+1个网格点的特征通量，uj-1为第j-1个网格点的特征通量。
4.如权利要求3所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，所述半网格点的光滑因子
Figure FDA0003640968990000015
的计算公式为/>
Figure FDA0003640968990000021
其中，Mj+1为第j+1个网格点上k次样条插值函数的k-1次导数，Mj为第j个网格点上k次样条插值函数的k-1次导数，k次样条插值函数为S(j)，a≤j≤b，b-a小于预设区间，h为相邻网格点之间的距离。
5.如权利要求1-4之一所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，步骤S30中，当所述网格点的间断性为光滑时，通过三次样条插值格式得到半网格点上通量。
6.如权利要求5所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，半网格点上通量
Figure FDA0003640968990000022
的计算公式为/>
Figure FDA0003640968990000023
其中，/>
Figure FDA0003640968990000024
为三次样条插值。
7.如权利要求6所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，步骤S30中，当所述网格点的间断性为间断时，半网格点上通量
Figure FDA0003640968990000025
的计算公式为
Figure FDA0003640968990000026
其中，/>
Figure FDA0003640968990000027
为第j个半网格点特征通量的最小值，/>
Figure FDA0003640968990000028
为第j个半网格点的特征通量的最大值。/>
8.如权利要求1所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，步骤S10中所述特征通量的计算方法为：根据控制方程计算各网格点上的初始通量，通过对所述初始通量分裂得到正负通量，对所述正负通量进行特征投影，得到特征通量。
9.如权利要求1所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，步骤S40中所述空间导数通过线性紧致差分格式得到，所述紧致差分格式为
Figure FDA0003640968990000031
其中，a为第一系数，b为第二系数，α为第三系数，f′j-1为第j-1个网格点上物理通量的导数的差分近似式，f′j为第j个网格点上物理通量的导数的差分近似式，f′j+1为第j+1个网格点上物理通量的导数的差分近似式，/>
Figure FDA0003640968990000032
为第j个网格点的右半网格点上通量/>
Figure FDA0003640968990000033
的反投影，/>
Figure FDA0003640968990000034
为第j个网格点的左半网格点上通量/>
Figure FDA0003640968990000035
的反投影，fj+1为第j+1个网格点上物理通量，fj-1为第j-1个网格点上物理通量；
当j为网格内点的编号时，第三系数
Figure FDA0003640968990000036
此时紧致差分格式为四阶显式格式；当j为网格边界点的编号时，第三系数α＝0，此时紧致差分格式为四阶紧致格式；
所述第一系数a的计算公式为
Figure FDA0003640968990000037
所述第二系数b的计算公式为/>
Figure FDA0003640968990000038
10.如权利要求1所述的一种高阶非线性激波捕捉空间离散方法，其特征在于，步骤S50中采用三阶龙格库塔方法对时间项进行离散，将时间推进至预设结束时间。
CN202210086558.0A 2022-01-25 2022-01-25 一种高阶非线性激波捕捉空间离散方法 Active CN114638173B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210086558.0A CN114638173B (zh) 2022-01-25 2022-01-25 一种高阶非线性激波捕捉空间离散方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210086558.0A CN114638173B (zh) 2022-01-25 2022-01-25 一种高阶非线性激波捕捉空间离散方法

Publications (2)

Publication Number Publication Date
CN114638173A CN114638173A (zh) 2022-06-17
CN114638173B true CN114638173B (zh) 2023-06-02

Family

ID=81946138

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210086558.0A Active CN114638173B (zh) 2022-01-25 2022-01-25 一种高阶非线性激波捕捉空间离散方法

Country Status (1)

Country Link
CN (1) CN114638173B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116384288B (zh) * 2023-06-05 2023-08-25 中国空气动力研究与发展中心计算空气动力研究所 一种可压缩流动高分辨率数值模拟方法、介质及设备
CN116542184B (zh) * 2023-07-05 2023-09-19 中国空气动力研究与发展中心计算空气动力研究所 粘性通量的计算方法、装置、终端设备和存储介质
CN116611369B (zh) * 2023-07-17 2023-09-29 中国空气动力研究与发展中心计算空气动力研究所 基于光滑度量量级与候选模板点数的插值方法及装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107391885A (zh) * 2017-08-29 2017-11-24 西北工业大学 基于有限体积法的剪切滑移动网格方法
CN108197367A (zh) * 2017-12-27 2018-06-22 中国空气动力研究与发展中心计算空气动力研究所 一种基于流场通量阶跃的高精度间断Galerkin人工粘性激波捕捉方法
CN110457806A (zh) * 2019-08-02 2019-11-15 南京航空航天大学 基于交错网格的中心五阶weno格式的全流场模拟方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2004061723A1 (ja) * 2002-12-27 2004-07-22 Riken Ｖ−ｃａｄデータを直接用いた非圧縮性粘性流体の流れ場の数値解析方法と装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107391885A (zh) * 2017-08-29 2017-11-24 西北工业大学 基于有限体积法的剪切滑移动网格方法
CN108197367A (zh) * 2017-12-27 2018-06-22 中国空气动力研究与发展中心计算空气动力研究所 一种基于流场通量阶跃的高精度间断Galerkin人工粘性激波捕捉方法
CN110457806A (zh) * 2019-08-02 2019-11-15 南京航空航天大学 基于交错网格的中心五阶weno格式的全流场模拟方法

Also Published As

Publication number Publication date
CN114638173A (zh) 2022-06-17

Similar Documents

Publication Publication Date Title
CN114638173B (zh) 一种高阶非线性激波捕捉空间离散方法
Fidkowski et al. Review of output-based error estimation and mesh adaptation in computational fluid dynamics
CN108763683B (zh) 一种三角函数框架下新weno格式构造方法
Manceau et al. Inhomogeneity and anisotropy effects on the redistribution term in Reynolds-averaged Navier–Stokes modelling
Shitov et al. Flutter of rectangular simply supported plates at low supersonic speeds
CN112214869B (zh) 一种求解欧拉方程的改进型高阶非线性空间离散方法
CN115238397B (zh) 高超声速飞行器热环境计算方法、装置和计算机设备
Ceze et al. Drag prediction using adaptive discontinuous finite elements
Pandya et al. Assessment of USM3D Hierarchical Adaptive Nonlinear Method Preconditioners for Three-Dimensional Cases
CN109726465B (zh) 基于非结构曲边网格的三维无粘低速绕流的数值模拟方法
Glushko et al. Computational method for turbulent supersonic flows
Michal et al. Comparison of fixed and adaptive unstructured grid results for drag prediction workshop 6
Daxini et al. Parametric shape optimization techniques based on Meshless methods: A review
CN111780949B (zh) 基于cfd分析的高速进气道前体风洞实验总压修正方法
CN113361173B (zh) 一种完全基于当地流场参数对转捩模型可压缩修正的方法
Fernández Entropy-stable hybridized discontinuous Galerkin methods for large-eddy simulation of transitional and turbulent flows
Huang et al. Piston problem for the isentropic Euler equations for a modified Chaplygin gas
CN112577657A (zh) 一种分离激波振荡产生的脉动载荷快速预测方法
CN113378440A (zh) 一种流固耦合数值模拟计算方法、装置及设备
Galbraith et al. Comparisons of HPCMP CREATETM-AV Kestrel-COFFE, SU2, and MIT SANS RANS Solutions using Output-Based Adapted Meshes for a Multi-Element Airfoil
CN111563314B (zh) 一种七点weno格式的构造方法
Righi et al. Uncertainties Quantification of CFD-Based Flutter Prediction
Bhatia et al. Higher-order transonic flutter predictions
CN109726454B (zh) 管路系统的流固耦合建模方法及装置
Arian On the coupling of aerodynamic and structural design

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