国内外学者的大量研究发现，在机床加工精度的影响因素中，其外部环境和内部热源引起的热变形占整个加工误差的40%~70%^[1-2]. 作为制约数控机床精度提高的最主要因素，机床主轴的性能对于机床的切削速度和加工精度至关重要 ^[3]，因此有必要在设计阶段通过热设计理论与技术来提升主轴的热态特性（如温度场、热变形、热应力等），达到减少机床热误差、提高加工精度、降低实验研究和样机制作成本的目的^[4-5].

连续体拓扑优化方法主要有变密度方法、渐进结构优化法、水平集法以及基于数学规划的移动渐进线法等，这些算法需要根据求得的梯度信息来获得最终解. 元胞自动机（CA）法在20世纪90年代开始运用在了结构优化上，例如，Bochenek等^[6]基于CA法的局部规则提出拓扑优化算法，该算法可以快速收敛且更加容易得到优化结果；Faramarzi等^[7]把CA法运用到桁架结构优化即拓扑优化中；Baetens等^[8]研究动力学特性对元胞自动机拓扑结构的影响. CA法无法解决场状态的全局分析问题.

混合元胞自动机法（hybrid cellular automaton method，HCAM）^[9]是由Tovar耦合了CA法和有限元法（FEM）提出的拓扑优化仿生设计方法，这是无梯度优化算法，计算过程中无需进行敏度分析，解决了CA法没有解决场状态的全局分析问题. HCAM采用FEM来评价场状态（应变能），用FEM计算应变能，作为CA的状态信息，该方法简单易懂、收敛性好、计算效率高、减少了数值不稳定性，可以用于多目标问题的求解，为有效进行结构的拓扑优化设计提供新的分析思路和技术手段^[10-11]. 近年来，Tovar等结合能量吸收控制法^[12]、谐波振动载荷法^[13]，将HCAM应用在了车辆防撞器的优化设计、能量收集器外壳和多材料防爆装甲板等的优化设计上^[14]. 郭连水等^[15]针对大变形非线性结构拓扑优化问题，提出基于HCAM的多空间域连续体结构拓扑优化方法，将基于应变、多域的拓扑优化算法引入HCAM中，进行大变形情况下的车辆防撞性研究^[16]. 此外，高云凯等^[17]采用HCAM对车辆的保险杠横梁和车身机构进行优化设计. Da等^[18]开展极端材料属性情况下的拓扑优化研究.

本文以数值传热学相关理论为基础，采用HCAM对温度−结构场耦合结构的连续体开展无梯度热拓扑设计研究，避免在设计过程中目标函数敏度信息的计算，减少数值的不稳定性；使用该方法对典型算例进行对比分析和优化计算，获得更好的结构构形. 使用HCAM对不同工况下的主轴结构算例进行热拓扑设计，获得更优的主轴热态特性，为减少主轴热致变形、提高机床加工精度提供参考.

目前，HCAM都是基于规则的结构元胞，其中最常见的选择是矩形元胞^[19]. 在许多情况下，矩形元胞不可能均匀地覆盖设计区域，在实际工程中分析和设计需要使用不规则的元胞来处理复杂的几何图形，在重点关注区域需要对元胞进行局部细化处理. 使用不规则元胞（如边长自适应变化的三角形元胞）来替代HCAM传统的矩形规则元胞，以更好、更灵活地拟合复杂的几何形状，使得求解拓扑构形时能够获得更优的解决方案. 引入局部网格细化的思想在拟合复杂模型的时候，在应力集中或应变急剧变化的区域进行局部元胞细化，从而不需要对整个结构都使用细元胞，以避免元胞数量过多而导致元胞单元与有限元网格之间映射的难度，加快计算速度，增加结果的真实性.

对于离散位置i和离散时间t，HCAM模型的元胞状态集合可以表示为

(1) ${a_i}(t) = \left\{ {a_i^1(t),a_i^2(t), \cdots, a_i^J(t)} \right\}.$

式中：J为元胞状态的个数， $a_i^J(t)$为对于离散时间t和离散位置i定义的第J种状态.

每个元胞的下一个时刻（t+1）状态由其时刻t状态及其邻域元胞的状态来共同决定. 元胞第J种的状态更新规则可以表示为

