一种多层组合梁自由振动分析的混合有限元法

蒋佳卿 王云 陈伟球 徐荣桥

摘要： 基于截面等效，提出了多层组合梁的二维分析模型。使用混合能变分原理，以包含频率函数的节点位移及其能量对偶的应力分量为单元节点未知量，引入混合元对梁长方向进行离散，建立动力学状态空间控制方程；接着引入微分求积法（Differential Quadrature Method， DQM）对控制方程进行梁高度方向的离散，求得组合梁在不同轴力与边界条件下的动力学方程。以混凝土?木材组合梁、波形钢腹板梁和钢?混组合梁为例进行验证。该方法基于二维理论，可为梁理论提供假设依据和误差分析的基准。

关键词： 自由振动； 组合梁； 混合有限元； DQM； 状态空间法

中图分类号： O327    文献标志码： A    文章编号： 1004-4523（2024）05-0856-08

DOI：10.16385/j.cnki.issn.1004-4523.2024.05.014

引  言

组合梁的动力特性是组合梁计算分析中的重要内容，利用自振频率可以识别结构可能出现的损伤。基于Timoshenko梁理论，Xu等 ［1］推导了考虑层间滑移的组合梁解析表达式，分析了其静力、动力与屈曲特性。考虑高阶剪切效应的影响，Vo等［2］提出了使用高阶多项式和非多项式函数的高阶剪切变形理论。Carrera等［3］比较多种不同高阶函数，包括多项式、三角级数、指数函数和Zig?zag函数，分析多层组合梁自振特性的异同。基于Zig?zag理论，胡霖远等［4］通过引入分段的位移线性分布假设与横向剪应力抛物线分布假设，建立了适用于分析波形钢腹板梁的自由振动分析理论。Shen等［5?6］基于状态空间架构，给出了考虑部分组合作用的组合梁自由振动控制方程，并分析了轴力的影响。通过引入“动力折减系数”与“频率折减系数”，侯忠明等［7］提出了钢?混组合梁动力问题的简便计算方法。考虑到钢?混组合梁在钢梁翼缘存在剪力滞效应，陈玉骥等［8］通过假设结构的纵向位移函数，提出了考虑剪力滞与层间滑移效应的组合梁一阶自振频率计算方法。

除了解析方法，有限元方法在组合梁动力分析中得到了较多应用。基于经典梁理论，Vo等［9］使用一维有限元单元分析了薄壁组合梁的动力特性。Chakrabarti等［10］提出了考虑层间滑移的高阶梁单元。Chalak等［11］提出了基于Zig?zag理论的有限元单元，并用其分析了包含柔软夹层的三明治结构组合梁的动力特性。

以上研究都基于一维梁理论，若将组合梁近似为二维结构，则无须假定截面上的变形或应力的分布，适用于截面材料特性变化较大、层数较多，或是边界条件沿梁高分布较为复杂的情况。Chen等［12］基于二维弹性理论，提出了使用状态空间法与微分求积法的组合梁自由振动分析方法。基于状态空间法，Xu等［13］给出了简支条件下考虑层间滑移的组合梁自由振动与屈曲荷载解析解，并使用微分求积法给出了在其他边界条件下的半解析表达式。因为将位移与其能量对偶的应力分量作为基本未知量，因此状态空间法在分析多层组合结构时具有独特优势。

本文在组合梁统一的二维模型基础上，推导了哈密顿体系下考虑固有频率与轴力作用的混合能量泛函。然后通过梁长方向的有限元离散，求得沿梁高度方向的状态空间方程。接着利用微分求积法（Differential Quadrature Method， DQM）对状态方程在梁高方向上进行分解。最后使用二分法解得组合梁的各阶自振频率。通过对混凝土?木材组合梁、波形钢腹板梁和钢?混组合梁的分析，验证了该方法的正确性。

1 组合梁的二维模型

不同形式的组合梁，包括钢?混组合梁与波形钢腹板梁，均可通过截面等效转化为宽度为的矩形截面组合梁。参考胡霖远等［4］、Johnson等 ［14］ 以及Xu等［15］的工作，通过把波形钢腹板等效为正交各向异性材料，可模拟波形钢腹板横向和纵向刚度的巨大差异。同时根据截面抗弯刚度与轴向刚度一致原则，对不同宽度的混凝土截面或钢梁截面进行材料参数转换，使其统一等效为矩形截面。

