CN117350116A - 一种考虑惯性载荷的热力耦合拓扑优化方法 - Google Patents

Abstract

本发明提供了一种考虑惯性载荷的结构热力耦合拓扑优化方法，并将其应用于航空发动机涡轮盘结构的拓扑优化中，验证了所提出方法的有效性。首先，分析模型的热力学边界条件，通过ABAQUS构建结构模型并导出模型的物理信息；其次，分别给出了热力学材料参数的插值惩罚格式，并基于此构建了考虑热载荷、机械载荷以及惯性载荷的热力耦合拓扑优化模型；然后，基于ABAQUS求解热传导方程及热力耦合方程，并推导了拓扑设计变量的灵敏度信息；最后，采用基于梯度的移动渐近线算法更新设计变量，获得热力耦合结构的最优拓扑构型。该方法可以处理非规则设计域的问题，为实际应用中的结构优化设计提供了有效的工具。

Description

一种考虑惯性载荷的热力耦合拓扑优化方法

技术领域

本发明属于热力耦合拓扑优化设计领域，涉及一种考虑惯性载荷的热力耦合拓扑优化方法。

背景技术

现代装备结构的优化设计在轻量化、强韧性、耐热性等方面的高要求对结构设计方法提出了越来越多的挑战。拓扑优化作为一种非直观的概念设计方法，突破了传统尺寸优化和形状优化的设计理念，可以通过对有限设计域内有限材料的充分利用及合理布局，对结构进行创新构型设计。现代装备结构通常会同时承受机械载荷和温度载荷的作用，热力场耦合作用将对结构性能产生重要影响，仅考虑单物理场的设计可能导致结构失效，进而损害装备结构的性能。研究如何在保障结构实现高承载、耐高温等功能的同时满足结构轻量化的要求至关重要。

目前，结构拓扑优化研究得到了广泛的关注与发展，并发展了一系列高效的方法。然而，在大部分研究中，优化通常在给定外载荷的条件下进行，即结构所受的外载荷与结构拓扑不相关。而结构自重、惯性力、热应力等工程问题中实际存在的体力是与结构拓扑直接相关的，其在优化过程中随着结构拓扑的变化而变化，对结构拓扑设计有着至关重要的影响。例如，航空发动机关键结构在极端服役环境下，往往同时承载着高速旋转所带来的离心力以及高温所带来的热应力。在拓扑优化设计的过程中，必须考虑这些较大载荷所带来的影响。然而，考虑惯性载荷的热力耦合拓扑优化问题中往往伴随着低密度区域寄生效应的问题，出现的大量灰度区域会极大影响结构优化后的性能，并且难以加工制造。因此，针对航空发动机关键结构例如涡轮盘开展考虑惯性载荷的热力耦合结构拓扑优化设计具有重要意义。

发明内容

本发明的目的是提供一种考虑惯性载荷的热力耦合拓扑优化方法，能够解决航空发动机涡轮盘结构旋转工况下的热力耦合拓扑优化问题。首先，分析模型的热力学边界条件，通过ABAQUS构建结构模型并导出模型的物理信息；其次，分别给出了热力学材料参数的插值惩罚格式，并基于此构建了考虑热载荷、机械载荷以及惯性载荷的热力耦合拓扑优化模型；然后，基于ABAQUS求解热传导方程及热力耦合方程，并推导了拓扑设计变量的灵敏度信息；最后，采用基于梯度的移动渐近线算法更新设计变量，获得热力耦合结构的最优拓扑构型。该方法可以在热力耦合拓扑优化过程中考虑惯性力的影响，较大程度上改善了低密度区域寄生效应的问题，为工程中多场耦合结构优化设计提供有效的工具。

本发明提供了一种考虑惯性载荷的热力耦合拓扑优化方法，可以求解考虑惯性载荷的热力耦合拓扑优化问题。本发明包括如下步骤：

步骤1：通过ABAQUS构建结构模型并获取模型的物理信息，包括网格数据，外载荷分布及约束边界、传热边界；

