基于二分法識別拉索索力與抗彎剛度

施 政，顏全勝

(華南理工大學土木與交通學院，510640，廣州)

基于二分法識別拉索索力與抗彎剛度

施 政，顏全勝*

(華南理工大學土木與交通學院，510640，廣州)

采用二分法對拉索的索力與抗彎剛度進行同時識別，來滿足橋梁工程中快速、準確的識別索力與抗彎剛度的要求。通過將頻率方程進行數(shù)學上的處理和簡化，成為簡化的頻率方程，并考慮初等函數(shù)曲線的性質(zhì)，采用二分法在指定區(qū)間內(nèi)迭代的方法，求解該簡化方程的根，可以實現(xiàn)由索力求解任意階的頻率，或由至少2階(頻階任意)的頻率識別索力和抗彎剛度。在Excel中建立數(shù)值拉索，利用其VBA平臺，編程實現(xiàn)了該方法的自動計算；與ANSYS結果對比，表明精度滿足工程要求。

頻率法；二分法；索力；抗彎剛度；參數(shù)識別；VBA

0 引言

頻率法測得各階頻率之后，如何識別未知的索力和抗彎剛度，是索力測試中最主要問題之一。索力的計算是基于抗彎剛度值的，因此精確的識別抗彎剛度，保證索力識別的準確性，是十分重要的。由于索的內(nèi)部平行鋼絲束之間相互擠壓而并非完全粘結，索彎曲過程中由于相對位移而產(chǎn)生摩擦，索外包聚乙烯保護層和索內(nèi)部其他構件等均有抗彎剛度等諸多原因，導致索的抗彎剛度是介于零和同截面特性的梁的抗彎剛度之間的數(shù)，無法精確確定。

對于物理剛度較大，兩端固支的短索，特征參數(shù)ξ較小的索而言，抗彎剛度的不精準，將導致較大的索力計算誤差。

識別抗彎剛度EI的前提是能夠計算對應于各階頻率fn的索力Tn。許多學者已研究并推導了由頻率求索力的經(jīng)驗公式，但其中多數(shù)不能由高階頻率計算索力，且有其特定的ξ適用范圍，有的對于ξ大的索識別準確，有的對于ξ小的索識別準確。

Hiroshi Zui[1]等通過繪出兩端固支索頻率方程的解曲線并進行曲線擬合，得到精度較高的經(jīng)驗公式，但公式基于1階和2階頻率求索力，采用高階頻率計算索力需要ξ>200；任偉新、陳剛[2]基于能量法與曲線擬合法得出由基頻計算索力的經(jīng)驗公式，只能通過基頻求索力；邵旭東[3]基于能量法和索的一階振型函數(shù)得到由剛性吊桿基頻計算索力的公式，但該公式的準確性與參數(shù)G有關，而G隨拉索參數(shù)不同而變化，引起精度的不穩(wěn)定；孟少平[4]通過一、二階頻率計算索力，避開了抗彎剛度的問題，但該式只能用第一、二階頻率計算索力，且不能識別抗彎剛度。魏金波[5]提出了對頻率方程進行簡化的思路，蘇成[7]提出了同時識別索力和抗彎剛度的方法，但使用的是3次樣條擬合技術且需要進行有限元分析，且試算索力與抗彎剛度均采用定步長法，其搜索效率低于二分法。

本文綜合文獻[5,7]的思路和方法，采用二分法[8]直接對非線性方程進行迭代求根，并利用Excel表格中的VBA功能編制程序，實現(xiàn)方便快捷的進行索力與抗彎剛度的識別。

1 索的特性

1.1索參數(shù)說明

本文中用到的主要參數(shù)包括：T=索力；EI=索的抗彎剛度；EIfac=抗彎剛度折減系數(shù)；l=索的長度；m=索的單位長度質(zhì)量；n=頻階；ωn=第n階圓頻率；fn=第n階頻率；Tn=由頻率fn計算的第n階索力；αl=文中構造函數(shù)的自變量；ξ=無量綱參數(shù)。

工程中的索各式各樣，索的特性各不相同，但頻率-索力曲線的樣式不同，其本質(zhì)在于特征參數(shù)ξ不同[6]。

(1)

2 不同邊界索的頻率與索力求解

本文研究的索為小垂度索在xy平面內(nèi)的振動。索的頻率方程寫為式(2)，其中y(x,t)為索的橫向位移坐標，x為沿索縱向坐標，φ(x)為振型函數(shù)，h(t)為索力隨索振動時的變化量，這里忽略它取h(t)=0，其他參數(shù)見1.1節(jié)。

(2)

φ(x)=C1sinαx+C2cosαx+C3sinhβx+C4coshβx

(3)

根據(jù)不同的邊界條件求出系數(shù)C1～C4，得到頻率方程。以下重點研究兩端固支邊界和一端簡支一端固支邊界的拉索的索力和抗彎剛度識別。

