(天津大學(xué) a.水利工程仿真與安全國家重點實驗室；b.建筑工程學(xué)院， 天津 300072)

0 引 言

渦激振動(Vortex-Induced Vibration， VIV)現(xiàn)象普遍存在于橋梁工程、土木工程和海洋工程等諸多領(lǐng)域，是橋梁斜拉索、工廠煙囪、海洋立管等大長徑比(長度與直徑的比值，L/D)圓柱結(jié)構(gòu)疲勞損傷的重要外界影響因素。近一個世紀(jì)以來，VIV一直是備受關(guān)注的核心科學(xué)問題之一，并取得了大量卓越的研究成果[1-4]。目前，人們對此復(fù)雜流-固耦合現(xiàn)象的理解已經(jīng)上升到了一個全新的高度。

彈性支撐剛性圓柱VIV的流體力特性研究開展較為廣泛，常將垂直于來流方向(即橫流向)以及平行于來流方向(即順流向)的流體力合力分解為與速度同相位的力——升力和脈動阻力以及與加速度同相位的力——附加質(zhì)量力。SARPKAYA[5]試驗表明在圓柱橫流向VIV“鎖定”現(xiàn)象發(fā)生時，響應(yīng)位移與橫流向流體力合力的相位角存在跳躍，該現(xiàn)象也被其他學(xué)者[6-7]通過數(shù)值仿真證實。GOPALKRISHNAN[8]繪制升力系數(shù)和附加質(zhì)量系數(shù)云圖，發(fā)現(xiàn)升力系數(shù)與響應(yīng)頻率和響應(yīng)位移密切相關(guān)。VIKESTAD等[9]通過能量轉(zhuǎn)化法分解橫流向流體力，獲得了附加質(zhì)量系數(shù)，并發(fā)現(xiàn)附加質(zhì)量系數(shù)與響應(yīng)頻率相關(guān)性強，隨約化速度增大而減小，在“鎖定”區(qū)域變化緩和。ARONSEN[10]研究圓柱順流向振動的流體力，并運用能量轉(zhuǎn)化法分解得到脈動阻力系數(shù)和順流向附加質(zhì)量系數(shù)。

與彈性支撐剛性圓柱相比，柔性圓柱VIV更復(fù)雜：結(jié)構(gòu)振動表現(xiàn)為多模態(tài)，模態(tài)競爭激烈，在響應(yīng)中出現(xiàn)高頻成分，順流向振動頻率一般近似為橫流向的2倍[11-13]。柔性圓柱VIV過程結(jié)構(gòu)振動變形，不同截面處受到的流體力不同。由于測量手段和試驗條件的限制，難以在不改變結(jié)構(gòu)外部流場的前提下測量柔性圓柱VIV的流體力分布。SANAATI等[14]測量柔性圓柱VIV總體流體力，研究表明：橫流向和順流向流體力合力在“鎖定”區(qū)域先增大后減小，在最大位移對應(yīng)的約化速度處達(dá)到最大值。SONI[15]開展模型試驗獲得柔性圓柱VIV截面處的運動軌跡，并迫使剛性圓柱以上述運動軌跡振動，假定剛性圓柱的水動力系數(shù)等價于柔性圓柱截面處的水動力系數(shù)，試驗發(fā)現(xiàn)：順流向脈動阻力系數(shù)和附加質(zhì)量系數(shù)受振動位移的影響較小，橫流向附加質(zhì)量系數(shù)出現(xiàn)負(fù)值，升力和脈動阻力出現(xiàn)倍頻成分。YIN等[16]采用與SONI[15]相似方法識別柔性圓柱VIV的水動力系數(shù)，結(jié)果表明：若柔性圓柱的振動位移呈現(xiàn)準(zhǔn)周期特性，則計算誤差較??；若響應(yīng)位移具有一定隨機性，則此方法的不確定性增大。