步骤2：构造与热力学相关的材料参数插值模型。对于传热分析方程中的材料参数，分别采用幂指数形式的变密度法(Solid Isotropic Material with Penalization,SIMP)插值形式与有理近似模型(rational approximation of material properties,RAMP)插值形式进行构建，结构第e个单元的热传导系数与热膨胀系数如下：

λ(x_e)＝x_e ^hλ₀ (15)

其中，h和m为惩罚因子，λ₀为实体单元的热传导系数，α₀为实体单元的热膨胀系数，λ(x_e)表示拓扑变量相关单元的热传导系数，α(x_e)表示拓扑变量相关单元的热膨胀系数；对于热力耦合分析方程中的材料参数，先采用SIMP法构建密度插值函数，结构第e个单元的密度：

ρ(x_e)＝x_e ^sρ₀ (17)

其中，s为惩罚因子，ρ₀为实体单元的密度，ρ(x_e)表示拓扑变量相关单元的密度；

由于热载荷和惯性载荷均与设计变量相关，其值大小会随着拓扑变量的更新而变化，故改进的SIMP插值模型计算单元的弹性模量：

其中，p，q均为惩罚因子，E₀为实体单元的弹性模量，γ为经验参数，E(x_e)表示拓扑变量相关单元的弹性模量；

步骤3：计算结构的惯性载荷，第e个单元的惯性载荷如下：

其中，v_e为单元体积，ω为角速度，r_e为转动半径，θ_e为转动角。

步骤4：构建考虑惯性载荷的热力耦合拓扑优化模型：

其中，T为结构温度场，U为结构位移场，H和K分别表示传热刚度矩阵和力学刚度矩阵，P为等效温度载荷，F为结构的机械载荷，F_ΔT表示结构的等效热载荷，S表示结构受到的总载荷，x为拓扑设计变量，C为结构柔顺度，F_I表示结构惯性载荷，V表示结构体积，V₀为结构初始体积，β为体积分数，避免出现计算奇异；

步骤5：进行传热有限元分析：

H(x)T＝P (21)

步骤6：进行热力耦合有限元分析：

K(x)U＝F+F_ΔT(x)+F_I(x)＝S(x) (22)

步骤7：对热力耦合拓扑优化模型进行灵敏度分析：

对式(8)两端同时对拓扑变量求导可得：

将由式(10)带入式(9)中，可求得目标函数关于拓扑变量的灵敏度：

步骤8：基于目标函数关于设计变量的灵敏度信息，采用基于梯度的MMA优化算法对设计变量进行更新迭代；

步骤9：判断收敛性，若不收敛，则回到步骤2，直至计算收敛后，得到热力耦合结构最优的拓扑构型。

进一步地，步骤7中求解目标函数关于设计变量的灵敏度信息，具体步骤包括：

基于对应的密度插值格式，可分别得出和/>如下：

其中表示单位弹性模量单元的刚度矩阵。等效热载荷关于拓扑变量的导数为：

其中，ε_e表示拓扑变量相关单元的热应变，ε₀表示实体单元的热应变，B表示形函数矩阵，D_e表示拓扑变量相关单元的弹性矩阵，D₀表示实体单元的弹性矩阵。将式(12)，(13)，(14)带入式(11)中，即可获得目标函数关于拓扑设计变量的灵敏度。

本发明的有益效果是：

(1)本发明构建了一种考虑惯性载荷的热力耦合拓扑优化模型，在体积约束条件下优化结构的柔顺度；

(2)本发明针对低密度区域寄生效应的问题，提出了一种光滑连续的材料插值模型，该模型能有效控制低密度区域灰度单元产生，并且易于求导。在此基础上，推导了拓扑设计变量的灵敏度信息，提出了一种高效的热力耦合拓扑优化求解算法；

(3)本发明针对考虑惯性载荷的热力耦合拓扑优化设计问题，联合ABAQUS进行了热力耦合有限元分析和拓扑优化设计，该方法可以处理非规则设计域的问题，为实际应用中的结构优化设计提供了有效的工具。

附图说明

图1为本发明的流程框图。

图2为结构设计域及边界条件示意图。

