基于兩種湍流模型的橋梁顫振導數識別研究及比較＊

祝志文，夏 昌

（湖南大學風工程試驗研究中心，湖南長沙 410082）

風洞試驗是研究大跨度橋梁顫振穩定性的主要手段.風洞試驗簡單直觀，但風洞試驗涉及模型制作、復雜的儀器設備，并且其試驗周期長，費用高，紊流風場的有效模擬比較困難.隨著計算機技術的快速發展，采用數值計算的方法來識別顫振導數成為可能，Larsen的離散渦法［1］，Zhu等開發的有限體積法［2］，曹豐產等開發的有限元方法［3］，都取得了重要的成果.其中部分方法還是在均勻流場中開展的，并沒有實現對湍流的模擬，雖然計算工作量減小，但網格等可能并沒有達到要求.大型CFD商業軟件FLUENT和CFX等以其通用性和便利的后處理功能受到數值風洞研究者的肯定，同時提供了多種湍流模型，因此它們在結構風工程領域的研究中越來越受到重視.

標準k－ε湍流模型是由Launder和Spalding［4］于1972年首先提出的，該湍流模型屬于高Re數模型，在近壁區并不直接求解，而是通過壁面函數將近壁區的變量與湍流核心區變量聯系起來求解，因此存在固有近邊界模擬的缺陷，同時在此模型中假設湍流黏性系數μt是相同的，而在流線彎曲和應變率大的情況下，湍流是各向異性的，因此，此模型用于強旋流、彎曲壁面流動或彎曲流線流動時，會產生一定的失真［5］.Menter在1993年對Wilcox提出的標準兩方程k－ω湍流模型進行修正后提出分區剪應力輸送（Shear Stress Transport）模型［6］，即SST湍流模型.該模型實際上是從邊界層內部的標準k－ω模型逐漸轉變到邊界層外部的高Re數的標準k－ε模型，在近壁區，它不再引入如壁面函數這樣的經驗公式，而是采用標準k－ω模型直接求解，能有效地使模型流線分離；在湍流核心區采用標準k－ε模型計算，同時，對湍流黏性系數μt進行修正，這樣既提高了數值模擬的精度，又有效減小了計算工作量.

過去由于受到計算條件的限制，標準k－ε湍流模型得到廣泛應用，同時人們根據需要提出了多種相關的修正模型，有力地推動了數值風洞的發展.近幾年來，隨著計算速度的大幅提高，研究者們才開始重視SST模型的使用，基于此模型利用數值計算的方法，成功地計算了橋梁斷面的三分力系數，數值計算值與試驗值吻合較好，體現了SST湍流模型的計算優勢.采用數值方法識別顫振導數時，研究者［7］常用標準k－ε湍流模型，采用SST湍流模型的相關報道卻很少見.

采用數值模擬橋梁斷面強迫振動，需要采用動網格技術.研究者［7］提出采用斷面外包矩形剛性網格的方法，以此提高數值模擬精度.但在矩形直角周圍的動網格質量卻不理想（疏密不均），由此可能引起在氣動力時程曲線中，出現大量“毛刺”.本文采用外包橢圓形剛性網格的形式，既能保證剛性網格區的質量，又能兼顧動網格的質量.

本文CFD計算基于大型商業軟件FLUENT平臺開展，湍流模型分別采用FLUENT提供的基于雷諾平均（RANS）的標準k－ε模型和SST模型，分析兩種湍流模型的使用條件，嘗試采用SST湍流模型識別橋梁斷面顫振導數，同時比較兩湍流模型在識別顫振導數上的計算精度.

1 控制方程及湍流模型分析

基本控制方程為連續方程與雷諾平均N－S方程（RANS）：

式中：ρ，μ分別表示流體的密度和動力粘度；ui，uj代表某個方向上的平均流速，u'i為速度分量的脈動量，對于二維問題，i與j的取值范圍為1，2；－項定義為Reynold應力［8］，這屬于新的未知量，為了使方程組可以封閉，故引入湍流模型（雷諾平均N－S方程），以便求解.