現(xiàn)階段柔性圓柱結(jié)構(gòu)VIV的預(yù)報主要依靠經(jīng)驗?zāi)Ｐ停A(yù)報結(jié)果與模型試驗結(jié)果存在一定差異[17]，主要原因是經(jīng)驗?zāi)Ｐ瓦x取的流體力系數(shù)來源于剛性圓柱試驗。如何準(zhǔn)確地識別柔性圓柱VIV流體力一直是工程領(lǐng)域和學(xué)術(shù)領(lǐng)域共同關(guān)注的問題。本文系統(tǒng)地總結(jié)現(xiàn)存VIV流體力計算和分解方法，闡述不同方法的計算過程，分析誤差來源和適用范圍，并對柔性圓柱VIV流體力系數(shù)識別方法的進(jìn)一步完善提出若干合理化建議。

1 VIV流體力計算方法

目前，廣泛采用根據(jù)結(jié)構(gòu)振動信息逆向求解流體力的方式確定柔性圓柱VIV流體力。工程中的大長徑比柔性圓柱結(jié)構(gòu)一般承受軸向預(yù)張力，其結(jié)構(gòu)振動控制方程可采用Euler-Bernoulli梁模型，公式可表示為

(1)

(2)

1.1 一般有限元法

柔性圓柱結(jié)構(gòu)可離散為若干個單元，每個單元包含2個節(jié)點，共計4個自由度。結(jié)合圓柱結(jié)構(gòu)兩端邊界條件，柔性圓柱振動的有限元方程形式為

(3)

C=αM+β(KE+KG)

(4)

式中：α和β為常系數(shù)，通過柔性圓柱在空氣中的自由衰減試驗測得。根據(jù)獲取的流體載荷矩陣Fy，確定柔性圓柱各節(jié)點處流體力合力fy。HUERA-HUARTE等[18]采用有限元法計算階梯來流下柔性圓柱VIV的流體力，并計算相應(yīng)的流體力系數(shù)。有限元法求解的阻力系數(shù)與測量值對比如圖1所示，計算結(jié)果與試驗結(jié)果吻合較好。圖1中橫坐標(biāo)為約化速度Vr(Vr=U/f1D, 其中U為來流速度，D為結(jié)構(gòu)外徑，f1為結(jié)構(gòu)靜水中一階固有頻率)，縱坐標(biāo)為平均阻力系數(shù)Cd(Cd=2Fd/(ρDLU2), 其中Fd為結(jié)構(gòu)順流向的平均阻力，ρ為外部流體密度，L為結(jié)構(gòu)長度)。WU[19]采用有限元法計算流體力并與VIV預(yù)報軟件VIVANA的結(jié)果進(jìn)行對比，發(fā)現(xiàn)兩種方法均可取得精度較好的計算結(jié)果。有限元法發(fā)展較成熟且穩(wěn)定可靠，是目前廣泛應(yīng)用的方法。除均質(zhì)性等材料本身性質(zhì)所造成的誤差之外，誤差的主要來源是計算過程中剛度、質(zhì)量、阻尼矩陣不會隨著時間變化而更新[18]，此外還受計算的空間分辨率即網(wǎng)格劃分的精細(xì)程度的影響。

圖1 有限元法求解的阻力系數(shù)與測量值對比

1.2 凝聚自由度有限元法

此方法基于有限元計算模型，將一般有限元方程中結(jié)構(gòu)轉(zhuǎn)角自由度凝聚。對于具有n個平動自由度和s個轉(zhuǎn)動自由度的系統(tǒng)，可通過矩陣H*轉(zhuǎn)換為僅有n個平動自由度的系統(tǒng)。

