利用延遲反饋控制柔性圓柱渦激振動數值研究

鄒 琳，王 程，徐勁力，陶 凡，左紅成

(武漢理工大學 機電工程學院，武漢 430070)

鈍體周圍非對稱的漩渦脫落在結構上產生周期性的脈動壓力，從而引發渦激振動(vortex induced vibration，VIV)現象[1-3]。當漩渦脫落頻率接近結構固有頻率時會發生大幅度的振動，導致結構疲勞損壞，也可以利用渦激振動產生新能源[4-5]。對于類似于海洋立管的柔性圓柱而言，由于展向具有柔性張力，渦激振動通常表現為更為復雜的駐波響應特性、行波響應特性和多模態振動特性等。

國內外學者對柔性圓柱渦激振動展開了大量數值及試驗研究。早期數值研究主要是利用直接數值模擬方法求解N-S方程，以獲得柔性圓柱的渦激振動響應，如Newman等[6]通過數值模擬方法研究了Re=100和Re=200下柔性圓柱的橫向流振動，發現振動模式為行波與駐波混合響應模式；Lucor等[7]通過直接數值模擬研究了長徑比超過500的細長圓柱在線型和指數剪切流流作用下的振動響應，研究認為指數剪切流速剖面會產生更大的振動頻率范圍，但和實際情況相比，圓柱單自由度振動具有一定的局限性。

Lie等[8]開展了剪切流作用下物理模型試驗，研究立管長細比高達3 000，分析發現立管的振動響應是不規則的，具有一定的帶寬，并且強度也隨著帶寬的不同有所變化。國內學者高云等[9]通過試驗研究，深入研究了柔性圓柱體在剪切來流下渦激振動響應，發現當折合速度較小時，響應頻率會全程參與渦激振動，隨著折合流速的上升，響應頻率則會間歇性參與渦激振動。宋磊建等[10]采用模型試驗的方法研究了均勻流下柔性立管的渦激振動響應特性及渦激力載荷特性，研究發現均勻流下柔性立管的渦激振動響應為位移和主導頻率不隨時間變化的穩態響應。

對于大長徑比的柔性圓柱體，數值模擬全尺度的振動響應存在一定的局限性[11-12]，而試驗方法研究成本相對較高[13-15]，因此有必要建立一種能夠快速預報響應細長圓柱渦激振動的經驗模型方法模擬細長柔性體振動響應主要特征。尾流振子模型是應用較為廣泛的經驗模型。最早Noack等[16]利用范德波爾振子方程推導出能夠預測細長鈍體渦脫頻率和角度的經驗模型，利用帶有擴散項的范德波爾方程能夠準確預測出尾流脫落的定性特征。之后Balasubramanian等[17]在此基礎上利用帶有擴散項的范德波爾振子方程，用來預測均勻流和剪切流下直圓柱和錐形圓柱的漩渦脫落。Facchinetti等[18]利用尾流振子方程和結構柔性振動方程，研究了三維柔性圓柱的渦激振動響應，并與試驗結果對比指出耦合模型能夠準確預測渦激振動響應的振動形態。高云等[19-20]通過尾流振子模型，對線型剪切來流下兩端鉸接的細長柔性圓柱的渦激振動響應做了詳細的探究。

綜上所述，當前研究主要是對細長柔性圓柱的渦激振動響應進行研究，如何控制三維柔性圓柱渦激振動的研究相對較少。因此，本文借助尾流振子模型，提出利用延遲反饋來實現柔性圓柱渦激振動主動控制，研究延遲時間τ和延遲增益kd對柔性圓柱振動幅值和振動頻率的影響，尋求控制柔性圓柱渦激振動的延遲反饋規律，為海洋立管等柔性鈍體的振動控制提供有益探索。

1 數學模型

柔性圓柱流固耦合作用下產生振動幅值主要受到柔性體張力的影響，本文考慮無限長張緊圓柱在受到均勻來流作用下的渦激振動響應，忽略結構的彎曲剛度，模型示意如圖1所示。

圖1 受到均勻來流的無限長張緊圓柱Fig.1 An infinite tension cylinder subjected to a stationary uniform flow

圖中,圓柱直徑為D，長度為L，兩端受到的軸向張力為N，在均勻來流Uf作用下產生振動，本文只考慮y方向的振動位移。參考文獻[21]，本文選取長徑比L/D=25進行討論。在延遲反饋控制的作用下，所選柔性圓柱部分的橫向位移Y可以表示為

