一种考虑载荷周期分解疲劳寿命计算的结构拓扑优化方法

刘宏亮，付 玲

（1.中联重科股份有限公司建设机械关键技术国家重点实验室，湖南 长沙 410013；2.湖南大学机械与运载工程学院，湖南 长沙 410082）

拓扑优化计算被广泛运用到工程机械动力学、静力学相关的众多结构计算和优化设计问题中，但困扰工程机械的结构疲劳损伤问题，拓扑优化却很难做出对应的求解。据评估，目前工程机械结构的损伤和破坏，80%以上由疲劳损伤引起，如挖掘机动臂的开裂、旋挖钻钻杆、起重机箱型主梁的疲劳损伤等。一方面，疲劳寿命评估是工程机械结构设计中必须考虑的一个环节，结构材料运用的改变如使用高强钢，结构设计形式的变化都需要经过疲劳试验测试、疲劳数值计算等技术手段来评估；另一方面，如将疲劳寿命约束计入结构拓扑优化，能在结构初始阶段得出疲劳性能较优的结构设计形式，可以避免后续的二次设计，增加设计效率。但疲劳计算需要对结构进行载荷时间历程下的瞬态计算，同时疲劳损伤是对应不同周期载荷下的分项疲劳损伤之和，因此，损伤函数是关于优化变量的离散函数（不可微函数），增加了优化计算中寻找全局最优解的难度。目前将疲劳寿命计入拓扑优化模型的研究和技术运用还比较少见，因此，如能通过相关算法将疲劳寿命计入结构拓扑优化计算中，解决离散函数问题，则能提升结构设计效率，提升智能化设计水平。

连续体结构拓扑优化目前已经大范围推广运用。Bendsøe等［1］提出将优化结构离散为多个微元体，通过微元体的删减来实现结构材料的优化布置，Mlejnek 等［2］通过引入中间密度惩罚函数将离散的0～1 问题转化为连续函数求解，后续Bendsøe 等［3］根据中间密度罚函数的思路提出了变密度插值方法。在此基础上，运用相关成熟的数学优化算法如优化准则法［4-6］、序列规划算法［7-8］进一步推广了连续体结构拓扑优化计算在工程结构设计中的运用，如泵车臂架［9］、飞机垂尾结构设计［10］、散热结构优化设计［11］等。变密度法在优化过程中会出现一些数值求解问题，如数值奇异［12］、局部最优问题［13］。为了避免变密度法相关的数值计算问题，Xie等［14］提出了一种新的拓扑优化方法，根据一定的评价标准如单元应力值等进行单元的删除和增加，称为单元生长进化；后续又改进了单元增加和删除策略，提出了双向单元生长计划拓扑优化算法［15］。

现阶段，考虑疲劳寿命的拓扑优化计算主要是将疲劳极限作为应力约束计入到传统的应力水平约束拓扑优化计算中［16-17］。为了进一步考虑疲劳相关特性对拓扑优化计算的影响，相关研究开始探索不同疲劳载荷（高周、低周、比例、非比例等）下疲劳损伤对结构材料布置上的影响。Holmberg 等［18］将计算分成疲劳计算和拓扑优化计算2 个部分，将疲劳许用应力p范数处理，然后将得到的疲劳许用应力作为约束进行拓扑优化计算。Jeong 等［19］考虑了常幅值动态和静态的载荷下疲劳约束拓扑优化计算方法，并给出了对应的结构动态频率响应的拓扑优化计算方法［20］。此外，Oest等［21］运用雨流计数法，并根据Basquin 方程得到疲劳破坏循环次数，给出了比例加载模型的疲劳约束拓扑优化计算。上述方法适应于比例载荷模型，对于更加复杂的非比例载荷问题的拓扑优化，需要大量的计算资源和计算时长，因此，需要通过相关的计算手段来简化疲劳计算。

一种是等效静态载荷方法，Jeong 等［22］给出了将瞬态计算转化为对应的等效静态载荷计算的伪灵敏度分析方法。另一种是线性摄动周期载荷分解，Zhang 等［23］将外部载荷表达成一系列单位载荷的线性叠加，然后对每个单位载荷进行疲劳计算之后叠加得到总的疲劳寿命值，并计入到拓扑优化数学模型中。上述计入疲劳约束的拓扑优化计算均根据变密度法的思路来进行结构材料的优化布置，计算中会产生离散函数不可导［23］等问题。除了运用变密度法，还可运用前述提到的单元生长进化算法来进行疲劳约束下的拓扑优化计算。魏居业［24］尝试了单元双向生长进化拓扑优化计算在疲劳约束方面的运用。根据Wang-Brown 多轴载荷循环计数法［25-26］，采用线性累计损伤计算疲劳寿命，最终成功运用单元生长进化算法得到了基于疲劳计算的拓扑优化设计构型，但拓扑构型存在较多离散结构，构型不清晰。

