基于子结构的周期性结构设计方法

吴紫俊 肖人彬

摘要 针对传统基于均匀化方法的周期性结构设计宏微尺度分离问题，提出了基于有限元子结构的周期性结构设计方法。在该方法中，假设宏观结构由构型相同的多个子结构构成。在宏观层面，利用子结构凝聚构建了超单元计算模型;在微观层面，采用SIMP优化方法建立了子结构内部单元的拓扑构型设计方法。通过反求超单元模型计算结果，实现了宏微观结构的耦合计算。在耦合计算过程中，由于超单元之间具有公共的部分边界，使得子结构构型间具有良好的材料连通性。同时，子结构划分源于宏观设计域的有限元网格，确保了宏微结构间具有明确的尺寸比例关系。最后，通过两类梁的计算实例，验证了所提方法的有效性。

关 键 词 子结构;周期性;微结构;拓扑优化

中图分类号 TB383.3     文献标志码 A

0 引言

周期性结构设计是现代结构设计的重要分支[1]，已被广泛应用于多个领域，如航空[2]、汽车[3]等的结构设计。周期性结构不仅承载能力强，重量轻[4-5]，还具有隔热[6]、隔音[7]等多功能物理特性，因其独特的结构组成形式和良好的设计加工性能，展现出了广阔的研究和应用前景[8-9]。周期性结构设计是一种多尺度的拓扑优化设计，方法上可分为2类：一类是基于均匀化的设计方法[10];另一类是基于宏观周期性约束的设计方法[11]。

均匀化方法设定材料的宏观结构由微观尺度的单胞周期性拓展而形成。周期性结构的单胞具有周期性应力和应变边界条件，从而获取单胞内材料的最优分布，实现不同结构性能的周期性单胞设计[12]。Vogiatzis等[13]基于水平集法提出了单相和多相负泊松比材料的设计方法，通过施加周期性的Dirichlet位移边界条件得到周期性微结构，并采用最速下降法构建微结构内部构型。Fantoni等[14]提出了一种多场渐近均匀化方法，用以分析具有周期性微结构的压电材料，利用周期摄动函数考虑具有代表性的微结构非均匀性的影响。Hunag等[15]将双向进化结构优化方法推广到具有最佳粘弹性特性的复合材料微观结构设计中，基于均匀化方法计算复合材料的性能參数，通过BESO方法实现单胞内材料的重新分配。贾娇等[16]研究了宏观传导条件对材料微结构的影响，并提出了基于宏观传导条件的周期性结构传热材料的研究，构建了基于变密度法周期性传热结构模型。杜义贤等[17]基于能量均匀化方法研究了抗剪性能更优的周期性点阵结构，并提出了基于经典胞元理论的周期性点阵结构的面内等效切变模量的解析表达式。

在均匀化方法中，微观单胞与宏观结构设计是相互耦合的设计过程，对微观结构的宏观有效性材料进行评价的同时，宏观结构的计算结果也影响微观结构的构型优化[18]。采用均匀化方法进行周期性结构设计，其尺度分离假设会带来2个问题：一个是微结构单胞间材料的连通性[19]，其微结构单胞优化是一个独立的优化过程，宏观位移或应力的变化往往引起微结构间材料不连通，使得所优化的结果只有理论构型，不具备制造性[20];另一个是均匀化方法假设微结构与宏观结构的尺寸比例约为10-9，当该比例接近于1时往往得不到合适的优化构型[21-22]。

宏观周期性约束的设计方法中，宏观结构被均匀地划分成若干个具有特定长度比例大小的单元胞体。每个单元胞体同一位置上的单元具有相同的材料密度值和灵敏度值[23]。优化变量为单元胞体中每个单元的材料密度。与均匀化方法相比，随着单元胞体的细化，其优化结果逐渐近似于均匀化方法设计结果。在基于BESO方法的周期性结构设计中[24]，需要在整个设计域中求解具有周期性边界约束的方程，存在计算量大的问题[25]。

为解决均匀化方法的尺度分离问题和宏观周期性约束设计方法中计算量过大的问题[26]，本文在有限元子结构方法的基础上，提出了利用子结构凝聚构建宏观结构的超单元计算模型，采用SIMP优化方法[27-28]建立子结构内部单元拓扑构型设计的周期性结构设计方法。

1 子结构凝聚及周期性结构优化问题定义

1.1 子结构凝聚

根据有限元子结构的定义，有限元网格模型可以重新划分为多个子结构[29]，如图1所示。

优化目标是宏观周期性结构的最小化柔度;设计变量为周期性子结构中包含的n个单元的相对密度[ρj];[K]、F分别为子结构对应超单元组装的整体刚度矩阵与外载荷向量，[Ubb，i]是第i个超单元的位移向量，[V]是当前周期性子结构所含的材料体积，N为宏观结构的子结构总数。

2 基于子结构的周期性结构优化模型

2.1 宏观尺度模型

