引 言

拓扑优化作为一种有效的自动化设计工具, 已广泛应用于各种结构的设计问题. 常用的拓扑优化方法主要包括SIMP (solid isotropic material with penalization)法^[1-2]、ESO (evolutionary structural optimization)法^[3]、水平集法^[4-5]、智能拓扑优化方法^[6-7]等. 此外, Guo等^[8-12]提出了一种基于移动可变形组件/孔洞的拓扑方法, 以可移动可变形组件为基本元件进行结构拓扑优化设计. 近期, Luo等^[13-16]提出了一种基于材料场级数展开 (material-field series-expansion, MFSE)的拓扑优化方法, 该方法能够有效避免结构中的棋盘格现象及网格依赖性问题, 并且能够大幅减少设计变量个数, 提高计算效率.

在结构的加工制造及使用过程中, 不可避免会存在各种不确定性, 比如加工误差导致的几何不确定性, 工作时的载荷不确定性等. 这些不确定性会使结构产生性能波动和关键指标降低, 严重影响结构的安全性. 现有处理不确定性的方法主要是基于概率框架进行的. 随着概率可靠性优化方法的发展, 一些实用的求解策略如功能度量法^[17-18], 序列优化方法^[19]及概率神经网络^[20]等都极大地促进了概率可靠性方法在实际工程中的应用. 此外, 实际工程中的许多不确定性, 如分布荷载、几何尺寸、材料属性等, 其在空间不同位置上的实现是随空间位置的变化而变化的, 这属于“不确定场”问题. 目前, 常用来描述不确定场的模型为随机场模型^[21-23]. 将制造误差引起的几何缺陷描述为随机阈值场模型, Kang等^[22]提出了一种考虑几何空间不确定性的概率可靠性拓扑优化方法. 通过将非侵入式PCE (the polynomial chaos expansion)方法与设计灵敏度分析相结合, Keshavarzzadeh等^[23]提出了一种随机场几何不确定性下的可靠性拓扑优化的系统方法.

针对工程实际中广泛存在的未知但有界不确定性, 由于缺少大量样本数据信息, 难以获得不确定性准确地概率分布特征, 因此概率不确定性理论不再适用. 作为概率可靠性理论的有效补充, 研究人员提出了许多描述参数不确定性的非概率模型, 比如模糊模型^[24-25]、证据理论^[26-27]、区间模型及凸模型^[28-34]等. Sofi等^[30]通过结合区间运算和安全系数方法, 将非概率可靠性指标描述为一个区间变量. Pantelides等^[33]提出了一种反优化技术进行考虑载荷不确定性的非概率可靠性的优化设计. 将结构边界长度变化描述为凸模型, Luo等^[34]完成了考虑结构几何不确定性的非概率可靠性优化设计, 优化后的结构具有更高的可靠性.

至于少样本的场不确定性问题, 由于样本数量有限, 因此随机场模型也不再适用. 通过引入EUI(external unit interval)变量, Muscolino等^[35-36]提出了一种区间场方法来量化未知但有界的不确定场. 基于空间相关性的数学定义和非概率级数展开方法, Luo等^[37]提出了一种处理有限样本下不确定场问题的有界场模型. 基于有界场模型, Zhan等^[38]提出了一种非概率可靠性指标进行结构场不确定性下的非概率可靠性评估, 并进一步完成了不确定载荷场作用下的结构非概率可靠性优化设计研究^[39].

文献调研显示, 尽管基于参数的非概率可靠性优化方法已应用于结构几何不确定性分析, 但针对工程实际中考虑少样本及空间变化特性的几何不确定性, 仍缺乏合理的可靠性优化模型. 因此, 本文将采用阈值技术^[40]进行具有空间变化特性的几何不确定性表征, 并将Heaviside过滤函数中的阈值$\eta $假定为有界不确定阈值场, 建立非概率有界场模型, 进而完成考虑结构几何不确定性的非概率可靠性拓扑优化设计.

本文的具体安排如下: 首先进行不确定场的非概率描述, 即将不确定场描述为非概率有界场模型; 第2节进行结构几何不确定性的描述, 即通过不确定阈值场来表示, 进而描述为有界场模型; 第3节为建立结构的非概率可靠性优化模型; 在第4节推导了模型的灵敏度信息, 并采用移动渐近线法 (method of moving asymptotes)^[41]求解优化问题; 最后通过2个数值算例验证了模型的有效性.

1.   不确定场的非概率描述

在本文中, 结构的几何不确定性将通过非概率有界场模型^[37]来描述. 假定设计域${\varOmega _{{\rm{dom}}}}$内存在一个有界不确定场$Z({\boldsymbol{x}}) \in \left[ {{Z^{{\rm{low}}}},{Z^{{\rm{upp}}}}} \right],{\boldsymbol{x}} \in {\varOmega _{{\rm{dom}}}}$,其中${Z^{{\rm{low}}}}$和${Z^{{\rm{upp}}}}$分别表示不确定场波动的上下界限. 将不确定场$Z({\boldsymbol{x}})$离散为$N$个观察点, 同时将设计域${\varOmega _{{\rm{dom}}}}$离散为$N$个单元. 在设计域内, 不确定场$Z({\boldsymbol{x}})$的观察点与单元的中心一一对应. 对不确定场$Z({\boldsymbol{x}})$进行标准化, 得

