技术领域
[0001] 本发明属于直升机旋翼系统动力学仿真分析技术领域，涉及一种基于三维实体有限元部件级共旋坐标法(C‑CR)的直升机旋翼系统动力学仿真高效求解方法。
相关背景技术
[0002] 直升机作为一种能够垂直起降和悬停的飞行器，在军事、救援、运输和勘探等领域具有广泛的应用。直升机的核心部件之一是旋翼系统，它是提供升力和推力的关键组件。旋翼的性能直接影响直升机的飞行稳定性、操纵性和安全性。因此，对旋翼系统的结构动力学进行深入研究具有重要意义。
[0003] 旋翼系统主要由旋翼桨叶、桨毂、旋翼桨杆和旋翼驱动系统等组成。在飞行过程中，旋翼系统承受着复杂的动态载荷，如离心力、惯性力和弹性变形等。这些载荷相互作用，使得旋翼系统的动力学特性非常复杂。为了确保直升机的安全性和可靠性，必须对主旋翼系统的结构动力学进行精确的仿真分析。传统的分析方法主要包括以下几方面：
[0004] (1)多体动力学分析：模拟旋翼系统中各部件的相对运动及其相互作用，分析其整体动力学特性。
[0005] (2)有限元分析：对主旋翼系统进行精细化建模，分析其在复杂载荷作用下的应力应变分布及疲劳寿命。
[0006] (3)实验测试：通过物理试验验证仿真模型的准确性和可靠性，但这种方法通常耗时且成本高。
[0007] 随着计算机技术和仿真技术的发展，采用计算机仿真技术进行主旋翼系统的结构动力学分析成为研究的重点方向。仿真技术能够在设计阶段对主旋翼系统的动态特性进行预测和优化，减少物理试验的次数和成本，提高设计效率和产品性能。当前的仿真分析方法虽然在一定程度上解决了旋翼系统动力学分析的问题，但在处理复杂的非线性动态响应和高精度数值模拟等方面仍面临诸多挑战。例如，旋翼桨叶在旋转过程中会产生复杂的离心力和惯性力，导致结构的非线性响应，这对仿真分析提出了更高的要求。
[0008] 在直升机旋翼系统仿真分析领域，通常采用梁单元对桨叶结构进行模拟。然而，随着直升机构形的不断更新迭代，出现了更多具有先进构形和复杂几何外形的桨叶。梁单元在面对这些复杂的直升机旋翼结构时，难以满足精确分析的需求，尤其在桨毂和轴承等附属构件的建模及分析方面存在较大困难。因此，开发一种高效、准确的旋翼系统结构动力学仿真分析方法，具有重要的理论意义和工程应用价值。
[0009] 为实现直升机旋翼系统的高保真动力学分析，基于三维实体单元的共旋坐标法(Full‑CR)是解决这一问题的最佳途径。三维实体单元可以满足任何复杂旋翼几何外形的全尺寸建模分析，并且其完备的后处理分析技术可以获得直升机旋翼结构的应力及应变分布状态。这对先进翼型的设计与优化具有重要意义。然而，采用三维实体单元进行计算面临的主要问题是：
[0010] (1)单元数量庞大：为了精确地描述复杂几何形状和捕捉细微结构特征，模型需要划分成大量的三维实体单元。这导致模型规模急剧增大。
[0011] (2)计算资源消耗大：由于单元数量多，计算所需的存储和内存需求显著增加，对计算资源要求高。
[0012] (3)计算时间长：庞大的计算规模导致计算时间大幅延长，尤其在非线性动力学分析中，迭代求解过程复杂且耗时。
[0013] (4)数值稳定性问题：在求解过程中，特别是在处理非线性和复杂接触问题时，容易出现数值稳定性问题。
具体实施方式
[0110] 以下结合具体实施例对本发明
[0111] 结合直升机旋翼动力学分析算例，展开具体实施方式叙述并进行该方法的计算精度和计算效率对比分析。本发明的具体实施方式如下：
[0112] 第一步、将直升机旋翼系统几何模型通过有限元网格划分离散成三维实体有限元模型
[0113] 在这个算例中，附图8给出了直升机倾转旋翼系统的几何模型和三维实体有限元模型。所述几何模型共包含三片桨叶，单片桨叶的长度为4.5m，由铝制成，并沿展向具有‑12°的预扭角。旋翼系统桨毂为无铰式刚性桨毂。所述旋翼系统的三维实体有限元模型包含
3450个单元和6090个节点，自由度总数为18270个。将以上有限元数据生成模型节点坐标和单元节点编号模型文件，并将单元材料属性赋予到的单元当中。
[0114] 第二步、划分三维实体有限元模型的共旋分析部件
[0115] 在第一步将直升机旋翼系统几何模型通过有限元网格划分离散成三维实体有限元模型后，将三维实体有限元模型分割成若干个部件。将所述算例模型划分为6个共旋分析部件，每个部件内的单元共享同一坐标系，该坐标系被称作部件共旋坐标系。所述部件共旋坐标系的生成方式是首先在每个部件上寻找到三个相距较远且不在同一条直线上的三个节点，称之为部件代表节点。通过三个部件代表节点的边平行方法快速构造出相应的部件共旋坐标系。
[0116] 所述部件共旋坐标系基矢量的构造方法如下：
[0117]
[0118] 其中，e1表示由部件上的第二个代表节点与第一个代表节点形成的单位向量；e3表示由部件上的第三、第一代表节点形成的向量与e1向量叉乘得到的单位向量；e2表示由单位向量e3和e1叉乘得到的单位向量；x1表示第一个代表节点的坐标列向量；x2表示第二个代表节点的坐标列向量；x3表示第三个节点的坐标列向量；||…||表示计算向量的模；
[0119] 部件的旋转矩阵可表示为 部件共旋坐标系的建立省去了在传统单元级共旋坐标法中对每个单元共旋坐标系的求解过程，在很大程度上降低了程序计算求解时间。
[0120] 第三步、采用时不变约束或时变约束实现各个部件的有效连接。
[0121] 在第二步将直升机旋翼系统三维实体有限元模型被分割成若干部件后，采用多体系统动力学思想，通过约束实现部件间的对接。所述的约束包括时不变约束和时变约束。在本算例中，同一片桨叶上的两个部件通过时不变约束实现再次对接。部件1、3和5分别代表每片桨叶中桨根处的部件编号，这些部件通过时变约束进行对接，形成了等效形式的刚性桨毂结构。
[0122] 所述时变约束与时不变约束具体如下：
[0123] 3.1时不变约束的约束方程以及约束雅可比矩阵计算
[0124] 所述时不变约束可实现各部件之间的有效对接。通过耦合节点自由度的方法对相邻两个部件之间的有限元节点进行耦合。耦合节点自由度的方法是通过限制两个部件之间每对重合节点在空间三个维度上的相对坐标来实现的。耦合节点自由度方法所产生的约束方程与多体动力学中球铰约束的约束方程相同，约束方程可表示为Φ0＝0，其中下标0是代表与时不变约束相关的项。在这个方程中以部件的节点坐标 作为部件的广义坐标在得到球铰约束的约束雅可比矩阵Φ0后，通过将约束方程对部件广义坐标求导可计算时不变约束的约束雅可比矩阵Φq_0。Φq_0是一个稀疏的常数矩阵，仅在与某个部件相关的广义坐标位置上呈现1或‑1的值，其他位置均为0。因此可预先计算并存储所有球铰约束的雅可比矩阵，在进行仿真的每次迭代求解过程中直接回代这些矩阵，从而提高计算效率。
[0125] 3.2时变约束的约束方程以及约束雅可比矩阵计算
[0126] 所述的时变约束包括角度驱动约束、垂直约束和距离约束等约束类型。通过这些约束的组合与叠加可实现直升机桨毂结构与部件之间的有效对接。所述角度驱动约束用于控制三维实体有限元模型的旋转运动。所述垂直约束和距离约束用于约束三维实体有限元模型的垂直角度以及相对距离。为实现所述三种时变约束，需要在三维实体有限元模型中引入外置辅助节点。引入外置辅助节点目的是将辅助节点与三维实体有限元模型中的有限元节点联合，共同完成时变约束约束方程Φ1的列写。所述辅助节点的广义坐标记为qaux。通过qaux和在3.1节中所述的部件的广义坐标 在辅助节点与部件的三维实体有限元模型节点之间生成向量A、B、C、D。向量A表示桨毂结构中的支臂，向量C为直升机旋翼系统转轴，向量B和D属于直升机机体上的两个固定向量，作用是固定整个旋翼系统。这些向量用于约束方程Φ1的列写。角度驱动约束通过控制向量A与B夹角θ的变化来实现。所述的时变约束的约束方程Φ1的表达形式如下：
[0127] Φ1A＝A·B‑|A||B|cosθ (23)
[0128] Φ1V＝C·D (24)
[0129] Φ1D＝|A|‑d (25)
[0130] 其中，下角标1代表时变约束；Φ1A表示角度驱动约束的约束方程；θ表示驱动角度；V D
Φ1表示垂直约束；Φ1表示距离约束；d表示距离约束的控制距离。
[0131] 将时变约束的约束方程Φ1对广义坐标qaux和 求导可得到时变约束的约束雅可比矩阵Φq_1。时变约束的约束雅可比矩阵Φq_1具有高度的稀疏性。
[0132] 第四步、建立旋翼系统的动力学方程。
[0133] 在经过以上三步的处理后，旋翼系统的广义坐标q是所有三维实体有限元节点的广义坐标和引入辅助节点广义坐标的合集：
[0134]
[0135] 其中， 表示第一个部件内三维实体有限元节点广义坐标列向量； 表示第二个部件内三维实体有限元节点广义坐标列向量； 表示第Nb个部件内三维实体有限元节点广义坐标列向量； 表示辅助节点广义坐标列向量。
[0136] 利用第一类拉格朗日方程，并考虑所有的约束方程(23～25)，直升机旋翼系统的动力学控制方程可写成以下指标3的微分代数方程(DAEs):
[0137]
[0138] 其中，表示直升机旋翼系统中所有的广义坐标对应的广义加速度；Φ(q,t)表示直升机旋翼系统总的约束方程，包含了3.1节中的时不变约束约束方程Φ0和3.2节中的时变约束约束方程Φ1；M是直升机旋翼系统的质量矩阵；Φq＝[Φq_0；Φq_1]是系统约束雅可比矩阵；λ是物理上反映约束力的拉格朗日乘子向量。fINT是所有节点的内力列向量，节点的内力列向量由部件的弹性内力计算得到；fEXT是直升机旋翼系统外力，包括重力和气动力。
[0139] 进一步的，内力向量fINT和直升机旋翼系统的质量矩阵M可分别表示成：
[0140]
[0141] 其中， 表示第一个部件的质量矩阵； 表示第二个部件的质量矩阵；表示第Nb个部件的质量矩阵；Maux表示辅助节点质量矩阵。
[0142] 第五步、高效求解直升机旋翼系统动力学方程，并获得结构动力学响应。
[0143] 5.1计算直升机旋翼系统雅可比矩阵
[0144] 为了计算直升机旋翼系统的雅可比矩阵，本算例进行了计算参数设置。仿真过程的物理时间为10s，共划分成10000个时间积分步，每步步长0.001s。牛顿‑拉夫逊迭代的收‑3敛公差设置为5×10 。Newmark积分参数分别为0.56和0.25。旋翼最大转速为20rad/s，并在第一秒内从静止匀加速到最大转速。在t＝1.2s时，旋翼开始倾转，旋翼的倾转角速度为
0.157rad/s。
[0145] 进一步采用典型的隐式Newmark直接积分方法在每个积分/时间步上将系统动力学方程离散为一组非线性代数方程。如果得到了前第n个时间步长的解，则对微分代数方程进行离散化，并在第n+1步引入Newmark假设产生以下非线性方程：
[0146]
[0147] 式中，h为积分步长；α和δ表示Newmark算法中的两个参数；tn+1表示时间离散化后的第n+1个时刻； 表示t＝n+1时刻的加速度列向量；qn+1表示t＝n+1时刻的广义坐标列向量； 表示t＝n+1时刻的速度列向量。
[0148] 经推导后得到以 为未知数的非线非线性方程组。采用Newton‑Raphson方法得到迭代步中的线性方程组为：
[0149] JΔX＝GRES (31)
[0150] 其中，J是雅可比矩阵； 为迭代增量；GRES＝[FT,ΨT]T为非线性方程的残差；
[0151] 进一步的，所述雅可比矩阵J为：
[0152]
[0153] 雅可比矩阵可被简写为：
[0154]
[0155] 5.2线性方程组的快速求解
[0156] 在公式(33)中，C矩阵是大规模块对角矩阵，表示为
[0157]
[0158] 其中，
[0159]
[0160] 其中 为第二节中所述的部件的旋转矩阵。
[0161] 为实现大规模线性方程组(31)的快速求解，可使用舒尔补方法。首先将线性方程组(31)进行分块处理，经分块后，方程组(31)可以被改写为：
[0162]
[0163] 由式(36)的第一个方程可解出：
[0164]
[0165] 将公式(37)联立公式(36)的第二个方程，进而可解出非线性代数方程的拉氏乘子增量：
[0166]
[0167] 最后将公式(38)回代到公式(37)可得到加速度增量 的最终结果。 这一项是矩阵J中的子矩阵C的舒尔补。因此求解大规模线性化方程组问题通过舒尔补方法转‑1化为在每个迭代步中求解 和C F这两个关键项。
[0168] 第六步、计算线性化方程组求解过程中的两个关键项
[0169] 6.1高效计算第一个关键项C‑1F
[0170] 由公式(34)可知，大规模矩阵C具有块对角结构，直接使用并行模式求解可得到以下表达式
[0171]
[0172] 式中 是C(i)矩阵在进行乔列什基分解后所得到的下三角矩阵。在公式(39)中涉及到两类矩阵计算，一类是由旋转矩阵为子矩阵形成的大规模块状对角矩阵与一个列向量进行相乘，如表达式(18)中的第1步与第4步。这步计算可统一表示为diag(R)V。另一类是由一大规模稀疏下三角矩阵的逆矩阵与一列向量相乘，如表达式(39)中的第2步与第3步，可记为 这两类矩阵相乘的高效计算是提升算法计算效率的关键。具体操作流程如附图6所示。首先，在处理第一种矩阵相乘问题时，需将V3m×1向量利用Reshape函数将其转化成V3×m矩阵，下标m代表有限元节点总数。接下来再利用旋转矩阵R与其完成相乘运算。最后将得到的矩阵再次利用Reshape函数将其转化成为3m×1的列向量。
[0173] 在计算 形式的矩阵相乘时，下三角矩阵 可预先计算并储存，在每次迭代计算过程中只需要回代即可，不需要重新组装和分解。
[0174] 6.2高效计算第二个关键项
[0175] 在关键项 中，C矩阵的逆矩阵可以表示为
[0176]
[0177] 其中 和D均为块对角矩阵，由此可以得到以下表达式
[0178]
[0179] 由表达式(41)可以看出，只要计算出项 再将其转置矩阵与其自身进行矩阵乘法计算便可得到关键项二的结果。
[0180] 首先计算这一项中的时不变部分。如附图7所示，对于每一个部件，存在以下表达式
[0181]
[0182] 式中 代表在角标为m,n位置处的元素值为1而其它位置元素值为0的矩阵。由3.1节可知，公式(42)中的 完全是常数矩阵，可在迭代运算之前将其进行预先计算并储存。当在迭代过程中计算得到当前时刻某部件的旋转矩阵 后，代入公式(21)便可得到部件i的 逐一遍历全部部件，最后可
组装得到 以上所述算法避免了反复进行大规模稀疏矩阵计算组装及相乘运算所带来的时间消耗。
[0183] 其次考虑公式(41)中关于时变部分项 的计算。由3.2节可知，这一项是关于D和Φq_1的函数。所述时变约束部分项若采用预先计算并回代的方式在程序实现上并不一定会带来较高的计算效率。实际上含有双时变项的计算无需对其进行公式(21)的预处理。因为在直升机旋翼系统中 通常为一个十分稀疏且具有较少列数的矩阵，将这一项直接放在迭代计算中进行反复计算并不会对计算耗时产生较大的影响。因此，这一项的计算放在每个迭代步中完成即可。
[0184] 计算上述所有步骤后，可完成直升机旋翼系统的动力学微分代数方程(27)的求解。进而得到动力学时程响应。所述动力学时程响应包括全局位移、弹性位移、速度和加速度等。至此，完成了本发明提出的基于C‑CR方法的直升机旋翼系统动力学高效仿真求解方法。
[0185] 第七步、计算结果分析
[0186] 为了对算例结果进行对比分析，我们在商用有限元软件ANSYS中建立了相同的模型。刚性桨毂结构采用MPC‑184刚性梁单元建模，而旋翼旋转驱动铰和倾转铰则使用MPC‑184销轴连接单元建模。同时，还采用了单元级共旋坐标(Full CR)方法对结果进行了验证。
[0187] 附图9展示了采用部件级共旋坐标法计算的倾转旋翼系统的动力学仿真效果图。所述C‑CR算法可以计算得到旋翼桨尖节点在整个倾转过程中的运动轨迹。我们将该节点的空间运动轨迹、运动速度以及加速度分别投影到三个坐标轴方向，并与ANSYS计算结果以及单元级共旋坐标(Full CR)方法计算结果进行对比，相关曲线如附图10、附图11和附图12所示。计算结果表明，三种不同方法得到的桨尖节点的位移、速度和加速度结果基本一致，误差均在允许范围内。由于弹性效应的影响，加速度曲线经历了高频震荡以及较大的整体运动。附图13、附图14和附图15分别展示了三种计算方法所得桨尖节点x、y、z三个方向的局部弹性位移的时程曲线。这些曲线的振动趋势基本一致，结果验证了C‑CR计算方法的正确性。
[0188] 在计算效率方面，ANSYS耗时32030秒，单元级共旋坐标法(Full‑CR)耗时24066秒，而C‑CR方法耗时4358秒，约为ANSYS耗时的1/7，是Full‑CR算法耗时的1/5。由此可见，C‑CR方法在处理直升机旋翼系统动力学问题时具有较高的计算效率。
[0189] 以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制，应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。
