渦激響應尾流雙振子模型的參數反演方法

顏凝雨，楊自豪*，張潔瓊，林 巍，劉孟源

(1.西北工業大學 數學與統計學院，西安 710072；2.中交懸浮隧道工程技術聯合研究組，珠海 519000)

1 引 言

任何非流線型物體在一定的恒定流速下，都會在物體兩側交替地產生脫離物理結構表面的周期性漩渦，從而誘發彈性支撐結構物的渦激振動現象VIV[1]。大量實際工程問題中，許多圓柱結構，如懸浮隧道、深水立管、高空桅桿和新型風力發電等，均可能出現渦激振動現象。雖然合理利用渦激振動會為人們帶來便利，但其誘發的共振也會對結構產生較大的疲勞損害。

研究圓柱結構的渦激振動主要有試驗法[2-4]、數值方法[5,6]和經驗模型[7,8]。試驗方法雖然直觀可靠，但成本較高；數值方法雖然在計算精度上占有一定優勢，但計算較為復雜，尤其是在高雷諾數的情況下。而經驗模型可較快計算渦激振動的主要特征，還可以在模型的指導下設計試驗或者研究某個參數對試驗結果的影響，因此經驗模型法正受到國內外諸多學者的關注，其中尾流振子模型是應用較為廣泛的一種模型。

Bishop等[9]認為流體的尾流對圓柱體結構的作用相當于一個非線性振子，首次提出了尾流振子模型。之后，Hartlen等[10]將Vander Pol方程與結構振動方程聯立，使尾流振子模型更加具體化。Skop等[11]進一步改進了Hartlen的模型，對模型經驗參數進行了研究；Landle[12]提出了一個擁有五階非線性氣動阻尼項的無量綱尾流振子模型；Facchinetti等[13]通過分析對比發現加速度耦合模型可以更準確地定量預測渦激振動的基本特征；Ogink等[14]對Facchinetti的尾流振子模型建立過程進行了具體分析，通過引入頻率耦合變量進一步改進了模型。但上述對尾流振子的研究大多針對橫向單自由度，而實際的渦激振動中橫向振動和順流向振動往往相互影響，許多學者也針對雙自由度的尾流振子模型做了研究。Krenk等[15]根據流固耦合的能量流平衡原則構建了雙振子模型；Srinil等[16]提出了一種改進的用于預測均勻流中柔性圓柱體的二維耦合尾流雙振子模型。康莊等[17]建立了非線性結構振子和流體振子的耦合模型，可以較準確預報低質量比圓柱結構雙自由度渦激振動的特征。Yang等[18]提出一種新的渦激振動非線性耦合尾流振子模型，其只使用一個尾流振蕩器來耦合橫向和順流向的運動。在實際應用中，不是所有模型的參數都可以通過物理方法測量得到，尾流振子模型的應用往往面臨參數選取的難題，參數選取的偏差可能導致計算結果的不準確，從而無法充分發揮經驗模型在研究渦激振動問題中的優勢。

解決復雜模型參數反演問題的方法有很多種，其中神經網絡方法以其高速尋優和自學習等優勢在參數反演問題中得到廣泛的應用，如Spichak等[19]運用BP神經網絡求解地電三維逆問題，證實神經網絡方法可以成功地用于地電三維逆問題參數反演；Cerdena等[20]在測定海洋層積云物理性質的研究中，運用遺傳算法優化的BP神經網絡對輻射傳輸模型中的參數進行了反演；Zhao等[21]提出一種人工蜂群算法優化的BP神經網絡，并對巖土工程中的模型參數進行了反演。國內外沒有針對尾流振子模型參數反演的研究工作。

本文研究尾流雙振子模型中參數的反演方法。首先利用降階法將模型變為一個常系數多變量的非線性一階常微分方程組，并利用龍格-庫塔方法進行求解。設計懸浮隧道物理模型的水槽試驗，在給定的多組物理參數和試驗數據下，構建最優化模型，利用遺傳算法優化的BP神經網絡智能算法反演未知參數。

2 渦激響應尾流雙振子模型

柔性圓柱結構雙自由度渦激振動模型[16]如 圖1 所示。假定直徑為D的圓柱結構無限長，水平和豎直方向，即X和Y兩個方向的彈簧剛度分別為Kx和Ky,水的流速為V,圓柱體受到雙向四個彈簧的作用力，在X和Y兩個空間方向上自由振動，無量綱形式的渦激響應尾流雙振子模型可表示為[16]

(1)

圖1 柔性圓柱結構雙自由度渦激振動

