考虑载荷不确定性的多材料结构稳健拓扑优化

2019-10-21赵清海张洪信蒋荣超华青松

振动与冲击 2019年19期

赵清海, 张洪信, 蒋荣超, 华青松, 袁林

(1.青岛大学电动汽车智能化动力集成技术国家地方联合工程研究中心，山东青岛 266071;2.青岛大学机电工程学院，山东青岛 266071)

以拓扑优化技术为代表的先进结构设计方法在结构优化领域起到关键的引领作用。相较于尺寸优化与形状优化，拓扑优化应用于概念设计阶段，层次更高且更为复杂，其优化结果为后续详细设计提供参考，对结构的性能以及成本等均起着决定性的影响[1]。

传统的拓扑优化是在确定性条件下进行结构设计，然而由于工艺差异、材料非均质等因素导致材料属性具有不确定性；同时由于制造误差、装配偏差等致使几何尺寸具有不确定性；且由于行驶工况时变性、运行环境多变性等引起载荷工况具有不确定性。在这些不确定性因素影响下，结构性能势必会出现较大波动，甚至发生破坏与失效[2-4]。因此考虑不确定性条件的拓扑优化设计方法具有重要的研究价值。

不确定性拓扑优化设计通常分为两种：可靠性拓扑优化设计[5-6]与稳健拓扑优化设计[7-8]。可靠性拓扑优化侧重安全性，获得满足约束条件失效概率下的最优设计方案；稳健拓扑优化侧重稳定性，降低性能对不确定性因素的敏感度。目前，稳健拓扑优化设计已成为国内外学者研究的热点之一。依据不确定性变量的数学描述，稳健拓扑优化设计又可分为两大类：非概率型与概率型[9-10]。Thore等[11]采用最不利工况法，研究非概率型载荷不确定性对应力约束下结构柔度最小化问题的影响。Zhao等[12]分别利用蒙特卡洛法与矩分解法，进行集中载荷与分布载荷不确定性的稳健拓扑优化设计。Martinez-Frutos等[13]结合GPU并行计算方法与稀疏网格技术进行满足载荷随机场分布特征的稳健拓扑优化研究。付志方等[14]提出了载荷不确定性条件下周期性结构稳健拓扑优化数学模型。Holmberg等[15]等提出广义纳什均衡博弈理论进行椭圆概率分布载荷不确定性的稳健拓扑优化设计。

上述文献研究大都基于单材料进行稳健拓扑优化设计，多材料结构较少涉及。然而，随着3D增材打印技术的出现，多材料结构加工制造成为可能[16-18]。目前多材料拓扑优化设计方法分为以下几类：均匀化/变密度法；相场法；水平集法以及组合优化法等。Sigmund等[19]探讨基于变密度法的三材料拓扑优化数学模型，但构建材料插值模型较为复杂。Wang等[20]建立彩色水平集数学模型，将多材料拓扑优化转为求解Hamilton-Jacobi偏微分方程组，但优化结果往往取决于初始设置。在此基础上，Wang等[21]提出一种多材料水平集拓扑优化模型，能够有效避免“冗余相”的产生，通过提高解的收敛性有效降低初值依赖性。Zhou等[22]提出相场模型，建立基于广义Cahn-Hilliard型偏微分方程的多材料拓扑优化数学模型，计算效率同样受到考验。Tavakoli等[23]结合变密度法与交替相激活算法，将多材料拓扑优化分解为内层两相拓扑优化子问题模型，外层进行各材料设计变量耦合，导致模型计算量较大。

基于此，本文针对多材料结构稳健拓扑优化方法进行研究，构建多材料插值模型，并考虑载荷工况的不确定性。采用有序各向同性微结构材料惩罚模型法(Ordered-Solid Isotropic Microstructures with Penalization, Ordered-SIMP)表征多材料插值模型。考虑载荷不确定性分别为随机变量与随机场分布，针对载荷满足随机场分布时，通过K-L展开将载荷随机场离散化为有限个随机变量加权和，进而采用稀疏网格数值积分方法进行稳健设计矩估计，将多材料稳健拓扑优化转化为求解多工况加权多目标确定性多材料拓扑优化设计问题。

