结构拓扑优化本着"物尽其用"的宗旨,以期最大限度地挖掘材料潜力,设计出满足苛刻要求的创新结构构型,被公认为结构设计领域最活跃、最具发展前景和挑战性的研究课题之一.

目前,拓扑优化理论的发展形成了各具特色的设计方法. 如变密度法^[1-2]、渐进结构法^[3-4]直接将网格单元作为拓扑优化的基本设计要素,通过确定单元的"有"、"无"得到材料的最优分布.由于物理概念清晰简洁且易于编程实现,此类方法已被广泛应用于高端机械产品设计和航空航天飞行器研制中,但所优化得到的结果边界一般呈非光滑的锯齿状.为解决这一问题,借助高维水平集函数隐式设计结构拓扑演化轮廓的水平集法^[5-6]被提出并迅速成长起来.然而,传统水平集法在优化过程中需不断进行重新初始化和速度场扩展等操作,还需限制步长以满足Courant-Friedrichs-Lewy条件,求解效率较低.

为消除传统变密度法的网格依赖性和棋盘格等数值不稳定现象,Guest等^[7]发展出节点密度法,并借助Heaviside阶跃函数得到清晰的0-1材料分布;为克服传统水平集法的数值计算困难以提升其优化效率和通用性,Wang等^[8]发明出参数化水平集法,无需直接求解Hamilton-Jacobi偏微分方程且可配合使用成熟的优化算法.节点密度法^[9-10]和参数化水平集法^[11-13]虽属于不同类型的拓扑优化方法,却有着一定的相似之处：结构模型由高一维的隐函数描述,此隐函数又由高阶连续的插值函数构造而成,它在节点密度法中就是具有^[0,1]界限的密度场,在参数化水平集法中就是水平集函数.Sigmund等^[14]认为近期不同类型拓扑优化方法的发展使其之间的差异越来越小.然而,该发展趋势下形成的新方法在减少设计变量方面的效果甚微.节点密度法和参数化水平集法所采用的插值函数的数量需能够保证密度场/水平集函数具有一定的光滑连续性和足够的拓扑描述能力,无疑在二维问题中就要使用上千甚至上万个设计变量.

Eschenauer等^[15]早期提出的泡泡法能够在保证拓扑优化设计精度的条件下有效降低变量数目：采用参数化B样条曲线精确表达孔洞边界,使用沿结构边界生成的有限元网格保证分析精度,通过逐步增加新孔丰富设计变量数目.但是,所采用的有限元拉格朗日网格给泡泡法的实施带来了极大不便,每次优化迭代后均需调整甚至重新划分网格.而近年来兴起的特征驱动类拓扑优化方法^[16-21]则采用了固定网格分析方式,并将具有一定造型能力且参数少的实体/空洞特征作为设计要素,通过特征的移动和变形,达到少变量高精度拓扑优化设计的目的.然而,这类方法需要在初始设计中布局数量充足的特征要素,难免会造成设计变量的浪费,并且存在初始布局依赖性问题(见图1).

本研究在借鉴特征驱动类方法的优点的基础上,对传统泡泡法进行改进,提出一种基于固定网格和拓扑导数的自适应泡泡拓扑优化方法.该方法采用有限胞元方法^[22-24]在固定网格下计算结构力学响应,避免了传统泡泡法中繁琐的网格更新和重划分操作;优化迭代时通过综合考虑拓扑导数信息^[25-27]和创新性地设置孔洞影响区域,不仅能自适应地在合理位置逐步引入新的孔洞,还能克服特征驱动类方法中存在的初始布局依赖性,从而有效提升了数值计算稳定性;同时采用参数少且变形能力强的光滑变形隐式曲线^[28-30]描述孔洞边界,有效减少了设计变量.本文首先介绍了自适应泡泡法的3个主要模块(固定网格分析、孔洞自适应引入和孔洞隐式描述),然后建立优化模型并给出解析灵敏度计算公式,最后通过经典悬臂梁和简支梁算例验证了方法的有 效性.

有限胞元方法(finite cellmethod,FCM)^[22]和扩展有限元方法(eXtended finite elementmethod,XFEM)^[31]是两种常见的高精度固定网格分析方法.前者在本质上是一种采用高阶形函数逼近待求物理场的虚拟区域法;后者的核心思想是通过改进形函数处理单元内材料属性存在的不连续性.由于本文工作仅涉及单材料结构的拓扑优化设计,而且考虑到高阶形函数可以带来较高的分析精度,因此选用FCM计算结构力学响应.

FCM的计算原理如图2所示,为求解结构区域$\varOmega$上的力学问题,先将$\varOmega$嵌入到规则的计算域$D$中,然后使用规则网格将$D$离散,并基于加权余量法或变分原理建立起离散后的平衡方程组,最后通过求解此方程组以得到待求物理场.对于线弹性力学问题,相应的线性方程组为\begin{equation} \label{eq1} {{\boldsymbol {KU}}} = {{\boldsymbol F}} (1) \end{equation} 式中,$\textbf{U}$为待求位移向量,$\textbf{K}$和$\textbf{F}$分别是刚度矩阵和载荷向量,计算公式如下

