载荷作用位置不确定条件下结构动态稳健性拓扑优化设计1)

王栋

(西北工业大学航空学院,西安 710072)

引言

拓扑优化设计的目的是在规定的区域内和边界条件下,以最优的方式配置可用的材料,使结构的力学性能及响应指标在给定静、动力载荷作用下达到最优[1-4].当结构的初始构型设计无法借助于设计者以往的经验可以确定时,拓扑优化设计能够在给定的载荷和一定的应力、位移等约束条件下搜寻出创新性的结构构型(即传力路径) 设计,充分发挥材料的性能潜力,从而显著提高工程设计的效率和质量[5-6].传统的结构优化设计问题中,结构所受的载荷环境、材料参数以及优化的目标函数、约束条件等均被认为是确定性的,即所获得的最优结构设计都是基于相对理想状态下的结果[6-10].然而在实际工程结构中,不确定性因素却是普遍存在和不可避免的,如材料的物理性能、模型几何尺寸、边界约束条件、外界施加的载荷等[11-17].这些系统参数即使是在小范围内随机性地变化也无疑会对结构的构型设计产生巨大的作用,最终将对结构的力学性能和响应造成一些不利的影响.若忽略设计参数的不确定性,依旧按照确定性情形对结构进行优化设计,所获得的优化结果抵抗结构参数(物理、几何)或载荷等不确定性因素的能力将会极大地减弱,甚至可能导致结构功能的提前失效[17].为了使所设计的结构具有较强的抵抗不确定性因素的能力,使其能正常、稳定地发挥预期的功能,设计者需要在结构的初步设计阶段就主动地考虑这些不确定性因素的影响[16-18].

航空航天、汽车、桥梁等结构在其服役期内经常受到各种时变载荷的作用,而一个时变的激励通常能近似分解为几个简谐激励.因此简谐激励是一种非常基础和典型的动载荷[7-8],研究简谐激励作用下结构的拓扑优化设计具有十分重要的理论和实际意义[19].然而,当前大多数关于结构动态拓扑优化的研究都是基于确定性参量来构建优化问题的数学模型,很少考虑工程实际中的不确定性因素,从而明显降低了优化设计结果的可靠性与可信度,难以满足现代设计的需求[5-9].实际工程中,由于外部环境的变化,使得作用于结构上的外载荷,包括其大小、方向与位置等都可能具有一定的随机不确定性[20-23].而且相较于材料属性、边界条件等诸多不确定性因素,载荷不确定性对结构拓扑构型的影响更大.因此研究载荷不确定性对连续体结构拓扑优化设计的影响具有非常现实的意义.

研究不确定性条件下结构拓扑优化设计可采用基于可靠性的拓扑优化方法(reliability-based topology optimization,RBTO)[11-12,19],或基于稳健性的拓扑优化方法(robust topology optimization,RTO)[14,20-22].若一个不确定性设计量的随机变化性无法或难以用概率模型表示,如在结构的初始设计阶段,准确地获得不确定性变量的概率分布函数或特征参数都是一件十分困难的事,但却可以简单地用一个有界的区间变量(bounded interval variable)来表示其随机性的变化.此时采用旨在减小结构响应对不确定性设计量敏感性的稳健性设计思想,完成对结构拓扑构型的优化设计,并使其具有较高的抵抗不确定性变量扰动的能力.在结构的初步设计阶段,通常很难准确获知外激励变化的分布模式,但其变化的界限却较容易确定.因此,载荷的不确定性更适于采用“未知但有界(uncertain-but-bounded)”区间变量来描述[12-14,17-18].传统考虑区间不确定性的拓扑优化设计采取的是嵌套双循环模式,即内循环采用区间分析的方法寻找不确定性变量的最大(差) 影响状态,外循环采用确定性的优化策略实现结构拓扑设计的更新与演化.本文采用灵敏度分析的方法研究外载荷作用位置的不确定性,按照确定性的优化策略,通过单循环模式完成结构在动载荷作用下的稳健性拓扑设计.保证在外激励作用位置出现扰动时,结构的动响应不发生较大的改变.