其中${Z_0}({\boldsymbol{x}})$为标准化的不确定场, 且满足$ - 1 \leqslant {Z_0}({\boldsymbol{x}}) \leqslant 1$.

考虑到在实际工程中结构的不确定边界是连续变化的, 因此假定有界不确定场$Z({\boldsymbol{x}})$的空间波动具有一定相关性. 在本文中, 设计域${\varOmega _{{\rm{dom}}}}$内任意两个观察点${{\boldsymbol{x}}_a}$和${{\boldsymbol{x}}_b}$对应的不确定场的相关性通过相关函数$R({{\boldsymbol{x}}_a},{{\boldsymbol{x}}_b})$来描述, $R({{\boldsymbol{x}}_a},{{\boldsymbol{x}}_b})$的表达式为

其中, 符号$\left\| \cdot \right\|$表示2范数. L是不确定场$Z({\boldsymbol{x}})$的相关长度, 它用来控制不确定场的空间波动程度. 当相关长度L较小时, 不确定场的空间波动较为剧烈, $L = 0$表示所有观察点处的变量都是不相关的; 当相关长度L较大时, 波动较为平缓, $L \to + \infty $表示不确定场完全相关, 即所有观察点处的不确定场值$Z({{\boldsymbol{x}}_i})$均相等.

设计域内各观察点之间的相关性构成了相关矩阵R, 因此R表示为

基于非概率级数扩展, 标准化的不确定场${Z_0}({\boldsymbol{x}})$可表示为

其中, $ {\lambda _j} $ 和 ${{\boldsymbol{\varphi }}_j}$ 为相关矩阵${\boldsymbol{R}}$的特征值和特征向量, 它们之间满足${\boldsymbol{R}}{{\boldsymbol{\varphi }}_j} = {\lambda _j}{{\boldsymbol{\varphi }}_j}{\text{ }}\left( {j = 1,2,\cdots,N} \right)$. 向量${{\boldsymbol{R}}_D}\left( {\boldsymbol{x}} \right) = \left[ {R\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_1}} \right),R\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_2}} \right),\cdots,R\left( {{\boldsymbol{x}},{{\boldsymbol{x}}_N}} \right)} \right]$. $ {\xi _j} $为非概率不确定系数, 向量${\boldsymbol{\xi }} = \left[ {{\xi _1},{\xi _2},\cdots,{\xi _N}} \right]$.

在式(4)中, 特征值$ {\lambda _j} $为降序排列, 特征值越小的项对不确定场的贡献值越小. 因此, 为提高模型的计算效率, 可对式(4)进行截断, 仅保留前$M$项, 即通过M个不确定系数($ {\xi _j}\left( {j = 1,2,\cdots,M} \right) $)来描述不确定场. 截断公式可表示为

其中, $\alpha $是一个很小的值. 在本文中, $\alpha $ = $1.0 \times {10^{ - 6}}$.

对式(4)进行截断后并代入式(1)中, 可得到不确定场的表达式为

已知标准化的不确定场$ {Z_0}({\boldsymbol{x}}) $的变化范围为${Z}_{0}({\boldsymbol{x}})\in \left[-1, 1\right]$, 得

对式(7)两边平方, 并写成向量的形式, 可以得到

在式(8)中, ${{\boldsymbol{\varLambda }}^{ - 1/2}}{\text{ = }}{\rm{diag}}\left( \dfrac{1}{{\sqrt {{\lambda _1}} }},\dfrac{1}{{\sqrt {{\lambda _2}} }},\cdots, \dfrac{1}{{\sqrt {{\lambda _M}} }} \right)$, ${\boldsymbol{\xi }} = \left[ {{\xi _1},{\xi _2},\cdots,{\xi _M}} \right]$, ${\boldsymbol{\varPhi }}$表示一个矩阵, 它包含了$M$个特征向量${{\boldsymbol{\varphi }}_j}$ .

引入${{\boldsymbol{W}}_i} = {{\boldsymbol{\varLambda }}^{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. } 2}}}{{\boldsymbol{\varPhi }}^{{\rm{T}}} }{{\boldsymbol{R}}_D}\left( {{{\boldsymbol{x}}_i}} \right){{\boldsymbol{R}}_D}{\left( {{{\boldsymbol{x}}_i}} \right)^{{\rm{T}}} }{\boldsymbol{\varPhi }}{{\boldsymbol{\varLambda }}^{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. } 2}}}$, 式(8)可进一步表示为

2.   几何不确定性的描述

在现有的拓扑优化设计中, Heaviside过滤技术已得到了广泛的应用. 与敏度法和密度法过滤相比, 优化的结构经Heaviside过滤可获得清晰的结构边界. 本文用到的Heaviside过滤函数定义为

其中, ${\rho _e}$为结构中第 $e$个单元的相对密度, $\rho _e^H$为${\rho _e}$经Heaviside过滤后的相对密度值; $\eta \in [0,1]$为该函数的阈值, $\delta $为光滑参数. 光滑参数$\delta = 20$时, 不同阈值$\eta $对应的Heaviside近似函数如图1所示.

图  1  光滑系数$\delta = 20$时, 不同阈值$\eta $下Heaviside函数

