连续体结构的静动态联合拓扑优化

秦东晨，赵鹤鸣，王然然

1 引言

拓扑优化以其能够进行结构创新型设计和结构大幅度减重设计而成为工程设计领域不可或缺的一种方法，对各种结构进行拓扑优化研究是现在比较热门的一项课题。结构动力学拓扑优化是结构优化中最富有挑战性的研究领域，结构拓扑的改进可大大改善结构的静、动态性能[1]。目前，结构的动力学拓扑优化比起静力学拓扑优化研究还很有限，几乎刚刚起步，由于连续体结构拓扑优化模型描述的困难和数值算法计算量的巨大，因而发展较慢。

在拓扑优化中，考虑动态性能的要求，指在考虑静态刚度的同时，兼顾结构的动态振动频率，使得拓扑结构更加合理。而针对频率的动力特性拓扑优化目标函数是在满足结构刚度约束的情况下改善结构的模态特性，使结构整体刚度提高、材料得到优化配置。关于动态拓扑优化设计的研究，国内外的文献还很有限，因此很有研究的必要性。文献[2]以结构最小化振动烈度值为目标，以模态频率为约束，建立结构振动烈度优化模型，采用变密度法固体各向同性微结构材料惩罚模型，实现发动机机体的低振动优化。文献[3]以静态柔顺度或节点位移最小化和动态特征值最大化加权函数为目标函数，提出并建立了静动态多目标连续体结构拓扑优化模型，并通过数值算例验证了模型及方法在静动态多目标下连续体结构拓扑优化设计中的可行性和有效性。

2 拓扑优化模型

静动态联合拓扑优化建立优化模型时，静态约束参数采用结构中关键节点的位移，节点位移基于有限元理论通过弹性力学方程求解，动态约束参数选择结构的基频，基频通过求解动力学特征方程获得，同时以结构的总重量最小作为目标函数，基于变密度法[4]建立优化模型如下：

式中：m—结构总重量；ρ—材料密度；xi—单元相对密度；K—结构刚度矩阵；λ—特征值；M—结构质量矩阵；F—施加的载荷；U—结构位移；Up—p点的位移；para1、para2用来表征对结构基频和p点位移的约束大小。

3 有限元理论

3.1 单元类型选择

对于一个力学或者物理问题，在建立其数学模型以后，用有限元方法对它进行分析的首要步骤是选择单元形式。对于平面问题，常见的有3节点三角形单元和4节点矩形单元，采用相同数目的节点时，矩形单元的精度要比三角形单元的精度高，而三角形单元简单灵活，对斜交的边界和曲线边界有很好的适应性[5]。由于研究的对象是梁结构，故而选用矩形单元，如图1所示。

图1 平面4节点矩形单元Fig.1 Plane Rectangular Element with Four Nodes

典型的4节点矩形单元位移函数的选择一般选取多项式，求解方便，也可以是三角函数或其它有理公式，多项式中的项数由单元节点数决定，多项式方程的个数取决于单元的自由度数，在多项式中考虑到材料的各向同性，需保证多项式中各项必须使x和y对位移有相同的影响，一般以杨辉三角形选择多项式，建立平面4节点矩形单元位移函数如下：

式中：δ—单元中任意一点的位移矢量；u—水平方向位移；v—竖

直方向位移；α1，α2，…，α8—待定系数。

将节点在局部坐标中的坐标值代入位移插值函数，可列出四个节点处的位移分量，即两组四元联立方程，由此可以求得位移模式中的8个未知参数α1，α2，…，α8，再把这些参数代入位移插值函数中，便可得到用节点位移表示的位移模式：

式中：a、b—矩形单元四个节点的坐标，由结构总尺寸和网格划分数量可求得。

3.2 单元刚度矩阵

确定了单元位移后，可以很方便的利用几何方程和物理方程求得单元的应变矩阵如下：

单元刚度矩阵具有对称性、奇异性、主元恒正等特性，可由式（5）求出，式中积分部分可以利用MATLAB中积分函数int或者采用两点Gauss型求积公式求解。

