CN112818573B - 一种用于非结构网格的获取边界层非当地变量信息的方法 - Google Patents

Abstract

本发明涉及一种用于非结构网格的获取边界层非当地变量信息的方法，属于流体力学中的边界层转捩预测技术领域。本发明方法，首先对网格进行了预处理，所有获取边界层信息的计算均在w(i,j)数组中完成，由于w(i,j)数组表征壁面网格单元的i彼此独立，因此可进行大规模并行计算，适用于现代的CFD计算方法。采用循环盒子法寻找最近壁面距离，在此基础上完成对网格的重排序，具有较高计算效率。由于w(i,j)数组基本沿着壁面法向方向，且进行了壁面距离修正、边界外缘以及边界层非当地变量参数修正，所获结果具有较高精度。

Description

一种用于非结构网格的获取边界层非当地变量信息的方法

技术领域

本发明涉及一种用于非结构网格的获取边界层非当地变量信息的方法，属于流体力学中的边界层转捩预测技术领域。

背景技术

长期以来，边界层转捩作为流体力学领域的一个前沿和难点问题，受到了广泛而持续的关注。尽管关于边界层转捩的基础理论和发生机理仍在研究发展中，但其在工程实践中的重要性已经日益突出。例如，在航空航天工程领域中，转捩的发生位置和发展过程对飞行器的升阻特性、边界层分离以及表面气动加热等具有显著的影响，因此，边界层转捩的准确预测，对于航空航天飞行器气动外形、控制系统等设计具有非常重要的意义。

边界层转捩表征了流体从简单的分层稳定状态到复杂的混沌湍流状态的过渡，是十分重要却又极其复杂的流动问题。由于边界层转捩会显著影响飞行器以及叶轮机叶片的流动分离位置、壁面摩擦阻力、壁面热流交换、进气道流场品质、发动机燃料混合以及燃烧等，因此，研究有效的边界层转捩预测方法，对飞行器气动性能的准确预测、热防护系统和推进系统的合理设计等，均具有十分重要的意义。

目前，国内外研究人员已经提出了多种边界层转捩预测方法，包括半经验的e^N方法、转捩模式、转捩准则、直接数值模拟和大涡模拟等。其中，转捩模式、转捩准则等方法，由于计算量小、可靠性强，经过风洞或飞行实验数据标定之后具有较高的预测精度，在工程中广泛使用。但是，这两种方法在构造过程中，会使用到一些边界层非当地变量信息，例如：边界层内横流速度w、最大横流速度w_max、边界层动量厚度θ、动量厚度雷诺数Re_θ、边界层外缘速度u_e、边界层外缘马赫数Ma_e以及边界层外缘温度T_e，等等，这些参数在现代的基于大规模并行和非结构网格的计算流体力学(CFD)方法中难以获取，极大限制了转捩模式、转捩准则方法在复杂外形边界层转捩预测中的应用，预测精度也难以提升。

发明内容

本发明的目的是为了克服现有技术的缺陷，为解决在基于大规模并行和非结构网格的计算流体力学方法中，难以并行获取边界层非当地变量信息的技术问题，导致所构造的转捩模式和转捩准则方法受几何外形约束，应用范围存在局限、预测精度难以进一步提升等，提出一种用于非结构网格的获取边界层非当地变量信息的方法。

一种用于非结构网格的获取边界层非当地变量信息的方法，包括如下步骤：

步骤1：进行网格预处理。

在进行转捩迭代计算之前，首先对读入并存储的原始非结构网格数据进行预处理，对其重新排序，并构造适用于大规模并行求解的数据存储方案。

具体如下：

步骤1.1：将非结构网格数据标记为壁面网格单元和流场网格单元。采用循环盒子法，求解每个流场网格单元的壁面距离，并标记与其距离最近的壁面网格单元。

具体地，采用循环盒子法求解的方法如下：

步骤1.1.1：统一编号所有壁面网格单元，获取所有壁面网格点三个方向(x,y,z)的最大最小坐标，形成一个包含所有壁面网格单元的外接长方形盒子。