(2) $a_i^J(t + 1) = R_i^J\left\{ {a_i(t),a_{i + \Delta 1}(t), \cdots a_{i + \Delta \mathop N\limits^ \wedge }(t)} \right\}.$

式中： $a_{i + \Delta 1}(t),\;\cdots,\; a_{i + \Delta \mathop N\limits^ \wedge }(t)$为元胞i的邻域元胞； $\mathop N\limits^ \wedge $为元胞的邻域元胞个数； $R_i^J$为元胞i的局部更新规则，对于所有位置的元胞更新规则都相同.

在元胞自动机中，局部规则是定义在局部搜索范围内的，某一元胞状态更新所要搜索的空间域是该元胞的邻胞，故在确定规则之前必须定义邻胞大小和类型. 将HCAM传统的规则大小的矩形元胞拓展为不规则大小的非矩形元胞，例如采用边长可变的三角形元胞，如图1所示. 如图1（a）所示为冯诺依曼型邻胞（邻胞数为3），只考虑了3个紧邻的邻胞，这些相邻元胞（灰色三角形）与中心元胞（黑色三角形）存在公共边；如图1（b）所示为摩尔型邻胞（邻胞数为12），与中心元胞存在公共边或公共顶点的元胞都可以认为是中心元胞的邻胞. 不规则元胞排列会造成邻胞的数量不同.

不规则元胞HCAM的邻胞查询规则如下：首先，取出所有元胞及其节点；然后对每个元胞的节点进行比较，若2个元胞有2个相同节点，则定义这2个元胞互为邻胞.

使用HCAM进行拓扑优化迭代的过程中，每个元胞的设计变量 ${x_i}(t)$可以定义为单元密度、几何尺寸、弹性模量、热传导系数等，每个元胞的状态场变量 ${S_i}(t)$可以定义为应变能、热能、热势耗散或者它们的函数，则元胞的状态可以表示为

(3) ${{{a}}_i}(t) = \left[ \begin{array}{c} {x_i}(t) \\ {S_i}(t) \\ \end{array} \right].$

设计变量 ${x_i}(t)$根据材料模型确定，将单元相对密度作为设计变量，场变量 ${S_i}(t)$的定义取决于要解决的优化问题，则在一离散的区域内，结构整体的应变能U为

(4) $U = \sum\limits_{i = 1}^N {{U_i}{v_i}} .$

式中： ${U_i}$为CA中元胞单元i内由于变形而储存的应变能密度， ${v_i}$为元胞单元i的体积.

(5) $\left. \begin{aligned} & {U_i} = \frac{1}{2}{{{\sigma}} ^{\rm{T}}}{{\varepsilon}}, \\ & {{\varepsilon}} = {[{\varepsilon _x},{\varepsilon _y},{\varepsilon _z},{\gamma _{xy}},{\gamma _{yz}},{\gamma _{xz}}]^{\rm{T}}}, \\ & {{\sigma}} = {[{\sigma _{x1}},{\sigma _{y1}},{\sigma _{z1}},{\tau _{xy}},{\tau _{yz}},{\tau _{xz}}]^{\rm{T}}}. \end{aligned} \right\}$

其中， ${\varepsilon _x}$、 ${\varepsilon _y}$、 ${\varepsilon _z}$、 ${\gamma _{xy}}$、 ${\gamma _{yz}}$、 ${\gamma _{xz}}$分别为节点的x向应变、y向应变、z向应变、切应变， ${\sigma _{x1}}$、 ${\sigma _{y1}}$、 ${\sigma _{z1}}$分别为x向、y向、z向正应力， ${\tau _{xy}}$、 ${\tau _{yz}}$、 ${\tau _{xz}}$为xy、yz、xz三个面上的切应力.

应变能密度是结构的局部刚性指标，因此使用应变能密度作为场变量，则最优化标准指应变能密度的均匀分布. 根据该标准，CA状态可以定义为

(6) ${{{a}}_i}(t) = \left[ \begin{array}{c} {x_i}(t) \\ {U_i}(t) \\ \end{array} \right].$

采用HCAM进行结构拓扑优化设计， ${U_i}(t)$的定义可以由具体的优化目标函数确定. 当优化的目的是使元胞的应变能密度均值 $\mathop {{U_i}}\limits^ - $与状态场变量设定值的差值最小化时，优化问题可以用下式表示：

