引 言

结构拓扑优化通常是在给定载荷和约束等条件下寻找材料的最佳分布,其在结构设计的初始阶段能高效的给出最优构型^[1].结构拓扑优化与尺寸和形状优化相比具有更大的效益,但同时也最具挑战性^[2].自从文献^[3]提出均匀化方法,结构拓扑优化在理论和应用上都取得了长足的发展^{[4, 5]}.目前比较常用的连续体结构拓扑优化方法有:均匀化法^[3],固体各向同性材料惩罚法^{[6, 7]},渐进结构优化法^{[8, 9, 10]},以及最近发展起来的水平集法^{[11, 12, 13]}.

通常情况下,结构拓扑优化是在确定条件下完成的,但是这种最优结构在抵抗外界扰动时有可能是比较脆弱的^[14],因此研究不确定条件下的结构拓扑优化问题是很有必要的.在实际的设计过程中存在各类不确定性,如载荷不确定^{[15, 16, 17]},材料不确定^{[14, 18, 19]},制造过程产生的几何尺寸和边界不确定^{[20, 21, 22]}.本文考虑载荷不确定性进行结构拓扑优化设计.

目前处理不确定性的结构拓扑优化方法主要包括两种:基于可靠度的拓扑优化和稳健性拓扑优化^{[23, 24, 25]}.基于可靠度的拓扑优化主要强调结构拓扑形式的高可靠性,该方法是在满足指定可靠度约束的基础上来最小化结构质量;稳健性拓扑优化主要考虑结构抵抗干扰的能力,该方法的目标函数通常是结构柔度的均值和方差^[26].由于其重要性,近几年许多学者在稳健性拓扑优化上进行了大量的研究工作.有人提出了以均值和标准差最小的多目标结构稳健性优化方法^[27],也有人用非概率的方法研究了几何非线性的结构拓扑优化问题^[28],文献^[19]基于水平集研究了在载荷和材料性能都是随机不确定性条件下的稳健拓扑优化方法,还有人将随机规划法和水平集法相结合分析了基于随机不确定载荷下的结构拓扑优化^[29],文献^{[25, 30, 31]}用线弹性叠加原理给出了单工况下结构柔度的均值和方差的计算方法,并以均值和标准差的加权和为目标进行结构拓扑优化.

在上述的研究中,大部分是在单一工况下^{[27, 28, 29, 30, 31]}进行的稳健性拓扑优化设计的,多工况条件下的研究较少^{[32, 33]}.然而,在实际工程应用中,工况往往比较复杂^{[34, 35]},因此考虑多工况下结构的拓扑优化更符合实际情况.基于此,本文将对多工况下的连续体结构稳健拓扑优化方法进行研究,其中不但考虑载荷大小的不确定性,还考虑了其方向的不确定性.

文章首先给出了多工况、不确定载荷下的结构稳健拓扑优化模型;然后,利用线弹性结构的位移叠加原理给出多工况下结构柔度的均值和方差的计算公式,并在此基础上给出了结构的灵敏度分析方法;最后,通过算例说明本方法的有效性以及优化结果的稳健性.本文主要讨论二维问题,但本文方法能很方便地推广到三维情形.

1 多工况线性结构稳健拓扑优化设计问题 1.1 优化模型

本文采用密度法进行拓扑优化问题的求解，设计变量为单元密度$\rho _e$，并采用约束最小刚度的插值模型将单元材料属性$E_e$ 与单元密度$\rho _e$联系起来^[36]

${{E}_{e}}({{\rho }_{e}})={{E}_{\min }}+\rho _{e}^{p}({{E}_{0}}-{{E}_{\min }}),{{\rho }_{e}}\in [0, 1]$

(1)

其中E₀为材料的弹性模量，E_min 是一个为了防止刚度矩阵奇异而赋予无材料单元的很小的刚度值. $p\ge 1$为惩罚因子，为了迫使单元密度在优化过程中取值趋向0或1，从而获得构型清晰的结构.

假设结构被离散为N个单元，${f}^\gamma (\gamma = 1,2,\cdots,M)$为M个工况，则多工况、载荷不确定条件下 的连续体结构稳健拓扑优化问题可以表示为