步骤1.1.2：采用二分方式，沿外接方盒子最长边将其等分，将相应的壁面网格单元分别存储于两个方盒子中。重复上述划分过程，当划分后的方盒子内含有壁面网格单元数小于50或者划分最大层次达到10层后，停止划分，形成N个方盒子，此时，每个方盒子包含有若干数量的壁面网格单元。

步骤1.1.3：对于任意流场网格单元，分别计算该点到N个方盒子表面的最近壁面距离，依照最近距离，采用快速排序法对N个方盒子进行排序，选出距离最近的方盒子。

此时，要对循环盒子法求解的壁面距离进行修正。

循环盒子法以流场网格单元格心点坐标到最近壁面网格单元格心点坐标的直线距离作为最小壁面距离。假设B₁点为任一空间流场网格单元的格心，A点为与其对应的最近壁面网格单元的格心，则A、B₁点间的距离d_AB1为原始的循环盒子法求得的壁面距离。

但是，在计算网格生成时，尤其是在复杂网格的情况下，很难保证AB₁线与对应壁面正交，因为AB₁线可能与壁面法向量间存在一个夹角α，求解的壁面距离d_AB1则会与真实的壁面距离d_AB2存在一定的偏差，因此，需要对壁面距离进行修正。修正形式为：d_AB2＝d_AB1×cosα。

步骤1.1.4：采用直接比较法，计算出步骤1.1.3中选取的任意网格单元到距离最近方盒子内所包含的壁面网格单元的距离，其中，最近的距离即为该点到最近方盒子的距离，同时记录所对应的壁面网格单元。

步骤1.1.5：对步骤1.1.3中所获得的距离第二近的方盒子，重复步骤1.1.4操作，如果找到相比于步骤1.1.4更加近的距离，即为点到壁面的最近距离，同时更新对应的壁面网格单元编号。

记录壁面网格单元，是为了对流场网格重新排序，从而建立以壁面网格为索引的全新数据结构，便于后续边界层参数的辨识。

步骤1.2：建立一个二维指针数组w(i,j)，其中，第一维i用于记录每个壁面网格单元序号，第二维j记录所有将该壁面网格单元标记为最近壁面网格单元的流场网格序号，并按壁面距离由大至小进行排序。

步骤1.3：二维指针数组w(i,j)中存储的是i个壁面单元对应的各自j个法向流场网格单元。计算每一个壁面网格对应的边界层外缘，求取边界层内横流速度w、最大横流速度w_max、边界层动量厚度θ、边界层动量厚度雷诺数Re_θ等参数。

其中，将边界层内所有网格单元的速度沿垂直于边界层外缘的无粘流速度的流向方向进行投影，获得边界层内所有网格单元的横流速度w，对该横流速度从1到j_edge进行比较判断，获得最大横流速度w_max。

边界层动量厚度θ计算如下：

其中，u表示边界层内的当地速度，u_e表示边界层外缘速度，dy表示沿壁面法向积分；

动量厚度雷诺数Re_θ计算如下：

其中，ρ_e表示边界层外缘密度，μ_e表示边界层外缘流体动力粘性系数。

步骤2：判断边界层外缘。

沿网格重排序获得的边界层查寻路径，即利用数组w(i,j)进行搜寻，当满足给定的边界层外缘判别准则时，即将其确定为边界层外缘。

具体地，边界层外缘判别的判据可以包括：

(1)速度判据。

可以设置为：当地流向速度达到自由来流速度的99％，即u/u_∞＝0.99；

(2)速度导数判据。

可以设置为：当地无量纲流向速度沿法向的梯度减小到0.01，即

其中，L表示参考长度，d表示壁面距离。

(3)质量通量判据。

可以设置为：当地质量通量达到来流质量通量的99％，即(ρu)/(ρu)_∞＝0.99。其中，ρ表示密度，u表示速度。

(4)总焓判据。

可以设置为：当地总焓h₀达到自由来流总焓h_0,∞的1.0005倍，即h₀/h_0,∞＝1.0005。

对于边界层内速度剖面存在拐点且边界层外缘存在激波时，可以选用总焓判据。

注意，上述4种判据的临界值，仅是本发明给出的参考值。本领域技术人员在具体使用时，可以根据实际情况进行调整。

