一种飞行器湍流绕流模拟的壁面距有效单元计算方法与流程

1.本发明属于计算流体力学技术领域，尤其涉及一种飞行器湍流绕流模拟的壁面距有效单元计算方法。


背景技术：

2.湍流是飞行器绕流中最为常见的流动。其运动随机不规则，具有不稳定、非线性、涡结构丰富、强耗散等特征。受计算机硬件技术的限制，利用计算流体力学数值模拟飞行器湍流绕流流场在目前以及未来数十年的主要途径均为雷诺平均方法（rans）或脱体涡模拟方法（des）。而在这两类湍流模拟方法中均存在一个对模拟精度至关重要的参数——壁面距，即计算网格单元格心到飞行器表面的最小距离。
3.当前，壁面距的确定方法主要包括搜索法和求解微分方程法两大类。求解微分方程法需迭代求解eikonal方程或hamilton
‑
jacobi方程，算法实现复杂；对于复杂外形，其求解效率与精度更是难以保证。搜索法中最为简单的是穷举搜索法，即计算每一流场单元对所有壁面单元的距离，取其中的最小值作为该单元到壁面的距离。不过，穷举搜索法效率低下，对于大规模动网格数值模拟这类需反复计算壁面距又无法提前存储且要求计算效率的情况难以适用。
4.对于提高搜索法的效率，目前主要通过引入高效搜索数据结构和壁面面元分组来减少搜索量实现。例如，现有文献提出一种名为“循环盒子法”的壁面距离计算方法，利用八叉树数据结构对壁面点按坐标排序，并通过循环二等分壁面点的外接长方体实现了对壁面单元的分组。然而，现有改进方法主要是从几何学或计算机学的角度出发。若从计算流体力学的角度出发，壁面距只有在近壁湍流区至关重要，而在其他区域该变量则完全不影响计算。显然，若在大规模动网格数值模拟中仍在每一物理时刻采用计算所有网格单元壁面距的方式，即便采用最高效的搜索数据结构与点
‑
面距离计算方法也仍会存在大量无效计算，从而严重影响数值模拟的计算效率。


技术实现要素：