式中：D—弹性矩阵可由广义胡克定律推导求出如下：

3.3 总刚度矩阵组装

在完成单元刚度矩阵的求解后，需要计算结构刚度矩阵从而可以根据弹性力学求解系统中所有节点的位移。结构刚度矩阵是单元刚度矩阵通过式（6）集成得到。

式中：K—结构总刚度矩阵；Ke—单元刚度矩阵；G—单元节点自由度与结构节点自由度的转化矩阵。

在实际编程计算的过程中这个集成过程是在计算得到单元刚度矩阵Ke后，按照单元的节点自由度编码，对号入座地叠加到结构刚度矩阵的响应位置上即可实现。具体是通过节点在网格中的位置，计算得到单元中所有节点的编号，然后依据节点自由度数计算单元节点对应系统自由度中的位置编号，并将所有单元的单元刚度矩阵依次加在对应位置上，即可获得总刚度矩阵。

4 特征值求解

固有频率是系统的固有属性，与外载荷无关，系统的自由振动方程如下：

理论分析与实践均表明，阻尼对结构的固有频率和振型影响不大，所以在求解结构的固有频率时，可以不计阻尼的影响，故而建立特征方程如下：

应用特征向量理论求解特征方程有着广泛的应用，特别是工程技术领域中的振动问题和稳定性问题，特征方程求解的方法有很多，如矢量迭代法、子空间法、阻尼法等[6]，获得特征值解后即可求得结构相对应的固有频率。主要采用MATLAB中eigs函数[7]和矢量迭代法求解特征值问题。

MATLAB提供了专门的函数来求解特征值问题，即函数eigs，一般用来获取大型特征值问题。在求得结构的整体刚度矩阵K和质量矩阵M后，可以直接调用函数eigs获得结构的某几个低阶特征值，其调用格式如下：

式中：K—结构总刚度矩阵，基于前面有限元推导过程求得；M—总质量矩阵，推导过程与K类似；num—返回特征值个数；SM—返回绝对值最小的特征值；LM—返回绝对值最大的特征值。

矢量迭代法是对结构进行动态特性分析时常用的一种方法，矢量迭代法可以分为正迭代法和逆迭代法，前者向结构的最高频率和振型收敛，而后者向结构的最低频率和振型收敛。因为在实际工程中对结构影响较大的是低阶的频率和振型，而且以结构的有限单元离散化模型为基础计算结构的高阶特征时精度很差，通常选用逆迭代法而不使用正迭代法[8]。在具体实施迭代求解计算的过程时，可以有各种不同的步骤，但是其基本思想及使用的迭代关系式是相同的，如式（9）所示。

式中：Φ—固有振型矩阵；Ω—固有频率矩阵。

逆迭代法是利用刚度矩阵窄带性质的首选方法，因为其是使用刚度矩阵的逆矩阵，所以此法向最低振型收敛。为保证刚度矩阵的窄带性质，不形成动力矩阵，而是把假定的位移向量和质量矩阵结合，然后求解基于刚度的平衡联立方程组后得到改进的位移向量，解方程组通常采用LU分解法。

5 灵敏度过滤

拓扑优化过程会产生数值不稳定现象，常见有棋盘格、网格依赖现象。为了抑制这些问题，常在求解完成后，在后处理中，采用滤波法对优化结果进行处理，常用的滤波方法有密度过滤、灵敏度过滤以及灰度过滤[9]这三种方法。灵敏度过滤这种方法较为常用，其基本思想是使用新的灵敏度代替旧的灵敏度实现数值不稳定的过滤。灵敏度过滤使用如下的公式进行：

式中：c—目标函数；γ—一个较小的数，一般取10-3，用来避免分母为零的情况。

6 优化准则算法

有限元法的求解精度受单元的数量与形状影响，所以为达到所要求的精度，需要单元的数量较多，求解数量较多的线性联立方程计算工作量大，采用优化准则方法是一种很好的选择。

