基于参数化建模和代理模型的大跨储煤结构气动外形优化研究

Study on Aerodynamic Shape Optimization of Large-span Coal-storage Structures Based on Parametric Modeling and Surrogate Models

摘要

为提升大跨储煤网壳结构的抗风性能，基于NURBS理论的参数化构形方法和代理模型技术，提出自由曲面网壳结构的气动外形优化方法. 利用Grasshopper参数化平台和Fortran语言进行程序的二次开发，实现了设计参数调整、模型自动更新和气动响应计算的一体化联动处理，显著提高了优化设计的自动化程度. 以NURBS曲线控制点为优化变量，位移响应极值为优化目标，建立了自由曲面网壳结构优化模型. 在此基础上，对三心圆柱面和球面储煤网壳结构进行气动外形优化. 结果表明：融合参数化建模和Kriging代理模型的气动外形优化方法，通过发掘结构合理的气动外形取得了显著优化效果；0°和30°风向下优化后柱壳结构的位移极值分别下降了32%和18%，球壳的位移极值下降了12%；优化后两种网壳结构的抗风性能均有所改善，且达到降低工程造价的目的，为大跨空间网壳结构的抗风设计选型提供新的途径.

Abstract

To improve the wind resistance performance of large-span coal-storage structures， an aerodynamic shape optimization method for free-form reticulated structures is proposed by utilizing the parametric configuration method based on non-uniform rational B-splines （NURBS） theory and surrogate modeling technique. A secondary development of the program is carried out based on Grasshopper and Fortran to achieve an integrated process for design parameter adjustment， automatic model updating， and aerodynamic response calculation， which significantly enhances the automation level of optimization design. The control points of NURBS curves serve as optimization variables， and the maximum displacement is set as the objective. An optimization model for free-form reticulated structures is established accordingly. On this basis， aerodynamic shape optimizations are performed on tri-cylindrical and spherical coal-storage reticulated structures. The results show that， the proposed method， which combines parametric modeling and a Kriging surrogate model， effectively identifies rational aerodynamic shapes and achieves effective optimization. The maximum displacements of optimal cylindrical shells under wind directions of 0° and 30°are reduced by 32% and 18%， respectively， while the extreme value for spherical shells is reduced by 12%. The wind resistance performance of both types of reticulated shells is improved， and the goal of reducing construction costs is also achieved. Thus， a new approach for wind-resistant design and form selection of large-span reticulated structures can be finally provided.

Graphical abstract

关键词

大跨储煤结构 / 气动外形优化 / 参数化建模 / 非均匀有理B样条 / 代理模型 / 气动响应

Key words

large-span coal-storage structure / aerodynamic shape optimization / parametric modeling / non-uniform rational B-splines （NURBS） / surrogate model / aerodynamic response

随着我国电力事业的蓬勃发展，钢网壳储煤结构的需求日益增加，且结构建设不断向大型化、大跨度方向发展. 近年国内完工的多项大型储煤结构，跨度普遍超过100 m，部分甚至超过200 m，例如国电宁夏方家庄电厂封闭煤棚跨度已达到229 m. 此类结构具有质量轻、刚度小及阻尼小等特点，风荷载已经成为其结构设计的主要控制荷载. 因此，研究和探讨如何改进设计方案，寻找提升结构抗风性能的有效措施变得迫切和重要.

建筑结构抗风优化的主要目的是在满足结构性能要求的前提下，通过优化设计手段，降低设计成本或改善结构抗风性能^［1］. 目前结构抗风优化研究多围绕构件尺寸优化和气动外形优化展开. 针对构件尺寸优化，通常在满足给定的约束条件和结构形状基础上，以总重量或造价为目标函数，对杆件截面或壳体厚度进行优化设计，并在高层建筑^［2-3］、大跨屋盖^［4］和工业厂房^［5］的结构方案设计中应用.

气动外形优化是以发掘结构最合理的气动外形为目标，能够在根源上提升结构的抗风性能. 近几十年来，国内外学者在这方面开展了广泛研究，并建立了一套优化设计方法，主要优化算法包括遗传算法、仿生学法和模拟退火算法等^［6-8］. 但对于多设计变量、多目标优化问题，上述方法仍需进行大规模数值计算来获取全局最优解，在优化效率方面并没有本质突破. Mooneghi等^［9］对建筑结构气动外形优化的发展状况进行了详细的梳理，指出代理模型优化算法具有高效的优化效率和良好的近似全局最优解的能力，已成为结构气动优化研究领域的一个热点. 从现有文献来看，学者们主要针对高层建筑进行抗风优化设计. Bernardini等^［10］和Elshaer等^［11］以高层建筑的平均阻力和脉动升力系数最小为优化目标，采用代理模型与计算流体动力学（computational fluid dynamics，CFD）相结合的优化方法，对高层的角部形状进行了优化，研究指出，代理模型优化算法所需的CFD计算量仅为常规遗传算法的0.75%；Ding等^［12］基于克里金（Kriging）代理模型、雷诺时均湍流模型和大涡模拟技术，提出了多保真的高层建筑气动外形优化设计方法；Jafari等^［13］以高层建筑及其幕墙的阻力系数最小为优化目标，采用响应面（response surface method，RSM）代理模型对幕墙体系的外形轮廓进行了气动优化设计. 针对大跨度屋盖结构，Qiu等^［14］建立了三心圆柱面屋盖几何形状的数学函数模型，以风荷载作用下结构位移响应最小为优化目标，采用Kriging代理模型确定了三心圆截面的最优设计参数.

