機動式高架棧橋在風浪耦合作用下的位移響應

李 霞，孫蘆忠，尹洪波，謝興坤

（1.國防科學技術大學，長沙 410072；2.解放軍理工大學 工程兵工程學院，南京 210007；3.二炮工程設計研究院，北京 100011）

機動式高架棧橋是一種可移動到任意位置架設的棧橋系統。一方面，由于橋面不與海面接觸，可以經受比一般浮游棧橋更復雜的海況。另一方面，高架棧橋與普通大型橋梁有顯著的不同，普通大型橋梁由于考慮通航，橋面下的凈空高度往往達數十米，正常海況下，波浪對橋下風場的影響可忽略。尹洪波研究了凈空高度對高架棧橋下風載的影響［1］，表明了凈空高度以及波浪對橋下風載的敏感性，即低凈空下波浪對風載的影響不可忽略。對于機動式高架棧橋，只需滿足高潮位時，橋面底部不受波浪抨擊即可，故橋下凈空較小。鑒于該特點，波浪對高架棧橋橋下風場影響顯著，即風、浪對高架棧橋動力響應的貢獻不是孤立的，必須同時考慮風、浪與結構的耦合受力關系。因此，高架棧橋在受波浪影響的強風作用下的受力分析及其位移響應是一個值得深入研究的問題。

由于問題的難度和特殊性，考慮風浪耦合作用下橋梁動力分析方面的文獻不多。已有文獻大多針對系泊系統、海洋平臺等海洋工程結構的波浪動力應［2］，或者先計算出風和波浪單獨作用時結構動力響應，然后疊加，以此作為風浪聯合作用下結構響應的近似值［3，4］。錢景峰對隨機風浪下的橋墩、橋塔結構的動力響應進行了計算，即同時將風力和波浪力施加于結構上，但還是沒有考慮風與波浪本身的耦合相關性［5］。目前亦有863課題資助的“臺風浪耦合作用下跨海峽橋梁動力模擬及防災減災技術”研究，劉高等在文獻［6］中對于跨海峽橋梁、臺風以及臺風掀起的巨浪之間的耦合進行了定性說明，但尚未見報道詳情的文獻。

風、浪與結構構成的流場是典型的湍流問題，其數值模擬方法主要有直接數值模擬（DNS）、雷諾平均N－S方程法（RANS）、大渦模擬（LES）和分離渦模擬方法（亦稱脫體渦方法，DES）等。其中 DES是 Spalart［7］于1997年提出的將RANS與LES相結合的混合方法，其基本思想是：在近壁區采用RANS湍流模型，充分利用了RANS可有效處理附體邊界層流動的優點；而在遠離物面的分離區內采用LES方法，小尺度渦采用亞格子模型來模擬，大尺度渦則對其直接進行模擬，故DES吸收了RANS和LES各自的優點。目前廣泛應用于基礎研究和工程領域［8－10］。

高架棧橋受到風、浪激勵而發生振動是一個復雜的流固耦合問題。目前流固耦合數值解法有兩大類：緊耦合和松耦合。緊耦合法方程規模龐大、計算量大，有學者用其計算了橋墩的地震反應［11］。松耦合中的交替法是將流體和固體分成兩個單獨的求解域，在數值求解過程中交替地求解這兩個區域，并通過耦合界面進行有關物理量的傳遞從而達到不同求解域的相互耦合。交替法求解的方程規模小，是目前解決流固耦合問題比較理想的方案，如方平治等對橢圓柱體渦激振動的數值模擬［12］，本文采用交替求解方法。

1 控制方程與數值方法

1.1 流體域控制方程—DES方法

本文采用基于Menterk－ωSST（Shear Stress Transport）兩方程湍流模型［12］的DES方法，模型方程為：

式中，ρ為空氣密度，k為湍流動能，ω為湍流比耗散率，u為來流速度，F1為混合函數，μt為渦粘性系數，Pk、Pω為湍流生成項，其具體定義及模型方程中的常數根據參考文獻［13］給出。

以上k－ωSST湍流模型k方程的耗散項中，湍流尺度參數lk－ω的表達式為：

