CN106682262A - 一种获取飞行器流场的数值模拟方法 - Google Patents

Abstract

本发明涉及一种获取飞行器流场的数值模拟方法，采用基于非结构直角网格的间断Galerkin方法来求解流场，采用四叉树结构的直角网格将待求解的流场区域离散成网格的集合；对于每个四边形网格，利用双线性坐标变换，将目前物理坐标系下每个网格映射到计算坐标系下一个规则网格上，获得标准正方形网格；通过数值计算获得计算坐标系下单元内守恒变量；确定物理坐标系下每个网格的上面、下面、左面和右面；根据单元交界面左右变量的差别，构造了间断探测器；将每个网格内的守恒变量在离散网格上进行流场显示，获得流场分布。本发明将基于非结构直角网格的间断Galerkin方法引入到流场计算中，提高了流场计算的精度，能够求解复杂流动现象。

Description

一种获取飞行器流场的数值模拟方法

技术领域

本发明涉及一种获取飞行器流场的数值模拟方法，属于计算流体动力学领域。

背景技术

计算流体动力学软件的结构一般由前处理、求解器和后处理三部分组成，前处理主要包括几何模型构建和划分网格；求解器包括控制方程的离散化及选用数值计算方法；后处理包括速度场、温度场、压力场及其它流场参数的计算机可视化及动画处理。

相比于传统的结构网格离散方法，非结构直角网格优越的几何灵活性可以适应任意复杂外形，并且其随机的数据结构非常利于进行网格自适应。但是，现有流场求解器以二阶精度的计算格式为基础，尽管其已经在复杂外形流动数值模拟中取得了较大的成功，但是二阶精度计算格式具有较大的数值耗散与数值色散，对一些复杂的流动现象，如分离流动、气动噪声、湍流等问题，仍然难以给出精细的流场结构，尤其对于湍流的大涡模拟、直接数值模拟，必须采用高阶精度格式。因此发展基于非结构网格的高精度计算方法十分重要。

在众多的高阶精度计算方法中，间断Galerkin(DG)有限元方法备受关注。DG方法保持了传统有限元方法和有限体积法的优点，融入了有限体积法中的数值通量、黎曼解、TVD和限制器等思想和概念，能够有效地解决含间断等问题；DG方法对计算网格没有特殊的限制，非常适用于复杂几何外形的数值模拟。并且通过选取适当的基函数，可以处理含有悬空节点的计算网格。DG方法具有良好的紧致特性，DG方法精度的提高可以通过提高单元内多项式的次数来实现，这克服了有限体积法通过扩大网格模板重构单元交接面处变量来提高精度的不足，此外，DG方法具有良好的守恒性、稳定性和收敛性。现有方法计算流场的精度较低，如何利用DG方法计算流场，提高计算精度是本领域亟待解决的技术问题。

发明内容

本发明的目的在于克服现有技术的不足，提出一种获取飞行器流场的数值模拟方法，该方法能够达到高于二阶的计算精度，有效地减小数值耗散和色散。

本发明目的通过如下技术方案予以实现：

提供一种获取飞行器流场的数值模拟方法，包括以下步骤：

(1)采用四叉树结构的直角网格将待求解的流场区域离散成网格的集合,包括外围直角网格和物面处的贴体网格，网格类型均为四边形；

(2)对于每个四边形网格，利用双线性坐标变换，将目前物理坐标系下每个网格映射到计算坐标系下一个规则网格上，获得标准正方形网格；

(3)利用间断Galerkin方法对流动方程进行离散求解，通过数值计算获得计算坐标系下单元内流动守恒变量 为权系数，w_i(ξ,η)为高斯积分的正交基；

(4)确定物理坐标系下每个交接面是网格单元的上面、下面、左面和右面中的哪一个面；

(5)通过每个面左右单元内守恒变量求得该面的左侧守恒变量u_l和该面的右侧守恒变量u_r；

(6)根据单元交接面左、右侧守恒变量的差别，构造间断探测器，第i个网格的探测因子：如果D_i<θ，则确定该网格的内守恒变量如果D_i>θ，则通过限制器对进行限制，重新计算该网格的内守恒变量h为网格的尺度，p为高斯多项式的次数，为左右单元守恒变量的平均值；

(7)将每个网格内守恒变量在离散网格上进行流场显示，获得流场分布。