$\left. \begin{array}{*{35}{l}} \begin{align} & Min:J=\alpha \mu (c)+\beta \sigma (c) \\ & s.t.:K{{u}^{\gamma }}(\omega )={{f}^{\gamma }}(\omega )(\omega \in \Theta )(\gamma =1,2,\cdots ,M) \\ & g=V(x)-{{V}_{\max }}=\sum\limits_{e=1}^{N}{{{v}_{e}}{{\rho }_{e}}-{{V}_{\max }}\le 0} \\ & 0\le \rho \le 1 \\ \end{align} & & & & \\ \end{array} \right\}$

(2)

其中$\mu (c)$与$\sigma (c)$分别为多工况、载荷不确定条件下结构柔度的整体均值和标准差，非负常数$\alpha$和$\beta$ 用 于权衡均值和标准差的重要性，且$\alpha + \beta = 1.K{u^\gamma }(\omega ) = {f^\gamma }(\omega )$为结构在第$\gamma $个工况下的有限元平衡方程，$\omega \in \Theta$表示载荷具有不确定性，$v_e$ 为单元e的体积，$\rho _e$ 为其密度. V_max 为许用材料体积，${\rho}$为由各单元密度组成的设计变量向量.

假设各工况之间相互独立，则整体柔度均值和方差用各工况的柔度均值和方差分别表示如下

$\hskip 15mm \mu (c) = f_\gamma = 1^{M} \varpi ^\gamma \mu (c^\gamma )$

(3)

$\hskip 15mm \sigma ^2(c) = f_\gamma = 1^{M} (\varpi^\gamma )^2\sigma ^2(c^\gamma )$

(4)

其中$\varpi ^\gamma $表示第$\gamma $个工况所占的比重，权衡该工况在整体中的重要性.

将式(3)和式(4)代入式(2)得到多工况、载荷不确定条件下的连续体结构稳健拓扑优化问题的最终表达形式

$\left. \begin{align} & Min:J=\alpha \sum\limits_{\gamma =1}^{M}{{{\varpi }^{\gamma }}\mu ({{c}^{\gamma }})}+\beta \sqrt{\sum\limits_{\gamma =1}^{M}{{{({{\varpi }^{\gamma }})}^{2}}{{\sigma }^{2}}({{c}^{\gamma }})}} \\ & s.t.:K{{u}^{\gamma }}(\omega )={{f}^{\gamma }}(\omega )(\omega \in \Theta )(\gamma =1,2,\cdots ,M) \\ & g=V(x)-{{V}_{\max }}=\sum\limits_{\gamma =1}^{M}{{{v}_{e}}{{\rho }_{e}}-{{V}_{\max }}}\le 0 \\ & g=V(x)-{{V}_{\max }}=\sum\limits_{\gamma =1}^{M}{{{v}_{e}}{{\rho }_{e}}-{{V}_{\max }}}\le 0 \\ \end{align} \right\}$

(5)

1.2 单工况下结构柔度均值和方差的计算^{[25, 30, 31]}

假设在某一工况$\gamma$下,结构承受$n^{\gamma}$个不确定载荷,并且每个载荷的大小$h_i^\gamma$和作用方向与x轴的夹角$\theta_i^\gamma$的概率分布均已知.设${f}_2i-1^\gamma和{\pmb f}_2i^\gamma$分别为单独在第i个不确定载荷作用点沿x方向和沿y方向施加单位载荷时的结构整体载荷向量,并且在${f}_2i-1^\gamma$和${f}_2i^\gamma$作用下结构的位移向量分别为${u}_2i-1^\gamma和{\pmb u}_2i^\gamma$,则根据位移叠加原理,对于任意的载荷

${F^\gamma } = \sum\limits_{i = 1}^{2{n^\gamma }} {\xi _i^\gamma f_i^\gamma {\text{ }}} $

(6)

相应的结构位移为

${U^\gamma } = {\text{ }}\sum\limits_{i = 1}^{2{n^\gamma }} {\xi _i^\gamma u_i^\gamma } $

(7)

其中

${{\xi }_{2}}i-{{1}^{\gamma }}=h_{i}^{\gamma }\cos (\theta _{i}^{\gamma })$

(8)

$\xi _2i^\gamma = h_i ^\gamma \sin (\theta _i ^\gamma )$

(9)

由式(6)和式(7)可得单工况下结构的柔度为

