潘旦光 1 ，付相球 1 ，韦杉杉 1 ，陈 钒 1,2 ，杨少平 3
(1．北京科技大学土木工程系，北京 100083；2．中电建路桥集团有限公司，北京 100048；3．中国水电建设集团十五工程局有限公司，西安 710068)
摘 要： 以等截面Euler梁的自由振动模态为Ritz基函数，提出了一种求解变高度箱梁剪力滞微分方程组的Ritz法。该方法首先进行与箱梁相同跨度相同边界条件等截面欧拉梁模态分析，然后将箱梁的竖向挠度和剪切转角用模态及其导数的线性组合来表达，从而将变分法所得箱梁剪力滞微分方程组转化为线性代数方程组进行求解。在此基础上，研究了参与计算模态阶数和截面高度变化率对计算误差的影响，算例分析结果表明：箱梁高度变化越大，Ritz法的收敛速度越慢；但随着参与计算模态阶数的增加，Ritz法将收敛到解析解。采用12阶以上模态进行计算，所得的剪力滞系数误差小于5%。
关键词： 箱梁；变高度；剪力滞；Euler梁；模态
箱型截面具有抗弯和抗扭刚度大的特点而成为梁桥中应用最广泛的一种截面。当腹板间距较大时，受剪应力的影响使翼缘板的正应力分布明显不均匀而存在剪力滞效应 [1―3] 。在剪力滞效应的各种分析方法中，利用有限元中的实体单元或壳单元可进行任何形式桥梁的力学分析 [4] ，但存在计算工作量大的特点。基于变分法的最小势能原理分析方法，既满足工程上的精度要求，又计算简便而成为最常用的分析方法 [5―6] 。对于等截面箱梁，剪力滞效应的微分方程是一个常微分方程，可以得到解析解。但实际的桥梁工程常由于力学和美学等的要求，而将桥梁设计为变高度形式。对于变高度箱梁的剪力滞效应问题是一个变系数微分方程组，难以得到解析解而常采用数值解。
张士铎和丁芸 [7] 采用中心差分法将变系数微分方程转变为代数方程组进行求解，从而得到变高度箱梁剪力滞效应的数值解。罗旗帜 [8] 、Zhang和Lin [9] 、周世军 [10] 提出了考虑竖向位移、转角和位移差自由度的梁单元，刘世忠等 [11] 还提出了可以考虑剪切变形自由度的梁单元，然后形成单元刚度矩阵和等效结点荷载的数值计算方法，并用模型实验验证了算法的有效性 [12] 。刘杰和成圣翱 [13] 采用样条配点法得到变截面箱梁微分方程的数值解。丁南宏等 [14] 提出等效刚度和等效刚度比法，将变截面箱梁近似为跨度相同的等截面箱梁，从而可利用等截面箱梁的解析解进行近似计算。
采用梁段有限单元法或有限差分法，计算的适应面广，其计算精度依赖于单元划分的数目。寻求简便且有很好精度的计算方法对箱梁的剪力滞分析有积极的意义。本文依据结构动力学中的Ritz展开定理，以相同跨度相同边界条件Euler梁的自由振动模态为Ritz函数基，建立了一种求解箱梁剪力滞效应的半解析解。然后研究箱梁截面高度变化率和参与模态阶数对计算精度的影响。
图1所示的箱梁，在竖向荷载作用下产生弯曲变形，设梁的竖向位移为w(x)，纵向位移为v(x, y)，则v(x, y)可表示为 [15] ：
式中：u(x)为剪切转角的最大差值；b为箱室净宽的一半；h i 为截面形心到上或下翼板中心距离，上翼板取为h u ，下翼板取为h b 。
在荷载作用下，体系总势能Π为：
式中：M(x)为截面的弯矩；E为弹性模量；G为剪切模量； 
为上翼板相对于截面形心的惯性矩， 
为下翼板相对于截面形心的惯性矩，t u 和t b 分别为上下翼板的厚度，ab为箱的外伸臂长度；I w 为腹板相对于截面形心的惯性矩。对于变高度箱梁而言，h、I sb 、I su 、I w 都是坐标x的函数。
图1 箱梁计算模型示意图
Fig.1 Sketch map of box girder calculation model
根据最小势能原理，在外力作用下处于平衡状态的结构，体系总势能的一阶变分为零，即：
由此可得箱梁的微分方程和边界条件 [13] ：
对于任意变高度箱梁，求解式(4)和式(5)耦合偏微分方程组的解析解通常是很困难的。而等截面Euler梁自由振动特征方程的解容易获得。本文利用等截面Euler梁自由振动的模态作为Ritz函数，形成一个求解式(4)和式(5)较简单的半解析解方法。
设有一根与变高度箱梁具有相同跨度和边界条件的等截面Euler悬臂梁，Euler梁单位长度的质量和截面抗弯刚度分别为ρA 0 和EI 0 。则 Euler梁第j阶自由振动的特征方程满足：
式中： j ω为第j阶自振频率； j φ为第j阶模态。 j φ可表达为：
式中α j 为超越方程1+cosα j lc oshα j l = 0 的解。由于ρA 0 和EI 0 的取值不影响振型函数，因此，计算中这两个参数可取任意大于零的实常数。显然式(8)的φ j 满足边界条件：
Ritz法进行微分方程求解的关键是寻找一组满足边界条件且线性无关的函数集。由于Euler悬臂梁在x=0处为固定边界， 
′(0)= 0。对于x=l处的自由边界， 
′(l )= 0 。同时， 
(x)和 
′(x )为线性无关函数，因此，以 
( x)为梁的位移 w(x)的 Ritz函数基，以 
′(x )为剪切转角的最大差值u(x)的Ritz函数基，将w(x)和u(x)分别用Ritz函数基的线性组合表示为：
因此，式(11)和式(12)所得的w(x)和u(x)必然满足式(6)的边界条件。从理论上讲，Euler梁有无穷多个主模态，即式(11)和式(12)中的 n应为∞，但在实际计算时，可取有限个低阶模态进行近似计算。将式(11)和式(12)代入式(4)和式(5)，并对式(4)方程求两阶导数，对式(5)方程求一阶导数，整理后可得：
式中 f ( x) = - d M( x)/dx。在方程两边同时乘以φ( j = 1 ,2,⋅⋅⋅,n)，然后沿全长积分，可得：
式中：
将 j 依次取值1,2,⋅⋅⋅,n ，重复利用式(15)和式(16)可得2n个代数方程，写成矩阵形式如下：
从式(17)推导过程来看，利用自由振动模态为Ritz函数基将复杂的微分方程组求解转化为一组线性代数方程的求解。
在求得未知向量 q 和 p 以后，由式(11)和式(12)不难得到变高度箱梁的位移w(x)和剪切转角的最大差值u(x)。由于箱梁翼板的应力可表示为：
将式(12)代入式(18)后可得：
因此，箱梁的剪力滞系数为：
从理论上讲，本文方法适用于任意变高度箱型梁剪力滞效应的计算。在这种方法中，式(11)和式(12)中的模态截断将产生误差。为分析模态截断的影响，下面首先对比本文方法和解析解所得等截面悬臂梁计算结果的差别。对于图1所示的等截面悬臂箱梁，在均布荷载下剪力滞系数的解析解为：
式中： 
为瑞斯纳常数。令翼板中点(y=0)的剪力滞系数为λ c ，翼板和腹板交角处(y=b)的剪力滞系数为λ e 。在均布荷载作用下本文方法和解析解所得的λ c 和λ e 沿x的变化规律及x=0和x=l/2处的剪力滞分布如图2所示。其中图2(c)和图2(d)只显示一半的剪力滞系数。
算例分析中截面的几何参数取为：h 1 =h 2 =3 m，b=2 m，t u =t b =0.4 m，t w =0.6 m，a=0.5，l=40 m。材料的弹性模量为E=50 GPa，泊松比为0.2。从理论上看，式(11)和式(12)可理解为函数的逼近问题，因此，随着n值的增加，本文方法所得λ c 和λ e 逐步逼近解析解，并收敛到解析解。
图2 等截面悬臂梁的剪力滞系数
Fig.2 Shear lag coefficient of a cantilever beam with uniform section
将解析解任一反应量记为 r * ，本文方法所得近似解记为r。两者之间的相对误差为：
则固定端(x=0)和跨中(x=l/2)λ c 和λ e 的相对误差随n的变化规律如图3所示。图3表明当n>6后，本文方法的计算误差即小于 5%。显示了良好的计算精度。
由于 
(x)和 
′(x )满足边界条件且线性无关，因此，当式(11)和式(12)中n→∞时，计算结果收敛到精确解。对于变截面箱梁的剪力滞微分方程难以得到解析解的情况，下面以 n=60的结果作为精确解，来分析本文计算方法的收敛特性。
图3 Ritz法的相对误差
Fig.3 The relative errors of Ritz method
对于截面高度线性衰减箱梁，梁高函数可表示为：
式中：β=（h 1 -h 2 )/h 1 为截面高度变化率；h 1 和 h 2 分别为固定端和自由度端截面的高度。
算例中箱梁截面的几何参数为：h 1 =80 mm，b=72 mm，t u =t b =4 mm，t w =6 mm，a=1，l=400 mm。材料的弹性模量为E=3 GPa，泊松比为0.385。荷载f(x)=0.3 N/m。不同截面高度变化率下固定端的 λ c 和λ e 的计算误差如图4所示。
λ c 和λ e 的计算误差表明：1)在相同n情况下，β越大，Ritz法的计算误差越大，这表明截面高度越均匀的箱梁，Ritz法的收敛速度越快；2)随着n的增加，Ritz法的计算误差逐步减少，当n≥12时，β≤ 0 .75的线性变截面梁的计算误差都小于5%；3)如果将静力问题看作是激振频率为 0的动力学问题。结合动力学中Ritz法的计算经验 [16] ，为得到第j阶模态的近似结果，通常采用j+8阶Ritz函数进行近似计算。在低频荷载作用下，悬臂结构主要是由前几阶模态的反应所控制，采用本文方法进行近似计算时，建议采用 12阶及以上的模态进行近似计算。
图4 线性衰减箱梁固定端剪力滞的相对误差
Fig.4The relative errors of shear lag coefficients at fixed end of tapered box girder
表1对比分析了均布荷载下线性衰减箱型梁等效刚度法 [14] 、本文方法和有限元方法的计算结果。其中等效刚度法以加权平均所得 I s /I按式(21)的计算结果，本文方法列出了n=12和n=60的计算结果。有限元分析采用ANSYS软件，β=0.5时的计算模型如图5所示。模型采用Solid185单元，模型共包括28400个单元和39996个结点。
表1 线性衰减箱梁固定端剪力滞系数
Table1 Shear lag coefficients at fixed end of tapered box girder
β =0.5时有限元解和本文方法所得固定端的剪力滞分布如图6所示。由图6可以看出：1)n=12和 n=60的剪力滞计算结果基本重合，这表明采用12阶模态进行变截面箱梁剪力滞计算是可行的；2)由于最小势能原理是在一定假定条件下的简化分析结果，主要表现在由剪力滞效应导致翼板的位移沿横向按3次函数变化和腹板采用平截面假定与实际情况有一定差别，因此，即使 n=60时，Ritz法的计算结果和有限元的计算结果也不完全相同，但不同β下，固定端的λ c 和λ e 的计算结果和有限元结果都比较接近，且变化规律相同。
图5 箱梁的有限元模型
Fig.5 Finite element model of box girder
图6 固定端的剪力滞分布(β=0.5)
Fig.6 Shear lag distribution at fixed end (β=0.5)
若梁高度按抛物线变化，即：
式中β=（h 1 -h 2 )/h 1 ，h 1 和h 2 分别为端部截面的高度。算例中箱梁截面的几何参数为：h 1 =80 mm，b=72 mm，t u =t b =4 mm，t w =6 mm，a=1，l=400 mm。材料的弹性模量为 E=3 GPa，泊松比为 0.385。以n=60的结果作为精确解，不同β下固定端的λ c 和λ e 的计算误差如图7所示。表2对比分析了均布荷载下截面高度抛物线变化箱型梁等效刚度法 [14] 、本文方法和有限元方法的计算结果。计算结果与线性衰减梁的变化规律基本相同，对于β≤0.75的抛物线梁，当n≥12时，Ritz法计算误差小于5%。
图7 抛物线箱梁固定端剪力滞的相对误差
Fig.7 The relative errors of shear lag coefficients at fixed end of parabolic box girder
表2 抛物线箱梁固定端剪力滞系数
Table2 Shear lag coefficients at fixed end of parabolic box girder
本文基于等截面Euler梁自由振动的模态，将变截面箱梁剪力滞微分方程组转化为线性代数方程组来求解。从方法介绍中可以看出，本文的计算工作主要是式(15)和式(16)中系数的积分计算。对于任意的函数，都可以采用数值积分的方法进行计算，从而使一些复杂箱梁的剪力滞问题的求解较为容易实现。这个方法的本质属于Ritz法，由于Euler梁自由振动的模态满足箱梁的位移边界条件和力的边界条件，因而计算结果具有较高的精度。由计算结果可以看出，随着参与计算模态数目的增加，本文方法的计算结果逐步收敛到解析解。对于截面衰减系数小于0.75的箱型梁，采用12阶以上模态进行剪力滞系数的计算，计算误差小于5%。
参考文献：
[1]Cambronero-Barrientos F, Diaz-Del-Valle J, Martinez-Martinez J A.Beam element for thin-walled beams with torsion, distortion and shear lag [J].Engineering Structures, 2017, 143(15): 571―588.
[2]Zhu L, Nie J G, Li F X, et al.Simplified analysis method accounting for shear-lag effect of steel-concrete composite decks [J].Journal of Constructional Steel Research, 2015, 115(7): 62―80.
[3]姚浩, 程进.基于变分原理的波形钢腹板箱梁挠度计算解析法[J].工程力学, 2016, 33(8): 177―184.>Yao Hao, Cheng Jin.Analytical method for calculatingdefection of corrugated box girders based on variational principle [J].Engineering Mechanics, 2016, 33(8): 177―184.(in Chinese)
[4]吴文清, 万水, 叶见曙, 等.波纹钢腹板组合箱梁剪力滞效应的空间有限元分析[J].土木工程学报, 2004,37(9): 31―36.Wu Wenqing, Wan Shui, Ye Jianshu, et al.3-D finite element analysis on shear lag effect in composite box girder with corrugated steel web [J].China Civil Engineering Journal, 2004, 37(9): 31―36.(in Chinese)
[5]张元海, 林丽霞.薄壁箱梁剪力滞效应分析的初参数法[J].工程力学, 2013, 30(8): 205―211.Zhang Yuanhai, Lin Lixia.Initial parameter method for analysis shear lag effect of thin-walled box girders [J].Engineering Mechanics, 2013, 30(8): 205―211.(in Chinese)
[6]甘亚南, 石飞停.梯形箱梁剪力滞后效应的精细化分析[J].计算力学学报, 2014, 31(3): 351―356.Gan Yanan, Shi Feiting.The delicate analysis of shear lag effect on trapezoidal box girders [J].Chinese Journal of Computational Mechanics, 2014, 31(3): 351―356.(in Chinese)
[7]张士铎, 丁芸.变截面悬臂箱梁负剪力滞差分解[J].重庆交通学院学报, 1984(4): 34―47.Zhang Shiduo, Ding Yun.Finite difference solution in shear lag effect on cantilever box gerder with linear varying depths [J].Journal of Chongqing Jiaotong University, 1984(4): 34―47.(in Chinese)
[8]罗旗帜.变截面多跨箱梁桥剪滞效应分析[J].中国公路学报, 1998, 11(1): 63―70.Luo Qizhi.Analysis of the shear lag effect on continuous box girder bridges with variable depth [J].China Journal of Highway and Transport, 1998, 11(1): 63―70.(in Chinese)
[9]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(1): 73―83.
[10]周世军.箱梁的剪力滞效应分析[J].工程力学, 2008,25(2): 204―208.Zhou Shijun.Shear lag analysis of box girders [J].Engineering Mechanics, 2008, 25(2): 204―208.(in Chinese)
[11]刘世忠, 欧阳永金, 吴亚平, 等.变截面薄壁箱梁剪力滞剪切变形效应分析[J].中国公路学报, 2002, 15(3):61―63, 67.Liu Shizhong, Ouyang Yongjin, Wu Yaping, et al.Non-uniform thin wall box analysis of considering both shear lag and shear deformation [J].China Journal of Highway and Transport, 2002, 15(3): 61―63, 67.(in Chinese)
[12]吴幼明, 罗旗帜, 岳珠峰.变高度连续箱梁剪力滞效应试验研究[J].实验力学, 2004, 19(1): 85―90.Wu Youming, Luo Qizhi, Yue Zhufeng.An experimental study on the shear lag effect of continuous box girder with varying depth [J].Journal of Experimental Mechanics, 2004, 19(1): 85―90.(in Chinese)
[13]刘杰, 成圣翱.变截面矩形箱梁剪力滞效应的样条配点法[J].中南公路工程, 1995(1): 35―38.Liu Jie, Cheng Shengao.Spline collocation method for the shear lag effect of non-uniform rectangle box girder[J].Central South Highway Engineering, 1995(1): 35―38.(in Chinese)
[14]丁南宏, 林丽霞, 钱永久.变截面箱梁剪力滞及剪切变形效应近似计算方法[J].铁道科学与工程学报,2011, 8(1): 14―18.Ding Nanhong, Lin Lixia, Qian Yongjiu.An approximate method to analyze the effect of shear lag and shear deformation of box beam with varying depth [J].Journal of Railway Science and Engineering, 2011, 8(1): 14―18.(in Chinese)
[15]项海帆.高等桥梁结构理论[M].北京: 人民交通出版社, 2001: 71―72.Xiang Haifan.Advanced theory of bridge structures [M].Beijing: China Communications Press, 2001: 71―72.(in Chinese)
[16]楼梦麟, 吴京宁.复杂梁动力问题的近似分析方法[J].上海力学, 1997, 18(3): 234―240.Lou Menglin, Wu Jingning.An approach to solve dynamic problems of complicated beams [J].Shanghai Journal of Mechanics, 1997, 18(3): 234―240.(in Chinese)
SEMI-ANALYTIC SOLUTION FOR SHEAR LAG EFFECT OF CANTILEVER BOX GIRDERS WITH VARYING DEPTH
PAN Dan-guang 1 , FU Xiang-qiu 1 , WEI Shan-shan 1 , CHEN Fan 1,2 , YANG Shao-ping 3
(1.Department of Civil Engineering, University of Science and Technology Beijing, Beijing 100083, China;
2.Power China Road Bridge Group Co., Beijing 100048, China; 3.Sinohydro Engineering Bureau 15 Co., Ltd, Xi’an 710068, China)
Abstract: By using the free vibration mode of an Euler beam as Ritz base function, a new Ritz method is proposed to solve the set of shear lag differential equations of a varying depth box girder.Firstly, the modal analysis of a uniform cross-sectional Euler beam which is the same length and boundary condition with the box girder is carried out.The vertical deflection and the shear rotation of the box girder are expressed by the linear combination of the model and its derivative.And then the set of shear-lag differential equations of the box girder obtained by the calculus of variations are transformed into a set of linear algebraic equations.Then, the influences of the number of modes and variation ratio of section height on the errors are investigated.The numerical examples show that: the more significant the height of the box girder varies, the slower the Ritz method converges; but the results by Ritz method would converge to the analytic solution with the increasing of the number of modes.The errors of shear lag coefficients are less than 5%, when more than 12 of modes are included.
Key words: box girder; varying depth; shear lag; Euler beam; modes of shape
杨少平(1990―)，男，陕西人，助工，学士，主要从事道路桥梁施工研究(E-mail: 505703213 @qq.com).
韦杉杉(1993―)，女，河北人，硕士生，主要从事工程结构与工程系统抗震方向的研究(E-mail: 814440337@qq.com);
陈 钒(1972―)，男，江西人，教授级高工，硕士，主要从事道路桥梁设计、施工研究(E-mail: cfsc_1037@163.com);
作者简介：
付相球(1994―)，男，江西人，硕士生，主要从事防灾减灾及防护工程研究(E-mail: 1376369310@qq.com);
基金项目： 北京市自然科学基金项目(8143037)；北京科技大学教改项目(JG2013M11).
通讯作者： 潘旦光(1974―)，男，浙江人，研究员，博士，博导，主要从事防灾减灾及防护工程研究(E-mail: pdg@ustb.edu.cn).
文章编号： 1000-4750(2018)09-0207-07
收稿日期： 2017-08-23；修改日期：2018-01-16
文献标志码： A
doi: 10.6052/j.issn.1000-4750.2017.08.0643
中图分类号： U441+5