基于多学科耦合仿真的螺旋桨飞机桨叶飞脱安全评估方法

技术领域

本发明涉及螺旋桨飞机桨叶飞脱安全评估技术领域，具体而言，涉及一种基于多学科耦合仿真的螺旋桨飞机桨叶飞脱安全评估方法。

背景技术

近年来，螺旋桨飞机由于各种原因导致桨叶飞脱而造成的飞行事故呈高发态势。2011年美军一架KC130T飞机在空中飞行时，2号发动机的螺旋桨叶片由于老锈问题断裂脱落，击中飞机一侧后将飞机切分成3块，导致飞机坠毁，机上16人全部遇难；2014年11月6日，加拿大航空一架庞巴迪Q400由于轮胎爆裂紧急迫降时，机身右倾着地导致螺旋桨叶片飞脱，从窗口插入机身，重伤一名女乘客；2017年3月17日，澳洲航空的一架双发萨博340客机在空中右发螺旋桨飞脱，幸好螺旋桨叶片没有砍到机身，飞行员成功迫降。针对此类事故，我国在最新的运输类飞机适航条款CCAR25.905(d)中已明确规定，必须采取设计预防措施，使得螺旋桨叶损坏或由于桨毂损坏而松脱时，对飞机的危害减至最小，该条款中提出可提供分析/计算文件来表明符合性，分析/计算中需要研究螺旋桨飞脱叶片的运动轨迹，主要内容就是建立叶片失效模型，安全评估螺旋桨叶片损坏或由于桨毂损坏而松脱时对飞机的危害性影响。因此建立螺旋桨飞机桨叶飞脱安全评估方法是非常必要且急需的。

此前关于螺旋桨飞机桨叶飞脱的评估分析方法，由于具有一定的危险性和代价昂贵的缺陷，导致适航取证难度极大。美国联邦航空管理局(FAA)对此类条款的飞机适航符合性曾有豁免，后来即使有取证要求，也只是在地面通过单螺旋桨疲劳台架试验，或通过数学模型以工程估算的方式进行即可，相比真实飞行条件，其准确度仍有所欠缺，目前国内外未见在非定常滑流流场中实现飞脱叶片运动轨迹仿真的高可信度评估方法的公开资料报道。

发明内容

本发明旨在提供一种基于多学科耦合仿真的螺旋桨飞机桨叶飞脱安全评估方法，以解决上述的问题。

本发明提供的一种基于多学科耦合仿真的螺旋桨飞机桨叶飞脱安全评估方法，包括如下步骤：

步骤S1，初始化阶段：输入动力学参数、机体-桨叶重叠网格、以及流场CFD计算参数；

步骤S2，滑流计算阶段：基于动力学参数、机体-桨叶重叠网格、以及流场CFD计算参数，在多学科耦合分析框架下完成桨叶滑流非定常流场计算，得到桨叶飞脱初始位置的桨叶滑流非定常瞬时流场以及飞脱桨叶动力学参数；

步骤S3，桨叶飞脱阶段：基于桨叶飞脱初始位置的桨叶滑流非定常瞬时流场以及飞脱桨叶动力学参数，在多学科耦合分析框架下完成桨叶飞脱过程的非定常流场计算和多体动力学计算的耦合计算，得到桨叶飞脱非定常流场，获得桨叶飞脱过程中不同时刻的动力学参数并判断飞脱的桨叶是否与飞机其他部件发生碰撞。

进一步的，步骤S1包括：

(1)向多体动力学计算模块输入桨叶旋转和桨叶飞脱的动力学参数；

(2)向动态并行重叠网格模块输入机体-桨叶重叠网格完成重叠网格插值关系装配；

(3)向URANS并行非定常CFD求解器输入流场CFD计算参数，包括定常初始流场及其流场初始化参数。

进一步的，步骤S2包括：

(1)多体动力学计算模块根据输入的桨叶旋转的动力学参数，在多学科耦合分析框架下进行多体动力学计算，实现桨叶的约束旋转运动，得到飞脱桨叶和正常旋转桨叶的运动参数；

(2)动态并行重叠网格模块根据输入的机体-桨叶重叠网格以及飞脱桨叶和正常旋转桨叶的的运动参数寻找桨叶网格与机体网格的动态插值关系，分别得到各正常旋转桨叶的动态插值边界；

(3)URANS并行非定常CFD求解器通过动态插值边界完成桨叶网格与机体网格的流场信息交换，并在多学科耦合分析框架下对定常初始流场完成桨叶滑流非定常流场计算，得到桨叶滑流非定常流场，并获得收敛状态下的桨叶滑流非定常瞬时流场以及动力学参数，此时的桨叶滑流非定常瞬时流场以及飞脱桨叶动力学参数，即为桨叶飞脱初始位置下的桨叶滑流非定常瞬时流场以及飞脱桨叶动力学参数，并将其作为桨叶飞脱阶段的非定常流场计算和多体动力学计算的耦合计算的初始条件。

