基于修正位移模式的宽箱梁剪力滞梁段有限元法

Beam-segment finite element method for shear lag of wide box girder based on modified displacement mode

摘要

以宽箱梁剪翘变形状态为研究对象，选取附加挠度为广义位移，在箱梁截面上引入翘曲中心和应力自平衡两个整体修正系数；从宽箱梁翼板横向位移与附加挠度的几何关系入手，综合运用广义坐标法、翘曲剪切变形协调条件和平衡微分方程建立了宽箱梁的剪力滞修正位移模式，利用能量变分法提出了一种适用于宽箱梁剪力滞分析的8自由度梁段有限元法，并通过等（变）截面连续箱梁算例来考证所提方法的正确性。研究表明：修正分析方法具有一般性，当忽略翘曲中心修正和翼板横向位移时可退化为传统分析方法；修正分析方法计算结果与有限元数值解吻合良好，传统分析方法与有限元结果的差值比随宽跨比的增大而增大，且变截面箱梁影响尤为显著；三跨变截面连续箱梁宽跨比取0.42时，传统分析方法计算的中跨跨中顶板肋处的剪力滞系数与有限元结果的差值比为20.65%，对应规范方法与修正分析方法计算结果的差值比为8.17%；以上充分说明传统分析方法不适用于宽箱梁的剪力滞效应分析，规范方法修正应考虑翼板横向位移和翘曲中心差异这两个主要因素。

Abstract

The shear warping deformation state of wide box girder is taken as the research object， the additional deflection is selected as the generalized displacement， and two overall correction coefficients for the warping center and stress self-balance are introduced on the box section. Starting from the geometric relationship between the lateral displacement and additional deflection for the box flange， the modified displacement mode for shear lag was established by comprehensively using the generalized coordinate method， the warping shear deformation coordination condition， and the equilibrium differential equation. A 8-degree-of-freedom beam-segment finite element method for shear lag is proposed using the energy variational method， and the correctness of the modified analysis method is verified through the examples of continuous box girder with equal or variable cross-section. Research has shown that the modified analysis method has generality， and when ignoring the correction of the warping center and the lateral displacement of the flange， it can degenerate into the traditional analysis methods. The calculation results of the modified analysis method are in good agreement with the finite element solutions， the difference ratio between the results of the traditional analysis method and the finite element increases with the increase of width-span ratio， and the influence of variable cross-section box girders is particularly significant. For a three span variable cross-section continuous box girder， the width-span ratio is 0.42， the difference ratio between the shear lag coefficient calculated by the traditional analysis method and the finite element result is 20.65% at the intersection of top and web plate for mid-span section， corresponding to a difference ratio of 8.17% between the calculation results of the standard method and the modified analysis method. Fully demonstrate that traditional analysis methods are not applicable to the analysis of shear lag effects in wide box girders， and that the revision of standard methods should consider the two main factors for lateral displacement of flanges and differences in warping centers.

Graphical abstract

关键词

宽箱梁 / 剪力滞 / 修正位移模式 / 梁段有限元 / 剪力滞系数

Key words

wide box girder / shear lag / modified displacement mode / beam-segment finite element / shear lag coefficient

0 引 言

薄壁箱梁剪力滞引起的非均匀应力分布和挠度增大等现象在桥梁设计中必须加以考虑^［1，2］。近年来，许多学者在窄箱梁剪力滞解析理论、试验和有限元仿真等方面已开展了大量研究^［3⁃6］，但已有研究表明，既有剪力滞变分解析法用于宽箱梁时出现解析精度低和分布规律不合理等现象^{［1，7，8］}，开展宽箱梁剪力滞研究对箱梁桥合理参数化设计及现行规范翼缘有效宽度规定的修订提供科学依据。

Reissner^［9］最早用最小势能驻值原理建立了双轴对称矩形截面箱梁的剪力滞解析解，自20世纪八十年代以来，被推广至带悬臂板的箱梁剪力滞分析中^{［4⁃6，10⁃12］}。然而，悬臂板的存在使箱梁翘曲正应力不再满足自平衡条件，钱寅泉等^［10］在剪力滞位移模式中引入翘曲幅值差异修正系数和均匀轴向位移的思路，使带悬臂板的斜腹板箱梁剪力滞变分解析法更趋于严密，至今被广泛应用于混凝土箱梁^{［11，12］}、组合箱梁^{［13，14］}及钢箱梁^［8，15］的剪力滞分析；但也存在一定的缺陷，当悬臂板宽与顶板半宽相等时，悬臂板翘曲位移修正系数失效。张元海等^［16］提出了考虑上翼板边界约束差异影响的悬臂板边界约束修正系数。此外，还有学者^［4，17］从箱梁弯曲剪流分布规律和翼板翘曲位移差方面做了剪力滞位移模式的修正研究。综上所述，既有文献在窄梁剪力滞方面做了许多研究工作，但将其运用于宽翼缘箱梁时，出现解析精度低及分布规律不客观等问题。从宽箱梁挠曲变形状态来看，上、下翼板除发生纵向翘曲外，还会在面内产生横向位移，宜看作四周支承的双向板来分析。因此，宽箱梁剪力滞位移模式中应补充翼板横向位移模式。

