低雷諾數變截面細長柔性結構流固耦合能量交換特征分析

鄧秀兵, 于曰旻, 龐璽源

(1. 浙江省水利水電勘測設計院有限責任公司,杭州 310000; 2. 海南大學 土木建筑工程學院, 海口 570228; 3. 上海交通大學 船舶海洋與建筑工程學院,上海 200240)

在土木工程領域,風、浪、流等環境動力荷載作用下的細長結構物受損甚至破壞屢見不鮮.比如,懸索橋中的主纜索作為典型的細長柔性結構,在來風作用下將產生渦激振動,甚至發生疲勞損壞;又如,在主纜索施工過程中,纜索截面形式發生非對稱性變化,當風速超過臨界值后,將產生空氣動力負阻尼,使得纜索振動逐漸增強,甚至超負荷振幅而破壞.因此,科學評估和有效抑制細長結構的流致動力響應是風、浪、流敏感細長柔性結構設計必須考慮的內容.然而,當前細長柔性結構的流固耦合效應方面研究不足,尤其精準模擬和分析細長結構流致振動響應的力學模型和算法方面仍較缺乏,且亟待深入理解其流固耦合機理并建立系統的分析方法與手段[1].

波浪形變截面圓柱體結構具有良好的流動減阻效果.近年來,國內外學者對該類型結構的繞流特性開展了較多研究,發現波浪形變截面能夠延后剪切層的相互作用,從而有效減小阻力.Ahmed等[2]采用實驗方法研究了波浪形圓柱的邊界層分離線和尾流結構,發現波浪形圓柱流動分離線呈現明顯的三維特性,在分離節點附近形成流向渦,且發生邊界層“上卷”現象,從而延遲或抑制剪切層中湍流的生成和發展.Lam等[3]利用多種流動顯示技術對波浪形圓柱的近尾跡進行了實驗研究,給出了平均速度和波動速度分量沿流向、展向和橫向的分布特征.實驗結果表明,波浪形柱渦旋平均形成長度比大于光滑圓柱.湍流統計分析也表明,光滑圓柱尾跡中的渦街更為規則,而由于波浪形圓柱后面的渦具有較強的三維效應,使得由湍流摻混增強而尾跡表現出更強的非相干流動結構.

目前變截面結構繞流分析主要集中在靜止或彈性支承的剛性波浪形結構的減阻和渦激振動特性方面.然而在實際土木工程中,纜索等細長結構的剛度系數通常較小,因此這種結構的柔性作用對流動減阻作用以及渦激振動響應的影響需要進一步研究,包括對初始動力激勵的敏感性方面也需要深入分析.在文獻[4-7]的基礎上,針對這一問題,本研究采用基于高精度譜單元方法的流固耦合分析方法,建立了細長變截面柔性圓柱體結構的渦激振動力學模型,對其在低速均勻流作用和駐波初始動力激勵下的流致振動機理進行了量化分析,獲得了包括波浪形結構的尾流特性、結構動力響應特性、能量傳遞規律、渦脫頻率展向變化特征,并對其在駐波擾動下的減阻、減振機理進行了深入探討,對波浪形變截面細長結構的工程設計和應用提供參考.

1 數值算法

細長結構的流致振動是一種典型的流固耦合問題.求解流固耦合的基本方法通常有兩大類,分別為界面匹配法和非界面匹配法,包括經典的任意拉格朗日-歐拉(ALE)方法和浸入邊界法(IBM).本文采用另一種由Dimas等[8]提出的隨體坐標系(Body-fitted Coordinates)方法,該方法在慣性坐標系中求解N-S方程,再通過坐標轉換至非慣性坐標系中,從而不需要動網格和浸泡邊界近似.

N-S方程和連續性方程在慣性坐標系(x′,y′,z′)內可描述為

(1)

(2)

由慣性坐標系至非慣性坐標系的坐標轉換采用如下關系式:

