一种直升机流场数值模拟的优化方法及系统

技术领域

本发明属于直升机气动力技术领域，涉及一种基于直升机流场数值模拟的优化方法及系统。

背景技术

直升机全机的气动力是检验直升机飞行能力的重要手段。研究气动力问题主要有三大手段：理论分析(基于流动简化假设，给出所研究问题的解析解)、实验模拟(以地面实验为研究手段)和数值计算。数值模拟方法将基本定律(质量守恒定律、牛顿第二定律、能量守恒定律、气体状态方程等)以及相关的本构方程和状态方程，形成描述空气状态(密度、速度、压强)随时间变化率的偏微分方程组，然后将方程中的积分、微分项近似地表示为离散的代数形式，使得积分形式的控制方程转化为代数方程组，最后通过计算机求解这些代数方程组，从而得到流场在时间、空间上的数值解。本发明主要通过对求解代数方程的优化提高了整个数值模拟的性能。代数方程系数矩阵可用一个稀疏矩阵表示，我们使用多波前法(Multifrontal Method)来计算。在多波前法中，波前矩阵的计算是一个稠密矩阵计算问题，并且由于计算任务划分为多个波前矩阵的计算，因此涉及波前矩阵的计算过程成为了核心步骤。大型线性方程组的直接求解器由于其具有较好的数值鲁棒性、可靠性以及高可用性而被广泛的应用。但是，它也具有高计算复杂度，高内存占用率的缺点。这也限制了直接求解方法在大规模问题(具有数亿未知数的矩阵)上的应用范围。本技术的目标就是在保证鲁棒性、计算性能和高可用性的前提下减少代数方程求解的时间消耗，能够提高整个直升机流场数值模拟的性能。

发明内容

本发明的一个目的是为了解决上述问题，提出一种直升机流场数值模拟的优化系统，能够加速对直升机飞行过程中气流状态进行数值模拟分析，为直升机的空气动力学设计提供数据参考。整体流程分为以下四个模块：数据获取模块、数据预处理模块、流场求解模块和网格边界数据更新模块。本发明主要优化了流场求解模块的性能，能够提高其求解速率。流场求解模块主要的计算量在代数方程的求解的部分。结合SRBT(Symmetric RandomButterfly Transformation)的随机化处理与低秩矩阵的运用，最大程度上减小了运算次数。在瓦片(tile)化的数据布局下，减小计算任务粒度，提升并行程度，引入低秩矩阵，并改进了已有的ACA(Adaptive Crossing Approximation)算法，解决了ACA算法在压缩一些低秩特性不好的瓦片块时计算迭代次数过多的缺点，大大降低了操作数，并且可在该应用所需最低要求的低秩阈值条件下完成压缩操作，再通过OpenMP task construct进行并行化，在共享内存的平台上大大提升了计算速度。最后运用迭代精化算法提升最后结果的精确度，在保证一定精确度的前提下大大提升了流场的求解速度，加速了流场数值模拟的更新速度。

本发明系统包括数据获取模块、数据预处理模块、流场求解模块、网格边界数据更新模块。

所述数据获取模块用于获取直升机机身几何参数及飞行过程中气流物理参数信息，气流物理参数信息包括构建直升机飞行过程中来流速度、迎角等；

所述数据预处理模块用于构建气流网格拓扑结构，其中气流网格拓扑结构按照重叠网格法进行网格划分，每个网格包含当前位置气流的物理参数信息；

所述流场求解模块用于对气流网格拓扑结构各网格内的流场进行数值求解，得到流场数据；

所述网格边界数据更新模块用于根据流场求解模块获得的流场数据，对气流网格拓扑结构的网格节点进行更新，达到流场数值模拟的目的。

本发明的另一个目的是提供一种直升机流场数值模拟的优化方法。

步骤(1)、获取直升机机身几何参数及飞行过程中气流物理参数信息，气流物理参数信息包括构建直升机飞行过程中来流速度、迎角等；

步骤(2)、根据步骤(1)获取的数据构建气流网格拓扑结构，其中气流网格拓扑结构按照重叠网格法进行网格划分，每个网格包含当前位置气流的物理参数信息；

步骤(3)、气流网格拓扑结构中各网格进行独立的流场计算，获得流场空间分布的网格边界数据，以提供流场更新需要的数值；

步骤(4)、根据步骤(3)流场空间分布的网格边界数据，对气流网格拓扑结构进行更新。

步骤(5)、重复步骤(2)、(3)、(4)直至流场收敛。

本发明的又一个目的是提供一种计算机可读存储介质，其上存储有计算机程序，当所述计算机程序在计算机中执行时，令计算机执行上述的方法。