三维有限元软件可精确分析箱梁的剪力滞效应，但设计人员更喜欢用与现行规范相适应的一维有限梁段法分析剪力滞引起的内力和变形。传统广义坐标法可准确建立剪力滞梁段单元刚度矩阵和荷载列阵^{［18，19］}，但求解过程冗杂，不便于工程应用。有学者以剪力滞齐次解为位移形函数，建立箱梁剪力滞梁段有限元法^{［20，21］}，一定程度上简化了求解过程，但剪力滞与初等梁变形耦合在一起分析，无法直观反映剪力滞对箱梁变形和内力的影响。Zhang^［22］和Arici^［23］等提出了以剪力滞变形状态为独立分析对象的梁段有限元法，可直接分析剪力滞引起的内力和变形，便于工程人员理解和使用。

本文以剪力滞引起的附加挠度为广义位移，综合运用薄壁杆件分析理论和有限元法，构建宽箱梁的翼板横向位移模式，提出一种基于修正位移模式的宽箱梁剪力滞梁段有限元法，并通过等（变）截面连续箱梁算例来考证所提方法的正确性。

1 修正剪翘位移模式

图1为任意竖向对称荷载q（z）作用下的箱梁挠曲变形分解示意，其中：u₀为初等梁纵向位移，u_ω 为剪力滞翘曲（剪翘）纵向位移，η和d分别为剪翘自平衡系数和翘曲中心修正系数；坐标系oxyz和\tilde{o}\tilde{x}\tilde{y}z分别用于分析初等梁和剪翘两种变形状态，坐标原点\tilde{o}和o位于翘曲中心和截面形心处；2b₁、2b₂、b₃为顶、底板宽和悬臂板宽；h_u 、h_b 为顶、底板中面至形心轴（即x轴）的距离，{\tilde{h}}_u和{\tilde{h}}_b为顶、底板中面至翘曲中性轴（即\tilde{x}轴）的距离，h为顶、底板中面之间的距离；t_u 、t_b 为顶、底板厚度，t_w 为腹板厚度；θ为腹板倾角。

基于弯曲剪力滞耦合变形状态，选取剪力滞附加挠度f（z）为广义位移^{［16，22］}，则箱梁截面任一点的纵向位移u（x，y，z）为：

u\left(x,y,z\right)=

-y\left(w^{\text{'}}(z)+f^{\text{'}}\left(z\right)\right)+\eta {\omega }_{\xi i}\left(x,y\right)f^{\text{'}}\left(z\right)=

-y{w}^{\text{'}}(z)-\left(y-{\eta }_i{\omega }_{\xi i}\left(x,y\right)\right)f^{\text{'}}\left(z\right)=

-y{w}^{\text{'}}(z)-{\omega }_i\left(x,y\right)f^{\text{'}}\left(z\right)（1）

式中：w'（z）为初等梁挠曲转角；f '（z）为剪力滞附加挠曲转角；ω_ξi （x，y）为剪翘基本位移函数，ω_i （x，y）为剪翘位移分布函数，i取c和t，分别代表修正和传统分析方法，下同。

考虑翘曲中心与截面形心差异影响的修正剪翘基本位移函数ω_ξc （x，y）为：

{\omega }_{\xi c}(x,y)=\left\{\begin{array}{l}-{\tilde{h}}_u\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }x}{\mathrm{2}{b}_{\mathrm{1}}}\right),\mathrm{顶板}\\ {\tilde{h}}_b{\tilde{\beta }}_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }x}{\mathrm{2}{b}_{\mathrm{2}}}\right),\mathrm{底板}\\ -{\tilde{h}}_u{\beta }_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }(b_{\mathrm{1}}+b_{\mathrm{3}}-x)}{\mathrm{2}{b}_{\mathrm{3}}}\right),\mathrm{悬臂}\mathrm{板}\end{array}\right.（2）

式中：{\tilde{\beta }}_{\mathrm{1}}和β₂分别为底板、悬臂板剪翘基本位移函数修正系数，即{\tilde{\beta }}_{\mathrm{1}}=（b₂∙{\tilde{h}}_b）/（b₁∙{\tilde{h}}_u），β₂=b₃/b₁。

若忽略翘曲中性轴与截面形心轴的差异（d=0），则退化后的剪翘基本位移函数ω_ξt （x，y）为：

{\omega }_{\xi t}(x,y)=\left\{\begin{array}{l}-h_u\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }x}{\mathrm{2}{b}_{\mathrm{1}}}\right),\mathrm{顶板}\\ h_b{\beta }_{\mathrm{1}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }x}{\mathrm{2}{b}_{\mathrm{2}}}\right),\mathrm{底板}\\ -h_u{\beta }_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }(b_{\mathrm{1}}+b_{\mathrm{3}}-x)}{\mathrm{2}{b}_{\mathrm{3}}}\right),\mathrm{悬臂}\mathrm{板}\end{array}\right.（3）