在 DES 方法中，lk－ω=min（lk－ω，CDESΔ），其中 Δ =max（Δx，Δy，Δz），為網格單元最大邊長，常數CDES=0.61。lk－ω控制著方程在 LES和 RANS之間的轉換。在靠近壁面的邊界層中，ω相當大，而湍流動能k值有限，此時lk－ω?Δ，湍流模型為k－ωSST模型，是雷諾平均RANS算法；當遠離邊界層時，ω迅速減小，這時lk－ω增長非常迅速，當lk－ω?Δ 時，Smagorinski大渦模擬模型（亞格子雷諾應力模型）起作用。因此，本方法可以在LES和RANS之間自適應，在網格密度達到要求的地方使用LES，在網格較粗的地方則使用RANS，從而可用較少的計算網格得到與大渦模擬相近的計算精度。

1.2 固體域控制方程

高架棧橋在隨時間變化的荷載作用下，滿足結構動力學的基本運動方程

其中，M為結構的質量矩陣，C為阻尼矩陣，K為剛度矩陣，為節點加速度向量，為節點速度向量，u為節點位移向量。

1.3 波浪與風的模擬

由于自由表面的存在，分析中采用VOF法來表征氣液兩相流。定義函數C，C和1－C分別代表計算區域內空氣和水占計算區域的體積分數。對于某個計算單元而言，存在下面三種情況：

（1）C=0，表示該單元完全被空氣充滿；

（2）C=1，表示該單元完全被水充滿；

（3）0＜C＜1，表示該單元部分是空氣，部分是水，有氣液交界面，自由表面屬于第三種情況。

體積函數C滿足

u為該單元的平均速度。對于兩相流，粘性系數μ和密度ρ均由體積函數C決定，ρ=ρ1C+ρ2（1－C），μ=μ1C+μ2（1－C），其中下標1和2分別代表空氣和水。

對于波浪的模擬，本文采用二階Stokes波浪理論，讓入口邊界處流入的水滿足該理論近似解析解如式（6）。

式中，u，w為速度分量，ζ為水面變化，θ為位相函數，K為波數，H為波高，λ為波長，h為水深。

對于風，采用Davenport提出的指數率模型來描述平均風速剖面，風速隨高度變化，即：

式中，U10為10 m高度處50年一遇10 min最大平均風速；Z為任意一點距地面高度，a為地面粗糙度。本文中取U10=15 m/s，a=0.12。

1.4 流固耦合求解方法

為全面描述棧橋結構的耦合位移響應，本文采用基于ALE的網格運動算法。ALE描述下的不可壓流體方程為：

式中，wj為網格運動的速度分量，ALE描述和Euler描述的區別僅在于將速度uj用相對速度uj－wj代替。若wj為零，即為Euler系；若wj等于質點速度，即為Lagrange系。

在交替求解法中，應用有限體積法和有限單元法分別求解流體域與固體域，在流體和固體交界面設置邊界單元，處理流體域和固體域的數據交換。流體域向固體域傳遞各種力，固體域向流體域傳遞網格位移和速度，從而實現兩個求解域的耦合。

2 數值模擬與分析

2.1 算法驗證

為了驗證本文流固耦合方法的有效性，將柔性方形鈍體結構橫風向振動的數值模擬結果與試驗結果、Tamura［14］的計算結果進行了對比。結構的三維尺度比值為：B∶D∶H=1∶2∶2，其中B、D和H分別為橫風向寬度、順風向長度和高度。簡化風速Vr定義為Vr=U0/（Bf0），其中，U0為來流平均風速，f0為結構自振頻率。采用DES分別模擬了Vr為6、13和20時結構的振動。圖1為三種工況下，無量綱化橫向位移的均方根值與其它結果的比較。圖中顯示，風速較小時，模擬結果與試驗結果接近；隨著風速增加，兩者相差較大，但本文結果比Tamura的模擬結果更好，能夠較好地模擬振幅隨風速增加而增加的趨勢。

圖1 橫風向位移隨簡化速度的變化Fig.1 Across-wind displacement changes with simplified velocity

2.2 計算模型

（1）固體域

為減少計算量，取圖2（a）所示結構為算例，橋由上部的箱型結構和下部的橋柱組成。箱型結構外部尺度寬4.5 m，長 12 m，高 1 m，鋼板厚 8 mm；內部橫向每1.5 m、縱向每3 m設有縱橫加肋板，板厚6 mm。橋柱高8 m，外徑25 cm，板厚8 mm，柱底為固支。整個模型采用同一種鋼材，其彈性模量E=2.0×1011Pa，泊松比μ=0.3，密度ρ=7 800 kg/m3。結構的固有頻率與模態見表1，其中，x，y，z軸正向分別對應高架棧橋沿橋長方向、沿橋寬方向和沿橋高方向。

