循环对称结构的多尺度拓扑优化方法

(西北工业大学 航宇材料结构一体化设计与增材制造装备技术国际联合研究中心，陕西省空天结构技术重点实验室，西安 710072)

1 引 言

循环对称结构是一类在极坐标系下沿圆周方向循环对称的特殊周期结构，如图1(a)所示的夹芯点阵薄壁圆筒结构，其轻量化高性能设计在航空航天与机械领域具有重要应用价值。同时，这类结构也广泛存在于骨骼和植物等自然界生物材料中[1-3]。然而，过去几十年建立的结构多尺度优化方法主要针对笛卡尔坐标系下周期结构，如图1(b)所示的夹芯点阵平板结构。其特点为所有点阵微结构单胞均为矩形单胞，可由某一单胞经几何平移操作获得。因此，所有单胞微结构均可假设具有相同的等效力学性能，并通过均匀化方法计算结构的宏观力学响应。由此可见，均匀化方法不仅是计算单胞微结构等效力学性能的关键步骤，也是关联微结构与宏观性能和实现多尺度优化的桥梁。目前，基于微结构均匀化方法，相继发展了固体各向同性惩罚方法(SIMP)[4,5]、渐进结构优化方法(ESO)[6,7]和水平集方法(LSM)[8]等微结构拓扑优化方法，同时也促进了功能梯度材料以及考虑尺度效应的微结构多目标优化与拓扑优化研究[9-12]。

图1 夹芯点阵

本文研究的循环对称结构多尺度拓扑优化问题具有一定的代表性。与平移周期结构相比，可以发现当循环对称结构的半径趋于无穷大时，该结构退化为平移周期结构。因此，如何建立循环对称结构的多尺度优化方法，从理论上实现两类周期结构的有机统一，不仅具有重要的工程应用价值，也是本文的科学意义所在。

然而，目前循环对称结构的拓扑优化仍停留在单一尺度上，即直接对周期微结构[13-18]进行优化。由于缺乏微结构力学性能的均匀化等效，计算效率受到极大挑战，难以实现大规模循环对称结构的多尺度优化。Zhang等[19]尝试将针对平移周期微结构等效的均匀化方法直接应用到循环对称微结构，并未考虑随半径位置以及θ方向变化时单胞微结构性能的不同。Chatzigorgiou等[20,21]提出的改进式均匀化方法可用于循环对称周期微结构均匀化等效，并通过矩阵变换得到笛卡尔坐标系下相应的均匀化等效性能。然而，该方法需要对不同半径的单胞进行均匀化，其高昂的计算成本极大限制了应用效果，导致相关工作仍停留在均匀化等效层面上，循环对称结构多尺度拓扑优化研究仍处于起始阶段。

2 循环对称结构均匀化方法

2.1 循环对称单胞均匀化等效弹性矩阵计算

图2 循环对称结构

(1)

(2)

表1比较了循环对称微结构单胞和平移周期微结构单胞的均匀化方法差异。

表1 循环对称和平移周期微结构的均匀化方法比较Tab.1 Homogenization methods of cyclic symmetrical structure and translational periodic structure

(3)

式中M为二阶张量的旋转矩阵。

(4)

2.2 循环对称单胞均匀化等效性能的不变性

由表1可知，对于循环对称单胞，由于微观算子包含宏观半径r，单胞微结构均匀化等效性能与其所在位置对应的宏观半径变量相关，理论上应对不同宏观半径的单胞分别进行均匀化等效计算。

图3 两个相似扇形单胞

(5)

(6)

(7)

比较式(6,7)可得

(8)

(9)

式中本文引入单胞无量纲特征参数λ用于描述扇形单胞中心处弧长与边长之比。其物理意义为，若两相似单胞的特征参数λ相同，即使宏观半径不同，其均匀化等效性能仍然相同。由此可知，图2所示的循环对称结构，若单胞微观边长与其宏观半径呈正比例关系，则其均匀化等效性能保持相同(推导忽略)。

2.3 循环对称单胞均匀化等效性能的可映射性