式中：β₁=（b₂∙h_b ）/（b₁∙h_u ）。

根据胡克定律，则箱梁截面上任一点的翘曲正应力σ_ωi （x，y，z）表达为：

{\sigma }_{\omega i}\left(x,y,z\right)=E\frac{\partial u_{\omega i}}{\partial z}=-E{\omega }_i\left(x,y\right)f^{″}\left(z\right)（4）

结合翘曲自平衡特性，将式（4）代入{\int }_A{\sigma }_{\omega c}\mathrm{d}A=\mathrm{0}和{\int }_A{\sigma }_{\omega c}y\mathrm{d}A=\mathrm{0}，可得修正翘曲位移模式中的两个参数η_c 和d分别为：

\left\{\begin{array}{l}{\eta }_c=\frac{\mathrm{\pi }{\tilde{I}}_x}{\mathrm{2}\left({\tilde{I}}_{\mathrm{1}}+{\tilde{\beta }}_{\mathrm{1}}{\tilde{I}}_{\mathrm{2}}+{\beta }_{\mathrm{2}}{\tilde{I}}_{\mathrm{3}}\right)}\\ d=\frac{{\tilde{\beta }}_{\mathrm{1}}{h}_b{A}_{\mathrm{2}}-h_u\left(A_{\mathrm{1}}+{\beta }_{\mathrm{2}}{A}_{\mathrm{3}}\right)}{A_{\mathrm{1}}+{\tilde{\beta }}_{\mathrm{1}}{A}_{\mathrm{2}}+{\beta }_{\mathrm{2}}{A}_{\mathrm{3}}}\end{array}\right.（5）

式中：{\tilde{I}}_i为各板绕\tilde{x}轴的惯性矩；{\tilde{I}}_x为横截面绕\tilde{x}轴的惯性矩；A₁、A₂、A₃为箱梁各板横截面积。

同理，将式（4）代入{\int }_A{\sigma }_{\omega t}y\mathrm{d}A=\mathrm{0}可得传统翘曲自平衡系数η_t ：

{\eta }_t=\frac{\mathrm{\pi }{I}_x}{\mathrm{2}\left(I_{\mathrm{1}}+{\beta }_{\mathrm{1}}{I}_{\mathrm{2}}+{\beta }_{\mathrm{2}}{I}_{\mathrm{3}}\right)}（6）

式中：I_i 为各板绕x轴的惯性矩；I_x 为截面绕x轴的惯性矩。

与初等梁弯矩M_x 对应的翘曲广义力矩M_ωi 为：

M_{\omega i}={\int }_A{\sigma }_{\omega i}{\omega }_{i}dA=-E{I}_{\omega i}{f}^{″}（7）

式中：I_ωi 为剪翘惯性矩，即I_{\omega i}={\int }_A{{\omega }_i}^{\mathrm{2}}\mathrm{d}A。

根据初等梁和剪翘纵向应力表达式，则箱梁截面任一点的剪力滞系数λ_i 为：

{\lambda }_i=\mathrm{1}+\frac{M_{\omega i}}{M_x}\cdot \frac{I_x{\omega }_i}{I_{\omega i}y}（8）

式中：M_x 为截面绕x轴的弯矩。

对于宽翼缘箱梁竖向挠曲分析时，宜当作双向受力板来处理，即考虑翼板横向变位对宽梁剪力滞效应的影响。图2为箱梁在竖向对称荷载作用下发生剪力滞附加挠曲位移f时的分析图示，利用三维空间有限元软件绘制了箱梁翼板的横向位移分布模式，根据广义坐标法，则箱梁翼板横向位移模式v_ω （x，y，z）定义为：

v_{\omega }(x,y,z)=\varphi (x,y)f(z)（9）

式中：ϕ（x，y）为翼板横向位移分布函数。

由于箱梁顶板、悬臂板和底板的横向位移受腹板侧向约束存在差异，故在横向分布模式中引入3个修正系数，结合图2翼板横向位移分布图，则箱梁翼板横向位移分布函数ϕ（x，y）为：

\varphi (x,y)=\left\{\begin{array}{l}{\alpha }_{\mathrm{1}}\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }x}{\mathrm{2}{b}_{\mathrm{1}}}\right)\right),\mathrm{顶板}\\ -{\alpha }_{\mathrm{2}}\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }x}{\mathrm{2}{b}_{\mathrm{2}}}\right)\right),\mathrm{底板}\\ {\alpha }_{\mathrm{1}}+{\alpha }_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\frac{\mathrm{\pi }(b_{\mathrm{1}}+b_{\mathrm{3}}-x)}{\mathrm{2}{b}_{\mathrm{3}}}\right),\mathrm{悬臂}\mathrm{板}\end{array}\right.（10）

式中：α₁、α₂、α₃为翼板横向位移修正系数。

根据翼板面内翘曲正应力σ_ω 与剪应力τ_ω 的平衡微分方程，可导出翘曲纵向位移分布函数ω_c 与横向位移分布函数ϕ之间的微分关系：