2.1兩端固支索

兩端固支索的頻率方程[1]如式(2)。

2(αl)(βl)[1-cos(αl)cosh(βl)]+[(βl)2-(αl)2]sin(αl)sinh(βl)=0

(4)

各參數(shù)表達式如式(3)。

(5)

(6)

(7)

(8)

ωn=2πfn

(9)

該方程為隱式超越方程，由于雙曲正弦、余弦函數(shù)帶來的強烈的非線性。直接使用N.R迭代法求解，如果給定的迭代初值x0不足夠靠近零點，初值處的切線斜率易趨于無窮大f′(x0)→∞，于是之后一系列迭代點xn始終在初值x0附近不下降，切線斜率始終很大，如此惡性循環(huán)以致收斂速度極為緩慢。然而，若采用不需要導數(shù)的二分法進行迭代即可避免上述問題。

根據(jù)文獻[5]的簡化法，當ξ較大時，結合雙曲正弦、余弦函數(shù)的性質(zhì)，和βl的表達式，可以近似認為式(10)～式(12)的假定成立。

(10)

cosh(βl)≈sinh(βl)

(11)

2αβ=0

(12)

于是，在式(10)～式(12)假定下，頻率方程可簡化為式(13)，并寫成等號一邊為正弦函數(shù)，另一邊為關于αl表達式的形式

(13)

2.1.1 求解頻率 已知索力T、抗彎剛度EI，求第n階頻率fn。

由式(5)得到

(14)

將式(14)代入式(13)化簡得到

(15)

構造2個關于αl的函數(shù)

h(αl)=tan(αl)

(16)

(17)

則方程(15)的求根的問題轉化為求函數(shù)h(αl)和g(αl)的交點問題，兩函數(shù)圖像見圖1。

圖1 h(αl)&g(αl)～αl圖像兩端固支索求頻率

根據(jù)表達式(17)，αl≥0時，g(αl)是αl=0

處略低于拋物線，而隨著αl的增大不斷逼近函數(shù)p(αl)的近似二次函數(shù)，在(0,+∞)上遞增。

由式(14)，αl與fn增減性相同，fn隨n增大而增大，且在αl≥0時，由圖1，隨著αl增大，兩函數(shù)的第1個交點在n=1的區(qū)間內(nèi)，n>1的每個區(qū)間內(nèi)有且僅有一個交點，可推得第n階頻率對應于第n個區(qū)間內(nèi)的交點。

因此求解流程如下：

f(αl)=h(αl)-g(αl)

(18)

2)將αl結果與已知索力T代入式(14)求fn。

2.1.2 求解索力 此時已知索的第n階頻率fn和抗彎剛度EI，要求對應的索力Tn。

由式(5)得

(19)

將式(19)代入式(13)化簡后得

(20)

構造2個關于αl的函數(shù)

h(αl)=tan(αl)

(21)

(22)

則求解方程(20)的根的問題轉化為求函數(shù)h(αl)和g(αl)的交點問題，圖2為兩函數(shù)的圖像。

圖2 h(αl)&g(αl)～αl圖像兩端固支索求索力

對于g(αl)，令分母為零時的αl取值為αlA，則有

(23)

即g(αl)在αlA處不連續(xù)。因αl≥0，對g(αl)求導易得g′(αl)>0，所以g(αl)為關于αl的在(0,αlA),(αlA,+∞)上的遞增函數(shù)，且在區(qū)間(0,αlA)上g(αl)>0，在區(qū)間(αlA,+∞)上g(αl)<0。

1)推導間斷點所在區(qū)間的區(qū)間號nA，它需要滿足式(24)，解出nA表達式(25)，Int表示取整；

(24)

(25)

2)由于兩函數(shù)在n=nA的區(qū)間無交點，在n>nA的每個區(qū)間只有一個交點，因此當n≥nA時，需要前移一個區(qū)間即令n=n+1，再求解。

3)根據(jù)式(29)，當αl>αlA時，會有Tn<0，其物理意義是，因為目標頻率太小，以至于當索受到軸向壓力時，頻率才能達到如此小。拉索是受拉的所以索力不能取負值。在迭代求EI過程中，一旦出現(xiàn)Tn<0，則EI肯定偏大，無需再計算，可直接二分法減小EI識別值，再重新迭代。

因此計算步驟為：

1)由式(25)計算nA；

f(αl)=h(αl)-g(αl)

(26)

3)將αl結果與已知頻率fn代入式(19)求Tn。

2.2一端簡支一端固支索

修改邊界條件，推導出頻率方程如式(27)，該式變形為式(28)，其中參數(shù)的表達式見式(5)～式(9)。

βlsin(αl)cosh(βl)-αlcos(αl)sinh(βl)=0

