唐 煜，鄭史雄，張龍奇，馬存明，谷紅強(qiáng)

(西南交通大學(xué)風(fēng)工程研究中心，四川成都 610031)

橋梁斷面氣動導(dǎo)納的數(shù)值識別方法研究

唐 煜，鄭史雄*，張龍奇，馬存明，谷紅強(qiáng)

(西南交通大學(xué)風(fēng)工程研究中心，四川成都 610031)

針對一種氣動導(dǎo)納的數(shù)值識別方法進(jìn)行研究。基于二維不可壓縮URANS方法，選用SST k-ω湍流模型，通過在來流中給定單一頻率的豎向諧波速度分量，計算相應(yīng)的橋梁斷面氣動力荷載時程，識別橋梁斷面的氣動導(dǎo)納。首先考查來流脈動特性在計算域內(nèi)的自保持能力，隨后再對平板和橋梁斷面的氣動導(dǎo)納進(jìn)行識別，所得結(jié)果與理論解和試驗(yàn)值相比較，并討論流場初始化條件的影響。結(jié)果表明:足夠小的網(wǎng)格尺寸和時間步長是來流脈動不發(fā)生明顯衰減的必要條件;平板的氣動導(dǎo)納識別結(jié)果與Sears函數(shù)高度吻合;數(shù)值識別的橋梁斷面升力氣動導(dǎo)納在低頻段與Sears函數(shù)一致，在高頻段略低，但與試驗(yàn)值較接近;力矩氣動導(dǎo)納與Sears函數(shù)有較大差異，但與試驗(yàn)值基本吻合;流場初始化條件對計算效率有影響。

橋梁斷面;氣動導(dǎo)納;數(shù)值方法;URANS;來流自保持;流場初始化條件

0 引言

氣動導(dǎo)納是描述自然風(fēng)脈動速度到結(jié)構(gòu)物所受氣動力的傳遞函數(shù)，是橋梁風(fēng)致抖振響應(yīng)分析中的關(guān)鍵因素之一。氣動導(dǎo)納最初誕生于航空領(lǐng)域，20世紀(jì)30年代Von Karman［1］和Sears［2］在解析機(jī)翼抖振力時引入了氣動導(dǎo)納函數(shù)，基于勢流理論推導(dǎo)得到薄翼氣動導(dǎo)納函數(shù)的理論解。隨后Liepmann［3］將Sears函數(shù)進(jìn)一步近似化簡為更易使用的形式。氣動導(dǎo)納與物體形狀及自然風(fēng)湍流特性有關(guān)，早期的研究對象為薄翼或平板等流線體，這使得研究者能夠從流體力學(xué)基本理論出發(fā)尋找氣動導(dǎo)納的解析解。對于復(fù)雜鈍體繞流問題，由于存在流體理論描述的困難，只能通過試驗(yàn)手段獲得其氣動導(dǎo)納。Davenport［4］給出了一定厚度矩形斷面的氣動導(dǎo)納經(jīng)驗(yàn)式。Holms［5］通過全橋氣彈模型試驗(yàn)研究了West Gate橋的抖振響應(yīng)，指出主梁斷面氣動導(dǎo)納大于Sears函數(shù)，并通過間接修正Sears函數(shù)中的部分系數(shù)，使其可用于主梁斷面的氣動導(dǎo)納描述。Walshe［6］基于Wye橋首次對橋梁斷面氣動導(dǎo)納開展試驗(yàn)間接測量。此后Larose、Scanlan、Diana、Costa、馬存明、王雄江等大量研究者［7-12］對橋梁斷面氣動導(dǎo)納進(jìn)行了卓有成效的探索，發(fā)展出各種各樣的識別方法，這些方法大致可歸為兩類，其一是同時測得脈動風(fēng)譜與抖振力譜，再依二者比值直接擬合出氣動導(dǎo)納經(jīng)驗(yàn)表達(dá)式;其二是不測試氣動力，測試脈動風(fēng)譜與氣彈模型隨機(jī)響應(yīng)譜，再以系統(tǒng)辨識方法間接識別氣動導(dǎo)納，而這些研究均建立在風(fēng)洞試驗(yàn)基礎(chǔ)之上。