\frac{\mathrm{2}{\mathrm{\pi }}^{\mathrm{2}}(\mathrm{1}+\mu )}{l^{\mathrm{2}}}{\omega }_c-{\omega }_c^{″}+{\varphi }^{\text{'}}=\mathrm{0}（11）

式中：μ为材料的泊松比；l为箱梁跨度。

将修正翘曲纵向位移分布函数式（2）、横向位移分布函数式（10）代入式（11），可得横向位移修正系数α_i 为：

\left\{\begin{array}{l}{\alpha }_{\mathrm{1}}=\frac{\mathrm{4}\mathrm{\pi }(\mathrm{1}+\mu ){\tilde{h}}_u{b}_{\mathrm{1}}}{l^{\mathrm{2}}}\\ {\alpha }_{\mathrm{2}}=\frac{\mathrm{4}\mathrm{\pi }(\mathrm{1}+\mu ){\tilde{h}}_b{b}_{\mathrm{2}}}{l^{\mathrm{2}}}\\ {\alpha }_{\mathrm{3}}=\frac{\mathrm{4}\mathrm{\pi }(\mathrm{1}+\mu ){\tilde{h}}_u{b}_{\mathrm{3}}}{l^{\mathrm{2}}}\end{array}\right.（12）

由式（12）可以看出，横向位移修正系数与泊松比、翼板宽度及其至翘曲中心的距离呈正比，与跨度的平方呈反比，即跨度越大，修正系数越小，从理论上证实了窄梁剪力滞分析时可忽略横向位移。

2 控制微分方程及梁段有限元

在任意竖向对称分布荷载q（z）作用下，箱梁弯曲总势能Π为：

\Pi =\frac{\mathrm{1}}{\mathrm{2}}\iiint \left(E{\epsilon }_{\omega i}^{\mathrm{2}}+G{\gamma }_{\omega i}^{\mathrm{2}}\right)\mathrm{d}x\mathrm{d}y\mathrm{d}z-{\int }_o^{l}qf\mathrm{d}z=

\frac{E{I}_{\omega i}}{\mathrm{2}}{\int }_{\mathrm{0}}^l{f^{″}}^{\mathrm{2}}\mathrm{d}z+\frac{G{A}_{\omega i}}{\mathrm{2}}{\int }_{\mathrm{0}}^l{f^{\text{'}}}^{\mathrm{2}}\mathrm{d}z-{\int }_{\mathrm{0}}^{l}qf\mathrm{d}z（13）

式中：E为弹性模量；G为剪切模量；ε_ωi 为翘曲正应变；γ_ωi 为翘曲剪应变；修正翘曲面积A_{\omega c}={\int }_A{\left(\varphi -{\omega }_c^{\text{'}}\right)}^{\mathrm{2}}\mathrm{d}A，传统翘曲面积A_{\omega t}={\int }_A{{\omega }_t^{\text{'}}}^{\mathrm{2}}\mathrm{d}A。

利用最小势能驻值原理对总势能式（13）求一阶变分运算，并令δΠ=0，可得剪力滞附加挠度控制微分方程为：

{f^{″}}^{″}-k_i^{\mathrm{2}}{f}^{″}=\frac{q}{E{I}_{\omega i}}（14）

式中：k_i 为Reissner参数^［9］，即k_i=\sqrt[]{\frac{G{A}_{\omega i}}{E{I}_{\omega i}}}。

根据总势能一阶变分运算的边界项可以得知，剪力滞附加挠度f对应的广义剪力Q_ω 为：

Q_{\omega }=-E{I}_{\omega }(f^{‴}-k_i^{\mathrm{2}}{f}^{\text{'}})（15）

求解微分方程式（14）所需的边界条件为：

固定端：f=\mathrm{0}，f^{\text{'}}=\mathrm{0}；简支端：f=\mathrm{0}，f^{″}=\mathrm{0}；自由端：f^{″}=\mathrm{0}，f^{‴}-k_i^{\mathrm{2}}{f}^{\text{'}}=\mathrm{0}。

式（14）对应的齐次线性微分方程的通解为：

f=C_{\mathrm{1}}+C_{\mathrm{2}}z+C_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)+C_{\mathrm{4}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)（16）

将式（16）代入剪力滞广义内力表达式，并令z=0，可导出用f₀、f₀'、M_ω0、Q_ω0四个初参数表示的剪力滞初参数解为：

f=f_{\mathrm{0}}+f_{\mathrm{0}}^{\text{'}}\frac{\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)}{k_i}+

\frac{M_{\omega \mathrm{0}}}{k_i^{\mathrm{2}}E{I}_{\omega i}}\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)\right)+\frac{Q_{\omega \mathrm{0}}}{k_i^{\mathrm{3}}E{I}_{\omega i}}\left(k_{i}z-\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)\right)（17）

f^{\text{'}}=f_{\mathrm{0}}^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)-\frac{M_{\omega \mathrm{0}}}{k_{i}E{I}_{\omega i}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)+

\frac{Q_{\omega \mathrm{0}}}{k_i^{\mathrm{2}}E{I}_{\omega i}}\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)\right)（18）

M_{\omega }=-f_{\mathrm{0}}^{\text{'}}{k}_{i}E{I}_{\omega i}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)+M_{\omega \mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)+

\frac{Q_{\omega \mathrm{0}}}{k_i}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)（19）

Q_{\omega }=Q_{\omega \mathrm{0}}（20）

式（17）~（20）中的4个初参数由对应边界条件确定。

为便于分析变截面宽箱梁的剪力滞效应，提出一种基于修正位移模式的8自由度梁段单元，其节点位移和荷载如图3所示，单元位移列向量 δ^e 为：

{\mathit{\delta }}^e={\left[\begin{array}{ccc}\begin{array}{cc}{w}_m& {\theta }_m\end{array}& f_m& \begin{array}{cccc}{\phi }_m& w_n& {\theta }_n& \begin{array}{cc}{f}_n& {\phi }_n\end{array}\end{array}\end{array}\right]}^{\mathrm{T}}（21）

式中：w_m 、θ_m 和w_n 、θ_n 分别为梁单元m端和n端对应的初等梁挠度和转角；f_m 、φ_m 和f_n 、φ_n 分别为梁单元m端和n端对应的剪力滞附加挠度和转角。

单元节点荷载列向量 F^e 为：

{\mathit{F}}^e={\left[\begin{array}{ccc}\begin{array}{cc}{Q}_m& M_m\end{array}& Q_{\omega m}& \begin{array}{cccc}{M}_{\omega m}& Q_n& M_n& \begin{array}{cc}{Q}_{\omega n}& M_{\omega n}\end{array}\end{array}\end{array}\right]}^{\mathrm{T}}（22）

式中：Q_m 、M_m 和Q_n 、M_n 分别为梁单元m端和n端对应的初等梁剪力和弯矩；Q_ωm 、M_ωm 和Q_ωn 、M_ωn 分别为梁单元m端和n端对应的剪力滞广义剪力和力矩。

该梁单元对应一个8×8的刚度矩阵，其中初等梁位移对应的矩阵元素可从文献［4，22，23］直接查阅，此处不再赘述。下面仅推导剪力滞位移对应的矩阵元素，以剪力滞附加挠度初参数解为单元位移分量，则剪力滞附加挠曲转角分量φ（z）为：

\phi (z)=-f^{\text{'}}=-f_{\mathrm{0}}^{\text{'}}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)+\frac{M_{\omega \mathrm{0}}}{k_{i}E{I}_{\omega i}}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)-

\frac{Q_{\omega \mathrm{0}}}{k_i^{\mathrm{2}}E{I}_{\omega i}}\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)\right)（23）

对于两端固定的梁单元，m端发生附加挠度f_m =1，其余节点位移均为零（即φ_m =0，f_n =0，φ_n =0），可得单元刚度矩阵第3列对应的剪力滞元素（即k₃₃、k₄₃、k₇₃和k₈₃）；同理，可求得第4列（即k₃₄、k₄₄、k₇₄和k₈₄）、第7列（即k₃₇、k₄₇、k₇₇和k₈₇）和第8列（即k₃₈、k₄₈、k₇₈和k₈₈）刚度矩阵剪力滞元素；则剪力滞对应的单元刚度矩阵上三角元素为：

k_{\mathrm{33}}=k_{\mathrm{77}}=-k_{\mathrm{37}}=\frac{k_i^{\mathrm{3}}E{I}_{\omega i}}{\xi }\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)；

k_{\mathrm{34}}=k_{\mathrm{38}}=-k_{\mathrm{47}}=-k_{\mathrm{78}}=-\frac{k_i^{\mathrm{2}}E{I}_{\omega i}}{\xi }\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)-\mathrm{1}\right)；k_{\mathrm{44}}=k_{\mathrm{88}}=\frac{k_{i}E{I}_{\omega i}}{\xi }\left(k_{i}l\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)-\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)\right)；

