CN115935740A - 一种考虑材料温度依赖性的热力耦合拓扑优化方法 - Google Patents

Abstract

本发明提出一种考虑材料温度依赖性的热力耦合拓扑优化方法。首先基于ABAQUS定义结构模型及热力学边界条件并导出相关物理信息；其次构造考虑温变特性的材料参数插值模型，并基于此给出考虑温度依赖性的热力耦合拓扑优化模型；然后通过求解温变相关的非线性热传导方程及热力耦合方程，计算热力耦合结构的性能响应，进而求解拓扑优化目标函数；随后基于伴随变量法推导了拓扑设计变量的灵敏度信息；最后采用基于梯度的移动渐近线(MMA)算法更新设计变量，最终获取热力耦合结构的拓扑构型。本发明研究了一种考虑材料温度依赖性的热力耦合拓扑优化方法，可以处理涡轮盘结构的材料特性参数随温度变化的问题，有效实现了大温度梯度下热力耦合结构的拓扑优化设计。

Description

一种考虑材料温度依赖性的热力耦合拓扑优化方法

技术领域

本发明属于热力耦合拓扑优化设计领域，涉及一种考虑材料温度依赖性的热力耦合拓扑优化方法。

背景技术

拓扑优化是一种在给定的约束条件下，在给定的设计域内寻找材料的最佳分布，从而使结构的性能达到最优的结构优化方法。拓扑优化方法通常是以结构材料的有无(连续体拓扑优化)或者杆系结构的节点布局(离散体拓扑优化)为优化对象，一般应用于结构的概念设计阶段。因此，拓扑优化方法相比于传统的尺寸优化和形状优化等结构优化方法具有更多的设计自由度和更大的设计空间，能够有效地缩短结构的设计周期、提高材料的使用效率、降低结构的开发成本以及有效提高结构的性能，是一个具有生命力和挑战性的研究方向。

在最近几十年，现代加工方法与技术的发展为复杂结构的加工制造提供了可能，

与此同时，拓扑优化理论及方法也得到了不断的完善与发展。因此，拓扑优化方法已经逐渐成为当前产品设计和开发的重要工具，并已经被广泛应用于车辆、船舶、桥梁以及航空航天等领域。较早的拓扑优化通常仅考虑结构的传热性能或机械性能，然而在现代高端装备的设计中，结构通常会同时承受机械载荷和温度载荷的作用。热力场耦合作用将对结构性能产生重要影响，仅考虑单物理场的设计可能导致结构失效，进而损害装备结构的性能。

航空发动机是工业皇冠上的明珠，其在使用中会面临很高的温度和载荷，且结构温度值总是在从较低值到较高值的较宽范围内变化。在这种情况下，机械和热材料特性(弹性模量、热传导系数、热膨胀系数等)将根据大的温度梯度而显着变化。目前，热力耦合结构拓扑优化通常是基于恒定的材料特性假设，然而当整个结构的温度梯度很大时，材料属性不变的假设不再合理，会导致分析和优化误差较大。因此，针对航空发动机关键结构例如涡轮盘开展考虑材料温度依赖性的热力耦合结构拓扑优化设计具有重要意义。

发明内容

为解决传统热力耦合拓扑优化中材料性能常数假设的局限性，本发明的目的是提供一种考虑材料温度依赖性的热力耦合拓扑优化方法，能够解决航空发动机涡轮盘结构大温度梯度下的热力耦合拓扑优化问题。首先，基于ABAQUS的前处理功能分析模型并导出其物理信息，并基于变密度法给出了考虑温变特性的材料插值模型；随后，构造了考虑温度依赖性的热力耦合拓扑优化模型，并通过求解非线性热传导方程及热力耦合方程，计算结构性能响应；最后，进行了灵敏度分析并采用基于梯度的移动渐近线(MMA)算法进行优化求解。该方法可以在热力耦合拓扑优化过程中考虑到材料热力学特性参数受温度影响的实际情况，为工程中多场耦合结构优化设计提供有效的工具。

本发明提供了一种考虑材料温度依赖性的热力耦合拓扑优化方法，可以求解考虑大温度梯度的热力耦合拓扑优化问题。本发明包括如下步骤：

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

步骤2：基于密度插值模型，构建结构的第e个单元的材料弹性模量E如下：

其中，E₀为实体材料的弹性模量，E_min为空相材料的弹性模量，p为惩罚因子，N为有限单元数，

表示通过密度过滤和Heaviside投影过滤后得到的单元物理密度。密度过滤函数可表示为：

其中N_e为单元e在过滤半径内的单元集合，v_i与ρ_i为前述集合中第i个单元的体积与密度，权函数H_ei可以表示为：

