引言

近20年来连续体结构拓扑优化获得了飞速的发展，成为当今工程结构概念性设计的强有力工具.各种拓扑优化方法，比如基于材料分布概念设计的均匀化方法^[1]、固体各项同性变密度法(solid isotropic material with penalization，SIMP) ^[2]和渐进结构优化方法(evolutionary structural optimization，ESO)^[3]，还有基于几何边界描述的水平集方法(level set method，LSM)^[4-5]和隐式拓扑描述函数法^[6-8]，还有国内隋允康等^[9-10]提出的独立连续映射法，郭旭等^[11-12]提出的可动变形构件法(moving morphable components，MMC)，都被用来解决各种结构优化设计问题.

考虑自重载荷的影响在许多大型土木工程结构设计问题中具有十分重要的意义，近年来激发了许多科研工作者的研究兴趣. Turteltaub等^[13]率先研究了自重载荷作用下的结构拓扑优化问题. Bruyneel和Duysinx^[14]也对考虑结构自重的拓扑优化问题做了详细的讨论，他们指出了求解此类问题所面临的3个主要困难，即低密度区域容易出现材料的寄生效应、目标函数非单调性问题及体积约束的无效性，为了解决低密度区域材料的寄生效应问题，他们提出了一种修正的不连续的SIMP模型. Yang等^[15]基于离散的ESO/BESO (bi-directional ESO)方法研究了考虑设计依赖自重载荷作用下的带有刚度约束的拓扑优化问题，但是他们的优化结果与Bruyneel和Duysinx的结果不同. Ansola等^[16]指出，对于考虑自重载荷的优化问题，传统ESO算法中单元敏度数的计算方法并不能使优化问题收敛到一个最优解，为此，他们提出了一种新的敏度计算策略来计算敏度，提高算法的收敛性，但是他们人为地增加体积约束上限的做法被认为是缺乏理论基础的^[17].高彤等^[18]提出了一种可变参数的“rational opproximation of material properties (RAMP)”模型，基于此模型分析了多种材料插值格式对惯性载荷作用下结构拓扑优化结果和迭代过程的影响，而且对“没有材料就没有惯性”的理论最优解问题进行了阐述和研究.张晖等^[19]针对自重载荷作用下连续体结构拓扑优化中遇到的常见问题，提出了一种将RAMP插值模型与平均敏度过滤技术相结合的求解策略，并详细讨论了不同模型参数对拓扑优化结果的影响. Huang等^[20]将BESO方法与RAMP插值模型相结合，详细讨论了考虑设计依赖自重载荷和固定载荷作用的连续体结构拓扑优化区别. Holmberg等^[21]通过对质量密度进行线性缩放以减小低密度区域载荷与刚度比值的大小，从而获得一个无寄生效应的优化结果.另外，导重法^[17]和修改的梯度投影法^[22]也被用来研究考虑自重载荷作用的结构拓扑优化问题.

总结上面处理考虑自重载荷的拓扑优化问题的方法主要可以分为3类：第1种是对单元刚度的插值格式进行修正^{[14, 22]}；第2种是对单元质量密度的插值格式进行修正^{[18, 21-22]}；第3种是采用RAMP插值模型^[17-20].这些研究表明，对于自重载荷作用的优化问题，RAMP插值模型比SIMP插值模型更有效.虽然上述方法在处理自重结构优化问题是有效的，但是需要对刚度或质量密度的插值格式进行修正，甚至有些方法都不能保证插值格式的连续性，而且上述工作更多的局限于二维平面的优化问题.

针对这些问题，本文将基于非均匀有理B样条(NURBS)基函数的等几何分析方法与拓扑描述函数法结合，提出一种有效求解三维自重结构优化问题的拓扑描述函数法.等几何分析是由Hughes等^[23]最先提出的，它是一种将传统的有限单元法与基于NURBS的计算机辅助设计工具相结合来求解偏微分方程的新数值方法.目前等几何分析已被广泛应用于结构分析^[23-26]、结构优化^[27-30]等各个领域.本文用等几何方法代替传统的有限元方法进行三维实体结构的位移场分析，以每个控制点所对应的拓扑描述函数值为设计变量，详细推导了基于NURBS插值的考虑自重载荷作用的三维结构拓扑优化模型及其灵敏度列式，最后通过多个数值算例验证了本文方法的稳定性和有效性.