在周期性结构的宏微结构优化中，宏观结构的刚度矩阵由式（5）建立。由于省略了内部节点，可提高其计算效率。而子结构的内部节点的位移可根据式（2）获得。此时，图1中的子结构模型的刚度矩阵可表示为

此时，所优化的结构将分布在整个宏观设计域中。在宏观结构计算过程中，其构建刚度矩阵的周期性子结构均有相同构型，对应的每个超单元凝聚矩阵均相同，即[K*bb，t=K*bb，s，t≠s] 。

2.2 子结构优化模型

在子结构构型的优化中，需要增加惩罚因子p对拓扑优化设计变量的密度进行惩罚，使其收敛于设定的密度上下界，从而获得子结构的“黑白”设计[30]：

式中，惩罚因子[p≥1]表明具有中间密度的子结构[0<ρj<1]的刚度矩阵会被惩罚。该惩罚因子值越大，其优化结构杆系特征越明显。式中，n表示该子结构内部的有限元单元总数，[Κ0]为实体材料下的有限元单元刚度矩阵，如图2所示。

与传统的SIMP方法不同，此时的惩罚因子p对子结构刚度矩阵的惩罚，不是凝聚后的超单元对应的刚度矩阵，而是子结构内所有单元对应的有限元单元刚度矩阵。当子结构构型优化结束后，需对其进行二次凝聚，从而在宏观结构中计算每个子结构边界的位移，实现宏微结构的下一次迭代计算。

2.3 灵敏度计算及变量迭代

根据所建立的子结构代理优化模型，上述基于子结构的周期性结构设计问题可重新定义为

式中：[κ]为拉格朗日因子，该因子可通过分半算法强制使设计域中的材料满足体积约束;体积导数[?V?ρj=vj]为每个子结构的体积约束，在此与子结构当前的材料体积相等。

在优化迭代计算中，通过OC方法进行变量更新，每个子结构的初始体积约束均相等，整个设计域在整个优化过程中的材料体积与初始设置的体积约束相等。即材料体积始终不变，优化过程只在周期性子结构中寻找材料的最优分配比例，达到最优设计。对应优化框架如图3所示。

3 数值计算实例

3.1 双端固支梁结构设计

以双端固支梁的拓扑优化设计为例[31]，设计域大小为[L1×L2=2×1]，其左右两侧固定，设计域中央承受[F=1]的垂直向压力如图4所示。所选用材料的弹性模量为[E0=1]，柏松比[υ0=0.3]。

把宏观设计域划分为3种不同数量的子结构：[Nx×Ny=2×1]，[Nx×Ny=4×2]和[Nx×Ny=8×4];子結构分别设置为[30×30]，[50×50]两类不同单元数的单元集合。根据宏观模型的定义，在优化过程中，设体积约束设置为0.3，由于宏观结构中没有设计变量，因此每个周期性子结构的体积约束也为0.3。在周期性子结构的优化中，当惩罚因子设置为1，所优化的周期性子结构拓扑结构中存在大量的中间密度单元，没有明显的杆系结构特征。为了得到更为精细的局部结构特征，在此把惩罚因子设置[p=3]。优化迭代步长设置为0.02，为避免周期性子结构中出现棋盘格现象，灵敏度的过滤半径设置为子结构内部单元长度的1.5倍，当相邻两次目标函数值相对误差小于0.000 1时所得拓扑构型即为最优结构。其优化结果如表1所示。

从优化结果中，可推断出以下结论：1）相同的宏观结构子结构划分，优化时所选取的子结构大小不同，其最终的优化拓扑构形也不同，其宏观结构的柔度值随着子结构单元数增多而减小;2）在同一大小的子结构下，宏观结构子结构划分不同，最终优化的周期性子结构构形基本相似，其宏观结构的柔度值随着子结构数量增多而逐渐增大。当子结构所含单元数增多时，意味着设计域的有限元网格更为细化，所得到的宏观优化结构特征越精细，其优化结构的柔度值越小。

从优化结果看，微结构边界间具有良好的材料连通性，改变了文献[18]中由于应变场突变引起的微结构间材料不连通的情况，为其结构加工制造提供了模型基础。同时，设计过程中宏观结构与微观构型间具有明确的尺寸比例关系，保证了设计出的结构性能与实际制造出的结构性能的一致性，避免了均匀化方法中因尺度分离引起的结构设计与制造的性能误差。

3.2 悬臂梁结构设计

在悬臂梁的拓扑优化设计中，设置设计域大小为[L1×L2=2×1]，其左侧固定，设计域右上方角点承受[F=1]的垂直向压力如图5所示。所选用材料的弹性模量为[E0=1]，柏松比[υ0=0.3]。

把宏观设计域划分为3种不同数量的子结构：[Nx×Ny=2×1]，[Nx×Ny=4×2]和[Nx×Ny=8×4];子结构设置为[40×40]的单元集合。为验证惩罚因子对优化构型的影响，分别设置惩罚因子[p=1]和[p=2]。灵敏度的过滤半径设置为子结构内部单元长度的1.2倍，优化迭代步长设置为0.02，当相邻两次目标函数值相对误差小于10-4时所得拓扑构型即为最优结构。其优化结构如表2所示。

