Smagorinsky常數(shù)對扁平箱梁氣動特性大渦模擬的影響研究

祝志文，劉帥永，劉震卿

（1.汕頭大學土木與環(huán)境工程系，廣東 汕頭 515063；2.華中科技大學土木工程與力學學院，湖北 武漢 430074）

引 言

扁平鋼箱梁不僅承載力大、結(jié)構(gòu)輕，而且因結(jié)構(gòu)高度小，導致風荷載小，是大跨度橋梁常用的主梁結(jié)構(gòu)形式，如廣州南沙大橋、丹麥大帶東橋、潤揚大橋和江陰大橋等均采用扁平鋼箱梁。在大跨度橋梁抗風中，扁平箱梁氣動特性的獲得目前主要通過節(jié)段模型風洞試驗，但風洞試驗涉及模型制作，試驗周期長、費用高，且模型縮尺比往往較小，使得風洞試驗的雷諾數(shù)往往遠低于實際橋梁繞流Re數(shù)。計算流體力學（CFD）基于計算機和數(shù)值方法求解流體動力學問題，在風工程研究中應用廣泛。隨著計算機特別是并行計算的發(fā)展，CFD 有望為風工程研究提供一種可能代替風洞試驗的手段。與風洞試驗相比，CFD 模擬無需模型制作和測量設(shè)備，能回避或解決風洞試驗中存在的諸多問題，且能形象而又細致地再現(xiàn)風的復雜流動，能借助其強大的后處理和流場可視化開展微觀分析，便于揭示流動機理。另外，CFD 模擬能夠生成滿足指定風場特性的入口來流，方便修改流場參數(shù)，可不受測力和測壓模型限制獲得氣動力，能夠靈活布置大量表面壓力采集點，也無采集管路系統(tǒng)頻響特性對風壓信號動態(tài)特性的影響。

與湍流模擬的其他CFD 方法相比，LES 采用過濾方法將流場變量分解成大尺度脈動和亞格子脈動兩部分，直接計算大尺度的脈動，亞格子脈動通過湍流模型模擬，其對流動的非定常特性有較高的捕捉能力。LES 最早由Smagorinsky［1］于1963年提出，隨著計算能力及數(shù)值方法的不斷改進，LES 在20世紀90年代后獲得了巨大發(fā)展，廣泛應用于復雜流動的模擬中，橋梁風工程方面的相關(guān)研究文獻也越來越多。

作為一種經(jīng)驗湍流模型，LES 目前存在的問題是Smagorinsky 常數(shù)CS的不確定：CS的取值并非定值，可能與流場類型、雷諾數(shù)及其他無量綱參數(shù)相關(guān)［2］。為此，學者對CS的合理取值進行了廣泛研究。對于截斷波數(shù)處于慣性子區(qū)范圍且過濾尺度與網(wǎng)格尺寸相等的均勻各向同性湍流，Lilly［3］基于能譜分析得到的CS取值為0.23。Deardorff［4］在槽道湍流的模擬中發(fā)現(xiàn)采用文獻［3］建議的CS值會導致亞格子應力黏性過大，而當CS取0.1 時，模擬結(jié)果與Laufer［5］試驗相近。Kwak 等［6］、Shaanan 等［7］、Ferziger 等［8］和Antonopoulos-Domis［9］基于Smagorinsky模式得到了與Comte-Benot 和Corrsin［10］試驗一致的能量衰減率，其CS取值范圍為0.19～0.24。Mason和Callen［11］研究認為CS取值與網(wǎng)格分辨率有關(guān)，當網(wǎng)格分辨率足夠高時，CS取0.2 能給出較好的結(jié)果，反之，如網(wǎng)格分辨率不高，CS的取值須小于0.2。Piomelli 等［12］研究表明即使網(wǎng)格分辨率比Mason 等［11］所使用網(wǎng)格高得多，CS的最佳取值仍為0.1 左右。Germano 等［13］提出了一種動態(tài)亞格子模型，將C2S作為時間和空間的函數(shù)，用不同特征長度的兩個濾波函數(shù)來確定?？梢娔壳斑€沒有解決CS的合理取值問題，而現(xiàn)有橋梁主梁LES 采用CS=0.1［14］，其合理性也是存疑的，因為繞扁平箱梁的流動，與自由剪切流和槽道湍流差別較大，包含鈍體繞流的分離、沖撞、自由剪切層和旋渦脫落等各種不同的流動結(jié)構(gòu)。

