三類氣動導納數值識別方法的適應性研究

張偉峰, 張志田, 張顯雄, 陳政清

(1. 湖南大學土 風工程與橋梁工程湖南省重點實驗室, 長沙 410082;2. 華北水利水電大學 土木與交通學院, 鄭州 450045; 3. 海南大學 土木與建筑工程學院, 海口 570228)

0 引 言

受來流湍流的影響，處于大氣邊界層中的橋梁都會受到抖振力的作用。抖振力與來流的脈動風特性、靜力三分力系數、氣動導納函數等有關。氣動導納作為聯系脈動風與抖振力的傳遞函數，其準確性對于橋梁抖振具有重要的意義。

針對于桁架橋，Davenport[1]用速度互相關來計算阻力氣動導納，而升力氣動導納采用基于勢流理論推導得到的Sears函數。對于流線型的斷面，在沒有試驗結果的前提下，Sears函數經常被采用。隨著試驗技術的進步，大量的研究者借助于風洞試驗進行氣動導納的研究。但是風洞試驗方法通常存在湍流積分尺度明顯偏小，低頻成份顯著不足[2-3]，難以得到任意目標風場，可重復性差等問題。

在風場效應可以線性疊加的情況下，氣動導納可以利用階躍響應函數得到。在機翼理論里，Wagner[4]和Küssner[5]通過考察機翼姿態的階躍變化、以及穿過半無限陣風場的氣動力行為，分別得到了Wagner函數和Küssner函數，用以描述機翼所受到的氣動力隨時間的演化。通過Fourier變化，Garrick[6]證明了Wagner函數和頻域里描述氣動自激力的Theodorsen函數互成Fourier變換對。同樣Küssner函數和Sears函數也成Fourier變換對。因此，如果能得到Wagner函數和Küssner函數，通過Fourier變換就可以得到Theodorsen函數和Sears函數。Caracoglia和Jones[7]最早設計了一套裝置識別得到不同斷面的Wagner函數。而Küssner函數的試驗識別，由于無法得到理想的階躍風場，至今還沒有在風洞試驗里實現過。

至今為止，CFD已經廣泛應用于橋梁抗風實踐中。但與CFD在靜風力系數[8-9]、渦激振動[10-11]、顫振導數識別[9,12]等方面的廣泛應用相比較，氣動導納的數值研究卻鮮有報道。Hejlesen等[13]利用離散渦方法，識別了四種橋梁斷面的氣動導納函數；唐煜等[14]基于雷諾平均方法，在簡諧來流中識別了平板斷面和箱梁斷面的氣動導納；Bruno等[15]利用階躍函數的方法計算了機翼斷面和箱梁斷面的氣動力隨時間演化的曲線并識別得出各自的氣動導納函數；張偉峰等[16-17]在簡諧來流與湍流中利用CFD研究了橋梁斷面氣動導納與湍流參數的關系，對數值計算結果與風洞試驗結果進行了比較。

風洞試驗中存在的問題，可在合理的CFD模擬中克服，從而有望提高氣動導納的識別精度。但氣動導納的CFD識別研究尚處于起步階段，不同識別方法的計算策略、適應性、計算精度等問題尚沒有研究者進行研究。本文基于前期研究成果，對三種氣動導納數值識別方法的計算策略、適應性等問題進行研究。

本文采用二維SSTk-ω湍流模型，首先研究了簡諧脈動流、湍流和豎向階躍流在計算域內的傳播特性和數值計算方法，然后在這三種來流下分別計算了平板斷面和箱梁斷面的非定常氣動力，并識別得出各自的氣動導納函數和階躍響應函數，最后分析比較了不同方法的適應性以及計算效率。本文三種氣動導納數值識別方法的研究，對提高氣動導納數值識別精度以及工程應用具有很好的參考價值。

1 氣動導納識別方法概述

1.1 簡諧脈動流識別方法

處于二維脈動風場中的橋梁斷面，受到的抖振力的頻域表達式為[18]：

當僅考慮豎向脈動風作用時，由式(1)可以得到升力的功率譜密度：

其中：Sw為豎向脈動風的功率譜密度，|χLw|2為氣動導納模的平方。

根據式(2)可以得到|χLw|2的表達式：

當考慮機翼斷面在豎向脈動風下的作用時，|χLw|2即Sears氣動導納函數，可以近似表示為[19]：

