一种基于模态法的并行化高精度颤振计算方法

技术领域

本发明属于航空气动弹性领域，涉及一种基于模态法的并行化高精度颤振计算方法。

背景技术

颤振问题是气动弹性力学中的经典问题，也是飞行器设计人员普遍关注的一个问题：因为它严重的影响了飞行安全性，因此多数飞行器在气动与结构设计完成后都要进行颤振校核。

工程中常用的颤振计算方法是基于面元法气动力的频域颤振计算方法，该方法对于结构和气动均为线性的航空器经典颤振问题具有普遍适用性；但是，对于实际工程中常出现的低雷诺数流动、强分离流动、含有化学反应及相变的流动以及跨声速流动等非线性流动现象，工程颤振计算方法预测不准确，需要基于高精度的计算流体力学(CFD)方法进行颤振计算。因此需要发展一种能够考虑非线性气动力的高效颤振计算方法。

一种比较常见的解决思路是CFD/CSD(计算固体力学)耦合策略，在流固交界面处进行力和位移等信息交换，这种方法拥有较高的计算精度，但是计算成本很高，工程适用性差。

另一种解决方法是利用等效降阶模型辨识任意工况下的非线性气动力，再利用辨识模型来预测气动力，进而求解颤振问题；这种方法需要大量的CFD数据样本来训练降阶模型，时间成本不可小觑。同时降阶模型在任意飞行工况下的范化能力有待考验。通常情况下，一个降阶模型只能对某些工况下的非线性气动力预测有较好的效果。

发明内容

本发明为了解决能够考虑气动力非线性，同时兼顾计算精度和计算效率的工程颤振计算的难题，提出了一种基于模态法的并行化高精度颤振计算方法。

所述的颤振计算方法，具体步骤如下：

步骤一、使用流场网格划分软件对飞行器的流场域进行网格划分，并在fluent软件输出翼面网格节点坐标、网格中心点坐标、面单元全局编号以及累加节点数到txt文件中。

每个面单元都包括至少三个面网格节点和一个网格中心点。

步骤二、在结构有限元软件中输出结构有限元节点坐标及各阶模态位移；

步骤三、以并行方式运行fluent软件，通过薄板样条插值(TPS)方法将读取到的各阶模态位移值利用几何关系，插值到飞行器表面的各网格节点上，得到每个流场网格节点分别对应的各阶模态位移。

步骤四、host节点读取每个网格中心点分别对应的各阶模态位移以及txt文件，并利用各分线程节点node计算各阶广义气动力之和，得到整个机翼的广义气动力。

具体过程为：

步骤401、针对并行方式进行的每个面单元，host节点读取各网格中心点对应的各阶模态位移并传递给各分线程节点node；

各分线程节点node采用全局编号匹配搜索算法得到自身包含的所有网格中心点的各阶模态位移，具体过程如下：

首先，将飞行器表面的所有面单元全局编号向量faceGlobalId，按照串行局部编号k