從試驗(yàn)研究所涉及的風(fēng)場來看，識別氣動導(dǎo)納所用風(fēng)場條件常見有兩種。一為頻率連續(xù)分布的局部各向同性紊流風(fēng)場(寬帶紊流)，湍流能譜滿足Kolmogorov“-5/3”定律，在該風(fēng)場中識別氣動導(dǎo)納可一次性獲得各個頻率下氣動導(dǎo)納值，此類風(fēng)場可由風(fēng)洞中固定的被動湍流模擬裝置(如尖劈、格柵、粗糙塊)制造產(chǎn)生，但目標(biāo)風(fēng)場的低頻和高頻成分湍動能往往不易滿足，使得一次識別全頻率氣動導(dǎo)納在頻率兩端特別是高頻段易導(dǎo)致較大誤差。另一類為頻率單一的諧波風(fēng)場(窄帶單頻紊流)，在該風(fēng)場中識別氣動導(dǎo)納一般采用分離頻率法，通過改變風(fēng)場卓越頻率進(jìn)行多次識別，能獨(dú)立保證高低頻段各自的識別精度，該類風(fēng)場只能在主動控制風(fēng)洞(多風(fēng)扇、主動葉柵)中獲得，是當(dāng)下研究的熱點(diǎn)。需指出的是，由于湍流雷諾應(yīng)力的能量再分配作用，初始時僅有單向速度脈動成分的一維湍流，會隨著流動向下游的推移演化逐漸發(fā)展成具有三個方向速度脈動的三維湍流，實(shí)驗(yàn)室內(nèi)生成的窄帶單頻紊流不可能保持嚴(yán)格意義上的單一頻率，脈動頻率必分布在一定頻段范圍內(nèi)，在此風(fēng)場中識別的氣動導(dǎo)納結(jié)果是有待商榷的。

風(fēng)洞試驗(yàn)中要測試獲得模型的氣動力有兩種途徑，一是用天平直接連接模型測力，二是先測量表面壓力再積分求合力。前者對模型制作要求十分苛刻，模型測力系統(tǒng)需具備足夠高的頻率(高于80 Hz)，模型在測試過程中不能發(fā)生變形和振動，否則將面臨模型慣性力難剔除的困難。后者僅適用于閉口箱梁斷面的氣動導(dǎo)納測試，無法考慮欄桿、軌道等橋梁附屬設(shè)施的影響，且測試所得氣動力還受測壓孔布置方式及壓力積分格式等多重因素影響嚴(yán)重，誤差較大。

現(xiàn)場實(shí)測是研究真實(shí)橋梁氣動導(dǎo)納的最佳手段。陶奇［13］采用測壓法現(xiàn)場實(shí)測研究蘇通大橋主梁斷面氣動導(dǎo)納，獲取第一手珍貴實(shí)測數(shù)據(jù)，研究結(jié)果表明該橋現(xiàn)場實(shí)測氣動導(dǎo)納結(jié)果與風(fēng)洞試驗(yàn)結(jié)果接近，與Sears函數(shù)有交叉。使用測壓法的現(xiàn)場實(shí)測結(jié)果同樣受測點(diǎn)布置和積分格式的影響，此外還要克服測壓管路過長引起的各測點(diǎn)信號不同步及脈動信號失真的困難。就本文作者收集的文獻(xiàn)資料來看，目前暫未見到其他有關(guān)橋梁斷面氣動導(dǎo)納現(xiàn)場實(shí)測的公開報道。

正因?yàn)殁g體擾流理論推導(dǎo)面臨的障礙難以逾越，試驗(yàn)研究亦存在不少問題，現(xiàn)場實(shí)測研究鳳毛麟角，近年來有少數(shù)學(xué)者開始嘗試借助CFD數(shù)值模擬技術(shù)識別氣動導(dǎo)納。與風(fēng)洞試驗(yàn)不同，數(shù)值識別中可直接生成完全單一頻率的諧波風(fēng)場，數(shù)值模型提取氣動力簡單便捷，研究周期短費(fèi)用低，潛力巨大。Uejima［14］基于二維雷諾平均的CFD數(shù)值模擬生成單頻諧波風(fēng)場，研究平板、長寬比5∶1的矩形和扁平六邊形斷面的氣動導(dǎo)納，平板模擬結(jié)果與Sears函數(shù)十分接近，另外兩個斷面模擬結(jié)果也與試驗(yàn)吻合較好，但Uejima對計算來流特性的自保持問題討論得不夠，使得其識別的折算頻率最大不超過0.7。韓艷［15］推導(dǎo)了6個復(fù)氣動導(dǎo)納定義式，并基于三維雷諾平均方法采用標(biāo)準(zhǔn)κ-ε模型研究平板的氣動導(dǎo)納，所得結(jié)果與試驗(yàn)值和Sears函數(shù)存在差別，其識別的折算頻率最大不超過0.6。Rasmussen［16］基于二維離散渦方法生成頻率連續(xù)的紊流風(fēng)場，研究平板的氣動導(dǎo)納結(jié)果同Liepmann近似解較吻合，但在高頻段(折算頻率大于2)吻合度不太理想，其還研究大貝爾特東橋主梁的氣動導(dǎo)納，數(shù)值所得氣動導(dǎo)納同試驗(yàn)結(jié)果整體基本吻合，而高頻段二者差異顯著。總的來說，橋梁斷面氣動導(dǎo)納的數(shù)值識別方法尚處于起步階段，而基于雷諾平均CFD數(shù)值模擬識別橋梁斷面氣動導(dǎo)納的研究則鮮有報道。

