張偉峰, 張志田, 張顯雄

(湖南大學(xué)土木工程學(xué)院 風(fēng)工程與橋梁工程湖南省重點(diǎn)實(shí)驗(yàn)室, 湖南 長沙 410082)

0 引 言

處于大氣邊界層中的大跨度橋梁，通常會受到靜風(fēng)力、周期性的渦激力和隨機(jī)氣動力的作用。其中抖振力是由空氣中的湍流引起，因此對于實(shí)際橋梁是不可避免的?；跅l帶假定的大跨度橋梁隨機(jī)抖振理論把來流湍流產(chǎn)生的氣動力效應(yīng)按單位長度的斷面氣動力和沿跨度方向的空間相關(guān)性來處理[1]。通常，沿橋梁跨度方向的相關(guān)效應(yīng)用聯(lián)合接收函數(shù)來表示，而不同波長的湍流對斷面的非定常效應(yīng)則用氣動導(dǎo)納函數(shù)來考慮。

氣動導(dǎo)納作為聯(lián)系隨機(jī)抖振力和脈動風(fēng)的傳遞函數(shù)，在橋梁隨機(jī)抖振分析中起著重要的作用。因此，準(zhǔn)確地識別橋梁斷面的氣動導(dǎo)納函數(shù)，是進(jìn)行橋梁抖振分析的關(guān)鍵。Sears[2]在勢流理論的基礎(chǔ)上通過計(jì)算豎向簡諧脈動風(fēng)作用下二維薄機(jī)翼斷面受到的氣動力，得到了氣動導(dǎo)納函數(shù)。Davenport[1]針對桁架梁橋，提出了順風(fēng)向脈動作用下桁架梁橋的阻力氣動導(dǎo)納計(jì)算公式。目前，對于鈍體的橋梁斷面，氣動導(dǎo)納的識別主要以風(fēng)洞試驗(yàn)為主。按照在風(fēng)洞中來流湍流的生成方式不同，氣動導(dǎo)納的識別可以分為被動格柵識別方法和主動控制格柵識別方法兩種。由于被動格柵的簡易性，大量研究者對此進(jìn)行了詳細(xì)的研究，如：Sarkar[3]、Larose和Tanaka[4]、Sankaran和Jancauskas[5]、Gu和Qin[6]等。相比于被動格柵，主動格柵需要特殊的裝置，因此目前對于主動格柵識別氣動導(dǎo)納的研究還較少[7-9]。與被動格柵相比，主動格柵利用運(yùn)動機(jī)構(gòu)生成指定的湍流度、湍流積分尺度和功率譜密度的風(fēng)場，有助于對氣動導(dǎo)納進(jìn)行更加深入的研究。

除了傳統(tǒng)的風(fēng)洞試驗(yàn)，近年來CFD數(shù)值模擬也被用來進(jìn)行氣動導(dǎo)納的研究。Bruno等[10]通過階躍函數(shù)，識別了一箱梁斷面在兩種湍流度下的氣動導(dǎo)納函數(shù)；Hejlesen等[11]利用二維的離散渦方法，計(jì)算了四種橋梁斷面的氣動導(dǎo)納函數(shù)；唐煜等[12]借助于商用計(jì)算流體軟件Fluent，計(jì)算了平板和箱梁斷面在簡諧來流下的氣動導(dǎo)納。相比于主動格柵，CFD數(shù)值模擬不需要專門的試驗(yàn)設(shè)備、易于控制來流的風(fēng)場特性、可重復(fù)性好等優(yōu)點(diǎn)。在氣動導(dǎo)納的識別上，CFD數(shù)值模擬具有巨大的潛力。目前CFD數(shù)值模擬雖然已經(jīng)廣泛用于橋梁斷面的靜氣動力[13]，渦激力[14]，氣動自激力[15]等計(jì)算，但是利用CFD識別氣動導(dǎo)納的研究尚處于起步階段。

本文基于二維不可壓URANS方法，分別計(jì)算了平板、箱梁、雙邊肋斷面在簡諧脈動來流和寬頻湍流中的非定常氣動力，并識別出各自的氣動導(dǎo)納函數(shù)。在此基礎(chǔ)上，通過考察斷面在三種脈動幅值來流下的氣動導(dǎo)納，研究了不同來流特性對各類斷面氣動導(dǎo)納的影響。

