薄壁箱梁具有良好的抗弯和抗扭性能在现代桥梁中得到广泛应用.当箱梁发生竖向位移时，由于翼板中剪力滞后的影响引起翼板纵向应力沿横向分布不均匀，而存在剪力滞效应^[1-5].

对于箱梁自由振动而言，剪力滞效应引起箱梁的动力特性发生显著变化.对于不同的跨高比，可忽略箱梁部分影响因素而建立相应的自由振动方程.当跨高比较大时，甘亚南等^[6]、吴有俊等^[7]以欧拉梁为基础，忽略了剪力滞引起纵向位移的惯性影响，分别研究了剪力滞引起的弹性势能对梁动力特性的影响.对于跨高比小的箱梁，此时截面的剪切变形和转动惯量影响不可忽略.张永健等^[8]分析了剪切变形和剪力滞效应对简支箱梁自振频率的影响.甘亚南等^[9-10]基于Timoshenko梁理论讨论了剪力滞效应对等截面箱梁自振频率的影响.周旺保等^[11]研究了剪力滞、剪切变形、转动惯量、滑移效应对钢-混凝土组合箱型梁动力特性的影响.

在上述箱梁自由振动的求解过程中，大部分是基于分离变量法得到箱梁自由振动的解析解.由于箱梁的自由振动控制方程是一个微分方程组，因此自振频率超越方程及模态函数很复杂.从梁的振动方程角度看，剪力滞效应对梁振动的影响可看作是欧拉梁振动方程修正后形成的系统，则可以利用摄动法求解箱梁自由振动.楼梦麟^[12]首先利用模态摄动法，将变系数微分方程的求解问题转化为代数方程组进行求解，从而简化了特征方程的求解.随后，楼梦麟等^[13]、潘旦光等^[14]将模态摄动法应用于变截面梁的振动.潘旦光等^[15-16]进一步将模态摄动法推广到变截面Timoshenko梁振动方程的求解.本文将利用模态摄动法的基本思想，研究箱梁的自振频率和振型的简化分析方法.在此基础上，根据箱梁位移和弯矩的关系，推导了箱梁自由振动的模态剪力滞系数，并分析跨宽比和翼缘板的刚度占梁总刚度百分比等对梁模态剪力滞系数的影响.

1 箱梁的自由振动方程

图 1所示的矩形薄壁箱型梁，若梁的竖向位移为w (x, t)，上下翼缘板的纵向位移函数为v_i(x, y, t)，且假定v_i(x, y, t)可表示为

$ {v_i}\left( {x, y, t} \right) = {h_i}\left[{\frac{{{\rm{d}}w}}{{{\rm{d}}x}} + \left( {1-\frac{{{y^3}}}{{{b^3}}}u\left( {x, t} \right)} \right)} \right]. $

(1)

式中：u(x, t)为翼缘板剪切转角的最大差值，x、y和z分别表示顺梁方向、垂直于梁方向和竖向，b为箱室净宽的一半，h_i截面形心到顶、底板的距离，i分别取顶、底板.

基于欧拉梁理论，梁的动能T为

$ T = \frac{1}{2}{\int\limits_0^l {\rho A\left( {\frac{{\partial w}}{{\partial t}}} \right)} ^2}{\rm{d}}x + \frac{1}{2}{\int\limits_0^l {\rho {I_s}\left( {\frac{{\partial u}}{{\partial t}}} \right)} ^2}{\rm{d}}x. $

(2)

梁的势能V为

$ \begin{array}{l} V = \frac{1}{2}\int\limits_0^l {EIw'{'^2}} {\rm{d}}x + \frac{1}{2}\int\limits_0^l {\frac{9}{5}\frac{G}{{{b^2}}}{I_s}{u^2}} {\rm{d}}x + \\ \;\;\;\;\;\frac{1}{2}\int\limits_0^l {E{I_s}\left( {\frac{3}{2}w''u' + \frac{9}{{14}}u{'^2}} \right)} {\rm{d}}x. \end{array} $

(3)

式中：ρ、E、G、A分别为箱梁的密度、弹性模量、剪切模量和截面面积.I为截面的转动惯量，I_s为顶板和底板的转动惯量.根据Hamilton原理可得梁的自由振动方程为：