为保持模型的一致性，组合梁的层间部分组合作用通过引入厚度极小的虚拟粘接层进行模拟。如图1中高度为（其值远小于其他层的高度）的薄层，其剪切模量为：

（1）

式中为剪力连接件的抗剪刚度，该式可方便地根据能量相等原则推导得到。如果需要考虑实际情况下剪力连接件在梁长方向的离散分布特性，那么虚拟粘接层的剪切模量也可在梁长方向作相应的设置。通过一个等效的薄层材料模拟由于剪力连接件刚度有限产生的层间相对滑移，使得所有的等效截面成为统一的层间完美连接的二维模型。

2 混合有限元与DQM结合解法

2.1 混合能泛函

考虑组合梁在轴力作用下的自由振动问题，假设梁的振动频率为，则二维系统的拉格朗日函数可表示为：

式中为位移振幅列阵，u和v分别对应x方向和y方向的位移振幅；若变量下标中包含逗号，则表示该变量对逗号后面坐标的偏导数；表示材料的密度；矩阵的取值可参考胡霖远等 ［4］的工作。

选取方向作为传递方向，通过勒让德变换［16］可得到如下自由振动问题的哈密顿函数：

式中  为法线方向与x轴平行的平面上的应力振幅矢量，以及：

求得哈密顿函数表达式后，混合能泛函可表示为：

式中表示平面区域。该混合能泛函以位移振幅和应力振幅矢量为自变函数，使得该泛函取极值的可能位移振幅和应力振幅是精确解。出于简便性，后面的文字描述中统一省略“振幅”。

2.2 混合有限元离散

对结构沿x方向进行有限元离散，y方向保持不变，并对自由振动下的单元内位移与应力做如下假设：

（6a）

（6b）

式中表示位移与应力随时间的变化函数，i和t分别表示虚数与时间变量；和分别表示单元节点上的位移和应力分量；和分别表示单元内部位移与应力的分布形函数。

有限元离散后，系统的混合能泛函也可以写为：

式中  H为整个梁高；m为x方向单元数量；和为第e个单元所占的区域。

将式（3）代入式（7），并结合式（6），得：

式中

求得单元内部的系数矩阵后，按节点号组装全部矩阵，可得系统的混合能泛函为：

式中和表示系统全体位移与应力分量；，，，，表示根据有限元节点号组装后的各系数矩阵。

对式（10）取变分为零，可得如下控制方程：

式中

考虑到位移与应力具有不同的量纲，为提高数值精度，可按下面方式进行无量纲化。首先引入：

式中为特征长度，可取为；E为具有弹性模量量纲的量。

按式（13）无量纲化后，式（11）可改写为：

式中

2.3 三节点单元

有限元离散采用三节点有限元单元，如图2所示，并取中间节点k为首节点i与末节点j之间的中点。形函数表达式如下：

（16）

式中

表示单元内的局部坐标，与全局坐标x的关系为：

由于形函数由参考坐标给出，因此根据复合函数求导规则：

（19）

将形函数表达式代入式（9），可得各单元系数矩阵的表达式如下：

2.4 状态方程的DQM求解

状态方程（14）为微分方程，虽可用传递矩阵直接求解，但是由于引入有限单元后维数较大，容易导致数值稳定问题。为此，可采用微分求积法［17?18］求解状态方程（14）。首先， 在组合梁第i层上布置k个离散点，其局部坐标为（）。根据DQM法，在离散点c上，状态向量的导数满足：

式中和表示第i层内和点处的状态量；为DQM权重系数，一种取法为［18］：

对于第i层上任一离散点上的状态量，式（14）可以改写为：

式中表示式（14）中矩阵A在第i层的取值。将式（21）代入式（23），可得在任一离散点c上满足：

组合第i层上所有离散点之间的关系，可得：

（25）

式中

式中  I表示与维数相同的单位矩阵。

考虑到原组合结构具有n层，将各层系数矩阵依次代入式（25），并组合可得：

（27）

式中

式中  系数矩阵表示组合梁的自身动力特性，还需针对不同边界条件进行相应修改。

2.5 边界条件引入

边界条件通过有限元分析方法中常用的“乘大数法”进行引入，并记大数为。设边界条件共约束f个状态量，即：

，（29）

式中表示给定状态量；表示给定值。

通过“乘大数法”引入边界条件式（29）后，式（27）可改写为：

（30）