优选的，步骤(4)中确定物理坐标系下每个网格的上面、下面、左面和右面的具体方法为：

(4.1)对于外围直角网格，定义每个交接面的指向，与水平轴垂直的面指向水平轴正方向，由该交接面的左单元指向右单元，对于右单元来说，该面为右单元的左面，对于左单元来说，该面为左单元的右面；与竖直垂直的面指向竖直轴正方向，由该单元面的左单元指向右单元，对左单元来说，该面为左单元的上面，对于右单元来说，该面为右单元的下面；

(4.2)当粘性层数为1时，对于与外围直角网格相接的粘性层网格，当该粘性层网格被确定为左单元时，与外围直角网格相接的面为该左单元的右面，按照右手法则确定其他三个面，即按照逆时针方向依次为右面、上面、左面和下面；当粘性层数大于1时，将母网格整体作为与外围直角网格相接的粘性层网格，该母网格内的各个子网格内各面的排序均与母网格相同。

优选的，步骤(5)中求得该面的左侧守恒变量u_l和该面的右侧守恒变量u_r具体方法为：左单元令ξ＝1,η＝0，右单元令ξ＝-1,η＝0，

优选的，步骤(6)中通过限制器对进行限制，重新计算该网格的内守恒变量的具体方法为：

(6.1)通过公式限制

如果同时成立，则停止；否则进入步骤(6.2)

(6.2)限制 如果则停止；否则进入步骤(6.3)；

(6.3)限制

其中上标l表示为其左单元，上标r表示为其右单元，上标u表示为其上单元，上标b表示为其下单元，α、β为控制耗散的参数。

本发明与现有技术相比具有如下优点：

(1)本发明将基于非结构直角网格的间断Galerkin方法引入到流场计算中，提高了流场计算的精度，能够求解复杂流动现象。

(2)采用间断Galerkin方法重构交接面上的流动变量时，设定了潜在的数据规则，快速确定各面对应关系。通过这种方法可以快速确定交接面的对应关系，有效地提高计算效率和减小存储量。

(3)引入限制器来抑制计算格式在间断附近的数值振荡，为了防止限制器在解的光滑区域仍起作用，从而降低了格式捕捉精细流场的能力，构造了间断探测器。

(4)采用本发明的方法获得的流场精度高于二阶精度，通过后处理进行高效的流场显示。

附图说明

图1(a)为锯齿形内面网格生成示意图；图1(b)为光顺后内面网格生成示意图；图1(c)粘性层网格生成示意图；

图2为物理坐标系向计算坐标系转换的示意图；

图3为物面粘性层单元交接面的对应关系示意图；

图4(a)为M6翼型算例的非结构直角网格分布；图4(b)为M6翼型算例计算所得的压力等值线图。

具体实施方式

本发明获取飞行器流场的数值模拟方法，采用基于非结构直角网格的间断Galerkin方法来求解流场，具体步骤如下：

1.采用四叉树结构的直角网格将求解域离散成单元的集合

首先根据选择的背景网格和设置参数，在计算域内划分出初始网格，之后对物体附近的单元进行多层四叉树加密，然后删除物体内部以及与物体相交的网格单元产生锯齿形内面，留下的网格称为外围的直角网格；然后对内面进行光顺并向物面投影，生成贴体的粘性层网格，对贴体粘性层网格(母网格)分割生成层数和间距相同的子网格，最终获得带有粘性层的直角网格，所有网格类型均为四边形单元，可以利用网格生成软件生成直角网格，生成过程参见图1。

2.对于每个四边形网格，利用双线性坐标变换，将目前物理坐标系下所有网格统一映射到计算坐标系下一个规则网格上，获得标准正方形网格。

在自然坐标系内积分限是规格化的，希望能在自然坐标系下按照规格化的数值积分方法求解积分表达式，为此需要建立两个坐标系之间的导数、体积微元和面积微元之间的变换关系。对于四边形单元，利用双线性坐标变换，可以将物理坐标系xoy下所有网格统一映射到计算坐标系ξO'η下一个规则闭区域S_i＝{(ξ,η)|-1≤ξ≤1,-1≤η≤1}，如图2所示。

两个坐标间导数变换关系式及雅可比行列式为|J|：

以左下角的点为1点，逆时针依次为2、3、4点，其中x_i，y_i为第i个点物理坐标系，N_i坐标变换的型函数为：