进一步的，步骤S3包括：

(1)多体动力学计算模块解除桨叶飞脱的旋转约束，并根据步骤S1输入的桨叶飞脱的动力学参数以及步骤S2计算得到的桨叶飞脱初始位置下的飞脱桨叶动力学参数进行多体动力学计算，实现桨叶飞脱的六自由度运动，得到飞脱桨叶运动参数，此过程中正常旋转桨叶仍然保持旋转约束；

(2)动态并行重叠网格模块根据输入的机体-桨叶重叠网格以及飞脱桨叶不同时刻运动参数和正常旋转桨叶此刻的旋转角度，来寻找桨叶网格与机体网格的动态插值关系，得到飞脱桨叶和正常旋转桨叶的动态插值边界；

(3)URANS并行非定常CFD求解器通过动态插值边界完成桨叶网格与机体网格的流场信息交换，并在多学科耦合分析框架下从桨叶飞脱初始位置下的滑流流场开始完成桨叶飞脱过程中不同时刻的流场计算，获得不同时刻桨叶飞脱的动力学参数并判断飞脱的桨叶是否与飞机其他部件发生碰撞。

进一步的，所述多学科耦合分析框架是指：

(1)建立多体约束系统的动力学微分-代数方程组进行非定常流场计算和多体动力学计算的耦合计算；

(2)在非定常流场计算和多体动力学计算的耦合计算中，在一个非定常流场计算时间步长内完成多步连续的多体动力学计算；

(3)在非定常流场计算和多体动力学计算的耦合计算中，进行四元数与欧拉角的相互转换。

进一步的，所述多体约束系统的动力学微分-代数方程组为：

其中，

进一步的，步骤S3中基于射线求交法实现飞脱的桨叶是否与飞机其他部件发生碰撞的自动检测。

具体地，所述基于射线求交法实现飞脱的桨叶是否与飞机其他部件发生碰撞的自动检测的方法为：基于几何光学原理，将飞脱的桨叶物面上的网格点P取做视点，通过模拟光的传播路径，对网格点P与飞机其它部件面C上的每个小面元进行求交运算，以此来确定网格点P与面C的交点个数，若交点为奇数，则网格点P位于面C内，说明飞脱的桨叶物面上的该网格点P与其它部件发生碰撞；遍历整个飞脱的桨叶物面，则能判断飞脱桨叶是否与飞机其它部件发生碰撞。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

本发明提出多学科耦合分析框架实现真实飞行条件下多片多轴、变转速螺旋桨桨叶在瞬时解除约束条件下的桨叶飞脱工况模拟和动力学仿真安全计算，能够准确计算桨叶的速度、轨迹、姿态、角速度等运动参数；并且发展飞脱桨叶与其他部件的自动碰撞检测能力，实现了碰撞瞬间部件参数的自动报告，能够为满足适航条款CCAR 25. 、905(d)的要求提供重要的技术支撑。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1为本发明的基于多学科耦合仿真的螺旋桨飞机桨叶飞脱安全评估方法的流程图。

图2为本发明的多学科耦合分析框架的原理图。

图3为本发明的基于射线求交法的碰撞自动检测示意图。

图4为本发明示例的机体-桨叶重叠网格的网格分布示意图。

图5为本发明示例的桨叶根据断裂示意图。

图6为本发明示例的质心位置示意图。

图7a为本发明示例的全机拉力系数计算收敛曲线图。

图7b为图7a的局部放大图。

图8为本发明示例的全机物面压力系数分布云图。

图9为本发明示例的全机空间涡量分布云图。

图10a为本发明示例的桨叶在45度位置开始飞脱时的轨迹图。

图10b为本发明示例的桨叶在45度位置飞脱0.00378秒后的轨迹图。

图10c为本发明示例的桨叶在45度位置飞脱0.00756秒后的轨迹图。

图10d为本发明示例的桨叶在45度位置飞脱0.01133秒后的轨迹图。

图10e为本发明示例的桨叶在45度位置飞脱0.01511秒后的轨迹图。

图10f为本发明示例的桨叶在45度位置飞脱0.01889秒后的轨迹图。

图10g为本发明示例的桨叶在45度位置飞脱0.02267秒后的轨迹图。

图10h为本发明示例的桨叶在45度位置飞脱0.02644秒后的轨迹图。

图10i为本发明示例的桨叶在45度位置飞脱0.03022秒后的轨迹图。