Figure  1.  The Heaviside function under different thresholds $\eta $ when smoothing parameter $\delta = 20$

下载: 全尺寸图片 幻灯片

几何结构的不确定性可通过阈值技术^[40]来表征, 即对结构设计域中的每个单元选择不同的阈值$\eta $进行Heaviside过滤. 在本文研究中, 将式(10)中的阈值$\eta $描述为一个有界不确定阈值场$\eta ({\boldsymbol{x}}) \in \left[ {\eta ^{{\rm{low}}}}, {\eta ^{{\rm{upp}}}} \right],{\boldsymbol{x}} \in {\varOmega _{{\rm{dom}}}}$,其中${\eta ^{{\rm{low}}}}$和${\eta ^{{\rm{upp}}}}$分别表示不确定阈值场波动的上下界限, 取值均在0 ~ 1之间. 基于第1节中的不确定场非概率描述(式(6)), 可得到不确定阈值场$\eta ({\boldsymbol{x}})$的表达式为

其中, 不确定系数${\xi _j}$的变化范围如式(9)所示.

为了进一步说明不确定阈值场$\eta ({\boldsymbol{x}})$对结构不确定性的影响, 图2以MBB梁为例分别给出了不确定阈值场$\eta ({\boldsymbol{x}})$和确定阈值$\eta ({\boldsymbol{x}}) \equiv 0.5$对结构边界的影响. 从图2中可以看出, 与阈值为恒定值($\eta ({\boldsymbol{x}}) \equiv 0.5$)相比, 当阈值描述为不确定场时, 结构的边界会存在一定的扰动(认为是几何不确定性).

图  2  阈值为定值($\eta ({\boldsymbol{x}}) \equiv 0.5$)及不确定阈值场$\eta ({\boldsymbol{x}})$下的MBB梁拓扑构型

Figure  2.  Topological structure of MBB beam considering the constant threshold $\eta ({\boldsymbol{x}}) \equiv 0.5$ and uncertain threshold field

下载: 全尺寸图片 幻灯片

3.   考虑几何不确定的结构非概率可靠性拓扑优化模型

3.1   考虑不确定阈值场的结构非概率可靠性指标

正如第2节描述的那样, 不确定阈值场$\eta ({\boldsymbol{x}})$ 会对结构的边界产生影响, 进而影响结构的性能. 因此, 结构性能可表示为不确定阈值场$\eta ({\boldsymbol{x}})$的函数. 在本文中, 结构的性能函数可表示为$C(\eta ) \leqslant {C^*}$,其中$C$为结构的柔顺性, ${C^*}$为给定的柔顺性约束. 令$g(\eta ) = {C^*} - C(\eta )$, 则$g(\eta )$称为极限状态函数.

由式(6)可知, 不确定阈值场$\eta ({\boldsymbol{x}})$是关于不确定系数${\boldsymbol{\xi }}$的函数. 因此, 在不确定系数${\boldsymbol{\xi }}$的空间中, 极限状态函数$g(\eta )$可进一步表示为$G({\boldsymbol{\xi }}) = g(\eta ({\boldsymbol{\xi }},{\boldsymbol{x}}))$.

依据有界不确定场的非概率可靠性指标的定义^[39], 极限状态函数$G({\boldsymbol{\xi }})$可将不确定系数${\boldsymbol{\xi }}$空间划分为可靠区和失效区, 如图3所示(为便于表示, 示意图中只考虑了2个不确定性系数${\xi _1}$和${\xi _2}$). 因此, 考虑几何不确定性的结构非概率可靠性指标定义为

图  3  有界场模型非概率可靠性指标示意图

Figure  3.  Schematic diagram of non-probabilistic reliability index for the bounded field model

下载: 全尺寸图片 幻灯片

即表示在可靠区域内所允许的最大不确定性. 其中, 式(11)的最优解${\beta ^*}$为非概率可靠性指标. ${\rm{jud}}(G({\boldsymbol{0}}))$ 用来判断可靠性指标${\beta ^*}$的正负, 具体的表达式为

3.2   基于MFSE模型的拓扑优化方法

为避免传统密度法拓扑优化中的棋盘格现象及网格依赖性问题, 本文将采用基于MFSE的拓扑优化方法^[13]. 在本方法中, 结构拓扑通过一个具有一定空间相关性的有界材料场函数$\varphi \left( {\boldsymbol{y}} \right) \in \left[ { - 1,1} \right],{\boldsymbol{y}} \in {\varOmega _{{\rm{dom}}}}$来描述. 经材料场级数展开, 并进行截断保留前${M_e}$(${M_e} \leqslant N$)项, 则场函数$\varphi \left( {\boldsymbol{y}} \right)$ 可表示为

其中, ${\kappa _k}$和${{\boldsymbol{\psi }}_k}$为相关矩阵${\boldsymbol{\varGamma }}$的特征值和特征向量. 相关矩阵${\boldsymbol{\varGamma }}$的表达式为

在式(13)中, ${\chi _k}$为不确定系数, ${\varGamma _D}\left( {\boldsymbol{y}} \right) = \left[ R\left( {\boldsymbol{y}}, {{\boldsymbol{y}}_1} \right),R\left( {{\boldsymbol{y}},{{\boldsymbol{y}}_2}} \right),\cdots,R\left( {{\boldsymbol{y}},{{\boldsymbol{y}}_N}} \right) \right]$. 在基于MFSE的拓扑优化方法中, ${\chi _k}$为设计变量.