λx=2ξxf*+γΩ/μx,λy=2ξy+γΩ/μy

(2)

式中ξx和ξy為結構水平和豎向運動阻尼比，取值范圍為[0.01,0.08]；f*為X和Y兩方向結構固有頻率之比，取值范圍為[0.5,1.5]，γ為Stall參數，取值為0.8；Ω=St·Vr為靜水中旋渦脫落頻率與圓柱橫流固有頻率之比，St為Strouhal數，取值為0.2，Vr=V/(Dfn)為約化速度,fn為Y方向結構固有頻率，μx和μy表達式為

(3)

此外，(αx,αy,βx，βy)為無量綱幾何非線性系數，表征位移非線性項對本身的影響。MD和ML為系統質量參數，

(4)

式中CD 0和CL0分別為靜置圓柱結構的阻力和升力系數，取值為0.2[16]和0.3；St為斯特魯哈爾數，取值為0.2；εx和εy為尾流系數；Λx和Λy為經驗耦合參數，均取為12。

對于質量分布均勻的圓柱體，有

(5)

且可取εx=0.3，εy由試驗數據擬合得到，其擬合表達式為

εy(m*)=b1eb2 m*

(6)

式中b1=0.00234，b2=0.228。

由式(1)可知尾流雙振子模型是一個多變量的二階非線性常微分方程組，現運用降階法將該模型化為一個常系數多變量的一階非線性常微分方程組。引入如下變量

u=(u1,u2,u3,u4,u5,u6,u7,u8)T=

(7)

式中T為矩陣的轉置。方程組(1)可轉化為一階常微分方程組(8)。

(8)

進一步，式(7，8)可以改寫為

(9)

式中M和g的定義如下

(10)

g(t,u,αx,αy,βx,βy)=(g1,g2,g3,g4,g5,g6,g7,g8)T

式中

g1=u5，g2=u6，g3=u7，g4=u8

MDΩ2u2-2πMLΩ2(u4u7/Vr)

MLΩ2u4+2πMDΩ2(u2u7/Vr)

(11)

給定初值u0,該一階非線性常微分方程系統可使用經典的龍格-庫塔方法進行求解。

3 圓柱結構渦激振動響應水槽試驗

3.1 試驗目的

尾流雙振子模型中(αx,αy,βx，βy)四個參數表征位移非線性項的影響，無法通過試驗直接測得，其他參數均可通過試驗或經驗獲得。若想得到上述四個參數的值，由模型(9)可知，一種可行的方案是通過圓柱體渦激振動時的位移數據來反演這四個未知參數。因此，本節設計圓柱結構的水槽試驗，測得圓柱結構在不同水流力(流速)作用下水平和豎直方向的位移時間序列及運動軌跡。

3.2 試驗設計

現有的圓柱結構渦激振動試驗不考慮懸浮隧道等實際工程結構試驗節段的真實剛度，試驗結果無法反映真實結構的剛度和自振頻率，試驗結果與實際值誤差較大。本節給出新型的圓柱結構水槽試驗設計方案。

首先運用等效固有振頻法計算懸浮隧道等實際工程結構某一節斷的真實剛度值k0,計算公式為

k0=(4π2m)/T2

(12)

式中m為某一節段的質量，T為該節段在水中的固有頻率。其次根據相似準則，將真實剛度k0按照比尺換算成試驗模型的剛度k1,并據此配置水平和豎向彈簧的剛度。最后將試驗模型放入水池中進行試驗，浮筒式試驗裝置如圖2所示。

圖2 試驗裝置設計

試驗采用直徑為8 cm和長度為70 cm的透明樹脂材料制成的空心圓柱體模擬深海懸浮隧道等實際工程結構的節段(等剛度模型)，模型重 3.3 kg，浮重比1.07，試驗比尺為 1∶157.5。圓柱體使用4根彈簧與透明固定板鏈接，其中橫向和豎向彈簧剛度分別取為30.6 N/m和21.6 N/m，彈簧原始長度均為100 mm，圓柱體、彈簧和透明固定板之間均采用L型鋼扣栓接，如圖3(a)所示。靜止狀態下，彈簧長度為192 mm，如圖3(b)所示。

