CN103413060B - 基于双控制体的格心网格数据三维激波特征定位方法 - Google Patents

Abstract

本发明公开了一种基于双控制体的格心网格数据三维激波特征提取方法，要解决的技术问题是解决目前三维计算流场格心网格数据激波特征提取准确性较差的问题。技术方案是采用Roe平均方法得到所有内部面中心处的压力值；采用高斯-格林公式得到各网格单元的压力梯度；对输入格点网格数据进行采样，得到一级采样点；判断一级采样点位置，若其所在面不是内部面，利用正则马赫数计算公式得到一级采样点处的最终正则马赫数，若一级采样点所在面是内部面则采用双控制体方法得到一级采样点处的最终正则马赫数。采用本发明可准确提取三维计算流场格心网格数据激波特征，提高复杂数值计算流场特征可视化的准确性。

Description

基于双控制体的格心网格数据三维激波特征定位方法

技术领域

本发明涉及三维激波特征定位方法，尤其是基于双控制体的格心网格数据三维激波特征定位方法。

背景技术

在计算流体力学领域，有限体积法(FiniteVolumeMethod，FVM)具有明确的物理含义，且计算效率高，是近年来发展十分迅速且广泛应用的数值计算方法。FVM常将数值解存储在网格单元中心，将单元看成是围绕中心的控制体，产生格心格式网格数据。

激波是一种十分重要的流动间断现象，其位置和强度对飞行器设计至关重要，只有成功定位激波特征区域，才能设计出安全性能满足要求的飞行器。激波是由于超音速气流遇到障碍物时产生的，在实际流场中，激波是非理想气流中厚度极薄(10^-5mm量级)的三维曲面。

目前已有的三维激波特征定位方法(如边缘检测法、梯度阈值法、兰金-于戈尼奥条件法、梯度方向法、两级采样法等)只适用于格点网格数据，无法直接定位格心网格流场数据中的激波特征。对于格心网格数据，采用现有的三维激波特征定位方法需要将其外推为格点网格数据，损失了数据精度。由于激波区域极薄，外推计算是一种平均运算，可能会抹平原始数据场中的间断和尖峰，使外推后的格点数据值域是原始格心数据的真子集，又由于激波特征区域通常极薄，这种外推的平均操作往往会抹平激波区域，使得格心网格数据的三维激波特征定位准确性较差。两级采样法(马千里；李思昆；曾亮.基于两级采样的非结构化网格流场多激波特征可视化方法.计算机研究与发展.2012,49(7):1450-1459.)是目前最佳的激波定位方法，该方法具体步骤如下。

步骤1：利用体积加权外推法或反转距离外推法将格心网格数据转化为格点网格数据；

步骤2：对格点网格数据进行采样，得到的采样点叫做一级采样点，基于网格顶点压力值确定一级采样点处的压力梯度，根据压力梯度确定一级采样点处的正则马赫数；

步骤3：根据正则马赫数对一级采样点进行初步筛选：满足正则马赫数大于等于1的采样点作为备选激波点(所述激波点指位于激波特征区域内的采样点)；

步骤4：通过二级采样点判断步骤3得到的备选激波点是否为真实激波点；

步骤4.1：确定备选激波点处的切平面及切平面法线，该法线及其反向延长线与一系列网格单元相交，产生一系列交点，选取在备选激波点两侧的两个交点(所述交点即二级采样点)；

步骤4.2：确定两个二级采样点处的速度；

步骤4.3：计算步骤4.2两个二级采样点处速度在切平面上的投影速度的大小；

步骤4.4：如果两个投影速度大小相同，则该备选激波点为真实激波点，否则为伪激波点(噪声)。

该方法得到一系列激波点，这些激波点组成的区域就是激波特征区域。

对于格点网格数据，两级采样法可以比较准确地定位激波特征。但是对于格心网格数据，两级采样法需要将格心网格数据外推为格点网格数据，然后进行激波定位。外推计算可能会抹平激波区域，导致格心网格数据的激波特征定位准确性较差，从而使得设计出的飞行器存在安全隐患。

技术方案

本发明要解决的技术问题是：针对流场数值模拟产生的格心网格数据，在两级采样法的基础上，提供一种基于双控制体的三维激波特征定位方法，提高三维激波定位的精度，解决格心网格三维流场数据的激波特征不能精确定位的问题。

本发明技术方案是：为实现基于格心网格数据的激波定位，首先直接基于计算流体力学模拟仿真得到的格心网格数据计算各网格单元内的压力梯度，然后采用一种基于双控制体的采样点正则马赫数计算方法得到采样点的正则马赫数，根据正则马赫数的值得到备选激波点；然后采用“两级采样法”对备选激波点进行过滤，得到真实激波点，这些激波点组成的区域就是激波特征区域。

具体技术方案为：