受文献[23，24]关于平移周期单胞存在的映射关系启发，本文基于特征参数λ建立了循环对称单胞等效性能计算的映射关系。其重要意义为，对于不同λ的扇形单胞，无需通过式(1)每次重新划分网格计算其等效性能，只需在标准单胞(λ=1)上计算求解等效性能。其中的基体材料属性以及最终单胞等效性能通过相应映射计算获得，由此极大简化了不同λ取值单胞的均匀化过程。

图4所示为标准单胞(λ=1)和任意λ取值单胞，参数分别使用上标^和～区分。存在如下几何关系

图4 单胞映射变化

(10)

通过理论推导，任意λ取值单胞的均匀化等效性能为

(11)

式中Ti j为映射张量的分量，映射矩阵T为

(i,j=1,2)(12)

(13)

(14)

(15)

2.4 循环对称单胞与平移周期单胞的统一性

考虑图5所示的循环对称单胞与平移周期单胞。

图5 不同周期单胞

平移周期单胞尺寸分别满足

(16)

(17)

(18)

即表1两坐标系下微观算子等价。因此，当r*→ ∞时，循环对称单胞与平移周期单胞均匀化等效性能等价。

3 单胞微结构参数化建模与等效

性能插值计算

3.1 单胞微结构的特征驱动参数化建模

基于文献[25,26]的前期研究工作，本文首次将特征驱动拓扑优化方法应用于循环对称单胞微结构的多尺度拓扑优化。图6所示为超椭圆特征，其水平集函数为

图6 超椭圆特征水平集函数及其相关设计变量

(19)

(20)

(21)

图7(左)所示为具有最大体积模量和最大剪切模量的两种平移周期单胞微结构[27]。对其分别采用相同数目的超椭圆特征进行参数化建模，则单胞微结构包含的超椭圆设计变量集合为s。

图7 两种典型的微结构

(22)

式中(θc i,rc i,ac i,wc i和lc i)T代表第i个超椭圆设计参数。生成的单胞微结构整体水平集函数可通过所有超椭圆特征水平集函数的布尔运算获得

(23)

采用以上最大体积模量和最大剪切模量单胞微结构分别作为初始和最终微结构构型，则任意中间态单胞微结构构型可通过设计变量插值确定。

sκ=sinitial+κ(sfinal-sinitial)

(24)

图8 由超椭圆特征组成的单胞微结构

(25)

极坐标系下，单胞微结构的体积分数计算式为

(26)

式中H(·)为海维赛德函数。

引入超椭圆宽度比例变量ω用于改变单胞体积分数，则有

sω=ωs,ω=diag(1,1,1,ω,1)

(27)

(28)

图9 两种不同演变形式的参数化微结构

(29)

(30)

3.2 参数化微结构等效力学性能的插值计算

单胞微结构的等效力学性能不仅取决于微结构构型和材料用量，还与其弧边比特征参数λ紧密相关。为了定量地描述参数化微结构的等效力学性能与微结构控制参数的依赖关系和减少重复计算成本，本文基于函数拟合插值模型[21]构造了等效力学性能与微结构控制参数的近似表达式，具体采用三元四次多项式插值方式描述三个控制参数与等效力学性能的关系。

(31)

式中n=4，βI J L为对应系数。

图11给出了均匀化等效性能C11的剖视图与λ=1透明截面处的截面视图。圆圈为原始样本点均匀化等效数据，均位于拟合曲面附近，具有良好的拟合精度。

图11 均匀化等效性能C11

4 密度与插值参数κ空间分布的B样条参数化建模与优化模型

4.1 基于B样条参数化的多尺度设计框架

(32)

[sinitial+κ(sfinal-sinitial)]

(33)

如图12所示，由于相邻单胞在公共边界上共享相同参数变量，B样条可以自动保证微结构的几何连续性，无需施加任何额外约束。而传统密度法由于不同单胞微结构的变量离散化定义，即使对微结构进行细分也无法从本质上改变相邻微结构之间连续性差的问题。