步骤3：在已获得边界层外缘信息基础之上，获取其他相关的非当地变量参数信息。

具体如下：

基于步骤2选定的边界层外缘判据，在w(i,j)数组中对每一个壁面网格单元i，逐个比较判断对应的j个壁面法向网格单元，将满足判据的网格单元标记为w(i,j_edge)，即为边界层外缘网格单元。指针w(i,j_edge)指向的网格单元存储的壁面距离信息即为边界层厚度δ，该网格单元存储的温度、马赫数和速度信息，即为边界层外缘温度T_e、边界层外缘马赫数Ma_e和边界层外缘速度u_e。

步骤4：对边界层外缘信息和非当地变量参数信息进行修正。

其中，对边界层外缘信息的修正，即对边界层厚度δ的修正，方法如下：

以无量纲总焓判断边界层外缘为例：A点位于边界层内，对应的最小壁面距离为d_A，总焓比为h_0A/h_0∞，B点位于边界层外，对应的最小壁面距离为d_B，总焓比为h_0B/h_0∞，精确的边界层外缘位于A、B点之间。对于非结构网格流场数据存储在计算网格的格心处，当采用边界层判据进行边界层外缘辨识时，A点不满足判据，B点满足判据，则认定B点处的壁面距离及对应的流场变量作为边界层外缘参数。当边界层附近网格相对稀疏时，边界层厚度分布会出现抖动，且比准确的边界层厚度值偏大。为了更加精确的辨识边界层参数，通过线性插值的方法对此进行修正，具体表达式如下：

δ＝d_A+(d_B-d_A)×k (3)

k＝(h_{0_crit}/h_0∞-h_0A/h_0∞)/(h_0B/h_0∞-h_0A/h_0∞) (4)

其中，h_{0_crit}/h_0∞为总焓判据中总焓比的临界值，k为定义的线性插值系数。

对其他非当地变量参数信息的修正，同样可以使用类似的插值处理方法，考虑到不同参数沿壁面法向变化趋势不同，具体构造略有差异。具体如下所示：

Ma_e＝Ma_A+(Ma_B-Ma_A)×k (5)

其中，Ma_A表示A点的边界层外缘马赫数，Ma_B表示B点的边界层外缘马赫数。

T_e＝T_A-(T_A-T_B)×k (6)

其中，T_A表示A点的边界层外缘温度，T_B表示B点的边界层外缘温度。

ρ_e＝ρ_A+(ρ_B-ρ_A)×k (7)

其中，ρ_A表示A点的边界层外缘密度，ρ_B表示B点的边界层外缘密度。

θ＝θ_A+(θ_B-θ_A)×k (8)

其中，θ_A表示A点的边界层动量厚度，θ_B表示B点的边界层动量厚度。

对横流速度w的修正方法如下：

横流速度w是边界层内出现的垂直于边界层外无粘流动的速度分量，其求解的关键在于无粘流速度矢量的定义。由于横流速度w在判别的边界层外缘之外仍存在一定的梯度变化，因此，直接采用边界层判据获取的边界层外缘速度作为无粘流速度并不适用于横流速度的求解，需要对无粘流速度进行重新定义。

考虑到横流速度w的分布特点，在边界层判据判定的边界层外缘的基础上，沿壁面法向继续搜寻，提出一种斜率控制的方法来定义无粘流速度。以总焓比判据为例，定义偏折角β_i如下：

其中，脚标i代表网格编号，i+1为数组w(i,j)中沿法向的下一个网格，d和h₀分别代表壁面距离和总焓，d_i、d_i+1分别表示网格编号为i、i+1网格单元的壁面距离。当β_i、β_i+1和β_i+2同时大于80°时，取网格编号为i+1的速度为无粘流速度。

步骤5：将修正后的边界层外缘信息和非当地变量参数信息，用于非结构网格的外形边界层转捩预测。

有益效果

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

1.本发明对网格进行了预处理，所有获取边界层信息的计算均在w(i,j)数组中完成，由于w(i,j)数组表征壁面网格单元的i彼此独立，因此可进行大规模并行计算，适用于现代的CFD计算方法。