3.利用间断Galerkin方法对流动控制方程进行离散求解，通过数值计算获得计算坐标系下单元内守恒变量

下面对二维流动控制方程Euler方程为例，利用间断Galerkin方法对控制方程进行推导：

其中U为守恒变量，ρ为网格中心点处的密度，u为x方向的速度，v为y方向的速度，E为能量密度，t为时间，F、G为数值通量， 其中p为压力。

通过坐标变换到计算坐标系可以得到：

在单元S_e上对控制体进行积分：

其中w_j为积分参数，表示求取梯度，S表示体积分单元，S_e为标准单元格的区域，为计算坐标系下的流动守恒变量，为待求量，

得到等效积分弱形式为：

其中Γ线积分单元，为积分单元S_e的边界。

在计算坐标系上，令：

将公式(2)带入公式(1)获得

通过数值方法求解公式(3)即可得到的值

选择一组正交基：

替换w_i

利用Gauss公式对公式(3)进行求解，采用三点积分公式，积分点为0,

对于线积分项采用进行计算，其中H_i为权系数，值为5/9,8/9,5/9；

对于面积分项采用进行计算，H_i、H_j均权系数，值为5/9,8/9,5/9；

求得的值，进而通过公式(2)计算获得计算坐标系下单元内守恒变量

4.为了确定物理坐标系下交接面是相邻单元的上面、下面、左面、右面中的哪一个面，设定了潜在的数据规则，方便计算搜索，如图3所示。

(a)对于外围直角网格，定义每个交接面单元的指向，与水平轴(X轴)垂直的面指向X轴正方向，由该单元面的左单元指向右单元，对于右单元来说，该面为右单元的左面，对于左单元来说，该面为左单元的右面；与竖直轴(Y轴)垂直的面指向Y轴正方向，由该单元面的左单元指向右单元，对左单元来说，该面为左单元的上面，对于右单元来说，该面为右单元的下面。

(b)对于贴体的粘性层网格，当该网格被确定为左单元时，与外围直角网格相接的面为该网格单元的右面，按照右手法则确定其他三个面，即按照逆时针方向依次为右面、上面、左面和下面；当该网格被确定为右单元时，与外围直角网格相接的面为该网格单元的左面，按照右手法则确定其他三个面，即按照逆时针方向依次为左面、下面、右面和上面；

当粘性层数为1时，按上述方法确定各面；当粘性层数大于1时，母网格内的各个子网格内各面的排序均与母网格相同。通过这种方法可以快速确定交接面的对应关系，有效地提高计算效率和减小存储量。

5.通过左右单元内守恒变量求得该面的左侧守恒变量u_l和该面的右侧守恒变量u_r，具体方法为：左单元令ξ＝1,η＝0，右单元令ξ＝-1,η＝0，

6.根据单元交界面左右变量的差别，构造了一种间断探测器。该间断探测器的构造原理简单，编程实现容易。该间断侦测器能够准确捕捉到激波的位置，从而只在“间断单元”引入限制器，能够有效地减少计算量、保证解在光滑区的计算精度。

定义第i个网格的探测因子：

h为控制体的特征长度,分段p次多项式,为左右单元守恒变量的平均值。

其中为左右单元守恒变量的平均值。以此来确定“问题单元”。

θ为侦测阀一般选取θ∈[0.5，1.0]，选用密度作为侦测的对象。

moment限制器保证格式的鲁棒性,避免产生数值振荡。moment限制器属于minmod类型的限制器。minmod限制器的工作原理是：如果一个单元的解超过了其相邻单元解的平均值，就减小这个单元解的斜率。moment限制器与此类似，当第k单元数值解的第i阶导数应落在第i-1阶导数与前后两单元i-1阶导数梯度范围之内。在不存在间断的区域,解的泰勒展开式说明了系数与解的阶导数是相关的,因而限制解的系数意味着限制了解的导数。

计算坐标下三阶精度解的表达式中存在ξ、η的交叉项,因此必须对系数进行两个方向的限制。网格单元(k,m)上解的限制器公式如下:

其中分别为基函数的系数。其中上标l标识为其左单元，上标r表示为其右单元，上标u表示为其上单元，上标b表示为其下单元，这里的α、β为控制耗散的参数，它们的值越小，限制的强度越大，耗散也就越大，计算则更稳定。数值实验表明取α∈[1，10]，β∈[1，2]是合适的。

Moment限制器的实行步骤为：