图3为第一个拓扑迭代步中的温度场求解结果。

图4为第一个拓扑迭代步中的位移场求解结果。

图5为经过第一步迭代更新的拓扑优化结果。

图6为基于惯性载荷的热力耦合拓扑优化结果。

图7为基于惯性载荷的热力耦合拓扑优化迭代曲线。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

该实施例提供了一种考虑惯性载荷的热力耦合拓扑优化方法，并将其应用于航空发动机涡轮盘的拓扑优化设计案例中。下面结合附图及涡轮盘拓扑优化设计实例对本发明作进一步详细说明：

图3给出了内孔半径0.6m，轮毂半径2.4m，厚度为0.01m的涡轮盘结构。图中虚线部分代表固定温度边界，外表面介质温度为T₂，孔洞内部介质温度为T₁。设计域边界的四个节点处施加如图所示的四个固定载荷F。本发明方法具体的计算步骤详细说明如下：

步骤1：通过ABAQUS构建模型并获取模型的物理信息，包括网格数据，受力情况及约束边界、传热边界。

该步骤中：通过ABAQUS进行前处理并导出涡轮盘模型的相关物理信息，获得6个txt文件，分别为elemnode文件，nodecoor文件，mat文件，disp文件，thermal文件，rotationset文件，依次包含了模型的网格单元(设计域网格划分为8240个平面应力单元)；节点坐标信息；材料属性信息(弹性模量E₀＝210Gpa，泊松比v＝0.3，单元厚度t＝0.01m)；力学边界信息(设计域圆心固定约束，加载力F＝10000N)；温度边界信息(外表面介质温度为T₂＝200℃，孔洞内部介质温度为T₁＝20℃)；旋转工况信息(涡轮盘转速为ω＝3000rpm)。

步骤2：构造与热力学相关的材料参数插值模型，对于传热分析方程中的材料参数，分别采用幂指数形式的变密度法(Solid Isotropic Material with Penalization,SIMP)插值形式与有理近似模型(Rational Approximation of Material Properties,RAMP)插值形式进行构建，构造涡轮盘结构第e个单元的热传导系数λ与热膨胀系数α如下：

λ(x_e)＝x_e ^hλ₀ (1)

其中，h和m为惩罚因子，λ₀为实体单元的热传导系数，α₀为实体单元的热膨胀系数，x_e为第e个单元的拓扑变量，取值范围在0到1之间。λ(x_e)表示拓扑变量相关单元的热传导系数，α(x_e)表示拓扑变量相关单元的热膨胀系数；

在该实施例中，实体材料的热膨胀系数α₀＝10.67×10^-6m/℃，导热系数λ₀＝50.37W·m^-1·℃^-1。h和m为惩罚因子，分别取h＝1，m＝5。所有单元的拓扑变量初始值均设置为1。

对于热力耦合分析方程中的材料参数，先采用SIMP法构建密度插值函数，结构第e个单元的密度：

ρ(x_e)＝x_e ^sρ₀ (29)

其中，s为惩罚因子，ρ₀为实体单元的密度，ρ(x_e)表示拓扑变量相关单元的密度。

此外，由于热载荷和惯性载荷均与设计变量相关，其值大小会随着拓扑变量的更新而变化，惯性载荷和热载荷的共同作用会对拓扑优化产生较大的影响。

主要表现为以下几个问题：

①在优化过程中柔顺度目标函数不是单调函数；

②低密度区域会产生寄生效应；③优化过程中大量中间密度单元的惯性力和热应力对拓扑优化迭代造成影响。SIMP插值模型不能有效解决这一类问题，该实施例提出了一种改进的SIMP插值模型计算单元的弹性模量：

其中，p，q均为惩罚因子，E₀为实体单元的弹性模量，γ为经验参数，E(x_e)表示拓扑变量相关单元的弹性模量；

该实施例中，实体材料的密度ρ₀＝8010g/cm³，弹性模量E₀＝210GPa。s，p，q，γ均为惩罚因子，分别取s＝1，p＝2，q＝5，γ＝8。