1 氣動導(dǎo)納識別方法

基于準(zhǔn)定常假定的抖振力理論把橋梁斷面上的隨機(jī)氣動力看作是靜氣動力、湍流脈動分量的函數(shù)，它忽略了不同頻率脈動風(fēng)對氣動力的影響。這在湍流積分尺度遠(yuǎn)遠(yuǎn)大于橋梁斷面的特征長度時(shí)是成立的。當(dāng)來流的湍流積分尺度與橋梁斷面的特征長度處于同一量級時(shí)，必須引入氣動導(dǎo)納函數(shù)來考慮氣動力的非定常效應(yīng)。在僅有豎向簡諧脈動來流作用下橋梁斷面的抖振力可以表示為[16]：

由式(1)和式(2)可以得到升力和扭矩的功率譜密度：

其中，SL和SM分別為升力和扭矩的功率譜密度；Sw為豎向脈動風(fēng)的功率譜密度；χLw2和χMw2為升力和扭矩關(guān)于豎向脈動風(fēng)的氣動導(dǎo)納。

據(jù)式(3)和式(4)可以得到χLw2和χMw2的表達(dá)式：

這里也給出Liepmann簡化的Sears函數(shù)以供參考，如下式所示：

其中k=fB/U，為無量綱頻率。

2 CFD數(shù)值模型

2.1 數(shù)值計(jì)算方法

本文數(shù)值計(jì)算采用三種不同流線程度的斷面，分別為：寬高比B/D=103∶1的平板斷面，流線型的箱梁斷面和鈍體的雙邊肋斷面。具體的斷面形式如圖1所示。

(a) 平板斷面

(b) 箱梁斷面

(c) 雙邊肋斷面

繞斷面的不可壓流動由Navier-Stokes方程控制，通過商用的流體計(jì)算軟件ANSYS Fluent求解。其中，壓強(qiáng)速度的耦合采用PISO方法，時(shí)間項(xiàng)的推進(jìn)采用二階的隱式格式，對流項(xiàng)的離散采用三階的MUSCL格式。湍流的模擬采用SSTk-ω湍流模型。在基于RANS方法的湍流模型中，SSTk-ω湍流模型由于對來流特性的不敏感性、對邊界層分離和再附的準(zhǔn)確預(yù)測已經(jīng)被廣泛應(yīng)用于工程實(shí)踐中[17-18]。需要注意的是，相比于直接數(shù)值模擬(DNS)和大渦模擬(LES)，兩方程的RANS湍流模型不能捕捉到高頻范圍的脈動，因此也不能模擬出高頻范圍內(nèi)的氣動力效應(yīng)。但是，現(xiàn)代大跨度橋梁關(guān)心的折算頻率一般位于0～1之間，在本文對應(yīng)于f=0～25.8 Hz，而RANS湍流模型能夠分辨的頻率范圍要遠(yuǎn)大于這個(gè)范圍，如圖6所示。

模型的計(jì)算域如圖2所示。對于簡諧脈動來流，計(jì)算域的左側(cè)采用速度入口邊界條件，水平向速度分量u不隨時(shí)間變化，豎向速度分量w為正弦脈動；計(jì)算域的上下邊界也采用速度入口，豎向速度除了隨時(shí)間變化以外也是空間的函數(shù)，以保持與左側(cè)的速度進(jìn)口相協(xié)調(diào)。對于寬頻湍流來流，計(jì)算域左側(cè)的速度入口，水平速度與豎向速度都采用寬頻的人工合成湍流；計(jì)算域的上下邊界采用對稱邊界條件。右側(cè)的出口都采用壓強(qiáng)出口條件。

圖2 計(jì)算域及邊界條件Fig.2 Computational domain and boundary conditions

此外，為了保證上下側(cè)邊界不對內(nèi)部的繞流場產(chǎn)生影響，上下側(cè)邊界之間的距離應(yīng)該取得充分大。通過對8B、12B、16B、20B四種工況的計(jì)算，結(jié)果表明12B的距離可以滿足本文的計(jì)算需求。

