一种流体问题中物理量求解方法

本发明属于流体力学计算，具体涉及一种流体问题中物理量求解方法。

背景技术：

1、高精度数值算法具备较低的数值耗散，相同条件下相比于低阶格式能够得到更为精细的流动结构，所以一直是计算流体力学、计算数学、计算电磁学、计算声学等领域的重点研究课题。在众多的高精度算法中，weno格式扮演着非常重要的角色。weno格式于1994年首次被liu等人提出，并在1996年被jiang和shu完善并给出了任意高阶weno格式的统一框架。基于该框架，众多学者展开了诸多的后续研究。比如：balsara等人构造了7-13阶精度的保单调weno格式；henrick等人设计了weno-m格式用于恢复传统weno格式在极值点附近的数值精度；dumbser等人发展了非结构网格体系下的有限体积weno方法，并通过tsoutsanis等人开发的ucns3d代码实现了非结构weno格式的工程化应用。上述weno格式构造过程相对比较复杂，涉及繁琐的线性权计算过程。因此，zhu等人于2016年设计了一类称之为非等距weno的格式。该格式的线性权值可以任取一组和为一的正数，极大地简化了实施过程(尤其是非结构网格体系下)，并且对于定常问题具备更加优异的定常收敛性。

2、毫无疑问的是weno格式是一类受欢迎的激波捕捉类高精度格式，具备极好的工程应用前景。但是，由于其构造过程中涉及特征分解以及复杂的非线性权值计算，从而导致其计算效率不高。因此，如何继承weno格式优良特性的同时降低其计算成本一直是weno格式的研究热点。杂交weno格式是一种可行的替代方案，其通过间断识别器对网格单元进行分类，然后在光滑区域使用廉价高效的线性方法，而在间断区域使用昂贵但激波捕捉能力优异的weno格式。这种结合使得数值算法能够充分发挥线性方法和非线性weno格式各自的优点。但是正如pirozzoli所言，在杂交weno格式的过程中，对网格单元进行分类的间断识别器这一角色显得尤为重要，如何能够快速、准确、高效且无需人为参与的区分光滑流场单元和间断流场单元一直是一个开放的研究课题。

技术实现思路

1、本发明所要解决的技术问题是针对上述现有技术的不足，提供一种流体问题中物理量求解方法，设计一种无需人为参与且准确高效的间断识别器，然后对计算区域内的网格单元进行分类，分别利用廉价高效的线性方法和激波捕捉能力优异的非等距weno格式对流体力学控制方程进行离散，最后利用具备tvd性质的三阶runge-kutta时间离散方法实现控制方程在时间方向上的推进，进而得到任意时刻计算区域内所有网格点处的物理量，并通过无粘和黏性可压缩流动问题验证方法的有效性。

2、为实现上述技术目的，本发明采取的技术方案为：

3、一种流体问题中物理量求解方法，包括以下步骤：

4、步骤一、在实际流动问题计算区域内，生成直角坐标系并进行网格剖分；

5、步骤二、基于剖分后的计算区域，得到描述流体运动的euler或navier-stokes方程组的空间离散化方程；

6、步骤三、基于所述空间离散化方程，对计算区域内的网格点进行分类并标记，得到光滑网格点和间断网格点；

7、步骤四、分别在光滑网格点处使用线性格式和间断网格点处使用非等距有限差分weno格式对其空间导数项进行空间离散；

8、步骤五、基于步骤四，采用经典的runge-kutta时间离散方法对流体力学控制方程组进行时间方向上的推进，从而得到任意时刻网格离散点处的物理量值，进而数值模拟出实际流动问题中各个时刻的状态。

9、为优化上述技术方案，采取的具体措施还包括：

10、上述的实际流动问题包括前台阶反射问题、rayleigh-taylor不稳定问题、粘性激波管问题、三维激波-剪切层相互作用问题和taylor-green涡演化问题。

11、上述的步骤一中，对计算区域构建笛卡尔直角坐标系，并进行网格剖分，相应的网格点坐标记为ii,j,k＝(xi,yj,zk)，其中x，y，z表示空间变量，下标i，j，k表示网格点序号。

12、上述的步骤二中，描述流体运动的euler或navier-stokes方程组为：

13、

14、其中，t表示时间变量，x＝(x,y,z)表示空间变量向量，u(x,t)＝(ρ,ρu,ρv,ρw,e)t表示守恒变量向量，t表示转置；

15、表示总能，ρ、u、v、w、e都是关于时间变量t和空间变量x，y，z的函数，分别表示流体密度、x方向速度、y方向速度、z方向速度以及内能；

16、p＝(γ-1)ρe表示压力，其中γ表示比热比常数；

17、▽表示散度算子，σ为euler方程和navier-stokes方程的转换符号，σ＝0表示euler方程，σ＝1表示navier-stokes方程；

18、re为reynolds数，表示对流通量向量，表示黏性通量向量，上标x，y，z表示对应的空间方向，下标c和v表示对流和黏性；

19、fc(u)和fv(u,▽u)中各分量和的具体表达式为

20、

21、

22、

23、

24、

25、

26、其中，τxx、τyy、τzz分别表示x，y，z方向的正应力，τxy、τxz、τyz分表表示xy、xz、yz平面的切应力，qx、qy、qz表示x，y，z方向的热通量；

27、基于空间直角坐标系下的剖分网格，可以得到如下空间离散方程

28、