图10j为本发明示例的桨叶在45度位置飞脱0.03412秒后的轨迹图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例

如图1所示，本实施例提出一种基于多学科耦合仿真的螺旋桨飞机桨叶飞脱安全评估方法，包括如下步骤：

步骤S1，初始化阶段：

(1)向多体动力学计算模块输入桨叶旋转和桨叶飞脱的动力学参数；其中：

桨叶旋转的动力学参数包括桨叶转轴、桨叶转速等；

桨叶飞脱的动力学参数包括桨叶重量、飞脱初始时刻转动惯量、桨叶质心位置等；

(2)向动态并行重叠网格模块输入机体-桨叶重叠网格完成重叠网格插值关系装配；其中，

机体-桨叶重叠网格是指包括螺旋桨桨叶的全机网格，在全机网格中，对全机除螺旋桨桨叶以外的部件及空间流场区域，均采用多块对接网格绘制；对螺旋桨桨叶采用重叠网格技术绘制，来降低局部部件的绘制难度。

(3)向URANS并行非定常CFD求解器输入流场CFD计算参数，包括定常初始流场及其流场初始化参数（来流速度、飞行高度、迎角、湍流模型等），通过将流场初始化参数赋值给定常初始流场能够得到具有流场初始化参数的定常初始流场并参与后续的非定常流场计算。

步骤S2，滑流计算阶段：

(1)多体动力学计算模块根据输入的桨叶旋转的动力学参数，在多学科耦合分析框架下进行多体动力学计算，实现桨叶的约束旋转运动，得到所有桨叶（包括飞脱桨叶和正常旋转桨叶）运动参数；

(2)动态并行重叠网格模块根据输入的机体-桨叶重叠网格以及飞脱桨叶和正常旋转桨叶的运动参数寻找桨叶网格与机体网格的动态插值关系，分别得到各正常旋转桨叶的动态插值边界；

(3)URANS并行非定常CFD求解器通过动态插值边界完成桨叶网格与机体网格的流场信息交换，并在多学科耦合分析框架下对定常初始流场完成桨叶滑流非定常流场计算，得到桨叶滑流非定常流场，并获得收敛状态下的桨叶滑流非定常瞬时流场以及动力学参数，此时的桨叶滑流非定常瞬时流场以及飞脱桨叶动力学参数（质量、质心和转动惯量由步骤S1输入，六自由度信息由此时刻的流场计算得到），即为桨叶飞脱初始位置下的桨叶滑流非定常瞬时流场以及飞脱桨叶动力学参数，并将其作为桨叶飞脱阶段的非定常流场计算和多体动力学计算的耦合计算（MBD/CFD耦合计算）的初始条件。

如图2所示，上述的多学科耦合分析框架是指：

1、建立多体约束系统的动力学微分-代数方程组进行非定常流场计算和多体动力学计算的耦合计算；

其中，多体动力学计算作为计算力学的一个分支，主要用于复杂动力学系统的建模与仿真。多体动力学计算的研究对象一般包括大范围空间运动的刚体和柔性体两种类型。本发明中多体动力学计算模块基于多体系统的笛卡尔描述方法，使用分析力学中的拉格朗日方程建立多体约束系统的动力学微分-代数方程组：

其中，

螺旋桨桨叶旋转和飞脱过程动力学仿真过程分为两个阶段，其中第1阶段为桨叶飞脱前，此时各个螺旋桨绕轴以给定转速旋转。每个螺旋桨桨叶包括两个约束：

约束1：螺旋桨桨叶参考点与固定的转轴约束点之间的球铰约束，约束方程为（3个自由度）：

其中

约束2：给定螺旋桨旋转角速度随时间的规律（运动学约束），相对于固定坐标系，3个自由度的角速度约束方程为：

其中

若第

上述多体约束系统的动力学微分-代数方程组分为两部分：由系统动力学微分方程组和系统约束代数方程组构成典型的微分-代数方程组（Differential-AlgebraicEquations，DAEs）。本发明使用了隐式多步向后差分法（Backward Differential Formula，BDF）将方程中的广义坐标的微分导数项转化为代数式：

其中，

2、在非定常流场计算和多体动力学计算的耦合计算中，在一个非定常流场计算时间步长内完成多步连续的多体动力学计算；

在非定常流场计算和多体动力学计算的耦合计算中，非定常流场计算的时间步长Dt

因此，在一个非定常流场计算时间步长内完成多步连续的多体动力学计算。为了提高计算精度，

式中，

3、在非定常流场计算和多体动力学计算的耦合计算中，进行四元数与欧拉角的相互转换；