本文基于材料属性的有理近似方法(rational approximation of material properties,RAMP)建立连续体结构简谐响应拓扑优化模型,可有效地避免虚假局部模态的产生[21,24].以结构动柔顺度最小化为目标[8-9],以整体材料用量为约束条件,采用移动渐近法(method of moving asymptotes,MMA)[25],以载荷扰动区间内最大绝对值灵敏度作为拓扑变量修改依据,研究在简谐激励作用下结构强迫振动稳定阶段响应的拓扑优化问题.在优化求解过程中采用启发式敏度过滤(sensitivity filter) 技术[1],保证优化求解过程的稳定性.最后用两个典型优化算例来证明本文所提方法的有效性,并与传统确定性条件下的优化结果进行对比来验证稳健性拓扑优化设计结果的可靠性.

1 稳健性拓扑优化模型

1.1 线性振动系统稳态响应幅值

当外激励频率远离结构的固有频率时,结构内的阻尼影响很小,可以暂时忽略不计.采用有限元法对连续体结构进行离散化处理.于是无阻尼结构受迫振动微分方程可表示为

式中,K和M分别是振动系统的刚度和质量矩阵,它们都是N×N阶的实对称矩阵,N表示结构的自由度数.u(t)和分别是结构的位移和加速度响应N维列阵.f(t)是系统所受外激励列阵.假设结构受到简谐外力的作用

其中,F(s)是外激励幅值N维列阵,s表示外激励的作用位置.由于s是一个不确定的变量,可在一定的区域内随机变化,因此外激励幅值列阵F(s) 也具有一定的不确定性,优化过程中将不断发生改变.ω 是外激励频率(弧度/秒),是一个给定的量.

在简谐载荷作用下,线性系统的稳态响应也是同一频率的简谐运动

式中,U是结构位移响应幅值列阵.由于不考虑结构内的阻尼,U是实数,其中各项的正、负值与结构的动态性能以及激振频率密切相关.将式(3)代入方程(1),消去时间项,则在频域内振动系统响应特征方程可表示为

其中,D(ω)=K-ω2M是振动系统的动刚度矩阵,是激振频率的函数.若激励频率不等于系统的固有频率,响应特征方程(4) 的系数矩阵不会奇异.本文采用直接方法精确计算结构在简谐激励作用下的稳态响应,不必考虑模态截尾带来的计算误差和虚假模态的干扰[7,24].求解方程(4)可得结构响应的幅值列阵U,进而可以得到时域内结构的动响应u(t).

1.2 结构拓扑优化数学模型

实际工程结构中,由于外载荷作用位置的不确定性,结构的响应幅值也具有一定的不确定性.本文以结构有限元模型中单元的相对密度ρe作为设计变量,在给定设计区域的边界与加载条件下,拓扑优化的目标函数设置为结构的动柔顺度最小[7-8].根据动态位移响应的特点,稳健性拓扑优化问题数学列式可表示为

式中,Cd是结构的动柔顺度,是一个正定的值.动柔顺度不仅与外激励的大小有关,而且也与载荷作用位置密切相关,即随着外激励位置s的变化而改变.sI是一个有界的区间变量,和分别代表区间的下限和上限[26-28].ρ是拓扑变量n维列阵,n代表结构的单元总数.ρmin是为了避免结构刚度矩阵奇异而设置的一个极小值,一般取ρmin=0.001.ve是单元的体积,fv是给定的材料体分比,即材料用量的上限.

以上模型与传统的拓扑优化模型最大区别在于:结构的动柔顺度Cd不再是一个确切的目标函数,而是变成一个与区间变量sI有关的区间不确定值,因而也具有一定的随机性.从而使得优化设计过程更加复杂,但也真实地反映了工程实际情况.

为了合理地描述结构的动力学性能,根据以往的研究结果,本文采用RAMP 模型对单元质量和刚度矩阵进行插值[21]

式中,q是插值模型的惩罚因子(q=5);m0和k0分别是实体材料单元的质量和刚度矩阵.

1.3 结构动柔顺度一阶灵敏度的计算

当外激励作用位置s在有界区间内随机变化时,根据文献[29]的定义,只需在外激励的名义(理想)作用点,如区间的中点