目前计入疲劳约束的拓扑优化计算研究还在不断完善中，上述的研究给出了一些可行的研究思路，包括比例、非比例载荷加载、单轴、多轴等，但仍然存在计算效率不高、三维结构拓扑优化构型不理想等一些问题，需要后续研究和完善。本文结合相关疲劳计算方法以及单元双向生长进化拓扑优化算法，运用二次程序开发等途径将疲劳寿命作为优化目标因素之一计入结构拓扑优化计算中，提出一种考虑载荷周期分解疲劳计算的拓扑优化计算方法。

1 单元双向生长进化拓扑优化算法

单元双向生长进化［15，27-28］拓扑优化计算在生长进化算法的基础上改进，使得单元在优化迭代过程中根据力学指标能同时增加和删除，增加了优化计算效率。单元生长进化算法的优化变量设定为每个单元体的材料密度值ρi，赋值1或0（1为生长保留单元，0为删除单元）为离散变量。这样能快速地计算结构的体积V(ρ) =∑ρi（下标i为单元编号变量）；同时，对应单元的刚度Ki计算为Eρi，E为材料的固有刚度，通常为弹性模量。优化目标通常是结构柔度：

在此基础上将结构柔度对优化变量进行求导，可得到每个单元体对应的优化计算灵敏度值［29］：

式中：Wi为单元i的结构计算应变能；Ki为单元的刚度矩阵。

在双向单元生长进化算法中，为了避免数值奇异，当单元需要删除时，ρi取一个极小值0.001。此外，通常会考虑距单元i中心点半径为r的子区域来过滤该单元的灵敏度值来确定单元的生长和删除［29］：

式中：T为半径r范围内其他单元的总数；rij为单元j中心点到单元i中心点的距离；权重因子ω(rij)计算为［29］

双向单元生长进化算法为了增加迭代收敛速度，通常会根据前后2 个迭代步的灵敏度值来修正上述过滤灵敏度值［29］：

2 载荷周期分解以及疲劳寿命计算

疲劳计算通常较为复杂，本文为了提升拓扑优化计算效率，给出一种高效疲劳计算方法。疲劳计算首先需要通过相应的计数方法得到载荷谱作用下对应的周期载荷计数。不管是单轴还是多轴，疲劳载荷作用时域曲线都能通过对应计数方法获得最终周期载荷分解，如图1所示，本文将分解得到的周期载荷称为单元基础载荷。常用的计数方法有雨流计数法和Wang-Brown多轴载荷循环计数法［25-27］。

图1 单元基础载荷分解示意图Fig.1 The schematic of basic unit load cycle decomposition

获得载荷谱对应的单元基础载荷后，对每个基础载荷单元运用有限元方法进行求解，获得对应单元载荷周期内的交变应力幅值σa和平均应力σm。根据Morrow 方法修正平均应力σm对交变应力幅值σa的影响［23］：

式中：σf为材料的疲劳强度修正。

最后，根据Basquin 方程计算对应基础载荷单元的疲劳损伤寿命［23］：

式中：b为材料疲劳强度指数，一般情况下为负值［23］。

从上述计算过程不难看出，应力幅值σa越大，平均应力σm越大，则疲劳寿命值越低。

在所有的基础单元载荷周期中，存在一个疲劳寿命最短的基础载荷周期，可称为最危险基础载荷周期。为了提高拓扑优化迭代计算效率，减少计算耗时同时提高迭代收敛速度，本文将从疲劳计数得到的单元基础载荷周期中提取出最危险载荷周期并以此载荷为计算载荷进行拓扑优化计算，即认为在最危险载荷周期下有最优疲劳寿命的拓扑构型同时也是全载荷谱范围内疲劳寿命最优的拓扑构型。根据上述思路，时域内的瞬态疲劳计算可以简化为单元周期载荷下的准静态计算，利于后续计入疲劳寿命的拓扑优化计算。

