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

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

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

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

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

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

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

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

0 引言

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

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

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

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

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

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

1 數(shù)值模型

1.1 基本控制方程

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

模型封閉系數(shù):

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

1.2 計(jì)算域和邊界條件

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

圖1 計(jì)算域Fig.1 Computational domain

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

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

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

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

2 入流自保持能力

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

在對(duì)橋梁斷面氣動(dòng)導(dǎo)納進(jìn)行數(shù)值識(shí)別之前，同樣需考察無障礙物計(jì)算模型入口處簡(jiǎn)諧速度分量的自保持能力，其中計(jì)算域網(wǎng)格分辨率和數(shù)值求解使用的非定常時(shí)間步長(zhǎng)是關(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)心研究對(duì)象的氣動(dòng)力或壓強(qiáng)分布，其給出的網(wǎng)格分辨率建議值可能不適合本文時(shí)變來流脈動(dòng)特征的自保持問題。在此從數(shù)值誤差角度，對(duì)網(wǎng)格尺寸進(jìn)行經(jīng)驗(yàn)估計(jì)并驗(yàn)證。

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

附加在主流上的豎向簡(jiǎn)諧脈動(dòng)速度對(duì)應(yīng)的空間波長(zhǎng)為:

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

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

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

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

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

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

圖2 豎向脈動(dòng)衰減示意Fig.2 Decay of the vertical velocity

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

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

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

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

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

2.2 時(shí)間步長(zhǎng)的影響

取單波網(wǎng)格數(shù)N為80，考查時(shí)間步長(zhǎng)對(duì)簡(jiǎn)諧脈動(dòng)速度幅值的影響?；趤砹魉俣群陀?jì)算域網(wǎng)格尺寸，定義一個(gè)無量綱時(shí)間步長(zhǎng)

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

表1 時(shí)間步長(zhǎng)對(duì)簡(jiǎn)諧脈動(dòng)速度幅值的影響Table 1 Decay of vertical velocity influenced by time step

3 氣動(dòng)導(dǎo)納識(shí)別

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3 流場(chǎng)初始化條件的討論

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

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

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

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

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

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

圖9 斷面升力時(shí)程Fig.9 Lift history of NJ4B deck

4 結(jié)論

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

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

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

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

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

建議進(jìn)一步開展順風(fēng)向簡(jiǎn)諧脈動(dòng)速度自保持能力研究，以及橋梁、高層建筑等結(jié)構(gòu)斷面阻力氣動(dòng)導(dǎo)納的數(shù)值識(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òng)導(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)王雄江，顧明，秦仙蓉.基于氣動(dòng)參數(shù)之間關(guān)系的橋梁斷面氣動(dòng)導(dǎo)納識(shí)別［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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òng)導(dǎo)納的實(shí)測(cè)研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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òng)導(dǎo)納函數(shù)的試驗(yàn)與數(shù)值模擬研究［J］.振動(dòng)工程學(xué)報(bào)，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-)，男，湖南常德人，博士研究生，從事計(jì)算風(fēng)工程研究.E-mail:tang0107@163.com

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

唐煜，鄭史雄，張龍奇，等.橋梁斷面氣動(dòng)導(dǎo)納的數(shù)值識(shí)別方法研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，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.