处,对振动系统进行分析,计算结构的动响应.并在单元层面上计算结构的动柔顺度Cd0

因此,其对单元的相对密度 ρe的一阶导数灵敏度可由下式计算得到

其中,ue是单元的动响应幅值列阵.de=ke-ω2me是单元的动刚度矩阵.由RAMP 模型材料插值函数可得

于是,通过一次有限元分析计算,即可获得振动结构的名义动柔顺度Cd0及其对拓扑设计变量的一阶导数灵敏度值.但当外激励的作用点发生偏移时,结构的动柔顺度将不再依靠有限元重复计算而得到.

2 外载荷位置不确定性的影响

2.1 载荷位置移动对结构动柔顺度的影响

在复杂的载荷环境下,虽然已经假设外激励的大小和方向保持不变,但当载荷作用位置发生改变时,结构的传力路径亦将发生改变.即外激励作用位置的变化必然引起结构动柔顺度的变化.若能够获得结构动柔顺度对载荷位置的一阶和二阶导数,则可以利用二阶泰勒级数展开公式,得到外激励作用位置发生改变以后,结构动柔顺度Cd(s)的显性表达式.即可避免有限元计算的重复运行,减少振动系统动响应方程的求解.

在载荷名义作用点周围一个较小的邻域内,利用二阶泰勒展开公式,得到载荷位置发生改变以后,结构动柔顺度Cd的近似表达式

式中,Δs=s-s0表示载荷作用位置的改变量.∂Cd0/∂s和∂2Cd0/∂s2分别是在名义作用点s0,动柔顺度对外载荷位置的一阶和二阶导数.

对于常用的四结点双线性平面应力单元,文献[30] 已经推导出了外载荷对其位置(坐标) 的一阶导数计算公式.于是在名义作用点,结构的动位移响应幅值列阵U对外激励位置的一阶导数,可由方程(4)两边对s求一阶微分得到

其中,∂F(s0)/∂s表示外激励对其作用位置的一阶导数.

将结构的名义动柔顺度改写成

则动柔顺度对载荷位置的一阶导数可表示为

据此,可进一步计算结构的动柔顺度Cd0对载荷作用位置的二阶导数

将式(10)、式(16)和式(17)代入式(13),即可利用该显性表达式快速计算由于外激励位置的移动而引起结构动柔顺度的改变值,以及对拓扑设计变量,即单元相对密度ρe的一阶导数

其中,式(18)右边第一项由式(11)计算得到,第二、三项可由式(11)分别对外载荷作用位置s连续微分得到

根据以上推导结果,即可利用式(18) 计算当外激励作用位置发生改变以后,结构的动柔顺度对所有设计变量的一阶灵敏度,这使得动态优化设计效率大为提高.

2.2 载荷位置移动对优化方法的影响

由于外激励作用点的不确定性,振动结构的动柔顺度及其一阶导数也都存在一定的不确定性,即动柔顺度对设计变量的灵敏度随载荷作用点的移动而改变.但是根据MMA 优化算法[25],每次迭代只能选取一个动柔顺度的灵敏度值进行优化设计,因此在运用拓扑优化算法时会遇到一定的困难.为了提高结构对外载荷位置扰动的抵抗能力,文献[29] 提出在载荷位置不确定区间内,按照绝对值最大的原则选取目标函数的一阶导数灵敏度值.从结构动力学分析角度来看,动柔顺度对设计变量一阶导数的绝对值越大,说明单元密度的改变引起动柔顺度值的变化也越大.以此为基础开展结构拓扑优化设计,可以充分考虑载荷作用位置随机变化对结构动柔顺度(或动刚度)的影响.在保证外激励作用区间完整的基础上,本文采用外激励F(s)是在3 个特殊作用点,即名义作用点s0和不确定区间两个端点和的灵敏度值,近似代替其在整个区间内的变化情况,即用这3 个作用位置覆盖整个外载荷作用点的不确定区域.每次迭代时选取三者中绝对值最大的灵敏度对设计变量进行更新,实现结构拓扑的演化