(7) $\left. \begin{array}{l} \min \; \left| {{e_i}} \right| = \left| {{{\mathop U\limits^ - }_i} - U^*} \right| ; \\ {\rm{s}}{\rm{.t}}{\rm{. }}\;\;0 < {x_{{\rm{min}}}} \leqslant {x_i} \leqslant 1,\;\;\;\;i = 1,2,\cdots,N. \end{array} \right\}$

式中： ${e_i}$为误差信号， $U^*$为应变能密度的设定值.

式（7）中有下限 ${x_{{\rm{min}}}}$的要求，可以避免有限元分析过程中产生奇异矩阵.

元胞的状态、结构整体的应变能、 ${U_i}(t)$与规则元胞模型一样，但是元胞单元i的应变能密度的均值 ${\bar U_i}(t)$有变化，可以表示为

(8) ${\bar U_i} = \frac{{({A_i}/{A}) {U_i} + \sum\limits_{j = 1}^{\hat N} {(A_j}/{A}){U_j}} }{{\hat N + 1}}.$

式中：U_j为与元胞单元i相邻的邻胞单元j的应变能密度；A为CA中平均每个元胞单元的面积；A_i为元胞单元i的面积；A_j为与元胞单元i相邻的邻胞单元j的面积.

为了扩大HCAM优化模型的适用范围，加快拓扑优化设计的结果收敛，以适应定载荷和变载荷的结构拓扑优化，计算拓扑优化过程中的应变能密度目标值：

(9) ${U^*} = {\rm{aif}} \times ({U_t}/{V_0}).$

式中：aif为权系数，V₀为初始结构的体积，U_t为迭代到t次时结构的总应变能.

收敛标准由不断更新的设计变量所用的设计规则类型来决定. 将结构质量 $M(t)$的改变 $\Delta M(t)$作为收敛准则，则结构质量为

(10) $ M(t) = \sum\limits_{i=1}^N {{x_i}(t)\rho {v_i}} . $

式中：M（t）为迭代到t次时结构的质量，ρ为实体材料密度.

当 $\Delta M(t)$接近于0时，即质量变化量没有进一步变化，优化的过程就收敛，该状态可以表示为

(11) $\Delta M(t) = M(t) - M(t - 1) \approx 0.$

为了避免收敛过快，收敛准则使用连续2次迭代的平均变化来确定：

(12) $({{\left| {\Delta M(t)} \right| + \left| {\Delta M(t - 1)} \right|}})/({{2{M_0}}}) \leqslant e.$

式中： ${M_0}$为初始结构的总质量.

设计域中的材料分布由局部控制规则决定，它寻求的是状态场变量平均值与状态场变量设定值的差值减小到最小. 采用反向控制策略，设计变量的变化用分段常数函数表示：