1 多材料结构拓扑优化问题描述

在给定位移、载荷边界条件和体积或质量约束条件下，以获得结构最佳性能(如：结构刚度最大化、特征值最大化或散热弱度最小化)为目标函数，寻求多种材料的用量、空洞构型以及在设计空间的最佳材料分布方案，即为多材料拓扑优化设计[24]。通过对设计区域进行网格划分，多材料拓扑优化转化为获取网格单元填充何种材料或空洞设计问题，多材料设计区域示意图如图1所示。

图1 表征3种材料及空洞的设计区域

2 基于Ordered-SIMP方法的多材料结构插值模型

常用的变密度法材料插值模型有SIMP法与RAMP法[25-26]。SIMP方法引入0-1区间连续变化的单元设计变量，建立材料属性(如：弹性模量、导热系数)与设计变量之间的非线性函数关系，通过控制设计变量取值决定单元取舍，实现材料最佳布局。SIMP法具有形式简单，易于实现敏度推导等优点。对于单材料设计问题，表达式为

E(xe)=(xe)pE0xe∈[0,1]

(1)

式中:xe为单元相对密度，即设计变量；p为惩罚因子；E为插值后材料弹性模量；E0为实体材料的弹性模量。

对于多材料结构设计，首先将各材料的弹性模量、密度进行归一化处理，使材料属性转化为无量纲的相对值，描述为

(2)

引入比例系数AE与平移系数BE，建立基于Ordered-SIMP的多材料插值模型为

(3)

图2 基于Ordered-SIMP的多材料插值模型

Fig.2 Multi-material interpolation model based on Ordered-SIMP

3 多材料结构稳健拓扑优化模型

对于稳健拓扑优化设计，所涉及变量一般分为确定性变量与不确定性变量。其中，确定性变量x，即为多材料结构拓扑优化设计变量；不确定性变量ξ为载荷工况，考虑服从一定的概率分布。基于概率论与数理统计理论，典型的度量结构性能稳健性的指标有：均值与标准差。因此，建立多材料结构稳健拓扑优化数学模型如下

(4)

式中:w为权重因子，w∈[0,1]；K、U与F分别为结构刚度矩阵、位移矢量与结构载荷；ve为单元体积；V0与f分别为设计区域总体积与体积比；Ne为单元总数；xmax与xmin分别为设计变量的上下限；μc(x,ξ)与σc(x,ξ)分别为结构柔度c(x,ξ)的均值和标准差，分别表示如下

(5)

(6)

式中，p(ξ)为ξ联合概率密度函数。结构柔度的均值和标准差相对于设计变量xe的敏度分别表示为

(7)

(8)

4 载荷工况随机场离散化

载荷工况的不确定性可采用随机变量或随机场来表征。对于存在空间关联性的分布载荷满足高斯随机场模型时，可通过K-L变换将随机场转化为有限个不相关的随机变量在相应权重下的累加和[27]。

定义连续空间域Ω的二维载荷随机场为ξ(θ,χ)，其中，θ为空间坐标，定义为θ=(θ1,θ2)；χ标记为随机坐标；则随机场的K-L展开可描述为

(9)

式中，e(θ)为均值；λi与fi(θ)分别为第i阶特征值和正交特征向量，满足Fredholm积分方程

(10)

式中，C(θ1,θ2)为随机场的协方差。ui(χ)为互不相关的随机变量，满足如下条件

(11)

式中，δij为Kronecker-delta函数，满足i=j时δij=1，其它为0。独立正交随机变量ui(χ)定义为

(12)

当K-L展开应用于随机场ξ(θ,χ)离散时，定义d维随机向量ξ(θ)，其元素映射于ξ(θ,χ)中d个观测值。则载荷随机场的K-L展开定义为

(13)

式中，e为d个观测点处随机场均值；λi与fi分别为相关矩阵C的第i阶特征值和正交特征函数，可由下式求解[28]

Cfi=λifi

(14)

式中，相关矩阵C定义为

(15)

在实际问题中，通常用从最大特征值依次降低的前几阶特征值相对应的随机变量来近似反映随机过程的主要概率特征，取前M项，设置M<

(16)