x=x′-ζx(z,t),y=y′-ζy(z,t)

(3)

式中:ζx(z,t)和ζy(z,t)分別表示結構在順流向和橫流向的位移.相應地,速度項和壓力項采用以下變換:

(4)

因此,將式(3)和(4)代入N-S方程式(1)和(2)可得:

A(u,p,ζ)

(5)

(6)

(7)

在此定義:

(8)

此外,本文采用線性張力梁模型描述細長柔性結構的動力學行為,該模型采用小變形假設,可寫成:

(9)

式中:ρc,k和T分別表示結構的單位長度的質量、阻尼比和張力.需要指出的是,張力T的大小將影響相速度c,c=(T/ρc)1/2;ζ(z,t)=(ζx,ζy)表示結構在順流向和橫流向的位移;F(z,t)是作用于結構上的流體力,通過對壓力和黏性力項沿結構固壁表面積分獲得,

(10)

n為指向結構外向的法線單位向量;s為結構表面微分.此外,假設結構的動力響應滿足沿展向的周期性條件,則對流場和結構變量可采用傅里葉級數(Fourier Expansion)表示,即

(11)

(12)

式中:β=2π/Lz表示展向的波數,Lz為結構長度;M為展開式中傅里葉模態個數;m為各階模態.將式(11)和(12)代入式(5)和(6)可獲得解耦后的二維模態方程:

(13)

(14)

(15)

對結構運動的解耦方程式(15),采用二階Newmark-β方法進行求解.對不定常流場模態方程式(13)和(14),采用Karniadakis等[9]提出的高階分步格式進行時間離散.具體地,對每一時間迭代步對速度和壓力進行解耦計算.

(16)

式中:αq和βq均為與強穩定化積分(Stiffly Stable Integration)有關的參數.

第2步,將考慮壓力梯度的作用修正速度場,并施加連續性約束條件和紐曼邊界條件:

(17)

第3步,考慮黏性項更新下一步的速度場:

(18)

式中:γ0為強穩定計算過程中向后差分系數(Backwards Differentiation Coefficient).

對上述時間離散后的模態方程,根據Karniadakis等[9-10]將(x,y)平面的二維空間域離散為四邊形有限元單元,并采用高斯-洛巴托-勒讓德(GLL)高階正交拉格朗日多項式作為形函數進行空間離散.

2 數值模型

2.1 物理模型

考慮在均勻來流作用下的波浪形變截面圓柱(見圖1),其直徑沿展向變化由下式確定:

圖1 均勻流下截面直徑沿展向余弦變化的波浪柱幾何示意圖

(19)

式中:Dz為展向相應位置的圓截面直徑.由于波浪形柱直徑在展向呈余弦變化,平均直徑為D=(Dmax+Dmin)/2.基于平均直徑和均勻來流速度的雷諾數取值為100,因此本文考慮的流動為層流流動.定義波浪形柱直徑最大截面所在的展向位置為幾何節點,波浪柱直徑最小截面所在的對應位置為幾何鞍點.A代表波狀表面的波高,取值范圍為0.1D～0.3D;展向波長設置為λ=4πD.A=0對應于光滑圓柱,也作為基準工況用于與波浪形柱工況進行對比分析.參照 Newman等[11]的研究,可以通過規定圓柱的初始振幅和速度來確定其初始條件,本研究中初始擾動為駐波:ζy(z,t)=acos(ωt)cos(2πz/Lz),其中a為振幅,ω為振動頻率,ω=2π/(cLz).

二維平面網格劃分如圖2所示.計算域幾何形狀為C形,由圓柱上游半徑為20D的半圓弧和下游矩形組成,其流向長度為30D,在橫流向長度設置為40D,圓柱放置在計算域橫流向的中心位置.在(x,y)平面內采用了877個基礎網格;由于對每一基礎單元再用6階正交多項式(P=6)對速度和壓力場進行近似,所以每個基礎網格再細分為5×5個高階網格.此外,根據文獻[11-12],沿展向采用64個傅里葉模態,滿足本文雷諾數下的精度要求.

