李威霖，牛華偉，華旭剛，陳政清

（湖南大學風工程與橋梁工程湖南省重點實驗室，湖南長沙 410082）

引 言

隨著橋梁跨徑增長至千米級和高聳建筑高度增至300 m 以上，由大氣湍流產生的抖振荷載可能成為結構的控制性動力荷載。復氣動導納函數(shù)是精確計算抖振時域荷載的關鍵參數(shù)［1?2］，是脈動來流與抖振力在幅值與相位上的傳遞函數(shù)。對于流線型斷面，氣動導納可用Sears 函數(shù)進行表示，但對于鈍體斷面，由于其復雜的流動分離現(xiàn)象，其氣動導納理論表達式仍是風工程領域尚未解決的難題。

氣動導納現(xiàn)有研究主要依賴于風洞試驗，根據(jù)來流特性不同可分為被動格柵產生的寬頻脈動湍流下識別［3?4］和主動風洞產生的諧波來流下識別［5］。被動格柵雖能模擬大氣邊界層風場的部分湍流特性，但較難產生大尺度的低頻漩渦［6］，而模擬的小尺度高頻漩渦又對鈍體斷面的流動分離影響很大。因此采用被動格柵風場識別的氣動導納結果常有一定的離散性［7］，識別結果需要進行重復驗證。

諧波脈動來流通常由多個振動機翼組成的主動格柵裝置產生，也可在多風扇主動風洞內產生［8］。在該風場下識別的氣動導納重復性較好，便于工程應用，如薄機翼斷面的氣動導納Sears 函數(shù)在正弦脈動流場下已被多名學者證實［9?11］。而對鈍體斷面，許多學者如Jancauskas 等［11］、Diana 等［12］、Ma 等［13］在主動風洞中識別了氣動導納，但仍有三個問題尚未得到有效解決：1）諧波風場參數(shù)如風攻角、諧波幅值對復氣動導納的影響；2）單頻率和多頻率組合諧波風場下識別的復氣動導納是否一致；3）結構運動狀態(tài)對復氣動導納的影響。

上述第一個問題涉及到鈍體斷面的氣動導納是否與脈動風幅值有關，這需在不同幅值下的脈動風進行氣動導納識別；第二個問題本質為不同頻率抖振力的可疊加性，Diana 等［12］通過主動風洞產生了多次諧波，驗證了流線型箱梁的豎風向復氣動導納滿足線性疊加性，但對鈍體斷面，升力與阻力氣動導納的可疊加性尚不清楚；第三個問題涉及到由脈動風引起的抖振力與結構運動引起自激力之間的可疊加性，Yan 等［14］發(fā)現(xiàn)在被動格柵湍流風場下結構運動狀態(tài)嚴重影響氣動導納，但在諧波風場下氣動導納是否受結構運動狀態(tài)影響尚不清楚。

由于需要復雜的機械裝置且難以準確控制所產生的諧波參數(shù)［15］，上述問題通過風洞試驗較難解決，相比而言，CFD 數(shù)值模擬在氣動導納識別方面有著巨大的潛力，具有改變來流風場特性的簡易性及流場可視性好等優(yōu)點。如Rasmussen 等［16］、Hejlesen等［17］通過離散渦方法識別了湍流場下橋梁斷面的氣動導納；唐煜等［9］、張偉峰等［10］、及Kavrakov 等［18］采用單頻率諧波識別了一系列橋梁斷面的氣動導納，但對涉及前述三個問題的研究鮮有報道。

鑒于此，本文基于2 維不可壓URANS 方法，首先介紹了在CFD 中如何產生單頻率、多頻率組合的豎風向和順風向脈動氣流，并針對寬高比為4 的矩形斷面，分別識別了其在不同風攻角、諧波幅值下的升力和阻力復氣動導納；然后識別了在多頻率組合諧波風場下的氣動導納來驗證抖振力的可疊加性；最后通過斷面在豎向和扭轉簡諧運動下的復氣動導納研究抖振力與自激力的可疊加性。

1 氣動導納及諧波風場產生

1.1 氣動導納定義