近年来，随着气动外形优化研究的不断深化，学者们陆续提出将非均匀有理B样条（non-uniform rational B-splines，NURBS）曲线或曲面为基础的形态创建方法应用于气动外形优化中，已成功应用于航空航天、汽车制造等领域. 但在建筑结构领域，将NURBS技术用于结构气动外形优化的研究较少. Wei等^［15］提出一种基于B样条插值算法的结构化网格自适应调整方法，能够有效模拟结构外形变化对其风效应的影响；He等^［16］基于Grasshopper平台提高了复杂外形建筑结构的建模效率，使用非支配遗传算法、Kriging代理模型和CFD技术对螺旋形超高层建筑进行气动外形优化；Ye等^［17］采用NURBS曲线对高层建筑架空层的外形和支撑构件截面进行了气动优化设计，可为该层风力发电机提供高风速、低湍流度的风力条件；Ogawa等^［18］和Yamashita等^［19］基于NURBS曲线创建了二维柱面网壳结构模型，以强风作用下结构位移响应最小为优化目标，利用多项式响应面代理模型对柱面网壳进行了自由曲面优化设计.

从工程应用角度来看，形态自由化已被证实为空间结构发展的主要趋势之一，同时NURBS曲线和曲面已成为结构设计领域的核心技术^［20-21］. 对于大跨度储煤结构，常见的曲面形状包括圆柱面、三心圆柱面和球面^［22］. 而基于NURBS法的曲面造型能够提供更灵活的形状定义和多样化的结构外形，这有助于通过气动外形优化进一步提升大跨储煤结构的抗风性能.

基于上述研究基础，本文借助Grasshopper参数化设计平台和Fortran程序设计语言，提出融合NURBS建模技术、代理模型优化算法、CFD数值模拟与有限元分析的自由曲面网壳结构气动外形优化方法，并通过三心圆柱面和球面两种典型储煤结构算例验证其优化效果，以期为大跨空间结构的抗风设计选型提供参考.

1 自由曲面网壳结构气动外形优化方法

1.1 基于NURBS方法的自由曲面参数化建模

自由曲面网壳结构气动外形优化的前提是对曲面进行参数化描述. NURBS方法，即非均匀有理B样条方法，通过调节若干个参数实现对曲面形状的控制，同时为标准的解析形状以及自由型曲线和曲面的精确表示与设计提供通用数学公式，从而提高气动外形优化设计的灵活性^［23］. 因此，本文基于NURBS方法生成储煤结构的曲面形状.

由分段有理多项式定义的NURBS曲线方程，按式（1）确定.

\mathit{C}(v)={\displaystyle \sum _{i=\mathrm{0}}^n}{R}_{i,p}(v){\mathit{P}}_i（1）

式中： C （v）为NURBS曲线上的点， P_i 为曲线特征控制点，两者的区别在于控制点的个数是有限的，而曲线上点的个数是无限的；R_i，p （v）为有理基函数，其计算式如下：

R_{i,p}(v)=\frac{N_{i,p}(v)w_i}{{\displaystyle \sum _{j=\mathrm{0}}^n}{N}_{j,p}(v)w_j}（2）

式中：N_i，p （v）和N_j，p （v）为p阶B样条基函数；w_i 和w_j 为相应控制点的权因子；n为V方向控制点总数；i和j为V方向的索引；v为自由曲线沿V方向的节点矢量，v∈［0， 1］. 以上参数定义如图1所示.

在NURBS曲线中，利用控制点位置和控制点权因子等参数可方便地对曲线形状进行调节. 为实现气动外形优化中对自由曲面的参数化构建，本文利用可视化节点式编程工具，即Rhino软件的Grasshopper参数化设计平台，对自由曲面网壳结构进行参数化构形. 首先，基于NURBS曲线构造一条母线和一条准线，以准线为轨道曲线将母线扫掠延伸，即可形成所需的曲面造型（图1），或者将母线绕Z轴单摆旋转360°，生成闭合的自由曲面；其次，采用等参线分割法将曲面划分为多个四边形网格^［24］，以此为基础，生成风致响应计算所需的单层或双层空间网壳结构.

1.2 气动外形优化模型

将气动外形优化问题转化为数学模型来求解，通常包括优化变量、目标函数和约束条件. 由于控制点是对曲面形状影响较大的参数，因此选取控制点坐标向量作为优化变量，自由曲面网壳结构的气动外形优化模型可表示为：

\left\{\begin{array}{l}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d}x=\left[x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n\right]\\ \mathrm{m}\mathrm{i}\mathrm{n}f(x)\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}:x_{\mathrm{m}\mathrm{i}\mathrm{n}}\le x\le x_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}\right.（3）

式中：x为外形优化变量的集合，n为优化变量个数； f（x）为优化目标函数，一般选取与结构气动性能相关的指标作为优化目标，如结构风致响应极值、整体风力系数和极值风速比等；x_max和x_min分别表示优化变量取值范围对应的上、下限值，一般根据设计要求或现行规范等确定.

1.3 Kriging代理模型