式中  I表示与维数相同的单位矩阵；表示维数为的列向量，其中n表示材料层数，k表示每层的DQM离散点个数，M表示有限元节点数。满足：

式中表示状态量在总体状态量矢量中的位置。

结合式（12），（15）和（26b），可知式（30）中的系数矩阵为频率项和轴力的函数。为使得方程（30）有非零解，系数矩阵应满足：

（32）

计算组合梁在无轴力作用下的自振频率，可通过设置式（12）中的轴力项，然后通过二分法求得满足式（32）的组合梁各阶自振频率。通过设定式（12）中的频率，可相应求得组合梁的屈曲荷载。

3 算  例

3.1 混凝土?木材组合梁

混凝土?木材组合梁由上层混凝土与下层木梁组成。上层混凝土的材料参数为，，下层木梁的材料参数为，，层间具有剪切刚度。表1展示了四种不同边界条件下，即两端简支（SS）、一端简支一端固支（SC）、两端固支（CC）、一段固支一段自由（CF），组合梁的前五阶自振频率，并将计算结果与Xu等 ［1］解析解的结果进行对比。

混凝土与木梁间通过剪力连接件相连接，连接件刚度决定了组合作用的强弱。当连接件刚度较小时，组合作用较弱，此时组合梁的固有频率较低；当连接件刚度逐渐增大时，组合作用随之增强，固有频率也相应提高。当连接件刚度很大时，层间近似完美连接，此时固有频率达到稳定的最大值。

图3反映了四种不同边界条件下，组合梁一阶频率随剪切刚度而变化的趋势图。不同支座条件下，一阶频率发生显著变化时对应的连接刚度区间不同。两端简支（SS）下，其频率敏感的连接刚度变化区间为，而在两端固支（CC）条件下，其频率敏感的连接刚度变化区间为。

3.2 波形钢腹板梁

基于二维分析模型，组合梁内部位移与应力的分布满足传递方程式（14），因而无须预先引入其分布的假设函数，适用于分析具有各向异性特征的波形钢腹板梁结构。参考胡霖远等［4］、Johnson等［14］ 以及Xu等［15］的工作，波形钢腹板可转换为等效正交各向异性材料，并将波形钢腹板梁（Corrugated Steel Web Beams， CSWB）等效为单位宽度下的矩形三层组合结构。各层的材料参数如下：混凝土顶板的等效弹性模量为，泊松比为，密度为。波形钢腹板的等效弹性模量为，，；泊松比为，；密度为。混凝土底板的等效弹性模量为，泊松比为，密度为。波形钢腹板梁长，高，其中混凝土顶板与底板厚度均为0.25 m。

基于式（32）求得的各阶自振频率结果，与Zig?zag理论（Zig?zag Beam Theory， ZBT）［4］和有限元数值方法计算所得的结果对比如表2～5所示。表中f表示自振频率，括号内表示该方法与有限元结果的相对误差（%）。有限元结果来自于ABAQUS软件使用八节点平面应力单元的计算结果，上、下混凝土板的单元尺寸均为0.025 m×0.450 m，波形钢腹板层的单元尺寸为0.200 m×0.450 m，在建立波形钢腹板梁的平面应力有限元模型后计算其动力特性。

不同边界条件下，基于此方法的计算结果与有限元结果吻合度较高，各阶频率值均与有限元结果接近。而Zig?zag理论的结果在低阶频率上精度较高；随着频率阶数的升高，其精度有所下降。这是由于随着振型阶数升高，组合梁内的剪切变形呈现出更显著的非线性特征，因此使用Zig?zag理论所得结果的误差将逐渐增大。基于二维模型的传递方程式（14），无须引入变形和应力的预先分布假设，因此对于各阶频率均能保持较低的误差水平。为了保证总体精度，本算例统一选用较大的DQM离散点数（k=20）并在材料界面附近加密布置离散点，偏多偏密的离散点布置使得低阶频率的误差反而偏大。

3.3 钢?混组合梁

钢?混组合梁在桥梁工程中应用广泛。通过在界面处布置抗剪连接件，钢?混组合梁将钢材与混凝土进行组合，发挥其各自优势。考虑一跨径的钢?混凝土试验组合梁。钢梁为的标准H型钢，容重，弹性模量和泊松比分别为和。混凝土厚，宽，容重，弹性模量和泊松比分别为和。混凝土与钢梁之间通过栓钉（直径）连接，并沿跨径方向每间隔0.2 m布置一排，每排2个，产生的界面抗剪强度。参考胡霖远等［4］和Xu等［15］的工作，可将该组合梁等效为单位宽度下的矩形五层组合结构。对于五层结构的分析，体现了该方法可方便地处理任意多层结构的特点。