橋梁斷面氣動導(dǎo)納數(shù)值識別的主要工作在于斷面繞流場CFD計算，故須對流場計算中常見的幾個問題予以足夠關(guān)注。本文在 Uejima［14］工作的基礎(chǔ)上，首先對無障礙物計算域的入流自保持能力進(jìn)行詳細(xì)討論，使其具備更廣的折算頻率識別范圍，隨后開展理想平板的氣動導(dǎo)納數(shù)值識別以驗(yàn)證方法，緊接著研究某橋梁斷面的氣動導(dǎo)納并與風(fēng)洞試驗(yàn)比較，最后還分析了流場初始化技術(shù)對數(shù)值計算效率的影響。與傳統(tǒng)與風(fēng)洞試驗(yàn)和現(xiàn)場實(shí)測相比，本文數(shù)值方法所生成的脈動風(fēng)場具有單一頻率，且在獲取氣動力的便捷性上亦有優(yōu)勢，為工程中識別橋梁斷面氣動導(dǎo)納建議了一種新的途徑。

1 數(shù)值模型

1.1 基本控制方程

假設(shè)橫向風(fēng)繞過橋梁斷面為二維不可壓縮的流動過程，以非定常雷諾平均納維-斯托克斯方程(URANS)為基本流動控制方程:

模型封閉系數(shù):

α=5/9，β=3/40，β*=9/100，σ=1/2，σ*= 1/2 時間離散采用二階隱式，對流項(xiàng)為二階迎風(fēng)格式，其它流動物理量的空間離散也為二階格式，以SIMPLEC算法處理壓強(qiáng)與速度耦合，數(shù)值求解在基于有限體積法的通用流體軟件Fluent中完成。

1.2 計算域和邊界條件

模型計算域見圖1，在計算域左側(cè)邊界為速度入口條件，x方向速度保持恒定不變，y方向速度隨時間簡諧變化，通過用戶自定義函數(shù)(UDF)編程實(shí)現(xiàn)。上下側(cè)邊界也為速度入口條件，y方向速度不僅隨時間簡諧變化，還隨著空間位置不同而變化，以保持與左側(cè)邊界的波動特征相協(xié)調(diào)，其余邊界參數(shù)設(shè)置與左側(cè)邊界類似。右側(cè)邊界為自由出流條件。橋梁斷面為無滑移壁面條件。

圖1 計算域Fig.1 Computational domain

1.3 基本參數(shù)和網(wǎng)格設(shè)計

數(shù)值模擬建模時未直接使用實(shí)際橋梁斷面的原始尺寸，而采用風(fēng)洞試驗(yàn)中常見的縮尺橋梁模型斷面尺寸(梁寬B=0.5 m～1 m)，并據(jù)此確定計算域大小，如圖1所示。該處理基于以下考慮:①氣動導(dǎo)納與橋梁斷面形狀有關(guān)，而受雷諾數(shù)影響較?。?3］;②現(xiàn)場實(shí)測極為少見，數(shù)值方法研究足尺斷面氣動導(dǎo)納缺乏實(shí)測數(shù)據(jù)支持;③足尺建模的計算域太大，二維網(wǎng)格數(shù)量隨計算域尺寸增大呈平方關(guān)系增長，計算量難以承受。

來流平均風(fēng)速U∞取為12 m/s，y向(豎向)速度簡諧脈動幅值A(chǔ)取0.336，相應(yīng)的湍流強(qiáng)度為2%。通常，來流的湍流強(qiáng)度會影響氣動導(dǎo)納的識別結(jié)果，但由于本文重點(diǎn)旨在介紹氣動導(dǎo)納的數(shù)值識別方法，對湍流強(qiáng)度的影響，未作詳細(xì)討論。

對無障礙物的計算域，按均勻方形結(jié)構(gòu)化網(wǎng)格對其進(jìn)行劃分。對有橋梁斷面的計算域網(wǎng)格劃分如圖3，保持外圍區(qū)域?yàn)榫鶆蚍叫谓Y(jié)構(gòu)化網(wǎng)格不變，僅在斷面附近區(qū)域內(nèi)使用四邊形非結(jié)構(gòu)化網(wǎng)格，斷面壁面網(wǎng)格局部加密以滿足湍流模型對壁面網(wǎng)格y+值的要求，并通過尺寸函數(shù)控制漸變網(wǎng)格大小從壁面到均勻區(qū)域的平緩漸變，其相應(yīng)的相鄰網(wǎng)格尺寸增長因子為1.06。

2 入流自保持能力

正確的數(shù)值方法必須和流動的物理性質(zhì)一致。通常，在外部繞流的CFD數(shù)值模擬中，隨時空變化的計算入流邊界條件須具備在無障礙物的整個計算域中自然保持其邊界特征的能力，應(yīng)避免發(fā)生明顯不符合計算者預(yù)期的幅值衰減(即數(shù)值耗散)和頻率畸變(即數(shù)值彌散/色散)，如雷諾平均(RANS)中大氣邊界層入流風(fēng)剖面自保持［18-19］，又如大渦模擬(LES)中充分發(fā)展均勻湍流脈動特性自保持［20］。

