马 晶 赵明宣 王浩淼 刘 湃,2) 亢 战

*(大连理工大学工业装备结构分析国家重点实验室,辽宁大连 116024)

†(上海宇航系统工程研究所,上海 201109)

引言

自20 世纪90 年代Bendsøe 和Kikuchi[1]基于数值均匀化的拓扑优化方法提出以来,结构拓扑优化技术经历了近30 年的快速发展,形成了较为成熟的理论体系,并被广泛应用于航空、航天[2-3]、汽车[4]等各领域结构的轻量化[5]、多功能材料[6]等性能设计.在概念设计阶段,结构拓扑优化技术在满足多种设计约束的前提下,同时优化结构的形状及拓扑以实现特定的优化目标.目前,被广泛采用的拓扑优化方法主要包括变密度法[7](solid isotropic material with penalization,SIMP)、渐进结构优化法[8-9]及水平集方法[10-13],并发展出为各种具体的理论[14]和方法[15].其中,水平集方法基于隐式边界描述,便于处理结构边界的演化、融合,并能够提供清晰的结构边界描述.

在工程领域,包含嵌入式组件的结构较为常见,如卫星支架与导航定位等功能性组件.通常,此类包含嵌入式组件的结构需要应用在空间紧俏的狭小空间(如航天工程[16-17](如图1 所示),微机电系统(micro-electro-mechanical systems,MEMS)以及汽车工程等),因而此类结构也常将部分功能性组件用于结构承载.

图1 典型的飞行器多组件结构系统[17]Fig.1 Typical aircraft multi-component structure system[17]

针对这类结构,功能性组件与支撑材料连接界面的失效,会导致结构的刚度下降甚至完全丧失承载力.因而,考虑连接界面的力学性能对多组件结构进行设计十分必要.此外,为实现多组件结构的最优性能,需同时对支撑材料的分布及功能性组件的形状及布局进行设计,并考虑组件最小距离控制等[16,18]设计约束避免组件间(如电磁等)功能的相互干涉.

目前,针对多材料及多组件结构设计,学者们提出了多种拓扑优化方法:Bendsøe 和Sigmund[19]提出基于变密度法的多相材料插值模型;Wang 等[10]提出color level set 方法,在水平集法框架下采用m个水平集函数描述2m相材料的分布;Wang 等[20]提出多材料水平集方法,利用m个水平集函数描述m+1 相材料的分布,有效避免了设计域内因多个水平集函数相交而产生的冗余区域.张卫红等[16-18]对多组件结构优化设计进行了一系列系统的研究,提出了基于有限圆形近似组件边界位置的多组件距离控制约束及多组件结构拓扑优化方法;此外还有学者对嵌入孔洞的布局优化[21]进行了研究.然而,上述多材料及多组件设计方法,并未考虑材料连接界面可能发生的张开及失效.由于结构拓扑随优化迭代实时演化,如何在优化过程中实时追踪连接界面位置,准确描述界面的张开及失效成为了学者们关注的热点问题之一.Lawry 等[22]针对material anchor 这一特殊双材料结构,考虑材料界面的滑动摩擦接触,对其刚度进行设计,并将此工作扩展到大变形[23]等现象;刘湃等[24-26]基于内聚力模型,结合扩展有限元方法在固定网格下提出考虑多材料及包含固定形状组件结构连接界面力学性能的拓扑优化方法,考虑了界面可能发生的张开与失效,并以刚度为目标对结构进行设计.此外,考虑梯度界面的双材料结构优化[27],考虑脆性材料界面的抗断裂结构拓扑[28]和复合材料[29]结构优化,以及以接触面分布力均匀化[30]为目标的结构拓扑优化等研究都探索了优化迭代时材料界面模型对最终设计的影响.