H_ei＝max(0,r_min-r_ei)  (12)

其中r_ei为单元e中心到单元i中心的距离，r_min为过滤半径。经过密度过滤之后，中间变量

还需经过Heaviside投影过滤，

可由以下表达式计算得到：

其中ρ₀过滤密度的阈值，β为投影参数；

步骤3：构造与温度相关的材料参数插值模型，结构第e个单元的热传导系数λ与热膨胀系数α如下：

其中，λ(T)表示实体材料随温度T变化的热传导系数，α(T)表示实体材料随温度T变化的热膨胀系数；

步骤4：构建考虑材料温度依赖性的热力耦合拓扑优化模型：

其中C表示结构柔度，U为节点位移，K为总体刚度矩阵，T为节点温度，P为根据给定的热边界条件计算得到的温度载荷。H为总传热系数矩阵，F_ΔT表示等效热载荷，拓扑设计变量ρ_e表示第e个单元的密度，V₀表示结构初始体积，V表示结构体积，f为体积分数；

步骤5：进行非线性传热有限元分析，总传热系数矩阵H可通过单元传热系数矩阵H_e组装得到：

式中

为单元热应变矩阵，λ为单元传热系数矩阵，λ₀为传热系数矩阵的常数矩阵，注意到传热系数矩阵的非线性部分，应用牛顿-拉夫森迭代方法求解非线性残差方程：

R＝H(T)T-P＝0  (18)

步骤6：对热力耦合有限元方程，等效热载荷

可以通过以下方法得到：

其中B_e和D_e为第e个单元的应变-位移矩阵和弹性矩阵，D₀是实体材料的弹性矩阵；

其中ε_e为热应变，α为热膨胀系数，T_ref为参考温度，φ^T为单元组装矩阵。由于材料的弹性性能和热膨胀系数都与物理密度有关，式(20)可以表示为：

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

KU＝F+F_ΔT                           (23)

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

对KU＝F+F_ΔT两端同时对拓扑变量求导可得：

将由式(31)带入式(24)中，可求得目标函数关于物理密度的灵敏度：

步骤9：考虑过滤的影响，由链式法则推导目标函数关于设计变量的灵敏度：

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

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

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

等效热载荷关于物理密度的导数为：

其展开可以进一步表示为：

为了求解方程中的未知项

首先需要利用步骤5中的非线性残差方程(9)得到：

于是可以得到：

其中κ₁为拉格朗日乘子，基于伴随变量法，令：

于是，式(28)可以化为：

将式(30)和(26)带入式(32)求得拓扑优化目标函数关于物理密度的灵敏度：

其中U_e表示单元位移矩阵。

本发明的有益效果是：

1、本发明将真实物理场中材料属性变化的特点考虑到结构拓扑优化设计中，构建了一种考虑温度依赖性的热力耦合拓扑优化模型；

2、本发明结合伴随方程法和牛顿-拉夫森迭代法，针对该非线性拓扑优化问题推导了拓扑设计变量的灵敏度信息，提出了一种高效的热力耦合拓扑优化求解算法；

3、本发明针对考虑材料温度依赖性的热力耦合拓扑优化设计问题，使用ABAQUS二次开发功能模块进行了热力耦合仿真计算和拓扑优化设计，该方法可以处理任意非规则设计域的问题，为实际应用中的结构优化设计提供了有效的工具。

附图说明

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

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

图3为材料温度依赖性的热力耦合拓扑优化结果示意图。

图4为材料温度依赖性的热力耦合拓扑优化迭代曲线示意图。

具体实施方式

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

该实施例提出了一种考虑材料温度依赖性的热力耦合拓扑优化方法。首先，基于材料的温变性定义相关材料系数的插值模型，并构建了考虑材料温度依赖性的拓扑优化模型；然后，基于牛顿-拉夫森迭代法求解非线性传热方程，并计算热力耦合结构的目标柔度值；最后，基于结合伴随方程法推导目标函数关于结构拓扑设计变量的灵敏度，并采用基于梯度的优化算法求解，最终获取热力耦合结构的拓扑构型。下面结合附图及航空发动机涡轮盘的拓扑优化设计具体实例对本发明作进一步详细说明：

图3给出了内孔半径0.6m，轮毂半径2.4m，厚度为0.01m的涡轮盘设计域，灰色区域为非设计域。外表面介质温度为300℃，孔洞内部介质温度为20℃。本发明方法具体的计算步骤详细说明如下：

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

该步骤中：通过ABAQUS进行前处理并导出涡轮盘模型的相关物理信息，获得5个txt文件，分别为elemnode文件，nodecoor文件，mat文件，disp文件，thermal文件，包含了模型的网格单元，设计域网格划分为8240个平面应力(CPS4)单元；节点坐标信息；材料属性(弹性模量E₀＝210Gpa，泊松比v＝0.3，单元厚度t＝0.01m)；边界约束信息，设计域圆心固定约束；加载力F的类型为集中力，大小为100kN；