$${{\boldsymbol K}} = \int_D {{{\boldsymbol B}}^{T}\beta {{\boldsymbol {DB}}}} \mbox{d}\varOmega (2)\\ {{\boldsymbol F}} = \int_D {{\boldsymbol N}^{T}\beta {{\boldsymbol f}}\mbox{d}\varOmega } + \int_{\varGamma _{N} } {{\boldsymbol N}^{T}{{{\boldsymbol t}}}\mbox{d}\varGamma } (3)$$

式中,$\textbf{B}$为应变矩阵,$\textbf{D}$为弹性矩阵, $\textbf{N}$为形函数矩阵(本工作采用双二次B样条基函数作为形函数^[24]),$\textbf{f}$为体积力向量,$\textbf{t}$为边界力向量,$\beta $为定义如下的标量因子 \begin{equation} \label{eq4} \beta = \left\{ {{\begin{array}{ll} 1, &\mbox{in }\varOmega \\ 0, & {\mbox{in }D/ \varOmega } \\ \end{array} }} \right.(4) \end{equation}

由式(1)~式(4)可知,FCM与传统有限元方法的求解方式相同,但需将由计算域$D$离散而成的规则胞元区分为三种类型分别处理：虚拟胞元(fictitiouscell)、物理胞元(physical cell)和边界胞元(boundary cell),见图2.其中,虚拟胞元内$\beta$为0,则无需积分计算;物理胞元内$\beta$为1,则按照有限元的积分方式进行计算;边界胞元内$\beta$取值不定,则在积分时一般采用四叉树(quadtree)细化技术生成沿结构边界高度聚集的高斯积分点(参见图2),从而实现对胞元内实体材料部分的高精度积分计算.

当结构边界发生变化时,FCM的网格无需更新或重划分,这是其与传统有限元方法的主要区别.

拓扑导数的概念源于泡泡法^[15],反映了无限小单一区域的扰动(引入孔洞、夹杂、源项甚至裂纹等)对给定函数的影响.Sokolowski等^[25]给出了一般性拓扑导数推导过程,Novotny等^[26]证明了拓扑导数是形状灵敏度的极限,使得拓扑导数的推导可以采用理论成熟的形状灵敏度分析方法.Novotny和Sokolowski^[27]详细介绍了拓扑导数的概念、基本理论以及在各类问题中的应用.

对于平面应力问题,若不考虑体积力等设计相关性载荷,则结构柔顺度($J$ = $F^{ T}$$\textbf{U}$)对结构内任意一点($\textbf{x}$$_{\varOmega } \in $$\varOmega $)处的圆孔引入的拓扑导数为^[27]\begin{equation} \label{eq5} \begin{array}{l} D_{T} ({{\boldsymbol x}}_\varOmega ) = \dfrac{4}{1 + \nu }{{\boldsymbol S}}({ {\boldsymbol x}}_\varOmega ):{{\boldsymbol E}}({{\boldsymbol x}}_\varOmega )-\\ \qquad \dfrac{1-3\nu }{1-\nu ^2}\mbox{tr}\left( {{{\boldsymbol S}}({{\boldsymbol x}}_\varOmega )} \right)\mbox{tr}\left( {{{\boldsymbol E}}({{\boldsymbol x}}_\varOmega )} \right) \\ \end{array} (5) \end{equation} 式中,$\textbf{S}$为应力张量,$\textbf{E}$为应变张量,$\nu $为泊松比,tr($\cdot )$表示张量的迹. 二阶张量$\textbf{S}$($\textbf{x}$$_{\varOmega })$和$\textbf{E}$($\textbf{x}$$_{\varOmega })$在笛卡尔坐标系下对应的矩阵为

$${{\boldsymbol S}}({{\boldsymbol x}}_\varOmega ) = \left( {\sigma _{ij} ({{\boldsymbol x}}_\varOmega )} \right) = {\left[ {{\begin{array}{*{20}c} {\sigma _x ({{\boldsymbol x}}_\varOmega )} & {\tau _{xy} ({{\boldsymbol x}}_\varOmega )} \\ {\tau _{yx} ({{\boldsymbol x}}_\varOmega )} & {\sigma _y ({{\boldsymbol x}}_\varOmega )} \end{array} }} \right]} (6)\\ {{\boldsymbol E}}({{\boldsymbol x}}_\varOmega ) = \left( {\varepsilon _{ij} ({{\boldsymbol x}}_\varOmega )} \right) = {\left[ {{\begin{array}{*{20}c} {\varepsilon _x ({{\boldsymbol x}}_\varOmega )} & {\dfrac{\gamma _{xy} ({ {\boldsymbol x}}_\varOmega )}{2}} \\ {\dfrac{\gamma _{yx} ({{\boldsymbol x}}_\varOmega )}{2}} & {\varepsilon _y ({{\boldsymbol x}}_\varOmega )} \end{array} }} \right] }(7) $$

式中,剪应力$\tau _{xy}(\textbf{x}_{\varOmega }) = {\tau }_{yx}(\textbf{x}_{\varOmega })$,剪应变$\gamma _{xy}(\textbf{x}_{\varOmega }) =\gamma_{yx}(\textbf{x}_{\varOmega })$,点$\textbf{x}_{\varOmega }$处的应力分量和应变分量计算如下

$$ {{\boldsymbol \sigma }}({{\boldsymbol x}}_\varOmega ) = (\sigma _x ({{\boldsymbol x}}_\varOmega ),\sigma _y ({{\boldsymbol x}}_\varOmega ),\tau _{xy} ({{\boldsymbol x}}_\varOmega ))^{T} = {{\boldsymbol D\varepsilon }}({{\boldsymbol x}}_\varOmega ) =\\ \qquad {{\boldsymbol D}}(\varepsilon _x ({{\boldsymbol x}}_\varOmega ),\varepsilon _y ({ {\boldsymbol x}}_\varOmega ),\gamma _{xy} ({{\boldsymbol x}}_\varOmega ))^{T} = {{\boldsymbol {DB}}}({{\boldsymbol x}}_\varOmega ){{\boldsymbol U}} (8) $$

式中参数定义同前.

由于孔洞的引入会削弱结构的刚度,也即使结构的柔顺度$J$增大, 因此由式(5)得到的拓扑导数满足$D_{ T}$($\textbf{x}$$_{\varOmega }) \ge $ 0. 拓扑导数值$D_{ T}$($\textbf{x}$$_{\varOmega })$反映了在点$\textbf{x}$$_{\varOmega }$处引入圆孔对结构柔顺度$J$的影响程度,若$D_{ T}$($\textbf{x}$$_{\varOmega })$等于或者接近于0,则说明在点$\textbf{x}$$_{\varOmega }$处开孔后结构柔顺度不增加或者增加量很小. 由此,在以柔顺度最小化为设计目标的拓扑优化设计中,开孔位置要优先选在$D_{ T}$($\textbf{x}$$_{\varOmega })$最小的点处.

结构拓扑优化是一个迭代过程,结构的物理场在每步优化迭代之后就会发生变化,这就需要依据变化后的物理场调整结构的形状和拓扑. 若每步因孔洞数量引入过多,使得形状/拓扑改变量过大,则极易导致优化设计收敛于局部最优解或者无法收敛. 此外,在引入孔洞之后还要对其进行形状优化设计,如果两个孔之间相距较近,则它们在形状优化过程中就会相互重叠甚至吞并,从而造成设计变量的 浪费.

为解决上述两个问题,本文定义了拓扑导数阈值,以限制每步优化迭代中潜在开孔位置的数量;并提出了"孔洞影响区域"概念,要求结构中已存在的孔洞都有各自独立的影响区域,在这些影响区域内不能引入新的孔洞.

本文所研究的拓扑优化设计问题是如何在矩形计算域$D$(即初始未引入孔洞的结构区域$\varOmega $)内优化材料分布.图3展示了单步优化迭代的孔洞自适应引入机制,具体步骤如下：

(1)在计算域$D$上设置均匀分布的采样点,并根据式(5)计算结构区域$\varOmega $内采样点的拓扑导数$D_{T}(\textbf{x}_{\varOmega })$. 本工作直接将有限胞元固定网格节点作为采样点,$\varOmega $外的采样点在图3中用黑色圆点标记.

(2)为限制孔洞引入数量,对采样点的拓扑导数值从小到大进行排序,并将第$m$个拓扑导数值作为阈值$\overline {D_{T} } $,要求孔洞只能在满足$D_{ T}(\textbf{x}_{\varOmega }) \le \overline {D_{T} } $的$m$个采样点处引入.图3已对$D_{ T}(\textbf{x}_{\varOmega })$最小的$m=10$个点从小到大标记了顺序,本工作使用的$m$一般选为采样点(也即固定网格节点)总数的3%~5%左右.

(3)设置已有孔洞的影响区域$A_{i}$,其可通过将孔洞边界整体向外偏移一定距离来实现,如图3中颜色较深的蓝色区域. 将结构中累计引入的孔洞数量记为$num$,则所有孔洞的影响区域为集合$ {U}_{i = 1}^{num} A_i $.

(4)为避免在影响区域$ {U}_{i = 1}^{num} A_i$内引入孔洞,将拓扑导数值不大于阈值$\overline {D_{T} }$的$m$个采样点,按照其是否属于$\varOmega $和$ {U} _{i= 1}^{num} A_i $进行分类. 对属于$\varOmega /{U} _{i = 1}^{num} A_i $的这类采样点(图3中的红色圆点),在其中$D_{T}(\textbf{x}_{\varOmega })$最小的点处引入圆孔.

(5)重复第3, 4步的操作,直到$\varOmega $内满足$D_{T}({{\boldsymbol x}}_\varOmega ) \le \overline {D_{T} }$的采样点均属于$ {U} _{i = 1}^{num} A_i $.

本文用于描述已引入孔洞的光滑变形隐式曲线,本质上是对超二次曲线的一种新型扩展,并因此继承了其性质与优点.超二次曲线能够用少量的参数表达形状各异的平面物体,而且在定义时既可采用隐式形式\begin{equation} \label{eq9} \left| {\frac{x}{a}} \right|^{2 / \varepsilon } + \left| {\frac{y}{b}} \right|^{2 / \varepsilon }-1 = 0 (9) \end{equation} 又可采用参数形式\begin{equation} \label{eq10} \left. {{\begin{array}{l} {x = a{sgn}\left( {\cos t} \right)\left| {\cos t} \right|^\varepsilon } \\ {y = b{sgn}\left( {\sin t} \right)\left| {\sin t} \right|^\varepsilon } \\ \end{array} }} \right\} (10) \end{equation}式中,$a$和$b$为直角坐标系中2个主轴方向的半轴长;$\varepsilon$为形状指数,$a,b$和$\varepsilon $恒为正,sgn($\cdot)$为符号函数,$t$的变化范围为$[0, 2\pi)$.