唐國強[20]闡述凝聚自由度的有限元法，識別立管VIV的流體力系數(shù)。LI等[21]采用此方法識別海底懸跨管道VIV流體力系數(shù)。通過模態(tài)分解法可將由柔性圓柱VIV試驗獲得的應(yīng)變信息轉(zhuǎn)化為位移，但轉(zhuǎn)角未計算求得。在凝聚自由度有限元法應(yīng)用過程中，可省去計算轉(zhuǎn)角的步驟，但凝聚自由度的過程存在一定誤差，當(dāng)結(jié)構(gòu)單元質(zhì)量較大即結(jié)構(gòu)慣性力較大時，此法的精度較差。當(dāng)結(jié)構(gòu)振動頻率較高時，凝聚自由度法也會引入較大的誤差。此方法僅在結(jié)構(gòu)單元質(zhì)量較小、振動頻率較低時才更合理適用[22-23]。

1.3 基于廣義積分變換的反解法(GITT)

結(jié)合兩端邊界條件，根據(jù)廣義積分變換法，柔性圓柱結(jié)構(gòu)橫流向振動頻率即特征值求解方程為

(5)

式中：λi為特征值；φi(z)為振型函數(shù)，i為第i階模態(tài)階次。根據(jù)特征方程的正交準(zhǔn)則和式(6)

(6)

(7)

(8)

GU等[24]采用廣義積分法計算柔性圓柱VIV流體力，計算結(jié)果與HUERA-HURATE等[18]的有限元法結(jié)果對比如圖2所示，結(jié)果吻合較好，進(jìn)一步表明廣義積分法對柔性圓柱VIV流體力的識別具有可行性。在圖2中，T1為HUERA-HURATE等的計算數(shù)據(jù)，T2為GU等的計算數(shù)據(jù)。該方法應(yīng)用還不甚廣泛，對其適用條件和誤差來源等有待進(jìn)一步研究[24]。

圖2 GITT法與有限元法求解的順流向阻力系數(shù)對比

1.4 基于卡爾曼濾波的反解法

式(3)采用狀態(tài)空間方程表示為

(9)

式中：W為系統(tǒng)狀態(tài)向量；Λ為系統(tǒng)輸出向量；A、B、G分別為具有一定維數(shù)的系統(tǒng)矩陣。

考慮測量噪聲v(t)和過程噪聲w(t)的影響，離散狀態(tài)空間平衡方程為

(10)

其中

(11)

式中：k為離散的時間步長；τ為時間。測量噪聲v(t)和過程噪聲w(t)的方差用R和Q表示。R和Q需在計算結(jié)構(gòu)VIV流體力前計算獲得，在試驗中測量噪聲方差R較易確定，但過程噪聲方差Q的計算較困難，可通過數(shù)值模擬確定。具體的流體力計算過程見文獻(xiàn)[19]。

在流體力計算中引入遺忘因子。遺忘因子的作用是加強當(dāng)前數(shù)據(jù)的影響并削弱歷史數(shù)據(jù)的影響，消除數(shù)據(jù)飽和現(xiàn)象。若遺忘因子等于1.0，則為基本最小二乘法，此時歷史數(shù)據(jù)的影響較強，噪聲的影響減弱。若遺忘因子較小，則為遞推最小二乘法，當(dāng)前數(shù)據(jù)的影響較強，計算收斂的速度較快。需權(quán)衡計算精度與計算效率選擇合理的遺忘因子。

上述方法基于卡爾曼濾波，核心思想是用前一步的歷史信息計算當(dāng)前時刻的流體力。MA等[25]最早將基于卡爾曼濾波的反解法用于懸臂梁結(jié)構(gòu)的輸入力識別，并通過與其他數(shù)值模擬法的對比驗證該方法的準(zhǔn)確性。WU[19]亦運用該方法識別柔性圓柱VIV流體力。在采用此法計算結(jié)構(gòu)流體力時，流體力與結(jié)構(gòu)響應(yīng)位移、響應(yīng)加速度之間存在不確定的時間延遲[19,25]。在流體力系數(shù)計算過程中需考慮流體力與響應(yīng)位移之間的相位差。結(jié)構(gòu)流體力計算時間延遲會在確定流體力系數(shù)時引入一定的誤差[19]。此外過程噪聲方差Q和測量噪聲方差R較大也會使測量結(jié)果產(chǎn)生誤差，但先進(jìn)測量方法已可將其控制在較低水平，在一定程度上保證反解的準(zhǔn)確性[25]。

