基于偶应力理论的双向进化结构优化法*

乔赫廷，高 群，熊 展

（沈阳工业大学机械工程学院，沈阳 110870）

0 引言

在微/纳尺度技术领域中，微机电系统（Micro-Electro-Mechanical System，MEMS）的结构设计逐渐成为研究的热点问题。拓扑优化方法由于其能够产生结构的创新式构型，可获得比传统的形状和尺寸优化更好的优化结果而被广泛应用于传统宏观机械结构概念设计中，成为机械领域中实现轻量化设计不可或缺的手段。但是在微观机械结构领域中的发展却较为迟缓。因此，为了在MEMS中实现结构的创新式构型设计，必须要研究在微观尺度下的结构拓扑优化方法。

在宏观机械领域，结构拓扑优化理论与方法逐步完善，拓扑优化方法被视为寻求创新结构设计的有力工具。S Mantovani等[1]针对方程式赛车转向柱的设计问题，提出了一种以质量最小化为目标，同时遵循严格的结构约束的拓扑优化方法。Quhao Li等[2]对于薄壁结构承重构件的加强筋布局和截面轻量化设计问题，提出了一种并行拓扑优化方法，显著提高了结构的性能，同时也降低了计算成本。毕政等[3]面向如何提升车辆底部防护组件的抗爆性能，以降低车身底板变形对车内乘员的威胁问题，提出了一种结合混合自动元胞法的拓扑优化设计方法，得到了防护组件中加强梁的最佳材料分布形式。同时，拓扑优化方法拓宽了其在微观尺度上的应用，以设计满足工作需求的微型结构。微机电系统中的结构设计也需要先进的方法来充分发挥其作用，因此进行微观结构的设计显得尤为重要。Abbas Homayouni-Amlashi等[4]为了设计一种从来自不同方向的面内力中获取能量的压电板能量采集器，采用拓扑优化方法寻找压电板的最佳布局，以最大化电输出并克服电荷抵消。Yi Gao等[5]考虑微观结构的不确定性，提出了一种基于可靠性的拓扑优化框架，避免了不确定性对所得拓扑体积分数和局部特征的影响，实现精确的目标可靠性。此外，随着材料微观结构尺寸变小，非均质材料的有效特性取决于结构尺寸和材料微观结构尺度的大小，这种现象称为尺度效应。众所周知，一些高阶介质理论，例如偶应力理论可以提供合理的方法来解决这个问题，捕获结构的力学性能和材料尺度参数，揭示优化设计中的尺度效应。Sergey P Pavlov等[6]在Timoshenko运动学假设和修正的偶应力理论的基础上，提出了一种基于拓扑优化的纳米梁微观结构优化方法，获得了纳米梁在任意静态和动态载荷以及不同边界条件下的最佳拓扑，以增加其刚度。Wenzheng Su等[7]为了描述和检验周期性蜂窝状固体微观结构拓扑设计的尺寸相关性，以实现结构基频最大化，在偶应力理论下建立优化模型，充分说明了周期性蜂窝状固体微观结构中优化结果的尺寸依赖性并改进了经典理论在尺度效应描述方面的弱点。

基于变密度（Solid isotropic material with penalization，SIMP）法的拓扑优化中，优化过程中会出现数值问题。如何解决数值问题已成为SIMP法优化求解中的重要组成部分。在这方面，许多研究者开展了非常有价值的工作。Quhao Li等[8]为了消除固有频率拓扑优化过程中经常出现的数值问题，结合“有界公式”和“鲁棒公式”建立了动态拓扑优化公式。Haitao Han等[9]提出了一种孔洞填充方法（HFM），用于控制最优结构中孔洞的存在，不仅限制了最优结构设计中的孔的数量，而且抑制了棋盘格和网格依赖现象。与SIMP法相比，双向进化结构优化（Bi-directional evolutionary structural optimization，BESO）法允许添加和移除元素，收敛速度更快，效率更高，得到的优化结构更清晰，现如今已成为一种寻找最优解更稳健且更有效的算法。Mohsen Teimouri等[10]采用BESO法对新的混合实体网格结构进行拓扑优化设计，显著提高了结构的机械性能。Yunkai Gao等[11]提出了一种改进的BESO法解决多约束拓扑优化问题，结果表明不同的约束组合产生不同的材料分布。Xiuyang Qian等[12]利用BESO法，对机翼截面进行拓扑优化设计，以去掉结构中的低效材料，证明了拓扑结构的最佳材料布局在降低应力集中和提高机翼壁板承载能力方面的表现更好。