从优化结果中可以看出：当惩罚因子设置为1，所优化的周期性子结构拓扑结构中存在大量的中间密度单元，没有明显的杆系结构特征，随着惩罚因子的增大，其杆系特征逐渐明显。同时，最优宏观构型的柔度值随着惩罚因子和宏观子结构划分数的增大而增大。从优化结构来看，所得到的宏观结构及其对应的微观构型，其材料具有很好的连通性，因此所得构型具有良好的可制造性。

需要注意的是由于SIMP方法中惩罚因子的影响，使得子结构内部有限元单元的相对密度不是离散的0或1。尽管通过子结构方法可获得具有良好材料连通性的优化结构，但是所产生的中间密度单元给制造带来了困难。

4 结论

本文提出了一种基于子结构的周期性微结构拓扑优化方法，解决了传统均匀化方法中尺度分离引起的微结构间材料不连通和结构设计与制造间性能误差的问题，为设计制造一体化提供了一种参考方法。本文研究结论可归纳如下。

1）通过子结构凝聚构建宏观结构的超单元计算模型，把宏观结构尺寸信息带入到了微结构中，同时，通过引入超单元间的共同部分边界确保了微结构设计过程中材料的连通性。

2）优化计算结果表明，宏观子结构划分与子结构内部有限元网格越多，其优化构型的杆系特征越明显。由于采用的变密度法SIMP，当惩罚因子取值较小时，得到大量的中间密度单元。

3）采用SIMP方法，不可避免存在中间密度单元，给实际加工制造带来困难。同时，所设计的宏微结构模型边界存在台阶现象，因此需要在可加工性方面做进一步的改进。

参考文献：

[1]    ZHU J，ZHANG W，BECKERS P. Integrated layout design of multi-component system[J]. International journal for numerical methods in engineering，2009，78（6）：631-651.

[2]    LIU S，HU R，LI Q，et al. Topology optimization-based lightweight primary mirror design of a large-aperture space telescope[J]. Applied optics，2014，53（35）：8318-8325.

[18]  XIA L，BREITKOPF P. Multiscale structural topology optimization with an approximate constitutive model for local material microstructure[J]. Computer Methods in Applied Mechanics and Engineering，2015，286：147-167.

[19]  SIVAPURAM R，DUNNING P D，Kim H A. Simultaneous material and structural optimization by multiscale topology optimization[J]. Structural and multidisciplinary optimization，2016，54（5）：1267-1281.

[20]  LI H，LUO Z，GAO L，et al. Topology optimization for functionally graded cellular composites with metamaterials by level sets[J]. Computer Methods in Applied Mechanics and Engineering，2018，328：340-364.

[21]  ZHU Y，LI S，DU Z，et al. A novel asymptotic-analysis-based homogenisation approach towards fast design of infill graded microstructures[J]. Journal of the Mechanics and Physics of Solids，2019，124：612-633.

[22]  LI Q，XU R，LIU J，et al. Topology optimization design of multi-scale structures with alterable microstructural length-width ratios[J]. Composite Structures，2019，230：111454.

[23]  XIE Y M，ZUO Z H，Huang X，et al. Convergence of topological patterns of optimal periodic structures under multiple scales[J]. Structural and Multidisciplinary Optimization，2012，46（1）：41-50.

[24]  XIA L，XIA Q，HUANG X，et al. Bi-directional evolutionary structural optimization on advanced structures and materials：a comprehensive review[J]. Archives of Computational Methods in Engineering，2018，25（2）：437-478.

[25]  FU J，XIA L，GAO L，et al. Topology optimization of periodic structures with substructuring[J]. Journal of Mechanical Design，2019，141（7）：071403.

[26]  HOU Y，ZHAO Q，Sapanathan T，et al. Parameter identifiability of ductile fracture criterion for DP steels using bi-level reduced surrogate model[J]. Engineering Failure Analysis，2019，100：300-311.

[27]  BENDS?E M P，SIGMUND O. Material interpolation schemes in topology optimization[J]. Archive of applied mechanics，1999，69（9-10）：635-654.

[28]  SIGMUND O. A 99 line topology optimization code written in Matlab[J]. Structural and multidisciplinary optimization，2001，21（2）：120-127.

[29]  WU Z，XIA L，Wang S，et al. Topology optimization of hierarchical lattice structures with substructuring[J]. Computer Methods in Applied Mechanics and Engineering，2019，345：602-617.

[30]  BENDS?E M P，SIGMUND O. Topology optimization：theory，methods and applications[M]. Berlin：Springer-Verlag，2003.

[31]  ALEXANDERSEN J，LAZAROV B S. Topology optimization of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner[J]. Computer Methods in Applied Mechanics and Engineering，2015，290：156-182.

[責任编辑    杨    屹]