（2）流體域

流體計算區域長100 m，寬40 m，高30 m，高架棧橋距風浪入口30 m，水深5 m，橋下凈空3 m，水面以上為風入口，如圖2。波浪模擬采用二階Stokes波浪模型（式（6）），波長取為λ=15 m，風場模擬采用式（7），取U10=15 m/s，α=0.12。由于三維流固耦合計算量巨大，采用 LES所需網格數極大，故流體域計算采用DES。網格為結構化的六面體網格，單元總量為3 261 356。在棧橋附近和水面附近加密網格以提高精度，如圖3。

表1 結構各階模態頻率和周期Tab.1 The frequencies and periods of structure’s different modes

圖2 計算模型示意圖Fig.2 Sketch of computing model

圖3 局部網格Fig.3 Partial mesh

2.3 算例及分析

根據上述模型，分別計算了波浪周期T為5 s、8 s和10 s時棧橋的運動響應，對于每種波周期分別考慮了三種工況：① 波高H=0?m，風速U10=15 m/s，即棧橋只受風載；② 波高H=1.2?m，U10=0 m/s，即棧橋只受波浪荷載；③ 波高H=1.2?m，U10=15 m/s，即棧橋受到風浪耦合作用。計算時間取100s，時間步長為0.01 s。為了跟蹤高架棧橋的位移響應，將其迎風面中心點A設為觀察點（如圖2（a）所示），記錄其位移時程，結果見圖4－圖8。

從圖4中可看出，棧橋僅受風載時，A點在偏離原平衡位置1.5 cm處做周期為8 s左右的緩慢衰減振蕩，該周期與結構一階模態的周期相當；當棧橋僅受到周期為5s的波浪荷載時，A點在x方向的振幅在零值左右振蕩，其響應周期明顯受到波浪周期的影響。而圖5表明，當棧橋受到風與周期為5 s的波浪耦合作用時，A點x方向的位移幅值明顯大于其受單個荷載時的位移幅值，而且最大位移分別是僅受波浪荷載時的2.06倍和僅受風載時的1.16倍，且相當于兩種荷載單獨作用下最大響應之和的74.2%。這說明當波浪周期為5 s時，風浪耦合的結構位移小于僅受風荷載和僅受波浪荷載時位移的簡單疊加，且幅度不小。造成這種現象的原因是：波浪擾亂了棧橋附近的風場，使得作用在棧橋上的平均風壓降低，這與對高架棧橋斷面的三分力系數分析的結果一致［15］。三種工況中A點y，z方向位移隨時間的變化與x方向規律一致，不過位移幅值的數量級為10－3m，比x方向小一個數量級。

對于棧橋的轉動，設其繞x，y，z軸的轉角分別為α，β，γ。從圖6可看出，高架棧橋的轉動很小，β的量級為10－3度，三種工況中α和γ的時程變化規律均與β一致，但小一個量級，這說明棧橋的運動主要為平動。不過，風浪耦合對于棧橋轉動的影響與平動不同，兩種荷載所致的角位移響應之和的幅值與風浪耦合時基本相當，只是相位有所差別，因此，風浪耦合主要是對線位移的影響較大。

當波浪周期為8 s時，結構發生了明顯共振，在100 s內A點在x方向的位移幅值迅速增加到60 cm以上，且受到風浪耦合作用時的最大位移要比僅受兩種荷載單獨作用時的最大位移之和大82.8%，圖7也非常清晰地看到，共振時風浪耦合響應遠大于簡單疊加的響應。同理，A點在y，z方向的響應也呈類似規律，不過幅值均比x方向小一個量級。結構的轉動隨時間也呈共振變化，風浪耦合作用下的最大轉角也比其它兩種工況下的轉角之和大，但幅值最大的β量級僅為10－2度。

波浪周期為10 s時，在波浪荷載或風浪耦合作用下，位移響應出現了明顯的“拍”現象。與共振時不同，風浪耦合作用下結構沿x方向的最大位移，比僅受兩種荷載單獨作用時的最大位移之和小16.7%。A點其它五個自由度方向的位移或角位移，也呈現類似的規律，只是幅值不同；同時，角位移也與波浪周期為5 s的特點類似，兩種荷載耦合對其幅值影響不大。

