[发明专利]一种用于非结构网格的获取边界层非当地变量信息的方法

权利要求书

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：将修正后的边界层外缘信息和非当地变量参数信息，用于非结构网格的外形边界层转捩预测。