圖2 二維網格拓撲結構

2.2 模型驗證

為了驗證力學模型和數值模擬方法的適用性,首先獲得了均勻來流作用下細長光滑圓柱渦激振動模擬結果,并與已有文獻結果進行對比分析.本模擬中質量比設置為m*=6.0,而其他模擬參數均與Newman等[11]相同.圖3所示為光滑圓柱在駐波和行波初始激勵下的橫向振幅(ξy)、阻力系數(Cd)和升力系數(Cl)模擬結果,與Newman等[11]結果吻合良好.尤其在駐波初始激勵下本模擬結果與Newman等[11]結果在展向節點(響應振幅最大的點)和反節點(響應振幅最小的點)位置坐標完全相同.在行波初始激勵下的行波響應幅值和行波速度等均與Newman等[11]結果一致,故說明本文力學模型和模擬方法的可靠性.

圖3 柔性光滑圓柱渦激振動擬結果

接著正交多項式階數對計算結果的影響進行驗證分析.通過變化多項式階數P改變高階網格的疏密;在展向上通過改變傅里葉模態個數,改變展向網格的疏密.由計算結果表明:P=6,M=64時模擬結果足夠收斂至精確解.

3 結果與討論

3.1 運動響應與水動力系數

圖4 ξy、Cd和Cl沿展向的均方根值分布

圖5為在A=0.10D,0.20D時的結構橫向振動響應時程模擬結果.由圖可知,兩種工況下節點和反節點的展向位置與光滑截面結構相同,即節點位于z/D=0, 2π, 4π,反節點位于z/D=π, 3π處,表明該工況下結構振動模式均與光滑表面結構相同.阻力系數在A=0.10D時其隨時間演化特征與光滑表面結構基本相同,均在節點處取得最大值,在反節點處取最小值.然而A=0.20D工況下的阻力系數時空分布則與光滑表面結構完全不同,在阻力系數分布中看不到明顯的駐波效應.在z/D=0, 4π時阻力系數并未取得最大值,且最大阻力系數僅有光滑表面結構的50%,表明A=0.20D擾動波高具有良好的減阻效果.進一步觀察升力系數時程可知,A=0.10D工況下波形柱在展向π

圖5 結構響應與水動力系數計算結果

3.2 尾流特性

為了進一步了解擾動表面波高對結構振動響應的影響機理,需要對柔性波形柱和柔性光滑柱近尾跡三維渦結構進行對比分析.圖6給出了不同擾動波高時展向渦量等值面(ωz=±1.0)計算結果.其中,圖6(a)～6(c)為柔性光滑柱展向渦量等值面透視圖,圖6(d)～6(f)為俯視圖,與Newman等[11]計算結果吻合良好.駐波初始激勵下柔性光滑柱展向渦量中存在明顯的交織結構.這種交織渦結構的展向特征是,在z/D=0, 2π截面處形成交錯脫落的卡門渦街結構;而在z/D=π, 3π截面處上、下表面同時脫落的完全對稱型渦結構.圖6中分別給出A=0.10D時模擬結果,顯然在z/D=π, 3π截面處的渦結構也呈現對稱分布.值得注意的是,在此工況下幾何節點處的展向渦量幅值較大從而升力系數幅值也隨之增大,而幾何鞍點處渦量幅值較小從而升力系數幅值也隨之減小,這與圖5所示結果一致.圖6中分別給出A=0.20D時模擬結果,在該工況下交織渦結構完全消失,彎曲渦管以交錯方式從波浪形表面上、下側分別脫離,并且在尾流中迅速消失.更重要的是,與柔性光滑柱相比,此時波形柱表面剪切層的卷曲和相互作用明顯較弱,使得渦形成區長度進一步增大.Lin 等[13]用數值方法研究了亞臨界雷諾數條件下具有相對較大展向波長的剛性波形柱周圍的流動特征,也得到了類似的結果.

