CN112765725A - 针对多维Euler方程的解析黎曼解算方法 - Google Patents

Abstract

本发明公开了一种针对多维Euler方程的解析黎曼解算方法，基于流动的物理规律，将多维问题化为一维激波管问题，多维速度分量分解为产生波系的运动和随波系的运动，基于速度分解，在法向方向上应用一维的解析黎曼解，而切向速度通过波系保持不变，进而利用该物理规律将解析黎曼解方法推广到多维的通量计算，实现多维Euler方程的解析黎曼解数值求解。本发明提出一种针对多维Euler方程的解析黎曼解算方法，其相对于一般的矢通量分裂法相比，格式耗散小，对于粘性边界层计算而言，可以采用较少的网格数，或者较稀疏的网格分布，即可得到好的结果。

Description

针对多维Euler方程的解析黎曼解算方法

技术领域

本发明涉及一种计算流体力学领域。更具体地说，本发明涉及一种用在高超声速流场的计算流体力学(Computational Fluid Dynamics，CFD)数值求解时采用的针对多维Euler方程的解析黎曼解算方法。

背景技术

在航空航天领域，飞行器设计首先为飞行器外形设计，即满足气动力热特性要求。飞行器的气动力热特性特性可以通过风洞实验获得，也可以通过计算流体力学(Computational Fluid Dynamics，CFD)软件开展数值仿真计算得到。

高超声速流场的计算流体力学(Computational Fluid Dynamics，CFD)数值求解时，如果不考虑粘性效应，则可以采用无粘性的Euler方程求解。计算中为了捕捉高超声速流场中的激波，一般采用有限体积方法和高精度高分辨率算法，如MUSCL类型格式。MUSCL格式一般分为三个步骤实现：首先，进行带限制器的插值计算，获得网格单元边界的左右值；其次，根据边界上左右值，求解黎曼问题，计算边界上的通量；最后，根据边界上流入流出的通量，计算网格单元的参数变化。

第二步中根据黎曼问题求通量的方法有很多种，大多采用近似黎曼解算法，如矢通量分裂算法(Flux-Vector-Splitting Scheme，FVS)和通量差分裂算法(Flux-Difference-Splitting Scheme，FDS)算法等。研究表明(黎作武，力学学报，2008)，基于近似黎曼解算法的数值格式具有较大的数值耗散，不适合于计算边界层问题，如高超声速气动热等。

发明内容

本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。

为了实现根据本发明的这些目的和其它优点，提供了一种针对多维Euler 方程的解析黎曼解算方法，其特征在于，基于流动的物理规律，将多维问题化为一维激波管问题，多维速度分量分解为产生波系的运动和随波系的运动，基于速度分解，在法向方向上应用一维的解析黎曼解，而切向速度通过波系保持不变，进而利用该物理规律将解析黎曼解方法推广到多维的通量计算，实现多维Euler方程的解析黎曼解数值求解。

优选的是，基于网格内部的流动参数U的变化率等于四周边界处通量的净值，则通量的计算过程被配置为包括：

S1、流动参数U＝(r，u，v，p)分布分别代表流体密度、x方向速度、y 方向速度和压力，将网格中心点处的U值，插值到网格边界处；

S2、基于公式一、公式二将左边和右边参数中的速度分解为切向速度u_n和法向速度u_τ：

u_n＝u_nx+vn_y；

u_τ＝-u_y+vn_y；

S3、基于法向速度，按照一维方式求解黎曼解问题，并根据不同的波系结构，计算得到边界上的流动参数值

S4、根据波系中心处的法向速度和切向速度，根据公式三、公式四反向推处x方向的速度u^*和y方向的速度v^*：

S5、根据界面处的流动参数(ρ*,u*,v*,p*)计算网格边界处的通量，并通过累加各边界上的通量得到网格中心处的参数变化。

本发明至少包括以下有益效果：其一，本发明提出的基于解析黎曼解的数值解算方法，其相对于一般的矢通量分裂法相比，格式耗散小，对于粘性边界层计算而言，可以采用较少的网格数，或者较稀疏的网格分布，即可得到好的结果。

