橋梁斷面顫振導數的分狀態多頻強迫振動識別

王林凱, 劉志文， 陳政清

(1. 安徽省交通規劃設計研究總院股份有限公司，合肥 230000；2. 湖南大學 風工程與橋梁工程湖南省重點實驗室，長沙 410082)

橋梁斷面顫振導數是研究橋梁氣動性能的重要參數，可通過風洞試驗或計算流體力學數值模擬獲得。風洞試驗識別橋梁顫振導數可以分為自由振動法與強迫振動法兩類。自由振動法試驗測試相對簡單，Scanlan等[1]最早提出分階段的顫振導數識別方法，即先通過豎彎和扭轉單自由度振動試驗識別出直接顫振導數，再用耦合振動試驗識別出耦合項顫振導數。Poulsen等[2]通過風洞試驗用自由振動方法識別了大帶東橋梁斷面的顫振導數。之后，國內外部分學者對自由振動識別方法進行了進一步研究與與應用[3-8]。強迫振動法則需要能夠驅使橋梁斷面以某一頻率和振幅做簡諧運動的專門裝置。早在1952年，Halfman[9]采用強迫振動法直接測量了機翼非定常氣動力，采用不同振動形式識別了不同攻角的顫振導數。Ukeguchi[10]首次將強迫振動測力法用于橋梁斷面顫振導數的識別。Li[11]采用強迫振動法在水洞中識別了橋梁斷面的顫振導數。陳政清等[12]率先在國內實現了用強迫振動法識別橋梁斷面顫振導數，并且研究了顫振導數識別的時域和頻域方法。郭震山[13]、牛華偉[14]在此基礎上進行了三自由度橋梁斷面的顫振導數識別與研究。相比自由振動識別方法，強迫振動法所需裝置復雜，但是具有程序穩健性高，數據重復性好，可控制折算風速范圍廣等優點而在近年來的工程以及研究應用中被廣泛地采用。

隨著計算流體力學的發展，橋梁斷面氣動導數的強迫振動數值模擬識別方法逐漸受到關注。Larsen[15]借助離散渦方法模擬了大帶東橋主梁截面的繞流場，并用分狀態強迫振動法識別了該橋主梁斷面的顫振導數。后來該方法得到了廣泛的研究與應用[16-19]。崔益華等[20]使用了耦合強迫振動數值模擬識別法，識別了橋梁斷面的顫振導數，并計算相應得橋梁顫振臨界風速，計算結果與試驗吻合較好，驗證了該方法的可行性。

綜上所述，目前關于橋梁斷面氣動導數試驗識別采用耦合強迫振動法或自由振動識別法，數值模擬識別則以分狀態單頻強迫振動識別法為主。隨著計算流體動力學的發展，從提高橋梁主梁斷面顫振導數識別精度和效率的角度考慮，進行橋梁斷面顫振導數快速識別方法的研究仍具有十分重要的價值和意義。本文在現有研究成果的基礎上，分別進行了薄平板斷面和流線型主梁斷面顫振導數彎扭耦合兩自由度強迫振動識別和分狀態多頻強迫振動識別數值模擬研究。

1 數值方法

1.1 流體力學控制方程

直角坐標系下，黏性不可壓縮流體基于雷諾平均(RANS)的連續性方程和N-S方程分別為：

(1)

(2)

式中：i,j=1，ρ為空氣密度，ρ=1.225 kg/m3，μ為動力黏性系數，μ=1.789 4×10-5kg/m·s。

1.2 湍流模型

(3)

渦黏模型包括零方程模型，一方程模型以及兩方程模型。常用的兩方程RANS渦黏模型有標準k-ε湍流，標準k-ω模型以及SSTk-ω模型等。

在標準k-ε模型和k-ω模型的基礎上，Menter[23]提出了SST(Shear Stress Transport)k-ω模型，湍流方程表示為：

(4)

(5)

式中：μeff, k與μeff, ω表示湍流項的等效黏度;Sk與Sω表示湍動能和比耗散率的源項。

該模型將k-ε模型中的ε方程寫成ω的形式，并且將其與標準k-ω模型進行線性組合，組合系數則表示為關于與壁面之間距離的函數。SSTk-ω模型湍流模型通過線性組合函數隨與壁面之間距離的變化現實了從近壁面的k-ω模型向遠離壁面的k-ε模型轉變，從而結合了兩種模型各自的特點，在黏性子層及遠離壁面湍流充分發展區都具有很好的計算性能。

1.3 動網格技術