選取天津水運工程科學研究院水槽廳的波流水槽為試驗場地，水槽寬70 cm，長30 m，水深 57 cm，造流范圍為0.02 m/s～0.5 m/s，水槽尾端采用斜坡式消波框防止波浪反射。將試驗模型放入0.5倍水深處進行試驗，并在模型截面貼上紅色方塊標志點，如圖3(b)所示。不斷改變水的流速，利用相機對標志點進行動態捕捉，然后利用 Matlab 編寫的視頻識別程序，抓取每一幀中標志點中心的位置，并將前后兩幀位置數據進行比較，從而得到圓柱體在不同物理流速下水平和垂向位移隨時間變化的數據。試驗共進行9組，每組試驗所對應的水的流速和約化速度列入表1。

圖3 試驗裝置

表1 9組試驗中水的流速和約化速度Tab.1 Water velocity and reduced velocity in 9 sets of tests

3.3 試驗結果

不同流速下圓柱結構的運動軌跡經過無量綱處理后如圖4所示。可以看出，運動軌跡呈8字形，且隨著流速的增大，模型在豎直方向的運動幅度增大，在第五組達到最大，考慮可能達到了鎖定流速區間，即當流速在此區間內時，圓柱結構振動頻率近乎不變，振幅維持較高水平；且隨著流速增加，模型所受橫向拖拽力增加，無量綱橫向位移一直增大，這與以往研究保持一致[8]。

4 參數反演的BP神經網絡模型

4.1 參數反演的最優化模型

根據第3節中每組試驗，可以得到對應的位移試驗數據