1)限制如果同时成立，则停止；否则进入步骤2)；

2)限制如果则停止；否则进入步骤3)；

3)限制

7.将每个网格内守恒变量在离散网格上进行流场显示，获得流场分布。。

参见图4(a)为M6翼型算例的非结构直角网格分布，图4(b)为M6翼型算例计算所得的压力等值线图，从图中可以看出高精度DG方法能够得到高于二次分布的流场信息，在网格内布置一定数量的节点，并输出其上的流场信息，经后处理可以得到高效的流场显示。

以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (4)

1.一种获取飞行器流场的数值模拟方法，其特征在于，包括以下步骤：

(1)采用四叉树结构的直角网格将待求解的流场区域离散成网格的集合,包括外围直角网格和物面处的贴体网格，网格类型均为四边形；

(2)对于每个四边形网格，利用双线性坐标变换，将目前物理坐标系下每个网格映射到计算坐标系下一个规则网格上，获得标准正方形网格；

(3)利用间断Galerkin方法对流动方程进行离散求解，通过数值计算获得计算坐标系下单元内流动守恒变量 为权系数，w_i(ξ,η)为高斯积分的正交基；

(4)确定物理坐标系下每个交接面是网格单元的上面、下面、左面和右面中的哪一个面；

(5)通过每个面左右单元内守恒变量求得该面的左侧守恒变量u_l和该面的右侧守恒变量u_r；

(6)根据单元交接面左、右侧守恒变量的差别，构造间断探测器，第i个网格的探测因子：如果D_i<θ，则确定该网格的内守恒变量如果D_i>θ，则通过限制器对进行限制，重新计算该网格的内守恒变量h为网格的尺度，p为高斯多项式的次数，为左右单元守恒变量的平均值；

(7)将每个网格内守恒变量在离散网格上进行流场显示，获得流场分布。

2.根据权利要求1所述的一种获取飞行器流场的数值模拟方法，其特征在于：步骤(4)中确定物理坐标系下每个网格的上面、下面、左面和右面的具体方法为：

(4.1)对于外围直角网格，定义每个交接面的指向，与水平轴垂直的面指向水平轴正方向，由该交接面的左单元指向右单元，对于右单元来说，该面为右单元的左面，对于左单元来说，该面为左单元的右面；与竖直垂直的面指向竖直轴正方向，由该单元面的左单元指向右单元，对左单元来说，该面为左单元的上面，对于右单元来说，该面为右单元的下面；

(4.2)当粘性层数为1时，对于与外围直角网格相接的粘性层网格，当该粘性层网格被确定为左单元时，与外围直角网格相接的面为该左单元的右面，按照右手法则确定其他三个面，即按照逆时针方向依次为右面、上面、左面和下面；当粘性层数大于1时，将母网格整体作为与外围直角网格相接的粘性层网格，该母网格内的各个子网格内各面的排序均与母网格相同。

3.根据权利要求1所述的一种获取飞行器流场的数值模拟方法，其特征在于，步骤(5)中求得该面的左侧守恒变量u_l和该面的右侧守恒变量u_r具体方法为：左单元令ξ＝1,η＝0，右单元令ξ＝-1,η＝0，

4.根据权利要求1所述的一种获取飞行器流场的数值模拟方法，其特征在于，步骤(6)中通过限制器对进行限制，重新计算该网格的内守恒变量的具体方法为：

(6.1)通过公式限制

${{\tilde{u}}_6}^{\&prime;}=\min \mathrm{mod}({\tilde{u}}_6,\mathrm{\&alpha;}({\tilde{u}}_3^u-{\tilde{u}}_3),\mathrm{\&alpha;}({\tilde{u}}_3-{\tilde{u}}_3^b\left)\right)$

${{\tilde{u}}_5}^{\&prime;}=\min \mathrm{mod}({\tilde{u}}_5,\mathrm{\&alpha;}({\tilde{u}}_2^r-{\tilde{u}}_2),\mathrm{\&alpha;}({\tilde{u}}_2-{\tilde{u}}_2^l\left)\right)$

如果同时成立，则停止；否则进入步骤(6.2)

(6.2)限制如果则停止；否则进入步骤(6.3)；

(6.3)限制

其中上标l表示为其左单元，上标r表示为其右单元，上标u表示为其上单元，上标b表示为其下单元，α、β为控制耗散的参数。