3 优化模型和优化计算流程

3.1 拓扑优化数学模型

本文拓扑优化计算基于有限元，将初始结构划分单元后赋予每个单元材料密度值ρi并设定其为拓扑优化计算设计变量，结合单元生长进化拓扑优化计算方法，ρi赋值1 或0（1 为单元生长，0 为删除单元）为离散变量。

单元双向生长进化拓扑优化计算通常以结构柔度最小为优化目标，本文通过引入单元基础载荷来定义不同的疲劳周期，并给出了对应的高效疲劳计算方法。为了计入疲劳寿命对拓扑优化计算结果的影响，本文将综合考虑结构柔度和疲劳寿命来设定拓扑优化目标，使得优化构型有较好结构强度的同时能有最优的疲劳寿命。结构一般要求有较大的疲劳寿命，因此为了和结构柔度最小化保持一致，采用单个基础单元载荷周期对应的最小疲劳损伤D=1/N来定义优化目标函数：

式中：C为结构的柔度D在周期载荷作用下的疲劳损伤，通常设定为所有单元中最大的损伤值max(Di)；i为单元编号。

此外，结构优化考虑重量限制，因此可以约束体积上限。设计变量ρi取值为0 或1，能快速地计算结构的体积V(ρ) =∑ρi。因此，以结构强度和疲劳性能为优化目标，计入体积约束以及有限元计算的拓扑优化数学模型建立如下：

3.2 疲劳损伤灵敏度分析

双向单元生长进化拓扑优化算法根据相关评价指标或迭代准则在优化迭代中自动进行单元增加和删除计算，单元生长和删除的评价指标根据灵敏度计算确定。结构柔度的灵敏度分析已在前文有详细叙述，这里给出疲劳损伤D=1/N的灵敏度分析计算过程。同样，对单个有限单元进行求解，单元i疲劳损伤的灵敏度可求解为

式中：为经过平均应力修正后的幅值应力，可等效计算为

式中：为修正后幅值应力对应的应变。

因此，单元i疲劳损伤的灵敏度最终可等效为

在本文后续的数值计算中，同样为避免数值奇异问题，当删除时，ρi取一个极小值0.001。至此，根据结构柔度和疲劳寿命的灵敏度计算，可设定合适的迭代准则来增加和删除单元，实现结构寻优。根据式（8），本文拓扑优化目标函数的灵敏度计算为

根据式（2）和式（12）即可求解。本文采用上述灵敏度值的绝对值进行计算即保证单元生长和删除迭代准则为正；确定灵敏度后根据式（5）进行修正，提升计算收敛效率。

3.3 优化流程

通过上述介绍，本文结合双向单元生长进化拓扑优化算法，寻求结构强度最优柔度同时也追求更长的结构疲劳寿命，将疲劳寿命作为目标元素计入到拓扑优化计算中。优化计算从建立基础构型以及设定载荷边界条件开始，通过有限元离散确定优化设计变量。建立疲劳计算模块，通过单元基础载荷周期分解将时域瞬态计算简化为周期载荷下的准静态计算，经过Morrow 修正幅值应力后经Basquin 方程求得各单元的疲劳寿命，并将单元最小疲劳寿命对应的疲劳损伤计入优化目标。同时考虑结构柔度的优化，将两者的乘积作为拓扑优化目标函数；将结构体积作为约束条件，在优化结构强度和疲劳寿命的前提下减轻结构重量，增加算法的实用性。通过对优化目标结构柔度和疲劳损伤进行灵敏度分析，可确定优化计算单元增加和删除的迭代准则。随后设定优化迭代计算终止条件，进行优化计算并获得最终拓扑优化构型。

优化过程如图2 所示。整个计算流程基于Python 进行程序模块开发而进行相关数值计算和数据提取。程序包含疲劳计算模块、有限元计算模块、单元生长和删除模块、优化迭代计算模块以及后处理模块。

图2 优化流程Fig.2 The schematic of optimization flow

4 算例分析

4.1 悬臂梁模型

首先，采用常见的三维悬臂梁模型来测试算法的有效性。基础构型和载荷边界条件如图3 所示，左端面刚性固定，在右端面下边线中点处施加载荷。材料采用钢材计算，弹性模量E为206 GPa，泊松比为0.3。根据文献［23］，疲劳计算参数b取值为-0.090 9，σf取值为749 MPa。设定疲劳随机载荷谱如图3 所示，根据疲劳计数法进行单元基础载荷周期分解，得到最危险基础载荷周期，对应最小作用载荷为-800 N，最大作用载荷为1 000 N。