类似于前面不确定场的处理方式, 依据有界不确定场函数$\varphi \left( {\boldsymbol{y}} \right) \in \left[ { - 1,1} \right]$的界限并引入符号${{\boldsymbol{H}}_m} = {{\boldsymbol{\kappa }}^{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. } 2}}}{{\boldsymbol{\varPsi }}^{{\rm{T}}} }{{\boldsymbol{\varGamma }}_D}\left( {{{\boldsymbol{y}}_m}} \right){{\boldsymbol{\varGamma }}_D}{\left( {{{\boldsymbol{y}}_m}} \right)^{{\rm{T}}} }{\boldsymbol{\varPsi }}{{\boldsymbol{\kappa }}^{{{ - 1} \mathord{\left/ {\vphantom {{ - 1} 2}} \right. } 2}}}$ , 则式(13)可转化为如下形式

其中, ${{ {\textit{χ}}}} = \left[ {{\chi _1},{\chi _2},\cdots,{\chi _{{M_e}}}} \right]$, ${\boldsymbol{\varPsi }}$为包含${M_e}$个特征向量的矩阵, ${{\boldsymbol{\kappa }}^{ - 1/2}}{\text{ = }}{\rm{diag}}\left( {\dfrac{1}{{\sqrt {{\kappa _1}} }},\dfrac{1}{{\sqrt {{\kappa _2}} }},\cdots,\dfrac{1}{{\sqrt {{\kappa _{{M_e}}}} }}} \right)$.

在本方法中, 假定结构的单元中心与不确定场函数$\varphi \left( {\boldsymbol{y}} \right)$的观察点${{\boldsymbol{y}}_m}$一一对应. 单元的插值函数可表示为

其中, ${E_0}$为实体结构的弹性模量, ${E_{\min }} = {10^{ - 9}}{E_0}$, 惩罚系数$q = 3$. $ {\rho ^H}\left( {{{\boldsymbol{y}}_m}} \right) $为$\rho \left( {{{\boldsymbol{y}}_m}} \right)$ 经Heaviside过滤后的值(式(10)), $\rho \left( {{{\boldsymbol{y}}_m}} \right) = \dfrac{{ {1 + \varphi \left( {{{\boldsymbol{y}}_m}} \right)} }}{2}$. 因此, 在本文中, $ {\rho ^H}\left( {{{\boldsymbol{y}}_m}} \right) $是关于阈值场$\eta \left( {\boldsymbol{x}} \right)$的函数,可表示为$ {\rho ^H}\left( {\eta ,{{\boldsymbol{y}}_m}} \right) $. 依据式(11)知, $ {\rho ^H}\left( {{{\boldsymbol{y}}_m}} \right) $也是关于不确定系数${\boldsymbol{\xi }}$的函数, $ {\rho ^H}\left( {{{\boldsymbol{y}}_m}} \right) $也可进一步表示为$ {\rho ^H}\left( {{\boldsymbol{\xi }},{{\boldsymbol{y}}_m}} \right) $ .

基于MFSE的结构柔顺性拓扑优化问题可表示为

其中, 设计变量为不确定系数向量${{ {\textit{χ}}}} = \left[ {{\chi _1},{\chi _2},\cdots,{\chi _{{M_e}}}} \right]$. ${\boldsymbol{F}}$为载荷向量, ${\boldsymbol{u}}$和${\boldsymbol{K}}$分别表示位移向量和结构总刚度阵. ${V_m}$表示第$m$个单元的体积, ${V^*}$ 为结构完全有材料时的体积, $f$为给定的结构体积分数约束.

3.3   结构非概率可靠性优化模型

在本文中, 结合MFSE优化模型(式(17))及考虑不确定场的可靠性指标(式(11)), 考虑结构几何不确定性的非概率可靠性拓扑优化问题可表示为: 在结构体积一定的情况下使结构的可靠性指标最大化, 进而来提高结构的可靠性, 优化列式如下所示

需要注意的是, 在式(18)的外层优化中, 很难准确地获得可靠性指标${\beta ^*}$对设计变量${{ {\textit{χ}}}}$的灵敏度信息. 因此, 可采用关心性能法^{[34, 39]}对式(18)进行等效变换. 基于关心性能法(详见文献[34,39]), 式(18)的非概率可靠性拓扑优化问题可以转化为

其中, 内层优化为获得在外层设计变量为${{ {\textit{χ}}}}$时的关心性能值$\sigma \left( {{{ {\textit{χ}}}},{{\boldsymbol{\xi }}^*}} \right)$, ${{\boldsymbol{\xi }}^*}$称为关心点. 确定了关心点${{\boldsymbol{\xi }}^*}$后便可得到此时对应的不确定阈值场$\eta ({\boldsymbol{x}})$(依据式(11))分布情况. 极限状态函数$ G\left( {{{ {\textit{χ}}}},{\boldsymbol{\xi }}} \right) $是结构柔顺性的函数, 表示为$ G\left( {{{ {\textit{χ}}}},{\boldsymbol{\xi }}} \right) = {C^*} - C\left( {{{ {\textit{χ}}}},{\boldsymbol{\xi }}} \right) $. 由于$ {C^*} $为一给定常数, 因此式(19)中外层优化的目标函数$ \max G\left( {{{ {\textit{χ}}}},{{\boldsymbol{\xi }}^*}} \right) $可等效表示为$ \min C\left( {{{ {\textit{χ}}}},{{\boldsymbol{\xi }}^*}} \right) $; 内层优化等效为$ \max C\left( {{{ {\textit{χ}}}},{\boldsymbol{\xi }}} \right) $. $\underline \beta $为给定的非概率可靠性指标下限值.