(13) ${x_i}(t) = \left\{ \begin{array}{l} + {c_1},\;\;\;\;{e_i}(t) > 0; \\ 0,\;\;\;\;\;\;\;\;{e_i}(t) = 0 ; \\ - {c_2},\;\;\;\;{e_i}(t) < 0 . \end{array} \right.$

式中： ${c_1}$、 ${c_2}$为常数.

为了避免数值计算中的不稳定性，每次迭代过程中设计变量的变化不得超过最大变化量.

为了验证基于HCAM的三角形不规则元胞拓扑优化算法的有效性和普适性，通过不同的图形形状和工况算例来进行拓扑优化设计.

如图2（a）所示的C型结构，基本结构是C型结构的平面板，材料的弹性模量E₀=200 GPa、泊松比μ=0.3、密度ρ=7 850 kg/m³. 约束和载荷如图2（a）所示，在两节点处施加固定约束，并在右侧上、下两角节点处施加方向相反、大小为1 000 N的集中力F. 在该案例研究过程中，采用三角形单元作为混合元胞自动机元胞形状，邻胞选择冯诺依曼型；有限元法计算时选择了PLANE183单元（6节点三角形单元），由于右边上、下两角处施加了载荷以及内圆区域会出现应力集中或应变急剧变化，对这三个区域进行局部细化，采用不规则元胞（网格）划分结构后有8 187个三角形元胞（网格），元胞单元尺寸与有限网格单元尺寸保持一致，以便相互映射，结果如图2（b）所示. 初始变量 ${x_i}(0) = 1$，收敛误差e=10⁻⁶，权系数为0.53，采用反向控制规则，采用的分段函数如下：

(14) ${x_i} = \left\{ \begin{aligned} & \;\;\;1,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{U_i} \geqslant {U^*}; \\ & \;0.75,\;\;\;\;\;\;\;\;\;\;{U^*} > {U_i} > 5{U^*}/{6}; \\ & \;\;0.5,\;\;\;\;\;\;\;\;5{U^*}/6 > {U_i} \geqslant 2{U^*}/3; \\ & \;0.25,\;\;\;\;\;\;\;2{U^*}/3 > {U_i} \geqslant {U^*}/2; \\ & \;\;\;0,\;\;\;\;\;\;\;\;\;\;\;\;\;{U^*}/2 > {U_i}. \end{aligned} \right.$

为了便于观察拓扑优化过程中的构形变化趋势，结合分段函数设置，定义了5种不同弹性模量的单元材料属性以作区分，即 ${E_0}$、0.75 ${E_0}$、0.5 ${E_0}$、0.25 ${E_0}$、0. 最终求解该结果的迭代次数为62次，得到的最终拓扑优化构形如图3（a）所示. 图3（a）中不同的深浅颜色单元对应的是不同材料的单元，其中黑色区域表示弹性模量为0的单元（去除区域），白色区域为弹性模量为 ${E_0}$的单元（保留区域），其他颜色介于两者之间，根据颜色深浅不同可以容易区分出拓扑优化后的构形材料变化趋势以及是否需要保留的程度. 将图3（a）通过HCAM获得的最终拓扑构形与图3（b）中的拓扑构形图^[20]相比，两者的构形结构非常相似，验证了提出算法的有效性和可靠性.

如图4（a）所示的半圆环结构，基本结构是半圆环结构的平面板，大圆半径R=10 m，小圆半径r=3 m，半圆环板的厚度为0.6 m. 材料的弹性模量E₀=1 000 Pa、泊松比μ=0.25、密度ρ=7 850 kg/m³. 约束和载荷如图4（a）所示，左侧圆环直线处加固定约束，右侧圆环直线的两端节点上分别施加方向相反、大小为1 N的集中力F，邻胞选择冯诺依曼型. 有限元法计算时，选择PLANE183单元，采用不规则元胞（网格）划分结构后有7 243个三角形元胞（网格），元胞尺寸与有限单元尺寸保持一致，对加载力的两角和内圆弧区域进行局部细化，结果如图4（b）所示. 初始变量 ${x_i}(0) = 1$，收敛误差e=10⁻⁶，权系数为1.33，采用的局部控制规则是反向控制规则.

最终求解该结果的迭代次数为41次，得到的最终拓扑优化构形如图5（a）所示. 将图5（a）通过HCAM获得的最终拓扑构形与图5（b）中的拓扑构形图^[20]相比，两者拓扑优化后得到的构形结构非常相似.

结构温度场的分布不均会引起结构的热应力，热应力的产生与温度变化和约束有关，在温度变化下，当结构发生自由变形时，不产生热应力；当自由变形受约束时，会产生热应力. 同一物体内部若温度分布不均匀，虽然物体不受外界约束，但由于各处温度不同，每部分受到不同温度相邻部分的约束不能自由伸缩，会产生热应力. 若温升过高，则产生的热应力将超过材料的屈服极限^[21]，因此温度−结构场耦合问题是进行结构分析与优化时需要考虑的重要因素. 考虑在热载荷影响下，通过热拓扑结构优化使得结构的应变能密度分布达到一个较优效果. 当获得了相应的应变和热应力结果时，可以基于HCAM通过寻找设置的应变能密度目标开展温度−结构场耦合情况下的结构拓扑优化.

基于能量守恒定律建立的对于结构模拟件中无内热源的稳定传热，热平衡方程式为

(15) $\frac{{{\partial ^2}\theta }}{{\partial {x^2}}} + \frac{{{\partial ^2}\theta }}{{\partial {y^2}}} + \frac{{{\partial ^2}\theta }}{{\partial {z^2}}} = 0.$

式中：x、y、z为空间坐标，θ为温度.

热问题的基本有限元方程可以由热平衡方程推导求得：

(16) ${{C}}\dot{{ \theta}} + {{{K}}_\theta }{{\theta }} = {{Q}}.$

式中： ${{C}}$为比热矩阵， ${{C}} = \int {_V\rho c{{N}}{{{N}}^{\rm{T}}}{\rm{d}}V} $，其中N为形函数矩阵； $\dot{{ \theta}} $为节点温度变化率向量； ${{{K}}_\theta }$为热传导矩阵， ${{{K}}_\theta } = \int {_V\lambda {{B}}{{{B}}^{\rm{T}}}{\rm{d}}V} $，其中B为几何矩阵；Q为热通量向量.

根据牛顿冷却定律和能量守恒定律，边界条件为

(17) $ {\left. { - \left(\frac{{\partial \theta }}{{\partial n}}\right)} \right|_S} = \frac{\alpha }{\lambda }({\theta _{\rm{s}}} - {\theta _{\rm{f}}}) + \frac{q}{\lambda }({\text{第}}3{\text{类边界条件}}). $

式中：n为结构表面法线方向（向物体外为正），λ为结构模拟件的热导率， $\alpha $为结构模型上的传热系数， ${\theta _{\rm{s}}}$为物体周界温度， ${\theta _{\rm{f}}}$为环境温度，S为结构模拟件边界， $q$为边界表面发热量.

通过线性热应力理论^[21]可知，微元体的总应变由2部分组成：一部分是由应力引起的，另一部分是由温度变化引起的. 用热应力及温差表示应变的平面应力下的广义虎克定律表示为

(18) $\left. \begin{aligned} &{\varepsilon _x} = [{\sigma _x} - \mu ({\sigma _y} + {\sigma _z})]/{E} + a \Delta \theta, \\ & {\varepsilon _y} = [{\sigma _y} - \mu ({\sigma _x} + {\sigma _z})] /{E}+ a \Delta \theta, \\ & {\varepsilon _z} = [{\sigma _z} - \mu ({\sigma _y} + {\sigma _x})]/{E} + a \Delta \theta, \\ & {\gamma _{xy}} = {\tau _{xy}}/G,\;{\gamma _{yz}} = {\tau _{yz}}/G,\;{\gamma _{xz}} = {\tau _{xz}}/G . \end{aligned} \right\}$

式中： ${\sigma _x}$、 ${\sigma _y}$、 ${\sigma _z}$分别为x向、y向、z向热应力，a、G分别为材料的热膨胀系数、切变模量， $\Delta \theta $为温度变化量.

引进符号：

(19) $ \varTheta = {\sigma _x} + {\sigma _y} + {\sigma _z},\;e = {\varepsilon _x} + {\varepsilon _y} + {\varepsilon _z} ,$

则式（18）可以化简为另一种形式的物理方程式：

(20) $\left. \begin{aligned} & {\varepsilon _x} = \frac{1}{{2G}}\left[ {{\sigma _x} - \frac{\mu }{{1 + \mu }}\varTheta } \right] + a\Delta \theta, \\ & {\varepsilon _y} = \frac{1}{{2G}}\left[ {{\sigma _y} - \frac{\mu }{{1 + \mu }}\varTheta } \right] + a\Delta \theta, \\ & {\varepsilon _z} = \frac{1}{{2G}}\left[ {{\sigma _z} - \frac{\mu }{{1 + \mu }}\varTheta } \right] + a\Delta \theta . \end{aligned} \right\}$