此外,还要利用灵敏度过滤技术,对动柔顺度灵敏度进行再分配,以便消除优化过程中的棋盘格现象,确保优化结果的可行和可靠性[28].

3 算例分析

3.1 简支平板结构

简支平板结构设计区域几何尺寸及边界约束如图1 所示,厚度为10 mm.理论上,结构在设计区域底边中点处承受一个集中简谐力f(t)=10ejωtkN 的作用,外激振频率150 Hz(ω=2π×150 rad/s).但是由于加载环境的复杂性,外激振实际作用在距离中点±10 mm 点某个不确定点上.将设计区域均匀划分成90×40 的有限元网格,采用四结点平面应力单元.假设材料的弹性模量E=200 GPa,密度ρ=7800 kg/m3,泊松比ν=0.3.材料体积约束系数fv=0.3.取过滤半径为单元尺寸的1.5 倍,优化过程收敛条件设置为相邻两次迭代各设计变量的最大改变量小于0.01[25].初始设计时假设材料均匀分布在设计域内,此时结构的第一阶固有频率为262.63 Hz.

图1 简支平面矩形板结构Fig.1 A simply supported rectangular panel

图2 分别示出了确定性拓扑优化设计(deterministic topology optimization,DTO)和稳健性设计所获得的结构最优构型.可见,虽然两种优化策略得到的结构构型基本类似,但细节设计还是有较大差异的.与确定性优化设计相比,稳健性拓扑优化设计主要有以下几个特点.

图2 两种设计策略结构拓扑优化结果对比Fig.2 Comparison of the topology optimizations on two different design strategies

(1)有更多的材料聚集在外载荷名义作用点附近(12:8 个单元宽度),这使得外激励作用区域得到显著加强,以便能更有效地抵抗外激励位置的扰动.毋庸置疑,这是考虑外载荷作用位置不确定性拓扑优化设计的必然结果[22,27,29].

(2) 内侧(中间) 两根斜撑杆(或构件) 的上端出现了分岔,导致与外缘构件连接处呈现圆滑过渡.内侧斜撑杆是竖向外力从作用点向约束支承点(基础)传递的主要路径.当外激励的作用点固定不动时,外激励将沿着斜撑杆先传到结构的外缘,然后再传递到结构的支承点,整个传力路径都是对称的.但是当外激励的作用点发生不确定性变化时,外激励沿内侧斜撑杆以及整个结构的传递都不再是均匀对称的了.在斜撑杆的上部将沿着不同的路径传递到外缘构件,以均衡外载荷的传递,抵消外激励作用位置的变动对结构动柔顺度(动刚度)影响.

(3) 内侧两根斜撑杆底部的截面尺寸明显增大,外侧两根斜撑杆的截面尺寸略有减小,而上、下水平构件的截面尺寸几乎未变.内侧两根斜撑杆与载荷作用点直接相连,主要抵御在竖直方向上的外力.外激励作用位置的变化会直接导致内侧斜撑杆截面上的内力分布不均匀,而较大截面尺寸的斜撑杆有助于外载荷更加有效地传递出去.因为可用的材料是一定的,因此这一部分增加的材料只能从外侧斜撑杆调配得到,这就导致了外侧两根斜撑杆截面尺寸的减小.

表1 列出了两种设计策略结构拓扑优化数值结果.首先,考虑外激励作用位置不确定性的稳健性拓扑设计,名义最小动柔顺度Cd0比确定性优化相应结果稍大一些(1.64%).这是因为有一部分材料用来加强载荷作用区域,以抵抗外激励作用位置的扰动;还要一部分材料用来在斜撑杆的上部构造额外的传力路径,增加载荷传递的可靠性.因此在相同材料体分比条件下,为了提高结构抵抗外载荷扰动的能力,结构的动刚度会有一定的损失,但这部分损失可以通过适当地增加结构材料得到有效补偿[27].

实际上,在相同的材料体积约束下,动柔顺度的名义最小值,与其对不确定性设计参量的敏感性分别是一对矛盾的两个方面.因此,在它们之间存在有一定的相互制约性,过度地强调一方必然会导致另一方性能的弱化或缺失.而优化设计过程就是要在相互妥协中找到矛盾双方最佳的解决方案.本例中,只要将优化结果的结构材料体分比增大到fv=0.306,就可以将稳健性拓扑优化设计的名义最小动柔顺度降到与确定性优化设计相同的结果,见表1 的最后一行所示.