1 基于NURBS插值的拓扑描述函数法 1.1 NURBS实体

三维NURBS实体定义如下^[31]

$ V(\xi, \eta, \zeta ) = \sum\limits_{i = 1}^n {\sum\limits_{j = 1}^m {\sum\limits_{k = 1}^l {P_{i, j, k} R_{i, j, k} } } } (\xi, \eta, \zeta ) $

(1)

其中，$R_{i, j, k} (\xi, \eta, \zeta )$为三变量NURBS基函数，表示为

$ R_{i, j, k}^{p, q, r} (\xi, \eta, \zeta ) =\\ \quad \dfrac{N_{i, p} (\xi )M_{j, q} (\eta )L_{k, r} (\zeta )\omega _{i, j, k} }{\sum\limits_{\bar {i} = 1}^n {\sum\limits_{\bar {j} = 1}^m {\sum\limits_{\bar {k} = 1}^l {N_{\bar {i}, p} (\xi )M_{\bar {j}, q} (\eta )L_{\bar {k}, r} (\zeta )\omega _{\bar {i}, \bar {j}, \bar {k}} } } } } $

(2)

式中，$N_{i, p}, M_{j, q}, L_{k, r} $分别为定义在开放性节点向量${\boldsymbol\varXi}_\xi = \left\{ {\xi _1, \xi _2, \cdots, \xi _{n + p + 1} } \right\}$，${\boldsymbol\varXi}_\eta = \left\{ {\eta _1, \eta _2, \cdots, \eta _{m + q + 1} } \right\}$，${\boldsymbol\varXi} _\zeta = \left\{ {\zeta _1, \zeta _2, \cdots, \zeta _{l + r + 1} } \right\}$上的$p$次、$q$次和$r$次B样条基函数，$n, m, l$分别代表$\xi $，$\eta $和$\zeta $方向上控制点个数，$P_{i, j, k} $为三维实体控制点，非负实数$\omega _{i, j, k} $为其对应的权重，当所有的权重都为1时，NURBS基函数就退化为B样条基函数.除非特殊声明，本文所有的算例均采用二次NURBS基函数进行插值.

1.2 基于NURBS插值的拓扑描述函数法

在拓扑描述函数法中，用来描述结构拓扑和形状的隐式函数可以定义为^[8]

$ \left.\begin{array}{ll} \phi (x) = 0, & \partial \varOmega \\ \phi (x) > 0, & \hbox{在} \ \varOmega \hbox{内} \\ \phi (x) < 0, & \hbox{在} \ D / \varOmega \hbox{内} \end{array} \!\!\right\} $

(3)

其中，${D}$是基准设计域，$\varOmega, \partial \varOmega $分别代表最优设计实心区域的边界.

结构拓扑优化的一个关键任务就是将不适定的离散的0-1优化问题转化为适定的连续优化问题.在TDF方法中，通过Heaviside函数很容易实现这一目的

$ H(\phi (x)) = \left\{ \!\!\begin{array}{ll} 0, & \phi (x) < 0 \\ 1, & \phi (x) \geqslant 0 \end{array} \right. $

(4)

目前有关拓扑描述函数法的大部分文献是基于有限元网格的，采用Q4单元形函数插值单元设计变量场^{[6, 8]}，本文通过NURBS基函数来插值单元的设计变量场，可表示如下

$ \phi ^h (x) = \sum\limits_{I = 1}^{n_{\rm cp}} R_I (\xi, \eta, \zeta )\phi _I $

(5)

其中$R(\xi, \eta, \zeta )$为式(2)确定的三变量NURBS基函数，$\phi _I $代表一个单元内第$I$个控制点所对应的拓扑描述函数值，$n_{cp} = (p + 1)\times (q + 1)\times (r + 1)$为一个单元相关联的控制点个数.

为了获得清晰的结构拓扑，本文对弹性模量和质量密度分别进行惩罚

$ \left.\begin{array}{l} E(x) = \Big [\sum\limits_{I = 1}^{n_{\rm cp}} R_I (\xi, \eta, \zeta )H(\phi _I ) \Big]^{g_1 }E_0 \\ \rho (x) = \Big [\sum\limits_{I = 1}^{n_{\rm cp}} R_I (\xi, \eta, \zeta )H(\phi _I ) \Big]^{g_2 }\rho _0 \end{array}\!\!\right\} $