2.本发明采用循环盒子法寻找最近壁面距离，在此基础上完成对网格的重排序，具有较高的计算效率。

3.本发明由于w(i,j)数组基本沿着壁面法向方向，且进行了壁面距离修正、边界外缘以及边界层非当地变量参数修正，因此结果具有较高的精度。

附图说明

图1为本发明方法的流程图；

图2为本发明方法的网格预处理示意图；

图3为本发明基于非结构网格对边界层外缘进行修正的示意图；

图4为本发明实施例1中，针对T3A平板，采用本发明方法得到的边界层厚度分布与理论值的对比。

图5为本发明实施例1中，针对T3A平板，采用本发明方法得到的边界层动量厚度分布与理论值的对比。

图6为本发明实施例2中，针对后掠圆柱，采用本发明方法得到的横流速度分布与文献结果的对比。

具体实施方式

下面结合附图及实施例对本发明做进一步详细说明。

实施例1

本实例描述了应用本发明方法在针对T3A平板时的具体实施方案。

图1为本发明方法的流程图。图2为本发明方法网格预处理示意图。图3为本发明基于非结构网格对边界层外缘进行修正的示意图。据此描述本发明在特定应用场景下的具体实施。

T3A平板计算域为长L＝1000mm、高H＝150mm的矩形，入口距离平板前缘80mm。计算采用混合网格，边界层内为结构网格，边界层外为非结构网格，壁面和平板前缘附近网格加密。边界层内沿流向为111个网格单元，沿法向为60个网格单元，壁面第一层网格高度保证y⁺＝1。壁面为无滑移固壁边条，入口和出口分别为速度入口和压力出口边界，上下远场设为对称面。自由来流雷诺数为Re_∞＝6.12×10⁵/m，来流速度为u_∞＝5.4m/s，壁面为绝热壁，采用全层流计算。

采用u/u_∞＝0.99和h₀/h_0,∞＝1.0005判据计算得到T3A平板在计算域0～1000mm范围内的边界层厚度δ和动量厚度θ分布，并与Blasius精确解和Karman近似解进行对比。Blasius精确解和Karman近似解分别为

θ＝(39/280)δ。图4和图5为本发明采用u/u_∞＝0.99和h₀/h_0,∞＝1.0005判据获得的边界层厚度δ分布和积分求得的边界层动量厚度θ分布与Blasius精确解和Karman近似解的对比，在整个计算域内，本发明结果均与理论结果十分接近，整体误差分别为Δδ＝(δ_Num-δ_theory)/L＜0.8％和Δθ＝(θ_Num-θ_theory)/L＜0.12％。

本发明方法，能准确判断边界层外缘并获得相关边界层非当地变量信息。

实施例2

本实例描述了应用本发明方法在针对后掠圆柱计算模型时的具体实施方案。

圆柱半径为30mm，后掠角为50°。自由来流和壁面条件为：Ma_∞＝3，Re_∞＝2.83×10⁷/m，T_∞＝121.42K，T_w＝300K。壁面第一层网格高度保证y⁺＝1，激波位置附近网格加密。壁面设置为无滑移固壁边条，两侧为周期性边条。

根据本发明给出的网格预处理、获取边界层信息的方法，得到了后掠圆柱边界层内的横流速度分布。图6为方位角

范围内得到的无量纲横流速度w/u_e与Balakumar等人文献结果的对比，可以看到，不同

站位的横流速度剖面基本一致，壁面和边界层外缘的值为0，在距离壁面0.04mm位置处横流速度达到最大值。往下游方向，随着

值增大，最大横流速度也逐渐增长，到

站位处，边界层内最大横流速度几乎达到边界层外缘速度的8％。

使用本发明的方法得到的横流速度剖面分布与Balakumar等人给出的结果几乎完全吻合，说明本发明方法能准确获取边界层内的横流速度信息。

Claims (7)

1.一种用于非结构网格的获取边界层非当地变量信息的方法，其特征在于，包括如下步骤：

步骤1：对原始非结构网格数据进行预处理；

步骤1.1：将非结构网格数据标记为壁面网格单元和流场网格单元，采用循环盒子法，求解每个流场网格单元的壁面距离，并标记与其距离最近的壁面网格单元；

