CN116305589A - 一种直升机桨叶结构降阶分析方法、系统、设备及介质 - Google Patents

一种直升机桨叶结构降阶分析方法、系统、设备及介质 Download PDF

Info

Publication number
CN116305589A
CN116305589A CN202310571708.1A CN202310571708A CN116305589A CN 116305589 A CN116305589 A CN 116305589A CN 202310571708 A CN202310571708 A CN 202310571708A CN 116305589 A CN116305589 A CN 116305589A
Authority
CN
China
Prior art keywords
finite element
load
displacement
node
unit
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202310571708.1A
Other languages
English (en)
Other versions
CN116305589B (zh
Inventor
董雷霆
卢志远
黄业增
李明净
李书
贺天鹏
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN202310571708.1A priority Critical patent/CN116305589B/zh
Publication of CN116305589A publication Critical patent/CN116305589A/zh
Application granted granted Critical
Publication of CN116305589B publication Critical patent/CN116305589B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/15—Vehicle, aircraft or watercraft design
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Geometry (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开一种直升机桨叶结构降阶分析方法、系统、设备及介质，涉及结构多尺度分析技术领域，先建立直升机桨叶三维单胞的有限元网格模型，并利用多载荷工况对有限元网格模型进行变形仿真，计算三维单胞的等效属性矩阵。然后建立直升机桨叶的简化模型，将直升机桨叶所受的实际载荷施加在简化模型上，基于有限元仿真计算得到简化模型的各个单元的广义单元应变，基于单元的广义单元应变和等效属性矩阵计算每一单元的广义单元应力，以进一步计算每一单元的应力场。本发明利用结构降阶技术，能更好地处理截面不连续的桨叶，适用于复合材料、含横向加强等多种形式的桨叶结构。

Description

一种直升机桨叶结构降阶分析方法、系统、设备及介质
技术领域
本发明涉及结构多尺度分析技术领域，特别是涉及一种基于三维单胞的直升机桨叶结构降阶分析方法、系统、设备及介质。
背景技术
直升机桨叶是直升机旋翼系统的重要组成部分，主要用于提供升力和推进力，桨叶的叶型通常为翼型，其具有在不同飞行状态下均可提供最大升力的优点。桨叶的长度和弯曲变形会对直升机的飞行性能产生影响，长的桨叶可提供更大的升力和更高的旋转速度，弯曲的桨叶可使直升机更加灵活。桨叶通常由铝合金、复合材料或钛合金制成，这些材料具有高强度、轻量化和耐久性的特点。随着制造加工工艺的不断进步，现代直升机桨叶结构呈现出材料构成复杂、内部结构形式复杂的发展趋势，这使得桨叶结构细节的变形响应情况愈发难以准确获得，预测桨叶承受复杂耦合载荷作用下的结构响应也变得愈发困难。桨叶整体的宏观响应（如位移、转角等）直接影响了直升机飞行中的升力和操纵特性，桨叶细节部位的局部变形（如材料交界处的应力、复合材料铺层的应变等）与桨叶的结构寿命、维护保障相关，更精准的仿真预测有助于合理制定维护策略、结构优化设计等工作，提升直升机桨叶的设计、使用效能。
针对上述典型的桨叶结构多尺度分析问题，降阶方法是一种通过建立低维模型，将结构细节响应与低维结构单胞变形之间的关系信息保留，以起到快速获取低维模型力学响应、进一步分析局部响应作用的方法。对于桨叶这种细长型的结构，可以用一维梁单元来等效桨叶的三维结构单胞。在桨叶结构降阶分析领域，由Hodges等学者开发的变分渐进梁截面分析(VABS)软件是业界广泛使用的分析工具，可以对复合材料的桨叶结构进行降阶，得到等效截面属性，进而通过CAMRAD等软件分析桨叶的气弹等动力学特性。该方法获得了业界的普遍应用，不过该方法主要适用于截面连续无突变的桨叶结构，对于含横向加强（如肋）的或是含空腔的桨叶结构处理起来有一定困难。
基于此，亟需一种适用于含横向加强（如肋）的或是含空腔的桨叶的结构降阶分析技术。
发明内容
本发明的目的是提供一种直升机桨叶结构降阶分析方法、系统、设备及介质，可适用于含横向加强（如肋）的或是含空腔的桨叶的结构降阶分析，处理简单、精度高。
为实现上述目的，本发明提供了如下方案：
一种直升机桨叶结构降阶分析方法，所述结构降阶分析方法包括：
建立直升机桨叶三维单胞的有限元网格模型；所述直升机桨叶由多个所述三维单胞组装得到；
对于用于三维单胞变形仿真的多载荷工况中的每一载荷工况，基于所述载荷工况的工况载荷列向量计算有限元载荷列向量；以所述有限元载荷列向量为输入，利用含约束的有限元矩阵表达式计算有限元位移列向量；所述工况载荷列向量包括施加在所述有限元网格模型的第一虚拟节点和第二虚拟节点的载荷；所述有限元载荷列向量包括施加在所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的载荷；所述有限元位移列向量包括所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的位移；
基于所有所述载荷工况的有限元位移列向量计算所述三维单胞的等效属性矩阵；
建立所述直升机桨叶的简化模型；将所述直升机桨叶所受的实际载荷施加在所述简化模型上，基于有限元仿真计算得到所述简化模型的各个单元的广义单元应变；对于每一所述单元，基于所述单元的广义单元应变和所述等效属性矩阵计算所述单元的广义单元应力；所述简化模型包括多个节点和位于相邻两个所述节点之间的单元，每一所述单元代表一所述三维单胞；
对于每一所述单元，基于所述单元的广义单元应力计算所述单元的应力场。
一种直升机桨叶结构降阶分析系统，所述结构降阶分析系统包括：
单胞模型构建模块，用于建立直升机桨叶三维单胞的有限元网格模型；所述直升机桨叶由多个所述三维单胞组装得到；
多工况变形仿真模块，用于对于用于三维单胞变形仿真的多载荷工况中的每一载荷工况，基于所述载荷工况的工况载荷列向量计算有限元载荷列向量；以所述有限元载荷列向量为输入，利用含约束的有限元矩阵表达式计算有限元位移列向量；所述工况载荷列向量包括施加在所述有限元网格模型的第一虚拟节点和第二虚拟节点的载荷；所述有限元载荷列向量包括施加在所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的载荷；所述有限元位移列向量包括所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的位移；
等效属性矩阵计算模块，用于基于所有所述载荷工况的有限元位移列向量计算所述三维单胞的等效属性矩阵；
桨叶模型构建模块，用于建立所述直升机桨叶的简化模型；将所述直升机桨叶所受的实际载荷施加在所述简化模型上，基于有限元仿真计算得到所述简化模型的各个单元的广义单元应变；对于每一所述单元，基于所述单元的广义单元应变和所述等效属性矩阵计算所述单元的广义单元应力；所述简化模型包括多个节点和位于相邻两个所述节点之间的单元，每一所述单元代表一所述三维单胞；
应力场计算模块，用于对于每一所述单元，基于所述单元的广义单元应力计算所述单元的应力场。
一种直升机桨叶结构降阶分析设备，包括：
处理器；以及
存储器，其中存储计算机可读程序指令，
其中，在所述计算机可读程序指令被所述处理器运行时执行上述的直升机桨叶结构降阶分析方法。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述直升机桨叶结构降阶分析方法的步骤。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明用于提供一种直升机桨叶结构降阶分析方法、系统、设备及介质，先建立直升机桨叶三维单胞的有限元网格模型，并利用多载荷工况对三维单胞的有限元网格模型进行变形仿真，计算每一载荷工况的有限元位移列向量，以进一步计算三维单胞的等效属性矩阵。然后建立直升机桨叶的简化模型，将直升机桨叶所受的实际载荷施加在简化模型上，基于有限元仿真计算得到简化模型的各个单元的广义单元应变，基于单元的广义单元应变和等效属性矩阵计算每一单元的广义单元应力，以进一步计算每一单元的应力场。本发明利用结构降阶技术，通过分析桨叶的一小段三维单胞，建立高效的低自由度降阶模型，同时保持了结构细节与降阶模型之间的变形映射关系，很好地解决了桨叶这种承受分布法向载荷作用的结构在发生变形时整体结构响应高效仿真和结构细节应力便捷求解的问题，能更好地处理截面不连续的桨叶，适用于复合材料、含横向加强等多种形式的桨叶结构。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例1所提供的结构降阶分析方法的方法流程图；
图2为本发明实施例1所提供的结构降阶分析方法的原理框图；
图3为本发明实施例1所提供的节点结组示意图；
图4为本发明实施例1所提供的示例1的桨叶截面示意图；
图5为本发明实施例1所提供的示例1的桨叶单胞承受载荷
Figure SMS_1
作用时的变形结果云图；
图6为本发明实施例1所提供的示例1的桨叶简化模型的挠曲变形仿真结果示意图；
图7为本发明实施例1所提供的示例2的桨叶的有限元模型示意图；
图8为本发明实施例1所提供的示例2的桨叶单胞的有限元网格模型；
图9为本发明实施例1所提供的示例2的桨叶简化模型挠曲位移仿真结果与精细化模型对比图；
图10为本发明实施例1所提供的示例2的桨叶细节应力恢复结果云图；
图11为本发明实施例2所提供的结构降阶分析系统的系统框图。
符号说明：
1-不锈钢前缘；2-玻璃纤维大梁；3-玻璃纤维Z型梁；4-碳纤维增强复合材料蒙皮；5-泡沫填充。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种直升机桨叶结构降阶分析方法、系统、设备及介质，可适用于含横向加强（如肋）的或是含空腔的桨叶的结构降阶分析，处理简单、精度高。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
实施例1：
本实施例用于提供一种基于三维单胞的直升机桨叶结构降阶分析方法，如图1和图2所示，所述结构降阶分析方法包括：
S1：建立直升机桨叶三维单胞的有限元网格模型；所述直升机桨叶由多个所述三维单胞组装得到；
S2：对于用于三维单胞变形仿真的多载荷工况中的每一载荷工况，基于所述载荷工况的工况载荷列向量计算有限元载荷列向量；以所述有限元载荷列向量为输入，利用含约束的有限元矩阵表达式计算有限元位移列向量；所述工况载荷列向量包括施加在所述有限元网格模型的第一虚拟节点和第二虚拟节点的载荷；所述有限元载荷列向量包括施加在所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的载荷；所述有限元位移列向量包括所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的位移；
S3：基于所有所述载荷工况的有限元位移列向量计算所述三维单胞的等效属性矩阵；
S4：建立所述直升机桨叶的简化模型；将所述直升机桨叶所受的实际载荷施加在所述简化模型上，基于有限元仿真计算得到所述简化模型的各个单元的广义单元应变；对于每一所述单元，基于所述单元的广义单元应变和所述等效属性矩阵计算所述单元的广义单元应力；所述简化模型包括多个节点和位于相邻两个所述节点之间的单元，每一所述单元代表一所述三维单胞；
S5：对于每一所述单元，基于所述单元的广义单元应力计算所述单元的应力场。
本实施例提出了一种基于三维单胞的直升机桨叶结构降阶分析方法，包括步骤：（一）直升机桨叶单胞建模（与S1相对应）；（二）桨叶单胞多工况变形仿真（与S2、S3相对应）；（三）桨叶简化模型仿真（与S4相对应）；（四）桨叶单胞细节应力恢复（与S5相对应）。以下，本实施例将按照上述四个步骤的顺序对本实施例的结构降阶分析方法进行详细介绍：
（一）直升机桨叶单胞建模
对直升机桨叶的结构进行分析，选取一适当长度的桨叶段作为直升机桨叶的三维单胞，直升机桨叶可由多个三维单胞组装得到，即多个三维单胞周期性排列后可组装为直升机桨叶。
根据直升机桨叶的尺寸和材料布局，对三维单胞进行有限元建模，具体的，首先根据三维单胞的实际组成和形状，在有限元分析软件中绘制出三维单胞的截面，截面如图4所示，该截面垂直于直升机桨叶的桨叶轴线方向。然后将该截面沿轴线方向进行拉伸，即可得到三维单胞的有限元模型，该有限元模型的尺寸和材料与三维单胞的实际尺寸和实际材料保持一致。最后根据多组分材料的分布对有限元模型中每一材料所对应的结构分别进行网格划分，建立沿轴线方向扫略型的网格，得到三维单胞的有限元网格模型，有限元网格模型如图5所示，网格的形状可选用任意形状，尽量使得两个端部截面（即图5的前端面和后端面）具有相同的网格分布。本实施例中，前端面和后端面是指垂直于轴线方向的最前端的截面和最后端的截面。
本实施例还对有限元网格模型前端面的网格节点和后端面的网格节点进行节点配对，如图3所示，网格节点是指网格的顶点，并设置虚拟节点的自由度，虚拟节点包括第一虚拟节点、第二虚拟节点和多个第三虚拟节点，第三虚拟节点的数量与前端面的网格节点的数量相同。
节点配对过程包括：根据网格节点的对应关系进行节点分组，定义前端面的网格节点自由度为
Figure SMS_3
，所属组别为/>
Figure SMS_9
，包括网格节点数量为/>
Figure SMS_13
，后端面的网格节点自由度为
Figure SMS_5
，所属组别为/>
Figure SMS_7
，包括网格节点数量与前端面包括网格节点数量相同。根据三维坐标的定义，将前端面的第/>
Figure SMS_11
个网格节点在/>
Figure SMS_15
三个方向上的位移自由度分别表示为
Figure SMS_2
，将后端面的第/>
Figure SMS_8
个网格节点在/>
Figure SMS_12
三个方向上的位移自由度分别表示为
Figure SMS_16
，其中，/>
Figure SMS_4
方向为桨叶轴线方向，y方向和z方向为桨叶截面内的两个垂直方向，上标中的正负号代表网格节点所属组别，上标中的数值/>
Figure SMS_6
代表网格节点在组内的序号，相同
Figure SMS_10
坐标的/>
Figure SMS_14
组内节点的序号保持相同，下标代表网格节点位移自由度的方向。
虚拟节点的自由度的设置过程包括：定义第一虚拟节点
Figure SMS_19
，第一虚拟节点/>
Figure SMS_22
在
Figure SMS_26
三个方向上的位移自由度分别表示为/>
Figure SMS_20
；定义第二虚拟节点/>
Figure SMS_24
，第二虚拟节点/>
Figure SMS_28
在/>
Figure SMS_31
三个方向上的位移自由度分别表示为/>
Figure SMS_18
；根据节点组/>
Figure SMS_21
的节点数量/>
Figure SMS_25
，定义/>
Figure SMS_29
个第三虚拟节点/>
Figure SMS_17
，其中第/>
Figure SMS_23
个第三虚拟节点在/>
Figure SMS_27
三个方向上的位移自由度分别表示为/>
Figure SMS_30
。
（二）桨叶单胞多工况变形仿真
对桨叶单胞在15种不同载荷工况下的变形进行仿真，每个载荷工况下的仿真均包括约束方程施加、载荷工况施加与计算两个步骤。15种载荷工况包括基础工况和由基础工况组合得到的组合工况，基础工况包括拉伸、弯曲、扭转和剪切，即本实施例的15种载荷工况是根据桨叶单胞的拉伸、弯曲、扭转、剪切等6组基础工况组合得到的。
每个载荷工况下的载荷施加情况由施加在上述第一虚拟节点
Figure SMS_32
和第二虚拟节点
Figure SMS_33
各方向的单位作用力所组成的工况载荷列向量表示，所有载荷工况的工况载荷列向量以矩阵形式/>
Figure SMS_34
组合，所得到的载荷矩阵可表示为：
Figure SMS_35
；
其中，载荷矩阵的第t列即为第t个载荷工况的工况载荷列向量
Figure SMS_36
，其值/>
Figure SMS_37
表示施加在虚拟节点自由度/>
Figure SMS_38
中的第/>
Figure SMS_39
个自由度的载荷（
Figure SMS_40
），比如第t个载荷工况的工况载荷列向量中第2行元素值代表施加在/>
Figure SMS_41
的载荷。
本实施例还需要施加约束方程和边界条件，在有限元公式所建立的矩阵关系基础上，构造含约束的有限元矩阵表达式，后续则基于该含约束的有限元矩阵表达式进行各个载荷工况下的单胞的变形仿真。具体的，构造含约束的有限元矩阵表达式的过程可以包括：
对已配对的节点和节点自由度构造约束方程表达式，约束方程如下：
Figure SMS_42
；
Figure SMS_43
；
Figure SMS_44
；
其中，
Figure SMS_55
为有限元网格模型前端面的第i个网格节点在x方向的位移自由度；/>
Figure SMS_46
为有限元网格模型后端面的第i个网格节点在x方向的位移自由度；/>
Figure SMS_51
为有限元网格模型的第一虚拟节点在x方向的位移自由度；/>
Figure SMS_59
为第i个网格节点的y坐标，由于前后端面节点的y、z坐标一致，故在此不区分前后端面；/>
Figure SMS_63
为参考轴线的y坐标；/>
Figure SMS_62
为有限元网格模型的第二虚拟节点在z方向的位移自由度；/>
Figure SMS_65
为第i个网格节点的z坐标；/>
Figure SMS_56
为参考轴线的z坐标；/>
Figure SMS_60
为有限元网格模型的第二虚拟节点在y方向的位移自由度；/>
Figure SMS_47
为有限元网格模型的第i个第三虚拟节点在x方向的位移自由度；/>
Figure SMS_52
为有限元网格模型前端面的第i个网格节点在y方向的位移自由度；/>
Figure SMS_48
为有限元网格模型后端面的第i个网格节点在y方向的位移自由度；/>
Figure SMS_50
为有限元网格模型的第一虚拟节点在y方向的位移自由度；/>
Figure SMS_54
为有限元网格模型的第二虚拟节点在x方向的位移自由度；/>
Figure SMS_58
为有限元网格模型的第i个第三虚拟节点在y方向的位移自由度；/>
Figure SMS_57
为有限元网格模型前端面的第i个网格节点在z方向的位移自由度；/>
Figure SMS_61
为有限元网格模型后端面的第i个网格节点在z方向的位移自由度；/>
Figure SMS_64
为有限元网格模型的第一虚拟节点在z方向的位移自由度；/>
Figure SMS_66
为有限元网格模型的第i个第三虚拟节点在z方向的位移自由度。/>
Figure SMS_45
是需要输入的两个参数，其定义为进行均匀化分析的参考轴线的/>
Figure SMS_49
坐标和/>
Figure SMS_53
坐标，参考轴线可任意选定，参考轴线的选取只会影响中间变量，而不会影响应力场的最终计算结果。网格节点的坐标可以根据有限元网格模型和坐标系确定，其为已知值。
在构建得到上述约束方程之后，利用约束方程施加方法施加约束方程，即可构建得到含约束的有限元矩阵表达式。约束方程施加方法包括但不限于罚函数法和拉格朗日乘子法等常见的含约束的线性代数方程组计算方法。
本实施例利用罚函数法施加约束方程，罚函数方法为有限元分析中的常规约束方程施加方法，其形式为：
Figure SMS_67
；
其中，
Figure SMS_68
表示在含约束条件下的能量泛函；/>
Figure SMS_71
表示在不含约束条件下的能量泛函；/>
Figure SMS_73
为罚函数因子；/>
Figure SMS_70
代表有限元网格模型的端面；/>
Figure SMS_72
代表以/>
Figure SMS_74
向量为变量的约束方程，该约束方程要求满足/>
Figure SMS_75
；/>
Figure SMS_69
为约束方程的非齐次项。
根据上述构建的约束方程，对于节点组
Figure SMS_76
中的第/>
Figure SMS_77
个网格节点，可以得到其约束方程为：
Figure SMS_78
；
其中，
Figure SMS_79
为约束方程所转化成的非齐次多项式的系数矩阵；/>
Figure SMS_80
为约束方程所涉及的位移自由度所组成的变量。
系数矩阵可表示为：
Figure SMS_81
。
Figure SMS_82
。
Figure SMS_83
。
根据罚函数方法，构造含约束的有限元矩阵表达式如下：
Figure SMS_84
；
其中，
Figure SMS_94
为三维单胞的增广后有限元刚度矩阵，其行和列的元素个数均等于有限元网格模型的各个网格节点的各个位移自由度、第一虚拟节点的各个位移自由度和第二虚拟节点的各个位移自由度的总个数；/>
Figure SMS_88
为罚函数因子；/>
Figure SMS_90
为对约束方程的系数矩阵进行数值积分所得到的计算矩阵，其行和列的元素个数均等于有限元网格模型的各个网格节点的各个位移自由度、第一虚拟节点的各个位移自由度和第二虚拟节点的各个位移自由度的总个数；/>
Figure SMS_96
为有限元位移列向量，包括有限元网格模型的各个网格节点的各个位移自由度、第一虚拟节点的各个位移自由度和第二虚拟节点的各个位移自由度的位移；/>
Figure SMS_100
为有限元载荷列向量，包括施加在有限元网格模型的各个网格节点的各个位移自由度、第一虚拟节点的各个位移自由度和第二虚拟节点的各个位移自由度的载荷；/>
Figure SMS_98
为根据约束方程的非齐次项得到的修正载荷列向量，其是对各个第三虚拟节点的各个位移自由度的数值进行数值积分计算得到的。需要指出的是，/>
Figure SMS_102
、/>
Figure SMS_95
均为/>
Figure SMS_99
维矩阵，/>
Figure SMS_85
，/>
Figure SMS_91
为节点自由度总个数，/>
Figure SMS_89
为有限元网格模型所包含的网格节点的数量，2代表第一虚拟节点和第二虚拟节点，每个节点均有3个位移自由度，故需要乘以3。/>
Figure SMS_93
、/>
Figure SMS_97
、/>
Figure SMS_101
均为含有/>
Figure SMS_87
个元素的列向量。当给定有限元载荷列向量/>
Figure SMS_92
时，经过上述含约束的有限元矩阵表达式可以计算得到有限元位移列向量/>
Figure SMS_86
。
Figure SMS_103
为在桨叶原始三维单胞的/>
Figure SMS_104
维有限元刚度矩阵增加6个空白行（即元素值为0）和空白列所形成的增广后有限元刚度矩阵，有限元刚度矩阵可以在有限元分析软件中计算得到，具体对网格的几何参数（尺寸和位置）和材料参数（材料）进行数值积分得到。
Figure SMS_105
是对约束方程的系数矩阵进行数值积分所得到的计算矩阵，在分析过程中保持不变，其计算公式为：
Figure SMS_106
；
在利用上式计算得到节点组
Figure SMS_107
中第/>
Figure SMS_108
个网格节点对应的单个矩阵/>
Figure SMS_109
后，按照节点自由度的对应关系即可将节点组/>
Figure SMS_110
中所有网格节点对应的单个矩阵组装成计算矩阵
Figure SMS_111
，该组装过程在有限元软件内可以直接实现。
罚函数因子
Figure SMS_112
需要取一个较大的正数，以保证使/>
Figure SMS_113
矩阵的主元比/>
Figure SMS_114
矩阵的主元大/>
Figure SMS_115
量级，通过/>
Figure SMS_116
的取值能够减小仿真误差。
由于每个载荷工况下
Figure SMS_117
的数值均会发生变化，故/>
Figure SMS_118
在每个载荷工况下需要分别进行数值积分求得，其计算公式为：
Figure SMS_119
；
在利用上式计算得到第
Figure SMS_120
个第三虚拟节点对应的单个列向量/>
Figure SMS_121
后，按照节点自由度的对应关系即可将所有第三虚拟节点对应的单个列向量组装成修正载荷列向量/>
Figure SMS_122
，该组装过程在有限元软件内可以直接实现。
Figure SMS_123
向量由第三虚拟节点自由度/>
Figure SMS_124
的取值决定，其在不同载荷工况下的取值如下：
当
Figure SMS_125
时，/>
Figure SMS_126
；
当
Figure SMS_127
时，
Figure SMS_128
；
Figure SMS_129
；
Figure SMS_130
；
其中，
Figure SMS_132
为第4个载荷工况下，前端面/>
Figure SMS_138
的第i个网格节点的位移在/>
Figure SMS_142
方向的分量，其等于第4个载荷工况下的有限元位移列向量/>
Figure SMS_133
中第i个网格节点的第j个位移自由度的相应数值；/>
Figure SMS_136
为第一位移修正系数，其取值根据/>
Figure SMS_140
的值求得，/>
Figure SMS_143
为第t个载荷工况的工况载荷列向量中第3行的元素值；/>
Figure SMS_131
为第3个载荷工况下，前端面/>
Figure SMS_137
的第i个网格节点的位移在/>
Figure SMS_141
方向的分量，其等于第3个载荷工况下的有限元位移列向量/>
Figure SMS_144
中第i个网格节点的第j个位移自由度的相应数值；/>
Figure SMS_134
为第二位移修正系数，其取值根据/>
Figure SMS_135
的值求得，/>
Figure SMS_139
为第t个载荷工况的工况载荷列向量中第2行的元素值。
在得到上述含约束的有限元矩阵表达式后，即可进行施加载荷计算变形的步骤，具体对桨叶单胞的有限元网格模型施加当前工况的工况载荷列向量，使用单胞均匀化分析方法计算单胞模型（即有限元网格模型）的等效属性矩阵。具体的，单胞均匀化分析方法包括如下步骤。
（1）生成有限元载荷列向量
根据当前载荷工况的工况载荷列向量
Figure SMS_145
的数值，计算需要施加在节点组/>
Figure SMS_146
内网格节点自由度上的载荷列向量/>
Figure SMS_147
，上标/>
Figure SMS_148
代表第/>
Figure SMS_149
个载荷工况，/>
Figure SMS_150
共包含/>
Figure SMS_151
个分量，表达式为：
Figure SMS_152
；
其中，
Figure SMS_154
表示/>
Figure SMS_156
组内第/>
Figure SMS_158
个网格节点在/>
Figure SMS_155
方向上所施加的载荷，/>
Figure SMS_157
表示/>
Figure SMS_159
组内第/>
Figure SMS_160
个网格节点在/>
Figure SMS_153
方向上所施加的载荷，其在不同载荷工况下的取值如下，即载荷列向量的计算公式为：
当
Figure SMS_161
时，/>
Figure SMS_162
；
其中，t代表第t个载荷工况；
Figure SMS_163
为施加在有限元网格模型前端面的第i个网格节点的j方向上的载荷；i=1，2，...，G+，G+为端面的网格节点的数量；j=x，y，z；/>
Figure SMS_164
为施加在有限元网格模型后端面的第i个网格节点的j方向上的载荷；施加在有限元网格模型其余网格节点的载荷为0。
当
Figure SMS_165
时，
Figure SMS_166
；
Figure SMS_167
；
Figure SMS_168
；
Figure SMS_169
；
其中，
Figure SMS_186
表示当/>
Figure SMS_190
时前端面的网格节点/>
Figure SMS_194
上的等效节点力在/>
Figure SMS_170
方向上的分量，可通过/>
Figure SMS_174
求得，/>
Figure SMS_178
表示第4个载荷工况求解得到的有限元位移列向量，取/>
Figure SMS_182
中与节点组/>
Figure SMS_198
中网格节点/>
Figure SMS_202
的/>
Figure SMS_205
方向自由度相对应的值作为/>
Figure SMS_207
，即/>
Figure SMS_203
为第4个载荷工况下，有限元网格模型前端面的第i个网格节点的j方向上的等效节点力；/>
Figure SMS_206
为第一载荷修正系数，其根据/>
Figure SMS_208
的值求得，/>
Figure SMS_209
为第t个载荷工况的工况载荷列向量中第3行的元素值；/>
Figure SMS_188
为第3个载荷工况下，有限元网格模型前端面的第i个网格节点的j方向上的等效节点力，可通过/>
Figure SMS_192
求得，/>
Figure SMS_196
表示第3个载荷工况求解得到的有限元位移列向量，取
Figure SMS_200
中与节点组/>
Figure SMS_171
中网格节点/>
Figure SMS_176
的/>
Figure SMS_180
方向自由度相对应的值作为/>
Figure SMS_184
；/>
Figure SMS_173
为第二载荷修正系数，其根据/>
Figure SMS_175
的值求得，/>
Figure SMS_179
为第t个载荷工况的工况载荷列向量中第2行的元素值；/>
Figure SMS_183
为第4个载荷工况下，有限元网格模型后端面的第i个网格节点的j方向上的等效节点力，可通过/>
Figure SMS_187
求得，取/>
Figure SMS_191
中与节点组/>
Figure SMS_195
中网格节点/>
Figure SMS_199
的/>
Figure SMS_177
方向自由度相对应的值作为/>
Figure SMS_181
；/>
Figure SMS_185
为第3个载荷工况下，有限元网格模型后端面的第i个网格节点的j方向上的等效节点力，可通过/>
Figure SMS_189
求得，取/>
Figure SMS_193
中与节点组/>
Figure SMS_197
中网格节点/>
Figure SMS_201
的/>
Figure SMS_204
方向自由度相对应的值作为/>
Figure SMS_172
。施加在有限元网格模型其余网格节点的载荷为0。
第一虚拟节点和第二虚拟节点的各个自由度的载荷值与载荷工况列向量的相应值相等。
基于上述公式，即可确定第t个载荷工况下的有限元载荷列向量。
（2）施加单胞边界条件
以通用的刚体位移、刚体转动约束方式限制
Figure SMS_210
节点组所在端面的位移和转角为0，具体选择三维单胞中的6个位移自由度分别进行限制，删去这6个位移自由度在/>
Figure SMS_211
、/>
Figure SMS_212
、/>
Figure SMS_213
、/>
Figure SMS_214
、/>
Figure SMS_215
中对应的行和列。上述6个位移自由度的具体选取方法为，选择后端面/>
Figure SMS_216
中的一个网格节点a，限制其三个方向的位移自由度，选择与网格节点a的y坐标相同的一个网格节点b，限制其x方向和z方向位移，选择与网格节点a的z坐标相同的一个网格节点c，限制其y方向位移。
（3）单胞变形求解
将载荷工况列向量和（1）中计算得到的载荷列向量按照自由度对应关系带入到含约束的有限元矩阵表达式等号右边的
Figure SMS_217
向量，再给罚函数因子/>
Figure SMS_218
赋值，进行矩阵求逆运算，即可计算得到第/>
Figure SMS_219
个载荷工况下的有限元位移列向量/>
Figure SMS_220
。需要指出的是，15个载荷工况依次进行，前4个载荷工况的计算结果需要在计算后11个载荷工况时调用。
（4）计算单胞等效属性
根据桨叶单胞在不同工况下的变形数据，计算得到单胞的等效属性矩阵，具体的，读取各个载荷工况下的第一虚拟节点
Figure SMS_221
和第二虚拟节点/>
Figure SMS_225
的位移/>
Figure SMS_228
，其中，下标/>
Figure SMS_222
代表方向，/>
Figure SMS_226
代表载荷工况序号，/>
Figure SMS_229
为第t个载荷工况下，第一虚拟节点的j方向的位移，/>
Figure SMS_231
为第t个载荷工况下，第二虚拟节点的j方向的位移。分别计算各个载荷工况下的桨叶单胞的总应变能/>
Figure SMS_223
，其中，上标/>
Figure SMS_224
代表载荷工况序号，第t个载荷工况下的桨叶单胞的总应变能/>
Figure SMS_227
的计算方法为：对有限元网格模型的所有网格节点的位移（基于第t个载荷工况的有限元位移列向量即可确定）和材料进行数值积分，即可得到第t个载荷工况下的桨叶单胞的总应变能/>
Figure SMS_230
。该过程可在有限元分析软件内完成。
单胞的等效属性矩阵的计算公式为：
Figure SMS_232
；
其中，
Figure SMS_233
为等效属性矩阵；/>
Figure SMS_234
为分块矩阵，其为6×6矩阵。
分块矩阵
Figure SMS_235
可通过以下方式求得，将/>
Figure SMS_236
表示为分块矩阵：
Figure SMS_237
；
其中，
Figure SMS_238
为第一分块，其为4×4矩阵；/>
Figure SMS_239
为第二分块；/>
Figure SMS_240
为第三分块。
Figure SMS_241
基于所有载荷工况的有限元位移列向量计算得到，具体基于/>
Figure SMS_242
求得，选取4个载荷工况，每个载荷工况下再选取4个值，其表达式为：
Figure SMS_243
。
Figure SMS_244
基于所有载荷工况下的三维单胞的总应变能计算得到，计算公式为：
Figure SMS_245
；
其中，
Figure SMS_246
为第t个载荷工况下的三维单胞的总应变能，根据第t个载荷工况的有限元位移列向量计算得到；/>
Figure SMS_250
为第t个载荷工况的工况载荷列向量；/>
Figure SMS_253
为第t个载荷工况的工况载荷列向量中第2行的元素值；/>
Figure SMS_248
为/>
Figure SMS_249
中第4行第4列的元素值；/>
Figure SMS_251
为第t个载荷工况的工况载荷列向量中第3行的元素值；/>
Figure SMS_252
为/>
Figure SMS_247
中第3行第3列的元素值。
将所有载荷工况的三维单胞的总应变能依次代入上述计算公式，即可最终确定
Figure SMS_254
中的每一个元素的值。
（三）桨叶简化模型仿真
建立桨叶的低自由度简化模型，施加其在待求工况时的气动力、离心力等分析载荷，进行有限元仿真，该仿真过程包括如下步骤：
（1）简化模型建立
建立直升机桨叶的简化模型，该简化模型可以包括多个节点和位于相邻两个节点之间的单元，本实施例的单元为线单元，即该单元可用直线表示，单元的数量与直升机桨叶实际包括的三维单胞的数量相同，每一单元代表一三维单胞。
（2）广义单元应变求解
由于有限元法中针对不同节点数量的单元具有不同的形函数
Figure SMS_255
，此处不再具体推导简化模型仿真的刚度矩阵，而是基于形函数来进行仿真，具体的，在有限元求解过程中，将直升机桨叶所受的气动力、离心力等实际载荷以等效节点力的方式施加在简化模型的每一节点上，在桨叶的根部（一般为第一个节点）施加固支边界条件，在有限元分析软件中计算得到简化模型的各个节点的位移，进一步利用形函数求得简化模型的各个单元的广义单元应变/>
Figure SMS_256
，下标/>
Figure SMS_257
代表单元编号。其中，各个单元的广义单元应变的求解过程具体可为：对于每一单元，对该单元对应的形函数N进行求导，再乘以该单元对应的两个节点的位移，即可得到该单元的广义单元应变，该过程可在有限元分析软件中实现。
（3）广义单元应力求解
使用上述的等效属性矩阵
Figure SMS_258
，建立简化模型的单元的本构关系：
Figure SMS_259
；
其中，
Figure SMS_260
为简化模型的单元的广义单元应力，表达式为
Figure SMS_264
，其分量表示作用在该单元上的6种内力，分别为拉力、扭矩、绕/>
Figure SMS_267
轴弯矩、绕/>
Figure SMS_263
轴弯矩、沿/>
Figure SMS_265
方向剪力、沿/>
Figure SMS_268
方向剪力；/>
Figure SMS_270
为简化模型的单元的广义单元应变，表达式为/>
Figure SMS_262
，其分量表示该单元上的6种广义应变，分别为拉伸应变、扭转应变、绕/>
Figure SMS_266
轴曲率、绕/>
Figure SMS_269
轴曲率、沿/>
Figure SMS_271
方向剪切应变、沿/>
Figure SMS_261
方向剪切应变。
则对于每一单元，基于单元的广义单元应变和等效属性矩阵计算单元的广义单元应力可以包括：计算单元的广义单元应变和等效属性矩阵的乘积，得到单元的广义单元应力。
（四）桨叶单胞细节应力恢复
利用有限元仿真结果，计算并输出桨叶结构单胞的细节应力场，本实施例计算应力场的方法包括但不限于线性叠加方法以及将广义单元应力以边界条件形式带入有限元方程中的计算求解方法。
本实施例以线性叠加方法为例，介绍本实施例的应力恢复过程，基于单元的广义单元应力计算单元的位移场，并根据位移场计算单元的应力场，具体包括：
（1）位移场求解
对于每一单元，基于单元的广义单元应力计算单元的位移场，具体的，根据简化模型的第i个单元的广义单元应力
Figure SMS_272
，使用线性叠加方法，计算第i个单元的位移场/>
Figure SMS_273
，计算公式为：
Figure SMS_274
；
其中，
Figure SMS_275
为第i个单元的位移场；/>
Figure SMS_276
为第i个单元的广义单元应力；/>
Figure SMS_277
为预设矩阵；/>
Figure SMS_278
为第t个载荷工况的工况载荷列向量；/>
Figure SMS_279
为第t个载荷工况下，约束方程所涉及的位移自由度所组成的变量对应的位移列向量，具体按照/>
Figure SMS_280
的表达式将第t个载荷工况的有限元位移列向量的部分元素组成该位移列向量。
其中，
Figure SMS_281
矩阵的表达式如下：
Figure SMS_282
；
（2）应力场求解
根据位移场，在有限元分析软件内即可计算单元的应力场。
以下，通过两个具体示例对本实施例的结构降阶分析方法进行进一步的介绍：
（一）示例一
在本示例中，考虑一个截面无变化的直升机桨叶，其由复合材料加工制成，组分包括不锈钢前缘1、玻璃纤维大梁2、玻璃纤维Z型梁3、碳纤维增强复合材料蒙皮4和泡沫填充5，其截面示意图如图4所示，采用NACA0012翼型，桨叶长度为5000mm，弦长为260mm。分析其在受均布升力作用时的变形和细节应力，升力的合力为1300N，分析时降阶过程的参考轴线为弦线的1/4处。
其有限元网格通过扫略方式生成，取桨叶单胞的长度为20mm，单胞的有限元网格数量为12751个，网格采用8节点六面体单元。根据有限元关系矩阵，生成桨叶结构单胞在各工况下的约束方程，通过罚函数方法施加，对各载荷工况进行分析仿真，并读取分析结果，在各载荷工况分析仿真后，得到桨叶单胞在承受不同载荷时的变形结果，如图5所示，其为桨叶单胞承受载荷
Figure SMS_283
作用时的变形结果云图。
利用计算单胞等效属性步骤，计算桨叶单胞的等效属性矩阵如下表1：
表1桨叶单胞等效属性
Figure SMS_284
根据等效属性，建立桨叶简化模型（50个单元），将升力转化为等效节点力施加在简化模型上，桨叶根部施加固支边界条件，计算得到的挠曲变形仿真结果如图6所示，最大挠度与精细化有限元结果相比误差为0.6%，在保持较高精确度的同时，将模型分析量大大降低，进而进一步分析各个单胞的细节应力场、应变场。
（二）示例二
在本示例中，考虑一个内部含加强肋结构的直升机桨叶，在轴线方向上不连续，如图7所示，其材料为玻璃纤维。桨叶采用NACA2412翼型，长度为3600mm，弦长为300mm。分析其在受均布升力作用时的变形和细节应力，分析时降阶过程的参考轴线为弦线的1/4处，升力在参考轴线上可以等效为法向的分布线载荷和绕轴向的扭矩，线载荷密度为1.53N/mm，扭矩为160.200Nm。
其有限元网格通过扫略方式生成，取桨叶单胞的长度为120mm，单胞的有限元网格数量为5090个，网格采用8节点六面体单元，如图8所示。根据有限元关系矩阵，生成桨叶结构单胞在各工况下的约束方程，通过罚函数方法施加，对各载荷工况进行分析仿真，并读取分析结果，在各载荷工况分析仿真后，得到桨叶单胞在承受不同载荷时的变形结果，进一步求得等效属性矩阵。将等效属性矩阵带入简化模型，计算挠曲位移，挠曲位移仿真结果与精细化模型的挠曲位移仿真结果对比如图9所示，位于桨叶中段的单胞的细节应力恢复结果云图如图10所示。
本实施例所提供的一种基于三维单胞的直升机桨叶结构降阶分析方法，包括步骤：直升机三维桨叶单胞建模，桨叶单胞多工况变形仿真，桨叶简化模型仿真，桨叶单胞细节应力恢复，利用了结构降阶技术，通过分析桨叶结构的一小段三维单胞，建立高效的低自由度降阶模型，同时保持了结构细节与降阶模型之间的变形映射关系，特别是受剪力作用下的变形结果具有较好精度，很好地解决了桨叶这种承受分布法向载荷作用的结构在发生变形时整体结构响应高效仿真和结构细节应力便捷求解的问题，将复杂拓扑、复杂材料组分的细长结构进行简化建模，处理成仅包含少量梁单元的仿真模型，显著提升了结构分析与仿真的效率。相较于目前普遍使用基于二维截面分析的桨叶结构降阶分析方法，对于沿展向无明显变化的桨叶结构，本实施例的结构降阶分析方法与当前方法具有相近的分析精度，而对于包含横向加强单元的复杂桨叶结构，当前方法不再适用，而本实施例的结构降阶分析方法仍能得到可靠的分析结构。即相比于现有的降阶分析工具，本实施例的结构降阶分析方法能更好地处理截面不连续的桨叶，具有良好的通用性和较高的保真度，适用于复合材料、含横向加强等多种形式的桨叶结构。
实施例2：
本实施例用于提供一种基于三维单胞的直升机桨叶结构降阶分析系统，如图11所示，所述结构降阶分析系统包括：
单胞模型构建模块M1，用于建立直升机桨叶三维单胞的有限元网格模型；所述直升机桨叶由多个所述三维单胞组装得到；
多工况变形仿真模块M2，用于对于用于三维单胞变形仿真的多载荷工况中的每一载荷工况，基于所述载荷工况的工况载荷列向量计算有限元载荷列向量；以所述有限元载荷列向量为输入，利用含约束的有限元矩阵表达式计算有限元位移列向量；所述工况载荷列向量包括施加在所述有限元网格模型的第一虚拟节点和第二虚拟节点的载荷；所述有限元载荷列向量包括施加在所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的载荷；所述有限元位移列向量包括所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的位移；
等效属性矩阵计算模块M3，用于基于所有所述载荷工况的有限元位移列向量计算所述三维单胞的等效属性矩阵；
桨叶模型构建模块M4，用于建立所述直升机桨叶的简化模型；将所述直升机桨叶所受的实际载荷施加在所述简化模型上，基于有限元仿真计算得到所述简化模型的各个单元的广义单元应变；对于每一所述单元，基于所述单元的广义单元应变和所述等效属性矩阵计算所述单元的广义单元应力；所述简化模型包括多个节点和位于相邻两个所述节点之间的单元，每一所述单元代表一所述三维单胞；
应力场计算模块M5，用于对于每一所述单元，基于所述单元的广义单元应力计算所述单元的应力场。
实施例3：
本实施例用于提供一种直升机桨叶结构降阶分析设备，包括：
处理器；以及
存储器，其中存储计算机可读程序指令，
其中，在所述计算机可读程序指令被所述处理器运行时执行如实施例1所述的直升机桨叶结构降阶分析方法。
实施例4：
本实施例用于提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现实施例1所述直升机桨叶结构降阶分析方法的步骤。
本说明书中每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种直升机桨叶结构降阶分析方法，其特征在于，所述结构降阶分析方法包括：
建立直升机桨叶三维单胞的有限元网格模型；所述直升机桨叶由多个所述三维单胞组装得到；
对于用于三维单胞变形仿真的多载荷工况中的每一载荷工况，基于所述载荷工况的工况载荷列向量计算有限元载荷列向量；以所述有限元载荷列向量为输入，利用含约束的有限元矩阵表达式计算有限元位移列向量；所述工况载荷列向量包括施加在所述有限元网格模型的第一虚拟节点和第二虚拟节点的载荷；所述有限元载荷列向量包括施加在所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的载荷；所述有限元位移列向量包括所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的位移；
基于所有所述载荷工况的有限元位移列向量计算所述三维单胞的等效属性矩阵；
建立所述直升机桨叶的简化模型；将所述直升机桨叶所受的实际载荷施加在所述简化模型上，基于有限元仿真计算得到所述简化模型的各个单元的广义单元应变；对于每一所述单元，基于所述单元的广义单元应变和所述等效属性矩阵计算所述单元的广义单元应力；所述简化模型包括多个节点和位于相邻两个所述节点之间的单元，每一所述单元代表一所述三维单胞；
对于每一所述单元，基于所述单元的广义单元应力计算所述单元的应力场。
2.根据权利要求1所述的直升机桨叶结构降阶分析方法，其特征在于，所述载荷工况包括基础工况和由所述基础工况组合得到的组合工况；所述基础工况包括拉伸、弯曲、扭转和剪切；
所有所述载荷工况的工况载荷列向量组成的载荷矩阵为：
Figure QLYQS_1
；
所述载荷矩阵的第t列即为第t个载荷工况的工况载荷列向量。
3.根据权利要求1所述的直升机桨叶结构降阶分析方法，其特征在于，所述有限元载荷列向量的计算公式包括：
当
Figure QLYQS_2
时，/>
Figure QLYQS_3
；
其中，t代表第t个载荷工况；
Figure QLYQS_4
为施加在有限元网格模型前端面的第i个网格节点的j方向上的载荷；i=1，2，...，G+，G+为端面的网格节点的数量；j=x，y，z；/>
Figure QLYQS_5
为施加在有限元网格模型后端面的第i个网格节点的j方向上的载荷；
当
Figure QLYQS_6
时，
Figure QLYQS_7
；
Figure QLYQS_8
；
其中，
Figure QLYQS_10
为第4个载荷工况下，有限元网格模型前端面的第i个网格节点的j方向上的等效节点力；/>
Figure QLYQS_12
为第一载荷修正系数；/>
Figure QLYQS_14
为第3个载荷工况下，有限元网格模型前端面的第i个网格节点的j方向上的等效节点力；/>
Figure QLYQS_11
为第二载荷修正系数；/>
Figure QLYQS_13
为第t个载荷工况的工况载荷列向量中第2行的元素值；/>
Figure QLYQS_15
为第t个载荷工况的工况载荷列向量中第3行的元素值；/>
Figure QLYQS_16
为第4个载荷工况下，有限元网格模型后端面的第i个网格节点的j方向上的等效节点力；/>
Figure QLYQS_9
为第3个载荷工况下，有限元网格模型后端面的第i个网格节点的j方向上的等效节点力。
4.根据权利要求1所述的直升机桨叶结构降阶分析方法，其特征在于，所述含约束的有限元矩阵表达式为：
Figure QLYQS_17
；
其中，
Figure QLYQS_18
为三维单胞的增广后有限元刚度矩阵；/>
Figure QLYQS_19
为罚函数因子；/>
Figure QLYQS_20
为对约束方程的系数矩阵进行数值积分所得到的计算矩阵；/>
Figure QLYQS_21
为有限元位移列向量；/>
Figure QLYQS_22
为有限元载荷列向量；
Figure QLYQS_23
为根据约束方程的非齐次项得到的修正载荷列向量；
所述约束方程为：
Figure QLYQS_24
；
Figure QLYQS_25
；
Figure QLYQS_26
；
其中，
Figure QLYQS_33
为有限元网格模型前端面的第i个网格节点在x方向的位移自由度；/>
Figure QLYQS_28
为有限元网格模型后端面的第i个网格节点在x方向的位移自由度；/>
Figure QLYQS_29
为有限元网格模型的第一虚拟节点在x方向的位移自由度；/>
Figure QLYQS_30
为第i个网格节点的y坐标；/>
Figure QLYQS_37
为参考轴线的y坐标；/>
Figure QLYQS_36
为有限元网格模型的第二虚拟节点在z方向的位移自由度；/>
Figure QLYQS_41
为第i个网格节点的z坐标；/>
Figure QLYQS_38
为参考轴线的z坐标；/>
Figure QLYQS_40
为有限元网格模型的第二虚拟节点在y方向的位移自由度；/>
Figure QLYQS_31
为有限元网格模型的第i个第三虚拟节点在x方向的位移自由度；/>
Figure QLYQS_34
为有限元网格模型前端面的第i个网格节点在y方向的位移自由度；/>
Figure QLYQS_42
为有限元网格模型后端面的第i个网格节点在y方向的位移自由度；/>
Figure QLYQS_44
为有限元网格模型的第一虚拟节点在y方向的位移自由度；/>
Figure QLYQS_43
为有限元网格模型的第二虚拟节点在x方向的位移自由度；/>
Figure QLYQS_45
为有限元网格模型的第i个第三虚拟节点在y方向的位移自由度；/>
Figure QLYQS_27
为有限元网格模型前端面的第i个网格节点在z方向的位移自由度；/>
Figure QLYQS_32
为有限元网格模型后端面的第i个网格节点在z方向的位移自由度；/>
Figure QLYQS_35
为有限元网格模型的第一虚拟节点在z方向的位移自由度；/>
Figure QLYQS_39
为有限元网格模型的第i个第三虚拟节点在z方向的位移自由度。
5.根据权利要求1所述的直升机桨叶结构降阶分析方法，其特征在于，所述等效属性矩阵的计算公式为：
Figure QLYQS_46
；
其中，
Figure QLYQS_47
为等效属性矩阵；/>
Figure QLYQS_48
为分块矩阵；
Figure QLYQS_49
；
其中，
Figure QLYQS_50
为第一分块；/>
Figure QLYQS_51
为第二分块；/>
Figure QLYQS_52
为第三分块；/>
Figure QLYQS_53
基于所有载荷工况的有限元位移列向量计算得到；/>
Figure QLYQS_54
基于所有载荷工况下的三维单胞的总应变能计算得到，计算公式为：
Figure QLYQS_55
；
其中，
Figure QLYQS_57
为第t个载荷工况下的三维单胞的总应变能，根据第t个载荷工况的有限元位移列向量计算得到；/>
Figure QLYQS_59
为第t个载荷工况的工况载荷列向量；/>
Figure QLYQS_60
为第t个载荷工况的工况载荷列向量中第2行的元素值；/>
Figure QLYQS_58
为/>
Figure QLYQS_61
中第4行第4列的元素值；/>
Figure QLYQS_62
为第t个载荷工况的工况载荷列向量中第3行的元素值；/>
Figure QLYQS_63
为/>
Figure QLYQS_56
中第3行第3列的元素值。
6.根据权利要求1所述的直升机桨叶结构降阶分析方法，其特征在于，所述基于所述单元的广义单元应变和所述等效属性矩阵计算所述单元的广义单元应力具体包括：计算所述单元的广义单元应变和所述等效属性矩阵的乘积，得到所述单元的广义单元应力。
7.根据权利要求1所述的直升机桨叶结构降阶分析方法，其特征在于，所述基于所述单元的广义单元应力计算所述单元的应力场具体包括：基于所述单元的广义单元应力计算所述单元的位移场，并根据所述位移场计算所述单元的应力场；
其中，所述位移场的计算公式为：
Figure QLYQS_64
；
其中，
Figure QLYQS_65
为第i个单元的位移场；/>
Figure QLYQS_66
为第i个单元的广义单元应力；/>
Figure QLYQS_67
为预设矩阵；/>
Figure QLYQS_68
为第t个载荷工况的工况载荷列向量；/>
Figure QLYQS_69
为由第t个载荷工况的有限元位移列向量的部分元素组成的位移列向量。
8.一种直升机桨叶结构降阶分析系统，其特征在于，所述结构降阶分析系统包括：
单胞模型构建模块，用于建立直升机桨叶三维单胞的有限元网格模型；所述直升机桨叶由多个所述三维单胞组装得到；
多工况变形仿真模块，用于对于用于三维单胞变形仿真的多载荷工况中的每一载荷工况，基于所述载荷工况的工况载荷列向量计算有限元载荷列向量；以所述有限元载荷列向量为输入，利用含约束的有限元矩阵表达式计算有限元位移列向量；所述工况载荷列向量包括施加在所述有限元网格模型的第一虚拟节点和第二虚拟节点的载荷；所述有限元载荷列向量包括施加在所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的载荷；所述有限元位移列向量包括所述有限元网格模型的各个网格节点、所述第一虚拟节点和所述第二虚拟节点的位移；
等效属性矩阵计算模块，用于基于所有所述载荷工况的有限元位移列向量计算所述三维单胞的等效属性矩阵；
桨叶模型构建模块，用于建立所述直升机桨叶的简化模型；将所述直升机桨叶所受的实际载荷施加在所述简化模型上，基于有限元仿真计算得到所述简化模型的各个单元的广义单元应变；对于每一所述单元，基于所述单元的广义单元应变和所述等效属性矩阵计算所述单元的广义单元应力；所述简化模型包括多个节点和位于相邻两个所述节点之间的单元，每一所述单元代表一所述三维单胞；
应力场计算模块，用于对于每一所述单元，基于所述单元的广义单元应力计算所述单元的应力场。
9.一种直升机桨叶结构降阶分析设备，其特征在于，包括：
处理器；以及
存储器，其中存储计算机可读程序指令，
其中，在所述计算机可读程序指令被所述处理器运行时执行如权利要求1-7任一项所述的直升机桨叶结构降阶分析方法。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-7任一项所述直升机桨叶结构降阶分析方法的步骤。
CN202310571708.1A 2023-05-22 2023-05-22 一种直升机桨叶结构降阶分析方法、系统、设备及介质 Active CN116305589B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310571708.1A CN116305589B (zh) 2023-05-22 2023-05-22 一种直升机桨叶结构降阶分析方法、系统、设备及介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310571708.1A CN116305589B (zh) 2023-05-22 2023-05-22 一种直升机桨叶结构降阶分析方法、系统、设备及介质

Publications (2)

Publication Number Publication Date
CN116305589A true CN116305589A (zh) 2023-06-23
CN116305589B CN116305589B (zh) 2023-08-11

Family

ID=86824335

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310571708.1A Active CN116305589B (zh) 2023-05-22 2023-05-22 一种直升机桨叶结构降阶分析方法、系统、设备及介质

Country Status (1)

Country Link
CN (1) CN116305589B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117540494A (zh) * 2024-01-10 2024-02-09 天目山实验室 一种航空舱段结构降阶分析方法、系统及设备

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109933837A (zh) * 2019-01-04 2019-06-25 中国电力科学研究院有限公司 一种风电机组强度快速评估方法和系统
US20200159879A1 (en) * 2018-04-24 2020-05-21 Blade Diagnostices Corporation Refinement of finite element model of integrally bladed disk
CN111507042A (zh) * 2020-04-29 2020-08-07 西安交通大学 基于叶端定时的旋转叶片动应力测量方法及其系统
CN112560167A (zh) * 2020-11-10 2021-03-26 北京航空航天大学 机翼结构力学高保真降阶仿真方法、电子设备、存储介质
CN113486552A (zh) * 2021-07-01 2021-10-08 江苏科技大学 一种单向复合材料弯曲刚度预测方法
CN114065394A (zh) * 2021-11-19 2022-02-18 中国直升机设计研究所 一种直升机机体主承力结构应力分析方法
CN114154379A (zh) * 2021-12-06 2022-03-08 中国船舶科学研究中心 一种极地船舶螺旋桨强度计算方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20200159879A1 (en) * 2018-04-24 2020-05-21 Blade Diagnostices Corporation Refinement of finite element model of integrally bladed disk
CN109933837A (zh) * 2019-01-04 2019-06-25 中国电力科学研究院有限公司 一种风电机组强度快速评估方法和系统
CN111507042A (zh) * 2020-04-29 2020-08-07 西安交通大学 基于叶端定时的旋转叶片动应力测量方法及其系统
CN112560167A (zh) * 2020-11-10 2021-03-26 北京航空航天大学 机翼结构力学高保真降阶仿真方法、电子设备、存储介质
CN113486552A (zh) * 2021-07-01 2021-10-08 江苏科技大学 一种单向复合材料弯曲刚度预测方法
CN114065394A (zh) * 2021-11-19 2022-02-18 中国直升机设计研究所 一种直升机机体主承力结构应力分析方法
CN114154379A (zh) * 2021-12-06 2022-03-08 中国船舶科学研究中心 一种极地船舶螺旋桨强度计算方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
刘涛;朱旭;杨成;孙会伟;: "水平轴风力机叶片的结构设计与有限元分析", 机械制造, no. 09, pages 8 - 10 *
周丞;林杰;刘勇;: "基于拓扑优化和形状优化的桨叶结构设计", 航空工程进展, no. 01, pages 113 - 119 *
许燕;孙文磊;黄艳华;: "基于逆向工程的风力机叶片动态特性分析", 机械设计与制造, no. 04, pages 214 - 216 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117540494A (zh) * 2024-01-10 2024-02-09 天目山实验室 一种航空舱段结构降阶分析方法、系统及设备
CN117540494B (zh) * 2024-01-10 2024-04-05 天目山实验室 一种航空舱段结构降阶分析方法、系统及设备

Also Published As

Publication number Publication date
CN116305589B (zh) 2023-08-11

Similar Documents

Publication Publication Date Title
Chimakurthi et al. Computational aeroelasticity framework for analyzing flapping wing micro air vehicles
CN106055764A (zh) 基于三维壳有限元‑梁模型的风力机叶片位移计算方法
CN112560167B (zh) 机翼结构力学高保真降阶仿真方法、电子设备、存储介质
Mian et al. Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach
CN116305589B (zh) 一种直升机桨叶结构降阶分析方法、系统、设备及介质
Nguyen et al. Modeling and analysis of insect-like flexible wings at low Reynolds number
CN112016167A (zh) 基于仿真和优化耦合的飞行器气动外形设计方法及系统
CN107391891A (zh) 一种基于模型融合方法的大展弦比机翼优化设计方法
CN113886967B (zh) 多巡航工况的大型飞机机翼气动弹性优化方法
Tang et al. Study on the aeroelastic responses of a wind turbine using a coupled multibody-FVW method
CN112580241A (zh) 一种基于结构降阶模型的非线性气动弹性动响应分析方法
Liu et al. An integrated high-fidelity approach for modeling flow-structure interaction in biological propulsion and its strong validation
Charvet et al. Numerical simulation of the flow over sails in real sailing conditions
CN117436322B (zh) 基于叶素理论的风力机叶片气动弹性仿真方法和介质
Díaz-Casás et al. Wind turbine design through evolutionary algorithms based on surrogate CFD methods
CN110147571B (zh) 一种组件结构的拓扑优化方法及装置
Aono et al. A numerical framework for isotropic and anisotropic flexible flapping wing aerodynamics and aeroelasticity
CN114139285B (zh) 一种柔性蒙皮机翼低雷诺数流固耦合分析方法
CN113486512B (zh) 一种功能梯度变厚度叶片模型的颤振分析方法
CN112214843B (zh) 一种用于风洞试验机翼模型有限元刚度修正方法及装置
US11983467B2 (en) Rapid aero modeling for computational experiments (RAM-C)
Lankford et al. Coupled aeroelastic study of a flexible micro air vehicle-scale flapping wing in hovering flight
Werthen et al. Comparison of different cross-sectional approaches for the structural design and optimization of composite wind turbine blades based on beam models
Ahmad et al. High-order accurate CFD/CSD simulation of the UH60 rotor in forward flight
Malhan et al. Cfd-csd coupled aeroelastic analysis of flexible flapping wings for mav applications: methodology validation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant