基于子結構中彎曲波的拉索索力識別方法

肖軍，代洋，張永水

(1.中交第二公路工程局有限公司，西安 710065；2.長安大學 公路學院，西安 710064；3.重慶交通大學 土木工程學院，重慶 400074)

拉索是斜拉橋的重要受力構件，其索力的準確識別在橋梁施工與運營階段十分重要。工程中使用最為廣泛的是動測法，又稱頻率法[1]。頻率法依賴于對拉索振動頻率的識別以及拉索頻率與索力的換算關系，其精度很大程度上取決于頻率與索力的換算關系。早期利用弦理論研究拉索振動，忽略了抗彎剛度，并且假設邊界條件為兩端鉸接，與實際情況存在較大偏差，因而識別結果并不理想[2-3]。由于頻率法本質是利用振動頻率對索力進行識別，抗彎剛度、邊界條件以及垂度對識別精度影響都比較顯著。由于對拉索抗彎剛度的計算缺乏精確的理論公式，通常只是給出合理的取值范圍，導致拉索索力識別時其抗彎剛度反而成為了待識別對象[4-5]；文獻[6-9]對不同邊界條件的情況做了大量分析，用彈簧剛度模擬邊界條件在靜力分析中可行，但對于動力分析，其動態邊界的剛度顯然依賴于其動力特性，因而，以截斷性邊界模擬拉索邊界不合理；同時，拉索的垂度使得拉索的索力是沿索長分布的函數，并非固定值，傳統索力識別方法無法只識別索力的平均值。鑒于頻率法依然存在的一些問題，部分學者轉向了利用行波進行索力識別的研究，McDaniel等[10]推導了梁的動力響應通解，不直接通過邊界條件建立特征方程求出待定系數，而是用不同測點的頻響反求待定系數，當測點數多于待定系數時，將有擬合誤差，通過擬合誤差最小化實現了梁的波數識別。Maes等[11]關注到梁的軸力與波數具有一一對應的關系，因此，將波數識別推廣應用于梁桿軸力識別，該方法在頻域中產生大量識別點集，提供了更大信息量，提升了穩定性。張松涵[12]提出了一種索力識別方法理論，利用選取的子索段的5個測點對波分量系數進行最小二乘求解，以波分量系數的擬合殘差最小作為索力識別的判定標準，由于其代入索力識別的波分量仍采用了Euler-Bernoulli梁模型的波數解，顯然在高頻響應中精度無法滿足。筆者對梁單元進行了修正，解決了Euler-Bernoulli梁模型對短粗梁以及高頻率段的不適用性，避免了Timoshenko梁模型存在的截斷頻率、兩個波速系的問題，對梁模型中4種波的特性進行探討，認為近場波由錨固端向梁中呈指數衰減只存在梁錨固局部位置處，并且隨著頻率的增加衰減得越快，忽略近場波后，通過3個測點的頻域響應采用通過最小二乘法擬合得到波分量系數，以擬合殘差最小為目標進行索力和抗彎剛度的識別，最后，通過拉索振動的數值模擬實驗驗證了方法的精確性。

1 拉索振動頻散關系

1.1 修正Timoshenko梁理論

Doyle[13]推導了Euler-Bernoulli梁理論的頻散關系，Lee等[14]在此基礎上考慮了剪切變形和軸向張力，推導了Timoshenko梁理論下的振動頻散關系，但Euler-Bernoulli梁模型由于忽略了剪切變形，以致在高頻段的誤差較大；而Timoshenko梁理論存在截止頻率，使得其具有兩個波速系，這不符合實際情況，其實，在Timoshenko梁理論的推導中，未引入剪切變形所引起的轉動慣量，一旦考慮之后，便可消除截止頻率，只留下一個波速系，并且增加了結構振動高頻響應的精確性[15]。

拉索微元段引入剪切變形所引起的轉動慣量后，平衡狀態如圖1所示。

圖1 微段平衡示意圖Fig.1 Cable balance diagram

根據圖1建立拉索微段平衡方程組

(1)

式中：y(x,t)為拉索的橫向位移；η(x,t)、λ(x,t)分別為拉索彎曲和剪切引起的截面轉角。

由Timoshenko梁[16]可知，剪力、彎矩、拉索橫向位移存在以下關系

(2)

(3)

式中：κ為截面的剪切變形系數，拉索的圓形截面可按式(4)計算；G為材料的剪切模量，對于各項同性的材料，G可以按式(5)計算。

(4)

(5)

將式(2)～式(5)代入式(1)，并進行Fourier變換，可得到拉索在頻域中的橫向振動微分方程