$\begin{align} & {{c}^{\gamma }}={{({{F}^{\gamma }})}^{T}}{{U}^{\gamma }}={{(\sum\limits_{i=1}^{2{{n}^{\gamma }}}{\xi _{i}^{\gamma }f_{i}^{\gamma }})}^{T}}(\sum\limits_{i=1}^{2{{n}^{\gamma }}}{\xi _{i}^{\gamma }u_{j}^{\gamma }})= \\ & \sum\limits_{i=1}^{2{{n}^{\gamma }}}{\sum\limits_{j=1}^{2{{n}^{\gamma }}}{\xi _{i}^{\gamma }\xi _{j}^{\gamma }{{(f_{i}^{\gamma })}^{T}}u_{j}^{\gamma }}}=\sum\limits_{i=1}^{2{{n}^{\gamma }}}{\sum\limits_{j=1}^{2{{n}^{\gamma }}}{\xi _{i}^{\gamma }\xi _{j}^{\gamma }{{c}_{i}}{{j}^{\gamma }}}} \\ \end{align}$

(10)

其中$c_ij^\gamma = ({f}_i^\gamma )^{T}{u}_j ^\gamma$表示载荷{f}_i ^\gamma 在载荷${f}_j ^\gamma$所引起的位移${u}_j^\gamma $上所作之功.

引入$\xi _i ^\gamma$ 的二阶和四阶中心矩

$\xi _ij^\gamma = E(\xi _i^\gamma \xi _j^\gamma ) \ \ \ \ (i,j = 1,2,\cdots ,2n^\gamma)$

(11)

$\xi _ijkl^\gamma = E(\xi _i^\gamma \xi _j^\gamma \xi _k^\gamma \xi _l^\gamma ) \ \ \ \ (i,j,k,l = 1,2,\cdots ,2n^\gamma)$

(12)

通常$\xi _ij^\gamma $与$\xi _ijkl^\gamma$ 的解析表达式很难给出，但可以通过蒙特卡罗方法对它们的值进行精确估计.

单一工况下结构柔度的均值和方差分别为

$\mu ({{c}^{\gamma }})=\sum\limits_{i,j-1}^{2{{n}^{\gamma }}}{\xi _{i}^{\gamma }\xi _{j}^{\gamma }{{c}_{i}}{{j}^{\gamma }}}=\sum\limits_{i,j-1}^{2{{n}^{\gamma }}}{{{\xi }_{i}}{{j}^{\gamma }}{{c}_{i}}{{j}^{\gamma }}}$

(13)

$\begin{align} & {{\sigma }^{2}}({{c}^{\gamma }})=E({{({{c}^{\gamma }})}^{2}})-{{\mu }^{2}}({{c}^{\gamma }})= \\ & \sum\limits_{i,j,k,l}^{2{{n}^{\gamma }}}{({{\xi }_{i}}jk{{l}^{\gamma }}-{{\xi }_{i}}{{j}^{\gamma }}{{\xi }_{k}}{{l}^{\gamma }}){{c}_{i}}{{j}^{\gamma }}{{c}_{k}}{{l}^{\gamma }}} \\ \end{align}$

(14)

由此可知,在$\xi_ij$与$\xi_ijkl$已知的情况下,不需要对所有工况都进行结构分析,只要获得结构在确定载荷工况$f_1^\gamma ,f_2^\gamma , \cdots ,{f_2}{n^\gamma }^\gamma $作用下的位移向量$u_1^\gamma ,u_2^\gamma , \cdots ,{u_2}{n^\gamma }^\gamma $,就可以求得结构柔度的均值与方差.注意到$c_ij^\gamma$为向量$f_i^\gamma $与${u}_j^\gamma$的内积,其计算量非常小,与有限元分析相比可以忽略不计.

1.3 多工况下结构柔度均值和方差的计算

根据以上部分的推导，分别将式(13)和式(14)代入式(3)和式(4)即得到M个工况下结构整体柔度的均值和方差

$\begin{align} & \mu (c)=\sum\limits_{\gamma =1}^{^{M}}{{{\varpi }^{\gamma }}\mu ({{c}^{\gamma }})}=\sum\limits_{\gamma =1}^{^{M}}{{{\varpi }^{\gamma }}}\sum\limits_{i,j=1}^{^{M}}{{{\xi }_{i}}{{_{j}}^{\gamma }}{{c}_{i}}{{_{j}}^{\gamma }}}= \\ & \sum\limits_{\gamma =1}^{^{M}}{{{\varpi }^{\gamma }}}\sum\limits_{i,j=1}^{^{M}}{{{\varpi }^{\gamma }}{{\xi }_{i}}{{j}^{\gamma }}{{c}_{i}}{{j}^{\gamma }}} \\ \end{align}$