图12 微结构连接光顺性的比较

4.2 循环对称结构多尺度拓扑优化

研究给定体积约束下结构柔顺度最小化问题。对应的数学模型为

min.C=UTKU

(34)

式中λm为中心线半径处的弧边比特征参数。

图13为基于B样条参数化的循环对称结构多尺度优化框架。

图13 循环对称结构多尺度优化流程

5 数值算例

5.1 曲梁多尺度拓扑优化

图14所示曲梁左端固支，右下端加载集中力F=1 N。轴线长l=600 mm，宽度h=200 mm，r0=1200/πmm，固体材料杨氏模量E=1 MPa,泊松比为0.3。假设材料用量体分比为50%，曲梁设计域在极坐标系下划分为20×60个单胞，并填充如图9(左)所示的参数化微结构，中心线半径处λm=1。当曲梁半径R=r0时，曲梁扇形角为90°。随着R的增大，曲梁的曲率逐渐变小。当R→ ∞时，曲梁的曲率趋近于零并退化为悬臂梁，其均匀化等效性能也与传统平移周期微结构一致。

图14 受集中力的曲梁

针对上述三种情况，将不同R的曲梁均置于相同B样条参数域中进行拓扑优化，其控制点位置不随R的变化而发生改变，如图15所示。

图15 不同半径曲梁置于同一个B样条参数域

为了验证该方法的正确性，说明平移周期结构是循环对称结构的特例，取R=1×104r0，使得曲梁接近于悬臂梁的形式。此时，曲梁的弧度约为1.57×10-4，几乎可忽略不计。多尺度优化结果和迭代曲线如图16和图17所示。前者与VCUT水平集方法[29]一致，有效验证了循环对称结构多尺度拓扑优化方法的正确性。

图16 R=1×104r0的曲梁多尺度优化结果

图17 R=1×104r0时曲梁柔顺度和体积分数的迭代曲线

表2给出了不同R下曲梁的多尺度优化结果。

可以发现，材料主要沿上下表面的加载路径分布，且局部微结构构型与其主应力方向相关，曲梁的上下表面构型参数趋向于0，以横竖杆为主。在上下主传力路径之间，存在一定的过渡结构。从 图18 可以看出，随着半径的增加，柔顺度逐渐增加；当半径趋于无穷大时，逐渐收敛到悬臂梁的柔顺度值。

图18 不同R值曲梁结构柔顺度

5.2 圆盘多尺度拓扑优化

该问题着重研究不同循环对称单胞排列方式、弧边比参数以及曲率半径对结构性能的影响。

如图19(a)所示，圆环内圈固定，圆环外圈施加切向载荷，外径与内径之差2(R0-r0)=50 mm。假设材料体积分数限制为40%，基体材料杨氏模量E=10 GPa，泊松比为0.3。图19(b，c)分别给出了等间距排列单胞和等比例排列单胞。最初将圆盘划分为6×54个单胞，优化过程中径向层数保持不变，圆周方向的数目与λm相关。本算例假定r0分别为15 mm,25 mm,50 mm,75 mm，λm的初始值分别为0.54,1.01,1.69,2.36。

图19 圆盘问题

对于等间距排列单胞，位于宏观坐标(r,θ)的微结构弧边比特征参数λ与中心线弧边比特征参数λm存在如下关系，

λ=2r/(r0+R0)λm

(35)

等间距排列单胞在不同内径下的优化设计结果如图20所示。可以看出，(a) 微结构插值构型参数κ始终接近1，形成X形的微结构。(b)λm总是收敛于1附近，表明具有相同弧长和边长的扇形单胞更有利于结构优化。(c) 无论内径长度如何，材料均在内径处堆积，微结构体积分数从内径到外径逐渐减小。 (d) 内径较小时，内径和外径的体积分数变化很大。随着内径的增加，内径与外径的体积分数差逐渐减小。

图20 等间距排列单胞在不同内径下的优化设计结果

图21 内径趋于无穷大时的结构演化