本文采用移动渐近线方法(MMA)^[41]求解该非概率可靠性优化问题. 优化模型的收敛准则定义为变量(内层为变量${\boldsymbol{\xi }}$, 外层为变量${{ {\textit{χ}}}}$)在相邻两步迭代中的最大变化值小于0.01.

4.   灵敏度分析

4.1   极限状态函数对不确定系数的灵敏度分析

本文考虑几何不确定性的非概率可靠性优化问题是基于梯度算法求解的, 因此灵敏度分析是必不可少的过程. 该优化模型(式(19))为嵌套优化, 首先求解内层优化的灵敏度信息, 即极限状态函数$ G\left( {{{ {\textit{χ}}}},{\boldsymbol{\xi }}} \right) $对不确定系数${\xi _j}$的灵敏度分析. 具体求解过程如下所示

平衡方程$ {\boldsymbol{Ku}} = {\boldsymbol{F}} $两边对不确定阈值场$\eta \left( {{{\boldsymbol{x}}_i}} \right)$求导, 并等式变换, 得

在本文中, 不确定场为阈值场, 其变化影响结构几何尺寸变换, 与结构所受载荷${\boldsymbol{F}}$无关. 因此, 式(21)中的$\dfrac{{{\text{d}}{\boldsymbol{F}}}}{{{\text{d}}\eta \left( {{{\boldsymbol{x}}_i}} \right)}} = 0$. 将式(21)代入式(20), 得

通过引入伴随向量${\boldsymbol{\gamma }}$(${\boldsymbol{K\gamma }} = \dfrac{{\partial G\left( {{{ {\textit{χ}}}},{\boldsymbol{\xi }}} \right)}}{{\partial {\boldsymbol{u}}}}$ ), 并在单元层级上表示, 式(22)可表示为

其中, 单元刚度阵对阈值的灵敏度$\dfrac{{{\text{d}}{{\boldsymbol{K}}_i}}}{{{\text{d}}\eta \left( {{{\boldsymbol{x}}_i}} \right)}}$可通过Heaviside函数(式(1))求得.

4.2   关心性能值对设计变量的灵敏度分析

优化模型中外层优化主要涉及的灵敏度为关心性能值$\sigma $对设计变量${\chi _k}$的灵敏度分析, 为实现优化模型的解耦, 假设关心性能点$ {\boldsymbol{\xi }} $对设计变量的灵敏度为0. 因此, 其表达式为

其中, 关心性能点${{\boldsymbol{\xi }}^*}$为优化模型内层的最优解.

依据链式法则, 式(24)表示为

将平衡方程$ {\boldsymbol{Ku}} = {\boldsymbol{F}} $对$ \rho \left( {{{\boldsymbol{y}}_m}} \right) $求导得$\dfrac{{{\text{d}}{\boldsymbol{u}}}}{{{\text{d}}\rho \left( {{{\boldsymbol{y}}_m}} \right)}}$的表达式为

已知结构拓扑构型与载荷无关, 因此有$\dfrac{{{\text{d}}{\boldsymbol{F}}}}{{{\text{d}}\rho \left( {{{\boldsymbol{y}}_m}} \right)}} = 0$. 将式(26)代入式(25)中, 可得

引入伴随向量${\boldsymbol{\gamma }}$, 式(27)表示为

依据式(16)的插值函数, 可得到

在单元层面上, 关心性能值$\sigma $对设计变量${\chi _k}$的灵敏度可表示为

4.3   非概率可靠性拓扑优化流程图

为便于理解式(19)的非概率可靠性拓扑优化模型, 图4给出了优化过程的流程图. 优化过程为嵌套优化, 其中右侧为内层优化, 目的是得到在几何不确定性下的关心性能值, 求解过程涉及极限状态函数$ G\left( {{{ {\textit{χ}}}},{\boldsymbol{\xi }}} \right) $对不确定系数 ${\xi _j}$的灵敏度分析(第4.1节内容); 左侧为外层优化, 即在给定体积约束下获得结构的最优拓扑, 求解过程涉及关心性能值$ G\left( {{{ {\textit{χ}}}},{{\boldsymbol{\xi }}^*}} \right) $对设计变量${\chi _k}$的灵敏度分析(第4.2节内容).

图  4  非概率可靠性拓扑优化流程图

Figure  4.  The flowchart of non-probabilistic reliability-based topology optimization

下载: 全尺寸图片 幻灯片

5.   数值算例