$ \rho A\ddot w + EI{w^{\left( 4 \right)}} + \frac{3}{4}E{I_s}{u^{\left( 3 \right)}} = 0, $

(4a)

$ \rho {I_s}\ddot u-\frac{3}{4}E{I_s}{w^{\left( 3 \right)}}-\frac{9}{{14}}E{I_s}u'' + \frac{9}{5}\frac{{G{I_s}}}{{{b^2}}}u = 0. $

(4b)

2 简支箱梁动力特性的模态摄动解

若式(4)所对应的第j阶模态的特征值为${{\bar \lambda }_j} = \bar \omega _j^2$，特征向量为${{\bar W}_j}\left( x \right)$和${{\bar U}_j}\left( x \right)$，则第j阶模态的特征方程为：

$ -{{\bar \lambda }_j}\rho A{{\bar W}_j} + EI{{\bar W}_j}^{\left( 4 \right)} + \frac{3}{4}E{I_s}{{\bar U}_j}^{\left( 3 \right)} = 0, $

(5a)

$ {{\bar \lambda }_j}\rho {I_s}{{\bar U}_j}' + \frac{3}{4}EI{{\bar W}_j}^{\left( 4 \right)} + \frac{9}{{14}}E{I_s}{{\bar U}_j}^{\left( 3 \right)}-\frac{{9G{I_s}}}{{5{b^2}}}{{\bar U}_j}' = 0. $

(5b)

忽略剪力滞效应后，具有相同跨度欧拉梁第j阶模态的控制方程为

$ -{\lambda _j}\rho A{W_j} + EI{W_j}^{\left( 4 \right)} = 0. $

(6)

当梁为简支梁时，与式(6)特征方程相对应的特征值和模态为

$ {\lambda _j} = {\left( {\frac{{j{\rm{\pi }}}}{l}} \right)^4}\frac{{EI}}{{\rho A}}, \;\;{W_j} = \sin \left( {\frac{{j{\rm{\pi }}}}{l}x} \right). $

(7)

直接模态摄动法的基本思想是把式(5)考虑剪力滞效应薄壁箱梁的特征方程看成式(6)所表示的等截面欧拉梁经过参数修改所得到的新系统，这个新系统的主模态函数以及特征值可以用等截面欧拉梁的模态特征经过简单的摄动分析求解.则可以假设：

$ {{\bar \lambda }_j} = {\lambda _j} + \Delta {\lambda _j}, $

(8a)

$ {{\bar W}_j} = {W_j} + \sum\limits_{k = 1, k \ne j}^n {{W_k}{p_{kj}}, } $

(8b)