为进一步提高超二次曲线的表达能力,使其能够描述非对称几何模型,周林等^[32-33]将超二次曲线的形状指数$\varepsilon$扩展成极角的函数,形成了扩展超二次曲线.本工作借鉴了此思想,将超二次曲线的半轴长$a$和$b$扩展成极角的函数,定义出表达能力(也即变形能力)更强的光滑变形隐式曲线^[30],用隐式形式表示如下\begin{equation} \label{eq11} \left| {\frac{x}{R\left( \theta \right)}} \right|^{2 / \varepsilon } + \left| {\frac{y}{R\left( \theta \right)}} \right|^{2 / \varepsilon }-1 = 0 (11) \end{equation}其中,$\theta $为极坐标系中的极角,半轴长函数$R$($\theta $)可借助B样条基函数定义为 \begin{equation} \label{eq12} R(\theta ) = \sum\limits_{j = 1}^n {N_{j,p} \left( \xi \right)P_j } = \sum\limits_{j = 1}^n {N_{j,p} \left( {\frac{\theta + \pi / 2}{2\pi }} \right)P_j } (12) \end{equation} 式中,$N_{j, p}(\xi )$为B样条基函数;$n$是B样条基函数的个数;$p$是B样条基函数的次数(本工作中的$p$取为2);$\xi $为参数区域坐标,$P_{j}$为控制参数,$P_{j}$恒为正且满足$P_{1}$ = $ P_{n}$ = ($P_{2}+P_{n-1})$/2,极角$\theta $可表示为坐标($x$, $y)$的函数 \begin{equation} \label{eq13} \theta = \arctan \frac{y}{x} + {sgn}x\left( {{ sgn}x-1} \right)\frac{\pi }{2} (13) \end{equation} 由上式可知$\theta $的值域为$[-\pi /2, 3\pi /2)$,则式(12)中参数区域坐标$\xi $的取值范围为$[0, 1)$.