對于通量φ，在任一控制體V內，其邊界運動時，N-S方程的積分形式表達式可以表示為：

(6)

式中：us為網格運動的速度。

為實現結構運動，采用“剛性網格區域+動網格區域”的方案，以薄平板為例的計算區域分區示意圖，如圖2所示。

圖1 薄平板網格區域劃分示意圖(mm)Fig.1 Schematic diagram of plate grid division(mm)

薄平板附近包圍一圈剛性運動區域網格，該部分網格隨薄平板一致做剛性運動。在剛性運動區域網格外設置一層三角形動網格區域，該區域網格會在運動過程中重新劃分網格。最外層為靜止區域網格，該區域內網格在計算過程中不會發生任何變化。計算域左側邊界為速度入口，計算域右側邊界為壓力出口，參考壓力設為零；計算域上、下側邊界為對稱邊界條件，模型表面為無滑移邊界。

動網格更新方法選擇“彈性光順+局部網格重構”的方式，既保證網格更新效率又可保證網格更新質量。動網格驅動宏選擇DEFINE_CG_MOTION，以指定剛性網格運動區域的相應運動速度，進而實現網格位移的更新。

2 主梁斷面氣動導數識別方法

2.1 氣動自激力描述

將Scanlan線性氣動自激力表達式表述為基于耦合強迫振動的形式：

(7)

(8)

(9)

式中：fh與fα分別表示強迫振動的豎向振動頻率和扭轉振動頻率。

2.2 彎扭耦合強迫振動氣動導數識別方法

文獻[20]中的彎扭耦合強迫振動識別方法為進行兩次彎扭同頻率的強迫振動識別，第一次彎扭位移同相位，第二次強迫振動彎扭位移設置為180°的相位差。為提高識別效率，本文選擇文獻[14]中的將彎扭頻率設為不一致的識別方法。假設結構做耦合運動形式位移表示為：

(10)

通過CFD計算獲取的氣動自激力L與M的時程曲線，按照式(4)與式(5)的形式進行最小二乘擬合獲取相應的顫振導數值。通過不同的彎扭頻率進行錯位組合識別出感興趣折算風速下的顫振導數。

以CFD耦合強迫振動獲取的升力系數和升力矩系數為數據來源，其三分力系數定義如下：

(11)

采用Matlab語言時域識別步驟可以表示為：

步驟1 賦予初值a=[0.1 0.1 0.1 0.1]；

步驟3 定義相關系數矩陣量為：

步驟4 定義目標擬合函數：

Fun=nlinfit(a,X)(a1N1X1+a2N2X2+a3N3X3+a4N4X4)

步驟5 進行非線性擬合：

H=nlinfit(X,CL,Fun,a)

A=nlinfit(X,CM,Fun,a)

2.3 分狀態多頻強迫振動氣動導數識別方法

在Scanlan線性自激力可疊加性的理論框架下，推出一種新的識別方式。假設模型分別進行單自由度多頻率強迫振動，即

(12)

每一振動項對顫振導數的影響系數取各自的振動頻率。當結構做豎向分狀態多頻強迫振動時，對應的升力和升力矩系數分別為：

(13)

(14)

當結構做扭轉分狀態多頻強迫振動時，對應的升力和升力矩系數分別為：

(15)

(16)

按式(12)分別驅動斷面做豎向單自由度多頻振動與扭轉單自由度多頻振動，通過兩次計算即可獲得主梁斷面的氣動力時程曲線，進而可識別出不同折算風速對應的橋梁斷面氣動導數。

3 典型斷面氣動導數識別

3.1 斷面參數

分別選取薄平板[21]和大帶東橋主梁斷面進行斷面氣動導數識別研究。薄平板斷面高D=20 mm，寬度B=450 mm，其外形如圖2(a)所示。大帶東橋主梁斷面高D=4.4 m，寬B=31 m，高寬比為7.05，斷面上有3%的橫坡。該斷面的剪切中心(S.C)位于距離橋面0.465倍截面高度處，其外形如圖2(b)所示。

(a) 薄平板尺寸(mm)

(b) 大帶東橋梁斷面尺寸(m)圖2 模型斷面參數Fig.2 The Section parameters of models

3.2 網格劃分與求解設置