由于栓钉沿梁长方向为间隔布置，采用两种方式对其进行模拟。方式（Ⅰ）首先计算每排栓钉处的界面抗剪强度：

式中表示梁长方向的栓钉排数；表示等效界面的宽度，即。求得抗剪强度后，可按式（1）求得栓钉位置处虚拟粘接层的剪切模量。而其他无栓钉处的剪切模量则取为一较小值（）。在进行有限元离散时，通过设定单元间不同的剪切模量，体现了该方法可处理材料特性沿梁长变化的优势。方式（Ⅱ）则不考虑栓钉的间隔布置，由式（1）直接求得对应于的全梁长等效剪切模量。

表6给出了两种计算方式下组合梁的固有频率计算结果，并与有限元数值结果以及张云龙等 ［19］的室内试验结果进行对比。有限元数值计算基于ANSYS软件，工字钢与混凝土使用实体单元（SOLID65）进行离散，栓钉使用梁单元 （BEAM188） 模拟。对比表明，采用方式（Ⅰ）的计算值与试验实测结果吻合度较高，误差仅为0.64%。而不考虑栓钉间隔布置的全梁长等效方式（Ⅱ），其频率计算结果较试验实测值偏低。这是由于将抗剪强度进行全梁长均匀化等效，降低了端部处的连接刚度，使得组合梁的固有频率下降。

4 结  论

基于组合梁统一的二维分析模型，导出了哈密顿体系下组合梁动力分析的混合能量泛函。以包含频率函数的节点位移及其能量对偶的应力分量为单元节点未知量，使用有限元离散建立了状态空间法下的混合有限元控制方程。状态空间法建立了状态量在传递方向上解析的传递关系，将复杂的多层结构界面连接关系转换为简单的矩阵相乘，对分析任意多层组合梁具有显著优势。而通过引入梁长方向的有限元离散，可方便地处理材料特性沿梁长变化的问题。最后，通过在传递方向上引入DQM方法，可有效防止由于引入有限元离散而导致的待解矩阵过大产生的数值稳定性问题。

部分连接的组合梁以及具有正交各向异性特征的波形钢腹板梁是工程中常见的组合梁形式，该方法准确分析了其在多种边界条件下的自振频率。对比分析波形钢腹板梁常用的Zig?zag理论，由于基于二维分析模型，该方法对高阶频率的计算十分准确，且随着频率的升高，其优势更加明显。而对于二层、三层及五层组合结构的分析，表明分析过程具有通用性，并可方便地适用于任意多层组合结构的自由振动分析。通过设置梁长方向上有限元单元之间不同的材料参数，可精细化处理材料特性沿梁长方向非均匀分布的情形，使计算模型及其结果与实际情况更加符合。

参考文献：

［1］Xu R， Wu Y. Static， dynamic， and buckling analysis of partial interaction composite members using Timoshenko's beam theory［J］. International Journal of Mechanical Sciences， 2007， 49（10）： 1139-1155.

［2］Vo T P， Thai H T. Vibration and buckling of composite beams using refined shear deformation theory ［J］. International Journal of Mechanical Sciences， 2012， 62（1）： 67-76.

［3］Carrera E， Filippi M， Zappino E. Free vibration analysis of laminated beam by polynomial， trigonometric， exponential and Zig-zag theories［J］. Journal of Composite Materials， 2014， 48（19）： 2299-2316.

［4］胡霖远，陈伟球，张治成，等. 基于 Zig-zag 理论的波形钢腹板梁自由振动分析［J］. 浙江大学学报（工学版），2019，53（3）： 503-511.

HU Lin-yuan， CHEN Wei-qiu， ZHANG Zhi-cheng， et al. Free vibration analysis of concrete beams with corrugated steel webs based on Zig-zag theory ［J］. Journal of Zhejiang University （Engineering Science）， 2019， 53（3）： 503-511.

［5］Shen X， Chen W， Wu Y， et al. Dynamic analysis of partial-interaction composite beams［J］. Composites Science and Technology， 2011， 71（10）： 1286-1294.