因此，研究LES 參數(shù)CS的合理取值，并提出相應建議，對采用LES 獲得扁平閉口箱梁的氣動特性有重要的工程意義，但目前無相關(guān)研究報道。

1 LES 的基本原理和控制方程

LES 的基本思想是，湍流由不同尺度的漩渦組成，大尺度的渦旋對湍流能量和雷諾應力的產(chǎn)生及流動量的擴散起主要作用。大渦的行為強烈依賴于邊界條件，隨流動的類型而異。小渦的上述貢獻較小，最小尺度的渦主要起耗散作用。高雷諾數(shù)下小渦近似于各向同性，受邊界條件影響較小，具有較大的共性［15］。雖然目前的計算機還不能計算到耗散尺度，但能夠小到慣性區(qū)尺度，所以可通過離散非定常Navier-Stokes（N-S）方程來確定大渦的行為，而用較通用的模型模擬小渦的作用。

LES 將N-S 方程中流場速度變量ui變成大尺度可直接求解的變量，加上過濾掉的亞格子尺度量u′i，即：

式中第一項表示大于過濾尺度的大渦量，是在整個計算域D上通過下述卷積積分獲得：

式中xj和x′j均為空間坐標向量；Gj為過濾函數(shù)，并滿足：

在有限差分法和有限體積法空間離散中，最常采用的是體積加權(quán)的盒函數(shù)過濾器，即：

式中Δj(j=1，2，3)為坐標軸方向的網(wǎng)格尺度，過濾得到的不可壓N-S 方程形式為：

式中ρ和υ分別為空氣密度和運動黏性；為過濾后的壓力。因無法同時求解和，需要將作分解，也即τij為亞網(wǎng)格剪應力張量，τij=為亞格子湍流黏性；δij為狄拉克函數(shù)；Sˉij為根據(jù)求解的大渦量確定的應變率張量，可表示為：

這樣，包含連續(xù)方程在內(nèi)的不可壓N-S 方程可表示為：

對LES 常采用的Smagorinsky 模型［3］，則亞格子湍流黏性可表示為：

與雷諾數(shù)均等時其他湍流模型需要多個經(jīng)驗化的模型參數(shù)相比，LES 僅需一個經(jīng)驗參數(shù)CS，因此，LES 的模型參數(shù)數(shù)量顯著減小。另外，從式（8）可見，當CS增加時，亞格子湍流黏性將增大，此時流動的耗散作用將增強，將導致非定常流動中亞格子尺度小渦的作用減弱，改變湍流中不同尺度渦之間能量級聯(lián)方式，降低LES 對流動非定常特性的捕捉能力，比如導致氣動力的脈動量減小、流動頻譜中高頻分量的衰減或缺失、頻譜變窄，也可能會改變漩渦脫落特性。此外，對流場的細部結(jié)構(gòu)，如主梁表面的流動分離、再附、剪切層運動和復雜尾跡流動均可能產(chǎn)生影響。

本文LES 基于Fluent 6.3.26 軟件，以大帶東橋主梁為對象［16］，根據(jù)LES 獲得的氣動特性與試驗結(jié)果的對比，探討扁平箱梁氣動特性LES 在CS上的合理取值。

2 數(shù)值方法和邊界條件

合理的CS取值需對比LES 結(jié)果與風洞試驗結(jié)果。本文以大帶東橋主跨加勁梁標準斷面施工階段為對象，不考慮橋面防撞欄等附屬設(shè)施。實橋主梁全寬31 m，高4.4 m，主梁扭轉(zhuǎn)中心（S.C.）位于橋軸線主梁底板以上2.35 m 處。CFD 模型縮尺比為1∶80，對應的模型截面寬B=0.3875 m，高D=0.055 m，順橋向模型長度取一倍模型截面寬度，也即L=B，主梁CFD 橫截面如圖1所示。