(6)

(7)

(8)

并且，

(9)

經Fourier變換后可得

(10)

(11)

要使式(11)存在非零解，則左端系數矩陣的行列式必然為0，可寫成

EI(κGA+Nx)k4+(-NxκGA+ω2ρAEI+

ω2ρIκGA)k2-ω2ρAκGA=0

(12)

對特征方程式(12)求解，便可得到修正Timoshenko梁的頻散關系

(13)

式中：α=EI(κGA+Nx)；β=-NxκGA+ω2ρAEI+ω2ρIκGA；γ=-ω2ρAκGA。

1.2 數值算例

基于梁理論的推導，結合一數值算例進一步探討修正Timoshenko梁、Euler-Bernoulli梁與Timoshenko梁理論的關系與區別。

對于一段連續均勻、等截面并且不考慮其長度的梁結構，假設其材料參數為：彈性模量E=200 GPa，泊松比μ=0.3，由式(4)算出截面剪切變形系數κ=0.86，由式(5)算出剪切模量G=76.92 GPa，密度ρ=7 800 kg/m3；幾何參數：圓形截面半徑為0.04 m，截面面積A=0.005 m2，截面慣性矩為Izz=2×10-6m4，對梁施加的初張力假設為Nx=600 MPa×0.005 m2=3 000 kN。

圖2給出了3種梁單元波數解與頻率的關系，其中，實部代表近場波，虛部代表行波。在頻率較低的情況，3種梁理論的頻散關系相差很小，但在較高頻段，Euler-Bernoulli梁的近場波數解存在無限增大的趨勢，這顯然是忽略了剪切變形和抗彎剛度造成的。而對于Timoshenko梁存在的截斷頻率，其實質是由于只考慮了彎曲變形產生的抗彎剛度，使得波數解中出現了頻率的四次方，而修正Timoshenko梁額外考慮了剪切變形引起的轉動慣量，恰好抵消掉了該項，因而避免了截斷頻率的產生[15]。

圖2 波數-頻率圖Fig.2 Wave number-frequence diagram

2 拉索動力響應的波分量分解

2.1 波分量理論

式(13)給出了修正Timoshenko梁的振動頻散關系，根據疊加原理可以得到忽略垂度的拉索頻域橫向自由振動通解(后文波數解均依賴修正Timoshenko梁理論)，即

Y(x,ω)=C1exp(k1x)+C2exp(k2x)+

C3exp(k3x)+C4exp(k4x)

(14)

拉索的頻域橫向自由振動為4項波分量exp(k1x)、exp(k2x)、exp(k3x)、exp(k4x)的疊加，對4種波分量做下列定義：

C1(ω)exp(k1x)[Re(k1)<0,Im(k1)=0]為近場波，沿x正方向衰減；

C2(ω)exp(k2x)[Re(k2)<0,Im(k2)=0]為近場波，沿x負方向衰減；

C3(ω)exp(k3x)[Re(k3)=0,Im(k3)<0]為行波，沿x正方向傳遞；

C4(ω)exp(k4x)[Re(k4)=0,Im(k4)<0]為行波，沿x負方向傳遞。

假設拉索兩端為固定約束，對式(14)代入邊界條件，可得到方程

H·C=0

(15)

式中

(16)

H=

(17)

同樣，要使拉索位移存在非0解，則系數矩陣行列式必為0，即ω={ωn||H(ωn)|=0}，因而，可以得到結構的模態分解表達式

(18)

2.2 數值算例

如圖3所示，一連續、均勻、長度為10 m的兩端固結梁，假設其材料參數：彈性模量E=200 GPa，泊松比μ=0.3，由式(4)算得截面剪切變形系數κ=0.86，由式(5)算得截面剪切模量G=76.92 GPa，密度ρ=7 800 kg/m3；幾何參數：圓形截面半徑0.04 m，截面積A=0.005 m2，截面慣性矩Izz=2×10-6m4。

圖3 結構示意圖Fig.3 structure diagram

將各參數代入式(17)，求得det(H)的值在0～1 000 Hz頻率段中的分布，如圖4所示，各個極小值點對應梁的固有頻率。

圖4 det(H)分布圖Fig.4 det(H)Distribution

式(14)表明拉索在各頻率點的響應由4種波疊加而成，為了進一步探討各種波的性質，將其分為兩組：

Y(x,ω)=C1exp(k1x)+C2exp(k2x)

(19)

為近場波分量；

Y(x,ω)=C3exp(k3x)+C4exp(k4x)

(20)

為行波分量。