将式（20）的三式相加，可以改写成

(21) $e - 3a\Delta \theta = \frac{{1 - 2\mu }}{E}\varTheta. $

式（21）表明，线应变之和与3倍温度应变之差正比于正应力之和. 将式（21）代入式（20），可得

(22) ${\sigma _x} = 2G{\varepsilon _x} + \lambda e - \beta \Delta \theta. $

式中：

$ \begin{split} &\lambda = \dfrac{{2G\mu }}{{1 - 2\mu }} = \dfrac{{\mu E}}{{(1 + \mu )(1 - 2\mu )}}, \\ & \beta = \dfrac{{aE}}{{1 - 2\mu }} =2Ga\dfrac{{1 + \mu }}{{1 - 2\mu }}.\end{split} $

可得由应变分量及温度分量表示的热应力分量公式：

(23) $ {\sigma _y} = 2G{\varepsilon _y} + \lambda e - \beta \Delta \theta, \; {\sigma _z} = 2G{\varepsilon _z} + \lambda e - \beta \Delta \theta. \\ $

对于切应力而言（按切应力互等定律），

(24) $ \left. \begin{aligned} & {\tau _{xy}} = {\tau _{yx}} = G{\gamma _{xy}} = G\left(\frac{{\partial v}}{{\partial y}} + \frac{{\partial u}}{{\partial x}}\right), \\ & {\tau _{yz}} = {\tau _{yz}} = G{\gamma _{yz}} = G\left(\frac{{\partial v}}{{\partial z}} + \frac{{\partial w}}{{\partial y}}\right) , \\ & {\tau _{zx}} = {\tau _{xz}} = G{\gamma _{zx}} = G\left(\frac{{\partial w}}{{\partial x}} + \frac{{\partial u}}{{\partial z}}\right) . \end{aligned} \right\} $