根據(jù)Sears 氣動導納經典理論，當正弦脈動氣流作用于橋梁斷面時，需采用氣動導納來考慮抖振力的非定常特性，單位長度下的抖振升力和阻力系數(shù)可表示為：

式中t為時間，f為脈動風頻率，CL和CD分別為平均升力與阻力系數(shù)，CL，α和CD，α為升力和阻力系數(shù)的導數(shù)，和為順風向和豎風向正弦氣流的復振幅，χLu和χLv分別為順風向和豎風向的升力復氣動導納，χDu和χDv分別為順風向和豎風向的阻力復氣動導納。

為與經典理論對比，將豎向氣動導納與Sears 函數(shù)對比，計算相位時風速點選矩形截面的中點位置處［19］；對順風向氣動導納的模常與Davenport［20］公式對比，而其相位由于仍未有相關的理論表達式，仍采用Sears 函數(shù)的相位，風速點選矩形的迎風側處，具體公式如下：

式中ω=πfB/U；J0，J1分別為第一類Bessel 函數(shù)，K0，K1為修正的第二類Bessel 函數(shù)；對阻力氣動導納，折減頻率k=fD/U。

1.2 諧波脈動流場產生

由以上定義式可知，氣動導納的識別需先產生正弦豎向與順風向脈動氣流。對理想二維不可壓縮流體，豎向和順風向脈動風場的理論表達式為：

式中u，v分別為脈動風速的順風向和豎風向分量；x為平均風速方向；U為平均風速；Au，Av為順風向和豎風向的脈動風幅值。

上述脈動氣流若想在計算域或風洞中實現(xiàn)，必須滿足兩個條件，第一為連續(xù)性歐拉方程，即：

將式（7）和（8）代入式（9）可知，對豎向脈動氣流的表達式已滿足歐拉方程，但對順風向脈動氣流卻不能滿足，從而在入口處會引起不可控的壓力，進而影響數(shù)值計算結果。

第二個條件需考慮上下邊界的影響，即在邊界處脈動風速分量v必為0，可知式（7）不滿足邊界條件，這樣會導致產生脈動氣流幅值會隨速度方向傳輸而折減，進而影響氣動導納的識別精度。

因此為保證所產生脈動氣流滿足上述兩個條件，引入渦旋數(shù)n和余弦函數(shù)以修正上述脈動風公式：

式中H為計算域高度，y為高度方向，將上述公式代入式（9）可得：

由式（10）和（11）可知：當渦旋數(shù)n＝1 時，即可在計算域中部產生僅有豎風向諧波脈動風；當n＝2時可產生順風向諧波脈動風。由于CFD 改變參數(shù)的靈活性，計算中可以控制所產生脈動風的幅值、風攻角，也可產生多頻率組合諧波以驗證氣動力的疊加性。

2 數(shù)值模型及驗證

2.1 4∶1 矩形網格及數(shù)值模型設置

基于OPENFOAM 開源軟件，將上述式（10）和（11）編程，并作用于計算域內入口以產生可控諧波風場。圖1為寬高比B/D=4∶1 的矩形斷面所采用的計算網格，整體尺寸為30B×20B，入口處采用速度入口，雷諾數(shù)為2×105，上部和下部采用對稱邊界。計算網格由4 個不同細化等級的矩形區(qū)域組成，靠近矩形斷面處的網格尺寸最小，細化等級最高，具體如圖2所示。

圖1 4∶1 矩形計算域示意圖（非真實比例）Fig.1 Overview of 4∶1 rectangular section in CFD（the figure does not respect real proportions）

圖2 近壁面網格尺寸圖Fig.2 Detailed meshes around solid boundary

離散方程組求解采用瞬態(tài)PISO 算法，湍流模型為k?ωSST，該設置廣泛應用于二維橋梁斷面三分力系數(shù)、顫振導數(shù)等氣動參數(shù)的識別。圖1的灰色區(qū)域可為動網格，使矩形做豎向或扭轉正弦運動，以研究運動狀態(tài)對氣動導納的影響。由于復氣動導納的識別為單一頻率法，需涉及多工況，因此所有計算均在湖南大學國家超級計算中心進行。

