CN108009370B - 一种结构应力敏度求解方法 - Google Patents

一种结构应力敏度求解方法 Download PDF

Info

Publication number
CN108009370B
CN108009370B CN201711330533.6A CN201711330533A CN108009370B CN 108009370 B CN108009370 B CN 108009370B CN 201711330533 A CN201711330533 A CN 201711330533A CN 108009370 B CN108009370 B CN 108009370B
Authority
CN
China
Prior art keywords
matrix
stress
unit
sensitivity
displacement
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.)
Active
Application number
CN201711330533.6A
Other languages
English (en)
Other versions
CN108009370A (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.)
AVIC Aircraft Strength Research Institute
Original Assignee
AVIC Aircraft Strength Research Institute
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 AVIC Aircraft Strength Research Institute filed Critical AVIC Aircraft Strength Research Institute
Priority to CN201711330533.6A priority Critical patent/CN108009370B/zh
Publication of CN108009370A publication Critical patent/CN108009370A/zh
Application granted granted Critical
Publication of CN108009370B publication Critical patent/CN108009370B/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/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
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/17—Mechanical parametric or variational design

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)

Abstract

本发明公开了一种结构应力敏度求解方法。所述结构应力敏度求解方法包括如下步骤：步骤1：计算结构有限元分析，从结构位移列阵中挑出当前单元节点对应的位移值；步骤2：进行总体刚度矩阵的广义逆阵计算；步骤3：建立设计变量，对金属杆板梁取其面积或厚度，对复合材料层合板取各个材料相同铺层角相同的单层层组合为单一设计变量；步骤4：挑选出当前作为约束的应力分量，根据当前单元类型和所求敏度的应力分量来确定构建挑选矩阵；步骤5：计算当前单元对应的微分刚度矩阵；步骤6：通过广义逆矩阵，应力挑选矩阵，应力矩阵，总体坐标到单元坐标的转换矩阵，位移选择矩阵，单元对设计的微分刚度矩阵及位移列阵的一系列矩阵连乘运算得到当前单元应力对设计变量的敏度计算；步骤7：按单元循环以上步骤4至所述步骤6完成所有单元应力约束的敏度求解计算。

Description