对于连续体结构空间的温度−结构场耦合分析问题而言，计算模型中有15个未知函数、6个应力分量、6个应变分量、3个位移分量，这15个未知函数应满足15个基本方程式，即3个平衡方程式、6个几何方程式、6个物理方程式. 将式（18）~（24）联立，当已知温度场分布时，可以求解热位移、热应力与热应变. 这是本文的温度-结构场耦合分析理论依据.

当采用HCAM进行温度−结构场耦合优化时，将设计空间内的材料分布情况通过CA控制局部变化规则得出，同时利用有限元软件获得因材料密度变化而引起的结果. 在进行有限元计算时，采用顺序耦合分析方法，先进行温度场分析求得结构中的温度分布，再进行结构场分析，将得到的温度场结果作为体载荷加到结构场分析中，求解该情况下的结构热应力分布，通过HCAM来获得热拓扑设计的优化构形.

采用如图6（a）所示的平板结构作为温度−结构场耦合拓扑优化测试用例. 该结构为左右两端和内部有固定约束的平面结构，基本尺寸为 $60\;{\rm{ mm}} \times $120 mm，板厚度为6.0 mm，材料的弹性模量E=100 GPa，泊松比μ=0.25，密度ρ=7 850 kg/m³，a=1.0×10⁻⁵，热导率λ₀=40 W/(m·K). 左、右两端边界上初始温度为 ${\theta _0} = 0$ °C，上、下边界上施加温度载荷 $\theta = 100$ °C，求重量最轻的拓扑构形. 在该案例计算过程中，设置了5种不同热传导率的单元材料，即λ₀、0.75λ₀、0.5λ₀、0.25λ₀、0. 结构邻胞选择冯诺依曼型，利用有限元法计算时选择Plane55单元（3节点三角形单元）来进行耦合场计算. 采用不规则元胞（网格）划分结构后有7 310个三角形元胞（网格），元胞尺寸与有限元网格尺寸保持一致，划分结果如图6（b）所示. 初始设计变量 ${x_i}(0) = 1$，收敛误差 $e = {10^{ - 5}}$，aif=9.9，局部控制规则是反向控制规则.

最终求解该结果的迭代次数为75次，迭代过程中得到的部分拓扑优化构形如图7所示. 图中，黑色区域代表的是须去除区域（热传导率为0），白色区域代表须保留区域（热传导率为λ₀），其他颜色介于两者之间.

如图8所示为优化前、后的结果对比. 从图8（a）、（b）可见，初始结构和最终优化构形的最大应力都出现在两端和中间约束处，黑色区域处的应力是最低的，这个黑色的区域与拓扑优化结果中得到的须去除区域保持一致. 初始结构的最大应力 ${\sigma _{\max }}$约为232 MPa，最终优化构形的 ${\sigma _{\max }}$约为180 MPa，在去除了多余材料的同时应力最大值有约22%的降幅；从图8（c）、（d）中初始结构和最终优化构形的温度场分布云图可以看出，最终优化构形比初始结构有更均匀的温度场分布，结构中的主体部分基本上处于同一温度范围内. 如图9所示为HCAM的性能指标P_i和质量比F随迭代次数t变化的历程曲线图. 可见，拓扑优化已经趋于收敛.