在多体动力学计算中采用四元数方法的优点是可以避免刚体转动过程中出现的奇异点，而在非定常流场计算的网格转动算法则基于欧拉角

步骤S3，桨叶飞脱阶段：

(1)多体动力学计算模块解除桨叶飞脱的旋转约束，并根据步骤S1输入的桨叶飞脱的动力学参数以及步骤S2计算得到的桨叶飞脱初始位置下的飞脱桨叶动力学参数进行多体动力学计算，实现桨叶飞脱的六自由度运动，得到飞脱桨叶运动参数，此过程中正常旋转桨叶仍然保持旋转约束；

(2)动态并行重叠网格模块根据输入的机体-桨叶重叠网格以及飞脱桨叶不同时刻运动参数和正常旋转桨叶此刻的旋转角度，来寻找桨叶网格与机体网格的动态插值关系，得到飞脱桨叶和正常旋转桨叶的动态插值边界；

(3)URANS并行非定常CFD求解器通过动态插值边界完成桨叶网格与机体网格的流场信息交换，并在多学科耦合分析框架下从桨叶飞脱初始位置下的滑流流场开始完成桨叶飞脱过程中不同时刻的流场计算，获得不同时刻桨叶飞脱的动力学参数并判断飞脱的桨叶是否与飞机其他部件发生碰撞。

当桨叶飞脱后需要重点关注它可能与飞机其它部件发生碰撞的位置，当发生碰撞时要能被有效检测同时终止后续的无效计算，所以需要在多体动力学计算和非定常流场计算（MBD/CFD）的多学科耦合分析框架下发展桨叶飞脱后与飞机其他部件发生碰撞的自动检测能力，本发明基于射线求交法实现飞脱的桨叶是否与飞机其他部件发生碰撞的自动检测。具体地：

如图3所示，基于几何光学原理，将飞脱的桨叶物面上的网格点P取做视点，通过模拟光的传播路径，对网格点P与飞机其它部件面C上的每个小面元(线段)进行求交运算，以此来确定网格点P与面C的交点个数，若交点为奇数，则网格点P位于面C内，说明飞脱的桨叶物面上的该网格点P与其它部件发生碰撞；通过该方法遍历整个飞脱的桨叶物面，则能判断飞脱桨叶是否与飞机其它部件发生碰撞。

示例：

S1，初始化阶段：飞机在0m高度的真空速484km/h，飞行迎角-3.2°，桨叶角39.31°，生成全机网格，在全机网格中，对全机除螺旋桨桨叶以外的部件及空间流场区域，均采用多块对接网格绘制；对螺旋桨桨叶采用重叠网格技术绘制，来降低局部部件的绘制难度，例如巡航构型共绘制了1108个网格，为了可以使任意的螺旋桨桨叶飞出，对每块桨叶都绘制了重叠网格，即24块重叠网格，重叠网格的网格单元总数约为1.0798亿，网格分布如图4。发动机转数1075r/min；桨叶断裂位置为桨叶根部断裂(见图5中画圈处)；桨叶质量：30.26kg；桨叶质心位置如图6所示。

S2，滑流计算阶段：进行桨叶未飞脱前的桨叶滑流非定常流场计算，通过求解非定常雷诺平均N-S方程，双时间的多重LU-SGS方法，空间格式采用Roe通量差分分裂方法，采用连续可微限制器，湍流模型为SST两方程模型。在非定常流场计算中，真实时间步长取为螺旋桨旋转5度对应的时间，即0.0007752s，真实时间步共计算720步，螺旋桨旋转10周，真实时间模拟0.5581s。图7a给出了对应的全机拉力系数（其中包括了桨叶拉力系数）计算桨叶滑流非定常流场的收敛曲线，可以看到螺旋桨旋转近两圈后全机气动特性计算趋于收敛，全机拉力系数计算收敛曲线局部放大图7b能看到明显的阻力周期性变化，说明桨叶滑流非定常流场计算可靠合理。图8显示了全机压力系数Cp的物面云图，可以明显看到机翼上表面被滑流影响的低压区。图9显示了以压力系数染色的螺旋桨空间涡量，考虑到桨叶飞出区域主要是在桨毂到机翼前缘这一区间，网格加密也集中于此。

S3，桨叶飞脱阶段：指定飞脱相位角从螺旋桨45º位置开始，这里给出每间隔两个计算步的桨叶飞行轨迹位置，见图10a～图10j，飞脱的桨叶在空中飞行约0.034s后桨尖与机体发生碰撞，当碰撞发生后，系统自动判断并终止计算机集群上的计算作业，并将碰撞瞬间飞脱桨叶的速度、姿态、角速度等运动参数信息以文本信息输出。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。