尽管作者在之前的工作[25]中对考虑连接界面力学性能的多组件结构拓扑优化进行了研究,但所考虑的内嵌组件形状固定不变,这直接导致组件与支撑结构之间连接界面的形状受限于给定组件的几何形状.然而,当考虑连接界面张开时,结构的性能不仅依赖于多组件的布局及支撑结构的拓扑,还显著依赖于连接界面的形状.因而,本文考虑连接界面力学性能,对形状可变的多组件与支撑结构拓扑的协同优化进行研究,采用超椭圆描述可变组件,显著扩展了设计空间,能够描述更为复杂的连接界面形状.本文针对包含多个内嵌组件的结构,考虑组件与支撑材料连接界面的力学性能,同时优化组件的形状、空间布局及支撑结构的拓扑,以实现结构的最优刚度设计(见图2).首先,采用超椭圆模型参数化组件的形状及布局(如空间位置、旋转角度等),并同时构造其水平集函数表达;进而,在固定网格下,以内聚力模型作为连接界面本构,结合水平集描述及扩展有限元方法,准确分析随优化迭代不断演化的连接界面的张开位移;进一步,在水平集方法框架下建立以刚度为目标的拓扑优化列式并推导解析的灵敏度表达,基于数学规划算法[31](method of moving asymptote,MMA)求解优化问题.

图2 包含嵌入组件的结构Fig.2 Structure with two embedded components

1 考虑界面力学性能的多组件结构分析

1.1 基于超椭圆描述的内嵌组件及其水平集函数构造

本文基于超椭圆模型参数化形状可变组件的形状及空间布局.超椭圆定义式为(|x/a|η+|y/b|η)=1,其中a,b分别为超椭圆的长轴、短轴,η 为超椭圆的阶次,且a,b,η ＞0.超椭圆上的点坐标(x,y)可由以t为参数的参数方程表达

超椭圆模型可以描述从菱形到矩形的连续形状变化.如图3 所示,不同超椭圆阶次η 控制了超椭圆的形状:η=2 时为椭圆,η →∞时为矩形.为保证组件边界的光滑性,避免可能出现的应力集中,本文限制参数η ∈[2,6].

图3 参数η 取不同值时超椭圆的形状Fig.3 Superellipse with different parameter η values

超椭圆与坐标轴的交点为(±a,0),(0,±b),其顶点为(±sa,±sb),是从圆点出发到矩形4 个顶点的连线与超椭圆的交点(图3 中的圆点),其中s=2-1/η.超椭圆的面积可以通过Beta 函数计算

基于超椭圆的显式参数表达,可构造其水平集函数为φ(γ)=-(|x′/a|η+|y′/b|η)+1,其中x′,y′为局部坐标,可由组件的中心坐标和角度(如图4 所示) 变换得到

图4 内嵌组件局部和整体坐标Fig.4 Local and global coordinate systems of the components

超椭圆组件的参数γ包括长短轴a,b以及旋转角度θ,以及组件的中心坐标(x0,y0).图5 给出了由参数γ构建的组件水平集函数.

图5 超椭圆组件的水平集函数Fig.5 Level set functions of superellipse components

1.2 结构拓扑及连接界面的水平集描述

本文中组件和支撑结构使用不同材料.本文采用多材料插值模型[20,25]描述设计域中内嵌组件,支撑材料的材料属性,并描述其连接界面位置.对于包含M个内嵌组件的结构,需要M+1 个水平集函数对其进行描述,其中φm(m=1,2,···,M)用于描述组件,φh用于描述支撑结构.利用Heaviside 函数可以表征材料分布如下

其中χm,χh,χv分别表征第m个组件、支撑结构和孔洞的材料分布.故任意点的材料属性可以由如下插值得到

组件q与支撑结构的连接界面可以通过插值表示为

考虑采用距离约束避免组件重叠的情况,式(6)可以简化为

设计域中的全部连接界面可表示为

图6 给出了M=2 时基于水平集的多相材料和连接界面的表征.

图6 基于水平集函数的多相材料及连接界面表征Fig.6 Multiphase material and interface characterization

1.3 连接界面内聚力模型

内聚力模型在断裂力学中被广泛应用于描述多相材料结构的连接界面,表征了界面之间的黏结力[32]和张开位移之间关系的.本文所采用的指数型内聚力模型[33](见图7) 基于光滑的非线性函数,可以描述界面的初始弹性阶段和软化阶段,并且考虑了界面法向张开位移和剪切方向位移的耦合.

图7 指数型内聚力模型示意图Fig.7 Schematic illustration of exponential cohesive zone model

如图8 所示,二维结构的材料连接界面从初始位置受力张开,界面上的P点在张开后移动到点P+和P-.

图8 连接界面张开示意图Fig.8 Schematic illustration of interface’s opening

界面间的黏结力t可表示为

其中,δs和δn分别表示界面剪切方向的位移向量以及法向张开位移向量,tequ和δequ分别表示等效界面黏结力和等效界面张开位移,其表达式如下

