流线型箱梁气动外形的数值优化

2021-09-26周志勇胡传新黄德睦

哈尔滨工业大学学报 2021年10期

孙强，周志勇，胡传新，秦鹏，黄德睦

(1.土木工程防灾国家重点实验室(同济大学)，上海200092； 2.中国路桥工程有限责任公司，北京100011； 3.武汉科技大学城市建设学院，武汉430081)

风致稳定性是大跨度桥梁设计中的控制性因素，气动外形优化是改善桥梁气动性能最基本有效的方法[1-2]。桥梁工程领域中气动外形优化的经典方法是基于风洞试验的迭代启发式过程，此方法依赖于风洞试验室的资源消耗和研究人员的精力投入，导致优化设计的范围有限，效果欠佳[3]。近年来CFD技术的快速改进使其在应用工业领域大放异彩，并且早在1997年，将其与数学优化算法结合用于挑选最优气动外形的想法就已经开始引起科学家的关注与研究[4]。目前，在机械、汽车和航空航天领域，基于CFD和数学优化算法的气动外形优化方法已被应用于配置最佳几何外形车辆、飞机机身和机翼、压缩机叶片、层流和湍流扩散器等，这种方法不仅能够严谨彻底地搜索整个设计域，并且消除了经典方法固有的风洞试验资源耗费和迭代试错过程效率低下的缺点[5-7]。由于CFD在土木工程领域中引入较晚[6-7]，导致这种结合数学策略的优化方法研究也有所滞后，因此关于将其移植到桥梁工程领域是否可行的研究显得需求迫切。

本文研究了基于数值计算与数学策略的流线型箱梁气动外形优化方法。数值计算是利用CFD模拟和颤振时域法进行计算，数学策略则是将试验设计、蚁群混合遗传算法与Kriging代理模型结合，通过求解某桥梁断面整个设计域中部分的解，构建和训练Kriging近似代理模型，得到满足一定精度的整个设计域所有的解，从而寻求最优断面。作为案例研究，以颤振性能为目标，优化苏通桥主梁截面外形。

1 优化准备

1.1 问题描述

在工程上，流线型箱梁气动外形优化是指在满足行车、静力、稳定等要求的前提下，通过主梁几何外形的修剪使其具有最佳的气动性能。相应地在数学上，流线型箱梁气动外形优化是将其几何形状表达为梁高、梁宽、风嘴角度上下腹板倾角几个设计变量的集合，称为设计参数；设计变量本身有一定的限制，一般六车道使得桥面宽度满足相应的不等式，流线型要求上下腹板倾角满足相应的关系式等等，这些关系式表达为约束；气动性能分为颤振性能、涡振性能等，在数学上对应目标函数值。鉴于以上，流线型箱梁气动外形优化可表达为多目标多约束多设计变量下的数学优化问题。

1.2 优化流程

由于主梁断面优化问题包含大变化范围的多设计变量、强敏感性的目标函数值，设计完全试验会产生巨量的试验次数；且风洞试验对人力、物力、时间成本消耗较大，因此在保证精度的前提下，需要寻找一种新的优化策略，流程见图1。策略分为3个组成部分：(a) 前期准备选定断面气动外形敏感参数及优化范围，在设计域中进行试验设计；(b) 优化过程以数值模拟为工具计算风致响应并结合数学策略(近似代理模型建立、验证与更新)得到整个设计域内的风致响应；(c) 结果分析是对各个参数进行统计分析，深入探究机理，最终选出气动外形良好断面。

1.3 设计域

将流线型箱梁断面气动敏感变量分为断面端部变量如风嘴角度、上下腹板倾角等；断面中部变量如梁高、梁宽等。目前已建造的大跨度桥梁中，断面端部总是采用较长的下斜腹板，且下腹板倾角与风嘴角度、风嘴长度等相互关联，因此在端部外形变量中，下腹板倾角对气动性能影响占比大；断面中部的梁高和梁宽两个参数一般不宜做出较大改变，因其不仅对气动性能影响较大，且会引起主梁抗压弯静力性能、质量和惯性矩动力性能等发生变化，但考虑到优化效果，认为可在一定范围内变化[8]；综上，选取主梁端部变量下腹板倾角A和中部变量梁高H作为设计变量。以苏通长江大桥主梁流线型箱梁断面为基准断面(A=13°,H=4.0 m)[9]，A和H两参数的优化范围取80%～120%，1%为梯度，设计域为{H:[3.2, 4.8],A:[10.4, 15.6] }。

1.4 试验设计