本节给出了2个数值算例来验证本文提出的考虑几何不确定性的非概率可靠性优化模型的有效性. 对于这2个算例, 结构的弹性模量和泊松比分别为${E_0} = 2.0 \times {10^5}{\text{ MPa}}$和$\upsilon = 0.3$. 在整个优化过程中, 结构始终处于线弹性阶段. 这2个算例都是Windows 10 操作系统性下基于MATLAB 软件计算的, 电脑的配置为AMD Ryzen Threadripper PRO 5965 WX 24-Cores 3.80 GHz, 128 GB RAM.

5.1   MBB梁结构的非概率可靠性优化

如图5所示为半个MBB梁结构, 结构的尺寸为$140{\text{ mm}} \times 70{\text{ mm}}$, 整个设计域离散为$140 \times 70\left( {9800} \right)$个平面应力单元. 结构右下角约束了Y方向的位移, 左边约束了$X$方向的位移. 集中载荷$F = 100{\text{ N}}$垂直作用于结构的左上角, 方向为$ - Y$方向. 不确定阈值场的范围为$ \eta (x)\in \left[0.25, 0.75\right] $, 在本算例中考虑了不确定阈值场2种不同的相关长度, 即$L = 80{\text{ mm}}$和$L = 30{\text{ mm}}$. 此外, 还讨论了2种不同的可靠性指标约束($\underline \beta = 1.0$和$\underline \beta = 1.5$)对拓扑结构的影响. 本算例给定的结构体积分数为$f = 50\% $.

图  5  MBB梁结构设计域

Figure  5.  Design domain for the MBB beam structure

下载: 全尺寸图片 幻灯片

作为对比, 首先给出了本算例的确定性优化结果, 即不确定阈值场$\eta ({\boldsymbol{x}}) \equiv 0.5$的情况(等价于相关长度$L \to + \infty $时的阈值场), 优化结果如图6所示. 考虑不确定阈值场的结构非概率可靠性拓扑优化结果如图7所示. 通过对比非概率可靠性优化结果(图7左列)与确定性优化结果(图6), 可以看出非概率可靠性优化的最优拓扑构型与确定性优化的不同, 并且非概率可靠性优化的最优解是通过增加更多的肋来保证结构在不确定阈值场下的可靠性.

图  6  MBB梁结构确定性拓扑优化结果(不确定阈值场$\eta ({\boldsymbol{x}}) \equiv 0.5$), 目标函数值$C = 34.22{\text{ N}} \cdot {\text{mm}}$

Figure  6.  Deterministic topology optimization solution for the MBB beam structure (uncertain threshold field $\eta ({\boldsymbol{x}}) \equiv 0.5$) and the objective function value $C = 34.22{\text{ N}} \cdot {\text{mm}}$

下载: 全尺寸图片 幻灯片

图  7  不同相关长度$L$及不同可靠性指标$\underline \beta $约束下的非概率可靠性拓扑优化结果(左列)和关心点处的不确定阈值场$\eta ({\boldsymbol{x}})$分布情况(右列)

Figure  7.  Non-probabilistic reliability-based topology optimization solution (left column) and the distribution of uncertain threshold fields at the concerned point (right column) with different correlation length $L$ and non-probability reliability index $\underline \beta $

下载: 全尺寸图片 幻灯片

此外, 图7(a)和图7(b)为考虑了不确定阈值场$\eta ({\boldsymbol{x}})$相同的相关长度($L = 80{\text{ mm}}$), 但优化模型(式(19))的可靠性指标约束值不同($\underline \beta = 1.0$和$\underline \beta = 1.5$). 对比这两个工况, 可以看出非概率可靠性拓扑优化结果存在一定差异(如图7(a)和图7(b)左列), 即结构中的“杆件”数量和位置存在明显不同. 另外, 当$\underline \beta $不同时, 尽管关心点处的不确定阈值场空间分布波动情况类似, 但不确定场$\eta ({\boldsymbol{x}})$的空间波动范围会随着可靠性指标约束值$\underline \beta $的增大而增大: (1) $\underline \beta = 1.0$时, 波动范围是$ \left[0.34, 0.75\right] $(如图7(a)右列); (2) $\underline \beta = 1.5$时, 波动范围是$ \left[0.19, 0.87\right] $(如图7(b)右列). 优化结果存在差异的原因在于, 随着可靠性指标约束值$\underline \beta $的增加, 进行结构的可靠性设计时考虑的几何不确定性范围更大, 关心点处不确定阈值场的波动范围也会更大(如图7(a)和图7(b)右列). 同时, 由于考虑了更多的不确定性, 那么设计的结构可靠性会更高, 最优的拓扑构型也会相应地发生变化(如图7(a)和图7(b)左列). 此外, 不确定阈值场$\eta ({\boldsymbol{x}})$的分布特点基本满足在结构左侧区域数值较大, 右侧区域数值较小的特点.

当非概率可靠性指标约束值为定值($\underline \beta = 1.0$)时, 不确定阈值场相关长度的变化也会对非概率可靠性优化的最优解产生影响(如图7(a)和图7(c)左列). 此外, 对比图7(a)和图7(c)右列图可以看出, 随着不确定阈值场相关长度L的减小, 不确定阈值场的空间波动程度也变得更为剧烈. 综上所述, 可靠性指标约束$\underline \beta $和不确定阈值场的相关长度L都会对结构的非概率可靠性优化结果产生影响.