然后，利用udf函数库中将局部编号转化为全局编号的函数F_ID，循环遍历各线程节点node内的所有面单元的串行局部编号，返回对应的全局编号F_ID(k

接着，将各面单元内的全局编号F_ID(k

最后，在进行广义气动力计算中，当计算出第k

步骤402、针对每个分线程节点node，该节点node调用自身包含的所有网格中心点的各阶模态值，分别乘以各网格中心点所在面单元的气动力得到该网格单元的各阶广义气动力，求和得到该节点node的各阶广义气动力之和；

在使用密度基求解各网格中心点所在面单元的气动力时，具体过程如下：

首先，飞行器表面的每个面单元对应一个体单元，获取各面单元的网格中心坐标以及体单元的网格中心坐标和压强。

然后，针对每个体单元，分别求出体单元中心到各面单元中心的压力梯度，以及体中心到各面中心的坐标向量。

最后，计算飞行器表面各面中心的压强值，乘以面单元的面积即为该面单元的气动力。

步骤403、分支节点node 0将并行方式进行的所有分线程节点node的气动力之和再次求和，得到各阶广义气动力之和。

步骤404、将各阶广义气动力之和经分支节点node 0，传输至host节点得到整个飞行器表面的广义气动力。

步骤五、利用整个飞行器表面的广义气动力，在host节点中求解用户自定义函数中编写的广义气动弹性运动方程组，得到各网格中心点对应的每阶模态的广义位移和广义速度，host节点将当前时间步下的各阶模态的广义位移和广义速度输出至txt文件，并将其传输至各分线程节点node；

步骤六、各分线程节点node执行户自定义函数：利用当前时间步下每阶模态的广义位移和广义速度，计算各网格中心点在飞行器表面的真实速度和位移；

步骤七、各分线程节点node利用各网格中心点的真实位移，控制各网格中心点所在面单元上的所有飞行器表面网格节点移动，并利用真实速度更新网格中心点的速度边界条件，完成该时间步的求解计算。

步骤八、判断求解计算是否到达终止时间步，如果是，进入步骤九，否则，返回步骤四；

步骤九、通过host节点输出的各个时间步的广义位移和广义速度的txt文件以及步骤二得到的结构有限元节点坐标、各阶模态数据计算各时间步的结构真实位移；

步骤十、监测不同来流速度下给定位置的结构真实位移-时间曲线，找到满足计算精度的颤振速度。

本发明与现有技术相比，具有以下优势：

(1)本发明一种基于模态法的并行化高精度颤振计算方法，通过在host节点向node节点传递数据的时候用到了全局编号匹配搜索技术，以及在用户自定义函数中编写广义气动弹性运动方程组；将结构方程求解过程写入用户自定义函数，避免了固体动力学方程求解，实现了在流体动力学软件中同时求解流场与结构动力学方程，进行时域颤振计算。大大提高了流固耦合计算速度，节约了计算成本。

(2)本发明一种基于模态法的并行化高精度颤振计算方法，通过对方法进行并行化处理，使得方法可以支持多核并行计算，进一步提升了求解的计算速度。

(3)本发明一种基于模态法的并行化高精度颤振计算方法，适用范围广，可以对不同速度(低速、跨音速以及超声速)范围下非线性气动力(跨声速、低雷诺数、大攻角等)参与的颤振问题进行求解。

附图说明

图1是本发明一种基于模态法的并行化高精度颤振计算方法的流程图；

图2是本发明在一个时间步内并行化的流固耦合计算过程图；

图3是本发明全局编号匹配搜索技术中不同编号方式下向量元素间的对应关系图；

图4是本发明体单元中心向面单元中心的外插值示意图；

图5是本发明在不同来流速度下AGARD 445.6翼尖处位移-时间曲线图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图对本发明作进一步的详细和深入描述。

本发明提出了一种针对非线性气动力的流固耦合颤振计算方法，适用于各种结构模型，如机翼，导弹和飞机，本发明利用fluent软件中的udf(用户自定义函数)功能，通过基于模态法的时域颤振求解算法，使得在非线性非定常气动力计算的同时能够求解结构变化参数，此外对程序进行了并行化处理，从而进行颤振分析，具有计算效率高、精度高和适用性广等优势。

本发明所述模态法的并行化高精度颤振计算方法，基于Ansys软件中的流体动力学计算和用户自定义函数功能，在用户自定义函数中开发颤振求解算法，具备位移插值模块、体压强到面压强转化模块、并行化处理模块、全局编号匹配搜索模块和结构动力学方程时域求解模块，能够应用于大攻角、分离流、低雷诺数、阻力舵等非线性气动力的工程颤振计算。

所述的颤振计算方法，如图1所示，具体步骤如下：

步骤一、使用流场网格划分软件对飞行器的流场域进行网格划分，并在fluent软件通过用户自定义函数输出翼面网格节点坐标、网格中心点坐标、面单元全局编号以及累加节点数到txt文件中。

每个面单元都包括至少三个面网格节点和一个网格中心点。

步骤二、在结构有限元软件中输出结构有限元节点坐标及各阶模态位移；

在基于Ansys用户自定义函数的位移插值模块中，读取结构有限元节点坐标、各阶模态位移、流体网格中的翼面网格节点坐标以及网格中心点坐标。

步骤三、以并行方式运行fluent软件，通过执行用户自定义函数中的薄板样条插值(TPS)方法将读取到的各阶模态位移值利用几何关系，插值到飞行器表面的各网格节点上，得到每个流场网格节点分别对应的各阶模态位移。

通过将薄板样条插值方法写入用户自定义函数，在fluent中实现了将结构有限元模型的模态位移插值到流场网格中，可以适用于任意有限元软件输出的模态格式；插值模态可以直接应用于Ansys或其他流体动力学软件的动网格边界。

步骤四、执行用户自定义函数使host节点读取每个网格中心点分别对应的各阶模态位移以及txt文件，并利用各分线程节点node计算各阶广义气动力之和，得到整个机翼的广义气动力。

主节点(host)：负责接收用户指令，并将其传递给第0分支节点node 0；同时，主节点也负责接收节点node 0从其余节点汇总的全部信息，处理并反馈给Cortex。在并行计算中，host节点中没有任何流场、网格相关信息，也不参与流体相关的求解计算。

分支节点node 0：接收来自主节点的指令，并且将其分配给自己以及其余分支节点；同时，节点node 0可以接受来自其他节点的数据，进行全局处理(全局和、差等)后反馈到主节点。在并行计算中，节点node 0除分配指令、汇总信息外，还承担流场的求解工作。该节点上存储着一个并行区块的流场、网格相关信息，并负责该区块流场的求解计算。

其他分支节点(node 1～n)：其他分支节点接收来自节点node 0的指令，同时也将自己的数据传输到节点node 0。在并行计算中，其他分支节点主要承担流场求解工作。每个分支节点上存储着一个对应并行区块的流场、网格相关信息，并负责该区块流场的求解计算。

在host节点向node节点传递数据时用到了全局编号匹配搜索技术以及在密度基求解计算气动力时用到的体单元压强与面单元压强转化技术；如图2所示，首先，将流体网格节点XYZ坐标写入文件1，将流体网格中心点XY坐标与对应面单元的全局编号、累加节点数写入文件2，host节点读取文件1、2中插值后的四阶模态位移传递给各分线程节点node，然后，每个分线程节点node，读取该节点上面单元的气动力，乘以所有网格中心点的各阶模态值，得到该节点的广义气动力；最后，将所有分线程节点node的广义气动力之和再次求和，得到广义气动力全局和，结合给定的广义气动弹性运动的初始广义位移、广义速度，或将上一步计算中得到的广义位移、广义速度结果设为初始值，用于求解广义气动弹性运动方程组，进一步计算飞行器表面各网格节点的真实速度，位移和面边界速度。

具体过程为：

步骤401、针对并行方式进行的每个面单元，执行用户自定义函数使host节点读取各网格中心点对应的各阶模态位移并传递给各分线程节点node；

各分线程节点node采用全局编号匹配搜索算法得到自身包含的所有网格中心点的各阶模态位移，具体过程如下：

首先，将飞行器表面的所有面单元全局编号向量faceGlobalId，按照串行局部编号k

向量faceGlobalId长度为Nc，按照串行局部编号k

然后，利用udf函数库中将局部编号转化为全局编号的函数F_ID，循环遍历各线程节点node内的所有面单元的串行局部编号，返回对应的全局编号F_ID(k

接着，将各面单元内的全局编号F_ID(k

最后，在进行广义气动力计算中，当计算出第k

全局编号匹配搜索算法由一个面单元循环构成，在每一个node i中执行：

其中，向量faceGlobalId是host节点读取的翼面的面单元全局编号向量，此信息包含在插值模块读取的翼面网格数据中。如图3所示，为各面网格单元编号之间的关系。由于UDF中全部节点循环是通过“面单元循环+面内节点循环”的嵌套实现的，因此只需用面单元的串行局部编号k

v

实际上，该方法适用于各种形式的混合网格，具有较好的通用性。

步骤402、针对每个分线程节点node，该节点node调用自身包含的所有网格中心点的各阶模态值，分别乘以各网格中心点所在面单元的气动力得到该网格单元的各阶广义气动力，求和得到该节点node的各阶广义气动力之和；

在使用密度基求解各网格中心点所在面单元的气动力时，具体过程如下：

首先，飞行器表面的每个面单元对应一个体单元，获取各面单元的网格中心坐标以及体单元的网格中心坐标和压强。

然后，针对每个体单元，分别求出体单元中心到各面单元中心的压力梯度，以及体中心到各面中心的坐标向量。

最后，计算飞行器表面各面中心的压强值，乘以面单元的面积即为该面单元的气动力。

体压强到面压强的转换是应用在计算飞行器表面的气动力上。

由于流体动力学求解器有两种求解方式：压力基与密度基。其中压力基求解适用于低速条件，密度基求解适用于压声速至超声速状态。压力基求解方式可以直接读取面网格中心的压强值，密度基求解无法实现。为了使本方法能够适用于密度基求解器从而进行高速状态下的颤振计算，网格面单元中心的压强值是通过其所在的体单元中心点的压强值做外插值获得的。其数学表达式如下所示：

其中，

其中，R为普适气体常量，M为气体摩尔质量，它们的比值

其中，

将该方程应用到体单元中心求解该点的压强梯度。同时，向量d

步骤403、执行用户自定义函数分使支节点node 0将并行方式进行的所有分线程节点node的气动力之和再次求和，得到各阶广义气动力之和。

步骤404、将各阶广义气动力之和经分支节点node 0，传输至host节点得到整个飞行器表面的广义气动力。

步骤五、利用整个飞行器表面的广义气动力，在host节点中求解用户自定义函数中编写的广义气动弹性运动方程组，得到各网格中心点对应的每阶模态的广义位移和广义速度，host节点将当前时间步下的各阶模态的广义位移和广义速度输出至txt文件，并将其传输至各分线程节点node；

在host节点向node节点传递数据的时候用到了全局编号匹配搜索技术，以及在用户自定义函数中编写广义气动弹性运动方程组，实现了在fluent流体计算的过程中同时求解结构运动方程。

在host中求解广义气动弹性运动方程组，得到各阶运动的广义位移和广义速度，并将其传输至各node i。其中host节点中广义气动弹性方程组的离散化处理采用了变步长的四阶龙格库塔方法，有较高的计算精度，亦可改为其他的时域数值积分方法。

步骤六、各分线程节点node执行用户自定义函数：利用当前时间步下每阶模态的广义位移和广义速度，计算飞行器表面各网格节点的真实速度和位移；

步骤七、各分线程节点node利用飞行器表面各网格节点的真实位移，控制所有飞行器表面网格节点移动，并利用真实速度更新网格节点的速度边界条件，完成该时间步的求解计算。

速度边界条件的作用是计算气动力。

步骤八、判断求解计算是否到达终止时间步，如果是，进入步骤九，否则，返回步骤四；

步骤九、通过host节点输出的各个时间步的广义位移和广义速度的txt文件以及步骤二得到的结构有限元节点坐标、各阶模态数据计算各时间步的结构真实位移；

步骤十、监测不同来流速度下给定位置的结构真实位移-时间曲线，找到满足计算精度的颤振速度。

如图5所示，横坐标表示时间，纵坐标表示后缘翼尖位移。从左到右位移幅值随时间分别为收敛、等幅、发散状态。由此可以判断该机翼的颤振速度约为175.3m/s。

表1计算的颤振速度与实验结果对比

本发明基于Ansys用户自定义函数功能的体压强到面压强转化方法，基于体网格中心数据的外插值方法，实现了由体压强数据得到面压强数据的功能，能够适用于高速流动的密度基气动力求解的颤振计算。基于Ansys用户自定义函数的程序并行化方法和全局编号匹配搜索方法。全局编号匹配搜索方法通过建立流场数据的全局编号与局部编号之间的联系，实现了并行化计算时各node节点与host节点间的数据提取与传递。在此基础上建立的程序并行化方法通过对颤振数据分发、分节点求解、收集的程序进行并行化处理，提高了运算效率，大大节省了运算时间。基于Ansys用户自定义函数功能的变步长龙格库塔算法的结构动力学方程求解。实现了在求解过程中按照整体截断误差对求解时间步进行调整，从而保证时域颤振计算精度。