流線型箱梁氣動外形的數值優化

孫 強，周志勇，胡傳新，秦 鵬，黃德睦

(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條曲線出現交叉，此時交互作用明顯。原因可能是梁高較大的斷面鈍體特性明顯，下腹板傾角產生的影響大。

圖8 優化參數統計分析結果Fig.8 Statistical analysis results of optimized parameters

圖9 平均方差分布Fig.9 Average variance distribution

圖10 典型交互作用示意Fig.10 Schematic of typical interaction

3.4 最優斷面

對上述統計分析總結并進行最值搜索，得到氣動外形最優斷面見圖11。虛線為原始蘇通橋斷面，實線為最優斷面，梁高增加20%，即4.8 m，下腹板傾角減小17%，即10.79°，顫振臨界風速為198 m/s，相比于原型斷面增加了15 m/s。

圖11 流線型箱梁最優化斷面與原斷面對比Fig.11 Comparison of optimized section and original section of streamlined box girder

4 結 論

以流線型箱梁氣動外形優化數值策略為研究對象，將蟻群混合遺傳算法、Kriging近似代理模型、CFD數值風洞模擬和ANSYS顫振時域法結合使用，得出一套完整的氣動外形優化方法，并以蘇通橋斷面優化為例檢驗了優化方法可行性。主要結論如下：

1) 采用CFD進行強迫振動識別顫振導數和ANSYS顫振時域法計算顫振臨界風速為182.5 m/s，與風洞試驗185 m/s相差較小，數值計算手段準確性較高。

2) 結合蟻群混合遺傳算法，Kriging數學代理模型更新13次后達到2%以內的預測誤差，數學代理模型預測顫振響應可行性較高。

3) 借助統計分析數學手段進行優化分析表明，下腹板傾角A相對于梁高H對流線型箱梁斷面顫振性能影響更大，兩個參數存在交互作用。

4) 優化結果顯示，基于蘇通長江大橋的流線型箱梁最優參數匹配斷面為梁高H=4.8 m，下腹板傾角A=10.79°。

5) 基于數值計算與數學策略的流線型箱梁氣動外形優化策略能夠應用于流線型箱梁斷面氣動外形優化，對今后橋梁斷面選型具有借鑒意義。