(15)

$\begin{align} & {{\sigma }^{2}}(c)=\sum\limits_{\gamma =1}^{M}{{{({{\varpi }^{\gamma }})}^{2}}{{\sigma }^{2}}({{c}^{\gamma }})}= \\ & \sum\limits_{\gamma =1}^{M}{{{({{\varpi }^{\gamma }})}^{2}}}\sum\limits_{i,j,k,l=1}^{M}{({{\xi }_{i}}{{_{jkl}}^{\gamma }}-{{\xi }_{i}}{{_{j}}^{\gamma }}{{\xi }_{k}}{{_{l}}^{\gamma }})=} \\ & \sum\limits_{\gamma =1}^{M}{\sum\limits_{i,j,k,l=1}^{M}{{{({{\varpi }^{\gamma }})}^{2}}({{\xi }_{i}}{{_{jkl}}^{\gamma }}-{{\xi }_{i}}{{_{j}}^{\gamma }}{{\xi }_{k}}{{_{l}}^{\gamma }}){{c}_{i}}{{_{j}}^{\gamma }}{{c}_{k}}{{_{l}}^{\gamma }}}} \\ \end{align}$

(16)

从而目标函数可以表示为

$\begin{align} & J=\alpha \mu (c)+\beta \sigma (c)=\alpha \sum\limits_{\gamma =1}^{M}{\sum\limits_{i,j=1}^{2{{n}^{\gamma }}}{{{\varpi }^{\gamma }}{{\xi }_{i}}{{j}^{\gamma }}{{c}_{i}}{{j}^{\gamma }}+}} \\ & \beta \sqrt{\sum\limits_{\gamma =1}^{M}{\sum\limits_{{{_{i}}_{,}}{{_{j,k,l}}_{=1}}}^{2{{n}^{\gamma }}}{{{({{\varpi }^{\gamma }})}^{2}}({{\xi }_{i}}jk{{l}^{\gamma }}-{{\xi }_{i}}{{j}^{\gamma }}{{\xi }_{k}}{{l}^{\gamma }}){{c}_{i}}{{j}^{\gamma }}{{c}_{k}}{{l}^{\gamma }}}}} \\ \end{align}$

(17)

1.4 结构灵敏度分析

本文采用移动渐近线方法 ^[36]进行拓扑优化问题的求解，因此需要计算结构的灵敏度信息.首先根据式(15)和式(16)可求 得结构柔度均值与标准差对密度$\rho_e$ 的偏导数表达式，分别为

$\frac{\partial \mu (c)}{\partial {{\rho }_{e}}}=\text{ }\sum\limits_{\gamma =1}^{M}{\sum\limits_{i,j}^{2{{n}^{\gamma }}}{{{\varpi }^{\gamma }}{{\xi }_{i}}{{j}^{\gamma }}\frac{\partial {{c}_{i}}{{_{j}}^{\gamma }}}{\partial {{\rho }_{e}}}}}$

(18)

$\begin{align} & \frac{\partial \sigma (c)}{\partial {{\rho }_{e}}}=\frac{1}{2\sigma (c)}\frac{\partial {{\sigma }^{2}}(c)}{\partial {{\rho }_{e}}}= \\ & \frac{1}{\sigma (c)}\sum\limits_{\gamma =1}^{^{2{{n}^{\gamma }}}}{\sum\limits_{i,j}^{^{2{{n}^{\gamma }}}}{{{({{\varpi }^{\gamma }})}^{2}}({{\xi }_{i}}jk{{l}^{\gamma }}-{{\xi }_{i}}{{_{j}}^{\gamma }}{{\xi }_{k}}{{_{l}}^{\gamma }}){{c}_{k}}{{l}^{\gamma }}]\frac{\partial {{c}_{i}}{{_{j}}^{\gamma }}}{\partial {{\rho }_{e}}}}1} \\ \end{align}$

(19)

引入变量