一种结构应力敏度求解方法
技术领域
本发明涉及结构应力分析技术领域，特别是涉及一种结构应力敏度求解方法。
背景技术
优化技术是满足日益提高的飞机综合性能要求及由继承设计、经验设计到创新设计的必由之路。
随着材料结构一体化等技术的发展，结构中可设计的变量规模日益膨胀，以十万变量以上的优化模型为基础开展精细优化设计已经提到日程，为对结构开展精细设计必须将应力作为约束进行考虑，由于工艺制备技术的发展使得按单元进行应力控制的设计结果可以不经过传统的圆整过程而直接进行生产制造，从程序的流程组织考虑利用刚度矩阵求逆对大变量的应力敏度求解是最为高效的程序实现方式，由于刚度矩阵通常是奇异的，按照常规的矩阵求逆算法将遇到数值求解困难使得对应力敏度的求解过程难以进行。
因此，希望有一种技术方案来克服或至少减轻现有技术的至少一个上述缺陷。
发明内容
本发明的目的在于提供一种结构应力敏度求解方法来克服或至少减轻现有技术的至少一个上述缺陷。
为实现上述目的，本发明提供一种结构应力敏度求解方法，所述结构应力敏度求解方法包括如下步骤：
步骤1：计算结构有限元分析，从结构位移列阵中挑出当前单元节点对应的位移值；
步骤2：进行总体刚度矩阵的广义逆阵计算；
步骤3：建立设计变量，对金属杆板梁取其面积或厚度，对复合材料层合板取各个材料相同铺层角相同的单层层组合为单一设计变量；
步骤4：挑选出当前作为约束的应力分量，根据当前单元类型和所求敏度的应力分量来确定构建挑选矩阵；
步骤5：计算当前单元对应的微分刚度矩阵；
步骤6：通过总体刚度矩阵的广义逆矩阵，应力挑选矩阵，应力矩阵，总体坐标到单元坐标的转换矩阵，位移选择矩阵，单元对设计的微分刚度矩阵及位移列阵的一系列矩阵连乘运算得到当前单元应力对设计变量的敏度计算；
步骤7：参加优化的设计变量所在的单元循环所述步骤4至所述步骤6完成所有单元应力约束的敏度求解计算。
本申请的结构应力敏度求解方法通过将总体刚度矩阵求逆运算从总体计算循环中分离出来并利用构造伪逆矩阵或moore-Penrose逆矩阵，利用A+的唯一确定性来解决奇异矩阵求逆问题，从而简化大规模应力约束对敏度计算的流程组织方式并极大的提高敏度求解效率。避免了反复进行矩阵聚缩运算的传统流程组织方式，以该技术方法结合有限元求解技术对敏度分析流程进行设计后，将使得应力敏度求解效率得到大幅提升。
附图说明
图1是本申请第一实施例的结构应力敏度求解方法的流程示意图。
图2是复合材料板壳元的分层示意图。
图3是板壳单元构型示意图。
图4以梁单元为示例的单元刚度矩阵示意图。
具体实施方式
为使本发明实施的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行更加详细的描述。在附图中，自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。所描述的实施例是本发明一部分实施例，而不是全部的实施例。下面通过参考附图描述的实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。下面结合附图对本发明的实施例进行详细说明。
在本发明的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明保护范围的限制。
图1是本申请第一实施例的结构应力敏度求解方法的流程示意图。图2是复合材料板壳元的分层示意图。图3是板壳单元构型示意图。图4以梁单元为示例的单元刚度矩阵示意图。
如图1所示的结构应力敏度求解方法包括如下步骤：
步骤1：获取结构有限元信息，计算结构位移响应，从结构位移列阵中挑出当前单元节点对应的位移值；在本实施例中，当前单元节点指的是要进行应力敏度分析的单元，其节点选取按照单元类型不同而不同，对杆、梁单元取两个连接点，对三角形、四边形单元取单元角点。
步骤2：进行总体刚度矩阵的广义逆阵计算；
步骤3：建立设计变量，对金属杆板梁取其面积或厚度，对复合材料层合板取各个材料相同铺层角相同的单层层组合为单一设计变量；对航空结构其有限元模型均以薄壁单元进行模拟，因此本方案的设计变量限定为杆、梁截面面积，板壳厚度，对复合材料按照层组设计。
步骤4：挑选出当前作为约束的应力分量，根据当前单元类型和所求敏度的应力分量来确定构建挑选矩阵；具体地，工程上每一个单元有六个应力分量，在开展应力约束优化前，设计员必须明确指定需要控制的各个结构部位的应力大小及方向，即应力约束是拉应力、压应力还是剪应力等，约束期望控制在多大水平，敏度求解时对每个单元仅仅挑选用户期望的应力分量进行分析。
步骤5：计算当前单元对应的微分刚度矩阵；
步骤6：通过广义逆矩阵，应力挑选矩阵，应力矩阵，总体坐标到单元坐标的转换矩阵，位移选择矩阵，单元对设计的微分刚度矩阵及位移列阵的一系列矩阵连乘运算得到当前单元应力对设计变量的敏度计算；
步骤7：按单元循环以上步骤4至所述步骤6完成所有单元应力约束的敏度求解计算。
本申请的结构应力敏度求解方法通过将总体刚度矩阵求逆运算从总体计算循环中分离出来并利用构造伪逆矩阵或moore-Penrose逆矩阵，利用A+的唯一确定性来解决奇异矩阵求逆问题，从而简化大规模应力约束对敏度计算的流程组织方式并极大的提高敏度求解效率。避免了反复进行矩阵聚缩运算的传统流程组织方式，以该技术方法结合有限元求解技术对敏度分析流程进行设计后，将使得应力敏度求解效率得到大幅提升。
下面以举例的方式对本申请进行进一步阐述。可以理解的是，该举例并不构成对本申请的任何限制。
步骤1：进行结构有限元分析，获得结构位移响应，从结构位移列阵中挑出当前单元节点对应的位移值[ue]n×1＝[T]n×N·[u]N×1；
步骤2：进行总体刚度矩阵的广义逆阵计算；称满足式下式(1)的矩阵G为A的伪逆矩阵，用符号A+表示。
AGA＝A
GAG＝G
(AG)T＝AG
(GA)T＝GA (1)
利用奇异值分解算法(SVD)，求得A+，对结构有限元分析来可以利用这一过程求解总体刚度矩阵的K的广义逆；
步骤3：建立设计变量(图1)，对金属杆板梁取其面积或厚度，对复合材料层合板取各个材料相同铺层角相同的单层层组合为单一设计变量(图2)；
步骤4：挑选出当前作为约束的应力分量，根据当前单元类型(图3为板壳单元拓扑形式)和所求敏度的应力分量来确定构建挑选矩阵；
步骤5：计算当前单元对应的微分刚度矩阵，(图4为梁单元的刚度矩阵形式)。对复合材料板壳元的优化设计，采用的是分层设计思想，即把材料相同和角度相同的单层集合定义为一个分层，分层的厚度等于各单层厚度相加，即
Figure BDA0001506541890000051
那么就需要把经典层合板理论中用单层坐标表示的刚度矩阵转换成用分层厚度Ti显式表达的形式，这样就比较容易得到各刚度矩阵相应的关于分层厚度的微分刚阵。拉伸刚度和剪切刚度是分层厚度Tk的显式表达式，因此其微分形式可通过直接对分层厚度Tk求导得到。故：
层合板刚度矩阵关于分层厚度的导数精确解为：
Figure BDA0001506541890000052
Figure BDA0001506541890000053
Figure BDA0001506541890000054
Figure BDA0001506541890000055
步骤6：通过广义逆矩阵，应力挑选矩阵，应力矩阵，总体坐标到单元坐标的转换矩阵，位移选择矩阵，单元对设计的微分刚度矩阵及位移列阵的一系列矩阵连乘运算得到当前单元应力对设计变量的敏度计算：
Figure BDA0001506541890000056
其中：
[K-1]N×N将在敏度正式计算之前利用步骤(2)直接求逆得到
[Tg]1×m根据当前单元类型和所求敏度的应力分量来确定
Figure BDA0001506541890000061
在开发位移敏度求解模块的时候已经解决
[u]N×1通过静力分析得到
[S]m×n应力矩阵
步骤7：重复步骤4到6直到对所用参与应力约束的单元完成计算计。
最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (1)

