引 言
轻量化设计如今是力学界、数学界与材料界的研究热点之一.在各种条件下寻求最优的结构拓扑形状和相应的微观材料分布，对研究者来说是非常有挑战性和吸引力的.
自20世纪70年代以来，均匀化理论开始被用来寻找材料的有效性质［1-3］.1988年Bendse等［4］利用均匀化理论对多孔材料进行结构优化设计，为均匀化理论应用于结构优化设计奠定了基础.1994年Sigmund［5］提出逆均匀化理论，根据逆均匀化理论，可以设计出不同性能的微观材料，为微观材料设计奠定了理论基础.2002年Rodrigues等［6］提出层级结构和材料的优化设计，即在确定性荷载作用下，寻求使结构刚度最大的宏观结构和微观材料分布，其中宏观结构由不同的微观材料构成，但其计算规模较大，工业制备较为困难.2007年阎军［7］的研究结果表明宏观结构上只有集中力作用时，最优微观单胞构型往往是各向同性的实体材料.2008年Liu等［8］对层级结构优化进行改进，提出了并发优化(concurrent optimization)理论，即宏观结构所用的微观材料均相同，既减小了计算规模又易于制备和应用.Coelho等［9］根据层级优化理论得到了三维的微观结构单胞.2011年Takezawa等［10］根据特征分析理论建立了不确定荷载作用下的结构优化模型，通过大量算例验证了此模型的合理性.2015年Guo等［11］提出了作用有不确定性荷载的超轻结构的多尺度结构优化，其结构优化为双层优化列式，上层求解宏观结构，下层转换为SDP(semi-definite programming)问题来求解，但是迭代收敛较慢.
本文针对存在的问题，基于特征值分析，利用并发结构优化理论，建立了双尺度超轻结构的多尺度优化模型.该优化方法被定义为一个最小柔度问题，即求使宏观结构柔度最小的宏观拓扑结构和微观材料分布.本文中作用的不确定性荷载被定义为其Euclid范数等于1，此时不确定性荷载可以视为名义值为零的特殊形式.此外，当荷载的定义域是线弹性定义域中一个有限小的子域时，可以合理利用荷载的子域而不是利用全局定义域来计算柔度.在这种情况下，全局荷载可以用局部荷载来代替，相应的全局位移同样用局部位移来代替.这样，用全局荷载定义的最小柔度问题可以利用“凝聚”方法转换为用维度相对较小的局部荷载来进行计算.“凝聚”方法是控制领域用于减少大型控制系统或状态方程的维度，这种思想在优化设计中也同样适用.如果原始问题是非常巨大的，利用“凝聚”概念重新定义荷载和位移向量，将不确定性荷载和位移组成的对称矩阵缩减到一个较小的维度，“凝聚”可以将原优化问题转化为小而简单的优化问题，相应地就减少了计算规模，提高了计算效率.
根据Rayleigh-Ritz定理，将最小柔度问题转换为求解一个对称矩阵的最大特征值最小问题.通过最大特征值对应的特征向量矩阵可以很容易地求得最不利荷载.特别是对点荷载问题，特征值表示结构的指定点上最不利荷载.
本文所提出的优化方法可以得到使宏观结构刚度最大的宏观拓扑结构和微观材料分布，从而得到相应的最不利荷载.本文拓扑优化中密度函数的更新是基于灵敏度分析［12］和移动渐近线方法.文中通过大量算例来说明宏观结构和微观材料分布的合理性，并得到最不利荷载作用下的最优宏观拓扑结构与微观材料分布.本文提出的优化方法不需要双层优化列式，数值迭代更稳定，对初始数值依赖不大，并且可以得到最不利荷载以验证所得结果的合理性.
1 双尺度结构优化模型
假设宏观结构上作用有Euclid范数为1的不确定性荷载，宏观设计域Ω是由一个个周期性基本单胞排列构成的孔洞材料组成(如图1所示)，基本单胞是由各向同性线弹性材料构成.其中微观设计域为1×1的正方形，记为Y.以刚度最大即柔度最小为目标函数，列式如下:
其中，ρ，γ分别为宏观和微观材料密度，F为全局荷载向量，‖·‖为Euclid范数，U为全局位移向量，c表示结构的柔度，P和P MI分别为宏观和微观体积约束.上述列式表示求作用最不利荷载时使结构柔度最小的宏观和微观材料分布.
图1 由周期性孔洞材料组成的设计域
Fig.1 A design domain composed of periodic cellular units
2 结构与灵敏度分析
2.1 结构分析
作用在结构上的局部荷载用F l表示，与其对应的位移向量记为U l，则柔度可表示为
局部荷载F l和U l可以表示为
其中，H为连接局部向量和全局向量的矩阵.
局部位移向量又可表示为如下形式:
其中
C为与局部荷载向量F l和局部位移向量U l维数相同的对称矩阵.因此柔度也可以表示为如下形式:
当‖F‖ =1时，根据Rayleigh-Ritz定理，可得
由上式可知，柔度c最大问题就等于矩阵C的最大特征值问题.则式(1)中的目标函数可以表示为如下形式:
即最不利荷载作用时的柔度最小问题转化为矩阵C的最大特征值最小问题.假设最不利荷载为F lworst，则式(9)表示如下:
式(13)为最大特征值和与之对应特征向量的特征方程，F lworst可由矩阵C的最大特征值对应的特征向量表示.
宏观结构采用PAMP(porous anisotropic material with penalization)方法表示，微观结构用SIMP(solid isotropic material with penalization)方法表示，假设微观结构所用材料的弹性模量为D B，微观结构所组成的微观材料的等效弹性模量为D H，则
其中，α，β为宏观尺度与微观尺度的惩罚因子，γ，ρ分别为微观尺度和宏观尺度的材料密度.已知全局节点荷载F，宏观结构有
其中，K为宏观结构的刚度矩阵，B为应变矩阵，U为全局节点位移.
D H为连接宏观和微观尺度的桥梁，可由均匀化方法获得
其中，u MI为微观结构节点位移，b为应变矩阵.在二维问题中，I为3×3的单位矩阵，Y 为单胞面积.
计算矩阵C的式(7)中包含刚度矩阵的逆矩阵，引入伴随矩阵Z，将式(7)改为如下形式:
其中
优化列式(1)转换为如下形式:
2.2 灵敏度分析
其中，r表示宏观尺寸的第r个单元，N表示宏观单元总个数.式(20)也可写为
上式两边同时对设计变量X进行求导，得
对于式(22)，两边同时对设计变量X进行求导，得
由式(24)可知，式(25)可表示为如下形式:
由于
将式(27)代入式(26)，得
则
同理可得
2.3 优化流程图
本文优化方法具体的流程如图2所示.
图2 优化流程图
Fig.2 The flow chart of the optimization procedure
3 数值算例
3.1 矩形设计域算例
如图3所示，悬臂梁设计域是长宽比为2∶1的矩形，左端固定端约束，右侧边中点处作用有不确定性荷载.微观单胞是1×1的正方形，本算例中的中间小正方形区域初始密度设为1，其他部分密度设为0.4.宏观设计域体积约束为0.2，微观设计域的体积约束为0.4，弹性模量E设定为1，微观材料Poisson比为0.3，惩罚因子在宏观尺度和微观尺度中分别设为3和4.
图3 宏观(悬臂梁)和微观设计域
Fig.3 Macro(cantilever)and micro design domains
本算例中的微观设计域划分为20×20和40×40的正方形网格，宏观设计域划分为50×25和80×40的正方形网格，优化结果如表1所示.
从表1中可以看出，宏观尺度的网格划分不同时，宏观结构拓扑形状基本相同，宏观结构拓扑形状与竖直方向荷载作用时的优化结果类似.在微观尺度网格划分不同时，微观材料分布情况类似，均为kagome形状单胞.从图4的收敛曲线可以看出，在迭代25步之后趋于收敛，各不同网格划分下的特征值收敛曲线基本接近，即收敛的特征值相近.也就是说本文中的模型宏观和微观尺度均不依赖网格划分.
悬臂梁右侧中点处作用有竖直向上或竖直向下的荷载时，悬臂梁产生的位移最大，也就是说使悬臂梁柔度最大的最不利荷载方向为竖直向上或竖直向下.而本算例中的最不利荷载向量为［0，－1］T，所得结果与现实情况相吻合，说明本文所建模型是合理的.
表1 不同网格划分下多尺度优化的最终结果
Table 1 Multi-scale structure optimization design results of different meshings
?
图4 不同网格划分下迭代收敛曲线
Fig.4 The iteration histories of different meshings
图5 宏观(反L形梁)和微观设计域
Fig.5 Macro(reverse L-shape beam)and micro design domains
3.2 反L形设计域
宏观反L形设计域结构尺寸如图5所示，宏观设计域用0.05×0.05的正方形网格来划分，上端固定端约束，左边中点处作用有单位不确定荷载.微观单胞是1×1的正方形，划分为20×20的正方形网格.
本算例中微观设计域的初始密度分为两种:微观设计域中间部分初始密度小于其余部分时定义为初始构型1，中间部分初始密度大于其余部分时定义为初始构型2.
宏观设计域体积约束为0.2，微观设计域的体积约束为0.4，弹性模量E设定为10，微观材料Poisson比为0.3，惩罚因子在宏观尺度和微观尺度中分别设为3和4.
微观设计域初始密度不同时，宏观尺度和微观尺度优化结果如表2所示.
表2 反L形梁优化结果
Table 2 The reverse L-shape beam optimization design results
?
通过表2可知，利用本文所提多尺度优化方法来建立模型，所得的宏观尺度优化结果与有名义值的不确定荷载作用时的宏观尺度优化结果［11］类似，微观结构拓扑形式类似.单位不确定荷载可以认为是名义值为零的一种特殊形式，无名义值的不确定荷载作用下，与有名义值的不确定优化结果差别不大.
微观尺度的初始密度不同时，利用本文模型计算的宏观结构拓扑结构基本相同，微观单胞相似，由微观单胞周期性组成的微观材料分布是类似的.如图6所示，两种初始构型下最终的收敛值接近.故本文所提计算方法对微观尺度的初始密度依赖不大.单尺度传统优化列式下得到的反L形梁的最不利荷载向量为［0.466 7，－0.884 4］T，本算例所得最不利荷载向量为［0.485 3，－0.874 3］T，所得结果相近，说明由本文所提出的多尺度优化方法是合理的.
图6 不同初始密度下的特征值收敛曲线
Fig.6 The iteration histories of different initial densities
图7 宏观桥形设计域和微观设计域
Fig.7 The macro bridge design domain and micro design domains
3.3 桥形设计域
如图7所示，宏观尺度设计域为一个类桥形的设计域，划分为50×25的正方形网格;微观尺度设计域为一个正方形设计域，划分为20×20的正方形网格.宏观设计域底部共有49个节点作用有竖向荷载，则矩阵C为一个49×49的对称矩阵.宏观总体积约束为0.2，微观体积约束为0.4，宏观和微观的惩罚因子均为3，微观所用材料的弹性模量设为10，Poisson比设为0.3.
如图7所示，微观设计域初始密度值分为两类:一类为初始构型1，颜色较浅的小正方形部分密度为0.4，其余部分密度为0.45;一类为初始构型2，颜色较深的小正方形部分密度为0.45，其余部分密度为0.4.不同初始构型下多尺度优化结果如表3所示.
表3 不同初始密度下多尺度优化的最终结果
Table 3 Multi-scale structure optimization design results of different initial densities
?
为了进一步验证该方法对微观单胞初始密度的依赖性，所设初始构型的密度差异较小，得到的微观单胞相同，相应的材料分布相同，宏观结构拓扑形状也相同，最优结构所对应的特征值相近(如图8所示)，最不利荷载分布也相似(如图9所示).故本文所提双尺度优化设计方法对微观单胞初始密度依赖不大.
图8 不同初始密度下的特征值收敛曲线
Fig.8 The iteration histories of different initial densities
图9 不同初始密度下的最不利荷载分布
Fig.9 The worst load distributions of different initial densities
4 结论与展望
本文提出了基于特征值分析的多尺度结构拓扑优化模型，通过大量算例可知，单位不确定荷载与有名义值的不确定荷载作用下的优化结果类似，单位不确定荷载可以认为是名义值为零的一种特殊形式.
本文所提出的双尺度优化设计方法在宏观和微观尺度上的网格划分不同时，宏观拓扑结构和微观材料分布相似，故本文所提优化设计方法对网格划分没有依赖性，且收敛快速，所得宏观结构和微观材料分布合理，微观单胞多数为kagome或类kagome形.
在单胞的初始密度不同时，如果初始密度相差较大，基本单胞的最优结果相差半个周期，但是由基本单胞组成的微观材料相同;如果初始密度相差较小，基本单胞的最优结果相同，因此微观设计域的初始密度对宏观拓扑结构和微观材料分布影响较小.
在确定性荷载作用下，以柔度为目标函数或者同时以柔度和热变形为目标函数［13］，均表明无法通过增加微观材料用量来提高宏观结构的承载能力，而且宏观结构与微观材料是有耦合效应的［14］.在不确定性荷载作用下，宏观拓扑结构和微观基本单胞的收敛速度有时也不同步，它们之间的耦合效应和材料用量关系还需进一步研究.
［1］ SANCHEZ-PALENCIA E.Non-Homogeneous Media and Vibration Theory［M］.Berlin:Springer-Verlag，1980.
［2］ BENSSOUSAN A，LIONS J L，PAPANICOULAU G.Asymptotic Analysis for Periodic Structures［M］.Amesterdam:North Holland Publishing Company，1978.
［3］ CIORANESCU D，PAULIN J SJ.Homogenization in open sets with holes［J］.Journal of Mathematical Analysis and Applications，1979，71(2):590-607.
［4］ BENDSE M P，KIKUCHI N.Generating optimal topologies in structural design using a homogenization method［J］.Computer Methods in Applied Mechanics and Engineering，1988，71(2):197-224.
［5］ SIGMUND O.Materials with prescribed constitutive parameters:an inverse homogenization problem［J］.Journals＆ Books Create Account Sign in International Journal of Solids and Structures，1994，31(17):2313-2329.
［6］ RODRIGUESH，GUEDESJ M，BENDSOE M P.Hierarchical optimization of material and structure［J］.Structural and Multidisciplinary Optimization，2002，24(1):1-10.
［7］ 阎军.超轻金属结构与材料性能多尺度分析与协同优化设计［D］.博士学位论文.大连:大连理工大学，2007.(YAN Jun.Multiscale analysis and concurrent optimization fo ultra-light metal structures and materials［D］.PhD Thesis.Dalian:Dalian University of Technology，2007.(in Chinese))
［8］ LIU L，YAN J，CHENG G D.Optimum structure with homogeneous optimum truss-like material［J］.Computers and Structures，2008，86(13/14):1417-1425.
［9］ COELHO P G，FERNANDES P R，GUEDES J M，et al.A hierarchical model for concurrent material and topology optimisation of three-dimensional structures［J］.Structural and Multidisciplinary Optimization，2008，35(2):107-115.
［10］ TAKEZAWA A，NII S，KITAMURA M，et al.Topology optimization for worst load conditions based on the eigenvalue analysis of an aggregated linear system［J］.Computer Methods in Applied Mechanics and Engineering，2011，200(25/28):2268-2281.
［11］ GUO X，ZHAO X F，ZHANG W S，et al.Multi-scale robust design and optimization considering load uncertainties［J］.Computer Methods in Applied Mechanics and Engineering，2015，283:994-1009.
［12］ ZHANG W S，ZHONG W L，GUO X.An explicit length scale control approach in SIMP-based topology optimization［J］.Computer Methods in Applied Mechanics and Engineering，2014，282:71-86.
［13］ 阎军，邓佳东，程耿东.基于柔顺性与热变形双目标的多孔材料与结构几何多尺度优化设计［J］.固体力学学报，2011，32(2):119-132.(YAN Jun，DENG Jiadong，CHENG Gengdong.Multi geometrical scale optimization for porous structure and material with multi-objective of structural compliance and thermal deformation［J］.Acta Mechanica Solida Sinica，2011，32(2):119-132.(in Chiniese))
［14］ YAN J，DUAN Z Y，ERIK L，et al.Concurrent multi-scale design optimization of composite frame structures using the Heaviside penalization of discrete material model［J］.Acta Mechanica Sinica，2016，32(3):430-441.