圖6 展向渦量瞬時等值面(ωz=±1.0)的透視圖和橫流向俯視圖

圖7所示為不同擾動波高條件下的渦量場順流向(ωx=±1.0)與橫流向(ωy=±1.0)分量瞬時等值面云圖.從該圖可清晰分辨出順流向渦結構和橫流向渦結構在尾跡附近的演化特征.顯然,與柔性光滑柱和A=0.10D擾動波高工況相比,A=0.20D時結構附近渦結構完全不同且渦量強度相對較低,并且沿下游耗散較快.Lin等[13]指出,附加流向渦在波形柱表面幾何鞍點和節點處分別誘導出上涌流和下涌流,使得在幾何鞍點處產生寬尾流,而在幾何節點處產生窄尾流.這與本文模擬結果非常一致,表明波形結構的渦激振動抑制和減阻效應的力學機理有相同之處.

圖7 3種工況下的順流向(ωx=±1.0)與橫流向(ωy=±1.0)瞬時等值面

為了進一步研究流向渦對結構響應的影響,圖8中給出了3種工況下(x,y)平面內流向渦量切片圖;其中圖7中分別為柔性光滑柱和在A=0.10D,0.20D時柔性波形柱的模擬結果.該結果表明,柔性光滑柱兩側和A=0.10D波形柱兩側均形成一對同向旋轉渦結構,而沿結構展向上、下側分布的兩對渦旋轉方向剛好相反.然而,在A=0.20D情況下,在同一展向位置處柔性波形柱兩側生成一對較強的反向旋轉渦,在這對反向旋轉渦的外側還分布著一對與各自強渦旋轉方向相反的弱渦.Lin 等[13]指出,最優控制波形柱附加產生的反向流向渦對剪切層穩定性具有重要作用,能夠防止剪切層與強渦結構的相互作用,使得渦旋形成長度增大.在A=0.20D的情況下觀察到的反向旋轉渦旋也穩定了剪切層,拉長了渦旋形成區長度,與本文之前觀察到的現象一致.

圖8 3種工況下(x, y)平面內圓柱周圍的瞬時等值面和相應的流向渦量示意圖

3.3 能量轉換特征與頻譜特性

結構的運動響應與流體與結構之間的能量轉化特征密切相關,根據 Newman 等[11]的研究,一個旋渦脫落周期的無量綱時均能量E(z)可定義為

(20)

式中:Γ代表無量綱脫落周期;Wl代表升力產生能量的功率,即Wl=?ζy/?t.當E(z)為正值時,代表能量從流體轉移至結構,E(z)為負值則轉移方向相反.Newman等[11]和Zhu等[12]研究指出對于駐波響應,E(z)取值與展向位置有關,而與無量綱時間無關.同時,E(z)在一個脫落周期內沿整數波長的積分應該為0,表明在一個旋渦脫落周期內,結構從流體中得到的能量和在流體中做功耗散的能量相同,在沒有外界能量輸入的前提下,整個流體和結構系統滿足能量守恒.因此,時均能量代表結構與流體進行能量交換的強度,其值越大則相應結構與流體之間的能量交換越頻繁,反映到實際工程中,就會引起結構的疲勞,加速結構的老化.

圖9所示為所有工況下時均能量隨展向分布計算結果.由圖9可知,A=0.10D柔性波形柱的時均能量比柔性光滑柱顯著增大,當擾動波面高度增至A=0.15D其時均能量相比于柔性光滑柱維持在同一水準.而當A≥0.20D時,E(z)始終在0附近小幅振蕩,說明在上述3種工況下波形柱與流體之間的能量交換接近于0,從能量角度驗證了波形柱與流體之間的動力響應被顯著抑制.同時,在z/D=0.37π～1.62π以及z/D=2.37π～3.62π之間柔性光滑柱E(z)為正值,在其余位置E(z)為負值;而A=0.10D波浪形結構在兩端E(z)為正值,而在跨中附近E(z)為負值,這與圖5中升力系數的時空演化模擬結果相互一致.