在對橋梁斷面氣動導(dǎo)納進(jìn)行數(shù)值識別之前，同樣需考察無障礙物計算模型入口處簡諧速度分量的自保持能力，其中計算域網(wǎng)格分辨率和數(shù)值求解使用的非定常時間步長是關(guān)鍵影響參數(shù)。

2.1 網(wǎng)格分辨率的影響

以往鈍體繞流的數(shù)值模擬中有研究者［21］結(jié)合自身經(jīng)驗(yàn)給出了鈍體網(wǎng)格分辨率的經(jīng)驗(yàn)值，且大多以壁面網(wǎng)格份數(shù)這種無量綱的形式表達(dá)。但他們重點(diǎn)關(guān)心研究對象的氣動力或壓強(qiáng)分布，其給出的網(wǎng)格分辨率建議值可能不適合本文時變來流脈動特征的自保持問題。在此從數(shù)值誤差角度，對網(wǎng)格尺寸進(jìn)行經(jīng)驗(yàn)估計并驗(yàn)證。

一般認(rèn)為數(shù)值誤差主要產(chǎn)生于把連續(xù)的微分方程轉(zhuǎn)換為各網(wǎng)格點(diǎn)處離散的代數(shù)方程的過程中，并隨著對流擴(kuò)散而存在于整個計算流場中，在時間步長和數(shù)值格式確定的情況下，離散誤差的大小主要取決于網(wǎng)格尺寸。

附加在主流上的豎向簡諧脈動速度對應(yīng)的空間波長為:

式中U∞為計算來流平均風(fēng)速，f為豎向脈動頻率。

要使該簡諧波長能夠在計算域內(nèi)有效傳播，則網(wǎng)格尺寸Δx應(yīng)至少比波長小一個量級:

其中，B為待識別斷面的寬度，N為單波長內(nèi)網(wǎng)格數(shù)，需進(jìn)一步研究確定。

在無障礙物空計算域內(nèi)試算，改變網(wǎng)格尺寸以考察y向簡諧脈動速度在計算域內(nèi)的自保持情況。待計算達(dá)到統(tǒng)計穩(wěn)定后，順主流方向沿直線y=0任意提取某時刻的脈動速度分量的空間分布，如圖2所示。圖中相鄰兩波峰間幅值衰減可用對數(shù)衰減率δ表示:

本文在此處約定，當(dāng)δ≤0.003時，認(rèn)為數(shù)值耗散可接受，來流脈動特征實(shí)現(xiàn)自保持。

兼顧計算效率的同時，應(yīng)盡可能地減小時間離散帶來的誤差。試算時采用的時間步長為0.000 1 s，采樣頻率10 000 Hz，而通常氣動導(dǎo)納研究考查的最大折算頻率約為10，若按斷面寬度0.6 m計，對應(yīng)的簡諧速度脈動頻率為31 Hz，采樣頻率約為最大簡諧速度脈動頻率的300倍，時間離散誤差的影響可忽略。

圖2 豎向脈動衰減示意Fig.2 Decay of the vertical velocity

在不同折算頻率下進(jìn)行數(shù)值試算，考查網(wǎng)格尺寸對來流脈動速度分量在計算域內(nèi)自保持能力的影響，結(jié)果見圖3。豎向簡諧脈動速度的衰減受網(wǎng)格尺寸影響顯著，隨著網(wǎng)格數(shù)量的倍增，脈動速度幅值的對數(shù)衰減率與之成平方關(guān)系減小，且兩個折算頻率下的衰減規(guī)律是一致的。

圖3 單波網(wǎng)格數(shù)對豎向脈動速度衰減的影響Fig.3 Decay of vertical velocity influenced by N

通過數(shù)據(jù)擬合，可得對數(shù)衰減率經(jīng)驗(yàn)表達(dá)式:

圖3中按該式繪制的曲線與兩個折算頻率下的試算結(jié)果均吻合良好，說明用該經(jīng)驗(yàn)式近似描述流場豎向簡諧速度脈動衰減規(guī)律是可行的。

當(dāng)單波網(wǎng)格數(shù)N達(dá)80時，對數(shù)衰減率不超過0.003。在進(jìn)行氣動導(dǎo)納數(shù)值識別的流場網(wǎng)格劃分時，建議按式(8)確定計算域最大網(wǎng)格尺寸，且式中單波網(wǎng)格數(shù)N應(yīng)不小于80，可基本保證來流脈動速度的自保持。

2.2 時間步長的影響

取單波網(wǎng)格數(shù)N為80，考查時間步長對簡諧脈動速度幅值的影響?；趤砹魉俣群陀嬎阌蚓W(wǎng)格尺寸，定義一個無量綱時間步長