$ {{\bar U}_j} = \sum\limits_{k = 1}^n {{W_k}'{q_{kj}}} . $

(8c)

从理论上来说，式(6)中梁有无穷多个主模态函数，即式(8)中n应该趋向于无穷.但是在实际计算时，通常只需要考虑有限个低阶模态进行近似的计算就可以满足要求.当解出Δλ_j，p_kj，q_kj这2n个未知数，即可计算箱梁的第j阶特征值及其对应的主模态函数.将(8)式代入(5)式可得：

$ \begin{array}{l} -\left( {{\lambda _j} + \Delta {\lambda _j}} \right)\rho A\left( {{W_j} + \sum\limits_{k = 1, k \ne j}^n {{W_k}{p_{kj}}} } \right) + EI\left( {{W_j}^{\left( 4 \right)} + } \right.\\ \left. {\sum\limits_{k = 1, k \ne j}^n {{W_k}^{\left( 4 \right)}{p_{kj}}} } \right) + \frac{3}{4}E{I_s}\sum\limits_{k = 1}^n {{W_k}^{\left( 4 \right)}{q_{kj}}} = 0, \end{array} $

(9a)

$ \begin{array}{c} \left( {{\lambda _j} + \Delta {\lambda _j}} \right)\rho {I_s}\sum\limits_{k = 1}^n {{W_k}''{q_{kj}}} + \frac{3}{4}E{I_s}\left( {{W_j}^{\left( 4 \right)} + \sum\limits_{k = 1, k \ne j}^n {{W_k}^{\left( 4 \right)}{p_{kj}}} } \right) + \\ \frac{9}{{14}}E{I_s}\sum\limits_{k = 1}^n {{W_k}^{\left( 4 \right)}{q_{kj}}}-\frac{{9G{I_s}}}{{5{b^2}}}\sum\limits_{k = 1}^n {{W_k}''{q_{kj}}} = 0. \end{array} $

(9b)

在方程(9)两边同时乘以W_i(x)(i=1, 2, ···n)，然后沿全长积分，根据欧拉梁的模态正交性，化简后可得：

$ \begin{array}{*{20}{l}} {-\left( {{\lambda _j} + \Delta {\lambda _j}} \right)\rho A\sum\limits_{k = 1, k \ne j}^n {{R_{ik}}{p_{kj}}} + EI\sum\limits_{k = 1, k \ne j}^n {{T_{ik}}{p_{kj}}-} }\\ {\Delta {\lambda _j}\rho A{R_{ij}} + \frac{3}{4}E{I_s}\sum\limits_{k = 1}^n {{T_{ik}}{q_{kj}}} = {\lambda _j}\rho A{R_{ij}}-EI{T_{ij}}, } \end{array} $

(10a)

$ \begin{array}{l} \frac{3}{4}E{I_s}\sum\limits_{k = 1, k \ne j}^n {{T_{ik}}{p_{kj}}} + \rho {I_s}\left( {{\lambda _j} + \Delta {\lambda _j}} \right)\sum\limits_{k = 1}^n {{S_{ik}}{q_{kj}}} + \\ \frac{9}{{14}}E{I_s}\sum\limits_{k = 1}^n {{T_{ik}}{q_{kj}}}-\frac{{9G{I_s}}}{{5{b^2}}}\sum\limits_{k = 1}^n {{S_{ik}}{q_{kj}}} =-\frac{3}{4}E{I_s}{T_{ij}}. \end{array} $

(10b)

式(10)中：

$ \begin{array}{c} {R_{ij}} = \int\limits_0^l {{W_i}{W_j}{\rm{d}}x} = \frac{l}{2}{\delta _{ij}}, \\ {S_{ij}} = \int\limits_0^l {{W_i}{W_j}{\rm{''d}}x} =-{\left( {\frac{{j{\rm{\pi }}}}{l}} \right)^2}\frac{1}{2}{\delta _{ij}}, \\ {T_{ij}} = \int\limits_0^l {{W_i}{W_j}^{\left( 4 \right)}{\rm{d}}x} = {\left( {\frac{{j{\rm{\pi }}}}{l}} \right)^4}\frac{1}{2}{\delta _{ij}}. \end{array} $

式中δ_ij为Kronecker符号.

依次取i为1, 2, …，n，重复利用式(10)可得2n个代数方程.将2n个代数方程写成矩阵形式：

$ \left( {\left[\begin{array}{l} {\mathit{\boldsymbol{C}}^{11}}{\mathit{\boldsymbol{C}}^{12}}\\ {\mathit{\boldsymbol{C}}^{21}}{\mathit{\boldsymbol{C}}^{22}} \end{array} \right] + \Delta {\lambda _j}\left[\begin{array}{l} {\mathit{\boldsymbol{D}}^{11}}0\\ 0{\mathit{\boldsymbol{D}}^{22}} \end{array} \right]} \right)\left\{ \begin{array}{l} \mathit{\boldsymbol{p}}\\ \mathit{\boldsymbol{q}} \end{array} \right\} = \left\{ \begin{array}{l} \mathit{\boldsymbol{M}}\\ \mathit{\boldsymbol{N}} \end{array} \right\}. $

(11)

式中：C¹¹、C¹²、C²¹、C²²、D¹¹、D²²都为n阶方阵，p、q、M、N为n阶向量.各个方阵和向量的元素分别为：

$ \begin{array}{c} \mathit{\boldsymbol{C}}_{ik}^{11} = \left\{ {\begin{array}{*{20}{c}} {-{\lambda _j}\rho A{R_{ik}} + EI{T_{ik}}, }&{k \ne j;}\\ {-{\lambda _j}\rho A{R_{ik}}, }&{k = j.} \end{array}} \right.\\ \mathit{\boldsymbol{C}}_{ik}^{12} = \frac{3}{4}E{I_s}{T_{ik}}.\\ \mathit{\boldsymbol{C}}_{ik}^{21} = \left\{ {\begin{array}{*{20}{c}} {\frac{3}{4}E{I_s}{T_{ik}}, }&{k \ne j;}\\ {0, }&{k = j.} \end{array}} \right.\\ \mathit{\boldsymbol{C}}_{ik}^{22} = {\lambda _j}\rho {I_s}{S_{ik}} + \frac{9}{{14}}E{I_s}{T_{ik}}-\frac{{9G{I_s}}}{{5{b^2}}}{S_{ik}}, \\ \mathit{\boldsymbol{D}}_{ik}^{11} = \left\{ {\begin{array}{*{20}{c}} { - \rho A{R_{ik}}}&{k \ne j;}\\ {0, }&{k = j.} \end{array}} \right.\\ \mathit{\boldsymbol{D}}_{ik}^{22} = \rho {I_s}{S_{ik}}.\\ \mathit{\boldsymbol{p}} = {\left\{ {{p_{1j}} \cdots {p_{\left( {j - 1} \right)j}}\;\;\Delta {\lambda _j}/{\lambda _j}\;\;{p_{\left( {j + 1} \right)j}} \cdots {p_{nj}}} \right\}^{\rm{T}}}, \\ \mathit{\boldsymbol{q}} = {\left\{ {{q_{1j}}\;\;{q_{2j}}\;\; \cdots \;\;{q_{nj}}} \right\}^{\rm{T}}}, \\ {\mathit{\boldsymbol{M}}_i} = - EI{T_{ij}} + {\lambda _j}\rho A{R_{ij}}, \;\;{\mathit{\boldsymbol{N}}_i} = - \frac{3}{4}E{I_s}{T_{ij}}. \end{array} $

显然C¹¹、C¹²、C²¹、C²²、D¹¹、D²²都为对角矩阵，M_i=0，N_i=0(i≠j)，因此，p、q只有第j个元素不为零，则式(11)只剩下两个方程.由此可得：

$ {{\bar \lambda }_j}-\frac{3}{4}\frac{{{I_s}}}{I}{q_{jj}}{\lambda _j} = {\lambda _j}, $

(12a)

$ {q_{jj}} = \frac{{-\frac{3}{4}E{I_s}{{\left( {\frac{{j{\rm{\pi }}}}{l}} \right)}^2}}}{{\frac{9}{{14}}E{I_s}{{\left( {\frac{{j{\rm{\pi }}}}{l}} \right)}^2} + \frac{{9G{I_s}}}{{5{b^2}}}-\rho {I_s}{{\bar \lambda }_j}}}. $

(12b)

将式(12b)代入式(12a)可得

$ \frac{{{{\bar \lambda }_j}}}{{{\lambda _j}}} = 1-\frac{{\frac{{{j^2}{{\rm{\pi }}^2}}}{{{l^2}}}\left( {\alpha-1} \right)}}{{\alpha \frac{{{j^2}{{\rm{\pi }}^2}}}{{{l^2}}} + {\beta ^2}-\frac{{14}}{9}\frac{{\alpha \rho {{\bar \lambda }_j}}}{E}}} $

(13)

式中$\alpha = \frac{1}{{1-7{I_{\rm{s}}}/8I}}, \beta = \frac{1}{b}\sqrt {\frac{{14G\alpha }}{{5E}}} $为瑞纳斯常数.当I_s=0，α=1，即梁不存在翼板而没有剪力滞效应，此时，梁的动力特性等于欧拉梁的动力特性.因此，式(13)中等式右边第二项为剪力滞效应所引起梁动力特性变化.当I_s>0时，α>1，由此使${{\bar \lambda }_j} < {\lambda _j}$，这表明剪力滞使梁的频率降低.对于等式右边第二项分母中的$\frac{{14\alpha \rho {{\bar \lambda }_j}}}{{9E}}$项是剪力滞所引起的惯性效应.若忽略剪力滞的惯性效应，即为吴有俊等^[7]所得箱形简支梁动力特性的计算公式.

由式(13)可得${{\bar \lambda }_j}$的两个根，其数值小的根为

$ {{\bar \lambda }_j} = 0.5\xi-\sqrt {0.25{\xi ^2}-\eta }, $

(14)

式中$\xi = \frac{{9E}}{{14\rho }}\left[{{{\left( {\frac{{j{\rm{\pi }}}}{l}} \right)}^2} + \frac{{{\beta ^2}}}{\alpha }} \right] + {\lambda _j}, \eta = \frac{{9E{\lambda _j}}}{{14\rho \alpha }}\left[{{{\left( {\frac{{j{\rm{\pi }}}}{l}} \right)}^2} + {\beta ^2}} \right]$.

将式(14)代入式(12b)可得q_jj的解为

$ {q_{jj}} = \frac{{4I}}{{3{I_s}}}\frac{{{{\bar \lambda }_j}-{\lambda _j}}}{{{\lambda _j}}}. $

(15)

由式(14)可知，采用直接模态摄动法计算等截面箱梁振动特性式的计算结果与式(8)中n的选取无关.当n趋近于无穷时，式(8)中代表了全模态的展开，是精确的坐标变换.所以式(14)和式(15)所得的特征值和特征向量是简支箱梁动力特性的精确解.

3 模态的剪力滞系数

由式(8)和(15)可知，箱形简支梁的模态为

$ {{\bar W}_j} = {W_j}, {{\bar U}_j} = {q_{jj}}{W_j}', $

(16)

则第j阶模态翼板中的正应力为

$ {\sigma _{ij}}\left( x \right) = E{h_i}\left[{\bar W'{'_j} + \left( {1-\frac{{{y^3}}}{{{b^3}}}} \right)\bar U{'_j}} \right], $

(17)

式(16)模态位移下，箱梁任意截面的弯矩为

$ {M_j}\left( x \right) =- \left[{EI\bar W'{'_j}\left( x \right) + \frac{3}{4}E{I_s}\bar U{'_j}\left( x \right)} \right]. $

(18)

在式(18)弯矩作用下，按欧拉梁理论所得翼板的正应力为

$ {{\bar \sigma }_{ij}}\left( x \right) =- \frac{{M\left( x \right){h_i}}}{I} = E{h_i}\left[{\bar W'{'_j}\left( x \right) + \frac{3}{4}\frac{{{I_s}}}{I}\bar U{'_j}\left( x \right)} \right]. $

(19)

根据剪力滞系数的定义：${\gamma _j} = \frac{{{\sigma _j}}}{{\bar \sigma }}$，由式(15)~(19)，可得到模态的剪力滞系数为

$ {\gamma _j}\left( x \right) = \frac{{{\lambda _j} + \left( {1-\frac{{{y^3}}}{{{b^3}}}} \right)\frac{{4I}}{{3{I_s}}}\left( {{{\bar \lambda }_j}-{\lambda _j}} \right)}}{{{{\bar \lambda }_j}}}. $

(20)

由式(20)可知，对简支梁而言，剪力滞系数和x无关，即纵向各截面的剪力系数相同.箱梁翼板和腹板交角处的剪力滞系数为

$ \gamma _j^e = \frac{{{\lambda _j}}}{{{{\bar \lambda }_j}}}. $

(21)

箱梁翼板和腹板交角处的剪力滞系数为最大剪力滞系数.最大剪力滞系数用于度量箱梁剪力滞影响的大小.式(21)表明单箱单室截面翼板和腹板交角处的剪力滞系数等于欧拉梁特征值和剪力滞影响下箱梁特征值的比.这表明对于箱梁自由振动而言，最大剪力滞系数既反映了箱梁剪力滞的大小，又反映了由于剪力滞引起梁自振频率的变化.最大剪力滞系数越大，梁的自振频率降低越多.

4 计算精度分析

式(14)所得的特征值和特征向量是在忽略剪切变形和转动惯量影响后所得的精确解.为验证本文方法计算结果的精度，下面对箱型截面梁分别进行欧拉梁理论，本文方法和有限元方法的自由振动计算.其中有限元方法是采用壳单元的分析结果.箱梁截面的材料参数为：E=35 GPa，G=15 GPa.箱梁截面形式和坐标如图(1)所示，各部位的尺寸为：t=0.25 m，b=3.55 m，h=2 m.当I_s/I=0.88时，腹板宽度为t_w=0.4 m；当I_s/I=0.94时，腹板宽度为t_w=0.2 m，同时选取梁的跨度为40、30和20 m.不同方法所得箱梁的前三阶的自振频率见表 1.

由计算结果可知本文的方法由于忽略了转动惯量、剪切变形以及翼缘板振动的影响而与有限元结果有一定差别，但是误差并不大，可满足工程需要.同时，本文方法所得各阶模态的频率都小于欧拉梁的频率，且模态阶数越高，两者的频率相差越大，这表明剪力滞对梁的自由振动有显著影响.

5 参数分析

以跨宽比(l/2b)和翼板相对转动惯量(I_s/I)为参数讨论箱梁自由振动时，模态剪力滞系数和梁自振频率的变化规律.参数分析时箱梁截面的跨度和腹板厚度为变量，其余参数同前.

5.1 模态剪力滞系数的横向分布

箱梁顶板的剪力滞分布沿箱梁中轴线左右对称，因此只画出一半的剪力滞系数.取梁的跨度为40 m，前4阶模态箱梁上翼板的剪力滞系数见图 2.

由计算结果可知：1)腹板附近剪力滞系数最大，且γ_j^e≥1，这表明简支梁自由振动时，梁的剪力滞效应为正剪力滞效应.事实上，由式(13)可知，剪力滞使梁的特征值降低，因此，γ_j^e≥1是剪力滞使梁振动频率降低的必然结果；2)前四阶模态的最大剪力滞系数分别为1.045、1.176、1.389和1.706.这表明随着振动阶数的增加，剪力滞后效应越来越明显.剪力滞使箱梁前四阶模态的频率分别降低了2.2%、7.8%、15.1%和23.4%.因此，对于有高阶模态参与振动的箱梁，必须考虑剪力滞效应的影响.

5.2 不同参数对最大剪力滞系数的影响

简支梁的最大剪力滞系数随l/2b和I_s/I的变化曲线见图 3、4.

由图可知:1)l/2b越大，γ^e越小.因此相同截面情况下，梁跨度越大，剪力滞影响越小.对于第一阶模态，当l/2b >4时，γ^e < 1.05.因此，对于仅需考虑第一阶模态振动的箱梁可忽略剪力滞效应的影响；2)I_s/I反映翼缘板刚度占总刚度的百分比.I_s/I显著地影响箱梁的剪力滞后效应，由此也显著影响箱梁自振频率.I_s/I越大，γ^e越大.这说明I_s/I比值越大，频率降低越多.但是，对于第一阶模态而言，I_s/I < 0.9时，γ^e < 1.05，此时剪力滞影响很小，可以忽略不计.但是对于2阶以上模态，当I_s/I >0.4时，则γ^e>1.05，此时，剪力滞影响不可忽略.

6 结论

本文基于欧拉梁的特征值和模态，利用模态摄动法将箱梁的自由振动方程组转化为非线性代数方程组来求解，从而简化计算.应用于等截面简支梁时，可得到箱梁主频率和模态的精确解.基于箱型梁的模态，进一步推导了模态剪力滞系数.由理论分析和数值计算可得如下结论：

1) 对于简支梁而言，模态剪力滞系数沿梁轴线方向不变.且单箱单室截面翼板和腹板交角处的剪力滞系数等于欧拉梁特征值和剪力滞影响下箱梁特征值的比.

2) 腹板附近剪力滞系数最大，且γ_j^e≥1，这表明简支箱梁自由振动时，梁的剪力滞效应为正剪力滞效应.

3) 随着模态阶数的增加，剪力滞效应越来越大，由此导致箱梁高阶模态的自振频率显著降低.

4) l/2b越小，I_s/I越大，γ^e越大.当l/2b >4或I_s/I < 0.9时，γ^e < 1.05.此时，对于仅需考虑第一阶模态振动的箱梁可忽略剪力滞效应的影响.除此以外，剪力滞效应对结构动力反应的影响不可忽略，因此，一旦激振荷载能激起箱梁高阶模态的振动，剪力滞效应将显著地影响结构的动力反应.