(i=1,…,n;ti∈[0,tmax](13)

而每給定一組參數(αx,αy,βx，βy),求解微分方程系統(9)，可得到對應的位移計算數據，記作

x(t,αx,αy,βx,βy),y(t,αx,αy,βx,βy)

(t∈[0,tmax])

(14)

下面構建參數反演的最優化模型，即求(αx,αy,βx,βy),使

(15)

由于本文關注的核心為振子的幅值和均方根，因此可引入向量函數

opt(x,y)=[rms(x),rms(y),max(x),max(y),

min(x),min(y)]T

式中max(·)和min(·)分別為向量的最大和最小分量，rms(·)為均方根。

N為向量x的長度。目標函數可以取為

(16)

求解最優化問題(15，16)即可得到未知參數的一組最優值。

圖4 不同流速下圓柱的運動軌跡

4.2 BP神經網絡的構建

人工神經網絡是一種由大量處理單元互聯組成的非線性、自適應信息處理系統，具有快速求解最優化解的能力，其最小最重要的單元叫神經元，每個神經元都有輸入連接和輸出連接，這些連接模擬了大腦中突觸的行為，每一個連接都有對應權重。BP神經網絡由輸入層、隱藏層和輸出層構成，是一種按誤差反向傳播訓練的多層前饋網絡。BP神經網絡求解優化問題的基本思想是利用BP算法，通過信息的正向傳播和誤差的反向傳播，不斷調整網絡中的權值和閾值，使輸出和樣本的實際標簽之間的誤差達到最小。渦激響應尾流雙振子模型的參數反演思路為，生成多組未知參數數據，通過求解尾流雙振子模型(9)得到每組未知參數下的位移時間序列數據；分別計算每組豎直和水平位移數據的最大值、最小值和均方根，并以這6個數據作為神經網絡的輸入，以未知的4個參數為輸出，優化目標為式(15)，即可訓練出對應的神經網絡模型；將試驗得到的一組位移數據代入已訓練好的神經網絡，即可輸出一組未知參數的值。

4.3 基于遺傳算法優化的BP神經網絡

BP神經網絡存在局部極值，誤差較大。現采用遺傳算法對BP神經網絡進行優化，即使用遺傳算法對BP神經網絡的初始權值和閾值進行優化，進一步提高參數反演的精度。遺傳算法以優勝劣汰的自然選擇和基因階段的遺傳變異原理作為基礎，是一種具有很強隨機搜索能力的智能算法，可以取得很好的搜索效果。遺傳算法對BP神經網絡的優化，實際上是對BP神經網絡的初始權值和閾值進行優化。首先確立BP神經網絡結構并得到一組初始權值和閾值，遺傳算法將對其編碼，并再次通過神經網絡訓練，得到的誤差作為適應度。然后不斷經過選擇、交叉和變異操作使適應度逐漸增加，最終滿足精度條件，得到最優權值和閾值。最后再代回BP神經網絡進行求解。具體流程如圖5所示。

4.4 參數設置及誤差定義

尾流雙振子模型的參數取值列入表2。

反演數據與試驗數據的誤差計算公式為

(17)

表2 參數取值

圖5 求解參數反演問題分析流程

5 參數反演結果分析

利用遺傳算法、BP神經網絡算法和基于遺傳算法優化的BP神經網絡算法反演參數。遺傳算法中初始種群數目取為20，進化代數取為80，交叉概率取為0.95，變異概率取為0.001。BP神經網絡由一個輸入層、一個隱藏層和一個輸出層構成，網絡模型參數選取如下，輸入層神經元個數取為6，隱藏層神經元個數取為8，輸出層神經元個數取為2，數據庫樣本個數取為8000，訓練集個數取為7200，驗證集個數取為800，訓練誤差目標取為10-4，學習效率取為0.4，傳遞函數取為{tansig,purelin}，未知參數αx和αy的取值范圍為(0,1)，未知參數βx和βy的取值范圍為(0,10)。此外，每組試驗數據分別訓練50次，取平均值作為參數的最終預測結果。需要指出的是，由于試驗誤差，試驗得到的位移數據沒有呈現穩定后的周期震蕩模型，但大致穩定在某個固定區間之內。為了減小參數反演誤差，現對神經網絡的輸入參數，即豎直和橫向位移的最大和最小實驗數據做以下處理，以第3組試驗的水平方向位移最大值為例，利用Matlab程序找到數據的峰值，再求出峰值的平均值作為神經網絡位移最大值的輸入，如圖6所示。基于遺傳算法優化的BP神經網絡算法流程如圖5所示，其中遺傳算法和BP神經網絡的初始參數如前所述。

圖6 試驗數據輸入最大值的處理

利用遺傳算法、BP神經網絡算法和基于遺傳算法優化的BP神經網絡算法反演參數的誤差Ret列入表3，誤差隨水的流速變化如圖7所示。由表3可知，三種方法所得結果的平均誤差分別為10.37%、7.53%和5.50%，基于遺傳算法優化的BP神經網絡算法的平均誤差相比遺傳算法和BP神經網絡算法分別降低了46.96%和26.96%。此外，由圖7可知，遺傳算法的穩定性最差，BP神經網絡次之，遺傳算法優化的BP神經網絡的穩定性最好。遺傳算法優化的BP神經網絡參數參演結果列入表4。

圖7 三種模型方法的誤差比較

圖8為利用遺傳算法優化的BP神經網絡的反演結果與第3、4和第6組試驗無量綱位移數據的對比，其中圖8左圖分別表示第3、4和第6組試驗方法與反演方法水平和豎直方向無量綱位移數據的對比，圖8右圖分別表示第3、4和第6組試驗數據與反演數據水平和豎直方向運動軌跡的對比。可以看出，基于遺傳算法優化的BP神經網絡算法反演結果與試驗結果吻合較好。因此，基于遺傳算法優化的BP神經網絡智能算法能夠精確穩定實現尾流雙振子模型的參數確定。

表4 遺傳算法優化的BP神經網絡參數反演結果Tab.4 Parameter inversion results of BP neural network optimized by genetic algorithm

最后以第3組參數反演問題為例進行靈敏度分析，研究反演參數的微小偏差對反演精度的影響。將兩個參數值分別在參數反演結果附近微小擾動，圖9給出了誤差Ret的值，其中灰色點代表遺傳算法優化的BP神經網絡參數反演結果。

可以看出，參數反演結果處于在局部最優解附近，進一步驗證了反演方法的正確性。此外，通過計算誤差沿兩個參數的變化率可發現，誤差Ret對參數βx和βy更加敏感。在反演結果(0.8214±0.2，10.0674±0.5)附近的參數取值(灰色實線區域)對應的誤差Ret約為6%，在反演結果(0.8214±0.5，10.0674±1)附近的參數取值(灰色虛線區域)對應的誤差Ret約為10%。該靈敏度分析結果可為類似問題的參數取值提供可靠的理論分析基礎。

圖8 第3、4和第6組試驗位移數據與反演數據對比

圖9 第3組參數反演問題中的參數取值域

6 結 論

本文主要研究了尾流雙振子模型中未知參數的反演方法，給出了基于物理試驗和BP神經網絡模型的參數反演方法。利用降階法將尾流雙振子模型轉換為一階常微分方程組，并使用龍格-庫塔方法實現模型的求解。給出了一種新型的圓柱結構水槽試驗設計方案，獲取了圓柱結構在水平和豎直方向的位移試驗數據。基于試驗數據，構建了最優化模型，并利用基于遺傳算法優化的BP神經網絡算法對尾流雙振子模型中四個未知參數進行了反演，誤差處于可接受范圍之內。本文的工作將為懸浮隧道等實際工程結構的渦激振動響應分析提供理論基礎。