表1所示試算結(jié)果反映時間步長對簡諧脈動速度幅值的影響。時間步長越小，脈動速度的衰減也越小，當(dāng)無量綱時間步長小至1.2時，計算來流脈動速度幅值的對數(shù)衰減率可控制在0.003以下，高、低折算頻率條件下試算得到的衰減規(guī)律基本一致。在進(jìn)行氣動導(dǎo)納識別時，建議無量綱時間步長取值應(yīng)不大于1.2，據(jù)此估算相應(yīng)計算工況的時間步長。另外，本研究試算中還觀察到網(wǎng)格分辨率和時間步長對脈動速度頻率無明顯影響，不再贅述。

表1 時間步長對簡諧脈動速度幅值的影響Table 1 Decay of vertical velocity influenced by time step

3 氣動導(dǎo)納識別

在數(shù)值識別氣動導(dǎo)納時，斷面升力的功率譜表示為:

其中，k為式(7)定義的折算頻率，|χL(k)|2為升力的氣動導(dǎo)納，|χM(k)|2為力矩的氣動導(dǎo)納，SL(fr)為升力頻譜，SM(fr)為力矩頻譜，CD、CL、CM分別為阻力、升力和力矩系數(shù)，C'L、C'M分別為升力系數(shù)和力矩系數(shù)對風(fēng)迎角的斜率，ρ為空氣密度，U∞為來流平均風(fēng)速，B為斷面寬度。

升力氣動導(dǎo)納的識別過程如下:①在均勻流條件下，通過CFD計算獲得C'L+CD值;②生成具有單一折算頻率簡諧脈動分速度的風(fēng)場，計算獲得該風(fēng)場中待研究斷面的升力時程;③對升力時程進(jìn)行快速傅里葉變換獲得頻譜數(shù)據(jù)，代入式(11)計算該折算頻率下的氣動導(dǎo)納值;④改變折算頻率，重復(fù)步驟②③，直至在所考察折算頻率范圍內(nèi)獲得足夠多的數(shù)據(jù)點(diǎn)，以擬合氣動導(dǎo)納函數(shù)。同理可識別力矩的氣動導(dǎo)納。

因本文重在介紹氣動導(dǎo)納的數(shù)值識別方法，故未對氣動導(dǎo)納理論如三維氣動導(dǎo)納［11］或六分量導(dǎo)納［15］等做深入研究。

3.1 平板氣動導(dǎo)納

對于理想平板升力氣動導(dǎo)納的Sear函數(shù)，Liepmann給出的理論近似解為:

數(shù)值計算中取寬厚比為100的矩形平板進(jìn)行建模，為減小鈍角引起的流動分離，對平板前后緣倒大圓角處理。對于平板，其阻力系數(shù)為零，升力系數(shù)斜率為2π。

圖4為平板升力氣動導(dǎo)納的識別結(jié)果，并與近似解析解相比較。不難看出，無論是在高折算頻率還是在低折算頻率，平板氣動導(dǎo)納的數(shù)值識別結(jié)果與按勢流理論推導(dǎo)得到的Sears函數(shù)都吻合得非常好。

圖4 平板的氣動導(dǎo)納Fig.4 Aerodynamic admittance of thin plate

3.2 箱梁氣動導(dǎo)納

南京三橋(NJ3B)主梁為一典型扁平鋼箱梁，寬高比11.6，本文作者之一曾對該橋主梁斷面氣動導(dǎo)納進(jìn)行了詳細(xì)的風(fēng)洞模型試驗(yàn)測試［11］，試驗(yàn)?zāi)Ｐ蛿嗝嫒鐖D5所示。在此以該斷面為例，驗(yàn)證數(shù)值方法對橋梁斷面氣動導(dǎo)納的識別能力。

圖5 南京三橋主梁模型斷面(mm)Fig.5 Section of NJ3B(unit:mm)

在均勻來流條件下，先識別模型斷面的三分力系數(shù)，結(jié)果見圖6，其中力矩方向以圖1中順時針為正。零度迎角來流時，計算得到的升力系數(shù)斜率為5.294，與文獻(xiàn)［11］中的試驗(yàn)值5.302十分吻合;力矩系數(shù)斜率為1.532，與試驗(yàn)值1.305比較接近;計算的阻力系數(shù)為0.213。這些數(shù)據(jù)作為基本參數(shù)，代入式(11)、式(12)用于斷面氣動導(dǎo)納的識別。

南京三橋主梁斷面的氣動導(dǎo)納識別結(jié)果見圖7，其中用于結(jié)果比較的試驗(yàn)測試風(fēng)場為固定格柵制造的紊流，豎向積分尺度為0.41 m。

圖6 南京三橋主梁靜力三分力系數(shù)Fig.6 Static coefficients of NJ3B

圖7 南京三橋斷面的氣動導(dǎo)納Fig.7 Aerodynamic admittance of NJ3B

就升力氣動導(dǎo)納而言，扁平箱梁斷面的數(shù)值識別結(jié)果與Sears函數(shù)相比，在低頻范圍內(nèi)非常接近，在高折算頻率段略低，但與已有試驗(yàn)結(jié)果在折算頻率(0.3，4)范圍內(nèi)亦吻合良好。