步骤3：计算结构的惯性载荷，第e个单元的惯性载荷如下：

其中，v_e为单元体积，ω为角速度，r_e为第e个单元转动半径，θ_e为第e个单元的转动角，θ_e＝[cosθ⁽ⁱ⁾sinθ⁽ⁱ⁾cosθ⁽ⁱ⁾sinθ⁽ⁱ⁾cosθ⁽ⁱ⁾sinθ⁽ⁱ⁾cosθ⁽ⁱ⁾sinθ⁽ⁱ⁾]^T，cosθ⁽ⁱ⁾为第i个节点的方向余弦。

在该实施例中，转速为ω＝3000rpm。对于该涡轮盘结构的8240个平面应力单元，每个有限单元都有相应的单元体积、转动半径以及转动角，以1号单元为例，其单元体积v₁＝2.5×10^-5m²，转动半径r₁＝2.3747m；

转动角θ₁＝[0.9493,0.3143,0.9493,0.3143,0.9493,0.3143,0.9493,0.3143]^T；

步骤4：针对涡轮盘结构构建考虑惯性载荷的热力耦合拓扑优化模型：

其中，T为结构温度场，U为结构位移场，H和K分别表示传热刚度矩阵和力学刚度矩阵，P为等效温度载荷，F为结构的机械载荷，F_ΔT表示结构的等效热载荷，S表示结构受到的总载荷，x为拓扑设计变量，C为结构柔顺度，F_I表示结构惯性载荷，V表示结构体积，V₀为结构初始体积，β为体积分数，设置为0.5，避免出现计算奇异；

步骤5、将插值后的材料参数λ(x_e)和α(x_e)传进ABAQUS，并进行传热有限元分析：

H(x)T＝P (7)

获得当前迭代步涡轮盘的温度场T，在第一个迭代步中，所求得的温度场如图3所示；

步骤6、再次将插值后的材料参数E(x_e)，ρ(x_e)传进ABAQUS，进行基于旋转工况的热力耦合有限元分析：

K(x)U＝F+F_ΔT(x)+F_I(x)＝S(x) (8)

可获得当前迭代步涡轮盘的位移场U，在第一个迭代步中，所求得的位移场如图4所示；

步骤7、对热力耦合拓扑优化模型进行灵敏度分析：

于对应的密度插值格式，可分别得出和/>如下：

等效热载荷关于拓扑变量的导数为：

同样的，对于该涡轮盘结构的8240个平面应力单元，每个有限单元都可求得相应的灵敏度值，以1号单元为例，所获得的目标函数灵敏度值为1328.9。

步骤8、基于目标函数对设计变量的灵敏度信息，采用基于梯度的MMA优化算法对设计变量进行更新迭代。第一个迭代步所求得的拓扑结果如图5所示；

步骤9、判断收敛性。若不收敛，则回到步骤3。直至计算收敛后，得到热力耦合结构最优的拓扑构型。

图6中给出了图2所示设计域的考虑惯性载荷的拓扑优化中间迭代过程和优化设计结果，其中图6中(d)是本方法得到的最终的拓扑结构构型。从图7所示的迭代历史可以看出，所提出的考虑惯性载荷的热力耦合拓扑优化方法具有良好的稳定性和收敛性。拓扑优化迭代过程表明，优化过程中低密度单元得到有效控制，形成的拓扑构型较为清晰，寄生效应显著减少，优化迭代过程也较为平稳。采用本方法得到优化结构更能满足航空发动机涡轮盘等在高速旋转工况下的性能需求。同时，本方法可以充分利用ABAQUS强大的前后置功能和有限元分析功能，进行大型复杂结构拓扑优化设计。

本发明针对惯性载荷的结构热力耦合拓扑优化问题，提出了一种针对弹性模量的改进材料插值惩罚模型，降低了由于拓扑相关载荷引起的寄生效应对优化构型的影响。首先，考虑结构热载荷、机械载荷以及惯性载荷，构造了相应的热力耦合拓扑优化模型；其次，研究了关于拓扑密度的材料插值惩罚格式组合方式，并基于此推导了结构柔顺度对拓扑设计变量的灵敏度；最后，采用基于梯度的优化算法进行拓扑变量迭代更新，最终获取热力耦合结构的拓扑构型。