圖9是結構線位移響應與波浪周期之間的關系，其中縱坐標為棧橋受風浪耦合作用時x方向的最大位移與僅單獨受兩種荷載時最大位移之和的比值δ，橫坐標為波浪周期與結構第一模態固有周期的差值。圖中表明，波浪周期越接近結構的一階固有周期，風浪耦合作用時最大位移與兩種荷載單獨作用時最大位移之和的比值越大，即風浪耦合效應越明顯。

3 結論

本文采用三維分離渦模擬方法，對單跨機動式高架棧橋在僅受風載、僅受波浪荷載及風浪耦合作用下的位移響應進行了數值模擬研究，得到以下結論：

（1）波浪荷載頻率若接近高架棧橋的固有頻率，可引起強烈的共振效應。此時，若風浪耦合作用，則棧橋附近海面的空氣流動受波浪影響，也呈現出與波浪頻率一致的振蕩。該特點在工程設計中應引起足夠重視。

（2）波浪荷載頻率若與高架棧橋固有頻率相差較大，由于波浪在水面附近激起的空氣波動擾亂了棧橋周圍的風場，引起作用在棧橋上的平均風載減小，造成結構在風浪耦合作用下的線位移響應幅值小于結構在兩種荷載單獨作用下的線位移響應之和的幅值。

（3）對于波浪頻率與棧橋固有頻率相差較大的情況，角位移的幅值很小，表明棧橋運動主要是平動；但角位移幅值與線位移規律不一致，風浪耦合對其影響不大。

以上結論說明，用簡單疊加的方法來計算風浪耦合作用下高架棧橋的位移響應將造成較大誤差。

［1］尹洪波.風浪耦合作用下高架棧橋動態響應研究［D］.南京：解放軍理工大學，2009.

［2］徐萬海，曾曉輝，吳應湘，等.深水張力腿平臺與系泊系統的耦合動力響應［J］.振動與沖擊，2009，28（2）：145－150.

［3］Bisht R S，Jain A K.Wind and wave induced behaviour of offshore guyed tower platforms［J］.Ocean Engineering，1998，25（7）：501－519.

［4］ Datta T K，Jain A K.Response of articulated tower platforms to random wind and wave forces［J］.Computers and Structures，1990，34（l）：137 －144.

［5］錢景峰.隨機風浪作用下結構的動力響應分析方法研究［D］.大連：大連理工大學，2010.

［6］劉 高，吳宏波，黃李驥，等.長大橋梁關鍵技術綜述［J］.公路，2009，5：53－64.

［7］Spalart P R，Jou W H，Strelets M，et al.Comments on the feasibility of LES for wings and on hybrid RANS/LES approach.∥ Advances in DNS/LES， First AFOSR International Conference on DNS and LES［C］.Ruston，LA，USA，1997.

［8］Saric S，et a1.A periodically perturbed backward-facing step flow by means of LES，DES and T-RANS：an example of flow separation control［J］.Journal of Fluids Engineering，2005，127（5）：879－887.

［9］ Aroon K V，Danesh K T.Detached eddy simulation of turbulent flow and heat transfer in a ribbed duct［J］.Journal of Fluids Engineering，2005，27：888－896.

［10］ Claus M P，Morton S A，Cummings R M，et al.DES turbulence modeling on the C－130 comparison between computational and experimental results［R］.AIAA 2005 －884，2005.

［11］陳國興，席仁強，王志華.考慮流固耦合的橋墩地震反應方法［J］.防災減災工程學報，2010，30（1）：1－9.

［12］方平治，杜明俠，顧 明.兩自由度橢圓柱體渦激振動的數值模擬［J］.振動與沖擊，2009，28（1）：51 －55.

［13］Menter F R.Zonal two equationk－ωturbulence models for aerodynamic flow［R］.AIAA，93－2906，1993.

［14］ Tamura T，Itoh Y.Three-dimensional vortical flows around a bluff cylinder in unstable oscillation［J］.J.Wind Eng.Indus.Aerodyn.，1997，67 －68：141 －154.

［15］孫蘆忠，尹洪波，李 霞，等.基于大渦模擬的高架棧橋斷面風載分析［J］.港工技術，2008，184（6）：5－9.