施 政，顏全勝

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

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

施 政，顏全勝*

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

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

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

0 引言

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

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

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

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

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

1 索的特性

1.1索參數(shù)說明

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

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

(1)

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

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

(2)

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

(3)

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

2.1兩端固支索

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

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

(4)

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

(5)

(6)

(7)

(8)

ωn=2πfn

(9)

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

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

(10)

cosh(βl)≈sinh(βl)

(11)

2αβ=0

(12)

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

(13)

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

由式(5)得到

(14)

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

(15)

構(gòu)造2個(gè)關(guān)于αl的函數(shù)

h(αl)=tan(αl)

(16)

(17)

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

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

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

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

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

因此求解流程如下：

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

(18)

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

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

由式(5)得

(19)

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

(20)

構(gòu)造2個(gè)關(guān)于αl的函數(shù)

h(αl)=tan(αl)

(21)

(22)

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

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

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

(23)

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

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

(24)

(25)

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

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

因此計(jì)算步驟為：

1)由式(25)計(jì)算nA；

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

(26)

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

2.2一端簡支一端固支索

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

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

(27)

(28)

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

h(αl)=tan(αl)

(29)

(30)

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

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

(31)

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

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

(32)

(33)

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

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

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

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

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

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

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

2)根據(jù)計(jì)算得到擬合直線斜率值k用二分法調(diào)整EIfac上下限：由于索力只有一個(gè)值，因此取得正確EIfac時(shí)，各階Tn應(yīng)該相同，即k=0。當(dāng)k>0說明EI偏小，則令EImin=EIfac；當(dāng)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，拉索應(yīng)力小于容許應(yīng)力σ<1.86E9Pa(1 860 MPa)。

4.2編制程序

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

4.3頻率計(jì)算

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

4.4索力和抗彎剛度識別

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

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

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

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

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

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

5 工程拉索測試

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

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

圖6 范和港斜拉橋

圖7 拾取振動(dòng)信號

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

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

6 結(jié)論

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

2)本文方法無需輸入索力預(yù)估值，即可識別索力，避免了索力預(yù)估值錯(cuò)誤導(dǎo)致索力識別錯(cuò)誤的潛在風(fēng)險(xiǎn)。

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

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

[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ì)算拉索拉力的實(shí)用公式[J].土木工程學(xué)報(bào),2005,38(11):26-31.

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

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

[5]魏金波.彈性支承鋼索受力動(dòng)力檢測理論和試驗(yàn)研究[D].上海:同濟(jì)大學(xué),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].杭州:浙江大學(xué),2010.

[11]魏金波,段欣.考慮抗彎剛度的索力動(dòng)力檢測法[J].陜西理工學(xué)院學(xué)報(bào)(自然科學(xué)版),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ù)識別研究。

國家自然科學(xué)基金項(xiàng)目(11202080)；廣東省交通運(yùn)輸廳科技項(xiàng)目(科技-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。