综上所述，本文的研究目的是提出一种在微观尺度下应用BESO法实施结构拓扑优化设计的方法。基于偶应力理论而不是经典连续介质理论构建了矩形8节点偶应力单元。在此基础上考察了各向同性偶应力介质刚度最大化的优化结果，并与相同工况下的经典连续介质结果进行了对比。数值算例表明，目前基于偶应力的优化模型详细说明了微型结构设计中优化结果的尺寸依赖性，并改进了经典优化模型在尺度效应描述方面的弱点。

1 偶应力理论

本文中，ui为线弹性偶应力模型中的位移矢量，κi为独立旋转矢量，χij=κj,i，为曲率张量。应变张量εij依赖于位移矢量和旋转矢量。应力张量σij和偶应力张量mij通过本构方程与εij和χij存在联系。本节阐述了在偶应力理论下的基本方程，其中包括几何方程、平衡方程和本构方程，以及广义有限元列式。与经典弹性理论相比，平面偶应力理论除了原有的σxx、σyy、σzz、τxy、τyx应力分量外，新增了mxz、myz两个偶应力分量，采用符号表示如图1所示。

图1 二维偶应力模型

1.1 基本方程

几何方程：

式中：εij为应变张量；uj,i为位移向量；ej i为置换符；κk为独立的旋转向量；χij为曲率张量；κj,i为旋转向量。

平衡方程：

式中：σij,i为应力张量；ψj为体积力；mij,i为偶应力张量；ei jk为置换张量；σik为独立的应力张量；φj为体力偶。

对于平面偶应力问题，剪应力是对称的，分别用τs和τa表示，即（τxy≠τyx）：

对称部分τs与剪切应变γx y相关：

式中：μ为剪切模量，μ=E/[ 2( 1+δ)]，E为杨氏模量，δ为泊松比。

反对称部分τa与刚体局部转角ω相关，并且它不是一个独立的量，具体可以表示成如下形式：

式中：u和v分别为面内的位移和转角。

与之对应的偶曲率为：

考虑各向同性材料，偶应力下的本构关系表示为：

式中：σi j、εk l分别为应力和应变；Dijkl为弹性本构矩阵。

边界条件：

对于平面应变条件下的二维连续体，本构方程可以用矩阵形式改写为：

式中：λ=D·δ/( 1+δ)( 1-2δ)，为拉梅常数；μg为类似于剪切模量的常数；lg为材料的特征长度。

1.2 有限元列式

本文研究了矩形8节点偶应力单元，其每个节点上分别有3个自由度，如图2所示，对于偶应力模型，根据势能原理，则单元总势能泛函为：

图2 平面矩形单元

位移场函数表述如下：

式中：ui为单元节点线位移；vi为节点角位移；ξ和η为局部坐标系中的自然坐标；Ni为形函数。

Ni在8节点单元中分别表示成：

将式（14）代入（13）取变分，并使其为0，得到的有限元列式可以简写成如下：

式中：Ke为单元刚度阵；Ue为节点位移列向量；Fe为节点载荷向量。

依据经典单元矩阵的表现形式，得出单元刚度矩阵，即：

2 双向进化结构优化算法

2.1 优化模型

刚度是影响结构能否稳定工作的关键参数，通常情况下，结构优化的最终目的是实现刚度最大化，也就是说，在给定的材料约束条件下，需要考虑的问题是如何使结构柔度最小，以提高抵抗弹性变形的能力。其优化问题可以定义为如下形式：

式中：c(x)为平均柔顺性；FT为力矢量的转置；U为位移矢量；K为总刚度矩阵；V p为单个元素的体积；V*为规定的总体积；N为网格中的元素总数；xp为设计变量，为了防止刚度矩阵产生奇异问题，使用较小的xmin值来表示空洞元素。

2.2 求解策略

在优化过程中，BESO法的工作原理与SIMP法有所不同，它的材料去除方案允许添加或移除材料，这主要取决于元素灵敏度数值的大小，当灵敏度为正时，添加材料；当灵敏度为负时，移除材料，同时要通过灵敏度过滤方案来避免数值不稳定现象的出现。因此，使用下式进行滤波：

式中：Q为设计域中单元的节点总数；b ho为元素中节点h和o中心的距离；ιo为灵敏度数；ς(b ho)为权重因子。

式中：R为滤波器半径。

式（20）表示要将在该滤波半径范围内的元素灵敏度值作为新的灵敏度数值，以获得更好的优化结果。

上述灵敏度滤波过程称为空间平滑，但是存在收敛问题，为了保证优化过程的绝对收敛，引入时间平滑的概念，其主要思想是通过取灵敏度的平均值来防止优化过程中结构性能发生振荡。它可以表示为：

其中t和t-1分别表示当前迭代次数和上一次迭代次数，在每次优化过程中，=φp将会被用到下一次迭代。

在BESO法中，为了确定下一次迭代的目标体积，采用如下方程计算：