2.2 網格獨立性驗證及三分力系數(shù)

為驗證上述數(shù)值模擬設置的正確性，首先進行了網格獨立性驗證，網格參數(shù)如表1所示，共采用了粗、中、細三種網格，y1為邊界層網格第一層高度，隨著網格細分逐漸降低，而邊界層數(shù)逐漸增大。Y+的平均值從7.69 到1.58 變化，網格總數(shù)量為5.5 萬到14 萬。

表1 不同網格的主要參數(shù)Tab.1 Main parameters for different meshes

圖3為4∶1 矩形斷面在均勻流下不同網格的壓力系數(shù)統(tǒng)計值分布，并與文獻［21］的風洞試驗對比，可知具有較好的網格獨立性，而且與試驗值吻合較好，證明了此數(shù)值模擬設置的正確性。為兼顧計算效率和準確性，以下運算均采用中等網格。

圖3 均勻流下不同網格壓力系數(shù)的統(tǒng)計值分布Fig.3 Distribution of statistical pressure coefficients for different meshes in uniform flow

圖4為±6°攻角范圍內4∶1 矩形斷面的阻力與升力系數(shù)，并分別與文獻［22?23］的CFD 和風洞試驗結果對比可知，雖然風攻角增大存在一定的偏差，但在可接受范圍內。

圖4 平均阻力系數(shù)與升力系數(shù)隨風攻角變化Fig.4 Drag and lift coefficients change as wind attack angles

2.3 諧波脈動流場

為檢查諧波來流效果，首先在空風場下圖1中的坐標原點處測試了風速時程曲線。圖5左側分別給出了順風向、豎風向單頻率及豎風向5 階頻率組合諧波的時程、FFT 曲線?？芍狢FD 流場內產生的諧波與公式（10）和（11）的輸入值在幅值和相位均符合較好；同時右側也給出了相應最后時間步下的流場圖，可知對豎向諧波流場渦旋數(shù)n=1，而順風向流場渦旋數(shù)n=2，驗證了上述理論公式（10）和（11）的正確性，也表明本文提出的方法能在CFD 中產生符合N?S 方程及邊界條件的可控脈動氣流。

圖5 坐標(x，y)=(0，0)處不同諧波風場的風速時程、FFT 曲線及對應的流場圖Fig.5 Time history and FFT of harmonic fields at the measured point(x，y)=(0，0)and the corresponding wind fields

3 氣動導納識別結果

基于前述模擬風場，首先研究風攻角、諧波幅值對豎風向與順風向復氣動導納的影響，然后通過作用多頻率諧波識別的氣動導納來驗證抖振力的線性疊加性，最后讓矩形斷面做豎向及扭轉簡諧運動來驗證抖振力與自激力的可疊加性。復氣動導納的識別由式（3）和（4）得到，其中抖振力的幅值和相位采用正弦函數(shù)擬合其時程曲線得到，較容易消除渦激力及高階分量的干擾。為保證識別精度、減少初始流場影響，每個單頻率下的模擬都保證至少有12 個周期的正弦波，而僅有最后8 個周期的正弦波用于計算氣動導納。

3.1 風攻角

由于4∶1 矩形為對稱截面，以下僅給出氣動導納在風攻角為0°，＋1°，＋3°下的識別結果，其中豎向單頻率諧波幅值Av/U均為0.0175，從而脈動風引起的最大風攻角為1°，對順風向氣動導納識別時，諧波幅值Au/U為0.05。

圖6為升力復氣動導納的絕對值與相位在不同風攻角下的識別結果，識別的折減頻率區(qū)間fB/U為0～0.4，可發(fā)現(xiàn)在0°和＋1°風攻角時絕對值與文獻［11］的實驗值吻合較好，表現(xiàn)為在低頻處大于Sears函數(shù)，并且數(shù)值保持為單位1，符合準定常準則，而后隨折減頻率增大而減少；當風攻角增大到＋3°，氣動導納與Sears 接近；升力氣動導納的相位在不同攻角下均與Sears 函數(shù)顯著不同，表現(xiàn)為隨折減頻率增大而逐漸降低，當風攻角增大到+3°時，在折減頻率為0.1 時發(fā)生了較大變化。

