两材料自由曲面壳形状-拓扑联合优化研究

Research on shape-topology joint optimization of free-form surface shells with two materials

摘要

针对由两种材料构成的自由曲面壳体结构，考虑形状和材料分布的耦合作用，通过形状-拓扑联合优化算法，将形状优化步和拓扑优化步交替迭代，获得两材料自由曲面壳的最优形状与材料最优分布形式.以四角点简支和四边简支两种常见壳体结构为研究对象，重点分析两种材料的弹性模量比值（E₁∶E₂）和体积比值（V₁∶V₂）对优化结果的影响.研究结果表明，材料弹性模量比值和体积比值对结构形态及材料分布具有显著影响.随着弹性模量比值减小，最优曲面曲率降低，高弹模材料向主要受力部位聚集；随着体积比值增加，最优曲面形态趋于平缓.当V₁∶V₂<5∶5时，高弹模材料集中于主要受力部位；当V₁∶V₂>5∶5时，高弹模材料沿主要受力部位的连接区域分布，形成有效传力路径.

Abstract

This paper focuses on the free-form surface shell structures composed of two materials and considers the coupling effect between shape and material distribution. By employing a shape-topology joint optimization algorithm， the shape optimization step and topology optimization step are iterated alternately to obtain the optimal shape and the optimal materials distribution for the two-material free-form surface shell. Taking the four-corner simply supported shell and the four-edge simply supported shell as the research objects， this study focuses on analyzing the influence of the elastic modulus ratio （E₁∶E₂） and volume ratio （V₁∶V₂） of the two materials on the optimization results. The research results show that the elastic modulus ratio and volume ratio of the two materials have a significant impact on the structural morphology and material distribution. As the elastic modulus ratio decreases， the curvature of the optimal surface decreases， and the high-modulus material gathers in the main load-bearing parts of the structure； as the volume ratio increases， the optimal surface tends to be gentle. When the volume ratio is less than 5∶5， the high-modulus material is mainly distributed in the main load-bearing parts of the structure； when the volume ratio is greater than 5∶5， the high-modulus material is also distributed along the connecting areas between the main load-bearing parts， forming an effective force transmission path.

Graphical abstract

关键词

自由曲面壳体 / 两种材料 / 形状-拓扑联合优化 / 材料参数影响

Key words

free-form surface shells / two materials / shape-topology joint optimization / influence of material parameters

自由曲面壳体结构因造型美观、空间利用率高以及力学性能卓越等优势，在建筑结构领域广受青睐.区别于规则曲面，自由曲面是指无法用一个或几个解析函数表达的曲面，例如悉尼歌剧院、北京大兴机场航站楼就是典型的自由曲面结构.相较于规则壳体结构，如何获得满足设计要求的最优曲面形状，即形状优化，是自由曲面壳体设计中的一个重要问题.随着自由曲面结构形式的推陈出新，自由曲面壳体不再限于由单一材料构成，例如位于德国德累斯顿工业大学校园内的CUBE壳结构，就是由玻璃纤维复合材料和混凝土两种材料构成.这就又将拓扑优化引入到结构设计中，即不同材料如何分布最优.本文的研究将针对两材料自由曲面壳体的形状和拓扑优化问题开展.

国内外学者已针对自由曲面形状优化问题开展了大量研究.例如，李欣等^［1］将非均匀有理B样条技术（non-uniform rational B-spline，NURBS）与寻求结构应变能最小的梯度算法相结合，提出了NURBS-GM法进行自由曲面形状优化.王磊等^［2］基于NURBS技术，通过直接调整型值点坐标来优化自由曲面形态，突破了传统控制点调整方法的局限性.San等^［3-4］考虑混凝土材料损伤效应，提出了融合NURBS参数化与Mazars损伤模型的自由曲面混凝土壳的形状优化框架，并在此基础上，进一步将研究扩展至材料徐变效应，提出基于广义Maxwell模型和NURBS参数化的约束优化方法，验证了徐变模型对于壳体长期性能分析的必要性.Jiang等^［5］通过NURBS技术描述曲面形状，并采用基于梯度的优化方法，探讨了水平投影面积约束对带开口自由曲面结构静力和动力性能的影响.

拓扑优化目前已发展出一些较成熟的方法，例如固体各向同性材料惩罚模型法（solid isotropic material with penalization， SIMP）、双向渐进结构优化法（bi-directional evolutionary structure optimization， BESO）等.为解决多材料的分布问题，一些学者开展了多种材料的拓扑优化研究.Zuo等^［6］提出了一种新型的多种材料SIMP模型，并解决了插值模型中的导数不连续问题.Ghabraie^［7］通过改进灵敏度过滤方案并引入逐步优化策略，提出了一种改进的BESO法，解决了传统方法的数值不稳定问题.刘继凯等^［8］提出了一种基于Ordered SIMP插值模型的点阵-实体复合结构拓扑优化方法，通过算法自动分配实体材料和点阵材料的比例.Li等^［9］通过引入材料利用效率的概念，提出了多材料双向进化结构优化方法.刘宏亮等^［10］通过建立离散变量下的多材料拓扑优化模型，研究了离散变量组合复杂性的问题.