其中M为进化体积比，当V＞V*时，则根据上式计算并更新得到下一次迭代的目标体积，当V＜V*时，则结构的体积将在之后迭代中保持不变。

图3所示为BESO法的循环流程，包括确定目标函数和目标体积约束的详细步骤。

3 数值算例

在这一节中，分别给出了几个例子来验证方法的可行性。所有算例均采用的是基于偶应力下的8节点矩形单元求解，并且假设所有固体材料的杨氏模量为1，泊松比为0.3，其进化体积比设为0.04。

优化过程的收敛准则表示如下：

其中ρc(x)设为0.5%，ϑV设为0.4%，Z为整数，取值为5。此外，如果迭代次数达到200次，则优化将被终止。

3.1 算例一

本例模型为长L=90 mm，宽H=30 mm的Michell桁架，底部两端固定约束，并在中间施加竖直向下的载荷，如图4所示。将该模型离散化为90×30个8节点偶应力单元，体积约束为0.5。首先从图5所示的变形示意图上可知，结构自身各部分的变形程度分配较为平均。符合Michell桁架的工作机理以及微型结构的工作机制。其次，研究了不同尺度下模型的优化结果，将偶应力材料的特征长度lg与梁高度H的比值称为尺度参数，记为ζ。如表1所示，依次考察了ζ为0.1，0.2，0.4，0.6，0.9以及1时的优化结果。随着ζ的增加，优化结果会有明显的差异。当ζ在0.1～0.4之间，拓扑结构几乎相同，具体地说，结构内部的孔洞变得越来越圆滑。当ζ=0.6时拓扑结构发生了显著的改变，孔洞由2个变为1个。当ζ=0.9时，孔洞彻底消失，最终得到稳定的结构。应该指出，当ζ增加到与域的维数相当时，拓扑结构会发生很大的变化。绘制桁架在ζ=0.6时的平均柔度、体积分数和拓扑优化的演变历史，以研究BESO法的质量，如图6所示。通过比较，随着体积分数的减小，平均柔度先增大，然后在达到目标体积后收敛到一个几乎恒定的值。经过48次迭代，最终收敛到稳定的拓扑结构。

图4 底部两端固支Michell桁架结构

图5 Michell桁架结构优化结果变形示意图

表1 Michell桁架的最佳拓扑结构

图6 Michell桁架的平均柔度和体积分数的演变历史

3.2 算例二

在这个例子中，采用长为L=80 mm、高为H=30 mm的悬臂梁作为研究对象，左端固定约束，上方承受均布载荷，如图7所示。宏观设计域离散化为160×60个8节点偶应力单元，材料的体积约束为0.4。最优构型的变形如图8所示，最大位移发生在结构上端，由此可以得出在微观尺度下悬臂梁的结构设计符合最大刚度要求。同时为了比较本文采用的BESO法和SIMP法的优化设计，在不同尺度下，应用BESO法求解MBB梁问题。从表2中可以看出，随着尺度参数的增大，MBB梁的柔度逐渐变小，迭代时间和次数减少，说明尺度的变化会对优化过程造成影响，并且宏观尺度和微观尺度下结构的性能有明显的差异。

图7 承受均布线载荷的悬臂梁

图8 均布线载荷悬臂梁优化结果变形示意图

表2 不同尺度下承受均布线载荷的悬臂梁的结果比较

表3所示为ζ取0.1、0.6和1尺度下的迭代历史图，其中2列和4列为算例模型同一尺度下不同算法的迭代过程，可以看出SIMP法在迭代过程中会出现灰度单元，这一点充分说明了本文采用的BESO法比SIMP法在数值求解中稳定性更好。并且该表的3～5列反映了材料微观尺度下结构的特性，因此偶应力介质拓扑优化具有明显的尺度效应。随着精密加工工艺的出现，往往需要设计微型结构的尺寸或与之配合的构件，而基于偶应力理论的拓扑优化方法在这一应用领域表现得更出色。

表3 迭代历史图

4 结束语

本文通过将偶应力理论与BESO法结合，针对微型结构实行最小柔度拓扑优化设计，事实证明了偶应力介质下结构存在显著的尺度效应，解决了经典介质理论无法得到微观结构的信息，提供了一种微观结构设计的新方法。最终的拓扑结果取决于尺度参数的大小，且不同尺度下结构的表现形式是不同的，当尺度参数相差较大时，宏观尺度下产生的是桁架结构，而微观尺度下会出现不规则的椭圆形孔。并且偶应力介质具有捕获微观结构信息的能力，因此尺度效应明显，这一点是经典连续介质无法做到的。此外，与SIMP法相比，本文所用的BESO法不会产生灰度单元，数值稳定性更好。