步骤1：基于计算流体力学模拟仿真得到的格心网格数据确定各网格单元内部面中心处的压力值：采用Roe平均方法，确定所有网格单元内部面(即被两个网格单元共享的面)的压力面通量即面中心处压力值，Roe平均方法见文献(QianliMa,HuaxunXu,LiangZeng,XunCai,SikunLi.Directraycastingofunstructuredcell-centereddatabydiscontinuityRoe-averagecomputation.TheVisualComputer,2010,26(6-8):1049-1059(马千里,徐华勋,曾亮,蔡勋,李思昆.基于不连续Roe平均计算的格心网格数据直接体可视化方法.TheVisualComputer,2010,26(6-8):1049-1059.))；

步骤2：采用高斯-格林公式(张筑生.数学分析新讲(第3册).北京大学出版社,1991年,第97页)，得到各网格单元的压力梯度，具体如下：

设当前要确定压力梯度的网格单元为C₀，它的第k个相邻网格单元记作C_k(k∈{1,2,…,n})，其中n表示C₀的相邻网格单元的总数。分别用p_k和n_k表示C₀与C_k邻接面的压力面通量和面法向，根据高斯-格林公式，网格单元C₀内的压力梯度▽p_C0，(▽p表示网格单元的压力梯度，▽p_x表示网格单元x的压力梯度，以下都用这种表示方法)满足公式1：

公式1其中，V和S分别是网格单元C₀的体积和表面积。若用S_k表示C₀与C_k邻接面的面积，则公式1可变换为公式2：

$V\·\▿p_{C_0}=(\underset{k=1}{\overset{n}{\Σ}}{p}_k\·n_k\·S_k)$ 公式2公式2中，只有C₀的压力梯度是未知数，求解公式2构成的方程即可得到网格单元C₀的压力梯度；

步骤3：采用背景技术所述的两级采样法步骤2获得一级采样点位置(根据两级采样法，采样点位置在各单元的面上)，即：对格点网格数据进行采样，得到一级采样点。

步骤4：判断一级采样点位置，根据其位置分别采用如下方法计算一级采样点处的正则马赫数：

步骤4.1：如果采样点所在面不是内部面(即一级采样点所在面只属于一个网格单元)，利用该网格单元的压力梯度▽p(步骤2得到)和文献(LovelyD.,R.Haimes.Shockdetectionfromcomputationalfluiddynamicsresults.InProceedingsofthe14thAIAAComputationalFluidDynamicsConference,1999.(LovelyD,RHaimes.从计算流体力学结果中检测激波.第14界美国航空航天学会计算流体力学会议,1999.))中的正则马赫数计算公式得到一级采样点处的正则马赫数，其中M_n表示一级采样点的正则马赫数，v表示采样点处的速度，w表示采样点处的局部音速，两者都可由输入数据(即格心网格数据)直接得到，▽p就是网格单元的压力梯度，||▽p||是▽p的2-范数，令最终正则马赫数M＝M_n，转步骤5；

步骤4.2：如果一级采样点所在面是内部面，采用双控制体方法得到一级采样点处的最终正则马赫数M。具体步骤为：

步骤4.2.1：确定一级采样点所在面所属的两个网格单元C_a和C_b(面和单元的所属关系可由输入数据直接得到)即一级采样点的双控制体；

步骤4.2.2：利用网格单元C_a的压力梯度(步骤2得到)和步骤4.1中的正则马赫数计算公式得到对应C_a的一级采样点的正则马赫数M_a，其中M_a表示利用C_a得到的一级采样点的正则马赫数，v表示该一级采样点处的速度，w表示一级采样点处的局部音速，表示C_a的压力梯度，是的2-范数；利用网格单元C_b的压力梯度(步骤2得到)和步骤4.1中的正则马赫数计算公式得到对应C_b的一级采样点的正则马赫数M_b，其中M_b表示利用C_b得到的一级采样点的正则马赫数，v表示该一级采样点处的速度，w表示一级采样点处的局部音速，表示C_b的压力梯度，是的2-范数；

步骤4.2.3：一级采样点处的最终正则马赫数M＝(M_a+M_b)/2；

步骤5：同“两级采样法”步骤3，根据一级采样点处的最终正则马赫数M对一级采样点进行初步筛选：若M大于等于1，则该一级采样点为备选激波点，若M小于1，则转步骤7；

步骤6：同“两级采样法”步骤4，通过二级采样点(二级采样点计算方法如背景技术所述)判断步骤5得到的备选激波点是否为真实激波点；

步骤7：结束。

与现有技术相比，采用本发明可达到以下技术效果：本发明直接基于原始格心数据进行激波特征定位，没有损失数据精度，有效解决了格心网格三维流场数据的激波特征定位问题，从而使得设计的飞行器具有更好的安全性能。

附图说明

图1是本发明总体流程图。

图2是当一级采样点所在面是内部面时，利用双控制体方法确定一级采样点处正则马赫数的流程图。

具体实施方式