計(jì)算域的網(wǎng)格采用混合網(wǎng)格的形式。在靠近模型的范圍內(nèi)采用結(jié)構(gòu)化的四邊形網(wǎng)格，近壁面的第一層網(wǎng)格高度yw=1.0×10-4B，可以保證壁面處的y+小于1。在結(jié)構(gòu)化網(wǎng)格之外是三角形的非結(jié)構(gòu)化過渡網(wǎng)格，網(wǎng)格的最大值為2.0×10-3B。此外，對尾流區(qū)內(nèi)的網(wǎng)格進(jìn)行加密。三角形網(wǎng)格的外面區(qū)域采用結(jié)構(gòu)化網(wǎng)格進(jìn)行填充。如圖3為箱梁斷面壁面附近的網(wǎng)格。

圖3 箱梁斷面附近的網(wǎng)格Fig.3 Computational mesh around the box section

2.2 正弦脈動流的數(shù)值計(jì)算

在無障礙物的流場中，保證來流的主要特征(如：湍流度、積分尺度、功率譜密度等)可以維持到計(jì)算域的出口不變，即來流的自保持特性，是CFD計(jì)算準(zhǔn)確性的基本要求。在基于有限體積法的CFD數(shù)值求解器中，誤差主要來源于微分方程的離散過程。一般認(rèn)為，計(jì)算域網(wǎng)格的分辨率和時(shí)間步長是來流自保持特性的關(guān)鍵影響因素。

對于豎向簡諧脈動來流，附加在主流上的豎向脈動速度對應(yīng)的波長為：

其中f為豎向脈動風(fēng)的頻率。

為了使簡諧脈動波在計(jì)算域內(nèi)有效的傳播，應(yīng)保證在一個(gè)波長內(nèi)有足夠多數(shù)量的網(wǎng)格，即：

λ=NΔx(9)

其中,Δx為網(wǎng)格的尺寸，N為一個(gè)波長內(nèi)的網(wǎng)格數(shù)量。

由折算頻率k的表達(dá)式及式(8)、式(9)，可得到用k表達(dá)的網(wǎng)格尺寸Δx的表達(dá)式：

根據(jù)唐煜等[12]的研究，當(dāng)N≥80時(shí)可以滿足流場中相鄰兩波峰間幅值的對數(shù)衰減率δ=In(Ai+1/Ai)≤ 0.003。其中Ai+1、Ai分別為相鄰兩波峰間的幅值。

另外，關(guān)于時(shí)間步長的取值，本文選擇在一個(gè)簡諧脈動周期內(nèi)保證至少有200個(gè)時(shí)間步，這比唐煜等[12]建議的無量綱時(shí)間步長Δt=dtU/Δx≤1.2要小。

圖4所示為k=0.03和k=1時(shí)，無障礙流場中坐標(biāo)(x,y)=(8B，0)處，簡諧脈動速度時(shí)程與目標(biāo)值的比較。從圖中可以看到，在低折算頻率與高折算頻率下模擬值與目標(biāo)值吻合良好，說明采用以上的數(shù)值設(shè)定可保證來流的自保持特性。

2.3 寬頻湍流的數(shù)值計(jì)算

入口邊界處的寬頻湍流，采用Davidson[19]等提出的人工譜合成方法。利用該方法可以生成指定湍流度、積分尺度和功率譜密度函數(shù)的隨機(jī)脈動風(fēng)速時(shí)程。本文以Karman譜作為目標(biāo)風(fēng)譜，湍流度與簡諧脈動來流的湍流度相接近。

為了降低湍流度、湍流積分尺度等湍流特性沿計(jì)算域的衰減，模型斷面到入口處的網(wǎng)格要足夠密。通常我們關(guān)心的折算頻率k≈0～1之間，參照正弦脈動來流的網(wǎng)格大小，并考慮計(jì)算效率的因素，模型斷面到入口處的網(wǎng)格取Δx=B/(50×1)。數(shù)值離散格式和時(shí)間步的取值同樣參照正弦簡諧脈動來流的設(shè)置，即，對流項(xiàng)的離散采用三階的MUSCL格式，無量綱時(shí)間步Δt取1。

