TOPOLOGY OPTIMIZATION METHOD FOR INTEGRATED THERMAL PROTECTION STRUCTURE CONSIDERING TRANSIENT TEMPERATURE AND STRESS CONSTRAINTS
-
摘要: 一体化热防护结构通常处于严酷的非稳态热环境, 热载荷作用的时间效应(即瞬态热效应)明显. 为了避免瞬态热分析的巨大计算消耗, 以往的一体化热防护结构优化设计研究通常将瞬态传热等效为相同热边界条件下的稳态传热, 将稳态传热分析的温度场作为设计热载荷. 然而, 已有的研究表明稳态传热无法准确等效瞬态传热的作用效果, 瞬态热效应对结构设计结果具有重要影响. 文章研究了考虑瞬态热效应的一体化热防护结构优化设计问题, 建立一种考虑瞬态温度和应力约束的一体化热防护结构拓扑优化方法. 该方法以SIMP (solid isotropic material with penalization) 法为基础, 构建两种针对一体化热防护结构的热弹性结构拓扑优化模型: (1)考虑材料体积分数、最大应力和底面最大温度约束, 以最小化结构应变能为目标的刚度设计模型; (2)考虑最大应力和底面最大温度约束, 以最小化材料体积分数为目标的轻量化设计模型. 通过求解瞬态热力耦合方程获得结构的热力耦合静力分析结果; 通过响应量在空间和时间域的凝聚积分函数表征结构响应在时域内的最大值, 并以此构建相应的约束和目标函数; 采用伴随法推导约束和目标函数的灵敏度表达式. 通过3个数值算例验证了本方法的有效性. 数值算例结果表明, 在瞬态传热条件下, 本方法能够准确反映瞬态热效应对一体化热防护结构设计结果的影响; 相比于基于稳态热分析的设计结果, 考虑瞬态热效应的设计结果具有更优的性能.Abstract: The integrated thermal protection structure is usually in a severe unstable thermal environment, and the time effect of thermal load, namely transient thermal effect, is obvious. In order to avoid huge calculation consumption of transient thermal analysis, previous optimization design studies of integrated thermal protection structures usually equivalent transient heat transfer to steady-state heat transfer under the same thermal boundary conditions, and take the temperature field of steady-state heat transfer analysis as the design thermal load. However, previous studies have shown that the steady-state heat transfer cannot accurately equivalent the effect of transient heat transfer, and the transient thermal effect has an important influence on the structural design results. In this paper, the optimization design problem of integrated thermal protection structure considering transient thermal effect is studied, and a topology optimization method of integrated thermal protection structure considering transient temperature and stress constraints is established. Based on the solid isotropic material with penalization (SIMP) method, two kinds of topology optimization models for integrated thermal protection structures are constructed: (1) the stiffness design model taking minimizing the structural strain energy as objective function, considering material volume fraction, maximum stress and maximum bottom-face temperature constraints; (2) the strength design model taking minimizing material volume fraction as objective function, considering maximum stress and maximum bottom-face temperature constraints. By solving the transient thermodynamic coupling equation, the thermodynamic coupling static analysis results of the structure are obtained. The maximum value of structural response in time domain is represented by the condensed integral function in space and time domains, which was taken as constraint and objective functions. The sensitivity expressions of objective function and constraint functions are derived by adjoint method. The effectiveness of the proposed method is verified by three numerical results. Numerical examples showed that the proposed method could accurately reflect the influence of transient thermal effects on the design results of integrated thermal protection structures under the condition of transient heat transfer. Compared with the design results based on steady-state thermal analysis, the design results considering transient thermal effects were significantly improved.
-
引言
高超声速飞行器在飞行过程中其表面要承受大量的气动压力和与空气摩擦产生的气动热. 兼具承载和热防护能力的一体化热防护结构成为了高超声速飞行器的理想选择[1-2]. 高超飞行器的飞行时间有限, 但面临着严酷的气动环境. 作用于高超声速飞行器表面的气动热载荷具有加热时间短、作用强度大的特点, 有很强的时间相关性, 为瞬态热载荷. 在瞬态热载荷的作用下, 结构的力学和热学性能会在工作期间随时间不断变化. 在设计过程中准确反应传热过程的瞬态热效应对于结构响应的影响、优化结果达到理想的承载和热防护效果至关重要. 然而现有的一体化热防护结构优化设计研究所采用的通常是将瞬态传热等效为稳态传热的设计方法. 该方法将相同热载荷边界条件下稳态传热分析的温度场作为热弹性响应的设计热载荷. 然而现有研究显示稳态传热无法准确等效瞬态传热的作用效果. 设计热载荷的简化会造成优化结构在实际热环境下无法满足设计要求. 随着武器装备的发展, 高超声速飞行器面临着更加严酷的气动环境和热环境, 对结构的强度、热防护和轻量化方面的需求更加苛刻. 迫切需要建立能够准确反应瞬态传热影响的一体化热防护结构优化方法, 通过更精准的设计来提升结构的性能和降低结构重量.
一体化热防护结构在工作时会同时受机械载荷和热载荷的作用, 是一个典型的热弹性结构. 热载荷产生的温度会以等效热载荷的形式影响热弹性结构的力学性能. 在实际工程中, 一体化热防护结构多应用于高超声速飞行器等航空航天设备. 因此, 瞬态传热成为一体化热防护结构设计工作必须考虑的因素. 在现有的一体化热防护结构设计研究中, 为了使设计结果能够在瞬态传热下达到理想的承载和热防护效果, 优化对象在瞬态传热作用下的相关结构响应通常作为约束条件被引入优化模型中. 例如, Zhao等[3]在利用基于模拟退火算法的多学科优化方法对瞬态热载荷下的一体化热防护结构进行轻量化设计的工作中将瞬态传热下的底面温度、屈曲特征值、von Mises应力和挠度约束作为优化模型的约束条件. Xie等[4]在所建立的以重量最小为目标的一体化热防护结构优化模型中, 将瞬态传热下的内部温度和局部应力作为约束条件. 除此之外, 在基于现有结构的新型一体化热防护结构设计工作中, 设计对象在瞬态传热下的传热和热弹性分析结果是新型结构几何参数调整的重要依据[5-8]. 以上研究表明, 对于一体化热防护结构尺寸和形状优化, 瞬态传热对于设计变量的取值具有重要影响.
拓扑优化是目前公认的工业装备结构创新设计工具[9], 能够在不需要初始设计的条件下, 通过合理设计给定区域内材料的分布, 自动获得创新结构构型. 由于具有比尺寸优化和形状优化更大的设计空间, 拓扑优化常被用于一体化热防护结构创新构型的优化设计. 瞬态传热分析由于在时间域内的迭代计算会产生大量的计算消耗. 因此, 现有的一体化热防护结构拓扑优化研究通常选择将瞬态传热转换为相同热载荷边界条件下的稳态传热来避免高昂计算消耗. 例如, Yang等[10]在所建立的针对热短路效应和热膨胀失配效应的新型综合热防护系统拓扑优化方法中, 将瞬态传热过程中的几个关键时间点的稳态温度场作为设计热载荷. 随后, Yang等[11]在基于最小净换热率和最小应变能原则所建立的降低等效导热率的同时保持结构刚度的一体化热防护结构拓扑优化方法中, 将瞬态传热转换为固定温度边界的稳态传热. Xu等[12]在所建立的以等效导热系数和弹性应变能最小为目标的一体化热防护结构单元拓扑优化方法中, 将瞬态传热过程中的几个关键时刻的温度场作为热屈曲模型的设计热载荷. 以上简化尽管会极大降低计算消耗, 但会造成优化所用热载荷与实际热载荷存在偏差. 原因是在瞬态热载荷的作用下, 一体化热防护结构的温度和热弹性响应会随着时间不断变化[13-22], 特定时间点和稳态分析的温度场无法准确等效瞬态热载荷的作用效果. 现有研究表明, 温度载荷对于一体化热防护结构的优化设计起到了决定性作用[23]. 除此之外, 在瞬态传热的作用下, 热载荷作用的时间效应, 即瞬态热效应, 对结构拓扑有重要影响[24-32]. 在优化过程中将瞬态传热过程转化为稳态温度场会使得优化结构的性能在实际热环境中下降, 无法达到预想的承载和热防护效果. 除此之外, 简化设计载荷会造成较高的安全设计裕度需求, 造成设计结果材料用量过多, 无法满足轻量化设计要求. 对于一体化热防护结构设计问题, 在设计过程中准确反映瞬态传热对于所关注结构响应的影响是实现结构精准设计, 在满足设计约束的条件下使设计结构性能最优的关键.
本文建立了一种考虑瞬态传热的一体化热防护结构拓扑优化方法. 该方法以SIMP法为基础, 通过最大光滑逼近函数在空间和时间域上的积分近似表示瞬态结构响应在整个时间区间内的最大值, 从而在优化过程中将结构响应在整个时域的变化都考虑在内, 实现对瞬态传热下一体化热防护结构的精准设计. 具体算例表明, 相比于以往研究常用的热载荷简化方法, 本文方法能够准确反应瞬态热效应对优化结构的影响. 优化结果不但达到了预期的效果, 而且性能有了明显提高. 本文接下来的结构如下: 第1节和第2节详细介绍考虑瞬态传热的一体化热防护结构拓扑优化设计方法的问题描述、优化模型和敏度推导过程; 第3节通过在3个不同算例中对不同方法的优化结果进行对比, 分析证明本文方法的必要性和有效性; 第4节总结本文的工作并给出最终结论.
1. 优化模型
1.1 问题描述
图1所示, 典型的一体化热防护结构由顶板(TFS)、腹板(web)、隔热层(insulation)和底板(BFS)组成. 其中顶面板受机械载荷P和瞬态热载荷f的作用, 底面板施加有位移边界条件$ \partial {\varOmega _m} $. 腹板将顶板和底板相连, 承担承载功能. 一般情况下, 顶板、底板和腹板均由弹性系数、导热系数和热容系数较高的金属材料制成. 腹板在承载的同时也会将顶板的热量向内传导. 为了阻止过多的热量传递到结构内部, 具有较低导热和热容系数的隔热材料被布置在腹板之间的空隙来承担隔热功能. 由此可见, 腹板是一体化热防护结构的关键. 理想的腹板既能够承载较大的机械载荷, 又能够使尽可能多的热被隔热层吸收和阻隔. 因此, 本文的设计目的是同时考虑机械载荷和瞬态热载荷的作用, 通过优化腹板的构型使一体化热防护结构的承载和热防护性能在给定的温度和应力约束下达到最优.
![]()
图 1 受机−热载荷作用的一体化热防护结构Figure 1. Integrated thermal protection structure under mechanical-thermal load瞬态热传导的控制方程为
$$ \rho c\frac{{\partial T}}{{\partial t}} - \nabla \cdot (k\nabla T) - \rho Q = 0 $$ (1) 其中$ T $表示结构的温度场, $ \rho $表示材料密度, $ c $表示材料的热容, $ t $表示瞬态传热过程的时间变量, $ k $表示材料的导热系数, $ Q $表示物体内部每单位质量产生的热能.
温度场会使受热物体产生初始热应变$ {\varepsilon ^{th}} $, 并以等效热载荷的形式作用结构的静力分析. 瞬态传热下热弹性结构的静力分析控制方程为
$$ \frac{{\partial {\sigma _x}}}{{\partial x}} + \frac{{\partial {\tau _{yx}}}}{{\partial y}} + {f_x} + {f_T} = 0 $$ (2) $$ \frac{{\partial {\sigma _y}}}{{\partial y}} + \frac{{\partial {\tau _{xy}}}}{{\partial x}} + {f_y} + {f_T} = 0 $$ (3) 其中, $ {\sigma _x} $为x方向正应力, $ {\sigma _y} $为y方向正应力, $ {\tau _{xy}} $为切应力, fx和fy分别为x方向和y方向的机械载荷, $ {f_T} $即为t时刻温度场引起的等效热载荷.
具有良好热防护效果的一体化热防护结构可以有效控制热量向系统内部传导. 本文选择将底板的时域最大温度作为设计约束. 在优化过程中, 结构的改变会影响温度的分布. 除此之外, 瞬态热载荷的时间相关性也会造成结构温度场的波动. 因此, 底板的最大温度通常不是一个连续函数. 这会给优化问题的求解带来困难. 为了解决这一问题, 文献[27]通过最大函数近似平滑理论, 将特定区域的时域最大温度近似表达为温度在整个时域和特定空间域内的凝聚积分函数
$$\left.\begin{split} &{T_{{\rm{Max}}}} = \frac{{\displaystyle\int_0^{{t_f}} {{\xi _T}(t){a_T}^{{\xi _T}(t)}{\rm{d}}t} }}{{\displaystyle\int_0^{{t_f}} {{a_T}^{{\xi _T}(t)}{\rm{d}}t} }} \\ &{\xi _T}(t) = \frac{{\displaystyle\int_{{\varOmega _T}} {T(x,t){a_T}^{T(x,t)}{\rm{d}}x} }}{{\displaystyle\int_{{\varOmega _T}} {{a_T}^{T(x,t)}{\rm{d}}x} }} \end{split}\right\} $$ (4) 其中, $ {t_f} $为瞬态热载荷的工作时间, ${\varOmega _T}$为温度控制区域. $ {a_T} $为温度凝聚积分的指数, 当$ {a_T} \to {\text{+}}\infty $时, ${T_{{\rm{Max}}}} \to {\rm{Max}}(T({\varOmega _T},t))$. 为了保证光滑逼近函数的准确性和优化算法的稳定性, 本文将凝聚积分指数设为在优化过程中随迭代过程逐步增加的变化参数.
一体化热防护结构是一个典型的热弹性结构. 已有的研究表明, 由于不能准确反映热弹性结构的变形与强度的关系, 柔顺性不适合作为热弹性结构强度设计的目标. 相比之下, 应变能更适合作为热弹性结构刚度设计的目标[33]. 与区域上的时域最大温度类似, 一体化热防护结构在瞬态传热下的时域最大应变能同样由结构应变能在时域的凝聚积分来近似表示
$$\left. \begin{split} & {\varPhi _{{\rm{Max}}}} = \frac{{\displaystyle\int_0^{{t_f}} {\varPhi (t){a_\varPhi }^{\varPhi (t)}{\rm{d}}t} }}{{\displaystyle\int_0^{{t_f}} {{a_\varPhi }^{\varPhi (t)}{\rm{d}}t} }} \\ & \varPhi (t) = \frac{{\text{1}}}{{\text{2}}}\displaystyle\int_V {{{({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))}^{\rm{T}}}{\boldsymbol{D}}({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))} {\rm{d}}x \end{split}\right\} $$ (5) 其中, ${a_\varPhi }$为应变能凝聚指数, 当${a_\varPhi } \to {\text{+}}\infty$时, ${\varPhi _{{\rm{Max}}}} \to {\rm{Max}}(\varPhi (t))$; ${\boldsymbol{\varepsilon}}$为结构在机械载荷和热载荷作用下的总应变; ${\boldsymbol{D}}$为弹性矩阵.
结构失效主要是由于工作过程中的最大应力超过强度极限引起的. 因此, 在优化设计过程中必须考虑最大应力约束来避免结构失效. 热弹性结构的应力计算公式为
$$ {\boldsymbol{\sigma}} = {\boldsymbol{D}}({\boldsymbol{\varepsilon}} - {{\boldsymbol{\varepsilon}} ^{th}}) $$ (6) 其最大值也可以由时域和空间域内的凝聚积分来表示
$$\left.\begin{split} & {\sigma _{{\rm{Max}}}} = \frac{{\displaystyle\int_0^{{t_f}} {{\xi _s}(t){a_s}^{{\xi _s}(t)}{\rm{d}}t} }}{{\displaystyle\int_0^{{t_f}} {{a_s}^{{\xi _s}(t)}{\rm{d}}t} }} \\ & {\xi _s}(t) = \frac{{\displaystyle\int_V {{\boldsymbol{\sigma}} (x,t){a_s}^{{\boldsymbol{\sigma}} (x,t)}{\rm{d}}x} }}{{\displaystyle\int_V {{a_s}^{{\boldsymbol{\sigma}} (x,t)}{\rm{d}}x} }} \end{split}\right\} $$ (7) 其中, $ {a_s} $为应力凝聚积分指数.
在以往的一体化热防护结构优化设计研究中, 为了避免在时域内迭代求解造成的高昂计算费用, 通常会选择将瞬态传热转化为相同热边界条件下的稳态传热. 在此情况下, 目标函数和约束函数的凝聚积分计算无需在时域上进行. 在稳态热分析中, 一体化热防护结构的底面通常被设为固定温度边界. 无论材料分布如何, 结构底面温度总是一定的, 在此情况下无法设置底面最高温度约束. 一般会用净传热率作为结构热防护能力优化的目标 [10]. 在稳态温度载荷下, 结构应变能、最大温度和净传热效率的表达式如下
$$ {\varPhi ^ * } = {\frac{{\text{1}}}{{\text{2}}}\int_V {{{({\boldsymbol{\varepsilon}} - {{\boldsymbol{\varepsilon}} ^{th}})}^{\rm{T}}}{\boldsymbol{D}}({\boldsymbol{\varepsilon}} - {{\boldsymbol{\varepsilon}} ^{th}}){\rm{d}}x} }$$ (8) $$ \sigma _{{\rm{Max}}}^ * = {\frac{{\displaystyle\int_V {{\boldsymbol{\sigma}} (x){a_s}^{{\boldsymbol{\sigma}} (x)}{\rm{d}}x} }}{{\displaystyle\int_V {{a_s}^{{\boldsymbol{\sigma}} (x)}{\rm{d}}x} }}}$$ (9) $$ Q = \int_V {k\nabla T{\rm{d}}x} $$ (10) 其中, $ \nabla T $为温度梯度场.
1.2 有限元列式
瞬态传热控制方程和热力耦合静力分析控制方程的有限元形式为
$$ {\boldsymbol{C}}\dot {\boldsymbol{T}}(t) + {\boldsymbol{KT}}(t) = {\boldsymbol{F}}(t) $$ (11) $$ {{\boldsymbol{K}}_{^m}}{\boldsymbol{U}}(t) = {{\boldsymbol{P}}_m}{\text{ + }}{{\boldsymbol{P}}_{th}}(t) $$ (12) 其中, ${\boldsymbol{C}}$为整体热容矩阵, ${\boldsymbol{K}}$为整体导热矩阵; ${\boldsymbol{F}}(t)$为施加瞬态热载荷矢量, ${\boldsymbol{T}}(t)$为t时刻的温度场向量, $\dot {\boldsymbol{T}}(t)$为t时刻的温度场向量关于时间的导数. ${{\boldsymbol{K}}_m}$是结构整体刚度矩阵, ${\boldsymbol{U}}(t)$是t时刻的位移向量, ${{\boldsymbol{P}}_m}$是机械载荷向量, ${{\boldsymbol{P}}_{th}}(t)$是由热应变引起的等效温度载荷向量. ${\boldsymbol{C}}$, ${\boldsymbol{K}}$和${{\boldsymbol{K}}_m}$均是由单元热容矩阵和单元导热矩阵组装得到
$$\qquad\qquad\qquad {\boldsymbol{C}}{\text{ = }}\sum\limits_{e = 1}^{{N_e}} {{{\boldsymbol{C}}_e}} $$ (13) $$\qquad\qquad\qquad {\boldsymbol{K}}{\text{ = }}\sum\limits_{e = 1}^{{N_e}} {{{\boldsymbol{K}}_e}} $$ (14) $$\qquad\qquad\qquad {{\boldsymbol{K}}_m}{\text{ = }}\sum\limits_{e = 1}^{{N_e}} {{{\boldsymbol{K}}_{me}}} $$ (15) 其中, $ {N_e} $为结构单元数量, ${{\boldsymbol{C}}_e}$, ${{\boldsymbol{K}}_e}$和${{\boldsymbol{K}}_{me}}$分别为单元热容矩阵、单元导热矩阵和单元刚度矩阵. 其表达式为
$$ {{\boldsymbol{C}}_e} = {c_e}\int_{{\varOmega _e}} {{N_i}} {N_j}{\rm{d}}{\varOmega _e} $$ (16) $$ {{\boldsymbol{K}}_e} = {k_e}\int_{{\varOmega _e}} {\frac{{\partial {N_i}}}{{\partial x}}\frac{{\partial {N_j}}}{{\partial x}} + \frac{{\partial {N_i}}}{{\partial y}}\frac{{\partial {N_j}}}{{\partial y}}} {\rm{d}}{\varOmega _e} $$ (17) $$ {{\boldsymbol{K}}_{me}} = \int_{{\varOmega _e}} {{{\boldsymbol{B}}_e}^{\rm{T}}} {{\boldsymbol{D}}_e}{{\boldsymbol{B}}_e}{\rm{d}}{\varOmega _e} $$ (18) $$ {{\boldsymbol{D}}_e} = {E_e} \cdot {{\boldsymbol{D}}_0} $$ (19) 其中, $ {N_i} $和$ {N_j} $为形函数, ${{\boldsymbol{B}}_e}$为单元应变−位移关系矩阵, ${{\boldsymbol{D}}_e}$为单元弹性矩阵, ${{\boldsymbol{D}}_0}$为单元弹性矩阵的常数项. $ {c_e} $, $ {k_e} $和$ {E_e} $分别为单元热容系数、单元导热系数和单元弹性模量. 本文采用的最常用的SIMP法来实现材料分布设计, 即每个有限单元内的材料密度$ {\rho _{{e}}} $作为设计变量, 通过优化设计域内材料密度的分布方式$ \chi (\rho ) $使结构性能达到最优. 单元材料属性参数热传导热容$ {c_e} $、系数$ {k_e} $和弹性模量$ {E_e} $的有限元插值模型分别表示为
$$\qquad\qquad\quad {c_e} = {c_2} + \rho _e^{{p_2}}({c_1} - {c_2}) $$ (20) $$\qquad\qquad\quad {k_e} = {k_2} + \rho _e^{{p_1}}({k_1} - {k_2}) $$ (21) $$\qquad\qquad\quad {E_e} = E_2^{} + \rho _e^{{p_3}}(E_1^{} - E_2^{}) $$ (22) 其中$ {\text{0}} \leqslant {\rho _e} \leqslant {\text{1}} $, 当$ {\rho _e}{\text{ = 1}} $时表示有限单元内仅有材料1, 当$ {\rho _e}{\text{ = 0}} $时表示有限元中单元内仅有材料2, 而$ {\text{0 < }}{\rho _e}{\text{ < 1}} $则假设为材料1和材料2的复合体. 下标1和2表示材料编号. $ {k_{\text{1}}} $, $ {c_{\text{1}}} $, $ {E_{\text{1}}} $, $ {k_2} $, $ {c_2} $ 和$ E_2^{} $分别为材料1的导热系数、材料1的热容系数、材料1的弹性模量、材料2的热容、材料2的导热系数和材料2的弹性模量. $ {p_1} $, $ {p_2} $, $ {p_3} $为惩罚因子, 其作用是使设计变量取值趋向于上限值或下限值. 在本文中$ {p_1} $, $ {p_2} $, $ {p_3} $均取值为3.
由温度场引起的等效温度载荷向量的表达式为
$$ {{\boldsymbol{P}}_{th}}{\text{(}}t{\text{) = }}{{\boldsymbol{S}}_P}({\boldsymbol{T}}(t) - {{\boldsymbol{T}}_0}) $$ (23) $$ {{\boldsymbol{S}}_p}{\text{ = }}\sum\limits_e {\int_{{\varOmega _e}} {{\boldsymbol{B}}_e^{\rm{T}}} {{\boldsymbol{D}}_e}} {{\boldsymbol{B}}_{the}}{\rm{d}}{\varOmega _e} $$ (24) $$ {{\boldsymbol{B}}_{the}} = \alpha [\begin{array}{*{20}{c}} 1&1&{0{]^{\rm{T}}}} \end{array} $$ (25) 其中, ${{\boldsymbol{T}}_0}$为初始温度场, $ \alpha $为单元热膨胀系数.
在分析过程中, 首先进行方程(11)瞬态传热分析得到结构的温度场. 时域和空域内的凝聚积分表示的区域温度控制函数可表达为
$$ \left.\begin{split} & {T_{{\rm{Max}}}} = \frac{{\displaystyle\int_{\text{0}}^{{t_f}} {{\xi _T}(t){a_T}^{{\xi _T}(t)}{\rm{d}}t} }}{{\displaystyle\int_{\text{0}}^{{t_f}} {{a_T}^{{\xi _T}(t)}{\rm{d}}t} }} \\ & {\xi _T}(t) = \frac{{\displaystyle\sum\limits_{{n_i} = 1}^{ {N_c}} {{T_{ni}}(t){a_T}^{T({n_i},t)}} }}{{\displaystyle\sum\limits_{{n_i} = 1}^{ {N_c}} {{a_T}^{{T_{ni}}(t)}} }} \end{split} \right\}$$ (26) 其中$ {N_c} $为温度控制区域内的单元节点数. 瞬态传热分析后, 将由温度场计算所得的等效温度载荷代入方程(12)中进行热力耦合分析计算, 最终得到结构的位移场与应变场. 结构应变能在时域内的凝聚积分可表达为
$$\left. \begin{split} & {\varPhi _{{\rm{Max}}}} = \frac{{\displaystyle\int_0^{{t_f}} {\varPhi (t){a_\varPhi }^{\varPhi (t)}{\rm{d}}t} }}{{\displaystyle\int_0^{{t_f}} {{a_\varPhi }^{\varPhi (t)}{\rm{d}}t} }} \\ & \varPhi (t) = {({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))^{\rm{T}}}{\boldsymbol{D}}({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t)) \end{split} \right\}$$ (27) 其中, ${\boldsymbol{\varepsilon}} (t)$为t时刻真实应变场, ${{\boldsymbol{\varepsilon}} ^{th}}(t)$为t时刻变温度场引起的单元初始应变场, 它们与t时刻的位移场向量和温度场向量的关系为
$$ {\boldsymbol{\varepsilon}} (t) = {\boldsymbol{BU}}(t),\qquad {{\boldsymbol{B}}{\text{ = }}\sum\limits_{e = 1}^{{N_e}} {{{\boldsymbol{B}}_e}} } $$ (28) $$ {{\boldsymbol{\varepsilon}} ^{th}}(t) = {{\boldsymbol{B}}_{th}}({\boldsymbol{T}}(t) - {{\boldsymbol{T}}_0}),\qquad {{{\boldsymbol{B}}_{th}}{\text{ = }}\sum\limits_{e = 1}^{{N_e}} {{{\boldsymbol{B}}_{the}}} } $$ (29) ${\boldsymbol{D}}$为总弹性矩阵, 由单元弹性矩阵组装得到
$$ {\boldsymbol{D}}{\text{ = }}\sum\limits_{e = 1}^{{N_e}} {{{\boldsymbol{D}}_e}} $$ (30) t时刻的单元应力计算公式为
$$ {{\boldsymbol{\sigma}} _e}(t) = {{\boldsymbol{D}}_e}({{\boldsymbol{\varepsilon}} _e}(t) - {{\boldsymbol{\varepsilon}} _e}^{th}(t)) $$ (31) $$ {{\boldsymbol{\sigma}} _e}(t) = [{{\boldsymbol{\sigma}} _x}(t);{{\boldsymbol{\sigma}} _y}(t);{{\boldsymbol{\tau}} _{xy}}(t)] $$ (32) von Mises应力可由单元应力计算得出
$$ {{\boldsymbol{\sigma}}^M _e}(t) = {\left( {{{\boldsymbol{\sigma}} ^{\rm{T}}_e}(t){\boldsymbol{V}}{{\boldsymbol{\sigma}} _e}(t)} \right)^{\tfrac{1}{2}}} $$ (33) 其中${\boldsymbol{V}}$为平面系数矩阵
$${\boldsymbol{V}}=\left[ \begin{matrix} 1 & -0.5 & 0 \\ -0.5 & 1 & 0 \\ 0 & 0 & 3 \\ \end{matrix} \right] $$ (34) 为了避免应力约束问题中的奇异现象, 要求应力约束松弛[34]
$$ {{\boldsymbol{\sigma}} ^{MR}_e}(t){\text{ = }}{\boldsymbol{x}}_e^q{{\boldsymbol{\sigma}}^M _e}(t) $$ (35) 其中, ${\boldsymbol{x}}_e^q$为设计变量向量, $ q $为释放参数. 结构 von Mises应力在时域和空间域内的凝聚积分可表示为
$$\left.\begin{split} & \sigma _{{\rm{Max}}}^M = \frac{{\displaystyle\int_{\text{0}}^{{t_f}} {{\xi _s}(t){a_s}^{{\xi _s}(t)}{\rm{d}}t} }}{{\displaystyle\int_{\text{0}}^{{t_f}} {{a_s}^{{\xi _s}(t)}{\rm{d}}t} }} \\ & {\xi _s}(t) = \frac{{\displaystyle\sum\limits_{e = 1}^{{N_e}} {{\boldsymbol{\sigma}} _e^{MR}(t){a_s}^{{\boldsymbol{\sigma}} _e^{MR}(t)}} }}{{\displaystyle\sum\limits_{e = 1}^{ {N_e}} {{a_s}^{{\boldsymbol{\sigma}} _e^{MR}(t)}} }} \end{split}\right\} $$ (36) 2. 拓扑优化方法
2.1 优化模型
这里主要构建两类优化模型, 一是以结构刚度为目标, 考虑材料体积分数、最大应力和底面温度约束; 二是结构轻量化为目标, 考虑最大应力和底面温度约束. 两类优化模型均由两种方法构建, 一种是本文提出的考虑瞬态传热影响的方法, 另一种是以往常用的将瞬态传热转换为稳态传热的方法.
一体化热防护结构同时受热载荷和机械载荷作用, 是典型的热弹性结构. 已有的研究表明, 对于热弹性结构, 应变能更适合作为热弹性结构刚度设计的目标. 因此, 对于刚度设计, 我们选择将应变能作为结构刚度设计的目标函数. 同时为了避免工作过程中的结构失效, 必须在优化设计时考虑最大应力约束. 为了使优化结构能够达到理想的热防护效果, 必须保证在结构底面的最大温度在工作过程中始终低于预设的可接受最大温度. 由此, 以瞬态热载荷作用下一体化热防护结构为设计对象, 采用本文方法所构建最小化结构时域最大应变能同时考虑最大应力约束、底面时域最大温度和材料体积分数为约束拓扑优化模型如下
$$\left.\begin{split} & {\rm{find}}:\qquad {\chi (\rho )} \\ & {\rm{minimize}}:\qquad {Obj = {\varPhi _{{\rm{Max}}}}} \\ & \begin{aligned} {{\rm{subject}}}\;{{\rm{to:}}} \qquad &{\frac{1}{A}\sum\limits_{e = 1}^{{N_e}} {{\rho _e}{A_e}} \leqslant \phi } \\ &{{\boldsymbol{C}}{{\dot {\boldsymbol{T}}}}(t) + {\boldsymbol{KT}}(t) = {\boldsymbol{F}}(t)} \\ & {{{\boldsymbol{K}}_m}{\boldsymbol{U}}(t) = {{\boldsymbol{P}}_m}{\text{ + }}{{\boldsymbol{P}}_{th}}(t)} \\ & \sigma _{{\rm{Max}}}^M \leqslant {{{\hat \sigma }_0}} \\ &{T_{{\rm{Max}}}} \leqslant {{{\hat T}_0}}\end{aligned} \end{split}\right\}$$ (37) 其中, A为设计域面积, Ae为设计域内的单元面积, $ \phi $为所设定的体积分数上限, $ {\hat \sigma _{\text{0}}} $为所设定的应力约束上限, $ {\hat T_{\text{0}}} $为所设定的温度约束上限. 以上优化模型以应力在时域和空间域的凝聚积分、底面温度在时域的凝聚积分和材料体积分数为约束函数, 以结构应变能在时域的凝聚积分为目标函数. 这样, 优化模型可以将所关注的结构响应在整个瞬态传热过程中的变化都考虑在内, 实现对瞬态传热下一体化热防护结构的精准设计.
对于轻量化设计, 以结构中金属材料的体积分数作为目标函数, 同时考虑最大应力约束和底面最大温度约束保证结构的承载和热防护能力. 由于本文的轻量化设计的目的是在满足应力和温度约束的基础上最小化结构金属材料的体积, 因此所建立的轻量化设计模型并未考虑结构的刚度. 由此, 以瞬态热载荷作用下一体化热防护结构为设计对象, 采用本文提出的方法构建以最小化结构金属材料体积分数为目标, 同时考虑最大应力约束、底面时域最大温度约束的拓扑优化模型如下
$$\left. \begin{split} & {\rm{find}}:\qquad {\chi (\rho )} \\ & {\rm{minimize}}:\qquad {Obj = \frac{1}{A}\sum\limits_{e = 1}^{{N_e}} {{\rho _e}{A_e}} } \\ & \begin{aligned} {{\rm{subject}}}\;{{\rm{to}}:} \qquad &{{\boldsymbol{C}}\dot {\boldsymbol{T}}(t) + {\boldsymbol{KT}}(t) = {\boldsymbol{F}}(t)} \\ & {{{\boldsymbol{K}}_m}{\boldsymbol{U}}(t) = {{\boldsymbol{P}}_m}{\text{ + }}{{\boldsymbol{P}}_{th}}(t)} \\ & \sigma _{{\rm{Max}}}^M \leqslant {{{\hat \sigma }_0}} \\ & {T_{{\rm{Max}}}} \leqslant {{{\hat T}_0}}\end{aligned} \end{split}\right\} $$ (38) 以上优化模型以应力在时域和空间域的凝聚积分、底面温度在时域的凝聚积分为约束函数, 以结构金属材料的体积分数为目标函数. 这样, 优化模型可以在保证结构的承载和热防护性能的条件下, 最大程度减轻一体化热防护结构的重量.
在以往的基于稳态热分析的一体化热防护结构拓扑优化研究中, 由于一体化热防护结构的底面通常被设为固定温度边界, 无法设置底面最高温度约束, 一般采用将所关注的承载和热传导目标的权重和作为目标函数的方法来实现承载和热防护的协同优化设计. 本文所建立的基于稳态传热分析的拓扑优化模型采用的是上述承载和热防护协同优化设计方法. 以稳态传热下一体化承载−热防护结构为设计对象, 考虑最大应力约束的刚度与热防护协同设计拓扑优化模型和轻量化与热防护协同设计拓扑优化模型分别为:
(1) 刚度与热防护协同设计拓扑优化模型
$$\left.\begin{split} & {\rm{find}}:\qquad {\chi (\rho )} \\ & {\rm{minimize}}:\qquad {Obj = w \cdot \frac{{{\varPhi ^ * }}}{{\varPhi _{{\rm{Max}}}^0}}} + (1 - w) \cdot \frac{Q}{{{Q^{\text{0}}}}} \\ & \begin{aligned} {{\rm{subject}}}\;{{\rm{to}}:}\qquad &{\frac{1}{A}\sum\limits_{e = 1}^{{N_e}} {{\rho _e}{A_e}} \leqslant \phi } \\ & {{\boldsymbol{KT}} = {\boldsymbol{F}}} \\ & {{{\boldsymbol{K}}_m}{\boldsymbol{U}} = {{\boldsymbol{P}}_m}{\text{ + }}{{\boldsymbol{P}}_{th}}} \\ & \sigma _{{\rm{Max}}}^ * \leqslant {{{\hat \sigma }_0}}\end{aligned} \end{split}\right\} $$ (39) (2) 轻量化与热防护协同设计拓扑优化模型
$$ \left.\begin{split} & {\rm{find}}:\qquad {\chi (\rho )} \\ & {\rm{minimize}}:\qquad {Obj = w \cdot \frac{1}{A}\sum\limits_{e = 1}^{{N_e}} {{\rho _e}{A_e}} } + (1 - w) \cdot \frac{Q}{{{Q^{\text{0}}}}} \\ & \begin{aligned} {{\rm{subject}}}\;{{\rm{to}}:} \qquad &{{\boldsymbol{KT}} = {\boldsymbol{F}}} \\ & {{{\boldsymbol{K}}_m}{\boldsymbol{U}} = {{\boldsymbol{P}}_m}{\text{ + }}{{\boldsymbol{P}}_{th}}} \\ & \sigma _{{\rm{Max}}}^ * \leqslant {{{\hat \sigma }_0}} \end{aligned}\end{split}\right\} $$ (40) 两种模型均以应力在空间和时间域的凝聚积分为约束. 刚度设计模型以结构应变能和净传热效率的加权之和为目标函数, 轻量化设计模型以金属材料体积分数和净传热效率的加权之和为目标函数. 因为应变能和净传热效率的数量级可能差别较大, 因此需要进行无量纲化处理. 采用的方法是将它们分别除以初始设计的应变能和净传热效率. $\varPhi _{{\rm{Max}}}^{\text{0}}$和$ {Q^{\text{0}}} $分别为在优化过程中初始迭代步的最大应变能和净传热效率. $ w $为权重系数, 当$ w{\text{ = 1}} $时, 为仅最小关注目标函数, 当$ w{\text{ = 0}} $时, 为仅最小化净传热效率, 当$ w \in (0,1) $时, 为最小化二者的加权之和. 为了平衡优化结构的承载和热防护功能, 本文权重系数的取值为$ w{\text{ = 0}}{\text{.5}} $. 其中, $ A $为设计域的总面积, $ {A_e} $为单元面积, $ \phi $为体积分数上限值.
2.2 敏度分析
本文采用基于梯度的优化求解算法对优化模型进行求解. 因此, 求解能够表达设计变量目标和约束关于设计变量梯度信息的灵敏度分析是非常重要的. 本文选择运用伴随法求解目标函数和约束函数关于设计变量的敏度解析式. 结构响应在时域的凝聚积分的一般形式为
$$ f = {{\int_0^{{t_f}} {\xi (t){a^{\xi (t)}}{\rm{d}}t} } \mathord{\left/ {\vphantom {{\int_0^{{t_f}} {\xi (t){a^{\xi (t)}}{\rm{d}}t} } {\int_0^{{t_f}} {{a^{\xi (t)}}{\rm{d}}t} }}} \right. } {\int_0^{{t_f}} {{a^{\xi (t)}}{\rm{d}}t} }} $$ (41) 其中, $ a $为凝聚指数, $ \xi (t) $为t时刻的响应. 当$ \xi (t) $为整体响应时, 其表达式为
$$ \xi (t) = \sum\limits_{ei = 1}^{N} {{\xi _{ei}}(t)} $$ (42) 其中, $ {\xi _{ei}}(t) $为t时刻的单个局部响应, $ N $为单个局部响应的总数. 当$ \xi (t) $为局部响应的最大值时, 其表达式为
$$ \xi (t) = {{\sum\limits_{ei = 1}^{N} {{\xi _{ei}}(t){a^{{\xi _{ei}}(t)}}} } \mathord{\left/ {\vphantom {{\sum\limits_{ei = 1}^{ N} {{\xi _{ei}}(t){a^{{\xi _{ei}}(t)}}} } {\sum\limits_{ei = 1}^{ N} {{a^{{\xi _{ei}}(t)}}} }}} \right. } {\sum\limits_{ei = 1}^{ N} {{a^{{\xi _{ei}}(t)}}} }} $$ (43) 引入拉格朗日乘子$ {{\boldsymbol{\lambda}} _1} $和$ {{\boldsymbol{\lambda}} _2} $后, 结构响应在时域的凝聚积分可表达为
$$\begin{split} & L{\text{ = }}f + \int_0^{{t_f}} {{\boldsymbol{\lambda}} _{\text{1}}^{\rm{T}}} \left( {{{\boldsymbol{K}}_m}{{{\boldsymbol{U}}}}(t) - {\boldsymbol{P}} - {{\boldsymbol{P}}_{th}}(t)} \right){\rm{d}}t + \\ & \qquad \int_0^{{t_f}} {{\boldsymbol{\lambda}} _2^{\rm{T}}} \left( {{\boldsymbol{C}}\dot {\boldsymbol{T}}(t) + {\boldsymbol{KT}}(t) - {\boldsymbol{F}}(t)} \right){\rm{d}}t \end{split} $$ (44) 在此情况下, 结构响应在时域的凝聚积分关于设计变量的显式解析式为
$$ \begin{split} & \frac{{\partial L}}{{\partial {\rho _e}}}{\text{ = }}\frac{{\partial f}}{{\partial {\rho _e}}} + \int_0^{{t_f}} {{\boldsymbol{\lambda}} _1^{\rm{T}}} \left(\frac{{\partial {{\boldsymbol{K}}_m}}}{{\partial {\rho _e}}}{\boldsymbol{U}}(t) + \right. \\ & \qquad \left.{{\boldsymbol{K}}_m}\frac{{\partial {\boldsymbol{U}}(t)}}{{\partial {\rho _e}}} - \frac{{\partial {{\boldsymbol{S}}_p}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t) - {{\boldsymbol{S}}_p}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}\right){\rm{d}}t + \\ & \qquad \int_0^{{t_f}} {{\boldsymbol{\lambda}} _2^{\rm{T}}} \left(\frac{{\partial {{{\boldsymbol{C}}}}}}{{\partial {\rho _e}}}\dot {\boldsymbol{T}}(t) + {{{\boldsymbol{C}}}}\frac{{\partial \dot {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}} + \frac{{\partial {\boldsymbol{K}}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t) + {\boldsymbol{K}}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}\right){\rm{d}}t \end{split} $$ (45) $$ \frac{{\partial f}}{{\partial {\rho _e}}}{\text{ = }}\int_0^{{t_f}} {\beta (t)} \frac{{\partial \xi (t)}}{{\partial {\rho _e}}}{\rm{d}}t $$ (46) $$ \beta (t) = \frac{1}{A}{a^{\xi (t)}}\left( {1 + \xi (t)\ln a} \right) - \frac{B}{{{A^2}}}{a^{\xi (t)}}\ln a $$ (47) $$ A{\text{ = }}\int_0^{{t_f}} {{a^{\xi (t)}}} {\rm{d}}t,\quad {B{\text{=}}} \int_0^{{t_f}} {\xi (t){a^{\xi (t)}}} {\rm{d}}t $$ (48) 根据分步积分法, $\displaystyle\int_0^{{t_f}} {{\boldsymbol{\lambda}} _2^{\rm{T}}} \dfrac{{\partial {{{\boldsymbol{C}}}}}}{{\partial {\rho _e}}}\dot {\boldsymbol{T}}(t){\rm{d}}t$可表征为
$$\begin{split} & \int_0^{{t_f}} {{\boldsymbol{\lambda}} _2^{\rm{T}}} {\boldsymbol{C}}\frac{{\partial \dot {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}{\rm{d}}t = {\boldsymbol{\lambda}} _2^{\rm{T}}({t_f}){\boldsymbol{C}}\frac{{\partial {\boldsymbol{T}}({t_f})}}{{\partial {\rho _e}}} - \\ & \qquad \int_0^{{t_f}} {\dot {\boldsymbol{\lambda}} _2^{\rm{T}}} (t){\boldsymbol{C}}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}{\rm{d}}t \end{split} $$ (49) 则方程(45)可表达为
$$\begin{split} & \frac{{\partial L}}{{\partial {\rho _e}}}{\text{ = }}\frac{{\partial f}}{{\partial {\rho _e}}} + \int_0^{{t_f}} {{\boldsymbol{\lambda}} _1^{\rm{T}}} \left(\frac{{\partial {{\boldsymbol{K}}_m}}}{{\partial {\rho _e}}}{\boldsymbol{U}}(t) +\right. \\ & \qquad \left.{{\boldsymbol{K}}_m}\frac{{\partial {\boldsymbol{U}}(t)}}{{\partial {\rho _e}}} - \frac{{\partial {{\boldsymbol{S}}_p}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t) - {{\boldsymbol{S}}_p}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}\right){\rm{d}}t + \\ & \qquad { {\boldsymbol{\lambda}} _2^{\rm{T}}({t_f}){{{\boldsymbol{C}}}}\frac{{\partial {\boldsymbol{T}}({t_f})}}{{\partial {\rho _e}}}} + \int_0^{{t_f}} {{\boldsymbol{\lambda}} _2^{\rm{T}}} \left(\frac{{\partial {{{\boldsymbol{C}}}}}}{{\partial {\rho _e}}}\dot {\boldsymbol{T}}(t) - \right. \\ & \qquad \left.\dot {\boldsymbol{\lambda}} _2^{\rm{T}}(t){{{\boldsymbol{C}}}}\frac{{\partial \dot {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}} + \frac{{\partial {\boldsymbol{K}}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t) + {\boldsymbol{K}}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}\right){\rm{d}}t \end{split} $$ (50) 当$ \xi (t) $为整体响应时, t时刻结构响应关于设计变量的显式解析式$\dfrac{{\partial \xi (t)}}{{\partial {\rho _e}}}$可表达为
$$ \frac{{\partial \xi (t)}}{{\partial {\rho _e}}} = \sum\limits_{ei = 1}^{ N} {\frac{{\partial {\xi _{ei}}(t)}}{{\partial {\rho _e}}}} $$ (51) 当$ \xi (t) $为局部响应最大值时, t时刻结构响应关于设计变量的显式解析式$\dfrac{{\partial \xi (t)}}{{\partial {\rho _e}}}$可表达为
$$ \frac{{\partial \xi (t)}}{{\partial {\rho _e}}} = \sum\limits_{ei = 1}^N {{\chi _{ei}}} {{{{{\boldsymbol{\varLambda}} }}_{ei}}} \frac{{\partial {{{\boldsymbol{\xi}} }}(t)}}{{\partial {\rho _e}}} = { {\textit{χ}}} {\boldsymbol{ \varLambda}} \frac{{\partial {{{\boldsymbol{\xi}} }}(t)}}{{\partial {\rho _e}}} $$ (52) ${{{\boldsymbol{\xi}} }}(t)$为t时刻结构响应向量. 在方程(52)中
$$ \begin{split} & {\chi _{ei}} = {\left(\sum\limits_{ei = 1}^N {{a^{{\xi _{ei}}(t)}}} \right)^{ - 1}}{a^{{\xi _{ei}}(t)}}\left( {1 + {\xi _{ei}}(t)\ln a} \right) - \\ & \qquad {{{\left(\sum\limits_{ei = 1}^N {{a^{{\xi _{ei}}(t)}}} \right)}^{ - 2}}\left(\sum\limits_{ei = 1}^N {{\xi _{ei}}(t)} {a^{{\xi _{ei}}(t)}}\right){a^{{\xi _{ei}}(t)}}\ln a} \end{split} $$ (53) $$ {{{{{\boldsymbol{\varLambda}} }}_{ei}}} = [\begin{array}{*{20}{c}} 0&0& \cdots &{\begin{array}{*{20}{c}} {\dfrac{{{\xi _{ei}}(t)}}{{\left| {{\xi _{ei}}(t)} \right|}}}& \cdots &0&0 \end{array}} \end{array}] $$ (54) 其中, $\dfrac{{{\xi _{ei}}(t)}}{{\left| {{\xi _{ei}}(t)} \right|}}$为${{{{{\boldsymbol{\varLambda}} }}_{ei}}}$ 的第$ ei $个元素. 随后, 将计算所得$\dfrac{{\partial \xi (t)}}{{\partial {\rho _e}}}$代入方程(50), 通过计算拉格朗日乘子的值即可求出敏度值.
结构应变能为整体响应, $ t $时刻结构应变能关于设计变量的显式解析式为
$$ \begin{split} & \frac{{\partial \varPhi (t)}}{{\partial {\rho _e}}} = {({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))^{\rm{T}}}{\boldsymbol{D}}\left(\frac{{\partial {\boldsymbol{\varepsilon}} (t)}}{{\partial {\rho _e}}} - \frac{{\partial {{\boldsymbol{\varepsilon}} ^{th}}(t)}}{{\partial {\rho _e}}}\right)+ \\ &\qquad {({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))^{\rm{T}}}\frac{{\partial {\boldsymbol{D}}}}{{\partial {\rho _e}}}({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))\end{split} $$ (55) 应变向量与总体弹性矩阵关于设计变量的显式解析式为
$$ \frac{{\partial {\boldsymbol{\varepsilon}} (t)}}{{\partial {\rho _e}}} = {\boldsymbol{B}}\frac{{\partial {\boldsymbol{U}}(t)}}{{\partial {\rho _e}}},\begin{array}{*{20}{c}} {}&{} \end{array}\frac{{\partial {{\boldsymbol{\varepsilon}} ^{th}}(t)}}{{\partial {\rho _e}}} = {{\boldsymbol{B}}_{th}}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}} $$ (56) 由此可得结构应变能的凝聚积分关于设计变量的显式解析式为
$$ \begin{split} & \frac{{\partial {L_\varPhi }(x)}}{{\partial {\rho _e}}} = \int_0^{{t_f}} {[2 {\beta _\varPhi }(t) \cdot {{({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))}^{\rm{T}}}{\boldsymbol{D}}} {\boldsymbol{B}} + \\ & \qquad {\boldsymbol{\lambda}} _{\varPhi 1}^{\rm{T}}(t){{\boldsymbol{K}}_m}]\frac{{\partial {\boldsymbol{U}}(t)}}{{\partial {\rho _e}}}{\rm{d}}t + \\ & \qquad \int_0^{{t_f}} {{\beta _\varPhi }(t){{\boldsymbol{\varepsilon}} ^{\rm{T}}}(t)\frac{{\partial {\boldsymbol{D}}}}{{\partial {\rho _e}}}{\boldsymbol{\varepsilon}} } (t) + {\boldsymbol{\lambda}} _{\varPhi 1}^{\rm{T}}(t)\left[\frac{{\partial {{\boldsymbol{K}}_m}}}{{\partial {\rho _e}}}{\boldsymbol{U}}(t)- \right. \\ & \qquad\left. \frac{{\partial {{\boldsymbol{S}}_p}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t)\right]{\rm{d}}t{\text{ + }}\int_{\text{0}}^{{t_f}} {[ - 2 {\beta _\varPhi }(t) \cdot } {({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t))^{\rm{T}}}{\boldsymbol{D}}{{\boldsymbol{B}}_{th}} - \\ & \qquad {\boldsymbol{\lambda}} _{\varPhi 1}^{\rm{T}}(t){{\boldsymbol{S}}_p} - \dot {\boldsymbol{\lambda}} _{\varPhi 2}^{\rm{T}}{\boldsymbol{C}} + {\boldsymbol{\lambda}} _{\varPhi 2}^{\rm{T}}{\boldsymbol{K}}]\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}{\rm{d}}t+ \\ & \qquad \int_0^{{t_f}} {{\boldsymbol{\lambda}} _{\varPhi 2}^{\rm{T}}(\frac{{\partial {\boldsymbol{C}}}}{{\partial {\rho _e}}}\dot {\boldsymbol{T}}(t) + \frac{{\partial {\boldsymbol{K}}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t)){\rm{d}}t}\\[-12pt] \end{split} $$ (57) 随后, 通过方程(58)和式(59)计算拉格朗日乘子的值
$$ {{\boldsymbol{K}}_m}{{\boldsymbol{\lambda}} _{\varPhi 1}}(t) = - 2{\beta _\varPhi }(t) \cdot {{\boldsymbol{B}}^{\rm{T}}}{\boldsymbol{D}}_{}^{\rm{T}}({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t)) $$ (58) $$\begin{split} & - {\boldsymbol{C}}{{\dot {\boldsymbol{\lambda}} }_{\varPhi 2}}(t) + {\boldsymbol{K}}{{\boldsymbol{\lambda}} _{\varPhi 2}}(t) = {\boldsymbol{S}}_p^{\rm{T}}{{\boldsymbol{\lambda}} _{\varPhi 1}}(t){\text{ + }} \\ &\qquad 2{\beta _\varPhi }(t) \cdot {{\boldsymbol{B}}^{\rm{T}}}{\boldsymbol{D}}_{}^{\rm{T}}({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t)),\quad {{{\boldsymbol{\lambda}} _{\varPhi 2}}({t_f})} = {\boldsymbol{0}} \end{split} $$ (59) 将计算所得拉格朗日乘子代入方程(57)中即可得到目标函数关于设计变量的敏度分析值.
结构最大应力为局部响应最大值, 因此采用局部最大响应的敏度求解方法计算t时刻结构最大应力凝聚积分关于设计变量的显式解析式为
$$ \begin{split} & \frac{{\partial {\xi _S}(t)}}{{\partial {\rho _e}}} = \sum\limits_{ni = 1}^{{N_e}} {{\chi _{ni}}} {{{{{\boldsymbol{\varLambda}} }}_{ni}}} \frac{{\partial {{\boldsymbol{\sigma}} ^{MR}}(t)}}{{\partial {\rho _e}}} = \\ &\qquad {{{ {\textit{χ}}}} _s}{ {\boldsymbol{\varLambda}} _s}\left[p \cdot x_e^{p - 1}{{\boldsymbol{\sigma}} ^M}(t) + x_e^p\left(\frac{{\partial {{\boldsymbol{\sigma}} ^M}(t)}}{{\partial {\boldsymbol{\sigma}} (t)}}\right)\frac{{\partial {\boldsymbol{\sigma}} (t)}}{{\partial {\rho _e}}}\right] \end{split} $$ (60) 其中, ${{\boldsymbol{\sigma}} ^M}(t)$和${\boldsymbol{\sigma }}(t)$分别为t时刻的von Mises应力向量和单元应力向量, $\left(\dfrac{{\partial {{\boldsymbol{\sigma}} ^M}(t)}}{{\partial {\boldsymbol{\sigma}} (t)}}\right)$为von Mises应力向量与单元应力向量的转换关系. 单元应力向量关于设计变量的敏度解析式为
$$ \frac{{\partial {\boldsymbol{\sigma}} (t)}}{{\partial {\rho _e}}} = \frac{{\partial {\boldsymbol{D}}}}{{\partial {\rho _e}}}({\boldsymbol{\varepsilon}} (t) - {{\boldsymbol{\varepsilon}} ^{th}}(t)) + {\boldsymbol{D}}\left({\boldsymbol{B}}\frac{{\partial {\boldsymbol{U}}(t)}}{{\partial {\rho _e}}} - {{\boldsymbol{B}}_{th}}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}\right) $$ (61) 由此可得最大应力的凝聚积分函数关于设计变量的显式解析式为
$$ \begin{split} & \frac{{\partial {L_s}(x)}}{{\partial {\rho _e}}} = \int_0^{{t_f}} {{\beta _s}(t)} { {{ {\textit{χ}}}}_s}{ {\boldsymbol{\varLambda}} _s}\Biggr[p \cdot x_e^{p - 1}{{\boldsymbol{\sigma}} ^M}(t)+ \\ & \qquad x_e^p\left(\frac{{\partial {{\boldsymbol{\sigma}} ^M}(t)}}{{\partial {\boldsymbol{\sigma}} (t)}}\right)\frac{{\partial {\boldsymbol{D}}}}{{\partial {\rho _e}}}({\boldsymbol{\varepsilon}} - {{\boldsymbol{\varepsilon}} ^{th}})\Biggr]{\rm{d}}t+ \\ & \qquad \int_{\text{0}}^{{t_f}} \left[ {\beta _s}(t){{{ {\textit{χ}}}} _s}{ {\boldsymbol{\varLambda}} _s}x_e^p\left(\frac{{\partial {{\boldsymbol{\sigma}} ^M}(t)}}{{\partial {\boldsymbol{\sigma}} (t)}}\right){\boldsymbol{DB}}+ {\boldsymbol{\lambda}} _{s{\text{1}}}^{\rm{T}}{{\boldsymbol{K}}_m}\right]\frac{{\partial {\boldsymbol{U}}(t)}}{{\partial {\rho _e}}}{\rm{d}}t {\text{ + }} \\ &\qquad \int_{\text{0}}^{{t_f}} {\Biggr[{\beta _s}(t){{ {{ {\textit{χ}}}} }_s}{{ {\boldsymbol{\varLambda}} }_s}x_e^p\left(\frac{{\partial {{\boldsymbol{\sigma}} ^M}(t)}}{{\partial {\boldsymbol{\sigma}} (t)}}\right){\boldsymbol{D}}{{\boldsymbol{B}}_{th}}}- \\ & \qquad {\boldsymbol{\lambda}} _{s{\text{1}}}^{\rm{T}}{{\boldsymbol{S}}_q} - \dot {\boldsymbol{\lambda}} _{s{\text{2}}}^{\rm{T}}{\boldsymbol{C}} + {\boldsymbol{\lambda}} _{s{\text{2}}}^{\rm{T}}{\boldsymbol{K}}\Biggr]\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}{\rm{d}}t{\text{ + }} \\ & \qquad \int_{\text{0}}^{{t_f}}\Biggr[ {{\boldsymbol{\lambda}} _{s{\text{1}}}^{\rm{T}}\left(\frac{{\partial {{\boldsymbol{K}}_m}}}{{\partial {\rho _e}}}{\boldsymbol{U}}(t) - \frac{{\partial {{\boldsymbol{S}}_q}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t)\right)}{\text{ + }} \\ & \qquad {\boldsymbol{\lambda}} _{s{\text{2}}}^{\rm{T}}\left(\frac{{\partial {\boldsymbol{C}}}}{{\partial {\rho _e}}}\dot {\boldsymbol{T}}(t){\text{ + }}\frac{{\partial {\boldsymbol{K}}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t)\right)\Biggr]{\rm{d}}t\\[-12pt] \end{split} $$ (62) 拉格朗日乘子由以下方程计算得出
$$ {{\boldsymbol{K}}_m}{\boldsymbol{\lambda}} _{s1}^{\rm{T}} = - {\beta _s}(t){{\boldsymbol{B}}^{\rm{T}}}{\boldsymbol{D}}_{}^{\rm{T}}{\left[x_e^p\left(\frac{{\partial {{\boldsymbol{\sigma}} ^M}(t)}}{{\partial {\boldsymbol{\sigma}} (t)}}\right)\right]^{\rm{T}}}{\boldsymbol{\varLambda}} _s^{\rm{T}}{{ {\textit{χ}}}}_s^{\rm{T}} $$ (63) $$\begin{split} &\qquad\quad {-\boldsymbol{C}}{{\dot {\boldsymbol{\lambda}} }_{s2}}(t) + {\boldsymbol{K}}{{\boldsymbol{\lambda}} _{s2}}(t) = \\ &\qquad\quad {\beta _s}(t){{\boldsymbol{B}}^{\rm{T}}_{th}}{\boldsymbol{D}}_{}^{\rm{T}}{\left[x_e^p\left(\frac{{\partial {{\boldsymbol{\sigma}} ^M}(t)}}{{\partial {\boldsymbol{\sigma}} (t)}}\right)\right]^{\rm{T}}}{\boldsymbol{\varLambda}} _s^{\rm{T}}{{ {\textit{χ}}}}_s^{\rm{T}}+ \\ &\qquad\quad {\boldsymbol{S}}_q^{\rm{T}}{{\boldsymbol{\lambda}} _{s1}}(t),\quad {{\boldsymbol{\lambda}} _{s2}}({t_f}) = {\boldsymbol{0}} \end{split} $$ (64) 将计算所得拉格朗日乘子代入方程(62)中, 即可得到应力时域最大逼近函数关于设计变量的敏度值.
最大温度凝聚积分也采用局部最大响应的敏度求解方法. t时刻结构最大温度凝聚积分关于设计变量的显式解析式为
$$ \frac{{\partial {\xi _T}(t)}}{{\partial {\rho _e}}} = \sum\limits_{ni = 1}^{{N_c}} {{\chi _{ni}}} {{{{{\boldsymbol{\varLambda}} }}_{ni}}} \frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}} = {{ {\textit{χ}}}_T}{ {\boldsymbol{\varLambda}}_T}\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}} $$ (65) 由此可得最大温度的凝聚积分函数关于设计变量的显式解析式为
$$ \begin{split} & \frac{{\partial {L_T}(x)}}{{\partial {\rho _e}}} = \int_0^{{t_f}} \Biggr[{({\beta _T}(t) {{{ {\textit{χ}}} _T}} - \dot {\boldsymbol{\lambda}} _T^{\rm{T}}(t){\boldsymbol{C}}} + {\boldsymbol{\lambda}} _T^{\rm{T}}(t){\boldsymbol{K}})\frac{{\partial {\boldsymbol{T}}(t)}}{{\partial {\rho _e}}}+ \\ &\qquad {\boldsymbol{\lambda}} _T^{\rm{T}}(t)\left(\frac{{\partial {\boldsymbol{C}}}}{{\partial {\rho _e}}}\dot {\boldsymbol{T}}(t) + \frac{{\partial {\boldsymbol{K}}}}{{\partial {\rho _e}}}{\boldsymbol{T}}(t)\right)\Biggr]{\rm{d}}t \end{split} $$ (66) 由此可得拉格朗日乘子的计算方程为
$$ - {\boldsymbol{C}}{{\dot {\boldsymbol{\lambda}} }_T}(t) + {\boldsymbol{K}}{{\boldsymbol{\lambda}}_T}(t) = - \beta (t) {{{\boldsymbol{\varLambda}} }} _T^{\rm{T}}{ {\textit{χ}}}_T^{\rm{T}},\quad {{\boldsymbol{\lambda}} _T}({t_f}) = {\boldsymbol{0}} $$ (67) 将计算所得拉格朗日乘子代入方程(66)即可得温度的时域最大逼近函数关于设计变量的敏度解析式.
3. 算例
本节的主要目的是通过具体的数学算例来证明本文方法的必要性和有效性. 所有算例的优化结果均由以下两种方法获得.
方法1: 以往常用的将瞬态传热转化为相同边界条件下稳态传热的方法. 在此情况下, 所关注的目标函数与净传热效率的加权之和被作为目标函数.
方法2: 本文提出的考虑瞬态传热影响的优化方法.
算例1通过二维一体化热防护结构的刚度和轻量化设计问题中两种方法所得优化结果的对比分析说明了本文所提出方法的必要性和有效性. 算例2通过三维一体化热防护结构设计问题证明本文方法可实现瞬态传热下一体化热防护结构的准确设计. 本节所有算例中的模型设计域均由双材料组成, 材料1为作为钛合金, 材料2为隔热石棉, 两种材料的材料属性如表1所示. 设计对象为材料1 的单元密度分布, 体积分数约束的对象也是材料1, 材料2默认分布于设计域内除材料1之外的区域.
表 1 所用材料的属性列表Table 1. Lists the properties of the materials usedDensity/
(kg·m−1)Young’s
modulus/
GPaPoisson’s
ratioThermal
conductivity/
(W·m−1·°C−1)Heat
capacity/
(J·°C−1·kg−1)CTE/
K−1mat-1 4620 96 0.36 21.9 522 ${\text{9} }{\text{.4} } \times {\text{1} }{ {\text{0} }^{ {{-6} } } }$ mat-2 50 0.0001 0.36 0.15 942 0 3.1 二维一体化热防护结构
如图2所示的二维一体化热防护结构模型尺寸为${\text{120}}\;{\rm{mm}} \times 80\;{\rm{mm}}$. 结构的上边和下边各有厚度为2 mm的非设计域, 即顶板和底板, 用黑色表示, 灰色部分为可设计域. 上边承受均布压力载荷, 大小为0.1 N/mm, 底边为固定支撑. 同时, 上边承受均布热流载荷, 大小为0.1 W/mm. 模型的热边界条件为: (1) 当采用本文方法进行瞬态热分析时, 其余边均为绝热边界; (2) 当采用方法1将瞬态传热转换为相同条件的稳态传热时, 底面为60 °C固定温度边界. 求在给定两种材料的体积分数(钛合金30%, 石棉70%)的条件下, 考虑时域最大应力约束(不超过46 MPa)和底面最大温度约束条件, 使结构的应变能或金属材料体积分数最小.
![]()
图 2 二维一体化承载−热防护结构优化模型Figure 2. Optimization model of two-dimensional integrated thermal protection structure3.1.1 刚度设计结果
在刚度设计中, 热载荷大小设定为0.1 W/mm. 在时域最大应力约束和底面最大温度约束(1800 s不超过60 °C, 2400 s不超过85 °C, 3600 s不超过140 °C)条件下通过两种设计方法得到了不同的优化结果. 图3给出了方法1的刚度最大化设计结果及对应的应力分析. 在实际计算过程中方法1平均每次迭代的时间为1.53 s. 由图3可知, 方法1所得优化结果在稳态传热下符合最大应力约束条件. 图4给出了当瞬态热载荷的工作时间为1800 s时, 本文方法所得优化结果及其应力和温度分析结果. 在实际计算过程中, 本文方法在工作时间为1800 s平均每次迭代的时间为6.19 s. 作为对照, 图4同样展示了方法1优化结果在瞬态传热下的应力和温度分析结果. 由图4可知, 在工作时间为1800 s时, 方法1优化结果的应变能为50.23, 本文方法优化结果的应变能为45.17, 比方法1低11.2%, 说明本文方法的优化结果具有更高的刚度. 本文方法所得优化结果的最大应力为46 MPa, 符合最大应力约束要求. 而方法1优化结果的最大应力为50.1 MPa, 超过了最大应力约束上限8.9%. 除此之外, 本文方法优化结果的底面最大温度为59.8 °C, 低于所设定的底面最大温度约束. 而方法1优化结果的底面最大温度为68.6 °C, 超过了底面最大温度约束上限14.3%. 为了进一步对两种方法进行比较分析, 图5给出了工作时间为2400 s的优化结果以及对应的瞬态热力耦合分析结果, 同时以方法1的结果作为对照. 在实际计算过程中, 本文方法在工作时间为2400 s平均每次迭代的时间为6.70 s. 由结果可知, 工作时间为2400 s时, 本文方法优化结果的应变能为46.95, 比方法1优化结果应变能51.86低10.46%, 说明本文方法优化得到了刚度更大的优化结构. 本文方法所得优化结果的最大应力为46 MPa, 符合最大应力约束要求. 而方法1优化结果的最大应力为50.5 MPa, 超过了最大应力约束上限14.2%. 本文方法所得优化结果的底面最大温度为84.4 °C, 符合底面最大温度约束. 而方法1优化结果的底面最大温度为98.2 °C, 超过了底面最大温度约束上限15.5%. 相似的情况也发生在图6所给出的工作时间为3600 s时本文方法的优化结果及其瞬态热力耦合分析结果. 在实际计算过程中, 本文方法在工作时间为3600 s平均每次迭代的时间为7.54 s. 本文方法优化结果的应变能为48.57, 比方法1优化结果的应变能52.32低7.72%, 具有更高的刚度. 本文方法优化结果的最大应力和底面最大温度分别为46 MPa和137.9 °C, 均符合所设定的约束条件. 而方法1优化结果的最大应力和底面最大温度分别为48.6 MPa和159.8 °C, 分别超过了约束上限5.6%和14.1%. 以上不同热载荷工作时间下优化结果的瞬态热力耦合分析结果在表2中给出. 由以上结果可知, 将瞬态传热转换为稳态传热的方法无法保证所得优化结果在瞬态传热下同样满足最大应力和底面温度约束要求; 而本文方法所得的设计结果不但在瞬态传热下具有比方法1优化结果更高的刚度, 且均能在各自的工作时间下满足应力和温度约束条件. 为了进一步对所得结果进行分析, 我们对比了所得优化结果在不同时间下的时域最大应力和底面最大温度. 由图7可知, 方法1的结果在不同时间下的最大应力和底面最大温度均不满足约束条件要求; 本文方法所得的结果都仅在其自身的优化时间下才同时满足最大应力约束和底面最大温度约束, 是当前瞬态热载荷工作时间下的最优解.
![]()
图 3 稳态传热拓扑优化结果Figure 3. Topology optimization results of steady-state heat transfer method表 2 优化结果瞬态热力耦合分析Table 2. Transient thermodynamic coupling analysis of optimization resultsPerformance Method tf = 1800 s tf = 2400 s tf = 3600 s strain energy method 1 50.23 51.86 53.52 proposed method 45.17 46.95 48.57 σMax/MPa method 1 50.1 50.5 48.6 proposed method 46.0 46.0 46.0 TBFSMax/°C method 1 68.6 98.2 159.8 proposed method 59.8 84.4 137.9 ![]()
图 7 不同工作时间所得优化结果的性能分析Figure 7. Analysis results of optimized structure obtained in different working times以上结果说明对于瞬态传热下的一体化热防护结构优化刚度设计问题, 以往常用的将瞬态传热等效为稳态传热的设计方法无法保证优化结果在瞬态传热下同样符合设计要求. 为了保证设计结果能够在瞬态传热下正常工作, 必须在优化设计过程中考虑瞬态热效应对结构响应的影响. 本文方法准确反映了瞬态热效应的影响, 不但得到了比以往将瞬态传热转化为稳态分析温度场方法刚度更高的优化结构, 且优化结构在当前热载荷作用时间下均满足应力和温度约束条件.
3.1.2 轻量化设计
在轻量化设计中, 瞬态热载荷大小设定为0.06 W/mm. 在时域最大应力约束和底面最大温度约束(1800 s不超过60 °C, 2400 s不超过85 °C, 3600 s不超过120 °C)条件下通过两种设计方法得到了不同的轻量化设计结果. 在实际计算过程中方法1平均每次迭代的时间为2.56 s, 本文方法在工作时间为1800 s, 2400 s和3600 s平均每次迭代的时间分别为6.42 s, 6.81 s和6.97 s. 方法1和本文方法的轻量化设计结果及其金属材料体积分数如图8所示. 由图可知, 方法1得到的设计结果的金属材料体积分数为0.232, 高于本文方法在1800 s和2400 s的优化结果, 0.179和0.216, 低于3600 s的优化结果, 0.251. 为了分析各设计结果的性能, 将图8中的设计方案分别进行了不同的工作时间下瞬态热弹性分析, 得到了它们的时域最大应力和底面最大温度如图9 ~ 图11所示. 由图9可知, 当tf = 1800 s时, 本文方法所得设计结果的最大应力和底面最大温度分别为46 MPa和56.64 °C, 均满足约束条件; 方法1设计结果的最大应力和底面最大温度分别为42.78 MPa和45.55 °C. 虽然最大应力和底面最大温度分别比本文方法结果低7.48%和24.35%, 但其金属材料体积分数比本文方法结果高29.61%. 由此可见在满足应力约束和温度约束的条件下, 本文方法得到了重量更轻的设计结果. 图10给出了 tf = 2400 s的分析结果. 结果显示本文方法设计结果的时域最大应力和底面最大温度分别为46 MPa和76.0 °C, 均满足约束条件; 相比之下, 虽然方法1设计结果的底面最大温度为66.2 °C, 比本文方法设计结果低14.8%, 但时域最大应力为48.146 MPa, 超过约束上限4.5%. 在此情况下, 方法1设计结果的金属材料体积分数虽然多于本文方法设计结果, 但最大应力还是超过最大应力约束上限. 在tf = 3600 s时, 方法1得到了比本文方法金属材料体积分数更小的设计结果. 然而图11的仿真结果显示, 相比于本文方法设计结果, 方法1设计结果底面最大温度由102.7 °C变为108.5 °C, 升高了5.65%, 而且其时域最大应力由46 MPa升高至52.6 MPa, 超过最大应力约束上限14.35%. 以上不同热载荷工作时间下优化结果的瞬态热力耦合分析结果在表3中给出. 由以上分析结果可知, 对于瞬态传热下一体化热防护结构的轻量化设计问题, 将瞬态传热转换为稳态传热的方法无法保证优化结果在瞬态传热下满足应力和温度约束条件; 而本文方法在不同热载荷工作时间下所得的设计结果不但均满足应力和温度约束条件, 而且在约束条件都满足的情况下具有更轻的重量.
![]()
图 8 轻量化设计结果及金属材料体积分数Figure 8. Lightweight design results and metal material volume fraction表 3 优化结果瞬态热力耦合分析Table 3. Transient thermodynamic coupling analysis of optimization resultsPerformance Method tf = 1800 s tf = 2400 s tf = 3600 s volume fraction method 1 0.232 0.232 0.232 proposed method 0.179 0.216 0.251 σMax/MPa method 1 42.78 48.1 52.6 proposed method 46.0 46.0 46.0 TBFSMax/°C method 1 45.55 66.2 108.5 proposed method 56.64 76.0 102.7 3.1.3 小结
本节通过二维瞬态传热下一体化热防护结构的刚度设计和轻量化设计问题对以往将瞬态传热转换为稳态传热设计方法和本文提出的考虑瞬态传热设计方法的优化结果进行了对比分析. 分析结果显示, 无论是刚度设计还是轻量化设计, 将瞬态传热转换为稳态传热设计方法无法确保设计结果在瞬态传热下满足应力和底面温度约束要求. 原因是在瞬态传热下, 热载荷的时间效应, 即瞬态热效应, 对结构响应有着重要影响. 由于稳态传热无法准确等效瞬态传热的作用效果, 优化结果在瞬态热环境下无法达到设计要求. 因此, 本文所建立的考虑瞬态传热影响的一体化热防护结构拓扑优化方法是必要的. 相比于将瞬态传热转化为稳态温度场方法, 本文方法由于在优化过程中考虑了瞬态热效应的影响, 在不同热载荷工作时间下所得的结果不但均满足应力和温度约束条件, 而且性能有了明显的提升, 证明了所提出的方法的有效性.
3.2 三维一体化热防护结构
如图12所示的三维一体化热防护结构模型尺寸为${\text{120}}\;{\rm{mm}} \times {\text{6}}0 \;{\rm{mm}}\times {\text{6}}0 \;{\rm{mm}}$. 结构的上边和下边各有厚度为2 mm的非设计域, 即TFS和BFS, 中间部分为设计域. 上边承受均布压力载荷, 大小为0.1 N/mm, 底边的四角为固定支撑. 同时, 上边承受均布热流载荷, 大小为0.1 W/mm, 热载荷作用时间为1000 s. 模型的热边界条件为: (1) 当采用本文方法进行瞬态热分析时, 其余边均为绝热边界; (2) 当采用方法1将瞬态传热转换为相同条件下的稳态传热时, 底面为60 °C固定温度边界. 求在给定两种材料的体积分数(钛合金30%, 石棉70%)的条件下, 考虑时域最大应力约束(不超过46 MPa)和底面最大温度约束(不超过60 °C)条件, 使结构的应变能最小.
![]()
图 12 三维一体化承载−热防护结构优化模型Figure 12. Optimization model of three-dimensional integrated thermal protection structure图13给出以往常用将瞬态传热转换为稳态传热方法的优化结果. 图14给出优化结果在瞬态传热下的热弹性仿真分析结果. 在实际计算过程中方法1平均每次迭代的时间为122.70 s, 本文方法平均每次迭代的时间为788.4 s. 仿真结果基于稳态传热分析的优化结果的应变能为26.45, 时域最大应力为50.50 MPa, 超过了最大应力约束上限9.78%. 优化结果的底面最大温度为55.1 °C, 符合底面最大温度约束条件. 图15给出本文方法的优化结果. 优化结果的应变能为25.31, 低于方法1的优化结果. 这表明本文方法的优化结果具有更高的刚度. 图16给出本文方法优化结果的瞬态热力耦合仿真结果. 结果显示, 优化结果的时域最大应力和底面时域最大温度分别为46 MPa和58.4 °C, 符合应力和温度约束要求. 由以上分析结果可知, 将瞬态传热转换为稳态传热问题方法的优化结果不但刚度低于本文方法的优化结果, 而且时域最大应力超过了应力约束上限. 原因是一体化热防护结构是一个热弹性结构, 准确反应传热过程的作用效果是实现对热弹性结构响应的准确控制的关键. 将瞬态传热转换为稳态传热问题的方法改变了热载荷的作用形式, 造成设计热载荷与实际热载荷存在偏差, 造成优化结构在实际热载荷下性能的下降. 以上不同热载荷工作时间下优化结果的瞬态热力耦合分析结果在表4中给出. 图17显示两种优化结果的最大应力在瞬态传热过程中均是随加热时间的延长而上升的. 将瞬态传热转换为稳态传热问题方法的优化结果的最大应力在地900 s达到了最大应力约束上限, 46 MPa, 并在最终时刻达到了50.5 MPa, 超过了最大应力约束上限. 这是所用设计方法无法准确控制结构响应造成的. 本文方法优化结果在最终时刻的最大应力为46 MPa, 刚好达到应力约束上限. 说明本文方法实现了对结构最大应力的准确控制. 以上结果表明, 将瞬态传热转化为稳态传热问题的方法无法实现对一体化热防护结构最大应力的准确控制. 本文方法准确反应了瞬态传热对结构响应的作用效果. 不但优化结果的刚度高于将瞬态传热转换为稳态传热问题方法的优化结果, 而且优化结果在整个瞬态传热时间区间内的最大应力始终不超过约束上限, 实现了对一体化热防护结构最大应力的准确控制. 以上分析说明本文所建立的考虑瞬态传热影响的一体化热防护结构拓扑优化方法可准确控制瞬态传热下的结构响应, 实现了对瞬态传热下一体化热防护结构的精准设计.
![]()
图 16 本文方法优化结构的分析结果Figure 16. Analysis of optimized structure obtained by proposed method表 4 优化结果瞬态热力耦合分析Table 4. Transient thermodynamic coupling analysis of optimization resultsStrain energy σMax/MPa TBFSMax/°C method 1 26.45 50.50 55.1 proposed method 25.31 46.0 58.4 4. 结论
针对如何在一体化热防护结构的拓扑优化设计过程中准确反应瞬态传热影响的问题, 本文建立了一种考虑瞬态传热的一体化热防护结构拓扑优化方法. 该方法以SIMP法为基础, 考虑结构响应在整个空间和时间域内的变化, 通过最大光滑逼近函数在空间和时间上的凝聚积分描述瞬态结构响应在时域的最大值, 实现对瞬态传热下一体化热防护结构的精准优化. 具体算例验证了方法的有效性和必要性. 通过优化结果的对比分析, 得到以下结论.
(1) 对于瞬态传热下一体化热防护结构的优化设计, 热载荷作用的时间效应, 即瞬态热效应, 对结构设计结果具有重要影响, 在设计过程中不能忽视. 无论是刚度设计还是轻量化设计, 将瞬态传热转换为稳态传热的设计方法无法确保设计结果在瞬态传热下满足应力和温度约束要求. 原因是传统基于稳态热分析的设计方法无法准确反映瞬态传热的作用效果, 导致优化结果无法达到设计要求.
(2) 本文方法可准确控制瞬态传热下的一体化热防护结构的结构响应, 能够保证所获得的设计结果满足应力和温度约束, 实现瞬态传热下一体化热防护结构的精准设计.
-
![]()
图 1 受机−热载荷作用的一体化热防护结构
Figure 1. Integrated thermal protection structure under mechanical-thermal load
![]()
图 2 二维一体化承载−热防护结构优化模型
Figure 2. Optimization model of two-dimensional integrated thermal protection structure
![]()
图 3 稳态传热拓扑优化结果
Figure 3. Topology optimization results of steady-state heat transfer method
![]()
图 7 不同工作时间所得优化结果的性能分析
Figure 7. Analysis results of optimized structure obtained in different working times
![]()
图 8 轻量化设计结果及金属材料体积分数
Figure 8. Lightweight design results and metal material volume fraction
![]()
图 12 三维一体化承载−热防护结构优化模型
Figure 12. Optimization model of three-dimensional integrated thermal protection structure
![]()
图 16 本文方法优化结构的分析结果
Figure 16. Analysis of optimized structure obtained by proposed method
表 1 所用材料的属性列表
Table 1 Lists the properties of the materials used
Density/
(kg·m−1)Young’s
modulus/
GPaPoisson’s
ratioThermal
conductivity/
(W·m−1·°C−1)Heat
capacity/
(J·°C−1·kg−1)CTE/
K−1mat-1 4620 96 0.36 21.9 522 ${\text{9} }{\text{.4} } \times {\text{1} }{ {\text{0} }^{ {{-6} } } }$ mat-2 50 0.0001 0.36 0.15 942 0 
下载: 导出CSV
表 2 优化结果瞬态热力耦合分析
Table 2 Transient thermodynamic coupling analysis of optimization results
Performance Method tf = 1800 s tf = 2400 s tf = 3600 s strain energy method 1 50.23 51.86 53.52 proposed method 45.17 46.95 48.57 σMax/MPa method 1 50.1 50.5 48.6 proposed method 46.0 46.0 46.0 TBFSMax/°C method 1 68.6 98.2 159.8 proposed method 59.8 84.4 137.9
下载: 导出CSV
表 3 优化结果瞬态热力耦合分析
Table 3 Transient thermodynamic coupling analysis of optimization results
Performance Method tf = 1800 s tf = 2400 s tf = 3600 s volume fraction method 1 0.232 0.232 0.232 proposed method 0.179 0.216 0.251 σMax/MPa method 1 42.78 48.1 52.6 proposed method 46.0 46.0 46.0 TBFSMax/°C method 1 45.55 66.2 108.5 proposed method 56.64 76.0 102.7
下载: 导出CSV
表 4 优化结果瞬态热力耦合分析
Table 4 Transient thermodynamic coupling analysis of optimization results
Strain energy σMax/MPa TBFSMax/°C method 1 26.45 50.50 55.1 proposed method 25.31 46.0 58.4
下载: 导出CSV
-
[1] 杨强, 解维华, 彭祖军等. 热防护设计分析技术发展中的新概念与新趋势. 航空学报, 2015, 36(9): 2981-2991 (Yang Qiang, Xie Weihua, Peng Zhujun, et al. New concepts and trends in development of thermal protection design and analysis technology. Acta Aeronautica et Astronautica Sinica, 2015, 36(9): 2981-2991 (in Chinese) doi: 10.7527/S1000-6893.2015.0137 Yang Q, Xie W H, Peng Z J, et al. New concepts and trends in development of thermal protection design and analysis technology [J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(9): 2981-2991. (in Chinese) doi: 10.7527/S1000-6893.2015.0137
[2] 解维华, 霍施宇, 杨强等. 新型一体化热防护系统热力分析与试验研究. 航空学报, 2013, 34(9): 2169-2176 (Xie Weihua, Huo Shiyu, Yang Qiang, et al. Thermal-mechanical analysis and test study of a new integrated thermal protection system. Acta Aeronautica et Astronautica Sinica, 2013, 34(9): 2169-2176 (in Chinese) Xie W H, Huo S Y, Yang Q, et al. Thermal-mechanical analysis and test study of a new integrated thermal protection system [J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(9): 2169-2176. (in Chinese)
[3] Zhao SY, Li JJ, Zhang CX, et al. Thermo-structural optimization of integrated thermal protection panels with one-layer and two-layer corrugated cores based on simulated annealing algorithm. Structural and Multidisciplinary Optimization, 2015, 51(2): 479-494 doi: 10.1007/s00158-014-1137-4
[4] Xie G, Wang Q, Sunden B, et al. Thermomechanical optimization of lightweight thermal protection system under aerodynamic heating. Applied Thermal Engineering, 2013, 59(1-2): 425-434 doi: 10.1016/j.applthermaleng.2013.06.002
[5] Wei K, Cheng X, Mo F, et al. Design and analysis of integrated thermal protection system based on lightweight C/SiC pyramidal lattice core sandwich panel. Materials & Design, 2016, 111: 435-444
[6] Chen Y, Tao Y, Xu B, et al. Assessment of thermal-mechanical performance with structural efficiency concept on design of lattice-core thermal protection system. Applied thermal engineering, 2018, 143: 200-208 doi: 10.1016/j.applthermaleng.2018.07.097
[7] Gogu C, Bapanapalli SK, Haftka RT, et al. Comparison of materials for integrated thermal protection systems for spacecraft reentry. Journal of Spacecraft and Rockets, 2009, 46(3): 501-513 doi: 10.2514/1.35669
[8] 吴书豪, 张永存, 刘书田. 一种考虑瞬态效应的散热结构导热路径设计的拓扑优化模型. 计算力学学报, 2018, 35(5): 547-551 (Wu Shuhao, Zhang Yongcun, Liu Shutian. A topology optimization model for conducting paths design of cooling structures considering transient effect. Chinese Journal of Computional Mechanics, 2018, 35(5): 547-551 (in Chinese) Wu S H, Zhang Y C, Liu S T. A topology optimization model for conducting paths design of cooling structures considering transient effect [J]. Chinese Journal of Computional Mechanics, 2018, 35(05): 547-551. (in Chinese))
[9] Bendsoe MP. Optimal shape design as a material distribution problem. Structural Optimization, 1989, 1(4): 193-202 doi: 10.1007/BF01650949
[10] Yang Q, Meng S, Xie W, et al. Effective mitigation of the thermal short and expansion mismatch effects of an integrated thermal protection system through topology optimization. Composites, Part B: Engineering, 2017, 118: 149-157
[11] Yang Q, Gao B, Xu Z, et al. Topology optimization for integrated thermal protection systems considering thermo-mechanical constraints. Applied Thermal Engineering, 2019, 150: 995-1001 doi: 10.1016/j.applthermaleng.2019.01.067
[12] Xu Q, Li S, Meng Y. Optimization and re-design of integrated thermal protection systems considering thermo-mechanical performance. Applied Sciences, 2021, 11(15): 6916 doi: 10.3390/app11156916
[13] Zhao SY, Li JJ, He XD. Comparison of thermo-structural responses for integrated thermal protection panels with different corrugated core configurations. Journal of Harbin Institute of Technology, 2013, 20(6): 21-28
[14] Xu JF. Thermal impact resistance of integrated thermal protection system of space vehicles. Advanced Materials Research, 2015, 1091: 103-108 doi: 10.4028/www.scientific.net/AMR.1091.103
[15] Meng S, Yang Q, Xie W, et al. Structure redesign of the integrated thermal protection system and fuzzy performance evaluation. AIAA Journal, 2016, 54(11): 1-10
[16] Wei K, Wang KY, Cheng XM, et al. Structural and thermal analysis of integrated thermal protection systems with C/SiC composite cellular core sandwich panels. Applied Thermal Engineering: Design, Processes, Equipment, Economics, 2018, 131: 209-220
[17] Xu Y, Xu N, Zhang W, et al. A multi-layer integrated thermal protection system with C/SiC composite and Ti alloy lattice sandwich. Composite Structures, 2019, 230: 111507 doi: 10.1016/j.compstruct.2019.111507
[18] Li Y, Zhang L, He R, et al. Integrated thermal protection system based on c/sic composite corrugated core sandwich plane structure. Aerospace Science and Technology, 2019, 91: 607-616 doi: 10.1016/j.ast.2019.05.048
[19] K Lin, K Hu, Gu D. Metallic integrated thermal protection structures inspired by the norway spruce stem: design, numerical simulation and selective laser melting fabrication. Optics and Laser Technology, 2019, 115: 9-19 doi: 10.1016/j.optlastec.2019.02.003
[20] Xu N, Xu Y, Zhang W, et al. Design and analysis of multi-layer integrated thermal protection system based on ceramic matrix composite and titanium alloy lattice sandwich. IOP Conference Series. Materials Science and Engineering, 2019, 531(1): 12059 doi: 10.1088/1757-899X/531/1/012059
[21] Shi S, Wang Y, Yan L, et al. Coupled ablation and thermal behavior of an all-composite structurally integrated thermal protection system: fabrication and modeling. Composite Structures, 2020, 251: 112623 doi: 10.1016/j.compstruct.2020.112623
[22] Cao C, Wang R, Xing X, et al. Performance improvement of integrated thermal protection system using shaped-stabilized composite phase change material. Applied Thermal Engineering, 2020, 164: 114529 doi: 10.1016/j.applthermaleng.2019.114529
[23] Xie G, Wang C, Ji T, et al. Investigation on thermal and thermomechanical performances of actively cooled corrugated sandwich structures. Applied Thermal Engineering, 2016, 103: 660-669 doi: 10.1016/j.applthermaleng.2016.04.117
[24] Zhuang C, Xiong Z. A global heat compliance measure based topology optimization for the transient heat conduction problem. Numerical Heat Transfer, Part B: Fundamentals, 2014, 65(5): 445-471 doi: 10.1080/10407790.2013.873309
[25] Zhuang CG, Xiong ZH. Temperature-constrained topology optimization of transient heat conduction problems. Numerical Heat Transfer, Part B:Fundamentals, 2015, 68(4): 366-385 doi: 10.1080/10407790.2015.1033306
[26] Long K, Wang X, Gu X. Multi-material topology optimization for the transient heat conduction problem using a sequential quadratic programming algorithm. Engineering Optimization, 2018, 50(12): 2091-2107 doi: 10.1080/0305215X.2017.1417401
[27] Wu SH, Zhang YC, Liu ST. Topology optimization for minimizing the maximum temperature of transient heat conduction structure. Structural and Multidisciplinary Optimization, 2019, 60(1): 69-82 doi: 10.1007/s00158-019-02196-9
[28] Hyun J, Kim HA. Level-set topology optimization for effective control of transient conductive heat response using eigenvalue. International Journal of Heat and Mass Transfer, 2019, 176: 121374
[29] Wu SH, Zhang YC, Liu ST. Transient thermal dissipation efficiency based method for topology optimization of transient heat conduction structures. International Journal of Heat and Mass Transfer, 2021, 170(3): 121004
[30] Zhang S, Yin J, Liu Y, et al. Multiobjective structure topology optimization of wind turbine brake pads considering thermal-structural coupling and brake vibration. Mathematical Problems in Engineering, 2018, 2018: 1-10
[31] Leader MK, Kennedy G. Thermoelastic topology optimization using steady-state and transient analysis for stress and thermal performance//AIAA Scitech, 2021 Forum
[32] Hu SB, Chen LP, Zhang Y, et al. Design 3D thermo-mechanical structures with multidisciplinary topology optimization. Advanced Materials Research, 2012, 466-467: 1212-1216 doi: 10.4028/www.scientific.net/AMR.466-467.1212
[33] Zhang W, Yang J, Xu Y, et al. Topology optimization of thermoelastic structures: mean compliance minimization or elastic strain energy minimization. Structural & Multidisciplinary Optimization, 2013, 49(3): 417-429
[34] Deaton JD, Grandhi RV. Topology optimization of thermal structures with stress constraints. Schizophrenia Research, 2013, 136(2): 259-263
计量
- 文章访问数: 973
- HTML全文浏览量: 247
- PDF下载量: 246