若将式(11)中的形状指数$\varepsilon$设为常数1,则光滑变形隐式曲线退化为文献^[28]中定义的闭合B样条曲线(closedB-splines),此时$R(\theta $)即为极角$\theta$对应的极径;若$\varepsilon$不等于1,则$R(\theta )$仍与极角$\theta$处的极径成正相关关系.图4(a)展示了光滑变形隐式曲线的构建原理,由此图和式(11)~式(13)可知,调整控制参数$P_{i}$会使$R(\theta)$发生变化,而极坐标系中隐式曲线上点的极径$\rho(\theta )$与$R(\theta)$正相关,则隐式曲线的形状也因此发生改变.

本工作使用光滑变形隐式曲线描述结构中已引入的孔洞,参数$\varepsilon$和$P_{j}$对孔洞形状都有直接影响,所需变量相对很少.为便于固定网格分析和孔洞合并操作等,本工作在优化过程中使用孔洞的隐函数,根据式(11)可写为\begin{equation} \label{eq14} \varphi = \left( {\left| {\frac{x}{R\left( \theta \right)}} \right|^{2 / \varepsilon } + \left| {\frac{y}{R\left( \theta \right)}} \right|^{2 / \varepsilon }} \right)-1 (14) \end{equation}其中,$\varphi =0$表示孔洞边界;$\varphi <0$则表示孔洞.若将上式中$R(\theta)$的参数$P_{j}$统一增加数值$r$,则相应的孔洞边界就会向外偏移.在不考虑计算域的情况下,可将此偏移区域设定为该孔洞的影响区域,如图4(b)所示.