对于由多种材料构成的自由曲面壳体结构，同时存在形状优化和拓扑优化问题.对此，学者们主要探索了三种优化策略，即单独优化、同时优化以及联合优化.单独优化是将形状和拓扑优化作为相互独立的优化流程，即先对结构进行形状或拓扑优化，再在此基础上进行另外一种优化^［11-12］.单独优化思路简单，易于实现，但是没有考虑结构形状和拓扑之间的耦合作用，导致优化结果并非真正的“最优解”.同时优化是指在一次优化计算中，同时进行形状和拓扑优化，即优化变量既包括形状参数又包括拓扑参数^［13-16］.同时优化考虑了结构形状和拓扑之间的相互影响，但是由于形状和拓扑变量同时在一个优化步内进行优化计算，优化变量较多，算法复杂，优化的时间比较长.联合优化是将拓扑优化步和形状优化步交替进行，使优化过程中结构的形状和拓扑能够相互影响并最终收敛.联合优化是目前应用较多的一种策略，例如Ansola等^［17-18］对自由曲面形状和拓扑的联合优化进行了探索，首先采用Coons patches插值技术建立合适的曲面边界并创建Coons曲面，对Coons曲面进行若干次形状优化，在此基础上，基于分层微结构，通过材料密度和方向参数优化并结合线性规划和Simplex算法对曲面进行若干步拓扑优化.这种方法既在一定程度上考虑了形状和拓扑之间的相互影响，又兼顾了计算效率.但该方法具有一定的局限性，即Coons曲面无法在不改变曲面边界形状的情况下改变曲面的内部形状，且在收敛条件方面，仅通过最大迭代次数和整体收敛性进行判断，未明确双收敛条件.总体而言，针对壳结构的联合优化已有一些可行方法，但在曲面参数化、优化策略具体实现方法等方面仍有改进空间.同时也发现，对于不同材料弹性模量、体积分配等因素对优化结果的影响也缺乏深入探讨.

为此，本文提出一种改进的形状-拓扑联合优化方法.其中形状优化采用基于NURBS曲面的梯度法，能够支持复杂自由曲面在边界固定情况下的局部形状调整；拓扑优化方法则引入两种材料SIMP变密度法，结合交换活性相法和优化准则法处理两种材料分布问题.此外，该方法通过设置双收敛条件，解决联合优化的收敛问题.基于改进后的形状-拓扑联合优化方法开展系统分析，着重探讨两种材料的弹性模量比值、体积比值对最优曲面形状和最优拓扑的影响.

1 自由曲面形状-拓扑联合优化方法

1.1 基于NURBS的曲面形状优化

本文采取的自由曲面形状优化方法以结构的最小应变能作为优化目标，通过调整结构的曲面形状来实现优化目的.首先使用非均匀有理B样条（NURBS）曲面进行几何建模（图1）.

NURBS曲面的表达式为：

\mathit{S}\left(u,v\right)=\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{\mathit{P}}_{i,j}\cdot w_{i,j}\cdot N_{i,k}\left(u\right)\cdot N_{j,l}\left(v\right)}{{\displaystyle \sum _{i=\mathrm{1}}^n}{\displaystyle \sum _{j=\mathrm{1}}^m}{w}_{i,j}\cdot N_{i,k}\left(u\right)\cdot N_{j,l}\left(v\right)}

\mathit{U}=\left\{\underset{k+\mathrm{1}}{\underbrace{\mathrm{0},\cdots ,\mathrm{0}}},u_{k+\mathrm{2}},\cdots ,u_n,\underset{k+\mathrm{1}}{\underbrace{\mathrm{1},\cdots ,\mathrm{1}}}\right\}

\mathit{V}=\left\{\underset{l+\mathrm{1}}{\underbrace{\mathrm{0},\cdots ,\mathrm{0}}},u_{l+\mathrm{2}},\cdots ,u_m,\underset{l+\mathrm{1}}{\underbrace{\mathrm{1},\cdots ,\mathrm{1}}}\right\}（1）

式中： P_i，j 表示u方向上编号为i、v方向上编号为j的控制点；w_i，j 是权因子；N_i，k （u）和N_j，l （v）分别是定义在矢量 U 和 V 上的B样条基函数.

NURBS曲面通过调整控制点坐标和权因子来改变曲面形状.考虑到权因子对曲面形状的影响较控制点坐标小^［1］，本文中仅以控制点坐标为优化变量，权因子设为定值.