(1)

式(1)的右端項P(Z,T)為作用在結構上的流體力，可以表示為P(Z,T)=ρU2DCL0q(Z,T)/4，其中CL0為固定圓柱的參考升力系數[22]，q(Z,T)是用來描述尾流運動的尾跡變量。φ(Z,T)=10-3sin(50πt)為初始位置函數，圖2計算了不同初始位置函數φ(Z,T)下柔性圓柱縱向振動位移的均方值，說明數值計算結果與φ(Z,T)的選取無關。

圖2 不同初始位置函數φ(Z,T)下柔性圓柱縱向振動位移的均方值Fig.2 RMS values of dimensionless displacements at different initial location function φ(Z,T)

尾流的運動可以用非線性振子來表示，常見的表示形式為范德波爾振子，寫作

(2)

式中，A和ε為經驗參數。引入無量綱參數y=Y/D，z=Z/D，t=TΩf，式(1)和式(2)可以寫作如下無量綱形式

(3)

(4)

2 數值方法及模型驗證

2.1 數值方法

(5)

式中，k_d為延遲項所對應的索引，將式(5)代入式(3)、式(4)可以得到如下迭代格式

(6)

(7)

(8)

2.2 模型驗證

首先對數值離散方法的時間步長和空間步長進行無關性驗證。如表1所示，時間步長Δt取較大值0.02和0.03時，對于較小的空間步長Δz，數值計算結果發散。當Δt=0.005，Δz=0.05時，中間節點處最大振幅為1.849，而Δt=0.01，Δz=0.1時，最大振幅為1.847，降低約0.1%；當Δt=0.01，Δz=0.4時，最大振動幅值為1.785，降低約3.4%。因此在后續的數值計算中，取Δt=0.01，Δz= 0.1，整個數值仿真時間Tend=1 400。

表1 不同時間和空間步長下，中間節點處結構振幅Tab.1 The vibration amplitude of cylinder at node under different space and time steps

為了驗證文中所選模型的正確性，取L/D=25，c=4,μ=1.785，St=0.16，與Newman等所用參數保持一致。不考慮延遲反饋作用下，結構的振動位移云圖如圖3(a)所示，橫軸的無量綱時間為t/2πSt。可以發現，在時間段200～240內，三維柔性圓柱的振動形式為駐波；當無量綱時間達到260時，振動形式由駐波逐漸向行波轉變，此后一直保持行波振動狀態，這一過程與Newman等DNS計算結果一致(如圖3(b)所示)，說明本文的計算模型是可行的。

圖3 渦激振動位移響應對比(L/D=25，c=4,μ=1.785)Fig.3 Comparison of VIV displacement response with time and space (L/D=25，c=4,μ=1.785)

3 結果和討論

3.1 延遲時間τ對振幅和頻率的影響

從圖4(a)～圖4(d)均可以看出，在特定的延遲增益kd下，振幅隨著延遲時間τ呈周期性變化，且每個周期內振幅最大值(定義為振幅峰值)保持不變，延遲時間τ并不能改變振幅峰值，其大小只與延遲增益kd有關，隨著kd的增加，振幅峰值要更加偏離未施加延遲反饋時的振動幅值0.647，當kd為0.05，0.10和0.15時，振動峰值分別為0.85，1.18和1.82。從圖4(a)可以看出，當延遲時間τ=0時，振動幅值分別為0.518，0.429和0.363，振動幅值隨延遲增益kd的增大而減小，均小于未受控時的振動幅值0.647，最大能夠降低43.9%。觀察圖4(b)～圖4(d)可以發現，當kd符號相反時，振動幅值變化曲線之間會相差半個周期的相位。隨著τ的增加，振幅最小值都大于τ=0時的振幅值，并且逐漸趨向于振幅峰值，而且隨著kd的增大，趨向的速度越來越快，如箭頭所示。綜上發現，不同的延遲時間τ耦合不同的延遲增益kd，可以實現不同柔性圓柱振動幅值的控制，這將有益于實現工程實踐中柔性鈍體的振動控制。

如圖5所示，振動頻率隨延遲時間也是呈周期性的變化，振動頻率隨著延遲時間τ的增加逐漸趨向于未施加延遲控制柔性圓柱振動頻率f=0.156 4，且通過圖5(b)～圖5(d)可以發現，延遲增益kd越大，趨向的速度逐漸越快，如箭頭所示。