1.一种结构应力敏度求解方法，其特征在于，所述结构应力敏度求解方法包括如下步骤：
步骤1：计算结构有限元分析，从结构位移列阵中挑出当前单元节点对应的位移值；
步骤2：进行总体刚度矩阵的广义逆阵计算；
具体包括：称满足式下式的矩阵G为A的伪逆矩阵，用符号A+表示：
AGA＝A
GAG＝G
(AG)T＝AG
(GA)T＝GA
利用奇异值分解算法(SVD)，求得A+，对结构有限元分析利用这一过程求解总体刚度矩阵的K的广义逆；
步骤3：建立设计变量，对金属杆板梁取其面积或厚度，对复合材料层合板取各个材料相同铺层角相同的单层组合为单一设计变量；
步骤4：挑选出当前作为约束的应力分量，根据当前单元类型和所求敏度的应力分量来确定构建挑选矩阵；
步骤5：计算当前单元对应的微分刚度矩阵，采用分层设计思想，分层的厚度等于各单层厚度相加，将单层坐标表示的刚度矩阵转换成用分层厚度显式表达的形式，得到各刚度矩阵相应的关于分层厚度的微分刚度矩阵；
步骤6：通过总体刚度矩阵的广义逆矩阵，应力挑选矩阵，应力矩阵，总体坐标到单元坐标的转换矩阵，位移选择矩阵，单元对设计的微分刚度矩阵及位移列阵的一系列矩阵连乘运算得到当前单元应力对设计变量的敏度计算；
具体为：
Figure FDA0003095478430000011
步骤7：参加优化的设计变量所在的单元循环所述步骤4至所述步骤6完成所有单元应力约束的敏度求解计算。
CN201711330533.6A 2017-12-13 2017-12-13 一种结构应力敏度求解方法 Active CN108009370B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201711330533.6A CN108009370B (zh) 2017-12-13 2017-12-13 一种结构应力敏度求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201711330533.6A CN108009370B (zh) 2017-12-13 2017-12-13 一种结构应力敏度求解方法

Publications (2)

Publication Number Publication Date
CN108009370A CN108009370A (zh) 2018-05-08
CN108009370B true CN108009370B (zh) 2021-08-17

Family