通過在來流中給定單一頻率的豎向簡諧脈動，在得到了斷面的氣動力時程后，就可以根據公式(3)識別出氣動導納。該法識別結果具有較好的平滑性，但一次只可以識別出一個頻率點的氣動導納。

1.2 湍流識別方法

對于湍流，由于水平脈動風的影響，無法直接從式(1)中得到氣動導納。這里我們采用互譜識別方法[20]。升力的自功率譜密度為：

升力關于水平向及豎向脈動風的互功率譜為：

由此可以得到升力的氣動導納：

其中：SLu、SLw為升力關于水平和豎向脈動風的互功率譜；Suw、Swu為水平脈動風和豎向脈動風的互譜。

相比于簡諧脈動流方法，湍流識別方法可以一次性得到所有頻率范圍內的氣動導納函數。但是，該法在某一頻率點的識別精度嚴重依賴于該頻率的信號強度，識別結果的平滑性很差，應用前需要進行處理。本文僅討論豎向脈動風關于升力的氣動導納，其它氣動導納可以采用相同的方法進行研究。

1.3 Küssner識別方法

以水平速度U飛行的機翼，穿過幅值為w0的豎向階躍陣風時，Küssner得到了其氣動力隨時間演變的公式[5]：

其中：s=2tU/B為無量綱時間，ψ(s)為Küssner函數，可以近似表示為[21]：

ψ(s)=1-0.5e-0.13s-0.5e-s(11)

在線性疊加原理成立的前提下，利用杜哈梅積分，任意分布的豎向脈動風w(s)引起的氣動升力可以表示為：

對上式進行變量替換，利用分部積分可以得到:

求出上式的功率譜密度，并與式(2)相比較，可以得到階躍函數ψ與氣動導納的關系：

2 自由來流在計算域中傳播的數值研究

數值計算在商業軟件ANSYS FLUENT 15.0中進行。湍流模型選用基于RANS的二維SSTk-ω湍流模型。

數值模擬斷面非定常氣動力的第一步是對計算域內來流的演變特性進行研究，確定合適的離散格式、網格尺寸、時間步等參數以保證來流的主要特征在計算域內準確傳播。對于湍流來說來流的主要特征有湍流度、積分尺度、功率譜密度等。對于簡諧脈動來流來說，需要保證來流的幅值在計算域內維持不變。而對于豎向階躍流來說，需要保證來流的階躍特性在計算域內維持不變。下面分別討論離散格式、網格尺寸、時間步等參數對來流在計算域中傳播的影響。

2.1 正弦脈動流在計算域中的傳播

對于簡諧脈動來流，計算域的左側入口采用速度入口邊界條件，u不隨時間變化，w為正弦脈動；計算域的上下邊界同樣采用速度入口，w既是時間的函數也是空間的函數；出口采用壓強出口邊界條件[16](如圖1所示)。

為保障簡諧波在計算域內有效傳播，應在一個波長內有足夠多數量的網格，即：

圖1 計算域及邊界條件Fig.1 Computational domain and boundary conditions

λ=NΔx(15)

其中:Δx為網格尺寸大小，N為一個波長內的網格數量。

由折算頻率k=fB/U及式(15)，可得到網格尺寸Δx的表達式：

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

為了研究對流項離散格式對簡諧脈動波在計算域內傳播的影響，分別采用三種不同的對流項離散格式，擴散項統一采用二階中心差分。從圖2可以看出，一階迎風格式由于截斷誤差的首項包含有二階導數，會導致較大的數值耗散，使得簡諧脈動波的幅值變小。相比較而言，二階迎風格式由于截斷誤差的首項為三階導數，因此它引起的數值耗散就遠遠小于一階迎風格式。三階MUSCL(Monotone Upstream-Centrered Schemes for Conservation Laws)格式的數值耗散最小，簡諧脈動在整個計算域內的衰減基本可以忽略不計。所以，對流項的離散選擇MUSCL格式。

時間項的離散統一采用二階隱式格式，圖3為無量綱時間步Δt=dt·U/B對簡諧脈動幅值的影響。可以看出，隨著無量綱時間步的減小，脈動幅值的衰減率也迅速減小。當無量綱時間取1時，基本可以滿足計算要求。