其中,β 是调节法向及切向位移耦合关系的权重系数,δne和tne分别表示连接界面上的有效张开位移和载荷向量

其中,n是界面的法向量.

等效界面张开位移和等效界面黏结力中只考虑了界面的法向拉伸以及剪切量的作用,因而排除了界面受压发生破坏的情况.tequ和δequ满足如下关系式

其中,σc是内聚力模型中的界面力的峰值,δc是与之对应的界面张开位移,e是欧拉数.Gf是断裂能,是tequ在δequ从零到无穷的积分,Gf=eσcδc就是图7中的阴影部分的面积.

1.4 扩展有限元方法及非线性问题求解

由于本文采用内聚力模型描述连接界面的力学性能,当计算结构响应时,需要对界面的张开位移(即界面处位移场的强间断)进行准确描述.扩展有限元法[34](extended finite element method,XFEM) 通过在相关节点引入附加的位移自由度并利用间断的扩充形函数进行位移插值,可以在固定网格下描述位移场的间断属性.

图9 中,连接界面穿过结构化有限元网格,在被切割的单元节点上引入的附加自由度,用于插值界面的张开位移.

图9 扩展有限元附加自由度示意图Fig.9 Schematic illustration of XFEM and additional degrees

结构中除界面以外的任意点x的位移可以由有限元的节点位移和附加自由度的位移插值得到

其中,D指设计域,Γcoh表示被界面切割单元的节点集合,Ni和Nj为标准有限元的形函数,ai和bj分别是标准有限元自由度和附加自由度,v(x)为移轴后的间断函数.界面上任意点x的张开位移向量udisc在全局坐标下可以表示为

而法向和切向的界面张开位移向量可表示为

在平衡状态下,考虑材料连接界面张开的多相材料结构的虚功原理表示为

其中,σ和ε为应力和应变向量,f为体力,～t是结构所受的边界力;δu,δudisc和δε分别为虚位移,界面张开虚位移和虚应变.

在数值计算中,材料连接界面为零水平集与离散网格的交点所连成的折线,进而将包含位移间断的有限单元三角化为若干子单元进行高斯积分.因为内聚力模型的非线性,本文采用增量形式的求解格式,将残差表示为R=fint-λfext(其中λ 是载荷因子) (参考文献[21-23])的形式.具体来说,第k+1 个增量步的残差力向量可以在第k增量步上展开为如下表达

其中,KT=∂R/∂d表示切向刚度矩阵,上标k+1 和k表示增量步编号,k+1Δd=k+1d-kd,k+1Δλ=k+1λ-kλ.在每个增量步中,采用Newton-Raphson(N-R)算法迭代搜索结构的平衡状态.

2 拓扑优化问题列式

本节基于速度场水平集方法和参数化的水平集函数建立拓扑优化列示,见方程(18)

对于支撑材料,在水平集网格上定义了一组对应其水平集函数φh的速度设计变量,即1,2,···,n,其中,n为水平集网格节点的数目).支撑材料边界的法向演化速度可以通过设计变量的双线性插值得到:

对于可变组件,超椭圆长短轴和阶数a,b,η,及其位置(x0,y0)和角度θ 为设计变量.每个组件对应一组设计变量γm,可以确定当前组件的水平集.将组件的参数和基体结构零水平集的演化速度作为协同优化问题的设计变量.

结构的外力功可以表示为材料中存储的应变能和材料连接界面力所做功的和,即

约束函数g1约束支撑材料的体积,其中表示给定的支撑材料体积分数值.约束函数g2和g3分别是组件的非重叠约束和最小距离约束.非重叠约束[12]g2将组件约束在设计域之内.g3为距离约束,使任意两个组件之间的距离(即两个超椭圆中心的距离与各自中心到顶点(见图3) 的距离之差) 大于对于不同的组件可以采用不同的值.K为总增量步,和表示支撑材料边界演化速度的下界和上界,其值受限于CFL 条件所规定的最大演化步长.本文使用MMA 优化求解算法[28]更新设计变量(水平集网格点上的速度设计变量)并通过插值获得结构边界的演化速度,然后使用快速行进方法(fast marching)将这些边界速度外推到整个域,进而基于Hamilton-Jacobi 方程

更新支撑材料拓扑.组件的水平集函数在每次参数更新后重新构建.

3 灵敏度分析