標準k－ε模型是高Re數的湍流模型，它只能求解湍流核心區的流動，直接來求解橋梁斷面近壁區的流動是不正確的.因此，使用此湍流模型時需要引入壁面函數，對湍流核心區的流動使用k－ε模型求解，對近壁區不進行求解，而是直接采用一組經驗公式，將近壁區的變量與湍流核心區的求解變量聯系起來求解.采用標準k－ε模型時，只要把第一層網格節點布置在湍流充分發展的區域就可以了，并不需要在橋梁斷面近壁區加密網格，因此減少了網格數量，節省了數值計算的時間.當計算機計算速度有限時，此模型顯示出明顯的優勢.為了合理地使用壁面函數，需要對斷面的Yplus值進行控制，一般要求滿足Yplus≥11.63［5］.

SST湍流模型是在k－ω模型基礎上發展而來的，且融合了k－ε模型的優點.k－ω湍流模型的優勢是在低Re數下的近壁計算，k－ε模型適合湍流核心區的計算.SST模型克服了標準k－ω湍流模型對自由流參數變化比較敏感的缺點，在近壁區采用k－ω湍流模型，在遠離壁面的流場中采用k－ε湍流模型.這樣充分利用了k－ε湍流模型對逆壓梯度流動具有較高的模擬精度和k－ε湍流模型對湍流自由流參數不敏感的優點.采用SST模型時，需要對近壁區進行數值計算，而不是采用類似的壁面函數，因此需要加密近壁區的網格，同時要合理控制第一層網格高度，大體上使得Yplus≤6［8］.

2 氣動自激力描述

氣動自激力的Scanlan表達式為：

式中：h（t），α（t）分別為豎彎與扭轉位移；Κ＝ωB／U為無量綱頻率；?（t），（t）分別為豎彎速度與扭轉速度；和（i＝1，2，3，4）即為橋梁斷面顫振導數，是橋梁斷面外形和折算風速Ur的函數（Ur＝U／fB）.

本文通過數值模擬的方法分別提取每個折算風速下的氣動自激力，然后采用最小二乘法識別相應折算風速下的8個顫振導數.

3 算 例

3.1 研究對象及計算域網格劃分

本文選取丹麥大帶東橋為研究對象，該橋為主跨1 624m的三跨連續鋼箱梁懸索橋，其主梁是上下游側帶風嘴的扁平箱梁斷面（如圖1所示）.箱梁截面全寬B＝31m，高寬比為7.05∶1.CFD計算的模型采用與風洞試驗相同的幾何縮尺比1∶80，數值模擬時不考慮欄桿和防撞墻等附屬物.計算域采用二維圓形區域，左側半圓弧為來流進口，到模型中心的距離為30B，右側半圓弧為來流出口，到模型中心的距離為30B.

圖1 大帶東橋主梁斷面（單位：m）Fig.1 Girder section of the Great Belt Bridge（unit：m）

由于CFD模擬時剛性橋梁斷面在每一時間步上運動，因而在每一時間步上需要重新對計算域網格進行劃分.為確保橋梁最大運動位移處有較好的網格質量，不至出現畸變網格甚至負體積網格，本文將計算域進行分區劃分，分成3個大小近似合理的區，如圖2所示，并采用不同的網格進行剖分.圍繞橋梁的稱為剛性網格區.橋梁斷面運動時，該區域網格與橋梁斷面剛性固定，并在每一時間步上與橋梁斷面同步運動.該區域外邊界為橢圓，通過對該橢圓適當分段，便于對該域進行四邊形結構網格劃分，以便能獲得較好的正交網格.計算域絕大部分區域稱靜止網格區，該區域外邊界是計算域外邊界，內邊界為離開剛性網格區外橢圓一定距離，且包圍剛性網格區的圓形.靜止網格區采用四邊形單元剖分，從內到外采用合適的網格放大比例.靜止網格區和剛性網格區在整個CFD模擬過程中一直使用計算開始時的網格系統，不再重新劃分網格.在靜止網格區和剛性網格區之間為動網格區，動網格區采用三角形單元剖分.在每一時間步上，該區域根據橋梁斷面的運動位置并由設定的網格系統質量要求重新進行網格劃分.緊靠橋梁斷面的區域流場變量變化劇烈，特別是斷面迎風側和斷面法向，因此網格劃分需要能適應流場變量的變化程度，并沿各個方向采用適度的網格放大率，實現與動網格區域網格的平順過渡，如圖3所示.