采用循环盒子法求解的方法如下：

步骤1.1.1：统一编号所有壁面网格单元，获取所有壁面网格点三个方向(x,y,z)的最大最小坐标，形成一个包含所有壁面网格单元的外接长方形盒子；

步骤1.1.2：采用二分方式，沿外接方盒子最长边将其等分，将相应的壁面网格单元分别存储于两个方盒子中；重复上述划分过程，当划分后的方盒子内含有壁面网格单元数小于50或者划分最大层次达到10层后，停止划分，形成N个方盒子，此时，每个方盒子包含有若干数量的壁面网格单元；

步骤1.1.3：对于任意流场网格单元，分别计算该网格单元到N个方盒子表面的最近壁面距离，依照最近距离，采用快速排序法对N个方盒子进行排序，选出距离最近的方盒子；

此时，对循环盒子法求解的壁面距离进行修正；

循环盒子法以流场网格单元格心点坐标到最近壁面网格单元格心点坐标的直线距离作为最小壁面距离，假设B₁点为任一空间流场网格单元的格心，A点为与其对应的最近壁面网格单元的格心，则A、B₁点间的距离d_AB1为原始的循环盒子法求得的壁面距离，对壁面距离进行修正：d_AB2＝d_AB1×cosα，其中，α为AB₁线与壁面法向量间的夹角，d_AB2为真实的壁面距离；

步骤1.1.4：采用直接比较法，计算出步骤1.1.3中选取的任意网格单元到距离最近方盒子内所包含的壁面网格单元的距离，其中，最近的距离即为该网格单元到最近方盒子的距离，同时记录所对应的壁面网格单元；

步骤1.1.5：对步骤1.1.3中所获得的距离第二近的方盒子，重复步骤1.1.4操作，如果找到相比于步骤1.1.4更加近的距离，即为网格单元到壁面的最近距离，同时更新对应的壁面网格单元编号；

步骤1.2：建立一个二维指针数组w(i,j)，其中，第一维i用于记录每个壁面网格单元序号，第二维j记录所有将该壁面网格单元标记为最近壁面网格单元的流场网格序号，并按壁面距离由大至小进行排序；

步骤1.3：二维指针数组w(i,j)中存储的是i个壁面单元对应的各自j个法向流场网格单元；计算每一个壁面网格对应的边界层外缘，求取边界层内横流速度w、最大横流速度w_max、边界层动量厚度θ、边界层动量厚度雷诺数Re_θ；

其中，将边界层内所有网格单元的速度沿垂直于边界层外缘的无粘流速度的流向方向进行投影，获得边界层内所有网格单元的横流速度w，对该横流速度从1到j_edge进行比较判断，其中下标edge表示边界层外缘网格单元，获得最大横流速度w_max；

边界层动量厚度θ计算如下：

其中，δ表示边界层厚度，u表示边界层内的当地速度，u_e表示边界层外缘速度，dy表示沿壁面法向积分；

动量厚度雷诺数Re_θ计算如下：

其中，ρ_e表示边界层外缘密度，μ_e表示边界层外缘流体动力粘性系数；

步骤2：判断边界层外缘；

沿网格重排序获得的边界层查寻路径，即，利用数组w(i,j)进行搜寻，当满足给定的边界层外缘判别准则时，即将其确定为边界层外缘；

其中，判别准则，内容包括速度判据、速度导数判据、质量通量判据和总焓判据；

步骤3：在已获得边界层外缘信息基础之上，获取其他相关的非当地变量参数信息；

基于步骤2选定的边界层外缘判据，在w(i,j)数组中对每一个壁面网格单元i，逐个比较判断对应的j个壁面法向网格单元，将满足判据的网格单元标记为w(i,j_edge)，即为边界层外缘网格单元；指针w(i,j_edge)指向的网格单元存储的壁面距离信息即为边界层厚度δ，该网格单元存储的温度、马赫数和速度信息，即为边界层外缘温度T_e、边界层外缘马赫数Ma_e和边界层外缘速度u_e；

步骤4：对边界层外缘信息和非当地变量参数信息进行修正；