代理模型是指用近似函数表示优化变量与目标函数值之间的关系，在气动外形优化中，可解决直接调用高可信度分析模型（CFD计算）所面临的优化效率过低的问题^［9，11］. 常用的代理模型有响应面（RSM）、Kriging、多项式和神经网络等. 其中，Kriging模型具有较好的预测非线性、多峰值函数的能力，较好地折中了优化效率和近似精度，具有很好的应用前景^［25］. 因此，本文选用Kriging模型建立目标函数和优化变量的近似数学模型.

Kriging模型是指在m维向量 X =［ x{}_{\mathrm{1}}^{(\mathrm{1})}， x{}_{\mathrm{2}}^{(\mathrm{2})}，…， x{}_m^{(n)}］^T上观测到的n个样本点的真实目标函数值，对其进行线性寻优、无偏内插估计的一种方法，未知函数 y（ X ）的数学表达式为：

y(\mathit{X})=\beta +Z(\mathit{X})（4）

式中：β为未知常数，是随机过程y（ X ）的均值；Z（ X ）假定为平稳高斯随机过程，其协方差矩阵为：

\mathrm{C}\mathrm{o}\mathrm{v}\left[Z({\mathit{X}}^{\left(i\right)}),Z({\mathit{X}}^{\left(j\right)})\right]={\sigma }^{\mathrm{2}}\mathit{R}\left[R({\mathit{X}}^{\left(i\right)},{\mathit{X}}^{\left(j\right)})\right]（5）

R({\mathit{X}}^{\left(i\right)},{\mathit{X}}^{\left(j\right)})=\mathrm{e}\mathrm{x}\mathrm{p}\left(-{\displaystyle \sum _{k=\mathrm{1}}^m}{\theta }_k{‖x_k^{\left(i\right)},x_k^{\left(j\right)}‖}^{\gamma }\right)（6）

式中：σ²为高斯随机过程的方差； R 为n×n阶相关矩阵；R（ X^（i）， X^（j））为任意两个样本点 X^（i）和 X^（j）的相关函数，i， j∈［1，n］；γ为指数，取值范围为0~2；q_k 为相关函数中未知向量 θ =［θ₁，θ₂， …， θ_m ］^T的第k个分量，该向量控制着最终代理模型的光滑性.

上述公式中存在三个待定参数β、σ²和 θ，其中β和σ²由最小二乘法估计得到：

\left\{\begin{array}{l}\widehat{\beta }=({\mathit{I}}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}{\mathit{I})}^{-\mathrm{1}}({\mathit{I}}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\mathit{Y})\\ {\widehat{\sigma }}^{\mathrm{2}}=\frac{\mathrm{1}}{n}{(\mathit{Y}-\mathit{I}\widehat{\beta })}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}(\mathit{Y}-\mathit{I}\widehat{\beta })\end{array}\right.（7）

式中： I 为m维单位向量； Y =［y（ X^（1）），y（ X^（2）），…， y（ X^（n））］. 待定参数 θ 则通过最大化似然函数来进行估计：

\begin{array}{l}L\left(\mathit{\theta }\right)=-\frac{\mathrm{1}}{\mathrm{2}}\left[n\mathrm{l}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }\right)+n\mathrm{l}\mathrm{n}{\sigma }^{\mathrm{2}}+\mathrm{l}\mathrm{n}\left|\mathit{R}\right|+\right.\\ \left.\frac{\mathrm{1}}{{\sigma }^{\mathrm{2}}}{\left(\mathit{Y}-\mathit{I}\beta \right)}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}\left(\mathit{Y}-\mathit{I}\beta \right)\right]\end{array}（8）

本文采用遗传算法搜索出最优 θ （θ_k ）值，使得似然函数值最大，确保所构造的代理模型具有最高精度. 综上，利用Kriging模型可获得未知样本点 X 对应的目标函数预测值，其计算公式为：

\widehat{y}(\mathit{X})=\widehat{\beta }+\mathit{r}{(\mathit{X})}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}(\mathit{Y}-\mathit{I}\widehat{\beta })（9）

式中： r （ X ）为n维相关向量， r （ X ）=［R（ X， X^（1））， R（ X， X^（2）），…， R（ X， X^（n））］^T.

Kriging模型本质上是一个插值函数模型，为提高模型的预测精度，需要对样本点进行动态更新，本文采用改善期望最大准则，即EI准则，来确定新增样本点. 对于最小化问题，EI公式如下：

E[I(x)]=(y_{\mathrm{m}\mathrm{i}\mathrm{n}}-\widehat{y})\Psi \left(\frac{y_{\mathrm{m}\mathrm{i}\mathrm{n}}-\widehat{y}}{\widehat{s}}\right)+

\widehat{s}\Phi \left(\frac{y_{\mathrm{m}\mathrm{i}\mathrm{n}}-\widehat{y}}{\widehat{s}}\right)（10）

式中：y_min为当前样本点确定的目标函数最小值； Ψ（·）和Φ（·）分别为标准正态分布函数及其概率密度函数；ŝ为均方根误差. 在实际应用中，将EI最大值处的点作为新增样本点添加至代理模型，可有效改善其预测精度.

1.4 基于代理模型的优化设计流程