本节在速度场水平集法框架下,推导目标及约束函数对于速度设计变量的灵敏度.首先将目标函数表示成各增量步对应值的叠加,即考虑平衡方程,建立优化问题的拉格朗日函数

通过将式(21)对状态变量求导,可得伴随方程.求解伴随方程可得伴随变量(见文献[23]).最终,目标函数的对支撑材料速度设计变量的灵敏度可表示为

其中,参数ξ3=∂D/∂φh为材料矩阵对于支撑材料水平集的导数,且利用其中的δ 函数,将域积分的灵敏度转换成沿零水平集的线积分形式.由于组件由超椭圆参数组γm直接控制,因此目标对设计变量的灵敏度可通过链式法则

得到,其中下标m=1,2,···,M,表示组件的相关参数.此外,由于组件的非重叠约束及距离约束的灵敏度只依赖于组件设计变量,同样可以通过链式法则获得.

4 数值算例

本章在平面应变假设下考虑包含内嵌组件的悬臂梁及MBB 梁拓扑优化问题.支撑材料(材料1)和组件材料(材料2) 的刚度分别为E1=30 000,E2=90 000,泊松比为ν1=ν2=0.15.为避免有限元分析中刚度矩阵的奇异,孔洞相的材料参数设置为E3=10-6,ν3=0.15.内聚力连接界面的参数设置(见方程(12))如下:Gf=0.05,界面力峰值σc=1,对应的界面张开位移峰值δc=Gf/eσc=0.018 4,内聚力模型中界面法向张开位移和切向张开位移耦合的权系数β=2,界面的受压刚度设置为大常数dn=108,以避免界面受压时发生穿透现象.所有算例都通过位移进行加载,在非线性有限元分析中,位移载荷通过20 个增量步逐步增加至dP=0.05.数值算例全部采用边长为10 的正方形双线性单元对设计域进行离散,并采用与有限元网格重合的水平集网格进行拓扑优化.设置演化步长0.5 倍水平集网格尺寸;为了保证水平集函数的符号距离特征,对支撑材料的水平集函数φh每3 个迭代步进行一次初始化;组件的水平集函数每次参数更新后都重新生成.

4.1 悬臂梁结构拓扑优化

图10 给出包含内嵌两个超椭圆组件的悬臂梁结构设计问题示意图.宽800、高400 的悬臂梁结构左边固定,在右边中点施加非零位移.

图10 悬臂梁结构示意图Fig.10 Schematic illustration of cantilever beam structure

在组件形状、布局及支撑材料拓扑协同优化过程中,可变组件和支撑材料的拓扑演化互相影响,不同的初始构型演化出不同的结构.初始时采用了如图11(a)的初始构型对悬臂梁结构进行优化(支撑材料引入稀疏分布的较大孔洞).

图11 拓扑优化不同的初始解Fig.11 Different initial designs

图12 为在此初始结构下协同优化迭代的两个中间结果及第一主应力分布图.如图12(a) 所示,由于孔洞的存在以及支撑材料的体积约束在优化的初始阶段占据主导,组件(图12 中超椭圆曲线表示的部分)被困在局部区域.并且随优化的进行,如图12(b)所示中间设计,结构左侧的组件移动到不受力的孔洞区域以避免组件的连接界面发生拉伸及剪切破坏造成结构的刚度损失.该局部解虽然不会发生界面破坏,但由于功能性组件并未承载,并未对结构的承载能力提供贡献.

图12 悬臂梁优化迭代中间设计的第一主应力分布Fig.12 The first principal stress distribution in intermediate designs of cantilever beam

为了避免这种情况的出现,本文采用了两个阶段的优化策略:首先在设计域布满支撑材料并对组件的布局和形状进行优化,为了避免组件过于靠近导致在孔洞演进时发生重叠,设定距离约束的值C0=60;经过一定步数的迭代,组件会移动到受压区域,此时在支撑材料域内引入加密分布的均匀小孔洞(见图11(b)),进行组件和支撑材料结构的协同优化.

表1 给出了3 种初始构型对应的优化结果,初始的组件均为半径50 处于不同位置的圆,限制a,b∈[40,50],且限定第一阶段的迭代步数为100 步.在3组优化结果中,两组件均演化为a,b=40,η=6 的超椭圆.但是由于组件参数和支撑结构的相互影响,使用上述优化策略得到的最终的拓扑结构和目标值有所不同,但是组件及其界面均处于受压区域(优化结果与图12(b)相似).