29、其中，ui,j,k表示守恒变量u(x,t)在网格点(xi,yj,zk)处的离散值，表示x方向对流通量在网格边界(xi±1/2,yj,zk)处的离散值，表示y方向对流通量fcy在网格边界(xi,yj±1/2,zk)处的离散值，表示z方向对流通量在网格边界(xi,yj,zk±1/2)处的离散值，表示x方向黏性通量在网格边界(xi±1/2,yj,zk)处的离散值，表示y方向黏性通量在网格边界(xi,yj±1/2,zk)处的离散值，表示z方向黏性通量在网格边界(xi,yj,zk±1/2)处的离散值。

30、上述的步骤三中，对计算区域内网格点进行分类并标记的具体过程如下：

31、步骤31、对于目标单元，即目标网格点ii,j,k，先固定y方向和z方向，即下标j和k不变，选取x方向左右相邻的五个单元组成模板得到一个四次多项式p1(x)，且满足如下条件：

32、

33、其中，δx表示x方向上的空间网格步长；

34、

35、其中，表示基函数，an,n＝0,1,2,3,4表示多项式系数且表达式如下

36、

37、

38、

39、

40、

41、步骤32、求多项式p1(x)的一次导数得到导数多项式p1′(x)：

42、

43、其中，上标'表示导数符号；

44、步骤33、求出导数多项式p1′(x)的最大值和最小值的弱形式；

45、步骤34，根据最大值和最小值的弱形式对所有的网格剖分点进行分类；

46、如果满足(max1±δx)与(min1±δx)异号或者(max2±δx)与(min2±δx)异号，则该网格点标记为含间断点的间断网格点，记为flag＝1；反之若同号则标记为光滑网格点，记为flag＝0；

47、步骤35，分别固定x方向和z方向以及x方向和y方向，重复步骤31-34，即可对所有单元进行分类和标记。

48、上述的步骤33中，导数多项式中的在区间[xi-5/2,xi]和[xi,xi+5/2]内分别具备单调性，所以其最大值和最小值在其区间端点处取得，具体的：

49、在区间[xi-5/2,xi]内有

50、

51、

52、在区间[xi,xi+5/2]有

53、

54、

55、上述的步骤四具体如下：

56、步骤41、固定y方向和z方向，即下标j和k不变，对于光滑网格点，采用线性格式得到网格边界处对流通量的五阶近似值，有

57、

58、

59、步骤42、对于间断网格点，采用非线性非等距weno格式得到流通量的五阶近似值；

60、步骤43，分别固定x方向和z方向以及x方向和y方向，重复步骤41-步骤42，即可实现所有网格单元对流项的空间离散。

61、上述的步骤42包括：

62、a)首先已知多项式p1(x)，对于目标单元ii,j,k及其左右两个网格点，另选取两个模板s(2)＝{ii-1,ii}和s(3)＝{ii,ii+1}，可以分别得到两个线性多项式p2(x)和p3(x)，满足如下条件

63、

64、

65、具体表达式分别如下：

66、

67、

68、b)任意选取三个和为1的正数作为线性权，记为γ1，γ2，γ3；

69、c)基于已有的三个多项式p1(x)，p2(x)和p3(x)，通过如下式子分别计算其对应的光滑因子：

70、

71、其中，n表示对应多项式的序号，α表示求和指标，r表示对应多项式的次数，表示多项式pn(x)对自变量x求α次导数；

72、d)基于线性权和光滑因子计算非线性权，如下

73、

74、其中，τ和为计算中的过渡值，ε＝10-6为防止分母为零的小量；

75、e)通过非等距weno格式得到网格边界处对流通量的五阶非线性近似值：

76、

77、

78、上述的步骤五中，假设时间步长为δt，tn表示第n时间层，表示在第n时间层网格点(xi,yj,zk)处u(xi,yj,zk,tn)的值，在tn时间层通过步骤四得到流体力学控制方程的空间导数离散，且记为然后利用满足tvd性质的runge-kutta时间离散方法离散半离散有限差分格式，从而得到时空全离散有限差分格式，如下：

79、

80、其中，和为计算中间的过渡值，空间方向和时间方向的离散完成后得到的时空全离散有限差分格式为关于时间层的迭代公式，已知初始流动状态的密度ρ、x方向速度u、y方向速度v、z方向速度w、内能e，然后根据迭代公式不断循环，即可得到流场在某一时刻或稳定状态时的物理量的数值近似值。

81、本发明具有以下有益效果：

82、针对已有的非等距weno格式能够实现可压缩流动问题的准确模拟，但其计算代价相对较大的问题，本发明的杂交非等距weno格式能够继承非等距weno格式本身所有的优良特性，同时大幅度提升其计算效率(粗略统计提升30％左右)。

83、此外，本发明通过新型杂交策略对计算网格点进行分类，从而在光滑区域使用高效的线性格式用于模拟流场的精细化结构，在间断区域使用具备基本无振荡特性的非等距weno格式模拟间断结构，进而实现可压缩流动问题的准确、高效、鲁棒求解。最后，本发明通过各种无粘和粘性可压缩流动问题验证其在流体力学问题模拟中的优越性。

84、本发明发展的间断识别器相比于已有的研究成果具备如下优点：

85、1)其减小了网格依赖性，其构造过程只跟目标单元的重构多项式有关，无需考虑网格类型；

86、2)其不包含任何人为参数，无需根据具体问题调整参数，极大的提升了其实用性；

87、3)其能够自动、准确、高效的区分出光滑单元和间断单元，进而很好地结合线性方法和非线性格式。