在圖4中，從低頻段中選取一個固有頻率，頻率值為19.387 Hz，為第3階固有頻率；同樣，從高頻段中選取一個固有頻率，頻率值為982.463 Hz，為第25階固有頻率。在以上這些頻率處，求得式(14)的基礎解系，再按照式(19)、式(20)將近場波與行波分別開來，結果如圖5所示。

圖5 波分量示意圖Fig.5 Wave component diagram

由圖5可知，近場波僅存在于邊界附近，以指數形式衰減，隨著頻率增大，衰減速度越快。

3 基于子結構彎曲波的索力識別方法

3.1 識別方法理論

由公式推導可知，拉索中任何位置的動力響應在每個頻率點均可寫為4個波分量的疊加，第2節對4個波分量的特性進行了研究，結果表明，近場波衰減迅速，一般只存在固結處相當小的范圍，且在高頻段內更忽略不計[17]，因而，可將拉索振動響應寫為

Y(x,ω)=C3exp(k3x)+C4exp(k4x)

(21)

于是，可以選擇拉索的某一小段作為研究對象，子索段中任意一點的響應依然滿足式(21)，與常規求解代入邊界條件不同，代入拉索內部測點的動力響應求解，因而避免了復雜的邊界條件[18]討論，現假設子索段上布置了M個測點，則

(22)

為了避免Fourier變換產生的譜泄露問題，Igawa等[19]提出利用Laplace變換來替代Fourier變換，取得了非常好的效果，將各測點的動力響應結果進行Laplace變換，轉換到頻域中，得到

(23)

式(21)為拉索振動方程的通解，各測點的動力響應均應滿足。于是，可以得到矩陣方程

(24)

為系數矩陣，取決于結構特征與外部激勵；

為觀測向量，為式(23)的第j列。

如果拉索的參數、測點布置以及各測點的響應結果已經得到，可以通過最小二乘法對波分量系數進行求解。

(25)

如果n<2，C(sj)存在無數解，無法進行參數識別；如果n=2，C(sj)只存在唯一解，依然無法進行參數修正；如果n>2，C(sj)有最小二乘解，存在擬合殘差：

(26)

若結構特征矩陣的各參數取值完全正確，并且觀測結果不存在噪音干擾，那么擬合殘差為0。實際工程中，索力作為待識別對象，無法正確估計，于是，可對索力值做參數修正，然后以標準化擬合殘差達到最小作為索力識別的判定標準，如式(27)所示。

P={N,EI}=

(27)

3.2 索力識別數值算例

以拉索模型為例，驗證子索段結構索力識別方法，假定拉索長度L=100 m，傾斜度sinα=0.6。材料參數：彈性模量E=200 GPa，泊松比μ=0.3，根據式(4)算得截面剪切變形系數κ=0.86，根據式(5)算得剪切模量G=76.92 GPa，密度ρ=7 800 kg/m3；拉索幾何參數：圓形截面半徑0.04 m，截面積A=0.005 m2，截面慣性矩Izz=2×10-6m4，在拉索的兩端均設置了0.5 m的硬索夾段，其抗彎剛度取值為拉索的20倍，近似模擬塔梁對其產生的影響，測點布置在離索夾外2 m位置處，測點間隔1 m，連續布置3個，如圖6所示。

圖6 拉索示意圖Fig.6 Cable diagram

采用ANSYS建立拉索動力模型[20]，假定初拉力為3 000 kN，錘擊位置設在離索夾外1 m處，假設錘擊力為三角形脈沖形式，取其幅值為500 N，作用的時長為t=0.01 s，采樣的頻率取100 Hz，采樣的點數為N=212=4 096，荷載的時域圖、Fourier系數譜見圖7。

圖7 荷載示意圖Fig.7 Force diagram

將求得的3個測點的時域動力響應經Laplas變換到頻域內(如圖8所示)，并按式(27)組裝形成測點觀測向量，按式(25)代入拉索的各個參數，測點的位置以第一個測點為0點，沿拉索方向建立x軸，形成結構特征矩陣。

圖8 頻域響應示意圖Fig.8 Frequency domain response diagram

通過式(30)在各頻率點計算擬合殘差，識別索力，為方便觀察，對各索力值的擬合殘差進行繪制，如圖9所示。識別值為3 026.8 kN，誤差僅為0.9%，具有相當高的精度。

圖9 索力識別結果示意圖Fig.9 Schematic diagram of cable force identification

4 索力識別影響因素分析

4.1 拉索抗彎剛度