采用梯度法对上述优化问题进行求解.梯度法通过计算目标函数在变量处的导数（梯度），并让变量沿着函数负梯度方向变化，使目标函数达到最小的极值.具体步骤包括：采用差分法求解变量梯度确定优化方向，通过黄金分割法确定优化步长，最后根据优化方向和优化步长更新曲面控制点坐标.

整体曲面结构的应变能SE可以表示为：

\mathrm{S}\mathrm{E}(\mathit{z})=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{D}}^{\mathrm{T}}\mathit{F},\begin{array}{c}\end{array}\mathit{z}=\left\{z_{\mathrm{1}},z_{\mathrm{2}},\cdots ,z_k\right\}（2）

式中： z 为曲面形状优化变量即控制点的竖向坐标矩阵； D^T为节点位移矩阵； F 为节点荷载矩阵.可以看出结构的应变能是一个标量，k为曲面控制点的个数.

梯度法的迭代格式如式（3）所示：

z_i^{\left(h+\mathrm{1}\right)}=z_i^{\left(h\right)}-{\alpha }^{\left(h\right)}\frac{\partial \mathrm{S}\mathrm{E}(\mathit{z})}{\partial z_i}（3）

式中：\alpha为优化步长，h为第h步优化步.梯度法的收敛条件可以表示为：

\left(\mathrm{S}{\mathrm{E}}^{(h)}-\mathrm{S}{\mathrm{E}}^{(h+\mathrm{1})}\right)/\mathrm{S}{\mathrm{E}}^{(h)}\le {\epsilon }_{\mathrm{s}}（4）

式中：SE^（h）为第h个优化步的结构应变能；ε_s为给定的收敛标准，本文取ε_s=0.001.当式（4）满足时，表示结构应变能随着优化的进行而逐渐稳定，结构形状优化收敛.

1.2 基于两种材料SIMP法的曲面拓扑优化

本文拓扑优化采用的两种材料SIMP变密度法可通过调整自由曲面上两种材料的分布，实现结构在给定荷载条件下的最小应变能，如图2所示.通过引入单元拓扑密度作为设计变量，将其与单元弹性模量建立对应关系，并采用SIMP插值模型来描述两种材料的最优分布问题.根据交换活性相法确定单元拓扑密度的变化域，再采用优化准则法更新材料的拓扑密度，最终获得材料的最优分布.

交换活性相法通过动态划分“活性相”（可改变材料状态的区域）和“非活性相”（材料状态固定的区域）以提升优化效率^［19］.其核心原理是通过设定材料密度阈值，仅对接近阈值的单元进行状态更新，避免全区域迭代导致的计算冗余，从而减少无效迭代并提升数值稳定性；采用基于Kuhn-Tucker最优条件的优化准则法，即通过建立目标函数与约束条件的灵敏度关系，迭代调整设计变量直至满足收敛准则^［20］.该方法计算效率高，适用于含体积约束的拓扑优化问题.根据优化准则法得出拓扑优化的收敛条件需满足拓扑密度变化率≤ε_T，其中ε_T为给定的收敛标准，本文取ε_T=0.001.

两种材料SIMP插值模型可以表示为：

E_i={\displaystyle \sum _{j=\mathrm{1}}^m}{E}_j{\left({\rho }_{i,j}\right)}^p,i=\mathrm{1,2},\cdots ,n（5）

式中：m为材料种数；n为单元个数；E_i 代表第i个单元的弹性模量；E_j 为当第j种材料存在时的单元弹性模量；\rho_i，j 为第i个单元第j种的材料密度值；p为惩罚因子.

本文拓扑优化采用基于Kuhn-Tucker条件的优化准则法，Kuhn-Tucker条件是一种约束极值的最优化条件，可以表示为：

-\frac{\partial \mathrm{S}\mathrm{E}}{\partial {\rho }_{i,j}^{\left(h\right)}}-\lambda \frac{\partial \left({\displaystyle \sum _{i=\mathrm{1}}^n}{v}_i{\rho }_{i,j}^{\left(h\right)}-V^{\mathrm{*}}\right)}{\partial {\rho }_{i,j}^{\left(h\right)}}=\mathrm{0}（6）

其中：\lambda为拉格朗日乘子；V^{\mathrm{*}}为材料的给定体积约束上限；V_i 为第i个单元的体积.结构应变能对密度的梯度\frac{\partial \mathrm{S}\mathrm{E}}{\partial {\rho }_{i,j}^{(h)}}表达式为：