確定了數值計算參數以后，圖4分別計算了無斷面計算域中，斷面中心位置處折算頻率k=0.03和k=1時速度時程與目標值的比較。可以看到，在低折算頻率和高折算頻率下模擬值與目標值均吻合良好，說明采用以上的數值設置可以保證來流的幅值特性。

圖3 無量綱時間步對脈動幅值的影響Fig.3 Influence of the time step on the amplitude of the fluctuating wind

(a) k=0.03

(b) k=1

2.2 湍流在計算域中的傳播

入口采用速度進口邊界條件。入口邊界處的脈動速度，采用Davidson[22]等人提出的人工譜合成方法。選取Karman譜作為目標風譜，豎向脈動風的湍流度與簡諧脈動來流的湍流度接近。上下邊界采用對稱邊界條件，出口采用壓強出口邊界條件。

為了減小湍流的主要特征在計算域內的衰減，需要適當加密模型斷面到入口處的網格。參考公式(15)，并綜合考慮計算效率的因素，模型到入口處網格的尺寸取Δx=B/(50×1)。對流項的離散采用三階MUSCL格式，Δt取1。

圖5為湍流度和積分尺度沿x軸的變化。可以看出，湍流度和積分尺度僅呈現輕微的衰減。

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

圖6為無障礙流場中(x，y)=(8B，0)處豎向脈動風速功率譜密度與目標值的比較。可見，除了較低頻和高頻以外，模擬值與目標值吻合良好。

圖6 (x，y)=(8B，0)處豎向脈動速度功率譜密度Fig.6 Power spectrum of the vertical gust at (x，y)=(8B，0) in the absence of body sections

2.3 豎向階躍流在計算域中的傳播

數值模擬豎向階躍流會遇到以下幾個問題：

(1) 由于控制流動的微分方程需要在空間和時間上進行離散，因此嚴格的豎向階躍流在CFD數值模擬時是不可能實現的。

(2) 由于豎向階躍流在空間和時間上的急劇變化，因此在求解流動方程時會導致非物理的數值振蕩。

(3) 由于分子黏性、湍動能黏度、數值耗散等影響會抹平來流的階躍特性。此外，由于流動在時間和空間上的演化，來流的階躍特性也會受到影響。

嚴格的階躍變化是不可能實現的，因此采用一個光滑變化的“階躍函數”H(s)來作為近似，如圖7所示。階躍函數H(s)的幅值為Hmax，假設當H(s0.99)=0.99Hmax時階躍函數達到穩定，其中s0.99為H(s)從零變化到0.99Hmax所用的時間。為了保證H(s)的階躍特性又考慮到數值穩定性，我們取s0.99<0.1。與Küssner函數的演化特性相比，這是一個相當小的值。s0.99與網格尺寸、時間步、空間和時間離散格式等有關。

圖7 階躍函數和H(s)Fig.7 Heaviside function and H(s)

在豎向階躍流的數值計算中，入口采用速度進口邊界且使豎向速度做階躍變化，出口采用壓強出口邊界條件，上、下邊界的流動條件一致所以采用周期性邊界條件。

網格尺寸對豎向階躍流的階躍特性具有顯著的影響。從圖9可以看出，隨著網格間距x的減小，階躍來流的斜率也隨之變大。當x=B/100時，s0.99<0.1，可以近似保證來流的階躍特性。

圖8 對流項離散格式對階躍特性的影響Fig.8 Tests of the interpolation schemes for the convection terms

圖9 網格大小對階躍特性的影響Fig.9 Effect of the grid spacing on the characteristics of the vertical indicial wind

圖10所示為時間步對豎向階躍來流階躍特性的影響。可以看出，當時間步較大時，由于計算域內速度的突變，會產生較大的越界現象。當無量綱時間t=0.001時，可以近似保證來流的階躍特性。

圖10 時間步對階躍特性的影響Fig.10 Effect of the time step on the characteristics of the vertical indicial wind

3 數值計算及結果

確定了不同來流條件的數值計算參數后，對接近于理想流線型的平板斷面和箱梁斷面的非定常氣動力進行CFD模擬。數值模型的尺寸如圖11所示。