ID=62058676

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201711330533.6A Active CN108009370B (zh) 2017-12-13 2017-12-13 一种结构应力敏度求解方法

Country Status (1)

Country Link
CN (1) CN108009370B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110096761B (zh) * 2019-04-11 2020-01-31 河海大学 一种针对自由曲面层合壳的形状与铺层顺序同步优化方法
CN113704888B (zh) * 2021-08-23 2024-02-23 中国飞机强度研究所 一种单元应力筛选方法
CN116182730B (zh) * 2023-04-14 2023-07-18 交通运输部公路科学研究所 一种基于光纤光栅阵列传感的桥梁变形监测系统及方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105808884A (zh) * 2016-03-30 2016-07-27 北京航空航天大学 一种基于分形理论的有界不确定性平面裂纹应力强度因子上下界的预测方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CH695515A5 (de) * 2002-03-08 2006-06-15 Alstom Technology Ltd Verfahren zur Bestimmung des elasto-plastischen Verhaltens von aus anisotropem Material bestehenden Bauteilen sowie Anwendung des Verfahrens.
CN102902845B (zh) * 2012-09-12 2015-01-28 北京航空航天大学 一种直升机旋翼桨叶剖面设计方法
CN105181510B (zh) * 2015-10-10 2018-06-05 中国飞机强度研究所 一种梯度材料宏观等效弹性模量的计算方法
CN106202699B (zh) * 2016-07-07 2019-07-19 中国飞机强度研究所 一种多位移约束下的敏度求解方法
CN106844965B (zh) * 2017-01-19 2020-05-22 山西省交通建设工程质量检测中心（有限公司） 一种基于静载试验识别连续梁桥实际刚度的方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105808884A (zh) * 2016-03-30 2016-07-27 北京航空航天大学 一种基于分形理论的有界不确定性平面裂纹应力强度因子上下界的预测方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于分层设计变量的复合材料层合板优化设计及灵敏度分析;赵国忠等;《工程力学》;20030430;第20卷(第2期);第60-65页 *

Also Published As

Publication number Publication date
CN108009370A (zh) 2018-05-08

Similar Documents

Publication Publication Date Title
CN108009370B (zh) 一种结构应力敏度求解方法
US11003807B2 (en) Techniques for generating materials to satisfy design criteria
Cinefra et al. MITC9 shell elements based on refined theories for the analysis of isotropic cylindrical structures
Schmit et al. Finite deflection structural analysis using plate and shell discreteelements.
CN109858133B (zh) 一种基于应力映射的点阵结构设计与优化方法
Stanford et al. Aeroelastic optimization of flapping wing venation: A cellular division approach
Oktay et al. Three-dimensional structural topology optimization of aerial vehicles under aerodynamic loads
Feil et al. A cross-sectional aeroelastic analysis and structural optimization tool for slender composite structures
Jaworski et al. Implementation features of composite materials effective mechanical characteristics finding method based on microlevel cellular structural models
CN106683169A (zh) 一种关节运动感知的稀疏局部分解及重构算法
CN103065015B (zh) 一种基于内力路径几何形态的承载结构低碳节材设计方法
CN102663152B (zh) 一种异型蜂窝蒙皮结构的有限元建模方法
CN112231863B (zh) 太阳翼电池阵基板建模方法、装置、设备及存储介质
CN112395746B (zh) 微结构族等效材料性质的计算方法、微结构、系统及介质
CN110321571A (zh) 一种蜂窝板壳结构的力学参数数值提取方法
CN109460590B (zh) 一种新型含柔性翅脉节点的昆虫翅膀有限元建模方法
CN103425830B (zh) 基于多点位移协调约束的结构拓扑优化方法
Carrera et al. Static analysis of reinforced thin-walled plates and shells by means of finite element models
Osman-Letelier et al. Shape optimization of concrete shells with ruled surface geometry using line geometry
Öztorun A rectangular finite element formulation
Yan et al. Multiscale analysis of thermal stress of lattice materials and its size effects
Querin et al. Layout optimization of multi-material continuum structures with the isolines topology design method
CN113836768B (zh) 一种基于曲线/曲面桥接节点的异质结构高精度仿真方法
Al-Sabah et al. Rotation-free finite element ‘yield line’analysis of non-isotropic slabs
Prohasky et al. SmartNodes:‘Light-weight’parametric structural design process with BESO

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