圖2 計算域分區Fig.2 Computational domain decomposition

值得注意的是，為了分別滿足本文所采用兩湍流模型的使用要求，需要劃分兩套網格，如圖3所示.這兩套網格計算域劃分形式完全一樣，只有剛性網格區的第一層網格高度和網格數量不一樣.因為需要通過控制第一層網格的高度來調節斷面的Yplus值，本文中圖3（a）第一層網格高度取為2mm（約為0.005B），而圖3（b）第一層網格高度取為0.25mm（約為0.000 6B）；兩套網格的網格數量分布見表1.兩套網格在數量上比較，只有剛性網格區的網格數量不同，這是由于使用SST模型時需要在近壁區加密網格.兩套網格中，網格數量都集中在動網格區，這是因為經過反復試算表明：適當增加動網格區的網格數量，可以有效消除氣動力系數時程曲線中的“毛刺”.

表1 計算域各區網格數量Tab.1 Meshes in each region for models

計算區域的網格劃分在Gambit中實現，然后分別導入到Fluent中進行數值計算.根據上面對兩湍流模型使用要求的分析，分別控制Yplus值.由圖4可知，斷面Yplus值符合兩湍流模型的使用要求（見第1部分分析）.

圖3 橋梁斷面近壁網格及局部Fig.3 Girds arrangement around bridge section and close view

圖4 橋梁斷面Yplus值分布Fig.4 Yplus distribution around bridge section

3.2 計算條件

3.2.1 橋梁斷面運動模式

采用單自由度單頻等幅正弦位移激勵橋梁斷面運動.對純豎彎運動，有扭轉自由度α（t）＝0，橋梁斷面豎彎運動位移是：

對于純扭轉運動，豎彎運動位移h（t）＝0，橋梁斷面扭轉運動位移是：

式中：fh，fα分別表示豎彎運動頻率和扭轉運動頻率；h0，α0分別表示豎彎運動和扭轉運動幅值，本文統一h0和α0的取值，h0取為0.025B，α0取為3°，小幅振動以滿足線性小擾動假設.

3.2.2 邊界條件和其他相關參數

計算域左側進口為模擬大氣邊界層速度和紊流度的速度入口，紊流度為5%；計算域右側為出流邊界條件，對應沿出口邊界法向速度梯度為零.識別不同折算風速下的顫振導數，本文通過改變強迫振動頻率fh和fα來實現無量風速Vr＝U／fB的改變，這樣保證了不同折算風速下的Re數相同，采用的非穩態計算時間步長為0.005s.

本文采用了FLUENT軟件中的動網格技術，主要參數設置有：采用Smoothing和Remeshing兩種動網格更新方法；網格光滑更新迭代次數設為200，彈性常數因子和邊界節點松弛都設為0.6；局部網格重劃分中網格最大畸變控制為0.4，網格尺寸重劃分迭代次數設為100.

4 速度場比較

為了更直觀地說明兩湍流模型的計算特點和比較兩湍流模型的計算精度，本文給出橋梁斷面近壁區在nT時刻處的速度矢量圖，T為模型強迫振動周期，nT時刻表示模型回到平衡位置的時刻.在模型強迫振動過程中，沒有看到明顯的漩渦脫落；同時對氣動升力進行頻譜分析，可以看到在主頻率中，只有模型進行周期運動的頻率.

如圖5（b）所示，使用SST湍流模型時，在橋梁斷面上下緣流線分離處成功地捕捉了渦，且在近壁區處，速度大小沿斷面法線方向呈現梯度變化，圖5（a）卻無法看到這些現象.這可以解釋為：標準k－ε模型在近壁區采用了壁面函數，在模擬近壁區繞流、流線分離和渦的形成時，出現了一定程度的失真.SST湍流模型在計算近壁區流動時，采用了標準kω模型，由于它屬于一種低Re模型，因此有效地提高了近壁區的計算精度.可見，在模擬橋梁斷面繞流和流線分離時，SST湍流模型的計算優勢明顯.

圖5 近壁區速度矢量分布前緣及局部放大Fig.5 Velocity vector around bridge section and close view

5 氣動導數識別結果