來流的水平速度U=8 m/s，對應的雷諾數Re=UB/ν≈ 1.6×105。對于簡諧脈動來流，豎向速度幅值取0.226 2 m，對應2%的湍流度。對于豎向階躍流，取Hmax/U=2%。所有斷面的網格劃分采用混合網格的形式。簡諧脈動來流下箱梁斷面的網格數為10.2萬～156萬，湍流的網格數為75萬，豎向階躍流的網格數為64萬。圖12所示為簡諧脈動來流下箱梁斷面壁面附近的網格。

(a) 平板斷面

(b) 箱梁斷面

圖12 箱梁斷面壁面附近的網格Fig.12 Computational grid close to the box girder model

計算域如圖1所示，圖中給出了簡諧脈動來流下的邊界條件。模型的前緣距入口的距離為8B，模型的后緣距出口的距離為25B。上下側邊界之間的距離應該取得充分大，以避免邊界對內部的流場產生影響，試算表明16B的距離可以滿足計算需求。

需要注意的是，根據2.3節分析的豎向階躍來流的網格要求，為了保證來流的階躍特性，會導致模型前緣的計算域中網格過密，給計算造成一定的負擔。此外，階躍來流傳播到斷面前緣也需要一定的時間。為此，對流場使用非均勻的初始化。步驟如下：

(1) 在均勻來流的情況下(U=U∞，w=0)，進行定常計算。

(2) 在x方向選擇一個坐標點x0，待流場穩定后，將x

3.1 平板斷面

理想平板斷面的氣動導納存在理論解，即Sears函數。穿越半無限空間的豎向脈動風的理想平板，其升力隨時間的變化滿足Küssner函數，因此平板斷面首先被用來驗證本文數值方法的準確性。

圖13 豎向階躍流計算示意圖Fig.13 Diagram of the set up of the vertical indicial wind

在簡諧脈動來流下，由于平板的氣動力呈簡諧變化，因此數值計算模擬了5 s的氣動力時程。湍流引起的氣動力具有隨機性，需要很長的采樣時間，本文數值模擬了42 s的氣動力時程。Küssner函數在無量綱時間s≈30時趨于定常，本文數值模擬了25 s的氣動力時程，對應于無量綱時間s=50。本文的計算在曙光W580I工作站上進行(16物理核，32G內存，2T硬盤)，完成上述計算需要的CPU核時和內存如表1所示。可以看到，湍流法計算花費的時間最長，內存占用也是最多的。簡諧脈動方法雖然一次計算時間最短，但是需要對所有感興趣的頻率點進行掃頻，所以總的計算花費反而是最大的。相比其它兩種方法，Küssner方法的計算效率最高。

表1 三種氣動導納識別方法的計算效率比較Table 1 Computing costs of the three numerical methods

圖14所示為平板斷面在三種來流下升力系數隨時間的變化圖。從圖中可以看出平板斷面由于沒有渦脫，因此升力系數時程僅包含來流脈動的頻率。而在寬頻來流下，升力系數時程呈現出隨機特性。在豎向階躍來流下，升力系數先急劇增長，再經歷了一段緩慢增長后，最終趨于穩定達到定常狀態。

(a) 湍流

(b) 簡諧脈動來流

(c) 豎向階躍來流

圖15為階躍響應函數及由其得到的升力氣動導納函數。從圖中可以看出識別得到的階躍響應在初始時刻和最終時刻吻合較好，在中間時刻較Sears函數要大。這一差距可能是由流體的粘流引起，而Sear函數是在無粘性的有勢流場中得到。

(a) 階躍響應

(b) 氣動導納

圖16給出了三種來流下識別得到的升力氣動導納函數。可以看到利用湍流法和簡諧脈動法識別的氣動導納吻合良好，而且與Sears函數較為一致。Küssner法識別的氣動導納整體上也與Sears函數較為吻合，但在低頻有一定的差距。這一差距的原因，主要是由于不滿足流動的有勢條件。此外，可能還由于Küssner法涉及到階躍響應函數的識別、求導與傅立葉變換等多個步驟，見公式(14)，在這個過程中存在著誤差的積累與傳遞。

圖16 三種不同方法平板斷面氣動導納函數的比較Fig.16 The comparison of the aerodynamic admittance between three different methods for plate section

3.2 箱梁斷面