鉴于在拓扑优化设计中要引入一定数量的孔洞,而且每个孔洞还要能够自由地移动,这就需要按照孔洞的引入顺序$i$将其隐函数具体写为\begin{equation} \label{eq15} \varphi _i = \left( {\left| {\frac{x-x_i }{R_i \left( {\theta _i } \right)}} \right|^{2 / \varepsilon _i } + \left| {\frac{y-y_i }{R_i \left( {\theta _i } \right)}} \right|^{2 / \varepsilon _i }} \right)-1 (15) \end{equation}其中

$$ R_i (\theta _i ) = \sum\limits_{j = 1}^n {N_{j,p} \left( {\xi _i } \right)P_{i,j} } = \sum\limits_{j = 1}^n {N_{j,p} \left( {\frac{\theta _i + \pi / 2}{{2\pi }}} \right)P_{i,j} } (16)\\ \theta _i = \arctan \frac{y-y_i }{x-x_i } + \frac{\pi }{2}{sgn}(x-x_i )\left[ {{sgn}(x-x_i )-1} \right] \quad(17) $$

式中,($x_{i}$, $y_{i})^{T}$是构建第$i$个孔洞所用的极坐标系的极点(参见图4(a)中的$O$点),$R_{i}(\theta_{i}),\varepsilon _{i},\theta _{i},P_{i, j}$和$\xi_{i}$分别是第$i$个孔洞的半轴长函数、形状指数、极角、控制参数和参数区域坐标.本工作中每个孔洞的B样条基函数$N_{j,p}$及其个数$n$都一样,因此在上式中没有添加标注$i$.如果已知计算域$D$的隐函数$\varPhi_{D}(\varPhi _{D}>0$时表示$D$,$\varPhi_{D}$可基于R函数理论精确得到^[34]),则结构区域$\varOmega$的隐函数为 \begin{equation} \label{eq18} \varPhi _\varOmega = \min \left\{ {\varPhi _D ,\left\{ {\varphi _i } \right\}_{i = 1}^{num} } \right\} (18) \end{equation}式中,$num$为已引入的孔洞数量,在优化迭代中会逐渐增加,$\varPhi_{\varOmega}$满足 \begin{equation} \label{eq19} \left. {{\begin{array}{l} {\varPhi _\varOmega ({{\boldsymbol x}}) = 0,} \\ {\varPhi _\varOmega ({{\boldsymbol x}}) > 0,} \\ {\varPhi _\varOmega ({{\boldsymbol x}}) < 0,} \\ \end{array} }\mbox{ }{\begin{array}{l} {\forall {{\boldsymbol x}} \in \partial \varOmega } \\ {\forall {{\boldsymbol x}} \in \varOmega / \partial \varOmega } \\ {\forall {{\boldsymbol x}} \in D / \varOmega } \\ \end{array} }} \right\} (19) \end{equation} 由上式可知,根据$\varPhi_{\varOmega }$的正负性可以直接判断出点$\textbf{x}$是否位于结构区域$\varOmega $之内,从而就能很好地与有限胞元方法进行结合以实施胞元类型识别、四叉树细化和高斯积分点区分等步骤. 而且结合Heaviside函数$H(\cdot)$还可将式(4)定义的$\beta $用$\varPhi $$_{\varOmega }$表示为 \begin{equation} \label{eq20} \beta \mbox{ = }H\left( {\varPhi _\varOmega \left( { {\boldsymbol x}} \right)} \right) = \left\{ {{\begin{array}{*{20}c} {1,} \\ {0,} \\ \end{array} }\mbox{ }{\begin{array}{*{20}c} {\varPhi _\varOmega \left( {{\boldsymbol x}} \right) \ge 0} \\ {\varPhi _\varOmega \left( {{\boldsymbol x}} \right) < 0} \\ \end{array} }} \right. (20) \end{equation}

按照图4(b)的方法可将第$i$个孔洞的边界向外偏移,若将增大后的孔洞的隐函数记为$\varphi_i^r $,则孔洞影响区域的集合$ \cup _{i = 1}^{num} A_i $(参见图3)的隐函数为\begin{equation}\label{eq21} \varPhi _A = \min \left\{ {\varPhi _\varOmega ,-\min \left\{ {\varphi _i^r } \right\}_{i = 1}^{num} } \right\}(21)\end{equation}根据{$\varPhi$}$_{A}$的正负性可以直接判断计算拓扑导数的采样点是否位于孔洞影响区域之内.

本工作旨在有限材料用量下实现结构刚度的最大(柔顺度最小)设计,这种典型拓扑优化问题的数学模型为 \begin{equation} \label{eq22}\left. \begin{array}{l} \mbox{find }\left\{ {x_i ,y_i ,\varepsilon _i ,\left\{ {P_{i,j} } \right\}_{j = 2}^{n-1} } \right\}_{i = 1}^{num} \\ \mbox{min }J = {\boldsymbol F}^{T}{{\boldsymbol U}} \\ \mbox{S. t. }{{\boldsymbol {KU}}} = {{\boldsymbol F}} \\ \mbox{ }V = \int_D {H\left( {\varPhi _\varOmega } \right)\mbox{d}\varOmega } \le V_{\lim } \\ \mbox{ }\varepsilon _i > 0,\mbox{ }P_{i,j} > 0\mbox{ }i = 1,2, \cdots ,{num}\mbox{ },j = 2,3, \cdots ,n-1 \\ \end{array} \right\} (22) \end{equation} 式中,$V$和$V_{lim}$分别为结构体积及其上限,其余参数定义同前. 每个孔洞有$n+1$个设计变量,现将设计变量统一记为 \begin{equation} \label{eq23} {{\boldsymbol v}} = \{v_k \}_{k = 1}^{{num}\times (n + 1)} = \mbox{ }\left\{ {x_i ,y_i ,\varepsilon _i ,\left\{ {P_{i,j} } \right\}_{j = 2}^{n-1} } \right\}_{i = 1}^{num} (23) \end{equation} 则柔顺度$J$对变量$v_{k}$的灵敏度为 \begin{equation} \label{eq24} \frac{\partial J}{\partial v_k } = \left( {\frac{\partial {{\boldsymbol F}}}{\partial v_k }} \right)^{ T}{{\boldsymbol U}} + {{\boldsymbol F}}^{T}\frac{\partial {{\boldsymbol U}}}{\partial v_k } (24) \end{equation} 如果$\textbf{F}$的大小、方向和施加位置在优化过程中不变,则式(24)的右端第一项为0,而且由式(1)可得 \begin{equation} \label{eq25} \frac{\partial {{\boldsymbol U}}}{\partial v_k } =-{{\boldsymbol K}}^{-1}\frac{\partial {{\boldsymbol K}}}{\partial v_k }{{\boldsymbol U}} (25) \end{equation} 于是式(24)可进一步简化为 \begin{equation} \label{eq26} \frac{\partial J}{\partial v_k } =-{\boldsymbol F}^{T}{{\boldsymbol K}}^{-1}\frac{\partial {{\boldsymbol K}}}{\partial v_k }{{\boldsymbol U}} =-{{\boldsymbol U}}^{ T}\frac{\partial {{\boldsymbol K}}}{\partial v_k }{ {\boldsymbol U}} (26) \end{equation} 其中,刚度矩阵灵敏度可根据式(2)和式(20)推导为