式中，当s足够接近于1时，随机场可有效的通过降维K-L展开来表征[29-30]。

5 稀疏网格数值积分方法

以Smolyak准则为基础的稀疏网格方法，其基本思想是利用一维配置点的特殊张量积操作进行线性组合来构建多维求积公式[31-34]。其优势体现在：配置点数目被限制在一定的范围之内，自动去除对计算精度贡献较小的节点。将稀疏网格应用于稳健设计进行统计矩估计，基于嵌套分层原理，定义一维层间差分格式为

(17)

对含有d维载荷随机变量的性能函数c，构建具有l-水平(l≥1)精度的稀疏网格数值积分格式为

(18)

式中:⊗代表张量积运算符；|k|为多维指标之和(|k|=k1+,…,+kd)。对应的稀疏网格配置点集合定义为

(19)

wi=

(20)

稀疏网格中配置点的数目为

(21)

通过调整水平精度l值，可有效提高稀疏网格积分精度。基于Newton-Cotes积分法则，构建Clenshaw-Curtis型稀疏网格HT∈[-1,1]，定义一维配置点为

(22)

配置点序列为

(23)

相应的权值计算为

(24)

采用稀疏网格数值积分方法进行多材料结构稳健拓扑优化设计的目标函数均值和标准差求解，计算表达式如下

(25)

(26)

目标函数均值与标准差相对于设计变量的灵敏度计算为

(27)

(28)

6 多材料结构设计变量更新

目前结构拓扑优化设计的数值求解算法分为：优化准则法(Optimality Criteria, OC)、数学规划算法(Mathematical Programming, MP)与智能优化算法。其中，OC算法具有收敛速度快，计算规模与设计变量的数目无关等优点[38-39]。本文选取OC算法进行多材料拓扑优化设计变量求解。

OC算法依据库恩-塔克(Kuhn-Tucker，K-T)条件作为优化设计准则[40]，引入拉格朗日乘积因子，建立多材料拓扑优化设计变量迭代格式

(29)

式中:k为迭代次数；η为阻尼系数，取值为0.5；move(k)为正的移动限值，定义为

move(k)=max(α(k)move0,mmin)

(30)

(31)

式中，λ1为拉格朗日乘子，通过半分法计算得到。

在优化迭代过程中，收敛准则根据设计变量相对变化率来判定，定义为

(32)

7 多材料稳健拓扑优化设计算法

综上所述，采用Ordered-SIMP方法进行多材料拓扑优化设计变量表征，同时借助K-L变换将载荷随机场进行离散化，通过稀疏网格数值积分方法进行稳健设计矩估计，进而将多材料稳健拓扑优化转化为求解一组多材料多工况加权多目标确定性拓扑优化设计。在计算过程中，采用OC算法进行设计变量更新，为了获得良好的优化结果，采用传统的灵敏度过滤方法，抑制数值不稳定性问题。多材料结构稳健拓扑优化设计，具体步骤描述如下：

步骤1 多材料设计变量x与相关设置初始化。

步骤2 采用K-L变换将载荷随机场ξ(θ)离散化为有限个随机变量ui(θ)(i=1,2,…,d)。

步骤4 优化循环开始(loop=1)

步骤5 输出多材料结构拓扑构型。

8 数值算例

借助二维数值算例验证所提方法的有效性，算例1与2分别考虑载荷为随机变量与随机场分布时，对多材料拓扑优化结果的影响。

算例1如图3所示L型平板结构，设计区域几何尺寸为L=60，顶端固定，右上端点位置处受垂直向下的集中载荷F作用，分析情况如下：①确定性载荷工况，幅值为1；②不确定性载荷工况，其中幅值与相位为两个独立的随机变量。相位θ满足连续均匀分布，分布区间为[-3π/4,-π/4]；幅值满足正态分布，均值和标准差分别为1和0.30。结构体积分数约束为20%。材料参数设置如表1所示。

图3 L型平板设计区域

表1 三种材料参数设置