圖1 CFD 主梁橫斷面尺寸和表面監(jiān)測點布置（單位：m）Fig.1 Girder cross section and surface monitoring points in CFD（Unit：m）

圖2（a）為計算域和主梁網(wǎng)格布置示意圖，其中入口、上下側(cè)邊界到主梁S.C.點的距離均為10B，下游出口到S.C.點的距離為20B，計算域展向長度等于模型長度L。模型在計算域內(nèi)的堵塞度為0.7%，遠低于3%的CFD 模型堵塞度要求。

圖2 計算域和網(wǎng)格布置Fig.2 Computational domain and mesh arrangement

采用結(jié)構(gòu)和非結(jié)構(gòu)六面體網(wǎng)格控制計算域網(wǎng)格數(shù)量和質(zhì)量，并控制同一網(wǎng)格方向相鄰網(wǎng)格的高度比不大于1.1，在流動變量變化梯度大的風嘴前緣和主梁表面折角附近加密網(wǎng)格，在物面法向采用精細的結(jié)構(gòu)化六面體網(wǎng)格劃分，并保證網(wǎng)格的正交性，提高模型尾流區(qū)網(wǎng)格分辨率。在流動變量變化平緩的區(qū)域采用較粗的網(wǎng)格，以降低網(wǎng)格總量。

來流風速U0= 5 m/s，風攻角為0°。計算域邊界條件，在入口設(shè)置為均勻流速度邊界；下游出口施加出流邊界；主梁表面采用無滑移壁面條件；計算域的頂面和底面，以及模型展向兩側(cè)面均設(shè)為對稱邊界，如圖2（a）所示。計算域初始場采用入口速度條件初始化。開展了LES 模擬的網(wǎng)格無關(guān)性和時間步無關(guān)性檢查，因篇幅限制，本文不再列出。綜合比較計算量和模擬精度要求，確定模型表面法向第一層網(wǎng)格高度為B/1300，時間步長為0.0002 s，模型表面Y+滿足LES<1 的要求。計算域網(wǎng)格總數(shù)為986544，對應的流動雷諾數(shù)為1.33×105。

3 CS取0.1 時主梁的氣動特性

CFD 計算在模型表面布設(shè)了160 個壓力監(jiān)測點，以捕捉不同時刻模型表面的壓力分布，如圖1所示。定義監(jiān)測點的壓力系數(shù)CP為：

式中P為監(jiān)測點壓力；CP的平均值和RMS 值分別加下標m 和RMS 表示。

定義氣動力作用點在圖1截面的S.C.點，分別定義模型斷面氣動阻力、升力和扭矩系數(shù)為：

式中FD，F(xiàn)L和M分別為作用在主梁上的阻力（順來流方向為正）、升力（垂直來流方向向上為正）和扭矩（順時針方向為正）；三分力系數(shù)的平均值和脈動RMS 值分別再加下標m 和RMS 表示。

定義主梁漩渦脫落St數(shù)：

式中fs為主梁漩渦脫落頻率（Hz）。

以后緣壓力監(jiān)測點90為參考，圖3給出了扁平箱梁緊靠6個棱角的壓力系數(shù)時程曲線?？梢娺@些時程曲線均表現(xiàn)出頻率相同但幅值不同的波動。除迎風側(cè)上斜腹板監(jiān)測點1風壓系數(shù)為正外，其他均為負；下表面風壓脈動明顯大于上表面，最大風壓脈動出現(xiàn)在迎風側(cè)下斜腹板與底板相交處，其次是迎風側(cè)下斜腹板緊靠前緣的點。另外，從時程峰值對比來看，監(jiān)測點壓力并不同步，部分同相、部分反相。

圖3 主梁表面典型監(jiān)測點風壓系數(shù)時程Fig.3 Pressure coefficient time histories at typical monitoring points on girder surface