${{w}_{i}}{{_{j}}^{\gamma }}=\alpha {{\varpi }^{\gamma }}{{\xi }_{i}}{{_{j}}^{\gamma }}+\frac{\beta }{\sigma (c)}\sum\limits_{\gamma =1}^{^{2{{n}^{\gamma }}}}{{{({{\varpi }^{\gamma }})}^{2}}({{\xi }_{i}}{{_{j}}_{kl}}-{{\xi }_{i}}_{j}{{\xi }_{k}}_{l}){{c}_{k}}{{_{l}}^{\gamma }}}$

(20)

则目标函数对$\rho _e$ 的偏导数，即目标函数灵敏度为

$\begin{align} & \frac{\partial J}{\partial {{\rho }_{e}}}=\alpha \frac{\partial \mu (c)}{\partial {{\rho }_{e}}}+\beta \frac{\partial \sigma (c)}{\partial {{\rho }_{e}}}= \\ & \sum\limits_{\gamma =1}^{M}{\sum\limits_{i,j=1}^{2{{n}^{\gamma }}}{{{w}_{i}}{{j}^{\gamma }}\frac{\partial {{c}_{i}}{{j}^{\gamma }}}{\partial {{\rho }_{e}}}}}= \\ & \sum\limits_{\gamma =1}^{M}{\sum\limits_{i,j=1}^{2{{n}^{\gamma }}}{{{w}_{i}}{{_{j}}^{\gamma }}\frac{\partial ({{(f_{i}^{\gamma })}^{T}}u_{j}^{\gamma })}{\partial {{\rho }_{e}}}}} \\ \end{align}$

(21)

由式(21)可知要求出该灵敏度需要求$\sum\limits_{\gamma =1}^{M}{{{(2{{n}^{\gamma }})}^{2}}}$次偏导数，为了减少运算量对灵敏度公式进行下面的推导.

由式(20)可知由$w_ij^\gamma$ 组成的矩阵${W}^\gamma$是一个对称矩阵，可以对其进行对角化分解，即存在$2n^\gamma$ 个实数$\{ \lambda _1^\gamma , \cdots ,{\lambda _2}{n^\gamma }^\gamma \} $以及一个正交矩阵${Q^\gamma }$使得

${({Q^\gamma })^T}{W^\gamma }{Q^\gamma } = diag\{ \lambda _1^\gamma ,\lambda _2^\gamma , \cdots ,{\lambda _2}{n^\gamma }^\gamma \} $

(22)

为了推导方便，分别取下面4个矩阵

$f_b^\gamma = [f_1^\gamma ,f_2^\gamma , \cdots ,{f_2}{n^\gamma }^\gamma ]$

(23)

$u_b^\gamma = [u_1^\gamma ,u_2^\gamma , \cdots ,{u_2}{n^\gamma }^\gamma ]$

(24)

$F_b^\gamma = f_b^\gamma {Q^\gamma } = [F_1^\gamma ,F_2^\gamma , \cdots ,{F_2}{n^\gamma }^\gamma ]$

(25)

$U_b^\gamma = u_b^\gamma {Q^\gamma } = [U_1^\gamma ,U_2^\gamma , \cdots ,{U_2}{n^\gamma }^\gamma ]$

(26)

假设$q_ij^\gamma$ 是${Q}^\gamma$ 第i行第j列的一个元素，${\pmb F}_i^\gamma 和{U}_j^\gamma$ 可以进一步表示为

$F_{i}^{\gamma }=\text{ }\sum\limits_{j=1}^{^{2{{n}^{\gamma }}}}{{{q}_{i}}{{_{j}}^{\gamma }}f_{j}^{\gamma }}$

(27)

$U_{i}^{\gamma }=\text{ }\sum\limits_{j=1}^{2{{n}^{\gamma }}}{{{q}_{i}}{{j}^{\gamma }}u_{j}^{\gamma }}$

(28)

构造函数

${{({{J}^{\gamma }})}^{*}}=\text{ }\sum\limits_{i,j=1}^{2{{n}^{\gamma }}}{\lambda _{i}^{\gamma }{{(F_{i}^{\gamma })}^{T}}U_{i}^{\gamma }}$

(29)

可以证明^[30]

$\frac{\partial {{({{J}^{\gamma }})}^{*}}}{\partial {{\rho }_{e}}}=\sum\limits_{i,j=1}^{2{{n}^{\gamma }}}{{{w}_{i}}{{j}^{\gamma }}\frac{\partial {{c}_{i}}{{j}^{\gamma }}}{\partial {{\rho }_{e}}}}$

(30)

将式(29)和式(30)代入式(21)可得

