基于子結(jié)構(gòu)中彎曲波的拉索索力識別方法

肖軍，代洋，張永水

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

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

1 拉索振動頻散關(guān)系

1.1 修正Timoshenko梁理論

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

拉索微元段引入剪切變形所引起的轉(zhuǎn)動慣量后，平衡狀態(tài)如圖1所示。

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

根據(jù)圖1建立拉索微段平衡方程組

(1)

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

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

(2)

(3)

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

(4)

(5)

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

(6)

(7)

(8)

并且，

(9)

經(jīng)Fourier變換后可得

(10)

(11)

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

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

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

(12)

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

(13)

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

1.2 數(shù)值算例

基于梁理論的推導(dǎo)，結(jié)合一數(shù)值算例進(jìn)一步探討修正Timoshenko梁、Euler-Bernoulli梁與Timoshenko梁理論的關(guān)系與區(qū)別。

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

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

圖2 波數(shù)-頻率圖Fig.2 Wave number-frequence diagram

2 拉索動力響應(yīng)的波分量分解

2.1 波分量理論

式(13)給出了修正Timoshenko梁的振動頻散關(guān)系，根據(jù)疊加原理可以得到忽略垂度的拉索頻域橫向自由振動通解(后文波數(shù)解均依賴修正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負(fù)方向衰減；

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

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

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

H·C=0

(15)

式中

(16)

H=

(17)

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

(18)

2.2 數(shù)值算例

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

圖3 結(jié)構(gòu)示意圖Fig.3 structure diagram

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

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

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

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

(19)

為近場波分量；

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

(20)

為行波分量。

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

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

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

3 基于子結(jié)構(gòu)彎曲波的索力識別方法

3.1 識別方法理論

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

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

(21)

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

(22)

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

(23)

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

(24)

為系數(shù)矩陣，取決于結(jié)構(gòu)特征與外部激勵；

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

如果拉索的參數(shù)、測點布置以及各測點的響應(yīng)結(jié)果已經(jīng)得到，可以通過最小二乘法對波分量系數(shù)進(jìn)行求解。

(25)

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

(26)

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

P={N,EI}=

(27)

3.2 索力識別數(shù)值算例

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

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

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

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

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

圖8 頻域響應(yīng)示意圖Fig.8 Frequency domain response diagram

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

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

4 索力識別影響因素分析

4.1 拉索抗彎剛度

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

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

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

4.2 數(shù)值算例

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

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

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

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

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

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

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

5)當(dāng)β小于預(yù)設(shè)值時，輸出索力結(jié)果。

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

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

4.3 噪音干擾

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

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

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

5 結(jié)論

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

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

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

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

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