不同情况下的考虑几何不确定性的非概率可靠性优化模型计算时间如表1所示. 其中, 内层优化和外层优化的计算时间均为迭代一次的平均时间, 内层优化的迭代步数为30步左右, 外层优化迭代为150步左右. 从表1中可以看出, 当不确定阈值场的相关长度较小($L = 30{\text{ mm}}$)时, 内层优化需要的时间更长(23.360 4 s). 这是因为当相关长度较小时, 则需要更多的截断项才能保证要求的不确定场截断精度(式(5)), 那么式(19)中内层优化中的变量数(${\xi _j}$)会更多, 计算时间也会相应地增加.

表  1  不同情况下非概率可靠性优化模型计算时间

Table  1.  The computing time of non-probability reliability-based topology for different cases

Correlationlength/mm	Reliability-based index	Average calculation timeT/s
		Inner-loop optimization	Outer-loop optimization
$L \to + \infty $ (deterministic)	—	—	0.720 5
$L = 80$	$\underline \beta = 1.0$	14.871 0	0.7207
$L = 80$	$\underline \beta = 1.5$	15.186 2	0.712 6
$L = 30$	$\underline \beta = 1.0$	23.360 4	0.723 1

下载: 导出CSV 

| 显示表格

为了说明优化结果的正确性, 我们将通过差分灵敏度分析法来验证第4节的解析灵敏度分析的正确性, 差分法采用的差分步长为0.001. 以$L = 80{\text{ mm}}$, $\underline \beta = 1.0$的工况为例, 内层优化共包含25个不确定系数${\xi _j}\left( {j = 1,2,\cdots,25} \right)$, 随机选取4个不确定系数 ($ {\xi _2},{\xi _5},{\xi _{10}},{\xi _{15}} $) 进行灵敏度对比, 对比结果如表2所示. 外层优化同样随机选取4个设计变量(${\chi _2},{\chi _{10}}, {\chi _{30}},{\chi _{70}}$), 对比结果如表3所示.

表  2  内层优化中不确定系数的灵敏度

Table  2.  Sensitivity of uncertainty coefficients in inner-loop optimization

j-th uncertainty coefficient	$ {\xi _2} $	$ {\xi _5} $	$ {\xi _{10}} $	$ {\xi _{15}} $
Difference method	25.223	7.009	1.772	0.251
Analytical method	25.211	7.022	1.787	0.257
Relative error/%	0.047	0.184	0.873	0.884

下载: 导出CSV 

| 显示表格

表  3  外层优化中设计变量的灵敏度

Table  3.  Sensitivity of design variables in outer-loop optimization

k-th design variable	$ {\chi _2} $	$ {\chi _{10}} $	$ {\chi _{30}} $	$ {\chi _{70}} $
Difference method	25.715	−45.458	24.742	−5.482
Analytical method	25.713	−45.459	24.757	−5.467
Relative error/%	0.008 4	0.002 5	0.061 5	0.257

下载: 导出CSV 

| 显示表格

从表2和表3的对比结果中可以看出, 第4节的解析灵敏度和差分法的灵敏度计算结果非常相近, 其相对误差均在 1%以内. 这证明了本文优化模型及解析灵敏度推导的正确性.

基于本文提出的考虑结构几何不确定性的非概率可靠性拓扑优化模型, 本算例中不同工况下的内层优化和外层的优化迭代历史分别如图8, 图9所示. 其中, 图8表示第一次的内层优化. 从图8中可以看出, 在整个内层优化的过程中, 目标函数都能够快速收敛到某一特定值(图8实线所示), 并且结构能够始终满足指定的可靠性指标约束(图8虚线所示). 此外, 图9的外层优化迭代历史也体现了优化模型稳定的收敛性. 因此, 这也进一步证明了采用梯度优化算法来求解本文提出的非概率可靠性优化模型是合理有效的.

图  8  内层优化迭代历史(实线表示目标函数, 虚线表示可靠性指标约束)

Figure  8.  Iteration history of inner-loop optimization (solid lines represent the objective function and dashed lines represent reliability index constraints)

下载: 全尺寸图片 幻灯片

图  9  确定性和非概率可靠性拓扑优化的迭代历史

Figure  9.  Iteration histories of the deterministic and non-probability reliability-based topology optimization

下载: 全尺寸图片 幻灯片

为了进一步验证非概率可靠性拓扑优化结果的合理性, 我们将非概率可靠性拓扑优化关心点 $ {{\boldsymbol{\xi }}^*} $(式19)处的几何不确定性(即表示为图7右列的不确定阈值场)施加到确定性拓扑优化结构(图6所示)和可靠性拓扑优化结构(图7左列)上, 来进行结构的柔顺性对比, 对比结果如表4所示. 以$L = 80{\text{ mm}}$; $\underline \beta = 1.0$的非概率可靠性拓扑优化结构为例(如表4第2, 3列所示), 当确定性优化结构中存在几何不确定性时(即关心点 $ {{\boldsymbol{\xi }}^*} $处的不确定阈值场), 结构柔顺性值为$C = 34.95{\text{ N}} \cdot {\text{mm}}$, 大于可靠性拓扑优化结构在该关心点$ {{\boldsymbol{\xi }}^*} $处的柔顺性值$C = 33.66{\text{ N}} \cdot {\text{mm}}$. 此外, 其他条件下的非概率可靠性拓扑优化结构同样满足这种情况: (1) $L = 80{\text{ mm}}$, $\underline \beta = 1.5$时, 如表4第4, 5列所示, $C = 35.57{\text{ N}} \cdot {\text{mm}}$和$C = 34.94{\text{ N}} \cdot {\text{mm}}$; (2) $L = 30{\text{ mm}}$, $\underline \beta = 1.0$时, 如表4第6,7列所示, $C = 38.78{\text{ N}} \cdot {\text{mm}}$和$C = 34.87{\text{ N}} \cdot {\text{mm}}$. 因此, 通过上述对比可以看出, 当存在几何不确定性时, 考虑几何不确定性的非概率可靠性优化结构具有更好的抵抗几何不确定性的能力, 即有更小的结构柔顺性.