拉索中同時存在幾何剛度與自身的抗彎剛度，最早的頻率法通過弦理論推導頻率、索力的關系，未考慮拉索自身的抗彎剛度，因而使得頻率法識別索力存在較大的誤差，大量研究證實，用動力方法來進行索力識別，必須考慮拉索自身的抗彎剛度，但對于拉索這種特殊結構，抗彎剛度如何取值，目前依然沒有準確的計算公式[1]。

Shimada[21]在研究中發現，拉索的抗彎剛度通常取計算抗彎剛度的0.5倍；Geier等[22]認為應該取計算抗彎剛度的2/3；謝曉峰[23]在研究中通過最小化頻率實測值與計算值之間的誤差對抗彎剛度進行修正，認為拉索的抗彎剛度通常取計算抗彎剛度的0.3～0.4倍。

因而，如何準確地利用抗彎剛度進行索力識別仍存在一定的研究價值，本文方法理論上可以進行多參數識別，但對索力和拉索自身的抗彎剛度同時識別，很容易出現錯誤。其實，通過理論推導不難發現，隨著頻率的增大，拉索抗彎剛度的影響也越來越大，若用于索力識別的抗彎剛度取值大于拉索的實際剛度，可以推斷，識別的索力必然隨著頻率點的增大而出現遞減的形式。

4.2 數值算例

在第3節的數值算例中，加入下述條件：截面的計算慣性矩為Izz=2×10-6m4，截面的實際慣性矩為Izz=0.6×10-6m4，將實際慣性矩代入模型中計算拉索動力響應，再用計算慣性矩進行索力識別，識別結果如圖10所示。

圖10 索力識別結果示意圖Fig.10 Schematic diagram of cable force identification

從圖10可以發現，由于進行索力識別的抗彎剛度取值比實際剛度大，因而出現了斜率為負的索力識別線，這與推斷一致。其實，要想識別出正確的索力結果，只需要修正抗彎剛度即可，可采用如下方法：

1)取α=1計算抗彎剛度αEI，在各頻率點計算索力，進行線性擬合，其斜率為β；

2)再次取值0<α<1，計算α對β的靈敏度k；

3)用α=α-βk更新抗彎剛度；

4)代入更新的抗彎剛度重新擬合索力，得到斜率β；

5)當β小于預設值時，輸出索力結果。

如圖11所示，經過數次抗彎剛度修正后，索力識別結果已經趨于一條水平線，此時α=0.303，與理論值0.3僅相差1%，索力識別結果為3 027.6 kN，誤差為0.92%，具有相當高的精度。

圖11 更新抗彎剛度后識別結果Fig.11 Recognition result after updating bending rigidity

4.3 噪音干擾

針對實際采集的振動信號通常存在某些干擾，并且通常是與頻率相關。在計算的拉索頻域振動響應中選取某個頻率段，加入隨機觀測誤差，再進行索力識別，結果如圖12所示。

圖12 包含干擾頻段索力識別圖Fig.12 Including interference band cable force identification diagram

由于是通過測點在每個頻率點的響應來識別索力，因而可以得到大量的索力識別值，若存在某個頻率段的干擾，通過識別圖很容易區分干擾頻域段，可以剔除掉該頻率段，以大量正確頻率點識別結果作為索力識別值。

5 結論

推導了修正Timoshenko梁振動模型的頻散關系，對拉索中的波分量進行了討論，提出了一種新的基于子結構中彎曲波的索力識別方法。該方法利用拉索中的行波，通過最小二乘法擬合波分量系數，以擬合殘差最小為目標進行索力識別；也對采用該方法進行索力識別的影響因素進行了探討。得到以下結論：

1)Euler-Bernoulli梁理論在粗短梁或者高頻率段存在較大誤差，Timoshenko梁存在截止頻率是由于忽略了剪切變形引起的轉動慣量，修正Timoshenko梁模型綜合全面地考慮了各種影響，是相對更加完善的梁理論，且在較高頻段具有更高的精度。

2)修正Timoshenko梁模型的拉索橫向振動的解由4個波系構成，可歸類為近場波與行波，距離梁端一定距離或較高的頻段可不考慮近場波的影響。

3)通過ANSYS建立拉索振動模型，用模態疊加法求解了拉索的動力響應，選取拉索的一個子結構，利用3個測點的動力響應識別了拉索子段的索力。此方法只需要拉索截面參數以及3個測點的相對位置即可進行索力識別，與邊界條件無關，在理論上具有十分高的精度。

4)對抗彎剛度以及激勵干擾對索力識別的影響進行了分析，提出了解決辦法，取得了良好的索力識別結果，理論偏差均不超過1%。