利用Gambit繪制相應的網格如圖3～4所示。薄平板邊界處使用邊界層4∶2網格能夠有效的減少網格數量，并且與薄平板兩端以結構化網格形式拓撲出去使得網格邊長增大導致相鄰網格長度比增大的現象相平衡，三角形網格數為23 840個，四邊形網格數為13 208個。大帶東橋主梁斷面采用外包橢圓加圓的方式向外拓撲，四邊形網格數量為39 941個，三角形網格數量為65 774，壁面處首層網格高度為1×10-5B。大帶東橋梁斷面網格文件導入Fluent中時需要進行80∶1縮尺，斷面寬度即為0.387 5 m。

圖3 薄平板網格示意圖Fig.3 Schematic diagram of flat plate grid

圖4 大帶東橋梁斷面網格圖Fig.4 Grid of the main girder section of the Great Belt East Bridge

采用商業流體力學計算軟件Ansys Fluent進行求解計算。計算選擇速度-壓力解耦的SIMPLEC算法，湍流模型選擇SSTk-ω模型，壓力方程采用二階格式離散，動量方程、湍動能方程及比耗散率方程均采用二階迎風格式離散，計算迭代殘差為10-6，給予初始均勻流場湍流特征，湍流強度為0.5%，湍流黏性比為2。薄平板流場相對簡單，湍流模型對計算結果的影響較小[24]。大帶東橋主梁斷面流場較為復雜，在強迫網格做簡諧運動前，先對大帶東橋梁斷面做靜止繞流計算以驗證網格與求解設置的可靠性。大帶東橋主梁斷面靜止繞流計算中，來流速度取12 m/s，計算時間步長取0.000 1 s，計算至流場充分穩定。St數定義如下：

(17)

式中：fs為渦脫頻率；D為斷面迎風面高度。計算結果如表1和圖5所示。由表1可知，本文的大帶東橋主梁斷面網格以及求解設置滿足精度要求，壁面處Y+分布符合湍流模型的要求。

圖5 大帶東橋主梁斷面靜止繞流計算結果Fig.5 The calculation results of flow of Great Belt East Bridge

3.3 彎扭耦合兩自由度強迫振動氣動導數識別

不同于文獻[14]中試驗過程中通過改變風速來改變折算風速，本文通過改變運動頻率來改變折算風速，計算工況如表2和表3所示，各個工況中的頻率的設置盡量使得折算風速為一整數值，為此薄平板算例中的來流速度值取10倍的薄平板寬度，即4.5 m/s，大帶東橋梁斷面的來流速度值取24倍的梁寬，即9.3 m/s。計算時，模型豎向振動的振幅為0.02B，扭轉振幅為2°，根據采樣定理，計算時間步長建議不大于模型驅動周期的0.02倍，本次統一取0.000 5 s，計算時間建議不小于各運動成分的10個周期時長。

薄平板顫振導數的識別結果如圖6所示，大帶東橋主梁斷面顫振導數識別結果如圖7所示。

表2 薄平板顫振導數計算工況Tab.2 Calculation conditions of flutter derivatives of the thin plate

表3 大帶東橋梁斷面顫振導數計算工況Tab.3 Calculation conditions of the flutter derivatives of the section of the Great Belt East Bridge

圖6 薄平板顫振導數識別結果Fig.6 Flutter derivative identification results

從圖6中可以看出，薄平板各個顫振導數均與平板理論解吻合較好，驗證了本文的耦合強迫振動識別方法的正確性與可靠性。從圖7中可以看出，大帶東橋主梁斷面顫振導數結果與相關文獻數值模擬的結果較為接近，并且與Poulsen等的試驗結果吻合較好。需要說明的是，祝志文等采用的是分狀態單自由度強迫振動法識別的顫振導數。

圖7 大帶東橋梁斷面顫振導數計算結果Fig.7 Calculation results of the flutter derivatives of girder section of the Great Belt East Bridge

3.4 氣動自激力可疊加性CFD驗證

(a) 升力時程曲線

(b) 力矩時程曲線圖8 氣動自激力對比圖Fig.8 Comparison of self-excited aerodynamic forces

從圖8中可以看出，由CFD獲取的升力、力矩時程曲線與用耦合振動法識別的顫振導數擬合表達式計算的時程曲線幾乎完全吻合，從側面為耦合強迫振動法提供了理論支持，驗證了該橋梁斷面的氣動自激力的可疊加性。同時也證明了CFD氣動力計算的準確性并不因為模型結構的運動復雜性提高而改，這與文獻[24]的結論是一致的。

3.5 分狀態多頻強迫振動氣動導數識別

在線性氣動自激力可疊加性的基礎上，按表4和表5中工況進行分狀態多頻強迫振動識別薄平板和大帶東橋主梁斷面的顫振導數。