圖6 升力復氣動導納隨風攻角的變化Fig.6 The lift AAF change as wind attack angles

圖7為不同風攻角下的阻力復氣動導納的絕對值和相位變化，折減頻率區(qū)間fD/U為0～0.1，可以發(fā)現(xiàn)阻力氣動導納的模在低頻處接近1，但在高頻處顯著低于Davenport 值，這與Kavrakov 識別箱梁斷面的結論很類似［18］。當攻角增大至+3°，與升力氣動導納類似，氣動導納值開始出現(xiàn)了偏差。阻力氣動導納的相位隨折減頻率先升后降，隨攻角增大也出現(xiàn)一定偏差。

圖7 阻力復氣動導納隨風攻角的變化Fig.7 The drag AAF change as wind attack angles

為探究較大風攻角引起升力與阻力氣動導納變化的原因，可參見如圖8所示的平均流線圖。在1°風攻角時，矩形頂部的負壓區(qū)氣泡先分離然后會發(fā)生再附著，但當風攻角增大至4°時，氣泡分離后不再附著，矩形后緣處的尾渦也發(fā)生顯著改變，從而引起了氣動特性的改變，這與Patruno 對5∶1 矩形截面的研究結論一致［24］。

3.2 脈動風幅值

為研究單頻率下諧波幅值對矩形氣動導納的影響，針對0°平均風攻角，識別了6 種不同幅值下的升力及阻力復氣動導納。如表2所示，其中前3 個工況為僅有豎向脈動風作用，引起的最大瞬時風攻角為1°～4°，后3 個工況為僅有順風向脈動風，幅值為（0.05～0.15）U。

表2 單頻率下不同諧波幅值工況Tab.2 Amplitude of the input harmonic

圖9為升力復氣動導納隨不同豎向諧波幅值的變化。當Av/U<0.035，即豎向脈動風引起的最大風攻角在2°以內時，升力氣動導納的模和相位基本保持不變；當Av/U增大至0.07，最大風攻角為4°時，氣動導納的模有一定偏離。對氣動導納的相位，不同的脈動風幅值對其影響較小。這表明，對矩形斷面的升力氣動導納應考慮脈動風幅值的影響，引起這種現(xiàn)象的原因也可由圖8解釋，當風攻角在±4°范圍內變化時，分離氣泡在附著與不再附著之間變化，導致氣動導納隨脈動幅值變化。

圖8 不同風攻角下平均風速流線圖Fig.8 Streamlines in different wind attack angles

圖9 升力復氣動導納隨豎向諧波幅值變化Fig.9 The lift AAF change as vertical harmonic amplitudes

圖10 為阻力復氣動導納隨順風向諧波幅值的變化曲線，順風向諧波幅值對阻力氣動導納的模和相位在低頻處影響較小，但在折減頻率大于0.05后，增大幅值會引起阻力氣動導納的模增大，對其相位也有一定影響。

圖10 阻力復氣動導納隨順風向諧波幅值變化Fig.10 The drag AAF change as along?wind harmonic am?plitudes

3.3 多頻率組合諧波

通過產生多頻率組合的豎向和順風向諧波風場，驗證4∶1 矩形斷面的抖振力疊加性。如表3所示，前3 個工況為僅有豎風向的包含2，3，5 階頻率的諧波，對順風向驗證了3 階頻率諧波，為避免諧波幅值的影響，設置多頻率組合諧波下的幅值與單諧波幅值保持一致。

表3 多頻率組合諧波工況Tab.3 Multi?frequency harmonic fields

圖11 為不同組合諧波下的抖振升力的傅里葉變換曲線。為方便對比，單諧波工況僅給出折減頻率為0.05 時的結果，可發(fā)現(xiàn)在該頻率處不同諧波工況的抖振升力幅值基本一致，2 階諧波與5 階諧波在折減頻率0.3 處幅值也基本一致。