$\frac{\partial J}{\partial {{\rho }_{e}}}=\sum\limits_{\gamma =1}^{M}{\sum\limits_{i=1}^{^{2{{n}^{\gamma }}}}{\lambda _{i}^{\gamma }\frac{\partial ({{(F_{i}^{\gamma })}^{T}}U_{i}^{\gamma })}{\partial {{\rho }_{e}}}}}$

(31)

由式(31)可知只需要求$\sum\limits_{\gamma =1}^{M}{2{{n}^{\gamma }}}$ 次偏导数，而矩阵${U}_i^\gamma$ 是向量$u_j^\gamma (j = 1,2, \cdots ,2{n^\gamma })$的线性组合，其计算量很小，这样由式(31)求解敏度信息运算量大大降低.

利用拉格朗日乘子法容易证明

$\begin{align} & \frac{\partial ({{(F_{i}^{\gamma })}^{T}}U_{i}^{\gamma })}{\partial {{\rho }_{e}}}=-{{(U_{i}^{\gamma })}^{T}}\frac{\partial K}{\partial {{\rho }_{e}}}U_{i}^{\gamma } \\ & =-p\rho _{e}^{p-1}({{E}_{0}}-{{E}_{\min }}){{({{U}_{i}}{{e}^{\gamma }})}^{T}}k_{e}^{0}{{U}_{i}}{{e}^{\gamma }} \\ \end{align}$

(32)

其中${k}_e^0$ 表示单位模量单元e充满材料时的刚度阵，${U_i}{_e^\gamma } $为结构在载荷${F}_i^\gamma$ 作用下单元e的位移向量. 将式(32)代入式(31)可得目标函数对设计变量的灵敏度

$\frac{\partial J}{\partial {{\rho }_{e}}}=-p\rho _{e}^{p-1}({{E}_{0}}-{{E}_{\min }})\sum\limits_{\gamma =1}^{M}{\sum\limits_{i=1}^{^{2{{n}^{\gamma }}}}{\lambda _{i}^{\gamma }{{({{U}_{i}}{{_{e}}^{\gamma }})}^{T}}k_{e}^{0}}}U_{ie}^{\gamma }$

(33)

体积约束对于设计变量的灵敏度为

$\frac{\partial g}{\partial {{\rho }_{e}}}={{v}_{e}}$

(34)

1.5 灵敏度过滤方法

为了获得计算数值稳定性^[37]，本文采用灵敏度过滤的方法，将目标函数的灵敏度${\partial J} / {\partial \rho _e }$修改为如下形式

$\frac{\partial {{J}^{*}}}{\partial {{\rho }_{e}}}=\frac{1}{\max (\upsilon ,{{\rho }_{e}})\sum\limits_{{}}{{{N}_{e}}{{H}_{e}}i}}\sum\limits_{{{_{i\in }}_{N}}}{{{N}_{e}}{{H}_{e}}i{{\rho }_{i}}\frac{\partial J}{\partial {{\rho }_{i}}}}$

(35)

其中$N_e$ 是与单元e的中心距不大于过滤半径$r_\min $ 的所有单元的集合. $\upsilon$ 是为了避免除以零而设的一个特别小的正数. 权重因子$H_ei$ 定义为

$H_ei = \max (0,r_\min - d(e,i))$

(36)

其中d(e,i)是单元e和i的中心的距离.

1.6 多工况线性结构稳健拓扑优化算法

综合上述推导公式可以求出多工况、载荷不确定条件下目标函数值及其灵敏度信息，在计算过程中，本文采用移动渐近线方法 算法对设计变量进行更新. 多工况线性结构稳健拓扑优化设计具体过程可以表示如下：

(1) 设计变量和参数初始化.

(2) 用蒙特卡罗法计算每个工况下不确定载荷的二阶和四阶中心距$\xi _ij^\gamma$与$\xi _ijkl^\gamma$ .

(3) 开始循环迭代，直到收敛为止：

① 解有限元方程

${K u}_i^\gamma = {f}_i^\gamma \ \ (i = 1,2,\cdots,2n^\gamma; \ \gamma = 1,2.\cdots,M)$

② 计算$c_ij^\gamma \ (i,j = 1,2,\cdots,2n^\gamma ; \ \gamma = 1,2,\cdots,M)$

③ 通过式(15)~式(17)求得柔度的均值和方差，并计算目标函数值