同样，当r0=15 mm,25 mm,50 mm,75 mm且单胞以等比例方式排列时，可以得到如图22所示与等间距排列单胞相似的优化结果。该排列方式确保参数λ在圆的所有位置均等于中心线处弧边比参数λm。

图22 等比例排列单胞在不同内径下的优化设计结果

表3比较了相同内径下不同单胞排列方式对柔顺度的影响。柔顺度结果表明，等比例排列单胞明显优于等间距排列单胞。当内径较小时，这两种排列结构的柔顺度差异很明显，并在r0=15 mm时二者柔顺度降幅达到了6.37%。随着内径的增大，二者柔顺度的差异逐渐减小。当r0→ ∞时，两种排列方式之间的差异将不再明显。因此，在环向剪切载荷作用下，等比例排列比等间距排列具有更大的优势，在工程中对于受切向载荷的圆筒结构可优先考虑使用等比例排列单胞设计。

表3 不同半径不同排列方式柔顺度差异Tab.3 Differences in compliance under different radii and different arrangements

(续表)

6 结 论

本文开展了循环对称结构多尺度拓扑优化方法研究，提出了用于表征等效性能映射计算的弧边比特征参数，阐明了循环对称单胞均匀化等效性能的不变性、可映射性以及与平移周期单胞的统一性，建立了微结构单胞等效性能的三元拟合插值模型与B样条参数化多尺度设计框架。

算例结果表明，对于受剪切载荷作用的循环对称结构，等比例单胞排列设计优于等间距单胞排列设计，且当结构半径趋于无穷大时，循环对称结构退化为平移周期性结构。

参考文献(References):

[1] Liu Z Q,Meyers M A,Zhang Z F,et al.Functional gradients and heterogeneities in biological materials:Design principles,functions,and bioinspired applications[J].ProgressinMaterialsScience,2017,88:467-498.

[2] Osorio L,Trujillo E,van Vuure A W,et al.Morphological aspects and mechanical properties of single bamboo fibers and flexural characterization of bamboo/epoxy composites[J].JournalofReinforcedPlasticsandComposites,2011,30(5):396-408.

[3] Chen W J,Zheng X N,Liu S T.Finite -element-mesh based method for modeling and optimization of lattice struc-tures for additive manufacturing[J].Materials,2018,11(11):2073.

[4] Bendsøe M P,Sigmund O.Material interpolation schemes in topology optimization [J].ArchiveofApp-liedMechanics,1999,69(9-10):635-654.

[5] Sigmund O.Tailoring materials with prescribed elastic properties[J].MechanicsofMaterials,1995,20(4):351-368.

[6] Xia L,Breitkopf P.Recent advances on topology optimization of multiscale nonlinear structures[J].ArchivesofComputationalMethodsinEngineering,2017,24(2):227-249.

[7] Radman A,Huang X,Xie Y M.Topology optimization of functionally graded cellular materials[J].JournalofMaterialsScience,2013,48(4):1503-1510.

[8] Li H,Luo Z,Gao L,et al.Topology optimization for functionally graded cellular composites with metamaterials by level sets [J].ComputerMethodsinApp-liedMechanicsandEngineering,2018,328:340-364.

[9] 邱克鹏,张卫红.功能梯度材料结构拓扑优化设计研究[J].西北工业大学学报,2010,28(6):851-857.(QIU Ke -peng,ZHANG Wei-hong.Improving topo -logical design of FGM (functionally graded material) structure[J].JournalofNorthwesternPolytechnicalUniversity,2010,28(6):851-857.(in Chinese))

[10] Wang Y Q,Zhang L,Daynes S,et al.Design of graded lattice structure with optimized mesostructures for additive manufacturing [J].Materials&Design,2018,142:114-123.

[11] 孙士平,张卫红.多相材料微结构多目标拓扑优化设计[J].力学学报,2006,38(5):633-638.(SUN Shi-ping,ZHANG Wei-hong.Multiple objective topology optimal design of multiphase microstructures[J].ChineseJournalofTheoreticalandAppliedMe-chanics,2006,38(5):633-638.(in Chinese))