表1 支撑材料体积分数约束为0.4 时,不同初始组件位置的悬臂梁结构优化Table 1 Optimization of the cantilever beam structures with different initial components’positions, =0.4

表1 支撑材料体积分数约束为0.4 时,不同初始组件位置的悬臂梁结构优化Table 1 Optimization of the cantilever beam structures with different initial components’positions, =0.4

进一步考察不同支撑材料体积分数对优化结果的影响,使用与表1 中Case 1 相同的初始构型,给定支撑材料体积分数为0.45 和0.5,得到如图13 所示的悬臂梁拓扑结构.从图13 中可见,不同支撑材料体积用量下,组件的最优布局及形状基本不变,且当支撑材料体积为0.5 时其拓扑展现出明显的非对称性.

图13 悬臂梁拓扑优化结果Fig.13 Optimized designs of the cantilever beam

4.2 MBB 梁结构优化

本小节以MBB 梁结构优化为例讨论了组件的形状对优化结果的影响.图14 给出包含内嵌组件的MBB 梁结构设计问题示意图.宽800,高400 的MBB梁结构左边约束水平方向位移,右下角点约束垂直方向位移,并在左下角施加非零位移.

图14 MBB 梁结构示意图Fig.14 Schematic illustration of MBB beam structure

首先,采用与表1 中Case 1 相同的初始构型(水平并列的半径50 的圆) 和优化策略及参数,得到如图15(a)的MBB 梁结构设计,目标值为-0.47,两组件演化为a,b=40,η=6 的超椭圆.经主应力分析,该优化设计中的组件及其界面主要承受压力作用(见图15(b)).

图15 MBB 梁的优化设计Fig.15 Optimized design of the MBB beam

进一步,在优化过程中将组件1 的形状限定为椭圆形,组件2 的形状限定为近似的矩形(即令η1∈[2,2.5],η2∈[4,6]),同样采用与表1 中Case 1 相同的初始构型和参数进行协同优化,得到如图16 所示的迭代历史和优化设计.

图16 MBB 梁结构的优化迭代过程Fig.16 Optimization iteration process of MBB beam structure

两阶段优化策略中,第一阶段依旧限制a,b∈[40,50],但是由于限制了两个超椭圆组件的不同的η范围,在此阶段内组件的形状及布局并未产生显著变化.在第二阶段对组件形状、布局和支撑材料拓扑进行协同优化,并设定组件的长短轴a,b∈[30,80].在协同优化的初始阶段,支撑材料逐渐减少至给定的体积约束值,在此过程中目标值-f随之增大.体积约束满足之后,通过结构拓扑和组件特征的优化,得到了MBB 梁优化设计.同样地,在优化设计中,经主应力分析,组件及其连接界面主要承受压应力作用.两个组件均演化为长轴80,短轴30 的超椭圆,且组件与支撑结构的界面趋于曲率较小的直线或曲线.

虽然本算例及4.1 节算例均以双组件布局与支撑结构拓扑的协同优化作为典型对象进行研究,但本文所采用的有限元分析及所提出的优化方法对包含更多组件的结构优化问题仍然适用.根据对已有优化结果的观察,当组件数目增多时,多组件的布局仍可能分布到结构的受压区域,且优化结果中连接界面趋向于曲率较小的曲线.另一方面,随着组件数目的增加,由于优化过程中组件距离约束等的作用,多个组件之间的相互影响更为复杂,设计问题可能具有更多的局部最优解.

5 结论

材料连接界面的力学性能对多相材料结构的力学性能及其优化设计影响显著.本文基于扩展有限元及内聚力模型考虑材料连接界面的力学性能,使用超椭圆描述了有一定设计自由度的、形状及布局可变的内嵌组件,提出了基于水平集描述的内嵌组件形状、布局及支撑材料拓扑的协同优化方法,本文给出了优化问题的解析灵敏度,并基于数学规划法在速度场水平集方法框架下对优化问题进行求解.通过嵌入组件的悬臂梁和MBB 梁的优化验证了优化策略的可行性.

数值算例显示,考虑连接界面性能的情况下,组件及界面主要处在压应力作用区域,且连接界面最优形状呈现为曲率较小的光滑曲线,该设计避免界面发生破坏进而有效提高了结构的承载力.基于此考虑连接界面的协同优化框架,任意形状的多组件内嵌结构的优化也可以通过引入B 样条的组件描述实现.