参考文献[5]的主轴系统结构，采用如图10所示的温度载荷作用下的主轴拓扑优化测试用例，材料的E=100 GPa，μ=0.25，ρ=7 850 kg/m³，a=1.0×10⁻⁵，λ₀=40 W/(m·K). 该结构Ⅰ−Ⅷ为安装轴承处，因此在Ⅰ、Ⅳ、Ⅴ和Ⅷ处设置径向固定约束，在Ⅱ、Ⅲ、Ⅵ和Ⅶ设置轴向固定约束，并且在Ⅰ−Ⅷ处施加θ₀=40 °C的温度载荷，主轴其余与环境接触的部位为θ₁=22 °C，求重量最轻的拓扑构形. 将图10所示的温度载荷作用下的主轴拓扑优化测试用例划分为8 792个不规则三角形元胞，元胞尺寸与有限元单元尺寸相同. 在该案例计算过程中，设置了5种不同热传导率的单元材料，即λ₀、0.75λ₀、0.5λ₀、0.25λ₀、0. 结构邻胞选择冯诺依曼型，有限元法计算时选择Plane55单元，并且在Ⅰ和Ⅱ、Ⅲ和Ⅳ、Ⅴ和Ⅵ、Ⅶ和Ⅷ相交处（将产生应力集中处）进行元胞（网格）局部细化处理，将结构划分为8 792个三角形元胞（网格），元胞尺寸与有限单元尺寸保持一致，划分结果如图11所示. ${x_i}(0) = 1$， $e = {10^{ - 4}}$，aif=0.75，局部控制规则是反向控制规则.

最终求解该结果的迭代次数为23次，迭代过程中得到的部分拓扑优化构形如图12所示. 如图13所示为优化前、后的结果对比. 从图13（a）、（b）可见，初始结构和最终优化构形的最大应力都出现在两端和中间约束处，黑色区域处的应力是最低的，这个黑色的区域与拓扑优化结果中得到的须去除区域保持一致. 此外，初始结构的 ${\sigma _{\max }}$约为170 MPa，最终优化构形的 ${\sigma _{\max }}$约为176 MPa，在去除了多余材料的同时，应力最大值只有3.5%的增幅；从图13（c）、（d）中结构优化前、后的温度场分布云图可以看出，最终优化构形比初始结构有更均匀的温度场分布，结构中的主体部分基本上处于同一温度范围内. 如图14所示为P_i和F随t变化的历程曲线图. 可见，拓扑优化已经趋于收敛.

采用如图10所示的主轴结构作为拓扑优化算例模型. 在Ⅰ、Ⅱ、Ⅶ和Ⅷ处施加q₀=90 W/m的热流密度载荷，在Ⅲ、Ⅳ、Ⅴ和Ⅵ处施加q₁=130 W/m的热流密度载荷，其余条件和计算设置与4.1节相同，如图15所示，求重量最轻的拓扑构形.

最终求解该结果的迭代次数为36次，迭代过程中得到的部分拓扑优化构形如图16所示. 如图17所示为优化前、后的结果对比. 从图17（a）、（b）可见，初始结构和最终优化构形的最大应力都出现在两端和中间约束处. 图中，黑色的区域与拓扑优化结果中得到的须去除区域保持一致. 初始结构的 ${\sigma _{\max }}$约为165 MPa，最终优化构形的 ${\sigma _{\max }}$约为186 MPa，在去除了多余材料的同时，应力最大值有11%的增幅；从图17（c）、（d）所示的结构优化前、后的温度场分布云图可以看出，最终优化构形比初始结构有更均匀的温度场分布，结构中的主体部分基本上处于同一温度范围内. 如图18所示为P_i和F随t变化的历程曲线图. 可见，拓扑优化已经趋于收敛.

本文提出基于不规则元胞的HCAM主轴耦合场热拓扑优化设计方法. 该方法以数值传热学相关理论为基础，将边长自适应变化的三角形元胞来替代HCAM传统的矩形规则元胞，更好、更灵活地适应了复杂的几何形状，避免了矩形元胞不能均匀地覆盖设计区域的问题，使得求解拓扑构形时能够获得更优的解决方案. 该方法引入局部网格细化的思想，可以实现在应力集中或应变急剧变化的区域进行元胞细化，不需要对整个结构都使用同一尺寸的规则元胞，从而有效地避免了元胞数量过多、导致元胞单元与有限元网格之间映射的难度. 本文利用不规则元胞的HCAM，对主轴结构进行温度−结构场耦合热拓扑优化设计. 通过温度载荷和热流载荷的算例发现，使用该方法可以有效地去除主轴结构多余材料，改善热态特性.