5.为了解决飞行器湍流绕流模拟中关键参数壁面距计算效率不足的问题，本发明提出一种飞行器湍流绕流模拟的壁面距有效单元计算方法，本发明的具体技术方案如下：一种飞行器湍流绕流模拟的壁面距有效单元计算方法，包括以下步骤：s1：读入并存储计算网格单元格点坐标与边界条件；s2：壁面单元预处理；s2
‑
1：收集壁面单元；s2
‑
2：确定壁面单元外接长方体的个数：根据计算网格壁面单元数，确定循环盒子法中使用的壁面单元外接长方体个数，壁面单元外接长方体是指包含其内壁面单元的最小长方体；s2
‑
3：根据所有壁面单元建立初始外接长方体，外接长方体所需记录的信息包括：
外接长方体的笛卡尔坐标范围、所含壁面单元数量、所含壁面单元编号；s2
‑
4：二分壁面点最多的外接长方体；将包含最多壁面单元的外接长方体沿其最长边二等分为两个外接长方体，并更新这两个外接长方体所需记录的信息，重复步骤s2
‑
4，直至外接长方体的总数达到步骤s2
‑
3确定的个数；s3：计算壁面距；s4：输出结果。
6.进一步地，所述步骤s2
‑
1包括以下步骤：s2
‑1‑
1：根据步骤s1读入的边界条件，确定计算网格单元中的壁面单元；s2
‑1‑
2：为壁面单元和壁面点编号，确保壁面单元编号与其某一顶点的编号相一致以便搜索；s2
‑1‑
3：分别按壁面点编号存储壁面点坐标、按壁面单元编号存储壁面单元与其顶点的对应关系、按壁面单元编号存储其相邻单元的编号；s2
‑1‑
4：遍历壁面单元，计算并存储各壁面单元的单位法向。
7.进一步地，所述步骤s3包括以下步骤：s3
‑
1：壁面距归零：分配所有计算网格单元存储壁面距的空间，并将其初值赋为0；遍历所有计算网格单元，对任一单元执行步骤s3
‑
2；s3
‑
2：确定壁面距计算方式；若计算网格单元紧邻单元的壁面距均小于设定阈值，则精确计算壁面距，继续步骤s3
‑
3；若计算网格单元紧邻单元的壁面距大于设定的阈值，当计算网格单元位于其紧邻单元的壁面距减小方向时，则精确计算壁面距，继续步骤s3
‑
3；当计算网格单元位于其紧邻单元的壁面距增大方向时，则粗略估计壁面距，跳转至步骤s3
‑
5；s3
‑
3：搜索距离最近的2个外接长方体；令x
i,c
为计算网格单元的格心坐标，x
i,min
、x
i,max
分别表示外接长方体在x
i
方向的最小值和最大值，则计算网格单元到外接长方体的距离d0为：（1）遍历所有外接长方体，根据式（1）筛选出距离最小的2个外接长方体，并对这2个外接长方体执行步骤s3
‑
4；s3
‑
4：精确计算壁面距；s3
‑
5：根据梯度估计壁面距：利用单侧差分计算已赋值紧邻单元处的壁面距梯度，根据该梯度、当前网格单元及其紧邻单元的坐标信息，估计当前网格单元的壁面距。
8.进一步地，所述步骤s3
‑
4包括以下步骤：依据计算网格单元与外接长方体的位置信息，分为两种情况：当计算网格单元所有格点均位于外接长方体内时，执行步骤s3
‑4‑
1；当计算网格单元存在不位于外接长方体内的格点时，执行步骤s3
‑4‑
2；
s3
‑4‑
1：遍历外接长方体内所有壁面单元，逐一计算当前网格单元格点到壁面单元的距离，并存储其中的最小值；s3
‑4‑
2：基于矢量法跳跃式访问外接长方体内壁面单元，在外接长方体内任选一壁面单元，计算并存储当前网格单元到其的距离；计算壁面单元面心指向当前网格单元格心的向量在该壁面单元上的投影方向，该投影方向所指紧邻壁面单元即为下一个计算的壁面单元；跳跃式访问当前外接长方体内壁面单元，直至投影方向指向当前外接长方体的边界。
9.本发明的有益效果在于：1.本发明的飞行器湍流绕流模拟的壁面距有效单元计算方法可显著提升飞行器大规模动网格湍流模拟中壁面距的计算效率，从而实现飞行器非定常气动特性预测的加速。
10.2.相比于现有基于求解微分方程的壁面距计算方法，本发明对于复杂外形具有更好的适用性，结果均基于几何解析关系计算，不引入数值误差。
11.3.相比于现有搜索法中的循环盒子法，本发明采用了更为高效合理的编码方式存储壁面单元，并从待计算网格单元和待搜索壁面单元两方面定义了有效单元，从而显著降低了总计算量。
12.4.相比于现有技术可在确保数值模拟同等精度的同时，显著降低壁面距的总计算量，从而有助于复杂飞行器非定常气动特性预测效率的提升；在步骤s2
‑
1中，本发明以壁面单元与其某一顶点共享编号为特征，为壁面单元与壁面点进行编号，确保了依据壁面点访问壁面单元的高效性；在步骤s3
‑
1中，本发明依据壁面距对湍流模型方程精度的影响，甄别了计算网格单元需要精确计算壁面距的有效单元：对于有效单元，将利用几何关系精确求解壁面距；而对于非有效单元，则利用梯度估计壁面距；在步骤s3
‑
3中，本发明对于计算网格单元不位于外接长方体中的情况建立了基于矢量法的跳跃式壁面单元访问方式，从而进一步降低了精确计算壁面距时访问壁面单元的数量。
附图说明
13.为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，可以根据这些附图获得其他的附图。其中：图1为本发明壁面距有效单元计算方法的流程图；图2为本发明计算rae2822翼型壁面距时生成外接长方体的示意图；图3为本发明与穷举搜索法计算m6机翼壁面距的结果对比。
具体实施方式
14.为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本发明的实施
例及实施例中的特征可以相互组合。
15.在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
16.如图1所示，本发明提出一种飞行器湍流绕流模拟的壁面距有效单元计算方法，本发明的有效单元存在两层含义：其一，待搜索壁面单元中的有效单元，即针对某一计算网格单元，可确定或辅助确定其壁面距的壁面单元；其二，待计算网格单元中的有效单元，即在数值模拟中壁面距对计算精度存在不可忽略影响的单元。对于第一类有效单元，先采用循环盒子法粗筛选壁面单元，再建立矢量法依据壁面单元紧邻关系跳跃式访问粗筛选后的壁面单元，从而降低访问壁面单元的次数。对于第二类有效单元，通过分析湍流模型方程中壁面距的影响范围，实现了影响域内搜索法精确计算、影响域外以梯度估算的壁面距计算策略。
17.具体地，一种飞行器湍流绕流模拟的壁面距有效单元计算方法，包括以下步骤：s1：读入并存储计算网格单元格点坐标与边界条件；s2：壁面单元预处理；s2
‑
1：收集壁面单元；s2
‑1‑
1：根据步骤s1读入的边界条件，确定计算网格单元中的壁面单元；s2
‑1‑
2：为壁面单元和壁面点编号，为便于搜索，编号时令壁面单元与其一尚未编号的顶点使用同一编号，具体选取顶点的方式不做限制；s2
‑1‑
3：分别按壁面点编号存储壁面点坐标、按壁面单元编号存储壁面单元与其顶点的对应关系、按壁面单元编号存储其相邻单元的编号；s2
‑1‑
4：遍历壁面单元，计算并存储各壁面单元的单位法向；对于二维情况，设(x1,y1)、(x2,y2)是壁面单元的2个顶点，则壁面单元的单位法向为：（1）对于三维情况，设(x1,y1,z1)、(x2,y2,z2)、(x3,y3,z3)、(x4,y4,z4)是壁面单元依序排列的4个顶点，则壁面单元的单位法向为：q1=(x3‑
x1y3‑
y1z3‑
z1)，q2=(x4‑
x2ꢀꢀ
y4‑
y2ꢀꢀ
z4‑
z2)（2）s2
‑
2：确定壁面单元外接长方体的个数n：根据计算网格壁面单元数，确定循环盒子法中使用的壁面单元外接长方体个数，壁面单元外接长方体是指包含其内壁面单元的最小长方体；在一些实施方式中，以30个壁面单元分组，计算效率最佳。因此，外接长方体的个数可取n=[壁面单元总数/30]，其中[
•
]代表四舍五入取整。
[0018]
s2
‑
3：根据所有壁面单元建立初始外接长方体，外接长方体所需记录的信息包括：外接长方体的笛卡尔坐标范围、所含壁面单元数量、所含壁面单元编号。
[0019]
s2
‑
4：二分壁面点最多的外接长方体；将包含最多壁面单元的外接长方体沿其最长边二等分为两个外接长方体，并更新这两个外接长方体所需记录的信息，重复步骤s2
‑
4，直至外接长方体的总数达到步骤s2
‑
3确定的个数。具体地：首先，遍历已有的外接长方体，根据其所含壁面单元数量确定包含最多壁面单元的外接长方体；其次，根据该外接长方体的笛卡尔坐标范围，确定其最长边；最长边的方向即为分割外接长方体的方向，该方向坐标范围的中点即为分割坐标；最后，根据分割方向及其坐标，将原外接长方体中的壁面单元分组，并更新新生成的两个外接长方体所需存储的信息。
[0020]
重复步骤s2
‑
4，直至外接长方体的总数达到步骤s2
‑
2确定的个数n。
[0021]
s3：计算壁面距；s3
‑
1：壁面距归零：分配所有计算网格单元存储壁面距的空间，并将其初值赋为0；遍历所有计算网格单元，对任一单元执行步骤s3
‑
2；s3
‑
2：确定壁面距计算方式；若计算网格单元紧邻单元的壁面距均小于设定阈值，则精确计算壁面距，继续步骤s3
‑
3；若计算网格单元紧邻单元的壁面距大于设定的阈值，当计算网格单元位于其紧邻单元的壁面距减小方向时，则精确计算壁面距，继续步骤s3
‑
3；当计算网格单元位于其紧邻单元的壁面距增大方向时，则粗略估计壁面距，跳转至步骤s3
‑
5；例如：分析无量纲化的s
‑
a湍流模型中包含壁面距项的量级表明，当无量纲壁面距比无量纲涡粘性系数大两个量级时，壁面距对湍流模型方程精度的影响可以忽略。由于涡粘性系数在远场处取0.1，则壁面距阈值可取10。
[0022]
若计算网格单元紧邻单元的壁面距均小于10，则精确计算壁面距，继续步骤s3
‑
3。若计算网格单元紧邻单元的壁面距大于10，当计算网格单元位于其紧邻单元的壁面距减小方向时，则精确计算壁面距，继续步骤s3
‑
3；当计算网格单元位于其紧邻单元的壁面距增大方向时，则粗略估计壁面距，跳转至步骤s3
‑
5。
[0023]
s3
‑
3：搜索距离最近的2个外接长方体；令x
i,c
为计算网格单元的格心坐标，x
i,min
、x
i,max
分别表示外接长方体在x
i
方向的最小值和最大值，则计算网格单元到外接长方体的距离d0为：（3）遍历所有外接长方体，根据式（1）筛选出距离最小的2个外接长方体，并对这2个外接长方体执行步骤s3
‑
4；s3
‑
4：精确计算壁面距；
依据计算网格单元与外接长方体的位置信息，分为两种情况：当计算网格单元所有格点均位于外接长方体内时，执行步骤s3
‑4‑
1；当计算网格单元存在不位于外接长方体内的格点时，执行步骤s3
‑4‑
2。
[0024]
s3
‑4‑
1：遍历外接长方体内所有壁面单元，逐一计算当前网格单元格点到壁面单元的距离，并存储其中的最小值；对于二维情况：设壁面单位法向为(n1, n2)，(x
c
, y
c
)为壁面单元的面心，壁面上一点坐标为(x1, y1)，则壁面所在直线为：n1(x
‑
x1) n2(y
‑
y1)=0
ꢀꢀꢀꢀꢀꢀ
（4）经过壁面单元面心，垂直于壁面的直线为：
‑
n2(x
‑
x
c
) n1(y
‑
y
c
)=0（5）则计算单元格心在壁面单元所在平面上的投影点为：x
0=
n
12
x
1 n
22
x
c
n1n
2 (y1‑
y
c
)y
0=
n
12
y
c n
22
y1 n1n
2 (x1‑
x
c
)
ꢀꢀꢀ
（6）根据投影点(x0, y0)指向壁面单元两个壁面点的向量可判断投影点是否在壁面单元内；若在壁面单元内，则两向量内积为负，此时壁面距为计算单元格心到壁面单元的距离；若不在，则两向量内积为正，此时壁面距为计算单元格心到壁面单元顶点的距离。
[0025]
计算网格单元格心到壁面单元的距离为d
w
=|n1(x
c
‑
x1) n2(y
c
‑
y1)|（7）计算网格单元格心到壁面点的距离为（8）对于三维情况：设待求单元的格心坐标为(x
c
,y
c
,z
c
)，壁面单元的单位法向为(n1,n2,n3)，壁面上一点的坐标为(x1,y1,z1)，则壁面单元所在平面为：n1x n2y n3z=dd=n1x1 n2y1 n3z1（9）设格心坐标在壁面单元所在平面上的投影为(x0, y0, z0)，则投影点满足平面方程：n1x0 n2y0 n3z0= d（10）投影点指向计算单元格心的方向即为壁面单元的法向，则（11）则投影点坐标为：x
0=
x
c
n1[d
–
(n1x
c
n2y
c
n3z
c
)]（12）y
0=
y
c
n2[d
–
(n1x
c
n2y
c
n3z
c
)]z
0=
z
c
n3[d
‑
(n1x
c
n2y
c
n3z
c
)]若所有壁面单元顶点指向投影点与下一个壁面点的向量的叉乘均同号，则投影点在壁面单元内，此时壁面距为计算单元格心到壁面单元的距离；若不在，则壁面距为计算单元格心到壁面单元顶点的距离。
[0026]
s3
‑4‑
2：基于矢量法跳跃式访问外接长方体内壁面单元，在外接长方体内任选一壁面单元，计算并存储当前网格单元到其的距离；计算壁面单元面心指向当前网格单元格心的向量在该壁面单元上的投影方向，该投影方向所指紧邻壁面单元即为下一个计算的壁面单元；跳跃式访问当前外接长方体内壁面单元，直至投影方向指向当前外接长方体的边界。
[0027]
s3
‑
5：根据梯度估计壁面距：利用单侧差分计算已赋值紧邻单元处的壁面距梯度，根据该梯度、当前网格单元及其紧邻单元的坐标信息，估计当前网格单元的壁面距。
[0028]
s4：输出结果。
[0029]
为了方便理解本发明的上述技术方案，以下通过具体实施例对本发明的上述技术方案进行详细说明。
[0030]
实施例1本实施例计算rae2822翼型，其中，图2示意了计算rae2822翼型问题壁面距时外接长方体的生成过程；图中，
“■”
标记壁面单元，数字为外接长方体所含的壁面单元数。左上子图即未分割图中的矩形即为初始外接长方体，包含304个壁面单元。根据步骤s2
‑
2，n取10，共需分割9次。最终生成的10个外接长方体如图2的最下图所示。
[0031]
实施例2本实施例计算m6机翼，图3对比了本发明与穷举搜索法确定的m6机翼的壁面距，图中云图为本发明结果，等值线为穷举搜索法结果。图3中，不同截面的云图与等值线均完全重合，因此本发明计算的壁面距结果正确可信。对于计算效率，穷举搜索法计算29.5万网格单元的壁面距需耗时215.2s，本发明计算相同网格仅需3.7s，耗时缩短98.3%。
[0032]
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。