针对不同样本点构造相应的自由曲面，并通过气动响应分析获得目标函数响应值，是开展基于代理模型气动外形优化设计的基础和前提. 为提高优化计算效率，本文搭建了参数化建模、流域网格划分、流体计算和结构有限元分析的一体化气动响应分析程序流程，如图2所示.

图2所示流程图中，建立基于NURBS曲线的自由曲面屋盖及其网壳模型均由Grasshopper平台二次开发完成；再采用Fortran语言，结合ICEM网格划分软件、Fluent流体计算软件和Ansys有限元软件的参数化程序设计语言（JRF、Journal和APDL等）进行二次开发，实现参数修改、信息输入、跨软件通信及计算结果返回等功能一体化联动处理. 该流程可实现Grasshopper模型变更时，流体计算模型和有限元模型即时动态更新，由计算实时返回目标函数计算结果.

基于Kriging代理模型的气动外形优化策略如 图3所示，具体步骤如下：

1）代理模型的构建依赖合理有效的设计空间抽样策略，即试验设计方法. 本文采用最优拉丁超立方法^［26］，从初始设计空间中选取初始样本点，保证样本点的填充均匀性.

2）利用一体化气动响应分析程序计算初始样本点的真实目标函数值，并将这些样本点和目标函数值保存至样本数据库中.

3）利用1.3节的数学方法构造初始Kriging代理模型.

4）根据EI准则，利用遗传算法对代理模型进行全局寻优，获得可能最优点及其近似最优目标函数值.

5）采用气动响应分析计算最优点所对应的真实目标函数值.

6）根据给定的收敛标准ε（ε=10^-5），判断优化是否达到收敛条件，计算式为：

（11）

式中：y_max为当前样本点确定的目标函数最大值. 若是，停止优化，否则将所得的最优点作为新增样本点添加至现有样本数据库，回到步骤3）重新构建代理模型，直至整个优化过程收敛.

2 三心圆柱面与球面储煤结构气动外形优化算例

2.1 优化对象

火电厂储煤结构要求跨度大、净空高，满足储煤与作业空间，柱面与球面网壳已成为储煤棚的主要结构形式. 本文以双层柱面网壳和单层球面网壳结构为优化对象，结构模型如图4所示. 柱面储煤棚采用两段封闭的正放四角锥三心圆柱面双层网壳结构，结构初始设计参数如下：网壳跨度L和纵向长度B均为108 m，矢高H=40.26 m，网格尺寸和网壳厚度均为4 m；截面主圆弧半径R=70.12 m，圆弧夹角θ=61.32°，小圆弧半径及其夹角分别为R₁=37.46 m和θ₁=54.34°. 上弦杆、下弦杆和腹杆截面分别选用ϕ168 mm×6 mm、ϕ159 mm×6 mm和ϕ89 mm×4 mm钢管，采用上弦节点的支承方式.

图4（b）为K8型单层球面网壳结构，跨度L= 60 m，矢跨比H/L=1/3，网格尺寸为3 m，周边三向铰接于6 m高的承台上，顶部节点标高为26 m. 肋杆截面选用ϕ114 mm×6 mm，纬杆和斜杆采用ϕ100 mm×6 mm型号钢管. 两种网壳模型的结构材料均为Q235B钢，采用Beam188梁单元模拟杆件，焊接球节点. 经校核，网壳结构的全部构件及节点设计均满足《空间网格结构技术规程》（JGJ 7—2010）^［27］要求.

2.2 数值模拟方法验证

准确模拟储煤结构表面的风荷载特性，是获得真实可靠目标函数值的必要条件. 图4所示柱面和球面网壳系列风洞试验在哈尔滨工业大学风洞与浪槽联合实验室中进行^{［14，22］}，本节将试验结果用以验证数值模拟结果的可信度.

计算域及网格划分如图5所示. 两种网壳模型的流场计算域尺寸取18L×10L×5H，模型中心置于计算域前1/3处，阻塞率均小于3%. 将计算域分为内外两部分：1）内计算域选用具有良好适应性的四面体非结构网格，对模型表面及其附近网格进行加密，通过网格独立性分析，最终确定最小网格高度为H/50，无量纲高度y⁺≤90；2）在远离模型的外围空间选用规则拓扑的六面体结构网格. 采用上述混合网格的离散方式，在外形优化过程中只需要重新划分内部流域网格即可，这样有利于提高网格划分效率以及模拟精度. 计算模型的总网格数控制在350万以内.

根据《建筑结构荷载规范》（GB 50009—2012）^［28］以及风洞模拟的风场，选用B类地貌，地面粗糙度指数α=0.15. 图6给出目标风场中平均风速和湍流度剖面的试验测量值、数值模拟值与规范理论值的对比情况. 试验平均风速为10 m/s，取模型跨度作为特征尺寸，计算得到柱面和球面网壳模型的雷诺数分别为3.73×10⁵和4.14×10⁵. 来流湍流特性通过给定的湍动能k、湍流耗散率ε或比耗散率ω公式计算：

k=\mathrm{1.5}(U{I}_{\mathrm{u}}{)}^{\mathrm{2}}（12）