本发明的再一个目的是提供一种计算设备，包括存储器和处理器，所述存储器中存储有可执行代码，所述处理器执行所述可执行代码时，实现上述的方法。

本发明提供的技术方案包含以下有益效果：

保证精确度需求的前提下，本方法可以大幅度提升流场的求解速度，提升流场的预测效率，为直升机的空气动力学设计提供数据参考。

附图说明

图1为本发明系统的模块架构图；

图2为LDLT瓦片算法进行矩阵分解的流程图；

图3(a)、(b)分别为列优先数据布局与瓦片块状数据布局；

图4为LDLT分解在求解器中各任务实现流程及并行时数据相关性；

图5(a)-(b)分别为压缩前后全秩矩阵数据格式与低秩矩阵数据格式；

图6为改进ACA算法的流程图；

图7为网格边界数据更新模块的流程图。

具体实施方式

下面结合附图，对本发明的具体实施方案作进一步详细描述：

本发明基于直升机流场数值模拟的优化方法得到的流场空间分布为飞行器气动外形优化设计提供数据支持；

如图1本发明基于直升机流场数值模拟的优化系统包括数据获取模块、数据预处理模块、流场求解模块、网格边界数据更新模块。

所述数据获取模块用于获取直升机机身几何参数及飞行过程中气流物理参数信息，气流物理参数信息包括构建直升机飞行过程中来流速度、迎角等；

所述数据预处理模块用于构建气流网格拓扑结构，其中气流网格拓扑结构按照重叠网格法进行网格划分，每个网格包含当前位置气流的物理参数信息；

所述流场求解模块用于对气流网格拓扑结构各网格内的流场进行数值求解，得到流场数据；

所述网格边界数据更新模块用于根据流场求解模块获得的流场数据，对气流网格拓扑结构的网格节点进行更新，达到流场数值模拟的目的。

本发明基于直升机流场数值模拟的优化方法具体是：

步骤(1)、获取直升机机身几何参数及飞行过程中气流物理参数信息，气流物理参数信息包括构建直升机飞行过程中来流速度、迎角等；

步骤(2)、根据步骤(1)获取的数据构建气流网格拓扑结构；其中气流网格拓扑结构按照重叠网格法进行网格划分，每个网格包含当前位置气流的物理参数信息；

所述的气流网格拓扑结构通过重叠网格法建立，具体是：

(1)将流场分为若干个具有重叠部分的区域；

(2)在各个区域上分别生成网格结构；

(3)各区域独立求解；

(4)在重叠部分处通过网格间插值进行流场信息交换。

步骤(3)、气流网格拓扑结构中各网格进行独立的流场计算，获得流场空间分布的网格边界数据，以提供流场更新需要的数值；

3-1空气流动主要遵循四个方面的定律：运动学方面，动力学方面，热力学方面和气体物理与化学属性方面，包括质量守恒定律，牛顿第二定律，能量守恒定律，熵方程，气体状态方程等。基于上面基本定律以及相关的本构方程和状态方程，可形成描述空气状态(密度、速度、压强)随时间变化率的偏微分方程组，即为Navier-Stocks方程，简称N-S方程。在气流网格拓扑结构中各网格上建立三维可压雷诺平均N-S方程：

其中W为守恒变量，F(W)和F

Ax＝b (1)

系数矩阵A矩阵维度为n，nb为分块大小,x为方程未知量；

系数矩阵A是由大小为nb*nb的数据块组成，共有nt*nt个数据块；

提取代数方程的系数矩阵A；

3-2通过SRBT方法对系数矩阵A进行随机化预处理,并将矩阵A

A

A

其中T表示转置，U表示随机矩阵，随机化处理后A

随机矩阵U为递归蝶形矩阵：

U＝U

其中d代表深度，通常小于或等于2，U

其中B

其中R

e

瓦片状指的是将矩阵划分为网格单位，并以一个瓦片块为单位进行顺序存储，细化了任务粒度，提升了cache的使用效率；

图3(a)为现有的矩阵列优先存储的数据布局，只是按照数据实际的结构存储，并无关联数据在实际使用过程中提取的关联性，没有很好的重用cache，使得数据存取时，cache命中率低，因而降低存取性能。图3(b)为瓦片状的数据布局，由于在矩阵操作的过程中是以单个瓦片块为单位，因此以瓦片状的数据布局能够提高cache的命中率，充分利用cache的快速存取的特性。将矩阵A

3-3用多波前法分解瓦片(tile)状矩阵A

多波前法是对与稀疏分解相关的运算做了重新的组织，这种重组织以消去树的结构为基础，通过这种重组织把分解大型稀疏矩阵的任务转化成了一系列的子任务，而每个子任务则只包括对一个更小的稠密波前矩阵Q的局部分解；