其中，对边界层外缘信息的修正，即对边界层厚度δ的修正，方法如下：

基于无量纲总焓判断边界层外缘，A点位于边界层内，对应的最小壁面距离为d_A，总焓比为h_0A/h_0∞，B点位于边界层外，对应的最小壁面距离为d_B，总焓比为h_0B/h_0∞，精确的边界层外缘位于A、B点之间；对于非结构网格流场数据存储在计算网格的格心处，当采用边界层判据进行边界层外缘辨识时，A点不满足判据，B点满足判据，则认定B点处的壁面距离及对应的流场变量作为边界层外缘参数；通过线性插值的方法对边界层厚度进行修正，具体表达式如下：

δ＝d_A+(d_B-d_A)×k (3)

k＝(h_{0_crit}/h_0∞-h_0A/h_0∞)/(h_0B/h_0∞-h_0A/h_0∞) (4)

其中，h_{0_crit}/h_0∞为总焓判据中总焓比的临界值，k为定义的线性插值系数；

对边界层外缘马赫数Ma_e的修正如下所示：

Ma_e＝Ma_A+(Ma_B-Ma_A)×k (5)

其中，Ma_A表示A点的边界层外缘马赫数，Ma_B表示B点的边界层外缘马赫数；

对边界层外缘温度T_e的修正如下所示：

T_e＝T_A-(T_A-T_B)×k (6)

其中，T_A表示A点的边界层外缘温度，T_B表示B点的边界层外缘温度；

对边界层外缘密度ρ_e表示边界层外缘密度的修正如下所示：

ρ_e＝ρ_A+(ρ_B-ρ_A)×k (7)

其中，ρ_A表示A点的边界层外缘密度，ρ_B表示B点的边界层外缘密度；

对边界层动量厚度θ的修正如下所示：

θ＝θ_A+(θ_B-θ_A)×k (8)

其中，θ_A表示A点的边界层动量厚度，θ_B表示B点的边界层动量厚度；

步骤5：将修正后的边界层外缘信息和非当地变量参数信息，用于非结构网格的外形边界层转捩预测。

2.如权利要求1所述的一种用于非结构网格的获取边界层非当地变量信息的方法，其特征在于，步骤2中，速度判据设置为：当地流向速度达到自由来流速度的99％，即u/u_∞＝0.99，其中u_∞表示来流速度。

3.如权利要求1所述的一种用于非结构网格的获取边界层非当地变量信息的方法，其特征在于，步骤2中，速度导数判据设置为：当地无量纲流向速度沿法向的梯度减小到0.01，即，

其中，L表示参考长度，d表示壁面距离。

4.如权利要求1所述的一种用于非结构网格的获取边界层非当地变量信息的方法，其特征在于，步骤2中，质量通量判据设置为：当地质量通量达到来流质量通量的99％，即(ρu)/(ρ_∞u_∞)＝0.99，其中，ρ表示密度，u表示速度，下标∞表示来流变量。

5.如权利要求1所述的一种用于非结构网格的获取边界层非当地变量信息的方法，其特征在于，步骤2中，总焓判据设置为：当地总焓达到来流总焓的1.0005倍，即h₀/h_0,∞＝1.0005，其中，h₀表示当地总焓，h_0,∞表示来流总焓。

6.如权利要求1所述的一种用于非结构网格的获取边界层非当地变量信息的方法，其特征在于，步骤2中，当边界层内速度剖面存在拐点且边界层外缘存在激波时，选用总焓判据。

7.如权利要求1所述的一种用于非结构网格的获取边界层非当地变量信息的方法，其特征在于，步骤3中，对横流速度w进行修正，方法如下：

在边界层判据判定的边界层外缘的基础上，沿壁面法向继续搜寻，提出一种斜率控制的方法来定义无粘流速度，根据总焓比判据，定义偏折角β_i如下：

其中，脚标i代表网格编号，i+1为数组w(i,j)中沿法向的下一个网格，d和h₀分别代表壁面距离和总焓，d_i、d_i+1分别表示网格编号为i、i+1网格单元的壁面距离；当β_i、β_i+1和β_i+2同时大于80°时，取网格编号为i+1的速度为无粘流速度。

Images