3.2 延遲增益kd對振幅和頻率的影響

圖6(a)和圖6(b)的延遲時間τ分別是圖4(a)振幅取最大和最小時所對應的延遲時間。從圖6(a)中可以看出，當kd由0逐漸增加時，振幅也逐漸增加，振幅只與kd有關。當kd從0逐漸減小，在τ取較小值(3.2，9.6和15.9)時，延遲增益對振幅的影響占主導，振幅隨kd減小而減小；當τ取值較大的情況下，延遲時間對振幅的影響占主導地位，振幅隨kd的減小呈現先減小后增大的趨勢。在相同負值kd下，延遲時間越大振幅越大，這與圖4所示振幅最小值逐漸趨向于峰值的變化趨勢是一致的，而且kd越大趨向速度越快。由于延遲增益kd由正變負時變化曲線會相差半個周期的相位，因此圖6(b)中振幅隨kd的變化規律與圖6(a)中所描述的規律剛好相反。

圖4 延遲時間(τ)對振動幅值的影響Fig.4 The effect of delay time (τ) on the vibration amplitude

圖5 延遲時間(τ)對振動頻率的影響Fig.5 The effect of delay time (τ) on the vibration frequency

圖6 延遲增益kd對振動幅值的影響Fig.6 The effect of delay gain (kd) on the vibration amplitude

為了探究柔性圓柱在延遲反饋作用下的能量轉化，根據高云等的研究引入瞬時能量轉化公式

(9)

式中，v(z,t)，y(z,t)和q(z,t)分別為瞬時速度、位移和尾跡變量。瞬時能量W為正值表示流場向結構傳遞能量，從而激發柔性圓柱振動；負值表示結構向流場傳遞能量。

圖7(a)所示為柔性圓柱中間節點處系統振動變量的值，可以發現未施加延遲反饋時變量q的最大值為4.821，當尾跡變量q增加時，振動位移y也相應增加。對于振動幅值較小的情況(τ=6.5)開始振動時間要小于振動幅值較大的情況(τ=3.2)，而且都要小于未施加延遲反饋時的開始振動時間，說明施加延遲反饋控制能夠降低系統的起振時間。圖7(a)中瞬時能量W變化曲線表明，W的最小值越小，對應產生的振動幅值越大。例如對于τ=3.2而言，此時W最小值為-0.07，要小于τ=6.5時的-0.002，但前者對應的振動幅值y更大。結合瞬時能量W的含義，負值越小表明結構傳遞到流場中的能量越多，流場中能量增加會導致尾跡變量q值增大，從而又會反過來增加結構的振動y，因此通過能量云圖7(b)可以看出，當τ=3.2時，由于能量在結構和流體中相互傳遞，駐波向行波演變的時間會增加(白色虛線框對應的時間)，對于抑制振動的情況(τ=6.5)，駐波向行波演變時間會減小。

上述可以發現瞬時能量W最小值對系統振動幅值增加或者減小有很大的影響。而結構振動增強的能量來源是外部延遲反饋控制系統，控制系統的延遲時間τ決定系統當前振動狀態與τ時間之前振動狀態之間的耦合關系，這種耦合關系會降低或者增加結構振動幅值。因此延遲增益kd決定了系統振動幅值的上限(振幅峰值)，延遲時間τ決定系統當前振動狀態與之前振動狀態的耦合關系，結合圖4(a)的結論，可以看出延遲時間τ并不能改變振幅峰值，峰值大小只與延遲增益kd有關。

圖7 當kd=0.15時，延遲時間(τ)對系統的影響Fig.7 The effect of delay time (τ) on VIV system when kd=0.15

圖8(a)和圖8(b)中的延遲時間分別是圖5(a)頻率取最大和最小時所對應的延遲時間。從圖8(a)中可以看出，當kd由-0.15逐漸增加到0.15時，對于延遲時間τ=1.4，τ=7.4和τ=13.4而言，振動頻率隨kd的增加而增加；當τ=19.4和τ=25.6時，振動頻率先保持不變，kd進入正值以后頻率逐漸減小，在kd=0.12時頻率又突然增大，最終值和其他延遲時間保持相同。對于圖8(b)而言，振動頻率隨著kd增大都呈現逐漸減小的趨勢。