其次,考虑外激励不确定性的稳健性拓扑设计收敛过程明显增长.由于每次优化迭代选取的动柔顺度对拓扑设计变量一阶导数都是绝对值最大的,单元之间灵敏度值分散性增大,利用基于梯度的优化算法获得稳定的结构拓扑优化结果更加困难.因此,稳健性优化迭代过程会显著增长.

从表1 的后二列可知,当外载荷作用位置发生改变以后,结构的动柔顺度较之其名义值均有所增大,即结构的动态刚度均有所降低.当外载荷移到其不确定位置区间的左(或右) 端点时,此时作用点的偏移量相对于板的宽度最大(1.11%).对于稳健性拓扑优化设计,结构的动柔顺度变化率远低于确定性优化的结果,特别是在保持结构的动柔顺度名义最小值相同的设计情形下.此外,如果由于某种原因,外激励作用点发生更大的偏移,确定性拓扑优化设计的结构动柔顺度将快速增大,远大于稳健性优化的结果,如图3 所示.这表明采用本文提出的稳健性拓扑优化设计后,结构抵抗外载荷作用位置扰动的能力有了显著地增强.即在不(或者少)增大结构材料用量的基础上,结构动态响应的稳健性明显提高.由此可见,在结构拓扑优化设计中考虑载荷不确定性是非常有必要的.

表1 结构拓扑优化结果以及动柔顺度随载荷位置变化情况Table 1 Optimized results and the dynamic compliance variations

图3 结构动柔顺度随外激励作用位置变化情况Fig.3 Variation of the structural dynamic compliance caused by the external loading point disturbances

比较结构的前三阶固有频率可以发现,结构的低阶固有频率,特别是其第一阶固有频率,经过拓扑优化后都得到明显提高,如表2 所示.这使得激励频率更加远离结构的共振区,结构的阻尼作用基本可以忽略.在当前激励频率下,结构的总体动刚度较初始设计有了显著增大,从而使结构的动响应幅值急剧减小.此外,在相同材料体分比约束条件下,稳健性拓扑优化设计的结构固有频率要略小于确定性优化的结果,这也是为了提高结构响应的稳健性,有效抵抗外激励作用位置扰动的必然结果.

表2 简支平板结构前三阶固有频率比较Table 2 Comparison of the first three natural frequencies of the rectangular panel structure

3.2 MBB 梁结构

图4 所示为MBB 梁结构设计区域及尺寸,左、右两端完全固定,假设其厚度为10 mm.理想状况下,一个集中简谐力f(t)=10ejωtkN 施加在梁上边缘的中点处,外激振频率是200 Hz.但实际上外载荷作用在中点周围±20 mm 区间内某个不确定的点上.将结构设计区域均匀划分成150×25 的有限元网格.体积约束系数fv=0.35,材料性能参数与上例相同.初始设计时,将所给的材料均匀分布在整个设计区域内.

图4 MBB 梁设计区域及外激励Fig.4 Design domain of MBB beam and the external force

图5 分别示出了确定性和稳健性优化设计得到的最优结构构型设计.可以看到,两种设计策略得到的结构拓扑优化结果依然存在较大的差别.与前例相同,稳健性拓扑优化设计驱使更多的材料聚集在外载荷作用点附近(14:6 个单元宽度).当外激励偶然偏离名义作用点时,较宽的受载区域能有效应对外激励位置移动的影响,以便更可靠地将外载荷传递到基础.在确定性优化设计时,中间两根斜杆张成90°的夹角.表明外激励先沿着相互正交的两个方向传递到下弦,然后再传递到固定基础上.而在稳健性优化设计时,中间两根斜杆张得角度更大,夹角超过了100°.此外,确定性优化设计中间两根斜杆的截面尺寸较小,外侧斜杆的截面尺寸较大.而稳健性优化设计中间斜杆的截面尺寸,特别是上部靠近外激励作用点处较大,外侧斜杆以及下弦中间的截面尺寸则变小,以满足对材料体积的约束.