$$ \begin{aligned} &\dfrac{\partial {{\boldsymbol K}}}{\partial v_k } = \int_D {{{\boldsymbol B}}^{T}{{\boldsymbol {DB}}}\dfrac{\partial H\left( {\varPhi _\varOmega } \right)}{\partial v_k }\mbox{d}\varOmega } =\\ & \qquad \int_D {{{\boldsymbol B}}^{T}{{\boldsymbol {DB}}}\dfrac{\mbox{d}H\left( {\varPhi _\varOmega } \right)}{\mbox{d}\varPhi _\varOmega }\dfrac{\partial \varPhi _\varOmega }{\partial v_k }\mbox{d}\varOmega } =\\ &\qquad \int_D {{{\boldsymbol B}}^{T}{{\boldsymbol {DB}}}\dfrac{\partial \varPhi _\varOmega }{\partial v_k }\dfrac{1}{\left\| {\nabla \varPhi _\varOmega } \right\|}\left( {\dfrac{\mbox{d}H\left( {\varPhi _\varOmega } \right)}{\mbox{d}\varPhi _\varOmega }\left\| {\nabla \varPhi _\varOmega } \right\|} \right)\mbox{d}\varOmega } =\\ & \qquad \int_{\partial \varOmega } {{{\boldsymbol B}}^{ T}{{\boldsymbol {DB}}}\dfrac{\partial \varPhi _\varOmega }{\partial v_k }\dfrac{1}{\left\| {\nabla \varPhi _\varOmega } \right\|}\hat {\delta }\mbox{d}\varOmega } =\\ &\qquad \int_{\partial \varOmega } {{{\boldsymbol B}}^{ T}{{\boldsymbol {DB}}}\dfrac{\partial \varPhi _\varOmega }{\partial v_k }\dfrac{1}{\left\| {\nabla \varPhi _\varOmega } \right\|}\mbox{d}\varGamma } (27) \end{aligned} $$