圖8 延遲增益kd對振動頻率的影響Fig.8 The effect of delay gain (kd) on vibration frequency

從圖9可以發現，未施加延遲反饋控制時，振動幅值為0.647，行波傳遞周期為6.41，振動頻率為0.156(1/6.41)；當τ=4.7時，振動幅值增大為0.899，行波傳遞周期為6.71，振動頻率減小為0.149；當τ=1.4時，振動幅值減小為0.521，行波傳遞周期為5.95，振動頻率增大為0.168。相較于未施加延遲的情況，振動幅值和振動頻率變化呈相反的趨勢，振動幅值增大，振動頻率就會減小。

圖9 振動位移和行波頻率時間的關系(kd =0.15)Fig.9 The relationship between vibration displacement and frequency of travelling wave (kd =0.15)

3.3 理論解釋

對于耦合方程式(3)和式(4)，設方程組的解為諧波行波的形式，寫作

y(z,t)=y0ei(kz-ωt-φ),q(z,t)=q0ei(kz-ωt)

(10)

式中：結構位移y和流體變量q都具有共同角頻率ω以及波數k；φ為變量之間的相位角，y0和q0為與時間無關的幅值。將式(10)代入式(3)和式(4)中，得到

(11)

式(11)對角項分別為結構和流體的色散關系，寫作

(12)

(13)

為了得到實數角頻率ω和波數k，式(11)的行列式為0，由此得到流固耦合系統的色散關系為

DFS(ω,k;τ,q0)=

DS(ω,k;τ)DF(ω,k;q0)+AMω2=0

(14)

令式(14)實部和虛部為0，分別得到

-(ω2+c2k2)(-ω2+1)-

(15)

(16)

分析式(15)和式(16)可以發現，振動角頻率ω、波數k以及尾跡變量q0與延遲參數kd和τ有關。首先不考慮延遲參數，得到角速度ω和q0隨波數k的變化，如圖10所示。可以發現，有無數個(ω，k，q0)組合滿足上述兩個方程，特定的組合與給定的邊界條件有關。對于無限長柔性圓柱而言，最大的振動幅值y0和尾跡變量q0會發生在ω=1處，此時波數k=0.25，q0=4.89。

確定的尾跡變量q0之后，根據式(15)和式(16)，ω和k與延遲參數相關。從圖11中可以發現，隨著τ的增加，頻率f=ω/2π逐漸趨向于0.159，這一變化規律與圖5中的變化規律一致，同時波數k也逐漸趨向于0.25。通過圖11(c)可以看出，當τ較大時，ω/k的值基本保持為4不變，這與Facchinetti等研究中的結論ω/k=c相一致。通過上述分析發現，當延遲時間τ不斷增加時，振動頻率f和波數k都趨向于未受控時的值，分別為0.159和0.25。

圖10 振動頻率ω和尾跡變量q0隨波數k的變化趨勢Fig.10 The variation of vibration frequency ω and wake variable q0 with wave number k

圖11 頻率ω、波數k以及ω/k隨延遲時間τ變化規律(kd=0.15)Fig.11 The variation of frequency ω,wave number k,and ω/k with delay time τ (kd=0.15)

4 結 論

本文以三維細長柔性圓柱為研究對象，基于結構梁振動模型和尾流振子模型，利用二階中心有限差分方法求解耦合方程，通過引入時滯反饋控制，討論了在不同的反饋參數下(延遲時間τ和延遲增益kd)柔性圓柱渦激振動的幅值和頻率變化情況。主要得出如下結論：

(1)探究了延遲時間τ對系統振幅和頻率的影響，發現三維柔性圓柱振動幅值和頻率隨τ的增加呈周期變化；相反的延遲增益會使振幅變化曲線相差半個周期的相位角，而且系統能夠達到的振幅峰值隨著kd的增大而增大。

(2)探究了在特定τ下kd對系統振幅和頻率的影響，發現延遲增益kd決定了系統振動幅值的上限；且振動頻率和幅值的變化呈相反的趨勢，振動幅值增大，振動頻率越小。

(3)通過理論分析得出流固耦合系統的色散關系，發現振動頻率ω和波數k隨著延遲時間τ的增加最終會保持不變；振幅峰值只與延遲增益kd有關，kd越大振幅峰值越大。