(6)

式中，$E_0$和$\rho _0 $分别为完全实体材料的杨氏模量和质量密度，$g_1$和$g_2 $均为惩罚指数，本文取$g_1 = 3$, $g_2 = 1$.

类似于有限元等参元的思想，在等几何分析中，NURBS基函数被用来同时对几何区域和位移场进行插值

$ \left. \begin{array}{l} \mathit{\boldsymbol{x}}(\xi, \eta, \zeta ) = \sum\limits_{I = 1}^{{n_{{\rm{cp}}}}} {{R_I}} (\xi, \eta, \zeta ){\mathit{\boldsymbol{\overline x}} _I} = \mathit{\boldsymbol{RX}}\\ \mathit{\boldsymbol{u}}(\xi, \eta, \zeta ) = \sum\limits_{I = 1}^{{n_{{\rm{cp}}}}} {{R_I}} (\xi, \eta, \zeta ){\mathit{\boldsymbol{\overline u}} _I} = \mathit{\boldsymbol{RU}} \end{array} \right\} $

(7)

其中，$\bar {\boldsymbol x}_I = \{ x_I, y_I, z_I\}$，$\bar {\boldsymbol u}_I = \{u_I, v_I, w_I \}$分别代表第$I$个控制点所对应的物理坐标向量和位移场向量.应变位移阵${\boldsymbol B}$可以由几何方程得到

$ {\boldsymbol \varepsilon} ={\boldsymbol L}{\boldsymbol u} \to {\boldsymbol \varepsilon } = {\boldsymbol B}{\boldsymbol U } $

(8)

然后通过虚位移原理推导单元刚度阵如下

$ {\boldsymbol K}_e = \int_{V_e } E(x){\boldsymbol B}^{\rm T} {\boldsymbol C}_0 {\boldsymbol B}{\rm{d}} V_e $

(9)

其中${\boldsymbol C}_0 $为$3\times 3$的弹性阵.

1.3 Heaviside函数的光滑性处理

在基于梯度的优化方法中，目标函数及约束函数关于设计变量的敏度信息对于数值求解优化问题是极其重要的.考虑到Heaviside函数的非光滑特征，本文采用文献[8]里面的方法对Heaviside函数进行光滑逼近

$ H_\varepsilon (\phi ) = \left\{ \begin{array}{ll} \delta, & {\phi \leqslant - \varepsilon } \\ {[1 + \sin (\pi \phi / 2\varepsilon )] / 2}, & { - \varepsilon < \phi < \varepsilon } \\ 1, & {\phi \geqslant \varepsilon } \end{array} \right. $

(10)

这里$\varepsilon $和$\delta $均是一个很小的正数，其中$\varepsilon $为Heaviside函数的磨光参数，$\delta $是为了避免刚度阵奇异而引进的参数，一般取$\varepsilon = 0.3 \sim 0.5, \delta = 10^{ - 5} \sim 10^{ - 3}$，本文取$\delta = 10^{ - 3}$.

正则化的Heaviside函数$H_\varepsilon (\phi )$的导函数可表示为

$ H'_\varepsilon (\phi ) = \left\{ \begin{array}{cl} {\dfrac{\pi }{4\varepsilon }\cos \dfrac{\pi \phi }{2\varepsilon } }, & {\left| \phi \right| < \varepsilon } \\ 0, & {\left| \phi \right| \geqslant \varepsilon } \end{array} \right. $

(11)

2 拓扑优化模型及灵敏度分析 2.1 考虑自重载荷的优化模型

结构拓扑优化的目的是寻求满足某些固定约束条件的结构的最优的材料分布，以获得某种最优的结构行为，如结构的重量最轻或刚度最大化.本文考虑自重载荷作用下的最小柔顺性问题.以许可材料体积为约束的最小柔顺性问题列式可表示如下

$ \left. \begin{array}{l} {\rm{Find}}:\;\;{\boldsymbol\phi} \\ {\rm{Min}}:\;\;c = {\mathit{\boldsymbol{F}}^{\rm{T}}}\mathit{\boldsymbol{U}}\\ {\rm{S}}.{\rm{t}}.:\;\;\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{KU}} = \mathit{\boldsymbol{F}}}\\ {V-{V_0}{f_{\rm{v}}} \le 0} \end{array}} \right. \end{array} \right\} $