步骤2、基于密度插值模型，构造涡轮盘结构第e个单元的材料弹性模量E如下：

其中，E₀为实体材料的弹性模量即210GPa，E_min为空相材料的弹性模量，E_min＝10^-5，p为惩罚因子p＝3，N为有限单元数即8240，

表示通过密度过滤和Heaviside投影过滤后得到的单元物理密度。密度过滤函数可表示为：

其中

为中间变量，N_e为单元e在过滤半径内的单元集合，v_i与ρ_i为前述集合中第i个单元的体积与密度，权函数H_ei可以表示为：

H_ei＝max(0,r_min-r_ei)  (3)

其中r_ei为单元e中心到单元i中心的距离，r_min为过滤半径，设置为3。经过密度过滤之后，中间变量

还需经过Heaviside投影过滤，

可由以下表达式计算得到：

其中ρ₀过滤密度的阈值，β为投影参数，其初始值设置为1；

步骤3、构造与温度相关的材料参数插值模型，涡轮盘结构第e个单元的热传导系数λ与热膨胀系数α如下：

其中热传导系数λ(T)＝2+0.1T W/(m·K)，热膨胀系数α(T)＝11.7×10^-6+11.7×10^-6T/℃；

步骤4、针对涡轮盘结构构建考虑材料温度依赖性的热力耦合拓扑优化模型：

其中C表示柔度，U为节点位移，K为总体刚度矩阵，T为节点温度，P为根据给定的热边界条件计算得到的温度载荷。H为总传热系数矩阵，F_ΔT表示等效热载荷，拓扑设计变量ρ_e表示第e个单元的密度，V₀表示结构初始体积，V表示结构体积，f为体积分数，设置为0.4；

步骤5、将插值后的材料参数

传进ABAQUS，并进行非线性传热有限元分析，总传热系数矩阵H可通过单元传热系数矩阵H_e组装得到：：

式中

为单元热应变矩阵，λ为单元传热系数矩阵，λ₀为传热系数矩阵的常数矩阵。注意到传热系数矩阵的非线性部分，应用牛顿-拉夫森迭代方法求解非线性残差方程：

R＝H(T)T-P＝0                          (9)

非线性方程可如下求解：

在第k个迭代步，解T^(k)可能不能精确拟合方程，即R(T^(k))≠0。为了得到第(k+1)次迭代步骤的解，残差用一阶泰勒级数展开线性逼近：

T(^k+1)＝T(^k)+ΔT(^k)                          (11)

其中

为切线刚度矩阵，可以用K_T表示。ΔT^(k)是可以通过R(T^(k+1))泰勒展开等于0来计算的修正：

为了获得切线刚度矩阵，对非线性残差方程求导可得：

为了求解

构造A矩阵：

A的第k行第l列元素a_kl可以表示为：

可以进一步写为：

其中温度T_e为第e个单元中心的温度，可以由节点温度得到：

其中

是第n个节点的形函数。利用T_e的表达式可以将a_kl的表达式改写为：

最终可以求得A矩阵：

结合以上求解内容，可获得当前迭代步涡轮盘的温度场；

步骤6、对热力耦合有限元方程，等效热载荷

可以通过以下方法得到：

其中B_e和D_e为第e个单元的应变-位移矩阵和弹性矩阵，D₀是实体材料的弹性矩阵，

其中ε_e为热应变，α为热膨胀系数，T_ref为参考温度，φ^T为单元组装矩阵，由于材料的弹性性能和热膨胀系数都与物理密度有关，等效热载荷可以进一步表示为：

步骤7、再次将插值后的材料参数

传进ABAQUS，进行热力耦合有限元分析：

KU＝F+F_ΔT                            (23)

可获得当前迭代步涡轮盘的位移场；

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

等效热载荷关于物理密度的导数为：

其展开可以进一步表示为：

利用式(5)中的非线性残差方程求解

简化求解获得：

其中κ₁为拉格朗日乘子，基于伴随变量法，令：

于是得到

的表达式：

对KU＝F+F_ΔT两端同时对拓扑变量求导可得：

将上式代入到目标函数求导的表达式中，可得目标函数灵敏度的表达式

再将

代入到目标函数求导的表达式中，可得目标函数关于物理密度的灵敏度：

其中U_e表示单元位移矩阵。

步骤9、根据多场过滤的链式法则推导完整的目标函数敏度：

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

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