就力矩氣動導(dǎo)納而言，扁平箱梁斷面的數(shù)值識別結(jié)果與Sears函數(shù)存在明顯差異。原因應(yīng)在于:Sears函數(shù)描述的對象為0°風(fēng)迎角下的類平板薄翼，當(dāng)迎角小范圍變化時，其氣動中心是基本固定的，而對于截面形式非對稱的扁平箱梁斷面，并不具備固定的氣動中心，造成了二者在結(jié)果比較上的困難。但數(shù)值識別的結(jié)果與已有試驗(yàn)在高折算頻率段卻有著很好的吻合度。

對于更高的折算頻率范圍(k＞4)，未做更進(jìn)一步的計算驗(yàn)證主要基于如下兩個原因:①本文識別氣動導(dǎo)納折算頻率范圍(0.01，4)已能滿足常規(guī)抖振計算需要。橋梁斷面氣動導(dǎo)納識別主要用于大跨度橋梁抖振計算，以橋面寬30 m為例，一般大跨度橋梁常見的設(shè)計風(fēng)速值約為35 m/s，主梁一階豎彎頻率約為0.2～0.4 Hz，則其對風(fēng)致振動敏感的折算頻率在0.17～0.34附近，本文所識別氣動導(dǎo)納的折算頻率范圍已將其包含在內(nèi)。②計算資源的限制。以本文識別的最高折算頻率k=4為例，為保證來流的自維持，二維模型的網(wǎng)格數(shù)量需達(dá)400萬，若要識別k= 10時的氣動導(dǎo)納值，單向網(wǎng)格尺寸需細(xì)化至該模型的2/5，由于數(shù)值模型的兩個維度，則計算模型網(wǎng)格總數(shù)將達(dá)2 500萬，普通工程計算難以承受。

3.3 流場初始化條件的討論

在數(shù)值模型建立后，需預(yù)先人為給定流場初始值作為迭代計算的初始化條件。合理的初始化流場能幫助數(shù)值計算過程順利收斂，提高計算效率，故流場的初始化條件是值得研究的。

對鈍體空氣繞流問題，大多數(shù)研究者習(xí)慣參考來流邊界條件的參數(shù)對所有流場內(nèi)點(diǎn)進(jìn)行賦值，所得初始流場物理量處處相等，為均勻常數(shù)場。當(dāng)入口邊界為定常條件時，這樣的常數(shù)初始場是非常合理的，而當(dāng)入口邊界為非定常時，由于入口的時變信息傳播至模型處需要一定的計算時間，常數(shù)初始場不一定為最優(yōu)，這時可考慮非均勻的初始化條件。

通過Fluent軟件預(yù)留的用戶自定義函數(shù)(UDF)接口自編程序，以計算域上下邊界處(圖1)的波動方程對全場速度進(jìn)行流場初始化。湍動能和比耗散率仍參考入流邊界取常數(shù)值均勻初始化。

對來流自保持問題中無障礙物計算域進(jìn)行數(shù)值模擬，監(jiān)視計算域(6B，0)點(diǎn)處y向速度時程，對比兩種初始化條件對計算過程的影響，如圖8所示。當(dāng)采用均勻條件初始化流場時，計算初級監(jiān)視點(diǎn)處y向速度基本為零，經(jīng)過一段時間后開始進(jìn)入等幅簡諧脈動?？蓪υ摃r間長度進(jìn)行估算:由邊界時變信息的傳播速度為U∞/f，入口邊界到計算關(guān)心區(qū)域的距離為6B，則傳播用時6Bf/U∞。這段計算時間對于氣動導(dǎo)納的識別沒有實(shí)際意義，僅是流場初始條件不夠好所造成的額外迭代過程耗費(fèi)。當(dāng)采用恰當(dāng)?shù)姆蔷鶆虺跏蓟瘲l件，可快速獲得流場穩(wěn)態(tài)振蕩解。

圖9所示為k=2時兩種初始條件下的南京三橋斷面升力曲線。從圖中曲線來看，隨著計算時間的推移，兩種初始條件下計算的斷面升力最終都將進(jìn)入周期振蕩，初始條件的影響主要表現(xiàn)在數(shù)值迭代開始時解的振蕩差異上。氣動導(dǎo)納識別所用升力時程須選取穩(wěn)定振蕩部分的數(shù)據(jù)，從這個意義上看，兩種初始化條件并無明顯的效率差別。造成該現(xiàn)象的原因在于:非均勻初始化技術(shù)主要是節(jié)約來流邊界時變信息傳播到數(shù)值計算關(guān)心區(qū)域的時間t1，而對由橋梁斷面附近因斷面存在所造成的有物理意義的流場收斂尚需時間t2演化，當(dāng)t1不顯著大于t2時，流場非均勻初始化的優(yōu)勢并不十分明顯。由此還可推斷，當(dāng)待識別斷面距離入口邊界較遠(yuǎn)或來流速度很小時，使用非均勻初始化條件將在一定程度上提高識別效率。