以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围之内。

Claims (2)

1.一种考虑惯性载荷的热力耦合拓扑优化方法，其特征在于，包括如下步骤：

步骤1：通过ABAQUS构建结构模型并获取模型的物理信息，包括网格数据，外载荷分布及约束边界、传热边界；

步骤2：构造与热力学相关的材料参数插值模型，对于传热分析方程中的材料参数，分别采用幂指数形式的变密度法(Solid Isotropic Material with Penalization,SIMP)插值形式与有理近似模型(Rational Approximation of Material Properties,RAMP)插值形式进行构建，结构第e个单元的热传导系数与热膨胀系数如下：

λ(x_e)＝x_e ^hλ₀ (1)

其中，h和m为惩罚因子，λ₀为实体单元的热传导系数，α₀为实体单元的热膨胀系数，λ(x_e)表示拓扑变量相关单元的热传导系数，α(x_e)表示拓扑变量相关单元的热膨胀系数；对于热力耦合分析方程中的材料参数，先采用SIMP法构建密度插值函数，结构第e个单元的密度：

ρ(x_e)＝x_e ^sρ₀ (3)

其中，s为惩罚因子，ρ₀为实体单元的密度，ρ(x_e)表示拓扑变量相关单元的密度；

由于热载荷和惯性载荷均与设计变量相关，其值大小会随着拓扑变量的更新而变化，改进的SIMP插值模型计算单元的弹性模量如下：

其中，p，q均为惩罚因子，E₀为实体单元的弹性模量，γ为经验参数，E(x_e)表示拓扑变量相关单元的弹性模量；

步骤3：计算结构的惯性载荷，第e个单元的惯性载荷F_I ^(e)如下：

F_I ^(e)＝v_eρ(x_e)ω²r_eθ_e (5)

其中，v_e为单元体积，ω为角速度，r_e为转动半径，θ_e为转动角；

步骤4：构建考虑惯性载荷的热力耦合拓扑优化模型：

其中，T为结构温度场，U为结构位移场，H和K分别表示传热刚度矩阵和力学刚度矩阵，P为等效温度载荷，F为结构的机械载荷，F_ΔT表示结构的等效热载荷，S表示结构受到的总载荷，x为拓扑设计变量，C为结构柔顺度，F_I表示结构惯性载荷，V表示结构体积，V₀为结构初始体积，β为体积分数，避免出现计算奇异；

步骤5：基于ABAQUS进行传热有限元分析：

H(x)T＝P (7)

步骤6：基于ABAQUS进行热力耦合有限元分析：

K(x)U＝F+F_ΔT(x)+F_I(x)＝S(x) (8)

步骤7：对热力耦合拓扑优化模型进行灵敏度分析：

对式(8)两端同时对拓扑变量求导可得：

将由式(10)带入式(9)中，可求得目标函数关于拓扑变量的灵敏度：

步骤8：基于目标函数关于设计变量的灵敏度信息，采用基于梯度的MMA优化算法对设计变量进行更新迭代；

步骤9：判断收敛性，若不收敛，则回到步骤2，直至计算收敛后，得到热力耦合结构最优的拓扑构型。

2.根据权利要求1所述的一种考虑惯性载荷的热力耦合拓扑优化方法，其特征在于：

步骤7中求解目标函数关于设计变量的灵敏度信息，具体步骤包括：

基于对应的密度插值格式，可分别得出和/>如下：

其中表示单位弹性模量单元的刚度矩阵。等效热载荷关于拓扑变量的导数为：

其中，ε_e表示拓扑变量相关单元的热应变，ε₀表示实体单元的热应变，B表示形函数矩阵，D_e表示拓扑变量相关单元的弹性矩阵，D₀表示实体单元的弹性矩阵，将式(12)，(13)，(14)带入式(11)中，即可获得目标函数关于拓扑设计变量的灵敏度。