图3 悬臂梁模型初始构型和载荷边界条件Fig.3 Basic configuration and load condition for the cantilever model

在此最危险载荷周期下按所建立的拓扑优化数学模型和优化计算流程进行拓扑优化计算，体积分数约束为0.1，即拓扑构型体积不大于初始构型的1/10。由图4 可知，体积分数逐渐减小到约束值0.1，同时结构疲劳寿命逐渐增加，最终优化构型的疲劳寿命值达到3.5×107。根据建立的Python 程序模块进行计算，优化迭代过程体积分数变化，疲劳寿命值的变化和优化构型的变化如图5 所示，最终拓扑优化构型为双层支撑结构。

图4 悬臂梁模型体积分数0.1 条件下的计入疲劳的拓扑优化计算结果Fig.4 Simulation result for the cantilever model with volume fraction

图5 计入疲劳的悬臂梁拓扑优化最终构型Fig.5 Topology configuration for the cantilever model considering the fatigue

此外，为了探索在优化目标中计入疲劳寿命的拓扑优化计算与常规只考虑结构柔度最小的拓扑优化计算的差异，在同样的载荷边界条件以及体积分数约束下，将优化目标函数调整为仅考虑结构柔度最小。最终拓扑构型以及应力云图如图6 所示，2 种计算最终拓扑构型有明显的结构形式差异，相较于考虑疲劳的双层板格支撑结构，不计入疲劳的拓扑构型表现为多层杆支撑结构。从应力上看，计入疲劳目标后结构的应力水平要小于未计入疲劳的拓扑计算构型。说明考虑疲劳后，拓扑构型的结构布置形式发生了明显改变，结构的应力分布更加均匀，同时应力峰值也相应降低。从疲劳寿命上看（见图7），计入疲劳后拓扑构型的最小疲劳寿命由未计入时的3.075×107次提高到3.457×107次，最小疲劳寿命提升12.42%。因此，将疲劳寿命值计入到拓扑优化目标后，最终的拓扑构型在考虑结构刚度足够大的同时能兼顾疲劳寿命的优化，使得结构应力分布更加均匀合理，结构疲劳寿命更优，疲劳性能提升。

图6 计入疲劳与未计入疲劳计算拓扑优化构型的应力云图对比Fig.6 Comparation between fatigue consideration and not

图7 计入疲劳与未计入疲劳计算拓扑优化构型的疲劳寿命对比Fig.7 Fatigue life comparation between fatigue consideration and not

4.2 L型支架模型

为了进一步验证算法的有效性，对经典L 型支架模型进行考虑疲劳计算的拓扑优化计算，并与文献［22］考虑疲劳计算得到的优化构型进行对比。拓扑优化初始构型以及边界条件设置如图8 左下所示，与文献［22］边界条件设置保持一致。同样考虑最优疲劳寿命优化，本文的拓扑构型与文献［22］中的拓扑构型相类似，都是桁架支撑结构，L 型支架底部都采用圆弧边界，同时用多个支撑结构于上边界相连，并合并与L 型结构折点处，总体上布置形式差异不大，验证了本文算法的有效性。此外，本文得出的优化构型，在L 型支架右端出现X 型支撑结构，可降低载荷加载处的应力集中以及疲劳损伤值，提升结构疲劳寿命。

图8 L型支架模型拓扑优化构型对比Fig.8 Comparation for L-shape bracket

同时为了验证本文基于载荷周期分解疲劳计算的拓扑优化求解的计算效率，将本文优化计算的迭代过程和文献［22］的优化计算迭代进行比较，为了直接显示优化寻优收敛速度，以优化迭代计算次数作为指标进行对比分析。本文基于载荷周期分解计算，最终经过123次优化迭代计算完成优化求解，迭代求解过程中，结构体积分数逐次减小，整体过程趋势光滑平缓，拓扑构型各主要框架支撑结构也逐步清晰，迭代40、80、123 次计算的拓扑构型如图9 所示。文献［22］采用载荷静态等效的方法进行疲劳损伤求解，求解效率较完整时域内瞬态计算的疲劳求解有较大提升，但整个计算过程需经过近1 000次优化迭代。优化过程分为2 个阶段：前期基本构型确定阶段和后期细节优化阶段，如图10所示。前期只需50次左右迭代即可确定基本构型，此阶段体积分数快速下降，但拓扑构型较模糊，支撑结构轮廓不清；后期细节优化阶段需经过950次左右迭代计算，此阶段体积分数平缓下降直到最后迭代终止。综合比较分析，本文基于载荷周期分解考虑疲劳计算的拓扑优化计算，迭代计算收敛速度始终处于较合理水平，总迭代次数较少，文献［22］前期体积分数快速下降，但后期需要较多次优化迭代才能使拓扑构型清晰。因此，本文基于载荷周期分解考虑疲劳计算的拓扑优化计算从迭代次数上反馈需要相对较少的迭代计算次数，有较高的计算效率。

