基于子结构中弯曲波的拉索索力识别方法

肖军，代洋，张永水

(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%。