④ 构造矩阵${W}^\gamma$ 并计算$\{\lambda _1 ,\lambda_2,\cdots ,\lambda _2n^\gamma \}$以及它们的相关特征向量来组成矩阵{Q}

⑤ 由式(28)计算${U}_i \ (i = 1,2,\cdots,2n^\gamma )$

⑥ 通过式(33)和式(34)计算目标函数和约束的敏度信息

⑦ 用式(35)来过滤目标函数的灵敏度

⑧ 用移动渐近线方法来更新设计变量并判断优化过程是否收敛. 若收敛，则迭代结束；否则转到①重新进行有限元分析等直至优化过程收敛.

从以上流程中可以看出，在M个工况、载荷不确定条件下进行结构稳健拓扑优化，如果每个工况承受$n^\gamma$ 个不确定载荷，则每次迭代需要对$\sum\limits_{\gamma =1}^{M}{2{{n}^{\gamma }}}$个工况进行有限元分析并进行$\sum\limits_{\gamma =1}^{M}{2{{n}^{\gamma }}}$ 次灵敏度分析，确定性分析时则仅需要对M个工况进行有限元分析并进行M次灵敏度分析. 因此载荷不确定条件下的结构稳健拓扑优化算法计算量要大于确定性条件下的算法.

2 数值算例

本节以两个算例来说明多工况、载荷不确定条件下的最优结构与确定条件下的最优结构不同,并且前者要比后者具有更好的稳健性.两个算例均采用单位双线性正方形单元对设计域进行离散.蒙特卡罗方法计算$\xi_ij^\gamma$与$\xi_ijkl^\gamma$时样本数目为1×10⁸,保证了其精度.实体与空单元的弹性模量分别为$E_0=1,E_\min=1\times 10^{-9}$,材料的泊松比均为\mu=0.3,惩罚因子p=3.优化结束的准则为每个单元的伪密度变化量均不超出0.01.

2.1 悬臂梁结构的拓扑优化

如图1所示,某悬臂梁结构的设计域为30×30的正方形区域,左边固定.工况一为结构右上角作用一个竖直向上的不确定载荷$P_1$,工况二为结构右下角作用一个竖直向下的不确定载荷$P_2$.载荷$P_1,P_2$的大小和方向均符合正态分布,大小的均值分别为$\mu _h1=1\,N,\mu _h2=1\,N,$标准差分别为$\sigma _h1=0.2\,N,\sigma _h2=0.05\,N;$方向的均值分别为${\mu _\theta }1 = \pi /2,{\mu _\theta }2 = - \pi /2$,标准差均为$\pi {\text{/8}}$.材料的许用体积为设计域的40%.

将整个设计区域离散为8 100(90×90)单元,将均值和方差的权重系数$\alpha$和$\beta$设置为0.5,两个工况的权重系数$\varpi^1$和$\varpi^2$分别设置为1.对多工况确定载荷情况、载荷大小不确定情况、载荷方向不确定情况、大小和方向均不确定情况进行优化,最优结构如图2.

由于结构的设计域、边界条件以及两个工况的确定载荷作用均上下对称,因此优化后的结构也均上下对称.其中只考虑载荷大小不确定性时,由于载荷大小的标准差分别为$\sigma$ _h1=0.2 N比$\sigma$ _h2=0.05 N要大,克服力P₁占结构的比例较大,所以获得构型的上边框要比下边框粗,承受力P₁的杆架要比承受力P₂的杆架粗.而只考虑载荷方向不确定性时获得的构型与确定载荷作用时相比上下边框加粗,这样能承受更大的水平力.同时考虑载荷大小与方向不确定性时获得的结构构型综合了载荷大小不确定性优化和载荷方向不确定性优化构型的特点,即上下边框同时加粗且上边框要比下边框粗.对于本例而言,由于两个工况载荷大小的标准差不同,导致结构构型的不对称,结构承受能力更偏向于标准差更大的工况,使得结构更加合理.而载荷方向不确定时,其可能含有水平分量,加粗的上下边框具有较好的承受水平载荷的能力,这样的结构虽然竖直承载能力有所下降,但是水平承载能力加强了,结构整体上更加稳健.

另外,为了说明工况对优化结果的影响,对载荷P₁和P₂同时作用的情况进行优化,优化结果如图3所示.由图2和图3对比可知,载荷P₁和P₂作为多工况作用时的最优结构和载荷P₁和P₂同时作用时的最优结构明显不同.