若在设计域中进行完全试验，则试验总次数为1 681次，故采取合理的试验设计方法制定抽样方案。典型试验设计方法有正交试验设计、均匀试验设计、(空间填充)拉丁超立方设计[10]等，示意见图2。均匀试验设计能避免采样随机性、易于选取，且可以提供设计全空间的信息；序列探索是基于特定模型， Lin[11]提出了一种序列探索试验设计(sequential exploratory experiment design, SEED)方法产生新的样本点来更新 Kriging 模型；综上，采用均匀试验+序列探索试验设计方法制定抽样方案，示意见图2(c)。

图1 流线型箱梁气动外形优化流程Fig.1 Aerodynamic shape optimization process of streamlined box girder

2 数值计算

颤振临界风速采用混合数值计算。根据设计变量A/H变化，先通过CFD模拟竖弯和扭转强迫振动得到三分力系数时程，再用最小二乘法拟合颤振导数并建立气动力模型，最后建立有限元结构模型以二维时域方法求解颤振临界风速。以苏通桥原断面作为算例对计算方法进行说明并检验准确性。

图2 典型试验设计Fig.2 Typical experimental designs

2.1 网格及模型

CFD数值模拟直接影响到近似模型响应面的精度和优化算法的效率, 是整个优化流程的基础。因此，合理的计算域划分与网格数量、湍流模型、雷诺数[12]等至关重要。二维模型计算域及网格划分见图3，计算域的左侧设置为速度入口，右侧设置为压力出口，上下壁面设置为滑移壁面，其中B表示梁宽，在距断面中心10B位置处设置左右交界面以及在距断面中心5B位置处设置环形交界面。流线型箱梁断面中心距入口20B、出口30B、上下边界10B，出入口边界及上下壁面边界已远离断面，数值模拟能够反映真实物理环境。数值风洞入口风速12 m/s，雷诺数与苏通桥试验报告相同，取为4.64×105，湍流模型采用学界广泛应用的SSTk-ω模型[2]。为合理划分网格并减少网格数量，整个计算域内采用分块划分网格的方式：ZONE1、ZONE2和ZONE4区域采用结构化网格，ZONE3区域采用结构化与非结构化混合网格，箱梁断面壁面边界层为结构化网格，网格厚度为2.5×10-5m。将ZONE3设置为竖向和扭转区域，箱梁断面运动时，交界面内网格随之运动，保证了边界层的正确模拟，设定ZONE2为变形网格区域，应用局部网格重划与弹簧光顺法协调控制网格的变形。

图3 计算域及网格划分示意Fig.3 Schematic of computational domain and meshes

2.2 无关性验证

计算域、网格划分方式及基本参数设定之后，在进行CFD 计算之前，还需对不同的计算参数进行必要的试算，以期得到兼顾效率和准确性的最佳计算参数，网格无关性与时间无关性验证便是选取网格数量与时间步长的必经步骤。网格选取遵守双倍网格尺寸变化规则进行，对图3所示网格，ZONE1、3、4分区网格数量不变，只对ZONE2进行网格重划分以得到不同总数的网格，ZONE2网格分辨率随着网格数量增加而增大。网格名称以Mesh为前缀、中间下划线、网格总数为后缀进行命名，不同网格及网格无关性验证结果(0°攻角)见表1，误差在可接受范围内。基于桥梁断面长度与入口边界风速考量，选取计算步长为5×10-4s、2×10-3s和1×10-3s，三者计算偏差在0.1%以内。最终选取计算网格总数为12万，时间步长为1×10-3s。

表1 不同网格数量及数值结果Tab.1 Grids of different quantities and numerical results

2.3 强迫振动与自激力

国内外许多学者在颤振分析的计算方法中进行过深入探索，综合考虑计算效率与准确性，选用二维时域法[13]进行颤振临界风速计算。采用文献[14]的脉冲响应函数表达自激力，并将其施加于模拟节段模型风洞试验的有限元模型以计算颤振临界风速。

图4 最小二乘拟合颤振导数Fig.4 Flutter derivatives by least squares fitting

2.4 颤振临界风速求解

在ANSYS中建立有限元结构分析模型，参考风洞试验节段模型设计计算，如图5所示，弹簧采用刚度相同杆单元BEAM4模拟，其余构件采用刚性杆单元模拟，每延米质量和转动惯量采用质量单元MASS21施加；箱梁断面几何缩尺比选用1∶70，风速比1∶8.367，各项参数见表2；给定风速即可进行计算。