图5 两种策略结构拓扑优化设计结果对比Fig.5 Comparison of the topology optimizations on two different design strategies

图5 两种策略结构拓扑优化设计结果对比(续)Fig.5 Comparison of the topology optimizations on two different design strategies(continued)

表3 给出了两种设计策略的优化数值结果.与前例一样,此时的稳健性拓扑设计动柔顺度名义最小值Cd0仍比确定性优化结果稍大一些(6.07%),即结构的动刚度有一些损失.但这同样可以通过适当地增加结构材料(fv=0.375)得到有效补偿.表3 的后四列给出了拓扑优化设计后,结构动柔顺度随外激励作用点移动的变化情况.可以看到,虽然外激励作用位置相对于梁结构的跨度最大扰动只有1.33%,但是确定性拓扑优化设计结构动柔顺度的变化率相当高,远大于稳健性拓扑优化设计相应的变化率.因此可知,运用本文所提方法获得的优化构型设计,极大地提高了结构的动响应对外激励作用位置扰动的抵抗能力.

表3 MBB 梁结构拓扑优化结果以及动柔顺度变化情况Table 3 Optimized results and the dynamic compliance variations

表4 给出了优化前后结构的前三阶固有频率.初始设计时,外激振频率与结构的第一阶固有频率比较接近.经过拓扑优化后,结构的第一阶固有频率得到显著地提高,使得激励频率远离结构的共振区.图6分别绘出了两种载荷作用状况下,结构的名义动柔顺度Cd0迭代收敛曲线,同时还示出了每隔30 迭代步的中间构型和最终优化设计结果.值得注意的是,在结构稳健性优化设计过程中,中间构型设计始终都有较多的材料聚积在外激励作用区间周围,以增强结构抵抗外激励作用点的随机扰动.

表4 MBB 梁结构前三阶固有频率比较Table 4 Comparison of the first three natural frequencies of the MBB beam structure

图6 两种设计策略拓扑优化收敛过程Fig.6 Convergence curves of the structural dynamic compliance on the two design strategies

4 结论

本文研究了在外激励作用位置不确定条件下,连续体结构动态稳健性拓扑优化设计问题.由于在结构的初步设计阶段,很难获知外载荷变化的分布规律,本文将外激励作用位置的不确定性用有界区间变量表示.通过灵敏度分析,借助二阶泰勒展开公式,将结构动柔顺度及其对拓扑变量的一阶灵敏度,均表示成激励作用位置变化量的显函数.运用有理近似材料性能模型和移动渐近算法,以绝对值最大的动柔顺度灵敏度作为变量修改设计的依据,充分考虑载荷位置的不确定性对结构构型设计的影响.本文用两个典型结构动态稳健性拓扑优化设计算例,验证了所提优化策略的可行性,分析了考虑外激励作用位置的不确定性对结构构型优化设计的影响,并比较了结构的动柔顺度对外激励位置扰动的抵抗能力.研究结果如下:

(1)考虑载荷作用位置不确定性条件下得到的稳健性结构拓扑优化设计,不仅载荷作用点周围的材料会增多,而且拓扑构型设计也与确定性条件下的拓扑优化结果明显不同,使所得结构在载荷位置随机扰动情况下具有更可靠的传力路径.

(2)当载荷作用位置发生扰动时,稳健性拓扑优化设计的结构动柔顺度变化率远低于相应的确定性优化结果.充分表明所得到的结构拓扑构型抵抗外激励位置扰动的能力更强,结构性能更加稳健.

(3)在相同体积约束情况下,稳健性拓扑优化设计的名义最小动柔顺度稍大于相应的确定性拓扑优化结果.这是因为需要用一部分材料加强外激励作用区域,并且还要构造额外的传力路径.因而结构稳健性拓扑优化设计会损失一部分动刚度,以换取其抵抗载荷位置扰动能力的提高.但这部分损失的结构动刚度,可以通过增加少许结构材料而得到有效弥补,并且其对外载荷位置扰动的敏感性会进一步降低.