[12] 张卫红,骆金威,戴高明,等.周期性多孔材料等效剪切模量与尺寸效应研究[J].力学学报,2011,43(1):144-153.(ZHANG Wei-hong,LUO Jin-wei,DAI Gao-ming,et al.Numerical predictions of effective shear modulus and size effect for periodic cellular materials [J].ChineseJournalofTheoreticalandApp-liedMechanics,2011,43(1):144-153.(in Chinese))

[13] Tang Y L,Dong G Y,Zhao Y F.A hybrid geometric modeling method for lattice structures fabricated by additive manufacturing[J].TheInternationalJournalofAdvancedManufacturingTechnology,2019,102(9-12):4011-4030.

[14] Dedè L,Borden M J,Hughes T J R.Isogeometric analysis for topology optimization with a phase field model[J].ArchivesofComputationalMethodsinEngineering,2012,19(3):427-465.

[15] Moses E,Fuchs M B,Ryvkin M.Topological design of modular structures under arbitrary loading[J].StructuralandMultidisciplinaryOptimization,2002,24(6):407-417.

[16] Gao T,Zhang W H.Topology optimization involving thermo -elastic stress loads[J].StructuralandMultidisciplinaryOptimization,2010,42(5):725-738.

[17] Zuo Z H,Xie Y M,Huang X D.Optimal topological design of periodic structures for natural frequencies[J].JournalofStructuralEngineering,2011,137(10):1229-1240.

[18] Maharaj Y,James K A.Metamaterial topology optimization of nonpneumatic tires with stress and buckling constraints[J].InternationalJournalforNumericalMethodsinEngineering,2020,121(7):1410-1439.

[19] Zhang H K,Wang Y G,Kang Z.Topology optimization for concurrent design of layer-wise graded lattice materials and structures[J].InternationalJournalofEngineeringScience,2019,138:26-49.

[20] Chatzigeorgiou G,Efendiev Y,Lagoudas D C.Homo -genization of aligned “fuzzy fiber” composites[J].InternationalJournalofSolidsandStructures,2011,48(19):2668-2680.

[21] Wang C,Zhu J H,Zhang W H,et al.Concurrent topology optimization design of structures and non-uniform parameterized lattice microstructures[J].StructuralandMultidisciplinaryOptimization,2018,58(1):35-50.

[22] Xu Z,Zhang W H,Gao T,et al.A B-spline multi-para-meterization method for multi-material topology optimization of thermoelastic structures[J].Struc-turalandMultidisciplinaryOptimization,2020,61(3):923-942.

[23] Chatzigeorgiou G,Efendiev Y,Charalambakis N,et al.Effective thermoelastic properties of composites with periodicity in cylindrical coordinates[J].InternationalJournalofSolidsandStructures,2012,49(18):2590-2603.

[24] Li Q H,Xu R,Liu J,et al.Topology optimization de -sign of multi-scale structures with alterable microstructural length-width ratios[J].CompositeStructures,2019,230:111454.

[25] Zhou Y,Zhang W H,Zhu J H,et al.Feature -driven topology optimization method with signed distance function[J].ComputerMethodsinAppliedMecha-nicsandEngineering,2016,310:1-32.

[26] Xu Z,Zhang W H,Zhou Y,et al.Multiscale topology optimization using feature -driven method[J].ChineseJournalofAeronautics,2020,33(2):621-633.

[27] Xia L,Breitkopf P.Design of materials using topology optimization and energy-based homogenization app -roach in matlab[J].StructuralandMultidisciplinaryOptimization,2015,52(6):1229-1241.

[28] Liu S T,Cheng G D,Gu Y,et al.Mapping method for sensitivity analysis of composite material property[J].StructuralandMultidisciplinaryOptimization,2002,24(3):212-217.

[29] Zong H M,Liu H,Ma Q P,et al.VCUT level set method for topology optimization of functionally graded cellular structures[J].ComputerMethodsinAppliedMechanicsandEngineering,2019,354:487-505.