(12)

其中$c$指结构的柔顺性目标函数，${\boldsymbol\phi} = [\phi _1, \phi _2, \cdots, \phi _N]$为优化问题设计变量的集合，在优化过程中不断更新. ${\boldsymbol K}$为结构的整体刚度阵，由单元刚度阵${\boldsymbol K}_e $组装形成，${\boldsymbol U}, {\boldsymbol F}$分别为位移向量和载荷列阵. $V$，$V_0 $分别指实际的材料体积和整个设计域的体积，$f_{\rm v} $为准许的材料体积分数.

对于考虑自重载荷作用的结构拓扑优化问题，结构载荷具有设计依赖性.受自重载荷作用的等效单元载荷向量可表示为

$ {\boldsymbol F}_e = \bar {\boldsymbol F}_e \int_{V_e } {\rho (x)g} {\rm{d}} V_e $

(13)

其中$g$为重力加速度，$\rho (x)$是由式(6)定义的质量密度.对于二维问题，$\bar {\boldsymbol F}_e $是一个包含$2\times n_{\rm cp}$个元素的列向量，可表示如下

$ {\boldsymbol F}_e^i = \left\{ \begin{array}{cl} 0, & {i = 2 j - 1} \\ - \dfrac{1}{n_{\rm cp}}, & {i = 2 j} \end{array} \right. $

(14)

对于三维问题，$\bar {\boldsymbol F}_e $是一个包含$3\times n_{\rm cp}$个元素的列向量

$ {\boldsymbol F}_e^i = \left\{\begin{array}{cl} 0 & {i = 3 j - 1, \ \ 3 j - 2} \\ - \dfrac{1}{n_{\rm cp}} & {i = 3j} \end{array} \right. $

(15)

其中$j = 1, 2 \cdots, n_{\rm cp}$.

在优化问题(12)中，每个单元的材料体积可表示为

$ V_e = \int_{V_e } \sum\limits_{I = 1}^{n_{\rm cp}} R_I H_\varepsilon (\phi _I ) {\rm{d}} V_e $

(16)

2.2 灵敏度分析

用式(10)中光滑后的Heaviside函数代替式(6)中的Heaviside函数，然后将弹性模量和质量密度分别对设计变量求导得

$ \left. \begin{array}{l} \frac{{\partial E(x)}}{{\partial {\phi _I}}} = {g_1}{[\sum\limits_{I = 1}^{{n_{cp}}} {{R_I}{H_\varepsilon }({\phi _I})}]^{{g_1} - 1}}{R_I}{H_{\varepsilon '}}({\phi _I}){E_0}\\ \frac{{\partial \rho (x)}}{{\partial {\phi _I}}} = {g_2}{[\sum\limits_{I = 1}^{{n_{cp}}} {{R_I}{H_\varepsilon }({\phi _I})}]^{{g_2} -1}}{R_I}{H_{\varepsilon '}}({\phi _I}){\rho _0} \end{array} \right\} $

(17)

则单元刚度阵、单元外载荷向量及单元的材料体积对设计变量的导数分别表示为

$ \left. \begin{array}{l} \frac{{\partial {\mathit{\boldsymbol{K}}_e}}}{{\partial {\phi _I}}} = \int_{{V_e}} {\frac{{\partial E(x)}}{{\partial {\varphi _I}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}} {\mathit{\boldsymbol{C}}_0}\mathit{\boldsymbol{B}}{\rm{d}}{V_e}\\ \frac{{\partial {{\bf{F}}_e}}}{{\partial {\phi _I}}} = {{\bar F}_e}\int_{{V_e}} {\frac{{\partial \rho (x)}}{{\partial {\phi _I}}}g} {\rm{d}}\\ {V_e}\frac{{\partial {V_e}}}{{\partial {\phi _I}}} = \int_{{V_e}} {{R_I}{H_{\varepsilon '}}({\phi _I})} {\rm{d}}V \end{array} \right\} $

(18)

利用伴随法可推导目标函数的灵敏度为

$ \dfrac{\partial c}{\partial \phi _I } = 2\dfrac{\partial {\boldsymbol F}^{\rm T}}{\partial \phi _I }{\boldsymbol U}-{\boldsymbol U}^{\rm T}\dfrac{\partial {\boldsymbol K}}{\partial \phi_I }{\boldsymbol U} $

(19)

简化为

$ \dfrac{\partial c}{\partial \phi _I } = 2\sum\limits_{Ie = 1}^{N_{Ie} } {\dfrac{\partial {\boldsymbol F}_{Ie}^{\rm T} }{\partial \phi _I }{\boldsymbol U}_{Ie} } - \sum\limits_{Ie = 1}^{N_{Ie} } {U_{Ie}^{\rm T} \dfrac{\partial {\boldsymbol K}_{Ie} }{\partial \phi _I } {\boldsymbol U}_{Ie} } $

(20)

其中$N_{Ie} $为与第$I$个控制点相关联的单元个数.

体积约束关于设计变量的灵敏度可表示为

$ \dfrac{\partial V}{\partial \phi _I } = \sum\limits_{Ie = 1}^{N_{Ie} } {\dfrac{\partial V_{Ie} }{\partial \phi _I }} $

(21)

式(20)和式(21)表明只有与第$I$个控制点相关联的单元才对目标函数和体积约束的灵敏度有贡献.

2.3 优化算法

由式(12)定义的基于拓扑描述函数理论的优化问题可以通过不同的优化算法来求解，如OC准则法^[6]、数学规划方法^[8]等，本文采用Svanberg教授提出的移动渐进线(method of moving asymptotes, MMA)方法^[32]求解.优化算法主要包括以下6个步骤：

(1)构建设计域：定义设计域几何参数，使用NURBS构建优化模型的几何区域.

(2)定义模型参数：定义实体材料杨氏模量和质量密度及材料容许的体积分数等模型参数.

(3)等几何分析：运用等几何分析方法求解位移场.

(4)灵敏度分析：计算目标函数值，及通过式(20)和式(21)分别计算目标函数和材料体积对设计变量的敏度值.

(5)优化模型求解：用MMA算法求解优化模型，进行设计变量更新.

(6)收敛性判断：本文以最大的迭代次数来判断优化过程是否停止；若收敛，则迭代停止，输出优化结果；若不收敛，则返回步骤(3)，继续迭代直到收敛为止.

3 数值算例与讨论

本节通过多个数值算例来验证本文算法的有效性和准确性.除非特殊声明，在所有算例中磨光参数取$\varepsilon =0.5$，设计变量$\phi_I$的上下限取为$[-0.1, 0.1]$，所有结构的初始设计变量均附值为0.01；完全实体材料的杨氏模量$E_0 $和泊松比分别设置为1 000.0和0.3、质量密度$\rho _0 = 1$，板的厚度均设置为0.1.

3.1 固定外载荷的最小柔顺性问题：L形梁

第一个算例考虑一个带有曲率的受固定单位载荷的L形梁优化问题，几何区域和边界条件如图 1所示.以占设计域$30\% $的体积作为约束.需要注意的是，在基于有限元网格的拓扑优化方法中，L形梁区域的离散通常需要分成两部分，对于带有曲率的L形梁，建模及网格划分更是困难，而在本文方法中，带有曲率的L形梁可以很方便地建模和分析.

图 2(a)和图 2(b)分别显示了本文方法基于3 332和12 804个控制点(设计变量)获得的L形梁的最优拓扑构型.可以看出本文方法在没有借助于灵敏度过滤情况下基于两种不同的离散网格所获得的拓扑构型本质上很相似，没有出现明显的网格依赖性问题，随着网格的细化，拓扑构型更加清晰，这说明本文方法是十分有效的.

图 3显示了本文方法基于12 804个控制点获得的柔顺性目标函数的迭代历程，其中还包括第10, 30和50个迭代步的拓扑构型.可以看出，本文方法迭代过程稳健，收敛性快，最终得到一个没有孤岛现象、边界清晰的拓扑构型.

3.2 2D自重拱设计

第二个算例考虑一个几何尺寸为$2\times 1$的受自重载荷作用的2D拱形结构设计问题.设计域和边界条件如图 4(a)所示，结构底边左右两个端点完全固支.以$50\% $的材料体积为约束. 图 4(b)显示了本文方法基于$80\times 40$个均匀分布的控制点获得的最优的拓扑构型. 图 5给出了目标函数和体积分数的迭代历史.从图 4和图 5看出，本文方法迭代过程稳健，优化结果没有出现文献[14]中报道的低密度区域材料的寄生效应，而且基于本文方法获得的目标函数在整个优化过程中始终保持单调递减的性质，体积分数是先减小后增大最后保持在$35.1\% $附近.

对于考虑自重载荷作用的结构拓扑优化问题来说，另外一个关键特征就是体积约束的无效性.为了更清楚地阐述这一特性，针对上面的优化问题，考虑体积分数约束从$10\% $到$100\% $变化，图 6显示了本文方法在不同的体积分数约束下获得的最优拓扑的实际的体积分数变化曲线.可以看出，对于2D自重载荷作用的拓扑优化问题，体积约束并不是一个有效的约束条件.

3.3 3D自重拱设计

第3个算例考虑一个几何尺寸为$2\times 1\times 1$的受自重载荷作用的3D拱形结构设计问题.设计域和边界条件如图 7(a)所示，结构底面左右两边完全固支.以$40 \% $的材料体积为约束.结构通过$26\times 11\times 13$个均匀分布的控制点来离散.

图 8显示了本文方法获得的最优的拓扑构型.可以看出优化过程产生了一个跨越两端支撑的拱形结构，以便最佳地承受其自身重量，拓扑构型与二维结果一致.

为了对比固定载荷和自重载荷作用对结构拓扑构型的影响，针对上面的优化问题考虑上表面受均布的单位载荷作用情形(上表面包含一层不可设计域)，设计域和边界条件如图 7(b)所示.以占设计域$40\% $的体积为约束，图 9显示了本文方法基于$26\times 11\times 13$个均匀分布的控制点获得的最优拓扑构型.可以看出优化拓扑类似于赵州桥结构，与自重拱形结构不同的是，它多了几根垂直支柱以支撑桥面上所受到的均布载荷.

图 10显示了3D自重拱(3D arch)和3D桥(3D bridge)结构设计过程中目标函数和体积分数的迭代历程.可以看出，本文方法迭代过程稳健，收敛快；对于3D自重拱问题，目标函数是单调递减的，最终的体积分数$(46.2\% )$要高于约束的体积分数$(40\% )$，说明对于3D考虑结构自重的优化问题，体积约束也是无效约束；对于3D桥结构设计问题，目标函数先增大后减小，在迭代25步后基本保持不变，体积分数先减小后增大最后保持为$40\% $的约束体积分数不变.

4 结论

本文提出一种基于NURBS插值的三维拓扑描述函数法，研究了考虑设计依赖自重载荷作用的2D/3D结构设计问题，详细推导了基于NURBS插值的考虑自重载荷作用的三维结构拓扑优化模型及其灵敏度列式，多个算例验证了本文方法的有效性和稳定性，主要得出以下结论：

(1)高阶NURBS基函数被同时用于三维NURBS实体片中的几何场、位移场及设计变量场插值，实现了几何模型、分析模型和优化模型的有效统一，确保了位移场及设计变量场的高阶连续性；

(2)无需对刚度和质量密度的插值格式进行修正，本文方法能够有效克服考虑自重载荷作用的优化问题中出现的材料的寄生效应及目标函数的非单调性问题；

(3)研究表明，本文方法优化迭代过程稳健，收敛性快.