1.5 基于最優(yōu)控制理論的反解法

基于最優(yōu)控制理論的反解法主要通過使?fàn)顟B(tài)向量W與輸出向量Λ無限接近，從而計算出結(jié)構(gòu)的流體力。此法本質(zhì)上是一種優(yōu)化算法，即得到誤差最小的最優(yōu)解。

假設(shè)向量Λ服從期望值為Λ0、標(biāo)準(zhǔn)差為σΛ的高斯分布。目標(biāo)函數(shù)可表示為

(12)

式中：QΛΛ為輸出向量的協(xié)方差矩陣；σΛ為標(biāo)準(zhǔn)差；I為單位矩陣。結(jié)構(gòu)流體力的目標(biāo)函數(shù)為

(13)

式中：QFyFy為結(jié)構(gòu)流體力的協(xié)方差矩陣；σFy為計算得到的流體力與真實流體力之間誤差的標(biāo)準(zhǔn)差。結(jié)合兩個目標(biāo)函數(shù)，基于最優(yōu)控制理論即可計算流體力[19]。

上述方法的計算結(jié)果受QΛΛ和QFyFy的影響較大。若QΛΛ較大，計算出的響應(yīng)位移與測量結(jié)果相關(guān)性強，計算結(jié)果受測量誤差的影響較大，造成結(jié)果的方差較大。若QFyFy較大，計算結(jié)果與測量結(jié)果的相關(guān)性弱，造成結(jié)果的偏差較大。若測量的數(shù)據(jù)準(zhǔn)確性較高，QΛΛ和QFyFy將處于相對穩(wěn)定的范圍，計算結(jié)果的方差和偏差均較小[19]。

WU[19]將該方法識別的柔性圓柱VIV流體力的結(jié)果與利用VIVANA模型和有限元法計算結(jié)果進(jìn)行對比，在低模態(tài)和高模態(tài)取較多測點的情況下，結(jié)果的相似性較好。如果在高模態(tài)下測點較少，結(jié)果會產(chǎn)生一定的誤差。該方法的主要難點在于數(shù)值計算，自身屬性不好的矩陣易在計算過程中出現(xiàn)不穩(wěn)定的問題。此外，上述方法的計算時間與自由度的四次方成正比，若自由度較大，應(yīng)用上述方法較為困難。在高模態(tài)下的計算準(zhǔn)確度仍有待驗證[26]。采用凝聚自由度的方法將柔性圓柱結(jié)構(gòu)的轉(zhuǎn)動自由度凝聚，能顯著改善數(shù)值計算條件[19]。

1.6 基于振動方程的直接反解法

SONG等[27]采用上述方法計算柔性圓柱結(jié)構(gòu)的流體力，并通過ABAQUS有限元計算軟件檢驗上述方法的計算精度。首先對柔性圓柱結(jié)構(gòu)施加已知載荷，運用ABAQUS求解結(jié)構(gòu)的振動響應(yīng)，根據(jù)振動響應(yīng)采用直接反解法計算柔性圓柱結(jié)構(gòu)受到的載荷，比較計算獲得載荷與施加載荷的差異。該方法的準(zhǔn)確程度直接取決于測點位移的精確程度，其計算精度亦受試驗測點個數(shù)的影響。若測點個數(shù)大于等于柔性圓柱結(jié)構(gòu)實際激發(fā)的模態(tài)數(shù)，直接反解法的計算精度較高；若測點個數(shù)小于柔性圓柱結(jié)構(gòu)實際激發(fā)的模態(tài)數(shù)，由直接反解法計算得到的結(jié)構(gòu)流體力與結(jié)構(gòu)實際受力差異較大，計算精度較低。

2 流體力分解方法

