橋梁斷面顫振導數(shù)的分狀態(tài)多頻強迫振動識別

王林凱, 劉志文， 陳政清

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

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

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

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

1 數(shù)值方法

1.1 流體力學控制方程

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

(1)

(2)

式中：i,j=1，ρ為空氣密度，ρ=1.225 kg/m3，μ為動力黏性系數(shù)，μ=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-ω模型進行線性組合，組合系數(shù)則表示為關于與壁面之間距離的函數(shù)。SSTk-ω模型湍流模型通過線性組合函數(shù)隨與壁面之間距離的變化現(xiàn)實了從近壁面的k-ω模型向遠離壁面的k-ε模型轉變，從而結合了兩種模型各自的特點，在黏性子層及遠離壁面湍流充分發(fā)展區(qū)都具有很好的計算性能。

1.3 動網(wǎng)格技術

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

(6)

式中：us為網(wǎng)格運動的速度。

為實現(xiàn)結構運動，采用“剛性網(wǎng)格區(qū)域+動網(wǎng)格區(qū)域”的方案，以薄平板為例的計算區(qū)域分區(qū)示意圖，如圖2所示。

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

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

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

2 主梁斷面氣動導數(shù)識別方法

2.1 氣動自激力描述

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

(7)

(8)

(9)

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

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

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

(10)

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

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

(11)

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

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

步驟3 定義相關系數(shù)矩陣量為：

步驟4 定義目標擬合函數(shù)：

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

步驟5 進行非線性擬合：

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

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

2.3 分狀態(tài)多頻強迫振動氣動導數(shù)識別方法

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

(12)

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

(13)

(14)

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

(15)

(16)

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

3 典型斷面氣動導數(shù)識別

3.1 斷面參數(shù)

分別選取薄平板[21]和大帶東橋主梁斷面進行斷面氣動導數(shù)識別研究。薄平板斷面高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 模型斷面參數(shù)Fig.2 The Section parameters of models

3.2 網(wǎng)格劃分與求解設置

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

圖3 薄平板網(wǎng)格示意圖Fig.3 Schematic diagram of flat plate grid

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

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

(17)

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

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

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

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

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

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

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

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

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

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

3.5 分狀態(tài)多頻強迫振動氣動導數(shù)識別

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

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

表5 大帶東橋主梁斷面分狀態(tài)多頻強迫振動顫振導數(shù)識別工況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°，計算獲取的氣動力系數(shù)時程曲線以及氣動力時程曲線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 大帶東橋主梁斷面扭轉分狀態(tài)多頻強迫振動自激力時程曲線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 大帶東橋主梁斷面扭轉分狀態(tài)多頻強迫振動自激力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中均可以看出，升力與力矩系數(shù)中，各個預設的振動頻率能量均較為集中，并且包含了所有預設的扭轉振動頻率。

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

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

圖18 大帶東橋斷面彎扭分狀態(tài)多頻識別法顫振導數(shù)識別結果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 大帶東橋梁斷面顫振導數(shù)不同計算方法識別效率對比Tab.6 Comparison of the efficiency of different methods for identifying flutter derivatives of the Great Belt East Bridge

4 結 論

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

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

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