表4 薄平板分狀態多頻強迫振動顫振導數識別工況Tab.4 Calculation conditions of the identification of flutter derivatives of the thin plate by step-by-step forced vibration with multi frequency

表5 大帶東橋主梁斷面分狀態多頻強迫振動顫振導數識別工況Tab.5 Calculation conditions of the identification of flutter derivatives of the girder section of Great Belt East Bridge by step-by-step forced vibration with multi frequency

驅動斷面運動的振幅仍然按耦合識別法設置，即豎彎0.02B，扭轉位2°，計算獲取的氣動力系數時程曲線以及氣動力時程曲線FFT幅值譜如圖9～16所示。

(a) 升力時程曲線

(b) 力矩時程曲線圖9 薄平板豎向多頻強迫振動自激力時程曲線Fig.9 The history curves of self-excited forces of thin flat plate by forced vertical vibration with multi frequency

(a)升力時程曲線FFT幅值譜

(b)力矩時程曲線FFT幅值譜圖10 薄平板豎向多頻強迫振動自激力FFT幅值譜Fig.10 The FFT amplitude spectrum of forced vertical vibration self-excited forces of thin flat plate with and multi frequency

(a) 升力時程曲線

(b) 力矩時程曲線圖11 薄平板扭轉多頻強迫振動自激力時程曲線Fig.11 The history curves of self-excited forces of thin flat plate by forced torsional vibration with multi frequency

(a) 升力時程曲線FFT幅值譜

(b) 力矩時程曲線FFT幅值譜圖12 薄平板扭轉多頻強迫振動自激力FFT幅值譜Fig.12 The FFT amplitude spectrum of forced torsional vibration self-excited forces of thin flat plate with multi frequency

(a) 升力時程曲線

(b) 力矩時程曲線圖13 大帶東橋主梁斷面豎向多頻強迫振動自激力時程曲線Fig.13 The history curves of self-excited forces of the Great Belt East Bridge girder section by forced vertical vibration with multi frequency

(a) 升力時程曲線FFT幅值譜

(b) 力矩時程曲線FFT幅值譜圖14 大帶東橋主梁斷面豎向多頻強迫振動自激力FFT幅值譜Fig.14 The FFT amplitude spectrum of forced vertical vibration self-excited forces of the Great Belt East Bridge girder section with multi frequency

(a) 升力時程曲線

(b) 力矩時程曲線圖15 大帶東橋主梁斷面扭轉分狀態多頻強迫振動自激力時程曲線Fig.15 The history curves of self-excited forces of Great Belt East Bridge girder section by forced torsional vibration with multi frequency

(a) 升力時程曲線FFT幅值譜

(b) 力矩時程曲線FFT幅值譜圖16 大帶東橋主梁斷面扭轉分狀態多頻強迫振動自激力FFT幅值譜Fig.16 The FFT amplitude spectrum of forced torsional vibration self-excited forces of the Great Belt East Bridge girder section with multi frequency

從圖10、12、14、16中均可以看出，升力與力矩系數中，各個預設的振動頻率能量均較為集中，并且包含了所有預設的扭轉振動頻率。

利用上述數據識別的薄平板顫振導數結果與Theodorson理論解進行對比，結果如圖17所示；大帶動橋主梁斷面顫振導數結果與耦合振動法識別的結果進行對比，如圖18所示。

圖17 薄平板彎扭分狀態識別法顫振導數識別結果Fig.17 Identification of flutter derivatives of a thin plate with a multi-frequency step-by-step identification method

圖18 大帶東橋斷面彎扭分狀態多頻識別法顫振導數識別結果Fig.18 Identification of flutter derivatives of the Great Belt East Bridge main girder section with a multi-frequency step-by-step identification method

表6 大帶東橋梁斷面顫振導數不同計算方法識別效率對比Tab.6 Comparison of the efficiency of different methods for identifying flutter derivatives of the Great Belt East Bridge

4 結 論

采用計算流體動力學方法對橋梁斷面氣動導數識別方法進行了研究，取得了如下主要研究成果：

(1)分別針對薄平板和大帶東橋主梁斷面，采用彎扭耦合強迫振動方法進行了主梁斷面氣動導數識別，識別結果分別于西奧多森理論解和已有文獻結果吻合較好，驗證了本文彎扭耦合強迫振動識別方法的精度。

(2)在驗證了氣動自激力可疊加性的基礎上，提出橋梁斷面氣動導數識別的分狀態多頻強迫振動識別法，計算結果表明該方法的精度與傳統彎扭耦合識別方法精度相當，而計算效率得到較大的提升。