［6］沈旭栋，陈伟球，徐荣桥. 有轴力的部分作用组合梁的动力分析［J］. 振动工程学报，2012，25（5）： 514-521.

SHEN Xu-dong， CHEN Wei-qiu， XU Rong-qiao. Dynamic analysis of partial-interactive composite beam with axial force［J］. Journal of Vibration Engineering， 2012， 25（5）： 514-521.

［7］侯忠明，夏禾，王元清，等. 钢-混凝土组合梁动力折减系数研究 ［J］. 振动与冲击，2015，34（4）： 74-81.

HOU Zhong-ming， XIA He， WANG Yuan-qing， et al.  Dynamic reduction coefficients for a steel-concrete composite beam ［J］. Journal of Vibration and Shock， 2015， 34（4）： 74-81.

［8］陈玉骥，罗旗帜. 钢-混凝土组合梁考虑剪力滞和滑移的一阶自振频率［J］. 烟台大学学报（自然科学与工程版），2015，28（1）： 49-53.

CHEN Yu-ji， LUO Qi-zhi. First order natural vibration frequency for steel-concrete composite beams considering influence of shear lags and slips［J］. Journal of Yantai University （Natural Science and Engineering Edition）， 2015， 28（1）： 49-53.

［9］Vo T P， Lee J. Free vibration of axially loaded thin-walled composite Timoshenko beams［J］. Archive of Applied Mechanics， 2011， 81（9）： 1165-1180.

［10］Chakrabarti A， Sheikh A H， Griffith M， et al. Dynamic response of composite beams with partial shear interaction using a higher-order beam theory［J］. Journal of Structural Engineering， 2013， 139（1）： 47-56.

［11］Chalak H D， Chakrabarti A， Iqbal M A， et al. Vibration of laminated sandwich beams having soft core［J］. Journal of Vibration and Control， 2012， 18（10）： 1422-1435.

［12］Chen W Q， Lü C F， Bian Z G. Elasticity solution for free vibration of laminated beams［J］. Composite Structures， 2003， 62（1）： 75-82.

［13］Xu R， Wu Y F. Free vibration and buckling of composite beams with interlayer slip by two-dimensional theory ［J］. Journal of Sound and Vibration， 2008， 313（3-5）： 875-890.

［14］Johnson R P， Cafolla J， Bernard C. Corrugated webs in plate girders for bridges ［J］. Proceedings of the Institution of Civil Engineers-Structures and Buildings， 1997， 122（2）： 157-164.

［15］Xu R， Wu Y F. Two-dimensional analytical solutions of simply supported composite beams with interlayer slips［J］. International Journal of Solids and Structures， 2007， 44（1）： 165-175.

［16］Courant R， Hilbert D. Methods of Mathematical Physics： Partial Differential Equations［M］. Singapore： John Wiley & Sons， 2008： 32-34.

［17］Bellman R E， Kashef B G， Casti J. Differential quadrature： A technique for the rapid solution of nonlinear partial differential equations［J］. Journal of Computational Physics， 1972， 10： 40-52.

［18］Shu C. Differential Quadrature and Its Application in Engineering［M］. London： Springer Science & Business Media， 2012： 29-31.

［19］张云龙，郭阳阳，王静，等. 钢-混凝土组合梁的固有频率及其振型［J］. 吉林大学学报（工学版），2020，50（2）： 581-588.

ZHANG Yun-long， GUO Yang-yang， WANG Jing， et al. Natural frequency and mode of vibration of steel?concrete composite beam［J］. Journal of Jilin University （Engineering and Technology Edition）， 2020， 50（2）： 581-588.

A mixed finite element method for free vibration analysis of multilayer composite beams

Abstract： In this paper， a two-dimensional analytical model for composite beams is first proposed through the equivalent transformation of the cross-section. Based on the mixed variational principle， the dynamic state equations are derived through finite element meshing and interpolation along the length of the beam， with frequency contained nodal displacements and their energy-conjugated stresses as element nodal variables. The differential quadrature method （DQM） is introduced to discretize the equations along the height of the beam， and natural frequencies of composite beams under different axial forces and boundary conditions are obtained. This method was verified by numerical examples about natural frequencies of three beams， i.e. a concrete-wood composite beam， a concrete beam with a corrugated steel web and steel-concrete composite beam. Since the proposed method is based on the two-dimensional theory， it can provide benchmarks for beam theories and error analyses.

Key words： free vibration；composite beam；mixed finite element；DQM；state space method