式中,$\hat {\delta }$为Heaviside函数的法向导数,称为Dirac delta函数,它能将区域积分转换为边界积分.体积$V$对变量$v_{k}$的灵敏度可参照上式进行推导,结果为\begin{equation} \label{eq28} \frac{\partial V}{\partial v_k } = \int_D {\frac{\partial H\left( {\varPhi _\varOmega } \right)}{\partial v_k }\mbox{d}\varOmega } = \int_{\partial \varOmega } {\frac{\partial \varPhi _\varOmega }{\partial v_k }\frac{1}{\left\| {\nabla \varPhi _\varOmega } \right\|}\mbox{d}\varGamma } (28) \end{equation}

由式(27)和式(28)可知,该优化问题灵敏度分析的关键在于求解结构边界$\partial\varOmega$上的$\partial \varPhi _\varOmega / \partial v_k$和$\nabla \varPhi _\varOmega $,而$\varPhi_{\varOmega }$又与$\varPhi $$_{D}$和$\varphi_{i}$相关,其中$\varPhi_{D}$属于设计无关量.此外,式(18)采用的是min函数,$\varPhi_{\varOmega}$在边界$\partial \varOmega$上某点附近的值及其变化,一般仅与在该点处值为0的孔洞隐函数$\varphi_{i}$相关,所以$\partial \varPhi _\varOmega/\partial v_k$和$\nabla\varPhi _\varOmega $的求解最终归为计算$\partial \varphi _i / \partial v_k$和$\nabla\varphi _i $. 其中$\partial \varphi _i /\partial v_k$按照$v_{k}$所包含的变量类型可分别计算如下

$$ \begin{aligned} &\frac{\partial \varphi _i }{\partial x_i } =-\frac{2}{\varepsilon _i R_i \left( {\theta _i } \right)}\frac{\partial R_i \left( {\theta _i } \right)}{\partial x_i }\left( {\left| {\frac{x-x_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }} + \left| {\frac{y-y_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }}} \right)-\\ &\qquad \frac{2}{\varepsilon _i \left( {x-x_i } \right)}\left| {\frac{x-x_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }} (29)\\ & \frac{\partial \varphi _i }{\partial y_i } =-\frac{2}{\varepsilon _i R_i \left( {\theta _i } \right)}\frac{\partial R_i \left( {\theta _i } \right)}{\partial y_i }\left( {\left| {\frac{x-x_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }} + \left| {\frac{y-y_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }}} \right)-\\ &\qquad \frac{2}{\varepsilon _i \left( {y-y_i } \right)}\left| {\frac{y-y_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }} (30)\\ & \frac{\partial \varphi _i }{\partial \varepsilon _i } \!=\! \frac{-2\left| {x - x_i } \right|^{\frac{2}{\varepsilon _i }}\ln \left| {\dfrac{x - x_i }{R_i \left( {\theta _i } \right)}} \right|-2\left| {y - y_i } \right|^{\frac{2}{\varepsilon _i }}\ln \left| {\dfrac{y - y_i }{R_i \left( {\theta _i } \right)}} \right|}{\varepsilon _i ^2\left| {R_i \left( {\theta _i } \right)} \right|^{\frac{2}{\varepsilon _i }}} (31)\\ & \frac{\partial \varphi _i }{\partial P_{i,j} } = \frac{-2N_{j,p} \left( {\dfrac{\theta _i + \pi / 2}{{2\pi }}} \right)}{\varepsilon _i R_i \left( {\theta _i } \right)}\left( {\left| {\dfrac{x-x_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }} + \left| {\frac{y-y_i }{R_i \left( {\theta _i } \right)}} \right|^{\frac{2}{\varepsilon _i }}} \right) \quad(32) \end{aligned} $$

在式(29)和式(30)中,

$$ \frac{\partial R_i \left( {\theta _i } \right)}{\partial x_i } = \frac{y-y_i }{\left( {x-x_i } \right)^2 + \left( {y-y_i } \right)^2}\sum\limits_{j = 1}^n {\dfrac{\mbox{d}N_{j,p} \left( {\frac{\theta _i + \pi / 2}{{2\pi }}} \right)}{\mbox{d}\theta _i }P_{i,j} }(33)\\ \frac{\partial R_i \left( {\theta _i } \right)}{\partial y_i } = \frac{-\left( {x-x_i } \right)}{\left( {x-x_i } \right)^2 + \left( {y-y_i } \right)^2}\sum\limits_{j = 1}^n {\dfrac{\mbox{d}N_{j,p} \left( {\frac{\theta _i + \pi / 2}{{2\pi }}} \right)}{\mbox{d}\theta _i }P_{i,j} } \quad (34) $$