以升力系數(shù)和橫流向附加質(zhì)量系數(shù)的計算方法為例介紹流體力分解的方法，脈動阻力系數(shù)和順流向附加質(zhì)量系數(shù)的計算方法與之類似，不再贅述。在柔性圓柱VIV過程中結(jié)構(gòu)實際受到的流體力可表示為三角函數(shù)形式，并分解為與速度同相位的力和與加速度同相位的力，公式為

fy=f0sin(ωt+φ)=f0cosφsin(ωt)+f0sinφcos(ωt)

(14)

式中：ω為橫流向的振動圓頻率；φ為橫流向流體力與振動位移相位差；f0為橫流向流體力幅值。

2.1 最小二乘法

柔性圓柱VIV橫流向流體力合力的真實值可用升力系數(shù)CL和橫流向附加質(zhì)量系數(shù)Cay表示為

(15)

(16)

通過最小二乘擬合使計算值與真實值的差最小，得到升力系數(shù)和橫流向附加質(zhì)量系數(shù)的計算公式為

(17)

其中

ARONSEN[10]介紹最小二乘法分解剛性圓柱的流體力合力。SONG等[27]根據(jù)直接反解法確定結(jié)構(gòu)流體力，采用最小二乘法獲得柔性圓柱VIV的升力系數(shù)、脈動阻力系數(shù)和附加質(zhì)量系數(shù)。最小二乘法是一種偏數(shù)值方法，在流體力分解過程中不涉及數(shù)值積分，通過Matlab等數(shù)學(xué)軟件采用最小二乘擬合的方式，計算得到與真實水動力系數(shù)誤差最小的升力系數(shù)和附加質(zhì)量系數(shù)。

2.2 能量轉(zhuǎn)化法

在柔性圓柱VIV過程中轉(zhuǎn)移能量的平均值可表示為

(18)

振動位移可表示y(t)=y0sin(ωyt)，能量轉(zhuǎn)移的平均值可表示為

(19)

由式(18)和式(19)可得

(20)

通過相同的推導(dǎo)方式計算與加速度同相位的項為

(21)

升力系數(shù)、橫流向附加質(zhì)量系數(shù)計算式為

(22)

能量轉(zhuǎn)化法是剛性圓柱流體力分解應(yīng)用較多的方法。ARONSEN[10]用該法分解剛性圓柱純順流向流體力。VIKESTAD等[9]用能量轉(zhuǎn)化法獲得剛性圓柱結(jié)構(gòu)的附加質(zhì)量系數(shù)，并與GOPALKRISHNAN[8]的結(jié)果進(jìn)行對比，結(jié)果吻合良好。TANG等[28]和WU[19]根據(jù)能量轉(zhuǎn)化法分解柔性圓柱流體力合力，與模型計算對比，結(jié)果亦比較準(zhǔn)確。

2.3 傅里葉平均法

傅里葉平均法是根據(jù)傅里葉變換的基本理論對流體力進(jìn)行分解的方法。將計算得到的流體力fy乘以cos(ωt)，在n個振動周期(TOSC)內(nèi)積分，可得與速度同相位的力為

(23)

將計算得到的流體力fy乘以sin(ωt)，在n個振動周期內(nèi)積分，可得與加速度同相位的力為

(24)

GOPALKRISHNAN[8]采用傅里葉平均法分解剛性圓柱橫流向流體力合力得到升力系數(shù)和附加質(zhì)量系數(shù)。傅里葉平均法在流體力分解過程中僅使用流體力信息，誤差主要來源于式(23)和式(24)的數(shù)值積分過程。

2.4 傳遞函數(shù)法

傳遞函數(shù)法為頻域方法，假定輸入函數(shù)(結(jié)構(gòu)VIV的頻率和位移)與輸出函數(shù)(結(jié)構(gòu)VIV的流體力)之間為線性關(guān)系。振動位移與流體力的相關(guān)函數(shù)H(ω)應(yīng)為位移、流體力的互功譜S(ω)與位移自功譜Syy(ω)的比值，即

(25)

當(dāng)已知振動位移的幅值y0時，f0和φ的計算式為