图3中给出了图2所示设计域的考虑材料温度依赖性的拓扑优化迭代过程和设计结果，其中图3中d是本方法得到的最终的拓扑结构构型。从迭代历史可以看出，所提出的考虑材料温度依赖性的热力耦合拓扑优化方法具有良好的稳定性和收敛性。拓扑优化迭代过程表明，结构整体构型会在较少的迭代步显现，后面的迭代计算会进一步完善构型的细节处理，说明了求解方法的计算高效性。采用本方法得到优化结构更能满足航空发动机涡轮盘等在大温度变化环境下工作的性能需求。并且，本方法可以充分利用ABAQUS强大的前后置功能和有限元分析功能，进行大型复杂结构拓扑优化设计。

本发明针对材料特性参数随温度变化的问题，提出了一种考虑材料温度依赖性的热力耦合拓扑优化方法。首先，基于ABAQUS定义结构模型及热力学边界条件并导出相关物理信息；其次，构造考虑温变特性的材料参数插值模型以及对应的热力耦合拓扑优化模型；然后，通过求解温变相关的非线性热传导方程及热力耦合方程，计算拓扑优化目标函数；随后，基于伴随变量法推导了拓扑设计变量的灵敏度信息；最后，采用基于梯度的优化算法更新设计变量，最终获取热力耦合结构的拓扑构型。

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

Claims (2)

1.一种考虑材料温度依赖性的热力耦合拓扑优化方法，其特征在于，包括如下步骤：

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

步骤2：基于密度插值模型，构建结构的第e个单元的材料弹性模量E如下：

其中，E₀为实体材料的弹性模量，E_min为空相材料的弹性模量，p为惩罚因子，N为有限单元数，

表示通过密度过滤和Heaviside投影过滤后得到的单元物理密度。密度过滤函数可表示为：

其中N_e为单元e在过滤半径内的单元集合，v_i与ρ_i为前述集合中第i个单元的体积与密度，权函数H_ei可以表示为：

H_ei＝max(0,r_min-r_ei)  (3)

其中r_ei为单元e中心到单元i中心的距离，r_min为过滤半径。经过密度过滤之后，中间变量

还需经过Heaviside投影过滤，

可由以下表达式计算得到：

其中ρ₀过滤密度的阈值，β为投影参数；

步骤3：构造与温度相关的材料参数插值模型，结构第e个单元的热传导系数λ与热膨胀系数α如下：

其中，λ(T)表示实体材料随温度T变化的热传导系数，α(T)表示实体材料随温度T变化的热膨胀系数；

步骤4：构建考虑材料温度依赖性的热力耦合拓扑优化模型：

其中C表示结构柔度，U为节点位移，K为总体刚度矩阵，T为节点温度，P为根据给定的热边界条件计算得到的温度载荷。H为总传热系数矩阵，F_ΔT表示等效热载荷，拓扑设计变量ρ_e表示第e个单元的密度，V₀表示结构初始体积，V表示结构体积，f为体积分数；

步骤5：进行非线性传热有限元分析，总传热系数矩阵H可通过单元传热系数矩阵H_e组装得到：

式中

为单元热应变矩阵，λ为单元传热系数矩阵，λ₀为传热系数矩阵的常数矩阵，注意到传热系数矩阵的非线性部分，应用牛顿-拉夫森迭代方法求解非线性残差方程：

R＝H(T)T-P＝0 (9)

步骤6：对热力耦合有限元方程，等效热载荷

可以通过以下方法得到：

其中B_e和D_e为第e个单元的应变-位移矩阵和弹性矩阵，D₀是实体材料的弹性矩阵；

其中ε_e为热应变，α为热膨胀系数，T_ref为参考温度，φ^T为单元组装矩阵。由于材料的弹性性能和热膨胀系数都与物理密度有关，式(20)可以表示为：

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

KU＝F+F_ΔT                           (23)

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

对KU＝F+F_ΔT两端同时对拓扑变量求导可得：

将由式(31)带入式(24)中，可求得目标函数关于物理密度的灵敏度：

步骤9：考虑过滤的影响，由链式法则推导目标函数关于设计变量的灵敏度：

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

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

2.根据权利要求1所述的一种考虑材料温度依赖性的热力耦合拓扑优化方法，其特征在于：步骤8中求解目标函数关于设计变量的灵敏度信息，具体步骤包括：

等效热载荷关于物理密度的导数为：

其展开可以进一步表示为：

为了求解方程中的未知项

首先需要利用步骤5中的非线性残差方程(9)得到：

于是可以得到：

其中κ₁为拉格朗日乘子，基于伴随变量法，令：

于是，式(28)可以化为：

将式(30)和(26)带入式(32)求得拓扑优化目标函数关于物理密度的灵敏度：

其中U_e表示单元位移矩阵。

Images