圖5所示為湍流度和湍流積分尺度沿計(jì)算域中心線y=0的變化。從圖中可以看出，湍流度和湍流積分尺度從入口進(jìn)入計(jì)算域后都呈現(xiàn)出一定程度的衰減，但是隨著入口距離的增大，湍流的演化逐漸平穩(wěn)，在距離入口1.5 m后湍流度和湍流積分尺度都趨于穩(wěn)定。

(a) k=0.03

(b) k=1

圖5 湍流度和積分尺度沿x軸的變化 (y=0)Fig.5 Turbulent intensity and turbulent length scale along the x-axis (y=0)

圖6所示為無障礙流場中坐標(biāo)(x,y)=(8B，0)處豎向脈動風(fēng)速的功率譜密度與目標(biāo)值的比較。由圖6可見，在我們關(guān)心的頻率范圍內(nèi)(f<30 Hz)，模擬的隨機(jī)豎向脈動風(fēng)速功率譜密度與目標(biāo)值吻合良好。

圖6 坐標(biāo)(x，y)=(8B，0)處豎向脈動速度功率譜密度Fig.6 Power spectrum of the vertical gust at (x，y)=(8B，0)

3 數(shù)值計(jì)算結(jié)果

3.1 單頻與寬頻脈動風(fēng)下的氣動導(dǎo)納

首先考慮兩種來流條件，即單頻簡諧脈動風(fēng)與寬頻湍流脈動風(fēng)。兩種情況的平均風(fēng)速均設(shè)為8 m/s。對于單頻簡諧脈動來流，其豎向簡諧速度的幅值為0.226，對應(yīng)的湍流度為2%，而順風(fēng)向的湍流度很小可忽略不計(jì)；對于寬頻湍流來流，豎向脈動速度的湍流度為2%，其順風(fēng)向的湍流度與豎向的湍流度基本一致。

在氣動導(dǎo)納識別前，首先在均勻來流的情況下計(jì)算出三種斷面的靜力三分力系數(shù)及其對迎角的導(dǎo)數(shù)，然后在無障礙物的流場中計(jì)算得到橋梁斷面中心點(diǎn)處的脈動風(fēng)速時(shí)程。最后在得到橋梁斷面受到的三分力時(shí)程后，根據(jù)本文第一節(jié)所述方法可以識別出氣動導(dǎo)納函數(shù)。

圖7分別為平板、流線型的箱梁和鈍體雙邊肋斷面在兩種不同來流下的氣動導(dǎo)納函數(shù)。

理想平板斷面的氣動導(dǎo)納，存在理論解，即Sears函數(shù)。所以，平板斷面可以用來驗(yàn)證本文數(shù)值方法的準(zhǔn)確性。由圖7可知，在豎向簡諧脈動來流下，平板斷面的氣動導(dǎo)納函數(shù)與Sears函數(shù)吻合良好，證明了本文數(shù)值識別氣動導(dǎo)納的可行性。相比之下，湍流場下的氣動導(dǎo)納函數(shù)與Sears函數(shù)具有相同的趨勢，但卻表現(xiàn)出很強(qiáng)的隨機(jī)跳躍性，與Sears函數(shù)之間有實(shí)質(zhì)性的差異。

對于流線型的箱梁斷面，數(shù)值識別結(jié)果與Sears函數(shù)相比，在低折算頻率范圍內(nèi)與Sears函數(shù)相近，而在高折算頻率范圍內(nèi)比Sears函數(shù)低。雙邊肋斷面的氣動導(dǎo)納函數(shù)與Sears函數(shù)之間差別明顯，這是符合鈍體斷面特征的合理現(xiàn)象。

(a) 平板斷面

(b) 箱梁斷面

(c) 雙邊肋斷面

3.2 單頻脈動風(fēng)幅值對氣動導(dǎo)納的影響