表2 结构分析建模参数Tab.2 Structural analysis modeling parameters

当风速取值为U=21～22 m/s时，扭转位移时程见图6：在来流风速为21.8 m/s时接近于等幅振动但仍收敛，22 m/s时则明显运动发散，故认为颤振临界风速为21.9 m/s，此时换算实桥风速为182.5 m/s，节段模型试验结果为185 m/s，计算误差为1.4%，认为计算非常准确，耗时5～10 min。综上，二维时域法可以作为数值计算工具完成优化过程中每一轮新断面的颤振临界风速求解。

图5 有限元结构分析模型Fig.5 Finite element model of structural analysis

图6 二维时域法计算扭转位移时程Fig.6 Torsional displacement time history calculated by two-dimensional time domain method

3 数学策略

数学策略优化思路为近似代理模型：首先用均匀试验设计结合数值计算构造元模型并验证精度，然后用序列探索试验设计方法搜索模型最优解和误差最大解并与真实数值计算值对比验证，若不满足精度要求则加入新计算的真实值更新模型，直到优化结束，更新模型转变为末模型。

3.1 元模型及验证

统计学习的目的是使学到的近似代理模型不仅对已有数据有很好的拟合，更重要的是对新数据有很好的预测，因此选取合适的代理模型、评价指标、验证更新方法等至关重要[16]。

Kriging模型具有随机特性，可提供预测位置的预测值和预测误差，适用于全局优化算法[17]，因此选取Kriging进行统计学习，得到元模型见图7(a)；进一步通过简洁而有效的留一交叉验证来合理评价模型精度[20]。选取相对误差、相对误差均值、标准差以及预测趋势作为评价指标。预测相对误差、相对误差均值与标准差为：

(1)

(2)

(3)

表3 留一交叉验证相对误差Tab.3 Relative error of leave-one-out cross-validation

3.2 代理模型更新

对于Kriging模型某些点(如采样点9)真实值与预测值偏差较大的问题，解决办法是采用序列探索试验设计方法增加采样点对模型进行更新。Kriging模型的预测值是一个随机变量，模型可以给出预测均值EI和均方误差MSE，基于这种随机预测特性，序列探索实验设计可以对模型的EI和MSE两个指标进行搜索以增加新的样本点，最终使得构建的模型能够捕捉到真正对象函数的趋势和变化[18]。

将初始样本点和所有的更新点用于建立Kriging模型，得到流线型箱梁-Kriging末模型见图7(c)。由EI值响应面可直观得出最大值出现在左上角，而不是元模型的右下角，这也说明了模型更新的必要性。

表4 流线型箱梁-Kriging模型更新过程Tab.4 Updating process of streamlined box girder-Kriging model

图7 流线型箱梁-Kriging建立及更新Fig.7 Modeling and updating process of streamlined box girder-Kriging model

3.3 参数统计分析

利用末模型得到设计域内每个采样点的颤振临界风速V预测值，根据文献[20]方法对优化参数A和H进行直观分析和方差分析，以探究它们对V的影响规律和程度，统计分析结果见图8。图8(a)给出ΔA对V的影响规律，V随A增加呈减小趋势，以分隔线为界趋势明显呈两个相似区间Ⅰ和Ⅱ，在各自区间内V类似呈正弦变化；图8(b)给出了ΔH对V的影响规律，V随H增加呈增大趋势，以分隔线为界趋势呈两个相似区间Ⅲ和Ⅳ，在各自区间内V类似呈余弦变化。根据V变化范围来看，A影响区间为[165, 190]，而H影响区间为[175, 185]，所以A影响相对较大。

图9给出了方差分析结果，扇形比例表示对应项的平均方差占平均方差总和的比例。根据方差分析理论，一个因素的平均方差所占比例越高，表明该因素对指标变化的影响越明显。A、H和A&H交互作用所占比例分别为63%、12%和25%，3个指标对颤振临界风速的影响依次递减，验证了直观分析所得结论。A和H两个参数不仅能各自独立对V产生影响，还会联合起来产生影响，这称为交互作用。如图10所示，以ΔH=(-20%, 0, 20%)作为研究对象，探究ΔA=[-20%, 20%]变化范围内V变化趋势。ΔH=-20%和0时，曲线趋势同增同减，无明显的交互作用；在ΔH=20%，ΔA=[0, 7%]时，曲线呈下降趋势，另外两条则呈上升趋势，且3条曲线出现交叉，此时交互作用明显。原因可能是梁高较大的断面钝体特性明显，下腹板倾角产生的影响大。