圖4（a）為氣動三分力系數(shù)時程。可見阻力系數(shù)時程呈現(xiàn)隨機波動特性，阻力時程對應的平均阻力系數(shù)為0.435，比基于剛性模型測壓積分獲得的阻力系數(shù)0.352 高18%［16］；LES 捕捉到了較大的升力系數(shù)波動，呈現(xiàn)出規(guī)律的波動特征，隱含繞流明顯的渦脫特征。圖4（b）升力時程的功率譜密度（PSD）分析表明，顯著峰值占優(yōu)的頻率對應的St數(shù)為0.29，也即為單一頻率渦脫，這與風洞試驗結(jié)果完全吻合［16］。

圖4 氣動力系數(shù)時程和渦脫St譜Fig.4 Aerodynamic coefficients time histories and vortex-shedding spectrum

圖5為從升力時程峰值時刻開始，在一個渦脫周期T內(nèi)主梁周圍壓力云圖。模型后緣上下表面出現(xiàn)明顯的交替脫落的漩渦，從而在主梁表面誘導出波動的氣動力，如圖4（a）所示。國內(nèi)外多座大跨度扁平鋼箱主梁橋梁出現(xiàn)了渦激振動，圖5的流動可視化有助于理解扁平箱梁渦脫產(chǎn)生的流動機理。

圖5 一個渦脫周期T 內(nèi)的壓力云圖Fig.5 Pressure contours around girder in one vortex-shedding cycle

4 不同CS值對主梁氣動特性的影響

開展了CS在0.032～0.70 范圍取值以及動態(tài)亞格子黏性模型的多個LES 計算，其他計算設(shè)置和條件均與上節(jié)相同。CS值從0.032 開始逐漸增大，為提高計算效率，后一個CS值計算工況在前一個CS值收斂結(jié)果上開展，由于CS值前后差別很小，前一次CS值計算的收斂結(jié)果是后一個CS值計算工況的良好初始場。每個CS值計算工況，在計算殘差收斂平穩(wěn)后，再計算65 個渦脫周期［17］，得到1.3×104個時間步上的主梁氣動力時程和監(jiān)測點風壓時程。

不同CS取值時LES 得到的主梁氣動力平均值和RMS 值如圖6所示?？梢娏ο禂?shù)平均值變化不明顯，阻力系數(shù)隨CS值的增大先有小幅增大，隨后小幅減小，最大變化量不大于10%；升力系數(shù)平均值也是先小幅增大再減小，最大差值約0.05。力系數(shù)RMS 值隨CS值的增大均呈現(xiàn)減小的趨勢，特別是升力系數(shù)RMS 值最顯著；當CS值大于0.50 后，力系數(shù)RMS 值為零，也即其時程將為一條水平直線。這說明，隨著CS值的增大，亞格子湍流黏性增大，此時流動的耗散作用增強，使得LES 對主梁繞流非定常特性的捕捉能力降低，導致氣動力脈動減小。

圖6 主梁力系數(shù)平均值和RMS 值隨CS的變化Fig.6 Mean and RMS values of aerodynamic coefficients of girder under various CS values

為探討CS取值對主梁漩渦脫落的影響，本文對不同CS取值獲得的升力時程，開展了PSD 分析，并將橫軸換算為渦脫St數(shù)，可得到相應CS取值工況下的主梁渦脫St譜。為區(qū)分主梁渦脫的單頻（頻譜只有一個峰值）和多頻（頻譜有多個明顯的峰值）特征，渦脫St譜也需反映單頻和多頻特征（風洞試驗St=0.28～0.29，為單頻［16］），本文對渦脫多頻的情況給出了前2 個峰值對應的St數(shù)。圖7給出了主梁渦脫St數(shù)隨CS值增大的變化，結(jié)合圖4（b），可見當0.064≤CS≤0.45，主梁渦脫均呈單峰特征，且渦脫頻率隨CS增大而小幅減小。圖8另外給出了多個CS值下的渦脫St譜，可見在上述范圍之外，主梁渦脫呈現(xiàn)多峰特征；參考試驗結(jié)果［16］，可見在0.064≤CS≤0.27 范圍內(nèi)，LES 均能給出主梁渦脫St數(shù)的合理估計。另外，采用動態(tài)亞格子黏性模型，LES 無法給出主梁渦脫的合理估計，如圖8（d）所示，這與Ferziger 等［18］的結(jié)論一致。