本文采用最小二乘法識別了基于兩湍流模型下不同折算風速的顫振導數，并將Poulsen［9］的試驗值同時列入以作參考，如圖6和圖7所示.由于該試驗值是在均勻流中通過橋梁模型自由振動實現識別獲得，而本文則考慮了紊流場，因而與本文結果沒有嚴格數量上的可比性，但可作為定性參考.采用數值模擬的識別結果與試驗值趨勢大體一致（除外），氣動導數，，，和的數值模擬值與試驗值相差不大，符合實際工程要求，因此，基于CFD商業軟件FLUENT，采用標準k－ε湍流模型和SST湍流模型來識別大帶東橋的氣動導數，都能取得較好的效果；本文采用SST湍流模型成功地預測了和曲線的趨勢，相關文獻中對這些顫振導數的識別卻不理想，因此相比而言，SST湍流模型的數值模擬精度優于標準k－ε湍流模型.

圖6 計算的大帶東橋顫振導數（1）Fig.6 Simulated flutter derivatives of the Great Belt Bridge（1）

圖7 計算的大帶東橋顫振導數（2）Fig.7 Simulated flutter derivatives of the Great Belt Bridge（2）

6 結 論

本文針對基于標準k－ε湍流模型和SST湍流模型的大跨橋梁顫振導數數值識別及比較研究，得到下述結論：

1）本文提出的計算域分區網格劃分方法，以及為保證內部剛性網格區網格正交性的橢圓外形邊界，提供了橋梁斷面繞流近流場模擬的較好計算網格，為標準k－ε湍流模型和SST湍流模型較準確地識別橋梁斷面的大部分顫振導數提供了保障.

2）雖然標準k－ε湍流模型和SST湍流模型均能較準確地識別橋梁斷面的大部分顫振導數，但在個別顫振導數識別上，標準k－ε湍流模型無法給出趨勢性的結果，表明SST湍流模型比標準k－ε湍流模型具有明顯的優勢.

3）不同的湍流模型有不同Yplus要求，網格尺度和布置必須滿足各個湍流模型的相應要求.

4）本文考慮了紊流場對顫振導數的影響，大部分導數數值模擬結果略大于試驗值，但紊流場中的各個參數（包括紊流強度、積分尺度等）對各個顫振導數的具體影響，還必須結合相關的風洞試驗和數值模擬進行精細化研究.

［1］ LARSEN A，WALTHER J H.Aeroelastic analysis of bridge girder sections based on discrete vortex simulations［J］.Journal of Wind Engineering and Industrial Aerodynamics，1997（67／68）：253－265.

［2］ ZHU Zhi－wen，GU Ming，CHEN Zheng－qing.Wind tunnel and CFD study on identification of flutter derivatives of a long－span self－anchored suspension bridge［J］.Computer－Aided Civil and Infrastructure Engineering，2007，22：541－554.

［3］ 曹豐產，項海帆，陳艾榮.橋梁斷面顫振導數和顫振臨界風速的數值計算［J］.空氣動力學學報，2000，18（1）：26－33.CAO Feng－chan，XIANG Hai－fan，CHEN Ai－rong.Numerical assessment of aerodynamic derivatives and critical wind speed of flutter of bridge decks［J］.ACTA Aerodynamic Sinica，2000，18（1）：26－33.（In Chinese）

［4］ LAUNDER B E，SPALDING D B.Lectures in mathematical models of turbulence［M］.London：Academic Press，1972.

［5］ VERSTEEG H K，MALALASEKERA W.An introduction to computational fluid dynamics：the finite volume method［M］.New York：Wiley，1995.

［6］ MENTER F R.Two－equation eddy－viscosity turbulence models for engineering applications［R］.New York：AIAA，1993.

［7］ HUANG Lin，LIAO Hai－li.Numerical simulation for aerodynamic derivatives of bridge deck［J］.Journal of Wind Engineering and Industrial Aerodynamics，2009，17（4）：719－729.

［8］ LIAW Kai－fan.Simulation of flow around bluff bodies and bridge deak sections using CFD［D］.Nottingham：University of Nottingham，School of Civil Engineering，2005.

［9］ POULSEN N K，DAMSGAARD A，REINHOLD T A.Determination of flutter derivatices for the great belt bridge［J］.Journal of Wind Engineering and Industrial Aerodynamics，1992，41（1／2／3）：153－164.