由于橋梁斷面的鈍體性質(zhì)，氣動導(dǎo)納一般不存在與理想平板斷面類似的理論解。在現(xiàn)代橋梁抖振分析中，某次試驗(yàn)識別的氣動導(dǎo)納經(jīng)常被應(yīng)用于任意風(fēng)場環(huán)境下的全橋抖振響應(yīng)分析。這種做法暗含了氣動導(dǎo)納不隨風(fēng)場特性改變，僅是無量綱折算頻率的單值函數(shù)。然而，對于鈍體斷面，一些研究者的試驗(yàn)表明，氣動導(dǎo)納函數(shù)除了與斷面的形式有關(guān)外，還與來流的湍流度、積分尺度等有關(guān)[7, 20]。

ANSYS Fluent，通過用戶自定義函數(shù)(UDF)，可以輕易地改變計(jì)算域入口邊界的來流特性。通過改變簡諧脈動來流的幅值，下面討論簡諧脈動風(fēng)幅值對不同流線型程度橋梁斷面氣動導(dǎo)納的影響。

三種簡諧脈動來流的特性如表1所示，來流的平均速度都為8 m/s。

表1 CFD計(jì)算用簡諧來流特性Table 1 Wind properties used in simulations

圖8、圖9分別為平板、流線型箱梁、鈍體雙邊肋斷面在不同來流下的升力氣動導(dǎo)納函數(shù)和扭矩氣動導(dǎo)納函數(shù)。

由圖8、圖9可以看到，平板與箱梁斷面的升力氣動導(dǎo)納、扭矩氣動導(dǎo)納除了在較低折算頻率處離散性稍大外，在其它折算頻率范圍內(nèi)基本一致，體現(xiàn)了流線型斷面對風(fēng)場不具有依賴性或者具有弱依賴性。相比較之下，鈍體雙邊肋斷面在三種來流下的升力氣

(a) 平板斷面

(b) 箱梁斷面

(c) 雙邊肋斷面

動導(dǎo)納和扭矩氣動導(dǎo)納差別較大，表現(xiàn)出較大的離散性，尤其是在高折算頻率范圍內(nèi)，體現(xiàn)了鈍體斷面對風(fēng)場的依賴性。表明鈍體的氣動導(dǎo)納并不能定性為斷面的函數(shù)，而是由斷面與風(fēng)場共同決定。這主要是由于鈍體斷面周圍會產(chǎn)生復(fù)雜的流動分離與再附現(xiàn)象，會對來流簡諧脈動引起的氣動力產(chǎn)生干擾。而又由于斷面周圍流動分離與再附現(xiàn)象與來流特性密切相關(guān)[19]，因此使得最終識別的氣動導(dǎo)納依賴與來流特性。

為了進(jìn)一步說明不同斷面形式氣動力的差異，作為例子，圖10、圖11、圖12分別給出了平板、箱梁和雙邊肋斷面當(dāng)k= 0.8時(shí)在均勻來流、wind field 1、wind field 2、wind field 3來流下的升力系數(shù)時(shí)程曲線。圖13、圖14、圖15分別為對應(yīng)的升力幅譜。其中平板斷面因升力為0，所以省略平板斷面在均勻來流下的升力系數(shù)和升力譜曲線。

(a) 平板斷面

(b) 箱梁斷面

(c) 雙邊肋斷面

(a) wind field 1

(b) wind field 2

(c) wind field 3

平板斷面因其流線性質(zhì)，升力時(shí)程只包含來流簡諧脈動引起的氣動力，對應(yīng)的卓越頻率為25.82Hz。箱梁斷面的尾部會產(chǎn)生周期性的旋渦脫落，從圖14可知，在均勻來流下因旋渦脫落引起的升力脈動的頻率為66.4 Hz。雙邊肋斷面的幾何外形復(fù)雜，在它周圍的流動現(xiàn)象也更為復(fù)雜。從圖15可知，雙邊肋斷面的升力時(shí)程信號包含有四個(gè)顯著的峰值，對應(yīng)的頻率分別為：f1=13.4Hz、f2=28.1Hz、f3=41.5Hz、f4=56.2 Hz。從圖15也可以看到，頻率f2因與來流的脈動頻率(25.82 Hz)接近，因此按頻率f2脫落的旋渦誘導(dǎo)的氣動力會對來流簡諧脈動引起的氣動力產(chǎn)生干擾，從而使該頻率處對應(yīng)的氣動力不再滿足疊加原理。

(a) uniform flow

(b) wind field 1

(c) wind field 2

