CN114912226A - 考虑离心载荷和应力约束对结构进行优化设计的方法 - Google Patents

考虑离心载荷和应力约束对结构进行优化设计的方法 Download PDF

Info

Publication number
CN114912226A
CN114912226A CN202210655155.3A CN202210655155A CN114912226A CN 114912226 A CN114912226 A CN 114912226A CN 202210655155 A CN202210655155 A CN 202210655155A CN 114912226 A CN114912226 A CN 114912226A
Authority
CN
China
Prior art keywords
stress
iteration
design
cell
value
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
CN202210655155.3A
Other languages
English (en)
Other versions
CN114912226B (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.)
Xiamen University
Original Assignee
Xiamen 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 Xiamen University filed Critical Xiamen University
Priority to CN202210655155.3A priority Critical patent/CN114912226B/zh
Publication of CN114912226A publication Critical patent/CN114912226A/zh
Priority to US18/169,144 priority patent/US20230409767A1/en
Application granted granted Critical
Publication of CN114912226B publication Critical patent/CN114912226B/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/13—Architectural design, e.g. computer-aided architectural design [CAAD] related to design of buildings, bridges, landscapes, production plants or roads
    • 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
    • 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
    • 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
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/06—Multi-objective optimisation, e.g. Pareto optimisation using simulated annealing [SA], ant colony algorithms or genetic algorithms [GA]
    • 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)
  • Geometry (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Architecture (AREA)
  • Civil Engineering (AREA)
  • Structural Engineering (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Feedback Control In General (AREA)
  • Testing Of Devices, Machine Parts, Or Other Structures Thereof (AREA)

Abstract

本申请公开了考虑离心载荷和应力约束对结构进行优化设计的方法，相比于现有技术，本申请的技术方案对放松系数c获得方式进行了改进，包括通过步骤S9.1至S9.8基于预测应力法对第二预测最大应力的计算，这减轻了优化过程中锯齿边界和模糊区域的存在对结构应力场计算的影响。其中最主要的是根据设计变量小于0.1和大于0.9的单元数量与单元总数量的比值，决定对各单元的弹性模量采用线性惩罚还是非线性惩罚。引入预测最大应力。能够更充分地利用材料许用应力，提高优化设计质量，获得较现有技术质量更轻的设计方案。

Description

考虑离心载荷和应力约束对结构进行优化设计的方法
技术领域
本申请涉及结构设计领域，具体涉及考虑离心载荷和应力约束对结构进行优化设计的方法。
背景技术
涡轮盘是航空发动机重要的承力部件。涡轮盘的重量直接关系到结构的效率，更轻的涡轮盘可以提高航空发动机的推重比。因此，需要对涡轮盘进行轻量化设计。
目前传统的单辐板涡轮盘的形状优化设计已经非常成熟，很难再提升涡轮盘结构的性能。拓扑优化作为一种强有力概念设计工具，广泛应用于工业制造领域。使用拓扑优化技术对涡轮盘进行优化，能够得到更加轻量化的新颖结构。增材制造技术的快速发展，为制造复杂构型的设计提供了技术支撑。拓扑优化能够在给定的约束条件下，在特定的区域内找到最优的材料分布，使给定的目标函数值达到最大值或最小值。
现有技术中，基于变密度法的拓扑优化已经应用于不复杂工业产品的设计。但在应用于受离心载荷和应力约束的涡轮盘的设计过程中，却经常出现优化设计出结构中最大应力值小于材料许用应力值70％的情况，即仍有很大的减轻涡轮盘质量的余量。
发明内容
本申请的目的在于克服背景技术中存在的上述缺陷或问题，提供一种考虑离心载荷和应力约束对结构进行优化设计的方法，其能够更充分地利用材料许用应力，提高优化设计质量，获得较现有技术质量更轻的设计方案。
为达成上述目的，采用如下技术方案：
考虑离心载荷和应力约束对结构进行优化设计的方法，其优化目标为在应力约束条件下结构质量最小，所述方法包括：对设计域内的结构进行离散的步骤、初始化离散后的各单元的设计变量xe形成初始设计变量集的步骤、利用变密度法进行迭代设计直至得到最优设计变量集的步骤和根据最优设计变量集对结构进行光滑化的步骤，其中，所有设计变量集中的所有单元的设计变量xe应满足0≤xe≤1；其特征是，每次迭代设计执行如下步骤：
S1：基于初始设计变量集或上次迭代形成的设计变量集，用密度过滤法对各单元的设计变量进行过滤，得到单元过滤密度
Figure BDA0003689171460000011
S2：用Heaviside函数对各单元的过滤密度
Figure BDA0003689171460000021
进行映射，得到第一密度ρe1：
Figure BDA0003689171460000022
其中，
Figure BDA0003689171460000023
用于控制映射曲线的光滑度，k为本次迭代的序号，k＝1时，
Figure BDA0003689171460000024
从第1次迭代起，每20次迭代更新一次，每次更新的
Figure BDA0003689171460000025
值为上一
Figure BDA0003689171460000026
值的1.19倍，但
Figure BDA0003689171460000027
值最大为8；
S3：用RAMP法计算各单元的第一弹性模量Ee1：
Figure BDA0003689171460000028
其中，Emax是设计材料的弹性模量值，Emin是为了避免奇异矩阵添加的小量；q1为第一惩罚因子，q1＝4；
S4：组装全局刚度矩阵，根据静力平衡方程计算结构位移；
S5：计算各单元的第一应力σe1；
S6：计算各单元的第一米塞斯应力σVM1；
S7：用RAMP法对各单元的第一米塞斯应力σVM1进行惩罚，得到各单元的第一惩罚应力
Figure BDA0003689171460000029
Figure BDA00036891714600000210
其中，q2为第二惩罚因子，q2＝-0.95；
S8：聚合各单元的第一惩罚应力
Figure BDA00036891714600000211
得到聚合应力
Figure BDA00036891714600000212
S9：对应力约束进行放松：
将应力约束放松为
Figure BDA00036891714600000213
其中，σl为材料的许用应力，c为本次迭代的放松系数，c的值每5次迭代更新一次，每次更新时
Figure BDA00036891714600000214
其中，
Figure BDA00036891714600000215
为本次迭代的第一最大预测应力，
Figure BDA00036891714600000216
其中，
Figure BDA00036891714600000217
为本次迭代的第二最大预测应力，第二最大预测应力
Figure BDA00036891714600000218
在本次迭代中通过如下步骤获得：
S9.1：用Heaviside函数对每个单元的过滤密度
Figure BDA00036891714600000219
进行映射，得到第二密度ρe2：
Figure BDA0003689171460000031
其中，
Figure BDA0003689171460000032
用于控制映射曲线的光滑度，＝1时，
Figure BDA0003689171460000033
从第1次迭代起，每20次迭代更新一次，每次更新的
Figure BDA0003689171460000034
值为上次
Figure BDA0003689171460000035
值加4，但
Figure BDA0003689171460000036
值最大为20；
S9.2：计算转变因子
Figure BDA0003689171460000037
Figure BDA0003689171460000038
其中，N为所有单元的总数量，Nb为第二密度ρe2大于第一阈值的单元的数量，Nw为第二密度ρe2小于第二阈值的单元的数量；所述第一阈值大于0.75，所述第二阈值小于0.25；
S9.3：计算各单元的第二弹性模量Ee2：
Figure BDA0003689171460000039
S9.4：组装全局刚度矩阵，根据静力平衡方程计算结构位移；
S9.5：计算各单元的第二应力σe2；
S9.6：计算各单元的第二米塞斯应力σVM2；
S9.7：用线性法对各单元的第二米塞斯应力σVM2进行惩罚，得到各单元的第二惩罚应力
Figure BDA00036891714600000310
Figure BDA00036891714600000311
S9.8：计算本次迭代的第二最大预测应力
Figure BDA00036891714600000312
Figure BDA00036891714600000313
S10：对目标函数Vf进行灵敏度分析；
其中，目标函数
Figure BDA00036891714600000314
其中，ve为第e个单元的体积；
S11：采用移动渐近线算法对优化问题进行求解计算，更新各单元设计变量，得到并储存本次迭代形成的设计变量集；
S12：判断本次迭代是否满足退出迭代条件，本次迭代如满足退出迭代条件则退出迭代，并将本次迭代形成的设计变量集记为最优设计变量集；如不满足退出迭代条件则进行下一次迭代；
以上步骤S9.1至S9.8的过程允许与步骤S2至S8的过程并行计算。
优选地，根据最优设计变量集对结构进行光滑化的步骤采用如下方法：在最优设计变量集中，对于任一单元，如果xe<0.5，则该单元所在空间无材料，如果xe>0.5，则该单元所在空间具有全部材料，如果xe＝0.5，则该单元位于交界面。
优选地，步骤S9.2中，第一阈值为0.9，第二阈值为0.1。
优选地，步骤S8中，采用P-norm法聚合各单元的第一惩罚应力
Figure BDA0003689171460000041
优选地，步骤S12中，判断本次迭代是否满足退出迭代条件的方法是对比本次迭代形成的设计变量集和上次迭代形成的设计变量集，或者对比本次迭代得到的目标函数值和上次迭代得到的目标函数值，并且第二最大预测应力小于等于材料许用应力。
优选地，退出迭代条件为本次迭代形成的设计变量集中的所有设计变量与上次迭代形成的设计变量集中对应的设计变量的差值的绝对值均小于第三阈值。
优选地，所述第三阈值为0.05。
相对于现有技术，上述方案具有的如下有益效果：
申请人发现，现有技术在将变密度法应用于离心载荷和应力约束下的结构分析时，由于在最后一步中需要在最优变量集基础上对结构进行光滑化，而选代优化形成的最优变量集由于存在锯齿边界和在边界处存在模糊区域，光滑化的过程使优化出的结果往往远离原来设定的应力约束条件。同时，离心载荷的存在，采用SIMP法在步骤S7中惩罚米塞斯应力(Von Mises应力)，会导致低密度单元具有较低的灵敏度，不利于拓扑优化。因此，现有技术中用变密度法应用于离心载荷和应力约束下的结构分析效果不能令人满意。
本发明相对于现有技术，主要的改进及其技术效果在于：
1、步骤S9中，对放松系数c获得方式的改进，包括通过步骤S9.1至S9.8基于预测应力法对
Figure BDA0003689171460000042
的计算。减轻了优化过程中锯齿边界和模糊区域的存在对结构应力场计算的影响。其中最主要的是根据设计变量小于0.1和大于0.9的单元数量与单元总数量的比值，决定对各单元的弹性模量采用线性惩罚还是非线性惩罚。引入最大预测应力，能够更充分地利用材料许用应力，提高优化设计质量，获得较现有技术质量更轻的设计方案。
2、步骤S7中，用RAMP法对各单元的的米塞斯应力进行非线性惩罚，且惩罚因子选用-0.95，能够在离心载荷和应力约束下很好地提高低密度单元的灵敏度，有利于拓扑优化。
附图说明
为了更清楚地说明实施例的技术方案，下面简要介绍所需要使用的附图：
图1为本发明的优化设计的方法流程图；
图2为典型涡轮盘结构示意图；
图3为实施例中涡轮盘设计域示意图及载荷条件；
图4为实施例中涡轮盘根据最优设计变量集构造的涡轮盘结构示意图；
图5为实施例中涡轮盘光滑化结构后的结构示意图；
图6为验证例的结构图。
具体实施方式
权利要求书和说明书中，除非另有限定，术语“第一”、“第二”或“第三”等，都是为了区别不同对象，而不是用于描述特定顺序。
权利要求书和说明书中，除非另有限定，术语“中心”、“横向”、“纵向”、“水平”、“垂直”、“顶”、“底”、“内”、“外”、“上”、“下”、“前”、“后”、“左”、“右”、“顺时针”、“逆时针”等指示的方位或位置关系乃基于附图所示的方位和位置关系，且仅是为了便于简化描述，而不是暗示所指的装置或元件必须具有特定的方位或以特定的方位构造和操作。
权利要求书和说明书中，除非另有限定，术语“固接”或“固定连接”，应作广义理解，即两者之间没有位移关系和相对转动关系的任何连接方式，也就是说包括不可拆卸地固定连接、可拆卸地固定连接、连为一体以及通过其他装置或元件固定连接。
权利要求书和说明书中，除非另有限定，术语“包括”、“具有”以及它们的变形，意为“包含但不限于”。
权利要求书和说明书中，各种术语在相关文献中有不同的名称，应当根据本领域技术人员对其过程的理解结合说明书全文的内容确定其真实的含义。
下面将结合附图，对实施例中的技术方案进行清楚、完整地描述。
本实施例涉及对涡轮盘的优化设计。涡轮盘盘体结构如图2所示。航空发动机工作时，涡轮盘处于高速旋转的状态。涡轮盘所受的载荷主要来自两个方面：一方面,涡轮盘自身高速旋转所产生的离心载荷；另一方面，涡轮叶片高速旋转所产生的对涡轮盘的拉力。涡轮自身产生的离心负荷与中心旋转轴对称；叶片高速旋转对涡轮盘的拉力可以近似均匀分布在盘缘面上，拉力对称于中心旋转轴。涡轮盘所受到的约束只有轴向位移约束。因此，涡轮盘是轴对称模型，可以将原三维模型简化为二维轴对称模型，可以大幅度节约计算资源，提高优化设计效率。同时，考虑到拓扑优化是概念设计，可以对现有的设计域进行扩展。
图3是扩展后涡轮盘轴对称模型示意图及载荷条件。其中阴影部分为非设计区域，空白区域为初始设计区域。涡轮盘的几何形状由以下参数定义：R1＝83mm为涡轮盘内缘半径，R2＝237mm为涡轮盘外缘半径，H1＝92mm为涡轮盘内缘高度，H2＝40mm为涡轮盘外缘高度，W＝2mm为非设计域宽度。涡轮盘的材料为线弹性材料：杨氏模量192GPa，密度8240kg/m3，泊松比0.3，材料的许用应力为1000MPa。涡轮盘的旋转速度为10000r/min；叶片产生的等效载荷为100MPa，作用在F所在的区域；在1点处轴向位移为0。
参见图1，图1示出了本实施例基于变密度法对受离心载荷和应力约束的上述涡轮盘结构进行优化设计的方法的各步骤：
第一步骤，通过有限元法，采用四节点矩形轴对称单元对设计域内的结构进行离散，结构被离散为N个单元，每个单元对应一个设计变量xe，设计变量xe应满足0≤xe≤1，其中e为单元的序号。所有设计变量组合成一个向量x。考虑应力约束，以结构的质量分数最小作为目标函数，优化目标可表示为：
Figure BDA0003689171460000061
st:σe-σl≤0
0≤xe≤1,e＝1,...,N
with:KU＝F
其中，Vf是结构质量分数函数(目标函数)。ρe1是第e个单元的第一物理密度，是关于设计变量xe的函数；ve是第e个单元的体积；σe是第e个单元的许用应力；σl是材料的屈服应力；K是组合的全局刚度矩阵；U是单元节点位移向量；F是单元节点等效载荷的向量。
第二步骤，初始化离散后的各单元的设计变量xe，形成初始设计变量集。本实施例中，设计域中各单元的设计变量xe均初始化为0.5，非设计域中各单元的设计变量xe始终为1。
第三步骤，利用变密度法进行迭代设计直至得到最优设计变量集；每次迭代设计执行如下步骤：
S1：基于初始设计变量集或上次迭代形成的设计变量集，用密度过滤法(过滤函数)对各单元的设计变量进行过滤，得到单元过滤密度
Figure BDA0003689171460000062
该步骤为现有技术，其实现方法为本领域技术人员所知悉。本实施例中：
Figure BDA0003689171460000063
其中，w(xi)为权重系数，即第i个单元对第e个单元影响的加权系数；
w(xi)＝max(0,rmin-||ci-ce||2；
其中，rmin为预设的过滤半径，本实施例中取值为3.5；ci为第i个单元的中心位
置；ce为第e个单元的中心位置。
S2：用Heaviside函数对各单元的过滤密度
Figure BDA0003689171460000071
进行映射，得到第一密度ρe1：
Figure BDA0003689171460000072
其中，
Figure BDA0003689171460000073
用于控制映射曲线的光滑度，k为本次迭代的序号，k＝1时，
Figure BDA0003689171460000074
从第1次迭代起，每20次迭代更新一次，每次更新的
Figure BDA0003689171460000075
值为上一
Figure BDA0003689171460000076
值的1.19倍，但
Figure BDA0003689171460000077
值最大为8。
S3：用RAMP法对单元的弹性模量进行惩罚，计算各单元的第一弹性模量Ee1：
Figure BDA0003689171460000078
其中，Emax是设计材料的弹性模量值，Emin是为了避免奇异矩阵添加的小量，本实施例中Emin＝1e-6Emax；q1为第一惩罚因子，q1＝4。
S4：组装全局刚度矩阵，根据静力平衡方程KU＝F计算结构位移；
S5：计算各单元的第一应力σe1：
σe1＝DBue；
其中，D是轴对称单元的弹性系数矩阵；B是单元的几何函数矩阵；ue是单元的位移矩阵。
S6：计算各单元的第一米塞斯应力(Von Mises应力)σVM1：
Figure BDA0003689171460000079
其中，V是常数矩阵，表达式为：
Figure BDA00036891714600000710
S7：用RAMP法对各单元的第一米塞斯应力σVM1进行惩罚，得到各单元的第一惩罚应力
Figure BDA00036891714600000711
Figure BDA00036891714600000712
其中，q2为第二惩罚因子，q2＝-0.95。
S8：采用P-norm法聚合各单元的第一惩罚应力
Figure BDA00036891714600000713
得到聚合应力
Figure BDA00036891714600000714
Figure BDA00036891714600000715
其中，P是聚合系数，本实施例中，P＝8。
S9：对应力约束进行放松：
将应力约束放松为
Figure BDA0003689171460000081
其中，σl为材料的许用应力，c为本次迭代的放松系数，c的值每5次迭代更新一次，每次更新时
Figure BDA0003689171460000082
其中，
Figure BDA0003689171460000083
为本次迭代的第一最大预测应力：
Figure BDA0003689171460000084
其中，
Figure BDA0003689171460000085
为本次迭代的第二最大预测应力，第二最大预测应力
Figure BDA0003689171460000086
在本次迭代中通过如下步骤获得：
S9.1：用Heaviside函数对每个单元的过滤密度
Figure BDA0003689171460000087
进行映射，得到第二密度ρe2：
Figure BDA0003689171460000088
其中，
Figure BDA0003689171460000089
用于控制映射曲线的光滑度，＝1时，
Figure BDA00036891714600000810
从第1次迭代起，每20次迭代更新一次，每次更新的
Figure BDA00036891714600000811
值为上次
Figure BDA00036891714600000812
值加4，但
Figure BDA00036891714600000813
值最大为20；
S9.2：计算转变因子
Figure BDA00036891714600000814
Figure BDA00036891714600000815
其中，N为所有单元的总数量，Nb为第二密度ρe2大于第一阈值的单元的数量，Nw为第二密度ρe2小于第二阈值的单元的数量；所述第一阈值大于0.75，所述第二阈值小于0.25；本实施例中，第一阈值为0.9，第二阈值为0.1；
S9.3：计算各单元的第二弹性模量Ee2：
Figure BDA00036891714600000816
S9.4：组装全局刚度矩阵，根据静力平衡方程KU＝F计算结构位移；
S9.5：计算各单元的第二应力σe2(方法与S5相同，不再赘述)；
S9.6：计算各单元的第二米塞斯应力σVM2(方法与S6相同，不再赘述)；
S9.7：用线性法对各单元的第二米塞斯应力σVM2进行惩罚，得到各单元的第二惩罚应力
Figure BDA00036891714600000817
S9.8：计算本次迭代的第二最大预测应力
Figure BDA0003689171460000091
Figure BDA0003689171460000092
以上步骤S9.1至S9.8的过程允许与步骤S2至S8的过程并行计算。本实施例中，两者并行计算。
S10：对目标函数Vf进行灵敏度分析：
目标函数的灵敏度
Figure BDA0003689171460000093
采用伴随法计算全局应力约束的灵敏度：
应力的增广矩阵：
Figure BDA0003689171460000094
其中，λ为伴随向量；
全局应力的灵敏度：
Figure BDA0003689171460000095
其中，
Figure BDA0003689171460000096
Figure BDA0003689171460000097
用链式法则计算物理密度ρe对设计变量xe的导数：
Figure BDA0003689171460000098
S11：采用移动渐近线算法(即MMA优化求解器)对优化问题进行求解计算，更新各单元设计变量，得到并储存本次迭代形成的设计变量集；
S12：判断本次迭代是否满足退出迭代条件，本次迭代如满足退出迭代条件则退出迭代，并将本次迭代形成的设计变量集记为最优设计变量集；如不满足退出迭代条件则进行下一次迭代；具体地，判断本次迭代是否满足退出迭代条件的方法是对比本次迭代形成的设计变量集和上次迭代形成的设计变量集，或者对比本次迭代得到的目标函数值和上次迭代得到的目标函数值，并且第二最大预测应力小于等于材料许用应力；本实施例中，退出迭代条件为本次迭代形成的设计变量集中的所有设计变量与上次迭代形成的设计变量集中对应的设计变量的差值的绝对值均小于第三阈值，第三阈值取0.05。本实施例中，还规定最小迭代步数大于200步，最大迭代步数小于400步。
图4示出了根据实施例中涡轮盘根据最优设计变量集构造的涡轮盘结构。优化结构的质量分数为0.101。与原始的涡轮盘相比，在满足应力的要求下，结构的质量下降了48％。
第四步骤：根据最优设计变量集对结构进行光滑化。本实施例中，具体方法为：在最优设计变量集中，对于任一单元，如果xe<0.5，则该单元所在空间无材料，如果xe>0.5，则该单元所在空间具有全部材料，如果xe＝0.5，则该单元位于交界面。
光滑处理后，得到的涡轮盘结构如图5所示。在ANSYS软件对重构后的结构进行有限元分析，网格尺寸为1mm,单元的总数为3012，如图5所示。采用相同的荷载、边界条件和材料性质，重构结构的最大von Mises应力是1004Mpa，应力相对误差为0.4％，满足使用要求。
除实施例外，申请人还对步骤S9及步骤S9.1至S9.8中获取第一最大预测应力的方法进行了验证，见图6的验证例中的结构，具体验证方法为：
对图6所示的轴对称结构的横截面图进行有限元分析计算。图中展现了带锯齿边界和模糊区域的结构的具体尺寸大小，其中void部分代表单元密度为0，solid部分代表单元密度为1，0.3所指的单元代表其密度为0.3，0.5和0.7有着相同的含义。结构的旋转角速度为1047.2rad/s，在右侧边长为6mm处施加的载荷为100Mpa。结构的左下角的轴向位移约束为0。连接单元密度为0.5的对角线和单位密度为1的边线，如图6所示的边线构成的区域，边线围成的图形即为重构的结构。采用贴体网格计算，结构的最大应力为366.95Mpa。从下面表1中可以看出，随着Heaviside函数中
Figure BDA0003689171460000101
值逐渐增大，可以看到，采用非线性惩罚得到的结构最大应力值始终大于线性惩罚得到的应力值；采用线性惩罚得到的结构最大应力值波动较小，并且与ANSYS中在同一单元尺寸下进行有限元分析得到的最大应力值接近。
表1：
Figure BDA0003689171460000102
值与计算得到的应力值的关系。
Figure BDA0003689171460000103
通过以上对实施例和验证例可知，本发明相对于现有技术，主要的改进及其技术效果在于：首先，步骤S9中，对放松系数c获得方式的改进，包括通过步骤S9.1至S9.8基于预测应力法对
Figure BDA0003689171460000104
的计算。减轻了优化过程中锯齿边界和模糊区域的存在对结构应力场计算的影响。其中最主要的是根据设计变量小于0.1和大于0.9的单元数量与单元总数量的比值，决定对各单元的弹性模量采用线性惩罚还是非线性惩罚。引入最大预测应力，能够更充分地利用材料许用应力，提高优化设计质量，获得较现有技术质量更轻的设计方案。其次，步骤S7中，用RAMP法对各单元的的米塞斯应力进行非线性惩罚，且惩罚因子选用-0.95，能够在离心载荷约束下很好地提高低密度单元的灵敏度，有利于拓扑优化。
上述说明书和实施例的描述，用于解释本申请的保护范围，但并不构成对本申请保护范围的限定。

Claims (7)

1.考虑离心载荷和应力约束对结构进行优化设计的方法，其优化目标为在应力约束条件下结构质量最小，所述方法包括：对设计域内的结构进行离散的步骤、初始化离散后的各单元的设计变量xe形成初始设计变量集的步骤、利用变密度法进行迭代设计直至得到最优设计变量集的步骤和根据最优设计变量集对结构进行光滑化的步骤，其中，所有设计变量集中的所有单元的设计变量xe应满足0≤xe≤1；其特征是，每次迭代设计执行如下步骤：
S1：基于初始设计变量集或上次迭代形成的设计变量集，用密度过滤法对各单元的设计变量进行过滤，得到单元过滤密度
Figure FDA0003689171450000011
S2：用Heaviside函数对各单元的过滤密度
Figure FDA0003689171450000012
进行映射，得到第一密度ρe1：
Figure FDA0003689171450000013
其中，
Figure FDA0003689171450000014
用于控制映射曲线的光滑度，k为本次迭代的序号，k＝1时，
Figure FDA0003689171450000015
从第1次迭代起，每20次迭代更新一次，每次更新的
Figure FDA0003689171450000016
值为上一
Figure FDA0003689171450000017
值的1.19倍，但
Figure FDA0003689171450000018
值最大为8；
S3：用RAMP法计算各单元的第一弹性模量Ee1：
Figure FDA0003689171450000019
其中，Emax是设计材料的弹性模量值，Emin是为了避免奇异矩阵添加的小量；q1为第一惩罚因子，q1＝4；
S4：组装全局刚度矩阵，根据静力平衡方程计算结构位移；
S5：计算各单元的第一应力σe1；
S6：计算各单元的第一米塞斯应力σVM1；
S7：用RAMP法对各单元的第一米塞斯应力σVM1进行惩罚，得到各单元的第一惩罚应力
Figure FDA00036891714500000110
Figure FDA00036891714500000111
其中，q2为第二惩罚因子，q2＝-0.95；
S8：聚合各单元的第一惩罚应力
Figure FDA00036891714500000112
得到聚合应力
Figure FDA00036891714500000113
S9：对应力约束进行放松：
将应力约束放松为
Figure FDA00036891714500000114
其中，σl为材料的许用应力，c为本次迭代的放松系数，c的值每5次迭代更新一次，每次更新时
Figure FDA00036891714500000115
其中，
Figure FDA00036891714500000116
为本次迭代的第一最大预测应力，
Figure FDA0003689171450000021
其中，
Figure FDA0003689171450000022
为本次迭代的第二最大预测应力，第二最大预测应力
Figure FDA0003689171450000023
在本次迭代中通过如下步骤获得：
S9.1：用Heaviside函数对每个单元的过滤密度
Figure FDA0003689171450000024
进行映射，得到第二密度ρe2：
Figure FDA0003689171450000025
其中，
Figure FDA0003689171450000026
用于控制映射曲线的光滑度，k＝1时，
Figure FDA0003689171450000027
从第1次迭代起，每20次迭代更新一次，每次更新的
Figure FDA0003689171450000028
值为上次
Figure FDA0003689171450000029
值加4，但
Figure FDA00036891714500000210
值最大为20；
S9.2：计算转变因子
Figure FDA00036891714500000211
Figure FDA00036891714500000212
其中，N为所有单元的总数量，Nb为第二密度ρe2大于第一阈值的单元的数量，Nw为第二密度ρe2小于第二阈值的单元的数量；所述第一阈值大于0.75，所述第二阈值小于0.25；
S9.3：计算各单元的第二弹性模量Ee2：
Figure FDA00036891714500000213
S9.4：组装全局刚度矩阵，根据静力平衡方程计算结构位移；
S9.5：计算各单元的第二应力σe2；
S9.6：计算各单元的第二米塞斯应力σVM2；
S9.7：用线性法对各单元的第二米塞斯应力σVM2进行惩罚，得到各单元的第二惩罚应力
Figure FDA00036891714500000214
Figure FDA00036891714500000215
S9.8：计算本次迭代的第二最大预测应力
Figure FDA00036891714500000216
Figure FDA00036891714500000217
S10：对目标函数Vf进行灵敏度分析；
其中，目标函数
Figure FDA0003689171450000031
其中，ve为第e个单元的体积；
S11：采用移动渐近线算法对优化问题进行求解计算，更新各单元设计变量，得到并储存本次迭代形成的设计变量集；
S12：判断本次迭代是否满足退出迭代条件，本次迭代如满足退出迭代条件则退出迭代，并将本次迭代形成的设计变量集记为最优设计变量集；如不满足退出迭代条件则进行下一次迭代；
以上步骤S9.1至S9.8的过程允许与步骤S2至S8的过程并行计算。
2.如权利要求1所述的考虑离心载荷和应力约束对结构进行优化设计的方法，其特征是，根据最优设计变量集对结构进行光滑化的步骤采用如下方法：
在最优设计变量集中，对于任一单元，如果xe＜0.5，则该单元所在空间无材料，如果xe＞0.5，则该单元所在空间具有全部材料，如果xe＝0.5，则该单元位于交界面。
3.如权利要求1所述的考虑离心载荷和应力约束对结构进行优化设计的方法，其特征是，步骤S9.2中，第一阈值为0.9，第二阈值为0.1。
4.如权利要求1所述的考虑离心载荷和应力约束对结构进行优化设计的方法，其特征是，步骤S8中，采用P-norm法聚合各单元的第一惩罚应力
Figure FDA0003689171450000032
5.如权利要求1所述的考虑离心载荷和应力约束对结构进行优化设计的方法，其特征是，步骤S12中，判断本次迭代是否满足退出迭代条件的方法是对比本次迭代形成的设计变量集和上次迭代形成的设计变量集，或者对比本次迭代得到的目标函数值和上次迭代得到的目标函数值，并且第二最大预测应力小于等于材料许用应力。
6.如权利要求5所述的考虑离心载荷和应力约束对结构进行优化设计的方法，其特征是，退出迭代条件为本次迭代形成的设计变量集中的所有设计变量与上次迭代形成的设计变量集中对应的设计变量的差值的绝对值均小于第三阈值。
7.如权利要求6所述的考虑离心载荷和应力约束对结构进行优化设计的方法，其特征是，所述第三阈值为0.05。
CN202210655155.3A 2022-06-10 2022-06-10 考虑离心载荷和应力约束对结构进行优化设计的方法 Active CN114912226B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202210655155.3A CN114912226B (zh) 2022-06-10 2022-06-10 考虑离心载荷和应力约束对结构进行优化设计的方法
US18/169,144 US20230409767A1 (en) 2022-06-10 2023-02-14 Method of optimizing a design of a structure by considering centrifugal loads and stress constraints

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210655155.3A CN114912226B (zh) 2022-06-10 2022-06-10 考虑离心载荷和应力约束对结构进行优化设计的方法

Publications (2)

Publication Number Publication Date
CN114912226A true CN114912226A (zh) 2022-08-16
CN114912226B CN114912226B (zh) 2024-07-30

Family

ID=82771579

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210655155.3A Active CN114912226B (zh) 2022-06-10 2022-06-10 考虑离心载荷和应力约束对结构进行优化设计的方法

Country Status (2)

Country Link
US (1) US20230409767A1 (zh)
CN (1) CN114912226B (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100274537A1 (en) * 2009-04-24 2010-10-28 Caterpillar, Inc. Stress-based Topology Optimization Method and Tool
CN102043883A (zh) * 2010-12-29 2011-05-04 长沙理工大学 基于材料破损约束的连续体结构拓扑设计建模及优化设计方法
CN107273613A (zh) * 2017-06-15 2017-10-20 华中科技大学 一种基于应力惩罚和自适应体积的结构拓扑优化设计方法
CN109508495A (zh) * 2018-11-12 2019-03-22 华东交通大学 一种基于k-s函数的柔顺机构全局应力约束拓扑优化方法
CN110414165A (zh) * 2019-08-01 2019-11-05 华东交通大学 一种基于全局应力约束的多相材料柔顺机构拓扑优化方法
CN112177678A (zh) * 2020-09-25 2021-01-05 厦门大学 带双内环空腔的涡轮盘结构及其设计方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100274537A1 (en) * 2009-04-24 2010-10-28 Caterpillar, Inc. Stress-based Topology Optimization Method and Tool
CN102043883A (zh) * 2010-12-29 2011-05-04 长沙理工大学 基于材料破损约束的连续体结构拓扑设计建模及优化设计方法
CN107273613A (zh) * 2017-06-15 2017-10-20 华中科技大学 一种基于应力惩罚和自适应体积的结构拓扑优化设计方法
CN109508495A (zh) * 2018-11-12 2019-03-22 华东交通大学 一种基于k-s函数的柔顺机构全局应力约束拓扑优化方法
CN110414165A (zh) * 2019-08-01 2019-11-05 华东交通大学 一种基于全局应力约束的多相材料柔顺机构拓扑优化方法
CN112177678A (zh) * 2020-09-25 2021-01-05 厦门大学 带双内环空腔的涡轮盘结构及其设计方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
罗阳军;: "基于D-P准则的压力相关材料结构拓扑优化", 力学学报, no. 05, 18 September 2011 (2011-09-18) *

Also Published As

Publication number Publication date
US20230409767A1 (en) 2023-12-21
CN114912226B (zh) 2024-07-30

Similar Documents

Publication Publication Date Title
Stanford et al. Optimal topology of aircraft rib and spar structures under aeroelastic loads
CN105183996B (zh) 面元修正与网格预先自适应计算方法
Oyama et al. Transonic axial-flow blade optimization: Evolutionary algorithms/three-dimensional Navier-Stokes solver
Kersken et al. Time-linearized and time-accurate 3D RANS methods for aeroelastic analysis in turbomachinery
CN110008512B (zh) 一种考虑承载特性的负泊松比点阵结构拓扑优化方法
Sleesongsom et al. Aircraft morphing wing design by using partial topology optimization
CN109508495A (zh) 一种基于k-s函数的柔顺机构全局应力约束拓扑优化方法
CN112765732B (zh) 一种基于选区激光熔化工艺的航空叶片拓扑优化设计方法
Krawczyk et al. Fluid structure interaction of a morphed wind turbine blade
Gray et al. Geometrically nonlinear high-fidelity aerostructural optimization for highly flexible wings
CN104750948B (zh) 一种处理飞行器设计中多极值多约束问题的优化方法
CN112417773B (zh) 多级轴流膨胀机的多学科优化设计方法、装置及设备
CN108268728A (zh) 基于两步式改进粒子群优化算法的汽车尾门结构优化方法
CN115392094A (zh) 一种基于热力耦合的涡轮盘结构优化方法
CN114970366B (zh) 一种功能梯度超材料结构优化设计方法及系统
CN109502017B (zh) 一种拓扑优化仿生无人机及其设计方法
Jacobson et al. Flutter-Constrained Optimization with the Linearized Frequency-Domain Approach
CN114912226A (zh) 考虑离心载荷和应力约束对结构进行优化设计的方法
Miller et al. The development of a flatback wind turbine airfoil family
CN109815518A (zh) 基于转动惯量约束的飞行器舵面设计方法
CN106611078B (zh) 产品的具有时间步长控制方案的高效显式有限元分析
JP6548532B2 (ja) 数値的導関数を用いた構造的トポロジー最適化
CN111597724A (zh) 考虑频带约束的结构动力学拓扑优化方法、系统及应用
CN109033661B (zh) 一种叶轮的设计方法
CN114065423B (zh) 快速评估航空发动机风扇叶片颤振的方法

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