圖8 監(jiān)視點(diǎn)的速度分量時程Fig.8 Time history of the vertical velocity for monitored point

圖9 斷面升力時程Fig.9 Lift history of NJ4B deck

4 結(jié)論

本文考查一種基于二維URANS的橋梁斷面氣動導(dǎo)納數(shù)值識別方法，得到如下結(jié)論:

(1)來流的豎向簡諧脈動速度在計算域內(nèi)的自保持能力受網(wǎng)格分辨率和時間步長的影響，單波長內(nèi)的網(wǎng)格數(shù)越大，對應(yīng)的無量綱時間步長越小，該簡諧脈動速度幅值衰減程度越小。

(2)提出了豎向簡諧脈動速度幅值衰減效應(yīng)的估算公式，據(jù)此確定數(shù)值模型最大網(wǎng)格尺寸與時間步長，可保證來流脈動速度幅值不發(fā)生明顯衰減，這也是氣動導(dǎo)納識別的前提條件。

(3)平板升力氣動導(dǎo)納的數(shù)值識別結(jié)果與Sears解析解高度吻合，證明了數(shù)值方法的正確性和可行性。橋梁斷面升力氣動導(dǎo)納數(shù)值識別結(jié)果在低頻段與Sears解吻合較好，在高頻段雖低于Sears解，但與既有文獻(xiàn)中的風(fēng)洞試驗(yàn)結(jié)果基本一致，為本文數(shù)值方法的工程應(yīng)用建立了信心。

(4)流場初始條件對計算效率確有影響?？傮w上看，數(shù)值識別時使用恰當(dāng)?shù)姆蔷鶆虺跏紬l件要比使用均勻初始條件有利于提高效率。

建議進(jìn)一步開展順風(fēng)向簡諧脈動速度自保持能力研究，以及橋梁、高層建筑等結(jié)構(gòu)斷面阻力氣動導(dǎo)納的數(shù)值識別工作。

［1］Karman T H V，Sears W R.Airfoil theory for non-uniform motion［J］.Journal of the Aeronautical Sciences，1938，5(10):379-390.

［2］Sears W R.Some aspects of non-stationary airfoil theory and its practical application［J］.Journal of the Aeronautical Sciences，1941，8(3):104-108.

［3］Liepmann H W.On the application of statistical concepts to the buffeting problem［J］.Journal of the Aeronautical Science，1952，19 (2):793-800.

［4］Davenport A G.Buffeting of a suspension bridge by storm winds［J］.Journal of Structural Engineering，1962，88(ST3):233-268.

［5］Holmes J D.Prediction of the response of a cable-stayed bridge to turbulence［C］.Proceedings of the Fourth International Conference on Wind Effects on Buildings and Structures，Heathrow，1975，187-197.

［6］Walshe D E，Wyatt T A.Measurement and application of the aerodynamic admittance function for a box-girder bridge［J］.Journal of Wind Engineering and Industrial Aerodynamics，1983，14(1): 211-222.doi:10.1016/0167-6105(83)90024-7

［7］Larose G L，Mann J.Gust loading on streamlined bridge decks［J］.Journal of Fluids and Structures，1998，12(5):511-536.doi:10.1006/jfls.1998.0161

［8］Scanlan R H，Jones N P.A form of aerodynamic admittance for use in bridge aeroelastic analysis［J］.Journal of Fluids and Structures，1999，13(7):1017-1027.doi:10.1006/jfls.1999.0243

［9］Diana G，Bruni S，Cigada A，et al.Complex aerodynamic admittance function role in buffeting response of a bridge deck［J］.Journal of Wind Engineering and Industrial Aerodynamics，2002，90 (12):2057-2072.doi:10.1016/S0167-6105(02)00321-5

［10］Costa c.Aerodynamic admittance functions and buffeting forces for bridges via indicial functions［J］.Journal of Fluids and Structures，2007，23(3):413-428.doi:10.1016/j.jfluidstructs.2006.10.002

［11］Ma Cunming.3D aerodynamic admittance research of streamlined box bridge decks［D］.［PhD Thesis］.Chengdu:Southwest Jiaotong University，2007.(in Chinese)馬存明.流線箱型橋梁斷面三維氣動導(dǎo)納研究［D］.［博士學(xué)位論文］.成都:西南交通大學(xué)，2007.

［12］Wang Xiongjiang，Gu Ming，Qin Xianrong.Estimation of aerodynamic admittance function of long-span bridges by the relationships among aerodynamic parameters［J］.Acta Aerodynamica Sinica，2009，27(6):707-712.(in Chinese)王雄江，顧明，秦仙蓉.基于氣動參數(shù)之間關(guān)系的橋梁斷面氣動導(dǎo)納識別［J］.空氣動力學(xué)學(xué)報，2009，27(6):707-712.