(d) wind field 3

表2為三種斷面在來流下升力幅譜及升力幅譜與來流脈動幅值的比值。從表2可以看出平板以及箱梁斷面升力時(shí)程的幅譜與來流脈動幅值的比值在三種來流中非常接近。這說明流線型斷面的氣動導(dǎo)納與來流的脈動幅值無關(guān)，僅是斷面及無量綱折算頻率的函數(shù)。相反，具有顯著分離流的鈍體雙邊肋斷面，斷面周圍不同尺度的渦誘導(dǎo)的氣動力會對來流簡諧脈動產(chǎn)生的氣動力產(chǎn)生干擾效應(yīng)。從表2可以看到雙邊肋斷面在不同來流下升力時(shí)程的幅譜與來流脈動幅值的比值差別非常大，分別為0.4217、0.3120和0.146。雙邊肋斷面對來流風(fēng)場的依耐性，表明鈍體斷面的氣動導(dǎo)納并不能定性為斷面的函數(shù)，而是由斷面與來流風(fēng)場特性共同決定的。

(a) uniform flow

(b) wind field 1

(c) wind field 2

(d) wind field 3

此外，值得一提的是，從圖12可以看出，對于雙邊肋斷面不同性質(zhì)的脈動風(fēng)對平均氣動力特性有著重要的影響，表現(xiàn)為隨著脈動風(fēng)幅值的增加，平均升力顯著降低。這一現(xiàn)象與Bearman對垂直與流向放置的平板的平均阻力的研究結(jié)果類似[22]。迄今為止，很多研究者對來流湍流度對鈍體斷面平均氣動力

(a) wind field 1

(b) wind field 2

(c) wind field 3

(a) uniform flow

(b) wind field 1

(c) wind field 2

(d) wind field 3

的影響進(jìn)行了較為深入的研究[21-23]。其中，Bearman認(rèn)為，這種現(xiàn)象本質(zhì)上來源于來流湍流對邊界層分離后的剪切層及其剪切層再附的影響。

(a) uniform flow

(b) wind field 1

(c) wind field 2

(d) wind field 3

圖15 雙邊肋斷面升力系數(shù)幅譜Fig.15 Amplitudespectra of lift coefficients for the double-girder section

4 結(jié) 論

基于二維非定常RANS模型，計(jì)算了三種不同流線程度的斷面在簡諧脈動來流和寬頻湍流下的氣動導(dǎo)納函數(shù)，在此基礎(chǔ)上研究了不同來流特性對給定斷面的氣動導(dǎo)納的影響。根據(jù)本文數(shù)值計(jì)算結(jié)果可得出以下結(jié)論：

1) 在豎向簡諧脈動來流情況下，平板斷面的氣動導(dǎo)納與Sears函數(shù)吻合良好，表明了在兩種不同形式的來流下，本文數(shù)值模型及數(shù)值方法的可行性。相比之下，湍流作用下，計(jì)算所得的氣動導(dǎo)納函數(shù)具有很強(qiáng)的隨機(jī)跳躍性，與薄平板斷面的理論解之間有實(shí)質(zhì)性的差異。

2) 本文計(jì)算的平板與箱形斷面的氣動導(dǎo)納函數(shù)均接近流線體斷面的特性，即對風(fēng)場不具有依賴性或弱依賴性，不同風(fēng)場形成的差異在數(shù)值模型的不確定性范圍內(nèi)；相比之下雙邊肋鈍體斷面的氣動導(dǎo)納表現(xiàn)出明顯的風(fēng)場依賴性，表明鈍體的氣動導(dǎo)納并不能定性為斷面的函數(shù)，而是由斷面與風(fēng)場共同決定。

3) 雙邊肋斷面的升力系數(shù)計(jì)算結(jié)果表明，不同性質(zhì)的脈動風(fēng)對平均氣動力特性有著重要的影響，在本文中表現(xiàn)為隨著脈動風(fēng)幅值的增加，平均升力顯著降低。這一現(xiàn)象與Bearman對垂直于流向放置的平板的平均阻力的研究結(jié)果類似[22]，本質(zhì)上來源于湍流特性對邊界層的發(fā)展以及平均風(fēng)場的影響。