3-4采用改进LDLT瓦片算法对波前矩阵Q进行LDLT分解，并利用OpenMP taskconstruct并行化，得到下三角矩阵L和对角矩阵D；

图2描述了采用改进LDLT瓦片算法分解的全过程，图4以nt＝3为例，即假定矩阵划分为3×3的块的情况下，在分解的全过程中，最外层一次循环过程中的任务关联性及执行次序。图中factor步骤即为以下步骤b)，solve步骤即为以下步骤c)，compress步骤即为以下步骤d)，update步骤即为以下步骤f)。图中的箭头表示各步骤直接的先后依赖性，同时也阐述了如何并行该分解过程。

所述LDLT分解是把矩阵Q分解为：

Q＝LDL

其中L、D矩阵中每个元素代表的是一个nb*nb的瓦片块。T是指对矩阵进行转置操作。

所述改进的LDLT瓦片算法具体如下：

a)初始化k＝1；

b)采用LDLT算法对波前矩阵Q对角块Q

其中1≤k≤nt,T表示转置；

c)根据上述得到的下三角矩阵L

W

其中k≤i≤nt。

d)采用改进ACA算法对瓦片块进行压缩处理；

图5中描述了压缩前后数据的格式变化，低秩压缩后的矩阵拥有存储空间占用率低，计算操作数低等特点，对于计算量占比最大的步骤e)来说，步骤e)前进行压缩处理，能大大提升步骤e)的计算性能。

图6中描述了改进ACA算法具体实现的流程。每一个循环都将根据ACA算法传统方式生成一组向量(一个行向量，一个列向量，分别并入矩阵X和Y中)来近似的表达整个矩阵的值，直到近似的精度达到低秩阈值，或者循环次数达到设定的循环阈值，则循环结束，压缩步骤结束。近似的精度可以通过低秩矩阵与原矩阵差的二范式来计算。传统的ACA算法在遇到低秩属性较差的矩阵时，会出现循环次数过多，操作时间过长，而且生成的向量组数过多的情况，对后续的计算优化也会不明显。因此我们设定了一个循环次数的阈值，来限制此情况的发生，当循环次数达到循环阈值时，我们将停止压缩，并输出结果，之后涉及该矩阵块的运算将采用原矩阵值而不是压缩后的矩阵。

因此，压缩算法步骤如下，将步骤c)第k列所有瓦片块的下三角矩阵L

其中，k+1≤i≤nt。

e)更新Q

判断步骤d)压缩后低秩矩阵X、Y的秩是否大于等于循环阈值，若是则使用原全秩矩阵L

若否则将公式(17)代入公式(18)，使用低秩矩阵更新Q

其中，k+1≤i≤nt，k+1≤j≤i

f)判断当前迭代下k是否为nt，若是则所有瓦片块分解结束，若否则更新k＝k+1，返回步骤b)。

3-5通过正向替代和反向替代得到最终方程组的解；

从上一步得到波前矩阵Q的分解结果

Q＝LDL

由一系列波前矩阵Q的局部分解结果综合得出原系数矩阵A的分解结果。

通过多波前法后续的计算得到原系数矩阵A的分解结果为

其中L

根据公式(4)可知

Ax＝U

结合公式(20)-(21)，可知：

求得中间变量y，然后将公式(1)带入(22)求得解x：

x＝Uy (23)

3-6使用一次迭代优化提升解的精确度。

3-6-1根据公式(23)得出的解x计算残差r：

r＝b-Ax (24)

3-6-2根据公式(24)得出的解x计算残差r:

Ad＝r (25)

其中d是修正值。

3-6-3更新修正值

x′＝x+d (26)

其中m代表迭代次数。

由于之前求得的解本身已有一定的精确度，并且考虑到执行迭代精化的时间成本比较高，因此只执行一次迭代精化操作，即上述步骤3-6-1至3-6-3只执行一次即可，因此x′即为所求解。

根据迭代优化更新后的解x′得到流场空间分布的网格边界数据。

步骤(4)、根据步骤(3)流场空间分布的网格边界数据，对气流网格拓扑结构进行更新；如图7具体是：

4-1动网格生成(网格节点更新)；

4-2更新几何数据(根据几何守恒定律计算控制面的法向速度)；

4-3重叠网格关系建立；

4-4外迭代子网格间插值。

步骤(5)、在非定场流体计算中，需要重复网格建立、网格内流场计算、网格更新的操作直至流场收敛。

因此，通过流场求解后的数据以及网格边界数据的更新可以达到直升机流场数值模拟的目的。