k_{\mathrm{48}}=\frac{k_{i}E{I}_{\omega i}}{\xi }\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)-k_{i}l\right)。

式中：\xi =k_{i}l\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)-\mathrm{2}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)+\mathrm{2}。

将梁单元节点位移代入式（18）和式（23），整理可得单元剪力滞位移的矩阵{\mathit{N}}_f表达式为：

f(z)={\mathit{N}}_f\delta（24）

{\mathit{N}}_f=\left[\begin{array}{ccc}\begin{array}{cc}\mathrm{0}& \mathrm{0}\end{array}& N_{\mathrm{1}}(z)& \begin{array}{cccc}{N}_{\mathrm{2}}(z)& \mathrm{0}& \mathrm{0}& \begin{array}{cc}{N}_{\mathrm{3}}(z)& N_{\mathrm{4}}(z)\end{array}\end{array}\end{array}\right]

其中：

N_{\mathrm{1}}(z)=\frac{\mathrm{1}}{\xi }\left[\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)+k_{i}l-k_{i}z\right)+\right.

\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)\right)\left.\left(\mathrm{1}+\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)\right)\right]

\begin{array}{l}{N}_{\mathrm{2}}(z)=\frac{\mathrm{1}}{k_i\xi }\left[\left(k_{i}l\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)-\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)\right)\right.\cdot \\ \left(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)-\mathrm{1}\right)+\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)-\mathrm{1}\right)k_{i}z+\\ \left.\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)-\mathrm{1}-k_{i}l\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)\right)\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)\right]\end{array}