图1是本发明的总体流程图。具体步骤说明如下：

步骤1：采用Roe平均方法得到所有内部面中心处的压力值；

步骤2：利用高斯-格林公式得到各网格单元的压力梯度；

步骤3：对格点网格数据进行采样，得到一级采样点位置；

步骤4：判断一级采样点位置，若一级采样点所在面不是内部面，则执行步骤4.1，若其所在面是内部面执行步骤4.2；

步骤4.1：利用正则马赫数计算公式得到一级采样点处的最终正则马赫数M，转步骤5；

步骤4.2：采用双控制体方法得到一级采样点处的最终正则马赫数M；

步骤5：如果一级采样点处的最终正则马赫数满足M>＝1，则该一级采样点为备选激波点，转步骤6，若M<1，转步骤7；

步骤6：根据二级采样点判断该备选激波点是否为真实激波点；

步骤7：结束。

图2是当一级采样点所在面是内部面时，利用双控制体方法计算一级采样点处正则马赫数的流程图。具体步骤说明如下：

步骤1：确定一级采样点所在面所属的两个网格单元C_a和C_b，即一级采样点的双控制体；

步骤2：分别用C_a和C_b计算一级采样点处的正则马赫数M_a和M_b；

步骤3：一级采样点处的最终正则马赫数M＝(M_a+M_b)/2。

Claims (1)

1.一种基于双控制体的格心网格数据三维激波特征定位方法，其特征在于包括以下步骤：

步骤1：基于计算流体力学模拟仿真得到的格心网格数据确定各网格单元内部面中心处的压力值，方法是：采用Roe平均方法，确定所有网格单元内部面的压力面通量即面中心处压力值；

步骤2：采用高斯-格林公式，得到各网格单元的压力梯度，方法是：

设当前要确定压力梯度的网格单元为C₀，它的第k个相邻网格单元记作C_k，k∈{1,2,…,n}，其中n表示C₀的相邻网格单元的总数；分别用p_k和n_k表示C₀与C_k邻接面的压力面通量和面法向，根据高斯-格林公式，网格单元C₀内的压力梯度满足公式1：

公式1

其中，V和S分别是网格单元C₀的体积和表面积，若用S_k表示C₀与C_k邻接面的面积，则公式1变换为公式2：

$V\&CenterDot;\&dtri;p_{C_0}=(\underset{k=1}{\overset{n}{\&Sigma;}}{p}_k\&CenterDot;n_k\&CenterDot;S_k)$ 公式2

公式2中，只有压力梯度是未知数，求解公式2构成的方程即得到网格单元C₀的压力梯度；用表示网格单元的压力梯度，表示网格单元x的压力梯度；

步骤3：对格点网格数据进行采样，得到一级采样点；

步骤4：判断一级采样点位置，根据一级采样点位置分别采用如下方法计算一级采样点处的正则马赫数：

步骤4.1：如果一级采样点所在面不是内部面即一级采样点所在面只属于一个网格单元，利用该网格单元的压力梯度和正则马赫数计算公式得到一级采样点处的正则马赫数M_n，v表示采样点处的速度，w表示采样点处的局部音速，v和w都由格心网格数据直接得到，是的2-范数，令最终正则马赫数M＝M_n，转步骤5；

步骤4.2：如果一级采样点所在面是内部面，采用双控制体方法得到一级采样点处的最终正则马赫数M，具体步骤为：

步骤4.2.1：确定一级采样点所在面所属的两个网格单元C_a和C_b，即一级采样点的双控制体；

步骤4.2.2：利用网格单元C_a的压力梯度和步骤4.1中的正则马赫数计算公式得到对应C_a的一级采样点的正则马赫数M_a， 是的2-范数；利用网格单元C_b的压力梯度和正则马赫数计算公式得到对应C_b的一级采样点的正则马赫数M_b， $M_b=\frac{v}{w}\&CenterDot;\frac{\&dtri;p_{C_b}}{\left|\right|\&dtri;p_{C_b}\left|\right|},$ 是的2-范数；

步骤4.2.3：一级采样点处的最终正则马赫数M＝(M_a+M_b)/2；

步骤5：根据一级采样点处的最终正则马赫数M对一级采样点进行初步筛选：若M大于等于1，则该一级采样点为备选激波点，若M小于1，转步骤7；

步骤6：通过二级采样点判断步骤5得到的备选激波点是否为真实激波点，方法是：

步骤6.1：确定备选激波点处的切平面及切平面法线，该法线及其反向延长线与一系列网格单元相交，产生一系列交点，选取在备选激波点两侧的两个交点作为二级采样点；

步骤6.2：确定两个二级采样点处的速度；

步骤6.3：计算步骤6.2两个二级采样点处速度在切平面上的投影速度的大小；

步骤6.4：如果两个投影速度大小相同，则该备选激波点为真实激波点，否则为伪激波点；

步骤7：结束。