其二，本发明提出的基于解析黎曼解的数值解算方法，降低目前普遍采用的近似黎曼解算法中的数值耗散，达到高精度求解粘性边界层流动的目的。该方法不同于从数学方程出发的特征向量方法，而是基于物理流动图像而导出的，形式上更简单，计算效果是更为优异。

本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。

附图说明

图1为典型的二维网格格心型有限体积方法的原理示意图；

图2为二维网格格心型有限体积方法插值原理示意图；

图3为几种典型激波管流动结构示意图；

图4为本发明解析黎曼解的计算结果示意图；

图5为近似黎曼解(FVS)格式的计算结果示意图；

图6为球锥外形在计算外形时的气动热计算对比图；

图7为近似黎曼解计算结果曲线图；

图8为解析黎曼解计算结果曲线图；

图9为迎风面与实验结果对比图；

图10为背风面与实验结果对比l6图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

本发明提出一种针对多维Euler方程的基于解析黎曼解的数值解算方法，其的计算过程实际上是求解一维激波管的流场。如果流动是二维或三维的，根据数学方程，应该计算多维方程的特征值和特征向量，然后计算波系分解，这种计算过程是复杂的。在实际计算中，我们根据流动的物理规律，将多维问题化为一维激波管问题，多维速度分量分解为产生波系的运动和随波系的运动，则可以将问题大为简化，从而实现多维Euler方程的解析黎曼解数值求解。其相对于目前普遍采用的近似黎曼解算法中的数值耗散，达到高精度求解粘性边界层流动的目的。该方法不同于从数学方程出发的特征向量方法，而是基于物理流动图像而导出的，形式上更为简单。

基于此，本发明的解析黎曼解可以大大提高气动热计算的精度和可信性。其原因在于解析黎曼解的计算过程实际上是求解一维激波管的流场。如果流动是二维或三维的，根据数学方程，应该计算多维方程的特征值和特征向量，然后计算波系分解，这种计算过程是复杂的。在实际计算中，我们根据流动的物理规律，将多维问题化为一维激波管问题，多维速度分量分解为产生波系的运动和随波系的运动，则可以将问题大为简化，从而实现多维Euler方程的解析黎曼解数值求解。

具体来说，有限体积方法分为格点型和格心型两种类型，从数值方法来看其二者没有实质差别。以二维流动和格心型有限体积方法来说明本发明的计算步骤。

图1给出了一个典型的二维网格格心型有限体积方法的原理。网格中心点以圆圈表示，四周有流动分辨从该网格的边界处流入或流出(称为通量)，该网格内部的流动参数U的变化率等于四周边界处通量的净值，即：

通量E和F的算法决定了数值方法的差别，称为数值格式。我们采用的是基于解析黎曼解的算法。一般来说，解析黎曼解是一维流动条件下求出的，对于多维流动，需要发明一种以解析黎曼解为基础的算法。具体步骤如下：

以通量E的计算过程为例。流动参数U＝(ρ,u,v,p)分布代表流体密度、 x方向速度、y方向速度和压力。

Step1：将网格中心点处的U值，插值到网格边界处，一般来说，边界处左边和右边通过插值得到的量值不同，分别标记为UL和UR，如图2所示。 Step2：分别将左边和右边参数中的速度分解为切向速度和法向速度，假设网格边界的法向和切向单位向量为：

那么法向速度的计算方法为：

u_n＝u_nx+vn_y；

切向速度为：

u_τ＝-u_y+vn_y；

Step3：只考虑法向速度，按照一维方式求解黎曼解问题，也就是一维激波管问题，图3为几种典型激波管流动结构。

根据不同的波系结构，可以计算得到边界上的流动参数值，以上标“*”号表示。在滑移线两侧，法向速度和压力是连续的，而密度和切向速度是间断的。密度值可以通过两侧的激波关系式或膨胀波关系式得到，而切向速度则分别等于左右参数值的切向投影。这是因为根据流动物理机理，切向速度穿过激波或膨胀波，其大小不变。