［13］Tao Qi，Liao Haili，Li Mingshui，et al.The actually test research of aerodynamic admittance of long-span cable-stayed bridge girder section［J］.Acta Aerodynamica Sinica，2010，28(006):655-660.(in Chinese)陶奇，廖海黎，李明水，等.大跨斜拉橋主梁斷面氣動導(dǎo)納的實(shí)測研究［J］.空氣動力學(xué)學(xué)報，2010，28(006):655-660.

［14］Uejima H，Kuroda S，Kobayashi H.Estimation of aerodynamic admittance by numerical computation［C］.BBAAⅥ International Colloquium on Bluff Bodies Aerodynamics and Applications，Milano，Italy，July，20-24 2008

［15］Han Yan，Chen Zhengqing.Experimental and numerical simulations studies on complex aerodynamic admittance functions of thin plate section［J］.Journal of Vibration Engineering，2009，22(2):200-206.(in Chinese)韓艷，陳政清.薄平板復(fù)氣動導(dǎo)納函數(shù)的試驗(yàn)與數(shù)值模擬研究［J］.振動工程學(xué)報，2009，22(2):200-206.

［16］Rasmussen J T，Hejlesen M M，Larsen A，et al.Discrete vortex method simulations of the aerodynamic admittance in bridge aerodynamics［J］.Journal of Wind Engineering and Industrial Aerodynamics，2010，98(12):754-766.doi:10.1016/j.jweia.2010.06.011

［17］Menter F R.Two-equation eddy-viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32(8):1598-1605.doi:10.2514/3.12149

［18］Richards P J，Younis B A.Comments on“Prediction of wind-generated pressure distribution around buildings”by E.H.Mathews［J］.Journal of Wind Engineering and Industrial Aerodynamics， 1990，34(1):107-110.doi:10.1016/0167-6105(90)90152-3

［19］Yang Y，Gu M，chen S，et al.New inflow boundary conditions for modelling the neutral equilibrium atmospheric boundary layer in computational wind engineering［J］.Journal of Wind Engineering and Industrial Aerodynamics，2009，97(2):88-95.doi:10.1016/j.jweia.2008.12.001

［20］Tamura T，Ono Y.LES analysis on aeroelastic instability of prisms in turbulent flow［J］.Journal of Wind Engineering and Industrial Aerodynamics，2003，91(12):1827-1846.doi:10.1016/j.jweia.2003.09.032

［21］Kravchenko A G，Moin P.Numerical studies of flow over a circular cylinder at Re=3900［J］.Physics of Fluids，2000，12(2):403-417.doi:10.1063/1.870318

Numerical method for identifying the aerodynamic admittance of bridge deck

Tang Yu，Zheng Shixiong*，Zhang Longqi，Ma Cunming，Gu Hongqiang
(Research Centre for Wind Engineering，Southwest Jiaotong University，Chengdu 610031，China)

A numerical method for identifying the aerodynamic admittance of bridge decks was proposed.A two-dimensional incompressible unsteady Reynolds average Navier-Stokes(URANS)approach with SST k-ω turbulence model was used.By generating an incoming flow with a smooth longitudinal component and a vertical single-frequency sinusoidal component velocity，the history of forces acting on the bridge deck was calculated through flow solution.Thus the aerodynamic admittance could be obtained by analyzing its spectrum.First of all，the self-sustaining capability of the generated inflow was studied and confirmed.Then a thin plate and a bridge deck section were taken to verify the method by comparing to the theory resolution and the experiment value.Finally，the initial condition for flow solution was discussed.The results showed that both the grid size and the time step should be small enough to keep the decay of the vertical sinusoidal component in an insignificant extent.The identified aerodynamic admittance of the thin plate got a quit perfect agreement with Sears function.The identified lift aerodynamic admittance of bridge deck was in accordance with Sears function in low reduced frequency，and less in high frequency，while it was closer to the experimental value.The identified moment aerodynamic admittance of bridge deck was also close to the experimental though it seemed different from Sears function.The initial condition of flow field had an effect on efficiency.

bridge deck;aerodynamic admittance;numerical method;URANS;self-sustaining inflow;initial conditions of flow field

TU973+.32;U441+.3

:Adoi:10.7638/kqdlxxb-2014.0079

0258-1825(2015)05-0706-08

2014-07-31;

:2014-10-10

國家自然科學(xué)基金(51378443)

唐煜(1987-)，男，湖南常德人，博士研究生，從事計算風(fēng)工程研究.E-mail:tang0107@163.com

鄭史雄*(1965-)，男，浙江江山人，教授，博導(dǎo)，從事橋梁結(jié)構(gòu)抗風(fēng)抗震研究.E-mail:zhengsx@swjtu.edu.cn

唐煜，鄭史雄，張龍奇，等.橋梁斷面氣動導(dǎo)納的數(shù)值識別方法研究［J］.空氣動力學(xué)學(xué)報，2015，33(5):706-713.

10.7638/kqdlxxb-2014.0079 Tang Y，Zheng S X，Zhang L Q，et al.Numerical method for identifying the aerodynamic admittance of bridge deck［J］.Acta Aerodynamica Sinica，2015，33(5):706-713.