结构稳健拓扑优化过程收敛曲线如图4所示，横轴表示迭代次数，纵轴表示目标函数.各种情况下本文算法都能很好的收敛，说明本文算法的稳定性.

从表1的优化结果中可以看出，确定性载荷条件下得到的结构在考虑载荷不确定性时的柔度均值和标准差一般高于不确定载荷条件下得到的结构的柔度均值和标准差，这证明了本文算法的有效性.注意到对于本算例，载荷大小不确定条件下与确定条件下的最优结构的目标函数值相差不超过0.5%，因此，载荷不确定性对本算例影响较小.

为了研究权重系数$\alpha$和$\beta$对结构拓扑优化结果的影响,本文将$\alpha$在[0,1]之间每隔0.1取值一次进行大小和方向均不确定的优化,所得均值和标准差如图5所示.由可知,随着$\alpha$的增加,方差在目标函数中的作用增强,均值的作用降低.

2.2 米歇尔结构的拓扑优化

如图6所示,问题的设计域为240×120的矩形区域,底部两端均简支,工况一为结构下端中心作用一个竖直向下的不确定载荷P₁,工况二为结构下端两侧各1/3处分别作用一个竖直向下的不确定载荷P₂和P₃.3个载荷的大小h_i和作用方向与水平方向夹角\theta_i均服从正态分布,它们的均值$\mu $_hi,$\mu$ _{$\theta$ i}和标准差$\sigma$_hi,$\sigma$ _{$\theta$ i}如表2所示.材料的许用体积为设计域的30%.

将整个设计区域离散为28 800(240×120)单元,将均值和方差的权重系数$\alpha$和$\beta$设置为0.5,两个工况的权重系数$\varpi^1$和$\varpi^2$分别设置为1.对多工况确定载荷情况、载荷大小不确定情况、载荷方向不确定情况、大小和方向同时不确定情况进行优化,结果如图7.

由于结构的设计域、边界条件以及两个工况的作用载荷均左右对称,因此确定载荷下优化后的结构也均左右对称.其中只考虑载荷大小不确定性时获得的构型与载荷确定条件下获得的结构构型基本一致.只考虑载荷方向不确定性时和同时考虑方向和大小不确定性时获得的构型与确定性时获得的结构构型相比发生了明显的变化,结构下方呈现出连接左右两个端点以及3个载荷作用点的杆件,这是由于载荷作用方向的不确定性产生了水平方向分量,而水平方向杆件具有较好的承受水平方向载荷的能力;这样的结构虽然竖直方向承载能力有所下降,但是水平方向承载能力加强了,结构整体上更加稳健.只考虑载荷方向不确定性时和同时考虑方向和大小不确定性时获得的构型也有区别,即中间立柱消失而两边立柱增粗.这是因为在载荷方向不确定的基础上加入了大小不确定性使得在水平方向力进一步增加,需要构型能承受更大的水平方向力.对本例而言,由于载荷方向不确定性直接导致了构型的变化,所以其对结构影响要比大小不确定性的影响更大.

结构稳健拓扑优化过程收敛曲线如图8所示,各种情况下本文算法都能很好的收敛,说明本文算法的稳定性.优化结果如表3所示,可以看出,确定性载荷条件下得到的结构在考虑载荷不确定性时的柔度均值和标准差一般高于不确定载荷条件下得到的结构的柔度均值和标准差,这证明了本文算法的有效性.

3 结论

本文研究了多工况、载荷不确定条件下的连续体结构拓扑优化问题，以最小化结构柔度均值与标准差的加权和为目标，给出了采用载荷大小和作用方向的概率分布表示载荷不确定性时的结构拓扑优化算法. 本文的主要结论：

(1)基于线性结构的位移叠加原理给出了多工况、载荷不确定条件下结构柔度均值与方差计算方法；

(2)给出了多工况、载荷不确定条件下结构目标函数灵敏度分析方法和结构拓扑优化算法；

(3)应用两个数值算例验证了算法的稳定性和有效性；

(4)研究表明，随着标准差权重$ \alpha$的增加，标准差在目标函数中的作用增强，均值的作用降低.

致 谢 衷心感谢瑞典皇家理工学院的Krister Svanberg教授提供移动渐近线法(MMA)的代码.