箱梁斷面在橋梁工程中應用十分普遍。本文選取的箱梁斷面寬高比B/D=10.2。圖17分別為箱梁斷面在三種來流下的升力系數時程。在寬頻來流下，同平板斷面一樣，升力系數時程呈現出隨機特性。在簡諧脈動來流下，由于斷面尾部的渦脫效應，升力系數時程除了包含來流的脈動頻率以外，還包含有渦脫頻率，對應的Strouhal數為0.2。豎向階躍來流下，升力系數時程與平板斷面的情況有本質的區別，呈現出很大的周期性振蕩，這也是由斷面尾部的渦脫造成的。從階躍響應函數和氣動導納的物理意義來看，它們并不包含渦脫的影響，這部分氣動力由渦激力模型考慮。因此為了得到階躍響應函數，首先使用FFT光滑濾波對升力系數進行處理，濾波光滑后的升力系數如圖17(c)所示。圖18為根據濾波光滑后的升力系數得到的階躍響應函數。可以看到箱梁斷面的階躍響應與Küssner函數相差較大，呈現出一個峰值。這可能是因為當階躍來流到達斷面前緣時，在斷面的前緣處引起瞬時的較大的邊界層分離，隨著邊界層的再附及向下游傳播，并最終與后緣脫落的漩渦合并，造成斷面表面壓強的變化引起的[23]。對于鈍體斷面，這種現象是Küssner方法所固有的。顯然Küssner方法并不適用于這種形式的斷面。

(a) 湍流

(b) 簡諧脈動來流

(c) 豎向階躍來流

圖18 箱梁斷面的階躍響應Fig.18 Indicial lift response function for box girder section

圖19所示為不同來流下箱梁斷面的升力氣動導納函數。從圖中可以看出，正弦來流和湍流下識別的氣動導納函數在低頻范圍內吻合良好，并且與Sears函數較為接近。在高頻范圍，正弦來流下識別的氣動導納函數略小于Sears函數，湍流下識別的氣動導納函數卻略大于Sears函數。Zhang等人在文獻[24]中討論了氣動導納的風場依賴性，指出對于具有顯著分離流的鈍體斷面，非定常氣動力的疊加原理不再成立，氣動導納應表現出對來流特性的依賴性。本文的計算表明，扁平箱梁斷面的氣動導納與來流風場特性僅僅是弱相關的。豎向階躍來流下識別的氣動導納整體上與其它兩種來流的結果相差較大，從前面的分析可見這種方法僅能應用于嚴格的流線型斷面，即與風場特性完全無關的斷面。

圖19 三種不同方法箱梁斷面氣動導納函數的比較Fig.19 The comparison of the aerodynamic admittance between three different methods for box girder section

4 結 論

本文借助于CFD方法，研究了適用于簡諧脈動來流、湍流和豎向階躍來流三種風場的計算方法。在此基礎上計算了平板和扁平箱梁斷面受到的氣動力，識別了階躍響應和氣動導納函數，得出以下結論：

(1) 要準確模擬三種來流風場需要采取不同的數值計算策略。需要按標準保證足夠的網格精度，而且對流項離散格式應采用高階離散格式，如三階的MUSCL格式。湍流的最小網格尺寸應根據最高截止頻率按照簡諧脈動來流的條件取值。豎向階躍來流由于在計算域內的局部位置會發生急劇變化，因此對流項的離散需要采用有界的格式，以防止越界現象產生。另外較大的時間步也會引起越界現象。

(2) 三種方法識別得到的平板斷面氣動導納與Sears函數吻合，說明了本文三種方法識別氣動導納的可行性。

(3) 簡諧來流和湍流下識別的箱梁斷面氣動導納差別不大，體現出了扁平箱梁斷面氣動導納對來流風場的弱相關性。在豎向階躍來流下識別的氣動導納函數與其它兩種來流下識別的氣動導納函數有較大的差距，體現了Küssner方法應用于鈍體斷面的局限性。

(4) 三種方法相比較而言，湍流的計算效率最高，可以一次性識別出所有頻率范圍的氣動導納函數，但識別結果平滑性差隨機跳躍性大；簡諧來流方法進行多次掃頻計算因而計算消耗大，但該法具有求解穩定、結果平滑可靠的優點；Küssner方法具有計算時間短的優勢，而且可以識別出階躍響應函數直接用于時域分析，但該法存在誤差積累與傳遞的問題，且僅能適用于氣動導納與風場嚴格無關的斷面。