(27)

(28)

類比2.2節(jié)的方法，構造兩函數(shù)

h(αl)=tan(αl)

(29)

(30)

圖3h(αl)&g(αl)～αl圖像一端固支一端鉸支求頻率

2.2.1 求解頻率 將式(10)代入式(30)得到式(31)

(31)

2)將得到αl代入式(14)計算頻率fn。

2.2.2 求解索力 由式(19)代入式(1)變形為式(32)，代入式(30)得

(32)

(33)

2)將1)中得到的αl值和已知頻率fn代入式(19)計算索力Tn。

圖4h(αl)&g(αl)～αl圖像一端固支一端鉸支求索力

3 二分法迭代的參數(shù)識別流程

二分法迭代的優(yōu)點是迭代次數(shù)較少、穩(wěn)定、收斂性較好，能夠保證迭代終值的精度，迭代次數(shù)由精度決定，達到1E-6精度平均迭代次數(shù)為log2106=19.93。因此，采用二分法迭代對簡化方程式(14)求根。識別索力與抗彎剛度的算法流程圖如圖5。

算法包括兩重循環(huán)：1)外循環(huán)，由fn求Tn；2)內(nèi)循環(huán)，由Tn求EI。

算法關鍵的二分法調(diào)整步驟用偽代碼描述：

1)EIfac的上下限初值：EImin=0,EImax=0。

2)根據(jù)計算得到擬合直線斜率值k用二分法調(diào)整EIfac上下限：由于索力只有一個值，因此取得正確EIfac時，各階Tn應該相同，即k=0。當k>0說明EI偏小，則令EImin=EIfac；當k<0說明EI偏大，則令EImax=EIfac。

3)αl的上下限初值：

4)根據(jù)f(αl)調(diào)整上下限：如果f(αl)>0，說明αl偏大，則令almax=αl；如果f(αl)<0，說明αl偏小，則令almin=αl。

5)為減少誤差，T取為各階Tn的平均值。

4 數(shù)值拉索測試

4.1建立拉索

圖5 二分法參數(shù)識別流程圖

建立9根數(shù)值拉索如表1所示，滿足ξ=1～500，邊界為固支和一端固支一端鉸支；小垂度索需滿足垂跨比δ=s/l0?1，注意控制等效彈模Eeq非常接近于彈模E=1.8E11，拉索應力小于容許應力σ<1.86E9Pa(1 860 MPa)。

4.2編制程序

在Excel的VBA平臺用Visual Basic語言編程識別邊界為固支、一端固支一端鉸支拉索的頻率計算程序，和索力、抗彎剛度的識別程序。

4.3頻率計算

在ANSYS中，對已經(jīng)建立的數(shù)值索，進行模態(tài)分析，得到模態(tài)頻率(由有限元計算)；在Excel中，根據(jù)輸入的拉索參數(shù)(已知索力與抗彎剛度)，計算各階頻率；對每一根索，將在ANSYS中計算的模態(tài)頻率與Excel中計算的頻率比較，如表2所示。

4.4索力和抗彎剛度識別

將ANSYS計算得到的n階頻率作為實測頻率，采用Excel中程序計算，得到索力與抗彎剛度的識別值，將該識別值與建模時的真實值比較，結果如表3所示。

分析表2中數(shù)據(jù)得到結論：

1)由表2知，本文方法的頻率計算誤差與ANSYS結果相比小于0.3%，精度滿足工程要求；

2)由表2知，索力識別誤差小于1%，抗彎剛度識別誤差小于0.3%，精度滿足工程要求。

表1 數(shù)值索的各項參數(shù)

表2 Excel與ANSYS的1～3階頻率結果比較、索力與抗彎剛度的識別結果

5 工程拉索測試

范和港斜拉橋位于廣東省惠州市惠東縣巽寮灣，為主跨300 m的雙塔單索面斜拉橋，兩塔的中跨和邊跨均設23對斜拉索，索的特征參數(shù)在之間，索力設計值在1 986～4 038 KN之間。

采用頻率法識別索力，取短索、中長索、長索各1根，邊界為兩端固支，分別用Zui公式[1](帶入抗彎剛度折減系數(shù)0.6)以及本文公式，根據(jù)實測頻率計算索力、抗彎剛度折減系數(shù)，并對比列于表3。

圖6 范和港斜拉橋

圖7 拾取振動信號

由表3知，本文的索力計算結果與Zui公式索力計算結果的誤差在0.1%以內(nèi)，識別得到的各索的抗彎剛度均在0.6EI0左右，且其誤差隨著索的特征參數(shù)的增大而增大，說明本文方法的計算結果較為精確，進一步說明識別抗彎剛度折減系數(shù)時，要使用特征參數(shù)小的索才能得到較高精度。這是由于特征參數(shù)小的索的幾何剛度所占比例偏小，物理剛度所占比例偏大，因此其頻率對抗彎剛度的靈敏度較高，識別較為精確；特征參數(shù)大的索的幾何剛度所占比例偏大，物理剛度所占比例偏小，因此其頻率對抗彎剛度的靈敏度很低，識別精度偏低。