Step4：根据波系中心处的法向速度和切向速度，可以反推处x方向的速度和y方向的速度，即

Step5：根据界面处的流动参数(ρ*,u*,v*,p*)计算通量即可。如

E＝(ρ*u*,ρ*u*u*+p*,ρ*u*v*,ρ*e*)；

通过上述5个步骤，即可计算通过网格边界处的通量。最终通过累加各边界上的通量即可得到网格中心处的参数变化。

本发明基于速度分解，在法向方向上应用一维的解析黎曼解黎曼解，而切向速度通过波系不变，利用该物理规律将解析黎曼解方法推广到多维的通量计算。故采用基于解析黎曼解的数值方法与一般的矢通量分裂法相比，格式耗散小，对于粘性边界层计算而言，可以采用较少的网格数，或者较稀疏的网格分布，即可得到好的结果。比如平板边界层速度计算结果，采用解析黎曼解的数值方法只需要5个网格点，即可得到与理论值一致的结果，而矢通量分裂法(近似黎曼解)需要的网格数超过20个。

以图4中本发明的解析黎曼解计算结果图5中的近似黎曼解(FVS)格式的计算结果进行比较，可知本发明的解析黎曼解计算的结果良好。

采用解析黎曼解方法求解高超声速飞行器的热环境具有更高的可信性。气动热计算非常依赖于网格分布的密度，不同的网格密度得到的气动热结果不同，具有很大的散布度，这称为计算结果具有网格相关性。如果采用解析黎曼解方法，则可以大大降低这种网格相关性。已有的计算实践表明，与近似黎曼解相比，气动热计算的网格密度可降低十倍。

图6是一个球锥外形的气动热计算对比。一般情况下，如果采用近似黎曼解，网格尺度需要小于0.00005，如果以网格雷诺数(即雷诺数乘以网格间距Re·Δy)衡量，网格雷诺数衡量在5以内才能得到如图7基本一致的计算结果，而采用解析黎曼解的计算方法，网格间距小于0.002就可以获得如图8 基本一致的结果，以行业术语而言，也就是网格收敛结果，或网格无关解。与近似黎曼解相比，网格间距可以大一个量级以上。

根据一个类似航天飞机外形的计算结果及计算与实验数据的比较(图9- 图10)。可以看出，计算与实验数据符合很好，充分说明本发明能够很好地解决工程实际问题，部分替代风洞实验结果。

以上方案只是一种较佳实例的说明，但并不局限于此。在实施本发明时，可以根据使用者需求进行适当的替换和/或修改。

这里说明的设备数量和处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言，可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (2)

1.一种针对多维Euler方程的解析黎曼解算方法，其特征在于，基于流动的物理规律，将多维问题化为一维激波管问题，多维速度分量分解为产生波系的运动和随波系的运动，基于速度分解，在法向方向上应用一维的解析黎曼解，而切向速度通过波系保持不变，进而利用该物理规律将解析黎曼解方法推广到多维的通量计算，实现多维Euler方程的解析黎曼解数值求解。

2.如权利要求1所述的针对多维Euler方程的解析黎曼解算方法，其特征在于，基于网格内部的流动参数U的变化率等于四周边界处通量的净值，则通量的计算过程被配置为包括：

S1、流动参数U＝(r，u，v，p)分布分别代表流体密度、x方向速度、y方向速度和压力，将网格中心点处的U值，插值到网格边界处；

S2、基于公式一、公式二将左边和右边参数中的速度分解为切向速度u_n和法向速度u_τ：

u_n＝u_nx+vn_y；

u_τ＝-u_y+vn_y；

S3、基于法向速度，按照一维方式求解黎曼解问题，并根据不同的波系结构，计算得到边界上的流动参数值

S4、根据波系中心处的法向速度和切向速度，根据公式三、公式四反向推处x方向的速度u^*和y方向的速度v^*：

S5、根据界面处的流动参数(ρ*,u*,v*,p*)计算网格边界处的通量，并通过累加各边界上的通量得到网格中心处的参数变化。

Images