圖7 主梁渦脫St數(shù)隨CS值的變化Fig.7 Girder vortex-shedding St versus various CS

圖8 不同CS值主梁渦脫St譜Fig.8 Girder vortex-shedding St spectrum under various CS

從圖7的渦脫St譜特征考慮，開展了0.032≤CS≤0.50 和動態(tài)亞格子黏性模型的LES 計算，將主梁表面壓力系數(shù)平均值和RMS 值與風洞實驗進行了對比。由圖9和10 可見，對全部CS值，表面壓力系數(shù)平均值與風洞試驗結(jié)果，除在靠近前緣附近2個監(jiān)測點外，其他吻合良好。在頂?shù)装逵L側(cè)折角位置處，LES 給出了比試驗更大的極值負壓估計，這可能與CFD 布置的監(jiān)測點更靠近折角有關(guān)。實際上，因風洞模型制作困難，測壓管無法像CFD 一樣能緊靠模型折角。

圖9 主梁上表面壓力系數(shù)平均值Fig.9 Mean pressure coefficients on girder upper surface

圖11和12 為0.032≤CS≤0.50 和動態(tài)亞格子黏性模型下，LES 獲得的主梁表面壓力系數(shù)RMS 值與風洞實驗的對比?？梢妼θ緾S值，壓力系數(shù)RMS 值均與風洞試驗結(jié)果存在差異，表現(xiàn)為LES低估了壓力脈動值，特別是在主梁的中后部，后緣附近的差別較大。上述結(jié)果與風洞試驗的差別，還隨著CS值的增大而進一步增大。動態(tài)亞格子黏性模型能得到較大的壓力系數(shù)RMS 值，但也未能獲得與試驗結(jié)果相一致的趨勢。顯然，CS值的增大不利于LES 捕捉流場的脈動量。

圖10 主梁下表面壓力系數(shù)平均值Fig.10 Mean pressure coefficients on girder lower surface

圖11 主梁上表面壓力系數(shù)RMS 值Fig.11 RMS pressure coefficients on girder top surface

5 結(jié) 論

圖12 主梁下表面壓力系數(shù)RMS 值Fig.12 RMS pressure coefficients on girder bottom surface

本文開展了CS在0.032～0.70 范圍以及動態(tài)亞格子黏性模型的LES 多工況模擬，獲得了主梁氣動力、漩渦脫落St數(shù)、表面壓力系數(shù)平均值和RMS 值，并與風洞試驗結(jié)果進行了對比，基于數(shù)據(jù)分析和討論，得到下述結(jié)論：

1）CS取0.1，LES 能捕捉到主梁氣動力和漩渦脫落的非定常特征，估計的St數(shù)與風洞試驗吻合。

2）CS值對主梁表面壓力系數(shù)和氣動力系數(shù)的平均值影響均不明顯，且其結(jié)果與風洞試驗吻合良好；CS值對壓力系數(shù)和氣動力系數(shù)的RMS 值影響明顯。隨著CS增大，亞格子湍流黏性增大，流動的耗散作用增強，使得LES 對流動的非定常特性捕捉能力降低，導致氣動力的脈動量減小。當CS值大于0.5 后，氣動力將不再出現(xiàn)脈動。

3）在頂?shù)装迩熬壵劢俏恢锰?，LES 給出了比試驗更大的平均負壓極值估計。壓力系數(shù)的RMS 值和分布均與風洞試驗結(jié)果存在差異，且差異隨CS值的增大而增大。動態(tài)亞格子黏性模型能給出較大的壓力系數(shù)RMS 值和與風洞試驗不同的分布，也無法給出主梁渦脫的合理估計，不建議采用。

4）在0.064≤CS≤0.27 范圍內(nèi)，LES 均能給出與風洞試驗一致的主梁渦脫單頻St值估計，從捕捉流動的非定常特性考慮，建議CS在上述范圍內(nèi)取小值。

另外，LES 獲得的壓力系數(shù)RMS 值和分布均與風洞試驗結(jié)果存在差異，是否與LES 計算域入口來流條件有關(guān)，值得進一步研究。