\begin{array}{c}{N}_{\mathrm{3}}(z)=\frac{\mathrm{1}}{\xi }\left[\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)\left(k_{i}z-\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)\right)+\right.\\ \left.\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)\right)\left(\mathrm{1}-\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)\right)\right]\end{array}

\begin{array}{c}{N}_{\mathrm{4}}(z)=\frac{\mathrm{1}}{k_i\xi }\left[\left(\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)-k_{i}l\right)\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}z)-\mathrm{1}\right)+\right.\\ \left.\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)-\mathrm{1}\right)\left(k_{i}z-\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}z)\right)\right]\end{array}

当单元上作用竖向均布荷载q时，运用虚功原理可求得对应于剪力滞位移的等效节点荷载向量{\mathit{P}}_q^f为：

{\mathit{P}}_q^f=q{\int }_l{\mathit{N}}_f^{\mathrm{T}}\mathrm{d}z=\frac{q}{\mathrm{2}}{\left[\begin{array}{cccccccc}\mathrm{0}& \mathrm{0}& l& -\zeta & \mathrm{0}& \mathrm{0}& l& \zeta \end{array}\right]}^{\mathrm{T}}（25）

式中：\zeta =\frac{\mathrm{1}}{k_i^{\mathrm{2}}\xi }\left[{\left(k_{i}l\right)}^{\mathrm{2}}\left(\mathrm{1}+\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)\right)+\right.\mathrm{4}\left(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(k_{i}l)-\right.

\left.\left.k_{i}l\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(k_{i}l)-\mathrm{1}\right)\right]。

3 算例分析

3.1　两跨等截面连续箱梁算例

以文献［20］中的两跨等截面有机玻璃连续箱梁为例，荷载作用位置及截面尺寸如图4（a）所示，单跨跨径l=800 mm，满跨作用均布荷载q=0.2 N/mm，材料弹性模量E=2.8 GPa，泊松比μ=0.37；利用ANSYS⁃Shell63单元和ABAQUS⁃C3D8R单元分别建立箱梁有限元模型，顶、底板划分单元均按照10 mm×10 mm，腹板沿高度方向平均划分，有限元模型如图4（b）（c）所示。

为验证本文修正分析方法的正确性和适用性，基于此算例，若仅改变闭口箱室宽度，即顶板半宽b₁=80 mm、120 mm和160 mm，对应宽跨比2b₁/l=0.20、0.30和0.40；若仅改变梁高，即梁高h=64 mm和96 mm，对应高宽比h/2b₁=0.20和0.30；运用本文梁段有限元法，将连续箱梁划分为41个节点，40个单元，按照修正分析方法（CAM）和传统分析方法（TAM）计算跨中截面上翼板剪力滞系数；同时，分析并提取两种有限元计算结果，绘制上翼板剪力滞系数分布图，如图5、图6所示，关键点剪力滞系数对比如表1所示。

由图5和表1可以看出：ANSYS和ABAQUS分析结果差值比较小，考证了有限元结果的可信性；修正分析方法计算得到的等截面连续箱梁跨中截面上翼板剪力滞系数与有限元数值解吻合良好，验证了修正方法的正确性；宽跨比越大，传统分析方法计算结果与有限元数值解差值比越大，与ANSYS和ABAQUS结果相比，跨中截面顶板肋处的剪力滞系数最大差值比为8.40%和7.58%；按照《公路钢筋混凝土及预应力混凝土桥涵设计规范》JTG3362-2018^［2］中翼缘有效宽度计算方法得到的剪力滞系数与修正分析方法差值比较大，其最大差值比为8.53%；以上说明，传统分析方法和规范计算方法放大了宽箱梁的剪力滞效应，故应考虑翼板横向位移和翘曲中心差异这两个主要因素。另外，图中给出剪力滞系数横向分布图，在悬臂板上理论计算结果与有限元数值解有一定的误差，其原因在于悬臂板和顶板横向边界约束存在较大差异，构建剪力滞翘曲位移模式时未考虑翼板横向边界约束差异影响。

由图6和表1可以看出：不同高宽比下，修正分析方法计算得到的剪力滞系数与有限元数值解吻合更好，进一步说明修正分析方法具有良好的适用性；宽高比对宽箱梁的剪力滞效应影响不显著；与ANSYS和ABAQUS数值解相比，传统分析方法得到的跨中截面顶板肋处的剪力滞系数最大差值比为6.87%和11.11%；规范方法计算的剪力滞系数与修正分析方法计算结果的最大差值比为8.66%。