\left\{\begin{array}{l}\epsilon =C_{\mathrm{\mu }}^{\mathrm{0.75}}{k}^{\mathrm{1.5}}/l\\ \omega =k^{\mathrm{0.5}}/C_{\mathrm{\mu }}^{\mathrm{0.25}}l\end{array}\right.（13）

式中：U为入口面的平均风速；I_u为湍流强度，模型顶点处的湍流度约为12%；C_μ为经验常数，常取0.09； l为湍流特征尺度.

流域顶面及两侧采用对称边界，模型表面及地面采用无滑移壁面条件，出口面采用完全发展边界条件. 压力和速度场的耦合采用SIMPLEC（semi-implicit method for pressure-linked equations-consistent）算法，对流项的空间离散采用二阶迎风格式，选用非平衡壁面函数来模拟壁面附近的流场. 当所有控制方程相对迭代残差降至10^-4以下，且模型表面风压基本不变时，判定流场达到稳态.

分别采用标准k-ε、RNG k-ε、SST k-ω和RSM湍流模型进行数值模拟. 将三心圆柱面与球面模型中心线条带的平均风压系数与风洞试验数据进行对比，如图7所示. 对于柱面屋盖，在本文特定的网格和设置参数下，除标准k-ε模型之外，其余三个湍流模型均有良好的模拟效果；RSM模型在球壳顶部及尾流区的模拟结果与试验数据吻合较好，整体误差控制在±12%以内，明显优于其他湍流模型. 综上，本文选用RSM模型作为储煤结构风场模拟的湍流模型.

2.3 优化数学模型

为寻求给定条件下抗风性能最佳的曲面造型，利用1.1节的NURBS方法和Grasshopper平台分别对三心圆柱面和球面两种典型储煤结构进行参数化建模. 表征自由曲面几何特征的轮廓线及其控制点，如图8所示. 两种储煤结构具有对称的曲面形状，控制点A和D分别决定了结构跨度和高度，B、C两个控制点则决定了曲面弧度. 图8所示的初始NURBS曲线对应柱面与球面储煤结构网壳的初始形状（如图4）. 基于控制点调整法，本文将其他参数固定，选取竖向坐标z_B 和水平向坐标x_C 作为形状参数优化变量，以实现自由曲面弧度的调节.

大跨网壳结构在强风作用下的位移响应是结构设计的主要控制目标之一，故本文以网壳结构位移响应为目标函数^{［14，18］}，优化问题可以描述为：

\left\{\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}f(z_B/H,x_C/L)=|w{|}_{\mathrm{m}\mathrm{a}\mathrm{x}}\\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}:{\overline{z}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le z_B/H\le {\overline{z}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\\ {\overline{x}}_{\mathrm{m}\mathrm{i}\mathrm{n}}\le x_C/L\le {\overline{x}}_{\mathrm{m}\mathrm{a}\mathrm{x}}\end{array}\right.（14）

式中：|w|_max为所有网壳节点的位移响应极值；z_B /H和x_C /L分别为无量纲竖向坐标和水平坐标调整系数，其中H为网壳矢高，L为网壳跨度. 综合考虑煤堆边缘堆高、斗轮机作业空间及安全距离等要求^［29］，最终确定三心圆柱壳优化变量取值范围为：z_B /H=［0.2， 0.8］，x_C /L=［0.35， 0.65］；球壳优化变量取值范围为：z_B /H=［0， 0.9］，x_C /L=［0， 0.4］.

为真实反映储煤结构外形及其刚度变化对气动特性的影响，在优化过程中通过调整网格节点数量，保证网格尺寸和网壳厚度与初始网壳参数取值相同. 取原型10 m高度处平均设计风速为25 m/s，采用本文方法对柱壳和球壳储煤结构进行气动外形优化，其中柱壳考虑典型风向0°和最不利风向30°两个风向角工况^［14］，每个工况重复3次优化以确保获得全局最优解，并利用设计空间内的所有样本点构建目标函数的响应曲面.

3 优化结果分析

3.1 柱面储煤结构

基于Kriging代理模型构建的柱面储煤结构形状参数与目标函数的响应曲面如图9所示. 由图可知，设计变量与目标函数之间呈现强烈的非线性关系，在容许挠度范围内存在多种外形可供选择. 其中，OS代表初始网壳模型，A1（B1）表示最优解，A2（B2）和A3（B3）作为结构气动性能对比分析的参考设计外形. 针对这几组形状设计参数，图10给出了优化前后储煤结构的形状对比.

由图9（a）和图10（a）可知，在0°风向角下初始三心圆柱面网壳的位移极值|w|_max为0.086 m. 具有最佳外形的网壳结构设计变量为z_B /H=0.24，x_C /L=0.39，两侧曲线相较于初始形状略呈下凹趋势，结构的最大位移为0.058 m，较初始结构下降了33%. 但随着NURBS曲线控制参数z_B /H和x_C /L的增大，例如参考点A2和A3，截面轮廓线逐渐扩张，结构位移极值会显著增大，|w|_max分别为0.141 m和0.185 m. 因此，0°风向角工况下的优化设计主要是利用NURBS方法对三心圆曲线弧度进行细部调整，使其更加圆滑，最终趋向于抗风性能更合理的气动外形.

对于30°斜风向优化工况［图9（b）］，寻优区内目标函数响应曲面的变化规律与0°风向角工况相似，仅波动幅度略有减小. 初始柱面网壳的位移极值为0.08 m，寻优最终所得的结构最大位移为0.066 m，优化后降低约18%. 由图10（b）可知，与初始结构相比，最佳外形网壳的曲线弧度明显减小，与圆柱面较为接近. 比较参考点B2和B3可知（|w|_max分别为0.111 m和0.177 m），结构外形会随着NURBS曲线参数z_B /H和x_C /L的增加逐渐向两侧扩张，但结构位移极值会显著增大.

为进一步分析优化前后柱面网壳的气动特性，图11、图12分别给出了结构表面平均风压系数和节点位移云图. 由图11、图12可知风向角对柱壳结构风荷载和风致响应的影响较大. 由图11可知，0°风向角下，迎风面靠近底部的区域为风压区，最大负压区出现在顶部区域；与初始形状相比，具有最佳外形的柱壳顶部风吸力有所减弱；随着柱壳曲线弧度的增加（A2和A3），顶部负压区的范围逐渐扩大，且最大风吸力由中间移向两侧外凸区域. 30°风向角下，最大风压力出现在柱壳迎风端部的上游区，风吸力极值则出现在端部的顶点及其下游区域；对比初始点OS和最优点B1可知，最优柱壳沿跨向的负压区及其最大风吸力呈减小趋势.

图12为柱面网壳结构的节点位移云图. 当风向角为0°时，优化前后柱壳的整体位移分布规律没有明显改变，最大位移始终出现在顶部，但随着柱壳两侧曲线曲率的减小，位移极值逐渐减小；30°斜风向下网壳的整体位移分布规律与0°风向工况基本一致，最大位移位于结构顶部. 该风向工况下，最优结构外形与1/3矢跨比圆柱面较为相近，其气动原理是通过降低结构顶部区域的风吸力值来控制其位移响应.

表1给出了不同形状设计参数下柱面网壳的其他气动参数统计结果. 由表可知，风向角为0°时，具有最佳外形柱壳结构的升力系数、阻力系数和极值内力均有下降，较初始结构分别下降了17%、32%和8%；当风向角为30°时，最优柱壳结构的风力系数和杆件内力略有降低，升力系数、阻力系数和极值内力分别较初始结构降低了9%、6%和11%.

上述优化结果表明，本文气动外形优化方法能够较好地改善柱面网壳的抗风性能. 由于结构总重量或总造价是结构设计中需要考虑的一种重要因素，图13给出了储煤结构总用钢量随优化变量的变化规律. 由图13可知，总用钢量随着z_B /H和x_C /L的减小逐渐降低. 0°和30°风向角下最优结构形状均趋向于类圆柱面，总用钢量分别为512 t和503 t，较初始结构（536 t）分别降低了4.5%和6.2%. 因此，综合考虑工程造价以及风向角的影响，建议采用NURBS曲线竖向控制点z_B /H变化范围为［0.2， 0.25］，水平控制点x_C /L变化范围为［0.35， 0.4］的结构设计方案，其带来的结构性能收益较为可观.

3.2 球面储煤结构

进一步选取常见的球面储煤结构为优化算例，通过代理模型构建的设计变量与目标函数响应曲面关系如图14所示. 图15对比了球面储煤结构的初始形状（OS）、最优形状（C1）和两组参照形状（C2和C3）. 由图14可知，结构初始形状为1/3矢跨比的球面网壳，节点位移极值为0.017 m，由优化锁定的最优解为|w|_max =0.015 m，较初始结构下降了约12%. 此外，最优结构的几何外形与球面较为接近，只是随着NURBS曲线参数的微调，其圆弧弧度略有减小.

对照参考形状C2和C3可知，随着NURBS曲线控制参数z_B /H和x_C /L的增大，结构两侧轮廓线逐渐向外扩张，此时风荷载作用下结构的位移响应会显著增大，参考形状C2对应的位移极值为0.058 m；反之，当NURBS曲线控制参数减小时，网壳结构的几何外形呈近似圆锥体（C3），节点的位移极值为0.028 m，但相对初始形状和最优形状分别增加了65%和87%.

图16给出了所选四种几何特征球面网壳0°风向角下平均风压云图和节点位移云图. 由图16（a）可知，对于OS和C1两组网壳，结构迎风面均表现为正压，最大风吸力出现在顶部区域，与初始球壳相较，具有最优形状的网壳顶部风吸力峰值明显减小，且尾流区开始出现正压，主要由气流在尾流区的再附引起；对比C2和C3网壳可知，当轮廓线逐渐向两侧上拱时，结构顶部的风吸区范围有增大趋势，而当网壳演变成近似圆锥体形状时，其顶部风吸力峰值将显著增大.

由图16（b）可知，初始球壳和具有最优形状网壳的节点位移具有相似的分布规律，只是节点位移极值从顶部附近沿着横风向逐渐移向两侧；对于上拱的C2网壳，位移极值相对集中在结构顶部区域，但其幅值会迅速增大；对于类锥形的C3网壳结构，节点位移分布较为对称，最大位移出现在横风向两侧区域.

表2对比了不同几何特征球面网壳的其他气动参数. 与初始球壳相比，具有最佳外形网壳结构的升力系数、阻力系数和极值内力均有降低，较初始结构分别降低了16%、20%和13%，表明本文形状优化方法能够改善结构的抗风性能. 另外两组形状参数下的球壳结构在风荷载作用下的受力特征相对较差，位移极值和杆件内力均有大幅增加.

球壳储煤结构总用钢量随形状设计参数的变化规律如图17所示. 与柱壳类似，随着NURBS曲线设计参数值的减小，结构总用钢量逐渐降低. 最优网壳的用钢量为251 t，较初始结构（268 t）降低了6.3%，且内部空间满足储煤和机械作业空间等要求.

4 结 论

1）针对自由曲面网壳结构气动外形优化问题，将参数化建模和代理模型技术引入结构气动外形优化领域，提出了融合NURBS曲面造型方法和Kriging代理模型的自由曲面网壳结构气动外形优化方法. 基于Grasshopper平台和Fortran语言，对相关软件进行二次开发，实现了参数化建模、流域网格划分、风场数值模拟和结构有限元分析的一体化联动，可显著提高优化设计的自动化程度.

2）采用本文优化方法进行三心圆柱面储煤结构气动外形优化. 通过代理模型技术构建NURBS形状控制参数与目标函数的响应曲面，从中选取最合适的解，保证网壳节点位移极值和其他气动参数指标均有所下降；由0°和30°风向工况可知，三心圆柱壳的最优形状均趋向类圆柱面结构，故在满足净高和工艺要求的前提下，建议结构设计时尽量减小落地圆弧的曲率，能够同时降低工程造价.

3）进一步开展球面储煤结构气动外形优化，其优化效果显著. 结果表明，以NURBS形状参数为优化设计变量，在保证网壳结构整体光顺的同时，发掘了合理的气动外形，并有效控制了结构的风致位移响应. 但考虑到单目标函数无法全面反映自由曲面网壳的气动性能，后续研究将以风致动力响应或稳定性为控制目标开展多目标优化，并进一步考虑结构大开口自然通风的影响，为大跨自由曲面网壳结构的抗风设计选型提供新思路.

参考文献

[1] 傅继阳，吴玖荣，徐安 ．高层建筑抗风优化设计和风振控制相关问题研究［J］．工程力学， 2022，3 9（5）：13-33．

[2] FU J Y， WU J R， XU A. Some issues on wind resistant optimization design and on wind-induced vibration control of tall buildings［J］. Engineering Mechanics， 2022， 39（5）： 13-33．（in Chinese）

[3] 黄铭枫 ．基于风振性能的高层建筑抗风设计优化［J］．工程力学，2013，30（2）： 240-246．

[4] HUANG M F. Wind-induced vibration performance-based optimal structural design of tall buildings［J］. Engineering Mechanics，2013， 30（2）： 240-246．（in Chinese）

[5] DENG T， FU J Y， ZHENG Q X，et al ．Performance-based wind-resistant optimization design for tall building structures［J］．Journal of Structural Engineering，2019，145（10）：04019103．

[6] 黄友钦，岳启哲，傅继阳，等 ．基于应力比法的双层网壳结构抗风优化［J］．结构工程师，2013， 29（4）：112-117．

[7] HUANG Y Q， YUE Q Z， FU J Y，et al. Wind resistant optimization of double-layer reticulated shells by the stress ratio method［J］. Structural Engineers， 2013， 29（4）： 112-117．（in Chinese）

[8] 吴玖荣， 董超超， 徐安， 等. 大跨度变截面门式刚架抗风优化设计研究［J］. 武汉理工大学学报， 2012， 34（2）： 69-74．

[9] WU J R， DONG C C， XU A，et al. Wind-resistant optimized design of long span portal-rigid frame with web-tapered section［J］. Journal of Wuhan University of Technology， 2012， 34（2）：69-74．（in Chinese）

[10] WANG Z Y， MULYANTO J A， ZHENG C R，et al ．Research on a surrogate model updating-based efficient multi-objective optimization framework for supertall buildings［J］．Journal of Building Engineering，2023，72：106702．

[11] TIWARI P K， VIDYARTHI D P ．Improved auto control ant colony optimization using lazy ant approach for grid scheduling problem［J］．Future Generation Computer Systems，2016，60：78-89．

[12] ABDEL-BASSET M， DING W P， EL-SHAHAT D ．A hybrid Harris Hawks optimization algorithm with simulated annealing for feature selection［J］．Artificial Intelligence Review，2021，54（1）：593-637．

[13] MOONEGHI M A， KARGARMOAKHAR R. Aerodynamic mitigation and shape optimization of buildings：review［J］. Journal of Building Engineering，2016，6：225-235．

[14] BERNARDINI E， SPENCE S M J， WEI D，et al ．Aerodynamic shape optimization of civil structures：a CFD-enabled Kriging-based approach［J］．Journal of Wind Engineering and Industrial Aerodynamics，2015，144：154-164．

[15] ELSHAER A， BITSUAMLAK G， DAMATTY A EL ．Enhancing wind performance of tall buildings using corner aerodynamic optimization［J］．Engineering Structures，2017，136：133-148．

[16] DING F， KAREEM A. A multi-fidelity shape optimization via surrogate modeling for civil structures［J］. Journal of Wind Engineering and Industrial Aerodynamics，2018，178：49-56．

[17] JAFARI M， ALIPOUR A. Aerodynamic shape optimization of rectangular and elliptical double-skin façades to mitigate wind-induced effects on tall buildings［J］．Journal of Wind Engineering and Industrial Aerodynamics， 2021， 213：104586．

[18] QIU Y， YU R， SAN B B，et al ．Aerodynamic shape optimization of large-span coal sheds for wind-induced effect mitigation using surrogate models［J］. Engineering Structures，2022，253：113818．

[19] WEI D， SPENCE S M J， KAREEM A，et al ．A structured mesh boundary motion approach for simulating wind effects on bluff bodies with changing boundaries［J］. Journal of Wind Engineering and Industrial Aerodynamics，2014，126：118-131．

[20] HE Z， LIANG T， LAI X，et al ．Vibration acceleration-integrated parameterized aerodynamic shape optimization of super high-rise buildings with spiral configurations［J］. Structural and Multidisciplinary Optimization，2022，65（10）：291．

[21] YE X L， ZHANG X L， WEERASURIYA A U，et al ．Optimum design parameters for a venturi-shaped roof to maximize the performance of building-integrated wind turbines［J］．Applied Energy，2024，355：122311．

[22] OGAWA T， OHSAKI M， MIYAMURA T， et al. Shape optimization of shell roofs subjected to strong wind by using a variable complexity model［J］. Journal of the International Association for Shell and Spatial Structures， 2005， 46（2）：108-115.

[23] YAMASHITA T， KUMAGAI T， OGAWA T， et al. Shape optimization of long-span translational free-form shell roofs in strong wind using multigrid method and variable complexity model［C］//2010 JMSE Proceedings of the 6th China-Japan-Korea Joint Symposium on Optimization of Structural and Mechanical Systems. Kyoto， Japan. JMSE， 2010.

[24] 石开荣， 林金龙， 姜正荣， 等 ．基于模拟植物生长算法的自由曲面空间网格结构形状优化方法［J］．建筑结构学报， 2022，43（12）： 276-289．

[25] SHI K R， LIN J L， JIANG Z R，et al ．Shape optimization method of free-form space grid structures based on plant growth simulation algorithm［J］．Journal of Building Structures， 2022， 43（12）： 276-289．（in Chinese）

[26] 王志成，曹正罡，赵林，等 ．基于混合算法的自由曲面网格结构多目标优化［J］．哈尔滨工业大学学报，2022，54（4）：111-117．

[27] WANG Z C， CAO Z G， ZHAO L，et al. Multi-objective optimization of free-form grid structures based on hybrid algorithm［J］. Journal of Harbin Institute of Technology，2022，54（4）：111-117．（in Chinese）

[28] 苏宁，彭士涛，孙瑛，等 ．大跨度煤棚主体结构设计风荷载的神经网络建模研究及应用［J］．建筑结构学报，2019，40（7）：34-41．

[29] SU N， PENG S T， SUN Y，et al ．Research and application on artificial neural network modeling of design wind load on large-span coal storage sheds［J］．Journal of Building Structures，2019，40（7）： 34-41．（in Chinese）

[30] 卢旦，李承铭，汪大绥， 等 ．自由曲面建筑一体化造型与优化设计研究［J］．建筑结构学报， 2010， 31（5）： 55-60．

[31] LU D， LI C M， WANG D S， et al ．Integrated modeling and shape optimization for free form surface architectures［J］．Journal of Building Structures， 2010， 31（5）： 55-60．（in Chinese）

[32] 陆金钰，姜昊瑞，廖杰 ．波浪形自由曲面网壳参数化建模及稳定性分析［J］．空间结构，2019，25（1）： 14-22．

[33] LU J Y， JIANG H R， LIAO J ．Parametric modeling and stability analysis of wave-shaped free-form grid shells［J］. Spatial Structures，2019，25（1）：14-22．（in Chinese）

[34] 张建侠，马义中，张延静，等 ．基于Kriging模型和两目标约束应对策略的代理优化算法［J］．计算机集成制造系统， 2018，24（12）：3001-3007．

[35] ZHANG J X， MA Y Z， ZHANG Y J，et al. Surrogate-based optimization algorithm based on Kriging models and bi-objective constraint-handling strategy［J］.Computer Integrated Manufacturing Systems， 2018， 24（12）： 3001-3007．（in Chinese）

[36] JIN R C， CHEN W， SUDJIANTO A ．An efficient algorithm for constructing optimal design of computer experiments［J］．Journal of Statistical Planning and Inference，2005，134（1）： 268-287．

[37] 空间网格结构技术规程：JGJ 7—2010 ［S］．北京：中国建筑工业出版社，2011．

[38] Technical specification for space frame structures：JGJ 7—2010 ［S］.Beijing：China Architecture & Building Press，2011．（in Chinese）

[39] 建筑结构荷载规范：GB 50009—2012 ［S］．北京： 中国建筑工业出版社，2012．

[40] Load code for the design of building structures：GB 50009—2012 ［S］. Beijing：China Architecture & Building Press，2012．（in Chinese）

[41] 杨勇，刘观发，方国华，等 ．四连跨干煤棚三心圆柱面网壳结构设计与研究［J］．建筑结构，2017，47（15）：6-10．

[42] YANG Y， LIU G F， FANG G H，et al ．Design and research on three-centered cylindrical reticulated shell structure of four-span dry-coal shed［J］．Building Structure，2017，47（15）：6-10．（in Chinese）

基金资助

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

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

National Natural Science