(26)

然后求出與速度同相位的力f0sinφ和與加速度同相位的力f0cosφ，ARONSEN[10]提出傳遞函數(shù)法用于剛性圓柱流體力合力分解。傳遞函數(shù)法的誤差相對較大，傳遞函數(shù)法在頻域內(nèi)進(jìn)行計算，誤差主要來源于頻率分辨率，計算的結(jié)果相對比較保守。

3 結(jié)論與展望

根據(jù)結(jié)構(gòu)響應(yīng)逆向求解結(jié)構(gòu)受力是識別柔性圓柱VIV流體力的常用方法。分析對比多種流體力計算方法，總結(jié)歸納了不同方法的適用范圍和誤差來源。研究發(fā)現(xiàn)：有限元法是較為傳統(tǒng)、成熟的一般方法，在柔性圓柱結(jié)構(gòu)VIV流體力計算方面應(yīng)用廣泛，即使因采用的矩陣在計算過程中不隨時間步長更新而產(chǎn)生一定誤差，總體計算結(jié)果仍相對穩(wěn)定[18]。凝聚自由度法在計算結(jié)構(gòu)流體力過程中雖可降低結(jié)構(gòu)自由度數(shù)目，提高計算效率，簡化計算步驟，但對于單元質(zhì)量較大、響應(yīng)頻率較高的結(jié)構(gòu)，凝聚自由度法將引入較大誤差[22-23]?；诳柭鼮V波和最小二乘法的反解法在計算流體力時，采用前一時刻響應(yīng)位移計算下一時刻流體力，存在不確定的時間延遲，在后續(xù)流體力分解過程中將引入不確定的誤差[19]?；谧顑?yōu)控制理論的反解法在數(shù)值計算過程中存在難點，對數(shù)據(jù)要求較高，質(zhì)量較差的數(shù)據(jù)易引起計算不穩(wěn)定的問題[19]；計算周期較長，計算時間與自由度數(shù)目的四次方成正比，受自由度數(shù)目的限制，雖然能通過凝聚自由度來提高計算效率，但凝聚自由度方法的應(yīng)用會引入新的誤差。直接反解法和廣義積分變換法在計算柔性結(jié)構(gòu)VIV流體力方面具有可行性，但應(yīng)用并不廣泛，其穩(wěn)定性需進(jìn)一步檢驗[24]。

為了深入研究柔性圓柱VIV流體力特性，通常將流體力合力分解為與速度同相位分量和與加速度同相位分量。流體力的分解方法亦有多種：頻域上的傳遞函數(shù)法誤差相對較大[10]，最小二乘法和能量轉(zhuǎn)移法應(yīng)用較為廣泛，其精確度與計算過程的積分計算關(guān)系較大[9-10,27]。應(yīng)用能量轉(zhuǎn)移法、傅里葉平均法以及最小二乘法的基本前提是VIV位移和流體力的時域信號應(yīng)近似為簡諧形式，時域信號越接近簡諧形式，誤差越小。

為了提高識別精度和結(jié)果的可靠性，未來關(guān)于圓柱VIV流體力系數(shù)識別的研究可著眼于以下幾方面：

(1) 目前的流體力計算方法僅考慮單個頻率，柔性圓柱VIV多個模態(tài)共存，響應(yīng)頻率存在倍頻成分，其影響有待深入研究。應(yīng)建立包含多頻成分的流體力計算模型，更全面獲取圓柱VIV流體力信息。

(2) 現(xiàn)有方法大多將模型試驗的應(yīng)變信息轉(zhuǎn)化為位移，通過位移求解結(jié)構(gòu)的受力，這一過程存在不可避免的誤差。為提高計算精度，可改進(jìn)計算方法，直接通過結(jié)構(gòu)應(yīng)變信息識別結(jié)構(gòu)受力。

(3) 改進(jìn)目前試驗測量手段，通過傳感器直接測量得到柔性圓柱VIV流體力分布，將大幅降低流體力識別誤差。