3.2　三跨变截面连续箱梁算例

以三跨变截面有机玻璃连续箱梁为例^［3］，荷载及截面尺寸如图7所示，各板厚度沿跨度不变，梁高按y=34+0.025x²变化，材料弹性模量E=2.6 GPa，泊松比μ=0.4；满跨作用均布荷载q=0.5 N/mm；利用ANSYS-Shell63单元和ABAQUS⁃C3D8R单元分别建立变截面箱梁有限元模型，顶、底板划分单元均按照20 mm×20 mm，腹板沿高度方向平均划分，有限元模型如图7（b）（c）所示。

基于三跨变截面连续箱梁算例，仅改变闭口箱室宽度，即顶板半宽b₁=100 mm、140 mm和180 mm时，对应中跨宽跨比2b₁/l₁=0.23、0.33和0.42；若仅改变梁高，即中跨跨中梁高h=54 mm和72 mm，梁高线形按y=h+0.025x²变化，对应高宽比h/2b₁=0.15和0.20；按照梁段有限元法，全梁共划分47个节点、46个单元，分析并计算中跨跨中截面上翼板剪力滞系数；同时分析提取两种有限元模型对应截面的剪力滞系数，绘制上翼板剪力滞系数分布图，如图8、图9所示，关键点剪力滞系数对比如表2所示。

由图8和表2可以看出：用ANSYS和ABAQUS分析得到的变截面连续箱梁数值解吻合良好，最大差值比为5.92%，考证了有限元结果的可靠性；修正分析方法计算的变截面连续箱梁中跨跨中截面剪力滞系数与有限元数值解吻合良好，进一步验证了修正分析方法的正确性；宽跨比越大，传统分析方法计算的剪力滞系数与有限元数值解差值比越大，宽跨比为0.42时，跨中截面顶板肋处的剪力滞系数与ANSYS和ABAQUS结果的差值比为20.65%和14.02%，说明传统计算方法不适用于变截面宽箱梁的剪力滞分析；按照现行《公路桥规》^［2］中翼缘有效宽度计算方法得到的剪力滞系数与修正方法计算结果的最大差值比为8.17%，显然，规范方法放大了宽箱梁的剪力滞效应。

由图9和表2可以看出：不同高宽比下，修正分析方法计算的变截面箱梁中跨跨中截面剪力滞系数与有限元数值解吻合良好；高宽比对宽箱梁剪力滞效应影响不明显；传统和修正分析方法计算的剪力滞系数与ANSYS数值解最大差值比分别为11.83%和4.14%，规范计算结果与修正分析方法的最大差值比为6.21%，说明传统分析方法和规范计算方法放大了变截面宽箱梁的剪力滞效应，规范翼缘有效宽度折减方法应充分考虑翼板横向位移和翘曲中心差异性这两个关键影响因子。

4 结 论

（1）本文建立的修正分析方法具有一般性，当忽略翘曲中心修正和翼板横向位移时，可退化为传统分析方法。

（2）参数分析表明：修正分析方法计算的截面剪力滞系数与有限元数值解吻合良好；高宽比对宽箱梁的剪力滞效应影响不明显；宽跨比越大，传统分析方法计算结果与有限元数值解差值比越大，宽跨比为0.40和0.42的等（变）截面连续箱梁跨中顶板肋处剪力滞系数差值比为8.40%和20.65%，表明传统方法不适用于宽箱梁的剪力滞效应分析。

（3）现行规范计算结果分析表明：宽跨比越大，规范方法计算的剪力滞系数与修正分析方法结果的差值比越大，宽跨比为0.40和0.42的等（变）截面连续箱梁跨中顶板肋处的剪力滞系数差值比为8.53%和8.17%；说明规范方法放大了宽箱梁的剪力滞效应，修正时应充分考虑宽箱梁翼板横向位移和翘曲中心差异这两个主要因素。

参考文献

[1] 项海帆. 高等桥梁结构理论[M]. 2版.北京: 人民交通出版社, 2018.

[2] JTG3362—2018. 公路钢筋混凝土及预应力混凝土桥涵设计规范 [S].

[3] Luo Q Z, Wu Y M, Tang J, et al. Experimental studies on shear lag of box girders[J]. Engineering Structures, 2002, 24(4): 469-477.

[4] Zhou M D, Zhang Y H, Lin P Z, et al. A new practical method for the flexural analysis of thin-walled symmetric cross-section box girders considering shear effect[J]. Thin-Walled Structures, 2022, 171: No.108710.

[5] 张玉元, 张元海, 张慧. 基于剪力滞附加挠度的箱梁弯曲自振特性分析[J]. 西南交通大学学报, 2024,59(6):1431-1439.

[6] Zhang Yu-yuan, Zhang Yuan-hai, Zhang Hui. Analysis of bending natural vibration characteristics of box girder based on shear lag additional deflection[J]. Journal of Southwest Jiaotong University, 2024,59(6):1431-1439.