在非线性问题的求解算法领域，OC准则法[10]较为古老，但是由于其简单性和高效性，目前还被广泛的应用于拓扑优化领域。OC准则是基于库伦塔克（K-T）条件的，当满足约束0≤X≤1，即该约束为不起作用约束时，当公式拉格朗日乘子λ时，则解为极值点。

将上述公式进行变形得到Bi＝1，其中：

将上述灵敏度分析结果代入上式中得到：

采用如下公式进行设计变量的迭代更新：

式中：m—移动步长，结构刚度最大化问题中通常取为0.2，η取为0.5，λ通过二分法求得。

7 算例分析

本节对平面简支梁模型，基于有限元理论进行固有频率求解和拓扑优化分析，两个算例均是在MATLAB软件中编程实现。

算例一：

平面结构，如图2所示。厚度为5mm，其材料特性为：弹性模量 E=69GPa，泊松比为 0.3，密度 ρ=1kg/cm3。

图2 平面简支梁结构Fig.2 Structure of a Simple Supported Beam

现分别采用MATLAB提供的函数eigs和矢量逆迭代法计算结构的前9阶特征值后求出对应的固有频率。将平面结构划分为（59×29）个四边形网格，有限元部分采用上文中的四节点矩形单元，对梁的左下和右下两个节点施加相应的位移约束，计算结果，如表1所示。

表1 eigs函数与矢量迭代法求得前9阶频率Tab.1 The First Nine Frequencies from the Eigs Function and Vector Iterative Method

由表1可以看出此算例结构的基频约为59.95Hz，采用矢量迭代法与eigs函数所得结果相差不大，但是在编写调试程序时，直接调用MATLAB中eigs函数明显便捷，所以后面在MATLAB中进行拓扑优化时，就采用直接调用函数的方法在每一步迭代中求解变化后结构的基频。

算例二：

简支梁结构和材料属性与算例一中一致，如图3所示。在梁的下边界中点P点处作用有一个竖直向下的集中载荷F=2100N，对此结构进行拓扑优化设计使其质量得到减轻，同时保证p点竖直位移不大于0.5mm，结构基频大于65Hz。

图3 平面简支梁受力示意图Fig.3 Schematic Diagram of the Force of a Simple Supported Beam

依据式（1）建立优化模型如下：

将平面结构依旧划分为（58×29）个四边形网格，为避免应力集中的影响，在程序中将集中载荷分散在P点及其附近几点上，经93步迭代计算后所得拓扑结构及结构基频、结构总质量和P点位移在迭代过程中的变化情况，如图4～图7所示。

图4 优化后所得拓扑结构Fig.4 Topology Structure After Optimization

图5 结构总质量在优化过程中的变化Fig.5 Changes of the Total Quality of Structure in the Process of Optimization

图6 结构基频在优化过程中的变化Fig.6 Changes of the Fundamental Frequency of Structure in the Process of Optimization

图7 P点位移在优化过程中的变化Fig.7 Changes of the Displacement of Point P in the Process of Optimization

从图中可以看出，优化后所得的拓扑结构相比原始结构784kg约减轻至445.1kg，同时拓扑优化后结构的基频约为65.73Hz，相比算例一中的初始结构有所提高，且优化后的结构在载荷F的作用下p点位移为-0.477mm，也满足位移约束要求。

8 结论

（1）以上算例是针对平面连续体结构的拓扑优化设计，显然，更换有限元中使用的单元类型为8节点长方体单元，并对拓扑优化过程中的迭代算法进行修改，亦可推广到三维连续体结构的拓扑优化设计。（2）基于有限元理论，推导了平面4节点矩形单元的刚度矩阵，用于拓扑优化迭代过程中节点位移的求解；同时分别采用MATLAB函数和矢量迭代法求解动力学特征方程，获得了结构的固有频率，完成了静动态联合拓扑优化约束条件的求解，从而保证优化后的结果既满足结构的刚度需求，又使得结构的动态性能得到改善。（3）基于变密度法实现了静动态联合拓扑优化，能够高效的求解计算获得相应的拓扑结构，实现了重量的减轻，提高了材料利用率，践行了结构轻量化设计理念，而且获得的拓扑结构在设计阶段具有重要的指导意义。