式中,B样条基函数$N_{j, p}$的导数计算方法可查看文献^[35].孔洞隐函数的梯度$\nabla \varphi _i $可直接计算如下\begin{equation}\label{eq35} \nabla \varphi _i = \left( {{\begin{array}{*{20}c} {\dfrac{\partial \varphi _i }{\partial x}} {\dfrac{\partial \varphi _i }{\partial y}} \\\end{array} }} \right) = \left( {{\begin{array}{*{20}c} {-\dfrac{\partial \varphi _i }{\partial x_i }} {-\dfrac{\partial \varphi _i }{\partial y_i }} \\\end{array} }} \right) (35)\end{equation}

本节将使用自适应泡泡法对经典悬臂梁结构和简支梁结构进行拓扑优化设计,以验证其有效性.由于下面两个算例的计算域和边界条件都是轴对称的,因此极点($x_{i}$,$y_{i})^{T}$不在对称轴上的孔洞应该被对称地引入.然而"孔洞影响区域"的设置可能会导致距离过近的对称孔洞无法被先后引入,为解决此问题,本节在优化过程中遵从"先同时引入位置对称的孔洞,再共同激活其影响区域"的原则.

算例中材料杨氏弹性模量和泊松比分别取为1和0.3,每个孔洞设有28个变量(参见式(23),其中$n=27$).体积上限$V_{lim}$设为初始结构体积的50%,优化算法选用Boss-Quattro$^{TM}$优化平台^[36]中的GCMMA(全局收敛移动渐近线算法),优化迭代终止条件为设计变量的最大变化率小于0.01%或迭代次数达到150.

悬臂梁结构的计算域和边界条件如图5所示,对此计算域采用了80$\times$40的有限胞元网格划分,网格节点即为计算拓扑导数的采样点,共81$\times$41个. 决定每步优化迭代中拓扑导数阈值$\overline {D_{T} }$的$m$取值为160,约为采样点总数的5%.孔洞的初始半径设为0.3(2个网格长度),图4(b)所示的设置孔洞影响区域的参数$r$取值为1.2(8个网格长度).一般而言,$m$值越大则每步引入的孔洞数量越多,$r$值越大则每步引入的孔洞数量越少.

图6绘制了悬臂梁结构的中间优化结果和第150步的最终优化结果,整个优化过程共自适应引入了17个孔洞,设计变量总数为476.由图6可见,在50步之后结构中没有引入新的孔洞,这说明在后面的迭代步中,拓扑导数值不大于阈值$\overline{D_{T} } $的$m$个采样点全部位于已有孔洞的影响区域之内.

图7反映了参数$m$和$r$对最终优化结果的影响规律.当只减小$m$到100时,则仅需引入9个孔洞就可得到与图6(h)相同的优化结果,此时变量数目为252,见图7(a);当只增大$r$到1.5时,则引入11个孔洞得到了与图6(h)相同的优化结果,见图7(b);如果在减小$m$的同时又增大$r$,则会因为孔洞引入数量过少而降低最终优化结果的拓扑复杂性,见图7(c);反之则会得到拓扑相对复杂的最终优化结果,见图7(d).

简支梁结构的计算域和边界条件如图8所示,对此计算域采用了160$\times$40的有限胞元网格划分,网格节点即为计算拓扑导数的采样点,共161$\times$41个. 参数$m$取值为200,约为采样点总数的3%.孔洞初始半径和参数$r$的取值同上个算例.图9绘制了简支梁结构的中间优化结果和第150步的最终优化结果,整个优化过程共自适应引入了17个孔洞,设计变量总数为476.

本文基于有限胞元方法、拓扑导数理论和光滑变形隐式曲线,提出了一种改进的自适应泡泡法.数值算例表明,该方法在继承传统泡泡法设计变量少、结果边界光滑清晰、无初始布局依赖性等优点的基础上,克服了其网格重划频繁、孔洞合并操作繁琐等不足.自适应泡泡法数值计算稳定性高,可通过调整拓扑导数阈值参数$m$和孔洞影响区域参数$r$,改变孔洞引入的频次和位置,从而达到控制最终优化结果拓扑复杂程度的目的.

值得一提的是,自适应泡泡法的三个主要模块——固定网格分析、孔洞自适应引入和孔洞隐式描述并不局限于采用前面所述的方式,也可分别使用传统有限元方法(边界单元刚度依据单元体分比进行打折)^[12]、渐进结构法的材料消除策略^[3]以及闭合B样条曲线(closedB-splines)^[28].这种替代方案能够以降低一定精度为代价,有效简化结构力学响应和灵敏度等的计算过程,从而大幅提升优化效率.相应的简化后的自适应泡泡法与渐进结构法有着较大的相似性,区别在于前者消除材料的方式是引入边界光滑可变的孔洞,而且孔洞引入的频次和位置可被自适应地调整.

进一步的研究将在自适应泡泡拓扑优化中考虑局部应力约束^[37-38]、精确变形约束^[39-40]、破损、安全设计原则^[41-42]以及增材制造工艺约束^[43-44]等,发展出一种面向工程设计与制造需求的高效率高精度拓扑优化方法.