[7] 李立峰, 李谋, 叶萌, 等. 变截面波形钢腹板-UHPC箱梁剪力滞效应试验与理论研究[J]. 铁道科学与工程学报, 2022, 19(6): 1675-1684.

[8] Li Li-feng, Li Mou, Ye Meng, et al. Experiment and throretical research on shear lag effect of variable section corrugated box girder with steel web and UHPC[J]. Journal of Railway Science and Engineering, 2022, 19(6): 1675-1684.

[9] Li Y F, Zhao J H, Bao L S, et al. Numerical and theoretical analysis of spatial shear lag effect in through wide box bowstring arch bridge main girder[J]. Structures, 2021, 33:141-151.

[10] Singh G. Shear flow capacity of stiffened flanges in steel box girders[J]. Thin-Walled Structures, 2023, 185: No.106957.

[11] Reissner E.Analysis of shear lag in box girders by the principle of the minimum potential energy[J]. Quarterly of Applied Mathematics, 1946, 4(3): 268-278.

[12] 钱寅泉, 倪元增. 单室箱桥的剪力滞分析[J]. 中国公路学报, 1989, 2(2): 28-38.

[13] Qian Yin-quan, Ni Yuan-zeng. Analysis of shear lag in single cell box girder bridges[J]. China Journal of Highway and Transport, 1989, 2(2): 28-38.

[14] Hu H, Zhou P M, Yang C H, et al. Effect of external prestressing on the special shear lag behavior of precast box girder bridges[J]. Structures, 2024, 60: No.105890.

[15] 张元海, 马云亮, 刘泽翔. 偏载作用下薄壁箱梁空间整体效应的综合分析[J]. 土木工程学报, 2024, 57(4): 48-58.

[16] Zhang Yuan-hai, Ma Yun-liang, Liu Ze-xiang. Comprehensive analysis on spatial integral effect of thin-walled box girders under eccentric load[J]. China Civil Engineering Journal, 2024, 57(4): 48-58.

[17] Chen Y Y, Dong J E, Xu T H, et al. The shear-lag effect of composite box girder bridges with corrugated steel webs and trusses[J]. Engineering Structures, 2019,181: 617-628.

[18] Wang C S, Zhang W T, Yang M Y. Shear lag effect of composite girders in cable-stayed bridges under dead loads[J]. Engineering Structures, 2023, 281:No. 115752.

[19] Li S Y, Huang Z W, Chen Z Q, et al. Investigation of the shear lag effect of an innovative steel box girder in cable-stayed bridges: experimental and bar simulation method analysis[J]. Journal of Bridge Engineering, 2024, 29 (3): No.04024006.

[20] 张元海, 胡玉茹, 林丽霞. 基于修正翘曲位移模式的薄壁箱梁剪力滞效应分析[J]. 土木工程学报, 2015, 48(6): 44-50.

[21] Zhang Yuan-hai, Hu Yu-ru, Lin Li-xia. Analysis on shear lag effect of thin-walled box girders based on a modified warping displacement mode[J]. China Civil Engineering Journal, 2015, 48(6): 44-50.

[22] LI X, Wan S, Mo Y, et al. An improved method for analyzing shear lag in thin-walled box-section beam with arbitrary width of cantilever flange[J]. Thin-Walled Structures, 2019, 140: 222-235.

[23] Vojnić-Purčar M, Prokić A, Bešević M. A numerical model for laminated composite thin-walled members with openings considering shear lag effect[J]. Engineering Structures, 2019, 185: 392-399.

[24] Cambronero B F, Díaz D V J, Martínez M J. Beam element for thin-walled beams with torsion, distortion, and shear lag[J]. Engineering Structures, 2017,143: 571-588.

[25] 韦成龙, 李斌, 曾庆元. 变截面连续箱梁桥剪力滞及剪切变形双重效应分析的传递矩阵法[J]. 工程力学, 2008(9): 111-117.

[26] Wei Cheng-long, Li Bin, Zeng Qing-yuan. Transfer matrix method considering both shear lag and shear deformation effects in non-uniform continuous box girder[J]. Engineering Mechanics, 2008(9): 111-117.

[27] 周世军, 周操, 李贤伟, 等. 变厚度翼缘板箱梁剪力滞分析的有限梁段法[J]. 铁道学报,2025,17(3):142-153.

[28] Zhou Shi-jun, Zhou Cao, Li Xian-wei, et al. Finite beam segment method for shear lag analysis on box girders with variable thickness flange plates[J]. Journal of the China Railway Society, 2025,17(3):142-153.

[29] Zhang Y H, Lin L X. Shear lag analysis of thin-walled box girders based on a new generalized displacement[J]. Engineering Structures, 2014, 61(3): 73-83.

[30] Arici M, Granata M, Longo G. Symplectic analysis of thin-walled curved box girders with torsion, distortion and shear lag warping effects[J]. Thin-Walled Structures, 2022, 175: No.109244.

基金资助

国家自然科学基金项目(52368020)

甘肃省高校青年博士支持项目(2024QB-047)

甘肃省博士后资助项目(386824)