\frac{\partial \mathrm{S}\mathrm{E}}{\partial {\rho }_{i,j}^{\left(h\right)}}=-{\displaystyle \sum _{i=\mathrm{1}}^n}p{\left({\rho }_{i,j}^{\left(h\right)}\right)}^{p-\mathrm{1}}{\mathit{d}}_i^{\mathrm{T}}{\mathit{k}}_i{\mathit{d}}_i（7）

式中： d_i 为单元位移向量； k_i 为单元刚度矩阵^［20］.

由式（6）可得出优化准则法中拓扑密度的迭代格式：

{\rho }_{i,j}^{\left(h+\mathrm{1}\right)}=\left\{\begin{array}{c}\mathrm{m}\mathrm{a}\mathrm{x}\left(l_{i,j}^{\left(h\right)},{\rho }_{i,j}^{\left(h\right)}-\mathrm{m}\mathrm{o}\mathrm{v}\mathrm{e}\right)\\ \mathrm{m}\mathrm{i}\mathrm{n}\left(\mathrm{m}\mathrm{i}\mathrm{n}\left(u_{i,j}^{\left(h\right)},r_{i,j}^{\left(h\right)}\right),{\rho }_{i,j}^{\left(h\right)}+\mathrm{m}\mathrm{o}\mathrm{v}\mathrm{e}\right)\end{array}\right.（8）

式中：l_{i,j}^{\left(h\right)}为拓扑密度变化的下限值；u_i，j^（h） 为变化域的上限值；r_i，j^（h） 表示第j种材料在上一步优化结束后的剩余拓扑密度；move表示每一个优化步中拓扑密度的变化限值，为保证算法稳定性，本文取move=0.2.

本文采用灵敏度过滤法来避免拓扑优化过程出现数值不稳定现象.该方法的核心原理是选取一个合适的过滤半径r_min，并对该半径范围内各单元的应变能梯度进行加权平均处理，从而提高拓扑优化过程的数值稳定性.以第j种材料为例，其加权平均公式可以表示为：

{\left(\frac{\partial \mathrm{S}\mathrm{E}}{\partial {\rho }_{i,j}^{\left(h\right)}}\right)}^{\mathrm{\text{'}}}=\frac{\mathrm{1}}{{\rho }_{i,j}^{\left(h\right)}\sum _{i\in N_{\mathrm{e}}}{H}_{\mathrm{e}i}}\sum _{i\in N_{\mathrm{e}}}{H}_{\mathrm{e}i}{\rho }_{i,j}^{\left(h\right)}\frac{\partial \mathrm{S}\mathrm{E}}{\partial {\rho }_{i,j}^{\left(h\right)}}（9）

式中：N_e表示过滤半径范围内的单元合集.H_ei 表示权重.

1.3 联合优化流程

借鉴联合优化的思路，将形状和拓扑优化方法相结合，自由曲面的形状优化步和拓扑优化步相互交替，以考虑形状和拓扑的相互影响，直至计算收敛.方法流程如图3所示，具体步骤如下.

1）输入曲面的材料属性（包括不同材料的弹性模量、泊松比等）、荷载、边界条件以及初始拓扑密度.

2）根据初始控制点坐标创建NURBS曲面作为自由曲面初始几何形状.

3）进入曲面的形状优化模块.该模块的主要任务是完成曲面形状优化的一个优化步：

①根据曲面划分合适的网格，建立曲面的有限元模型；

②根据梯度法确定的优化方向和步长来更新曲面控制点坐标；

③根据形状优化产生的新的控制点坐标创建新的NURBS曲面作为下一步优化的曲面几何形状模型.

4）进入曲面的拓扑优化模块.拓扑优化模块的主要任务是完成曲面拓扑优化的一个优化步：

①根据形状优化更新的曲面几何形状模型来划分合适的网格，建立曲面的有限元模型；

②采用两种材料SIMP插值模型建立拓扑密度和材料弹性模量之间的相互关系；

③采用交换活性相法和优化准则法来更新材料的拓扑密度.

5）判断此时结构是否符合收敛条件，若是满足则输出优化结果，若不满足则返回第2）步继续计算，直到满足收敛条件，输出优化结果.由于联合优化方法中结构的形状和拓扑都会发生变化，所以在判断收敛时，结构的形状和拓扑需要同时达到各自的收敛条件，即需要同时满足应变能变化率≤0.1%，拓扑密度变化率≤0.1%，联合优化方法才收敛.

1.4 算例验证

为了验证本文提出的自由曲面联合优化方法的有效性，选取文献［17］中的模型进行优化和对比验证.文献［17］在形状优化方面，利用Coons patches插值技术构建曲面，通过调整边界曲线间接改变曲面形态，采用四节点Reissner-Mindlin平面壳单元进行有限元分析.拓扑优化则是将设计补丁材料状态参数化，引入可变密度微观结构多孔材料处理离散问题.联合优化过程中，采用顺序线性规划、单纯形法及最优准则法求解，且形状和拓扑优化交替迭代.

本文算例模型及网格划分和控制点布置如图4所示，模型为一长宽均为6 m的自由曲面壳体，壳体厚度为0.01 m，四角点简支约束，曲面中点施加一个竖直向下的集中荷载F，大小为100 N.考虑到结构的对称性，取整个模型的1/4结构作为分析对象，控制点数为15个，以控制点的竖向坐标作为优化变量.所有控制点的加权系数都设为0.5，在优化过程中能够保持各个控制点对曲面影响的均衡，简化优化过程，减少计算量，以及避免因权重变化带来的曲面不规则性，有助于在保证曲面平滑的同时，集中优化资源在关键的竖向坐标调整上.有限元仍然采用三角形薄壳单元，单元数为1 800，节点数为961.拓扑优化以结构的最小应变能为目标，为与文献［17］一致.这里取两种材料弹性模量比值为1∶0，即进行单一材料拓扑优化.实体材料弹性模量为210 GPa，材料的泊松比为0.3，以实体材料体积占比50%作为体积约束，SIMP插值模型的惩罚因子p=3，过滤半径r_min=0.16 m.

采用本文的联合优化方法对模型进行优化，优化过程如图5所示，整体结构的应变能呈现不断下降趋势，并在第330个优化步收敛，结构的应变能相比初始结构下降了96.5%，实体材料的体积占比为49.5%.从图中可以看出，整个结构起初呈现为相对平坦的曲面状态，随着优化过程推进，呈现为集中荷载位置的曲面高度持续增加，壳体矢高由初始的0.6 m增加到1.56 m，说明曲面在集中力的作用下需要一个较大的曲率以更好地承受外荷载.

将本文优化结果与文献［17］优化结果进行对比，如图6所示.可以看出，两者的优化结果基本一致，都会呈现出集中力荷载位置隆起的形状，曲面矢高分别为1.56 m和1.52 m，仅相差2.6%；拓扑分布也较为相近，实体都分布在角点到中心位置的连线上，与结构主要传力路径一致，但本文模型的“X”形传力路径更为清晰，这是因为本文采用两种材料SIMP模型，结合交换活性相法划分优化域，并利用优化准则法迭代更新密度，有效减少了无效迭代，使材料分布更聚焦于主传力路径.对比结果表明，本文提出的自由曲面结构形状-拓扑联合优化方法更准确且有效.

2 参数影响分析

2.1 分析方案

本文以两种典型的壳体模型为分析对象，即集中荷载作用下四角点简支模型和均布荷载作用下四边简支模型（见图7和图8），通过参数分析研究材料弹性模量比值和体积比值对优化结果的影响.

集中荷载作用下四角点简支模型如图7所示，模型为一长宽均为10 m的平面为矩形的自由曲面壳体，A、B、C、D四角点设置简支支承，曲面中心点E处受竖直向下大小为1×10³ kN的集中力.控制点布置如图7（b）所示，控制点的加权系数设置为0.5，在优化过程中保持不变.曲面的厚度为0.14 m，材料的泊松比设为0.2.有限元计算采用三角形薄壳单元，单元数为3 200，节点数为1 681.SIMP法惩罚因子取为p=3，过滤半径取为r_min=0.25 m.

均布荷载作用下四边简支模型如图8所示，平面尺寸为10 m×10 m，曲面四边设置简支支承，受竖直向下均布力，大小为q=1×10³ kN/m²，控制点布置如图8（b）所示，控制点加权系数设置为0.5.曲面的厚度为0.14 m，材料的泊松比设为0.2.有限元计算采用三角形薄壳单元，经网格敏感度分析后，单元数取为3 200，节点数为1 681.SIMP法惩罚因子取为p=3，过滤半径取为r_min=0.25 m.

为了分析材料弹性模量比值对曲面优化结果的影响，将高弹模材料的弹性模量设定为定值E₁=30 GPa，并调整低弹模材料的弹性模量E₂=0 GPa、3 GPa、 6 GPa、24 GPa，即高弹模材料和低弹模材料的弹性模量比值为E₁∶E₂=1∶0、1∶0.1、1∶0.2、1∶0.8.在分析弹性模量比值的影响时，将材料的体积比值设定为定值，取V₁∶V₂=2∶8.在分析材料体积比值对曲面优化结果的影响时，将高弹模材料体积V₁与低弹模材料体积V₂分别设定为V₁∶V₂=1∶9、2∶8、5∶5、8∶2、9∶1，并将材料的弹性模量比值设为定值，取E₁∶E₂=1∶0.2.

2.2 弹性模量比值的影响

2.2.1 集中荷载作用下四角点简支模型

采用本文的联合优化方法，分别对不同弹性模量比值的情况进行了优化，得到最优形状-拓扑形式，如图9所示，图中黑色代表高弹模材料，白色代表低弹模材料，后文均采用该表示方式.

可以看到，弹性模量比值的变化对曲面的整体形状影响较小，最优形状均为中心点高、四角点低的正高斯曲面.但弹性模量比值的降低会导致曲面最大高斯曲率整体减小（图10），当E₁∶E₂由1∶0降至1∶0.8时，曲面最大高斯曲率从0.15降至0.125，表明低弹模材料弹性模量的增加提升了曲面抵抗外荷载的能力，仅需较小的曲面曲率就可获得最优的受力性能.

由图9的拓扑优化情况可知，当弹性模量比值E₁∶E₂=1∶0时，即仅有一种材料的情况，材料沿中心至角点连线形成连续“X”形传力路径，中心点局部呈网格状.当弹性模量比值为E₁∶E₂=1∶0.1时，高弹模材料主要集中分布在中心点和简支角点处，而并没有在对角线上形成连线.随着弹性模量比值的减小（E₁∶E₂=1∶0.2，E₁∶E₂=1∶0.8），高弹模材料逐渐向集中荷载作用点（中心）和简支角点聚集，拓扑由连续路径转向局部集中.这是因为当E₁∶E₂=1∶0时，结构仅有一种材料，其拓扑需要形成一个连续的传力路径以将集中荷载力传至支座；当E₁∶E₂=1∶0.1、1∶0.2、1∶0.8，曲面由两种材料构成且为连续体，这时较强的材料不需要形成连续拓扑，仅分布在结构主要受力部位即可，随着E₁∶E₂比值的减小，这一趋势变得更为明显.

进一步比较不同弹性模量比值对应的最优曲面结构的位移和Mises应力，如图11所示，其中位移为曲面的竖向位移，且竖直向下为正.可以看出，当弹性模量比值取不同数值时，曲面的最大位移均出现在集中荷载处及其附近区域，但随着弹性模量比值的减小，高位移区域的范围逐渐减小；Mises应力的最大位置均出现在集中荷载以及简支角点处，这也说明了集中荷载和简支角点位置为曲面的主要受力部位.

图12给出了最优结构的最大位移及最大Mises应力随弹性模量比值的变化趋势.随着E₁∶E₂的降低，最大位移下降55.1%.这是由于低弹模材料的弹性模量逐渐增大，曲面整体高度变大，刚度增大.直到当弹性模量比值小于1∶0.2之后，最大位移的变化趋势才变得平稳.另外由于曲面的大致形状以及最大Mises应力出现的位置没有发生变化，弹性模量比值对最大Mises应力的影响比较小，在不同弹性模量比值下结果的波动范围仅为4.5%.

2.2.2 均布荷载作用下四边简支模型

此处给出体积比值取为V₁∶V₂=2∶8时，弹性模量比值的参数分析结果，如图13所示.

在均布荷载作用下的四边简支模型中，弹性模量比值的影响机理与集中荷载作用下四角点简支模型一致（低弹模材料弹性模量的提升会降低曲面曲率）（图14）.但材料分布模式呈现显著差异.由图13的拓扑优化情况可知，当弹性模量比值E₁∶E₂=1∶0时即仅有一种材料的情况，材料以网格状分布覆盖全曲面，形成均匀传力路径.而随着弹性模量比值的减小，高弹模材料越来越向曲面四边的中点位置聚集，表明均布荷载作用下高弹模材料更倾向于强化边界支承区域以分散荷载.

图15给出了最优结构的最大位移及最大Mises应力随弹性模量比值的变化趋势.当E₁∶E₂=1∶0.8时，最大位移较单材料模型降低96.8%，最大Mises应力变化可忽略（波动范围仅为11.3%）.

2.3 材料体积比值的影响

2.3.1 集中荷载作用下四角点简支模型

本节考虑两种材料的不同体积比值对集中荷载作用下四角点简支模型进行了优化分析，最优结果如 图16所示.可以看到，体积比值的变化也未影响最优曲面的整体几何特征，仅是改变了曲面的曲率，体积比值越大，最大高斯曲率越小（图17）.但体积比值的变化对最优拓扑结果有较大影响.

由图16的拓扑优化情况可知，当体积比值比较小时（V₁∶V₂≤2∶8），高弹模材料集中在四个简支角点以及集中荷载位置；而当体积比值逐渐增大（V₁∶V₂≥ 5∶5），材料沿荷载传递路径（中心至角点连线）形成连续“X”形分布，并逐渐填充次要区域；体积比值V₁∶V₂=8∶2和V₁∶V₂=9∶1时，除了曲面四边中点处外均为高弹模材料.

图18给出了最优结构的最大位移及最大Mises应力随体积比值的变化趋势.如图所示，随着体积比值的增加，位移逐渐减小，体积比值增至9∶1时，最大位移降低33.0%，而应力波动不足3.5%，表明高弹模材料体积占比对刚度提升具有主导作用.

2.3.2 均布荷载作用下四边简支模型

此处给出弹性模量比值E₁∶E₂=1∶0.2时，均布荷载作用下四边简支模型不同体积比值的参数分析结果，如图19所示.

体积比值对曲面形态的影响趋势与集中荷载作用下角点简支模型相似，表现为曲面曲率降低以及最大位移减小（图20、图21），且不同体积比值对最大Mises应力结果的影响较小，但材料分布呈现有所不同.

由图19的拓扑优化情况可知，当体积比值V₁∶V₂≤2∶8时，高弹模材料主要分布在四条简支边的中部位置；而随着体积比值的增大，当体积比值V₁∶V₂≥ 5∶5时，高弹模材料覆盖范围扩大，相邻边中点间形成连续连接；在V₁∶V₂=9∶1的优化结果中，除在曲面的四个角点以及曲面中心位置外，均分布有高弹模材料.

3 结 论

本文采用改进的形状-拓扑联合优化方法对不同边界条件、不同荷载形式下的两材料自由曲面模型进行形状-拓扑联合优化研究，系统探讨了材料弹性模量比值和体积比值对优化结果的形状、拓扑及结构力学性能的影响，得出如下结论：

1）材料弹性模量比值对结构的最优形状和拓扑均有影响.对于本文研究的两种模型，随着弹性模量比值E₁∶E₂（E₁>E₂）的减小，最优曲面的整体曲率降低.这是因为低弹模材料刚度（E₂）的提升增强了结构整体刚度，从而使得曲面能够在满足结构力学性能要求的前提下，以更小的曲率实现荷载的高效传递.

2）当E₁∶E₂=1∶0时，即相当于单一材料的情况，对于集中荷载作用下的四角点简支模型，高弹模材料形成“X”形连续传力路径，随着弹性模量比值的降低，高弹模材料逐渐向荷载作用点（壳体中心）和支座（壳体角点）聚集，拓扑趋于局部集中；对于均布荷载作用下的四边简支模型，高弹模材料呈网格状分布以传递均布荷载，随着弹性模量比值的降低，材料向四边中点集中.也就是说，两种材料的材料弹性模量差异越大，高弹模材料越倾向于分布在主应力路径关键区域，而非形成连续结构.

3）两种材料体积比值的变化同样会对曲面的形态特征和拓扑情况产生影响.当V₁∶V₂从1∶9增至9∶1时，集中荷载作用的四角点简支模型及均布荷载作用的四边简支模型曲面最大高斯曲率均逐步下降，表明高弹模材料体积占比增加可显著提升结构刚度，使得曲面优化形态趋于平坦.

4）从最优拓扑结果来看，当V₁∶V₂较小时（如1∶9）集中荷载作用下四角点简支模型的高弹模材料仅分布于荷载点与支座附近，而均布荷载作用下四边简支模型的高弹模材料集中于四边中点；当比值增大至5∶5及以上时，前者的高弹模材料沿荷载传递路径（中心至角点连线）形成连续“X”形分布，后者高弹模材料覆盖范围同样扩大，甚至填充次要受力区域.这表明随着高弹模材料体积占比提高，高弹模材料除了分布在结构的主要受力部位外，还会沿着主要受力部位之间的连接区域进行分布.

5）弹性模量比值和体积比值对最优结构的位移影响显著，对应力影响较小.例如，对于集中荷载作用的四角点简支模型，当E₁∶E₂=1∶0.8时的最大位移较单材料（E₁∶E₂=1∶0）的情况降低约55.1%，而最大Mises应力变化较小（波动范围仅为4.5%），表明刚度提升主要通过优化材料分布而非应力重分布实现.当V₁∶V₂从1∶9增至9∶1时，集中荷载作用的四角点简支模型最大位移降低约33.0%，而最大Mises应力在不同体积比值下的波动不足3.5%，说明通过增加高弹模材料体积占比，可有效提升刚度.

参考文献

[1] 李欣， 武岳， 崔昌禹 ．自由曲面结构形态创建的NURBS-GM方法［J］．土木工程学报，2011，44（10）： 60-66．

[2] LI X， WU Y， CUI C Y. NURBS-GM method for computational morphogenesis of free form structures［J］. China Civil Engineering Journal， 2011，44（10）： 60-66．（in Chinese）

[3] 王磊， 杨彬， 张其林 ．非均匀有理B样条曲面形状优化方法［J］．湖南大学学报（自然科学版），2012，39（7）： 14-19．

[4] WANG L， YANG B， ZHANG Q L. Shape optimization of non-uniform rational B-spline surface［J］. Journal of Hunan University （Natural Sciences）， 2012，39（7）： 14-19．（in Chinese）

[5] SAN B B， FENG D M， QIU Y ．Shape optimization of concrete free-form shells considering material damage［J］. Engineering Optimization，2022，54（12）：1981-1998．

[6] SAN B B， HE H Y， FENG D M， et al. Constrained shape optimization of free-form shells considering material creep［J］. Engineering Optimization，2022，54（10）：1787-1800.

[7] JIANG B S， SU C H， HOU B W ．The influence of invariant horizontal projection area constraint on shape optimization of free-form surface structure with openings［J］. Structures，2024，66：106818．

[8] ZUO W J， SAITOU K ．Multi-material topology optimization using ordered SIMP interpolation［J］．Structural and Multidisciplinary Optimization，2017，55（2）：477-491．

[9] GHABRAIE K ．An improved soft-kill BESO algorithm for optimal distribution of single or multiple material phases［J］．Structural and Multidisciplinary Optimization，2015，52（4）：773-790．

[10] 刘继凯， 张乘虎， 袁志玲，等 ．基于Ordered SIMP插值模型的点阵-实体复合结构拓扑优化设计方法［J］．湖南大学学报（自然科学版）， 2022， 49（2）： 13-19．

[11] LIU J K， ZHANG C H， YUAN Z L，et al ．Topology optimization method for lattice-solid structure design based on Ordered SIMP interpolation［J］. Journal of Hunan University （Natural Sciences），2022， 49（2）： 13-19．（in Chinese）

[12] LI Y， XIE Y M. Evolutionary topology optimization for structures made of multiple materials with different properties in tension and compression［J］. Composite Structures，2021，259：113497.

[13] 刘宏亮， 王成， 梁缘 ．基于离散变量的多材料结构拓扑优化设计［J］．计算力学学报，2024，41（4）：611-617．

[14] LIU H L， WANG C， LIANG Y. Topology optimization design of multi-material structure based on discrete variables［J］. Chinese Journal of Computational Mechanics，2024，41（4）：611-617．（in Chinese）

[15] WANG L， ZHANG Q L， YANG B ．Combined shape and topology optimization of free form shells［J］. Advanced Materials Research，2010，163-167： 2356-2360．

[16] 冯若强，葛金明， 胡理鹏， 等 ．基于B样条曲线的自由曲面索支撑空间网格结构多目标形态优化［J］．土木工程学报， 2015， 48（6）： 17-24．

[17] FENG R Q， GE J M， HU L P，et al. Multi-objective shape optimization of free-form cable-braced grid shells with B-spline method［J］. China Civil Engineering Journal，2015，48（6）：17-24．（in Chinese）

[18] HASSANI B， TAVAKKOLI S M， GHASEMNEJAD H ．Simultaneous shape and topology optimization of shell structures［J］．Structural and Multidisciplinary Optimization，2013，48（1）：221-233．

[19] MENG X C， XIONG Y L， XIE Y M，et al ．Shape-thickness-topology coupled optimization of free-form shells［J］．Automation in Construction， 2022， 142： 104476．

[20] 赵兴忠， 苗闯， 高博青， 等 ．基于鲁棒性的自由曲面结构优化设计研究［J］．建筑结构学报， 2014， 35（6）： 153-158．

[21] ZHAO X Z， MIAO C， GAO B Q，et al ．Optimization design research of free-form structures based on robustness［J］．Journal of Building Structures， 2014， 35（6）： 153-158．（in Chinese）

[22] 高天贺，田阔，黄蕾， 等 ．数据驱动的曲面构件形状-拓扑协同优化方法［J］．航空学报，2024，45（2）：428806．

[23] GAO T H， TIAN K， HUANG L，et al ．Data-driven shape-topology optimization method for curved shells［J］. Acta Aeronautica et Astronautica Sinica， 2024， 45（2）： 428806．（in Chinese）

[24] ANSOLA R， CANALES J， TÁRRAGO J A，et al ．An integrated approach for shape and topology optimization of shell structures［J］. Computers & Structures，2002，80（5/6）：449-458．

[25] ANSOLA R， CANALES J， TÁRRAGO J A，et al. Combined shape and reinforcement layout optimization of shell structures［J］. Structural and Multidisciplinary Optimization，2004，27（4）：219-227．

[26] TAVAKOLI R， MOHSENI S M. Alternating active-phase algorithm for multimaterial topology optimization problems：a 115-line MATLAB implementation［J］. Structural and Multidisciplinary Optimization， 2014， 49（4）： 621-642．

[27] SIGMUND O. Design of material structures using topology optimization［D］. Kongens Lyngby：Technical University of Denmark， 1994.

基金资助

国家自然科学基金资助项目(51978240)

National Natural Science