圖9 細長柔性結構渦激振動時的無量綱時均能量分布

完全消除初始瞬態效應后,對3種典型工況下橫向位移和升力系數時程進行了頻域分析.

圖10中分別為柔性光滑柱橫向位移和升力系數功率譜密度(PSD)展向分布計算結果,可以據此確定旋渦脫落頻率為fD/U∞=0.163,此外,在PSD的展向分布云圖中可觀察到,譜峰幅值與圓柱展向波長密切相關.Bourguet等[14]研究了長徑比為200的細長柔性結構渦激振動問題,并也報道了類似現象.由此可以進一步證實,結構的振動響應與波長相關.此外,阻力頻率(這里沒有顯示)是升力頻率的兩倍,說明尾跡是經典的2S模式.

圖10 光滑圓柱無量綱時均能量分布

圖11中分別為A=0.10D柔性波浪柱橫向位移和升力系數PSD展向分布模擬結果,可以得知其主頻為fD/U∞=0.161,略小于光滑柱時相應頻率.從圖11(b)可以看到,與光滑柱相比,A=0.10D柔性波形柱振動時程在跨中對應的PSD值明顯減小,表明該位置附近升力的波動強度顯著減小,與圖5顯示的現象一致.

圖11 A=0.1D波浪柱無量綱時均能量分布

圖12中分別為A=0.20D柔性波形柱橫向位移和升力系數PSD的展向分布計算結果.對于A=0.20D柔性波形柱,除了主頻fD/U∞=0.161外,在fD/U∞=0.126處也可識別出次峰頻率,在圖12用白色虛線表示.這種低頻次生頻率,在PSD-Cl中比在PSD-ζy中更強.次生頻率產生的原因涉及流固耦合的非線性問題,其機理有待進一步探討.

圖12 A=0.2D波浪柱無量綱時均能量分布

4 結論

采用基于高精度譜單元方法的流固耦合分析方法,建立了細長柔性波形柱體結構的渦激振動力學模型,對其在低速均勻流作用和駐波初始動力激勵下的流致振動機理進行了量化分析.結果顯示,當擾動波高A≥0.20D時,柔性波形柱橫向振幅相比于光滑柱明顯減小,升阻力系數顯著降低,水動力特性得到明顯改善,表明在合適波高下,柔性波形柱具有良好的流致振動抑制作用.進一步比較擾動波高分別為A=0.10D,0.20D的結構橫向振幅和水動力系數的時空演變特征,大致確定出波形柱抑制渦激振動的控制失效范圍(A<0.10D)以及控制有效范圍(0.20D≤A<0.30D).

為了探究擾動表面波高對結構振動響應的影響機理,對擾動表面結構和無擾動表面結構近尾跡三維渦結構進行了對比分析.結果表明,與柔性光滑柱相比,A=0.20D工況下交織渦結構完全消失,使得波形擾動截面剪切層的卷曲和相互作用明顯較弱,渦形成區長度進一步增大.對柔性光滑柱在A=0.10D,0.20D時波形柱的流向渦結構瞬時等值面進行了模擬比較,結果表明在A=0.20D的情況下觀察到的反向旋轉渦旋也穩定了剪切層,拉長了渦旋形成區長度.比較了所有工況下時均能量的展向分布計算結果,發現當A≥0.20D時,E(z)始終在0附近小幅振蕩,從能量轉移角度進一步驗證柔性波形柱對渦激振動的抑制作用.在完全消除初始瞬態效應后,對3種典型工況下橫向位移和升力系數時程進行了頻域分析,確定了3種工況下的漩渦脫落頻率,對于A=0.20D波形柱工況下觀察到的次生頻率,其產生的原因涉及流固耦合的非線性問題,有待對其機理進行進一步探討.