图9 L型支架本文优化迭代求解过程Fig.9 The iteration process for the L-shape bracket

图10 L型支架优化迭代求解过程（改编自文献［22］）Fig.10 The iteration process for the L-shape bracket（adapted from［22］）

5 计算参数的影响分析

本文计入疲劳寿命的单元生长进化拓扑优化模型中设定了一些计算参数，为了考察算法的稳定性，需分析参数取值不同时，拓扑优化构型的变动情况。疲劳计算过程中设定的参数b和σf是疲劳计算参数，主要影响疲劳寿命值的计算精度，对拓扑优化构型影响有限；并且本文关于b和σf的取值根据参考文献选取，已经过多方测试和验证对钢制构件的疲劳计算有较高的计算精度［23］，因此不作比较计算。其他参数还有灵敏度计算过滤半径以及初始构型的单元网格大小。

现以悬臂梁模型为例分析上述2个参数对拓扑构型的影响。在前述计算中，初始构型网格尺寸为3，灵敏度过滤半径为9。首先，测试不同过滤半径对最终拓扑优化构型的影响。如图11所示，灵敏度值过滤半径在一定范围内对最终拓扑构型基本没有影响。相对地，单元网格尺寸对拓扑构型的影响则比较明显，如图12 所示。单元网格尺寸越小，优化构型细微结构越多，反之网格尺寸越大则拓扑构型更简洁。因此，本文所建立的将疲劳寿命计入目标函数的拓扑优化计算方法，对灵敏度计算过滤半径的变动不敏感，主要的影响因素是单元网格尺寸。对于过小的单元网格尺寸，优化构型的棋盘格现象比较明显，主要由算法对每个单元网格单独求解并进行离散生长和删除引起。后续的改进措施不在本文的讨论范围内，将在后面的工作中提及。因此，计算时需设定合适的网格尺寸，可尝试几种网格尺寸大小，取计算构型最好的构型进行结构设计。

图11 灵敏度计算过滤半径变化对拓扑构型的影响Fig.11 The influence of the sensitivity calculation radius on the topology configuration

图12 单元网格尺寸对拓扑构型的影响Fig.12 The influence of the mesh size on the topology configuration

6 结语

本文结合双向单元生长进化拓扑优化算法，将疲劳计算考虑到结构拓扑优化中，给出了一种计入疲劳寿命的结构拓扑优化计算方法。通过相关对比分析，本文建立的疲劳寿命计入目标函数的拓扑优化算法能给出精度较高的拓扑优化构型，同时相较于有关文献，优化有较高的计算效率。同时本文还对算法中相关参数的影响进行分析，给出了对应的参数选择方案。

本文给出的考虑疲劳计算的拓扑优化方法能用较少的迭代次数给出合理的拓扑构型，但仍然有如下几点需后续进一步研究进行改进和验证：

（1）虽然计入了过滤算法改善拓扑构型成型结构边界的光滑度，但受到有限元单元形状的限制，拓扑构型没能给出光顺的边界轮廓，后续还需算法改进，完成拓扑构型边界光顺化处理。

（2）在悬臂梁模型算例中，对比未考虑疲劳计算的拓扑优化构型，考虑疲劳后的拓扑构型采用同样的疲劳寿命计算方法所得的疲劳寿命更优，但提升程度没达到数量级提升。因此，考虑疲劳计算的复杂性和计算精度，后续还需通过样件制造和试验测试来进一步验证。

（3）本文只考虑最危险载荷下的疲劳计算，其他载荷的影响尤其是对于多轴疲劳问题暂未考虑，上述简化的影响还需后续的研究进一步分析，算法总体有效性还需后续实机产品测试进行验证。