圖11 不同組合諧波下抖振升力的FFT 曲線Fig.11 The FFT of buffeting lift coefficient under different vertical harmonics

圖12 給出了上述多頻率組合豎向諧波下升力氣動導納結果?？芍?，3 及5 階諧波作用下，該矩形斷面的升力氣動導納在幅值和相位上均具有很好的疊加性，驗證了在諧波脈動風場下，不同頻率的抖振升力具有很好的線性疊加性，這與Diana 等［12］對Humber 橋升力氣動導納的結論一致。

圖12 豎向組合諧波下的升力復氣動導納Fig.12 The lift AAF in different vertical multi?frequency harmonics

圖13 為順風向3 階諧波組合作用下的阻力復氣動導納結果?？芍? 階諧波與單諧波下結果基本一致，表明該斷面的抖振阻力在幅值和相位上也具有很好的疊加性。

圖13 順風向組合諧波下的阻力復氣動導納Fig.13 ThedragAAFindifferentalong?wind multi?frequency harmonics

3.4 運動狀態(tài)

自激力由結構的運動產生，而由脈動風引起的抖振力與自激力之間的可疊加性尚未得到驗證，這也是廣泛應用于橋梁顫抖振分析的理論基礎。大跨橋梁抗風設計關心的折減頻率常在0～1 之間，因此選取豎向和扭轉簡諧運動的折減頻率fB/U=0.25，由于斷面運動和輸入諧波的頻率不同，抖振力和自激力能很精確分離。為避免結構運動引起的瞬時風攻角過大而引起氣動力的非線性，控制豎向和扭轉的運動幅度，以保證由斷面運動引起的最大風攻角在±1°以內變化［25］。

圖14 為升力復氣動導納的模和相位隨不同運動狀態(tài)的變化曲線，可發(fā)現(xiàn)豎向和扭轉簡諧運動對升力復氣動導納影響不大，除在折減頻率為0.2 時有一定偏差外，其他基本一致。這表明在此運動狀態(tài)下，將該斷面的抖振升力與自激力分開處理在諧波風場下是可行的。

圖14 運動狀態(tài)對升力復氣動導納的影響Fig.14 The lift AAF change as different motions

圖15 為阻力復氣動導納隨不同運動狀態(tài)的變化，可知豎向簡諧運動對阻力氣動導納影響較小；但當斷面做扭轉簡諧運動時，阻力氣動導納在高頻處增大，且接近Davenport 理論值。這表明在諧波風場下，由高頻脈動風引起的抖振阻力受扭轉運動的影響較大。

圖15 運動狀態(tài)對阻力復氣動導納的影響Fig.15 The drag AAF change as different motions

4 結 論

基于2?D URANS，研究了4∶1 矩形斷面在不同諧波脈動氣流參數(shù)及結構運動狀態(tài)下的順風向和豎風向復氣動導納函數(shù)，根據(jù)數(shù)值計算結果可得出以下結論：

（1）提出的模擬方法能有效在CFD 內產生幅值可控的豎向和順風向諧波脈動風場，測試值與目標值基本一致。

（2）在0°風攻角下，識別的4∶1 矩形斷面升力復氣動導納與風洞試驗值吻合好，表明該數(shù)值方法可應用于鈍體斷面復氣動導納的識別；較大的風攻角和諧波幅值均會對復氣動導納產生影響，根本的原因是過大的瞬時風攻角使氣流在矩形頂部發(fā)生分離后不再附著，從而引起氣動特性的改變。

（3）在多頻率組合和單頻率諧波風場下識別的升力與阻力氣動導納基本一致，驗證了不同頻率的抖振力具有較好的可疊加性。

（4）當豎向和扭轉簡諧運動引起的風攻角在±1°時，升力氣動導納與斷面靜止時基本一致，在顫抖振分析中將該斷面的抖振升力與自激力分開處理在諧波風場下是可行的。

（5）當簡諧運動引起的風攻角在±1°范圍內，豎向簡諧運動對阻力復氣動導納影響不大，但當折減頻率fD/U>0.04 后，扭轉簡諧運動會增大阻力氣動導納的絕對值。