设计区域离散为6 400(80×80)个平面四边形单元，对于稳健设计目标柔度的均值与标准差权系数分别设置为0.5。将载荷幅值与相位以及相应的权值通过Clenshaw-Curtis型稀疏网格进行计算，配置点与权值信息如图4所示。设置不同材料组合方案，(1)方案I：材料A、材料B、材料C与孔洞；(2)方案II：材料A、材料C与孔洞；(3)方案III：材料B、材料C与孔洞；(4)方案IV：材料C与孔洞。确定性设计与稳健设计最优拓扑材料分布方案如图5所示。

由计算结果可知，针对不同材料组合方案，所提方法能有效获得多材料结构确定性设计与稳健设计材料分布构型，表明所提方法的有效性。与确定性设计相比较，稳健设计获得不同拓扑构型，特别是在结构左上区域，出现材料填充增添载荷传递路径，这是由于稳健设计所考虑载荷工况的相位发生变化，结构需承受水平方向的载荷分力作用。因此，针对载荷工况不确定性问题，稳健设计所得结果具有更加良好的稳定性。

(a) 配置点分布

(b) 权值分布

图4 稀疏网格配置点与权值信息(d=2；l=4；No=65)

Fig.4 Sparse grid points and weights information (d=2;l=4; No=65)

(a) 确定性设计

(b) 稳健设计

(1) 方案I:材料A、材料B、材料C与孔洞

(a) 确定性设计

(b) 稳健设计

(2) 方案II：材料A、材料C与孔洞

(a) 确定性设计

(b) 稳健设计

(3) 方案III：材料B、材料C与孔洞

(a) 确定性设计

(b) 稳健设计

(4) 方案IV：材料C与孔洞

图5 多材料稳健拓扑优化设计结果

Fig.5 Results of robust topology optimization for multi-materials

对于多材料结构，不同材料组合条件下各材料分布方案具有差异性，这是由于各材料的EN与ρN值不同导致。以Case I为例，在多材料拓扑优化设计过程中，在弹性变形由大到小设计区域，分别由材料A、B、C与Void依次进行填充，同时结构受到材料体积约束的限制。

多材料结构拓扑优化设计目标函数均值与标准差迭代曲线如图6所示(以Case I为例)。结果表明，所提确定性拓扑优化设计与稳健拓扑优化设计方法具有良好的收敛性，且收敛速度与材料的组合数目无关，这是由于算法的收敛设置。同时，确定性设计得到的结构目标柔度的均值与标准差均高于稳健设计所获得结构柔度的均值与标准差，所提多材料结构稳健设计方法具有良好的鲁棒性。充分证明了所提方法的有效性。

图6 柔度均值与标准差随迭代次数的变化曲线(Case I)

Fig.6 The iterative history of the mean and standard deviation of compliance(Case I)

算例2如图7所示简支梁结构，设计空间为120×40的平面四边形区域，底部两端进行位移约束，顶端施加均布载荷，针对以下两种载荷工况进行探讨：①载荷工况满足确定性条件，幅值为1；②载荷工况满足不确定性条件，服从随机场分布，均值为1，标准差为0.3。随机场的协方差满足关系式

(32)

图7 简支梁设计空间

将设计空间离散为4 800(120×40)个平面四边形单元，材料的体积约束设定为0.3。设置三种不同的三材料组合方式，密度与弹性模量的归一化关系曲线如图8所示。

对于稳健设计目标函数均值与标准差的权重系数分别设置为0.5，对于载荷工况随机场，采用K-L展开离散化为三个随机变量，表示为

(33)

式中，μF为载荷均值。借助稀疏网格数值积分方法，获得三个随机变量离散配置点空间分布，如图9所示。确定性拓扑优化设计与稳健拓扑优化设计结果如图10所示。

由于结构的设计区域、垂向位移约束以及载荷工况均具有左右对称特征，因此确定性设计与稳健设计拓扑构型具有对称性。但随机场载荷条件下的材料分布较确定性载荷条件发生明显的变化，主要集中在结构中间下部填充杆状材料，这是由于稳健设计时不同幅值情况下多载荷工况影响相互叠加造成的。当载荷发生变动时，结构整体承载更加稳健，表明稳健设计的必要性。同时可得到所提出的多材料拓扑优化设计方法能有效获得多材料分布构型。针对不同的材料组合方案，该模型均具有良好的适用性。