表  4  可靠性拓扑优化结构和确定性优化结构在关心点${{{\boldsymbol{\xi}} }^{*}}$处的柔顺性对比

Table  4.  Comparison of compliance between reliability-based and deterministic topology optimization structures at concerned point ${{{\boldsymbol{\xi}} }^{*}}$

Case	$L=80\text{ mm, }\underline{\beta }=1.0$		$L=80\text{ mm, }\underline{\beta }=1.5$		$L=30\text{ mm, }\underline{\beta }=1.0$
Topological structure
Uncertain threshold field in ${{\boldsymbol{\xi }}^{*}} $
Compliance C/(N·mm)	34.95	33.66	35.57	34.94	38.78	34.87

下载: 导出CSV 

| 显示表格

5.2   悬臂梁结构的非概率可靠性优化

在本算例中, 以悬臂梁结构为例, 如图10所示, 悬臂梁结构设计域的尺寸为$112{\text{ mm}} \times 70{\text{ mm}}$, 离散为$7840\;\left( {112 \times 70} \right)$个单元. 悬臂梁结构的左侧完全约束, 集中力$F = 100{\text{ N}}$作用于结构的右下角, 方向为$ - Y$方向. 不确定阈值场的变化范围为$ \eta (x)\in \left[0.25, 0.75\right] $. 结构的体积分数约束同样设定为$f = 50\% $. 本算例考虑了一种不确定阈值场的相关长度$L = 30{\text{ mm}}$和一种可靠性指标约束值$\underline \beta = 1.0$来进行考虑几何不确定性的结构非概率可靠性优化设计.

图  10  悬臂梁结构设计域

Figure  10.  Design domain for the cantilever structure

下载: 全尺寸图片 幻灯片

作为对比, 本算例中同样进行了悬臂梁结构的确定性拓扑优化设计, 在整个优化过程中阈值场在设计域内各点处的阈值均为恒定值, 即$\eta ({\boldsymbol{x}}) \equiv 0.5$, 悬臂梁结构的确定性拓扑优化结果如图11所示. 考虑不确定阈值场的结构非概率可靠性拓扑优化结果及关心点对应的不确定阈值场分布情况如图12所示($L = 40{\text{ mm}},$ $\underline \beta = 1.0$). 对比2个优化结果可以看出, 在材料体积一定的情况下, 可靠性优化设计的结构产生了更多的肋来抵抗结构中的几何不确定性, 以提高结构的可靠性. 本算例进一步证明了本文提出的考虑几何不确定性的非概率可靠性优化模型的有效性.

图  11  悬臂梁结构确定性拓扑优化结果(不确定阈值场$\eta ({\boldsymbol{x}}) \equiv 0.5$), 目标函数值$C = 19.62{\text{ N}} \cdot {\text{mm}}$

Figure  11.  Deterministic topology optimization solution for the cantilever structure (uncertain threshold field $\eta ({\boldsymbol{x}}) \equiv 0.5$) and the objective function value $C = 19.62{\text{ N}} \cdot {\text{mm}}$

下载: 全尺寸图片 幻灯片

图  12  相关长度$L = 40{\text{ mm}}$及可靠性指标约束值$\underline \beta = 1.0$时的悬臂梁结构非概率可靠性拓扑优化结果和关心点处的阈值场分布, 目标函数值$C = 19.79{\text{ N}} \cdot {\text{mm}}$

Figure  12.  Non-probabilistic reliability-based topology optimization solution and the distribution of uncertain threshold fields at the concerned point with correlation length $L = 40{\text{ mm}}$ and non-probability reliability index$\underline \beta = 1.0$ and the objective function value $C = 19.79{\text{ N}} \cdot {\text{mm}}$

下载: 全尺寸图片 幻灯片

6.   结 论

本文提出了一种考虑结构几何不确定性的非概率可靠性优化模型, 其中结构的几何不确定性通过有界不确定阈值场函数来表示. 基于材料场级数展开, 不确定阈值场可通过一组不确定系数来描述. 优化模型为嵌套优化问题, 内层是进行结构的可靠性评估, 外层为确定结构的最优布局. 本文采用了基于梯度的优化算法来求解该非概率可靠性优化模型. 数值算例表明, 采用本文提出的非概率可靠性优化方法, 可以得到分布更为合理的拓扑构型来提高结构在考虑几何不确定性下的可靠性.