表3 模態(tài)頻率計算結果對比

6 結論

1)本文方法只適用于求解兩端固支或鉸支的小垂度拉索，理論上可求解任意指定的頻階n的索力Tn和頻率fn，與有限元索力和頻率結果的誤差普遍小于1%。

2)本文方法無需輸入索力預估值，即可識別索力，避免了索力預估值錯誤導致索力識別錯誤的潛在風險。

3)輸入的實測頻率，頻階可以是任意的(至少2個，互不重復)，優(yōu)于文獻[2]、[3]中只能根據(jù)基頻或二階頻率進行索力識別。

4)采用二分法迭代的收斂速度快且穩(wěn)定，計算精度高于多數(shù)擬合公式的精度，且適用于特征參數(shù)ξ>2的拉索，適用范圍廣。數(shù)據(jù)的輸入和輸出過程集成在Excel環(huán)境中，拉索參數(shù)，實測頻率與索力、抗彎剛度計算結果方便于查看與修改，工程實用性強。

[1] Hiroshi Zui,Tohru Shinke,Yoshio Namita.Practical Formulas for Estimation of Cable Tension by Vibration Method[J].Joumal of Structural Engineering,ASCE,1996,122(6):651-656.

[2]任偉新,陳剛.由基頻計算拉索拉力的實用公式[J].土木工程學報,2005,38(11):26-31.

[3]邵旭東, 李國峰, 李立峰. 吊桿振動分析與力的測量[J].中外公路,2004,24(6):29-31.

[4]孟少平,楊睿,王景全.一類精確考慮抗彎剛度影響的系桿拱橋索力測量新公式[J].公路交通科技,2008,25(6):87-98.

[5]魏金波.彈性支承鋼索受力動力檢測理論和試驗研究[D].上海:同濟大學,2009.

[6]Irvine H M.Cable Structures[M].Cambridge,MA:The MIT Press,1981.

[7]蘇成,徐郁峰,韓大建.頻率法測量索力中的參數(shù)分析與索抗彎剛度的識別[J].公路交通科技,2005,22(5):75-78.

[8]Li F C,Zhang L N,Tian S Z,etal.Analysis of Identified Cable Force of Cable-Stayed Bridge[C].International Conference on Transportation Engineering,2009:3020-3025.

[9]Song Y,Cui X P,Lei Y,etal.Intelligent Sensors with Application to Identify the Steel Cable Forces[C].Earth and Space,ASCE:2010,1626-1633.

[10]朱鐵誠.拉索參數(shù)識別的多頻率法研究與索力無線監(jiān)測系統(tǒng)開發(fā)[D].杭州:浙江大學,2010.

[11]魏金波,段欣.考慮抗彎剛度的索力動力檢測法[J].陜西理工學院學報(自然科學版),2009,25(3):30-33.

[12]Ren W X,Hu W H.Cable modal parameters identification[J].Journal of Engineering Mechanics,2009,135,(1):51-61.

IdentificationofCableTensionandFlexuralRigiditybyDichotomyMethod

SHI Zheng,YAN Quansheng*

(School of Civil Engineering and Transportation,South China University of Technology,510640,Guangzhou,PRC)

Dichotomy method is used in cable parameter identification,in order to satisfy the need of identifying cable tension and flexural rigidity quickly and accurately.A simplified frequency equation can be acquired by adopting mathematical treatment and simplification to the frequency equation.Taking the property of elementary function curve into consideration,dichotomy method is adopted to find the solution of simplified frequency equation in specified interval according to frequency mode,to realize the identification of cable tension and flexural rigidity according to at least 2(or more) arbitrary modes of frequency.Cable frequency of arbitrary mode can also be calculated according to given cable tension and flexural rigidity.Numeric cables are created in Excel and a program is developed base on VBA platform to implement the parameter identification automatically.The calculation result compared to ANSYS result,shows that the identification error satisfies the engineering demand.

vibration method;dichotomy method;cable tension;flexural rigidity;parameter identification;VBA

2014-08-26;

2014-09-28

施 政(1991-)，男，江西南昌人，碩士研究生，研究方向：拉索參數(shù)識別研究。

國家自然科學基金項目(11202080)；廣東省交通運輸廳科技項目(科技-2012-02-024)。

10.13990/j.issn1001-3679.2014.05.020

U433.38

A

1001-3679(2014)05-0667-07

*通訊作者：顏全勝(1968-), 男, 博士, 教授; 主要從事大跨度橋梁研究， E-mail:cvqshyan@scut.edu.cn。