一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法与流程

本发明涉及桥梁数据计算技术领域，尤其涉及一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法。

背景技术：

大跨度桥梁结构通常采用薄壁箱梁的截面形式，薄壁箱梁宽跨比较大，剪力滞后效应显著，同时由于薄壁箱梁结构自重较轻，承受的活载所占的荷载比重较大，因此其动力特性复杂。目前用于考虑剪力滞后效应的薄壁箱梁结构动力特性的分析方法主要有弹性理论法、比拟杆法、能量变分法以及有限单元法。弹性理论解法以经典的弹性理论为基础，对于简单的力学模型问题，能够获得较为精确的结果。对于大型复杂结构的动力分析，会面临计算过程繁琐、计算工作量大的问题，这使得弹性理论解法在解决工程实际问题中受到了一定的限制。比拟杆法由于进行了若干基本假设，将力学模型简化，因此仅适用于一般等截面箱梁，对某些复杂结构考虑剪力滞后效应的动力分析仍存在一定的困难。能量变分法在箱梁结构的初步设计阶段发挥了良好的效果，但该法只适应于等截面箱梁，对于变截面箱梁剪力滞后效应的求解还有待于进一步研究。另外，该法将翼板结构按照平面应力问题进行分析，这一假定对于翼板自由端受力形态的求解存在较大的误差，同时不同的位移模式假定对分析结果也具有一定的影响，如何合理地选取翼板轴向位移分布函数还有待进一步研究。伴随着计算机性能的不断提升，有限单元法逐渐成为薄壁结构动力特性分析的主要方法。有限单元法的主要优点是求解速度快，计算精度较高，但也存在一定的问题：首先，当节点自由度和单元位移参数较多时，有限单元法的计算效率会大大降低；其次，通常借助于有限单元程序进行考虑剪力滞后效应的结构动力特性分析，程序中节点的位移插值函数是基于拉格朗日(lagrange)或者赫米特(hermite)插值多项式进行构建，随着对剪力滞后问题的不断深入研究，实测结果显示上述插值多项式并不能精确描述薄壁钢箱梁横截面的轴向位移分布规律。第三，有限元法需要计算机具有强大的存储空间且分析计算会耗费较长的时间，不便于工程问题计算采用。

技术实现要素：

本发明用以克服现有技术的不足，提供一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法。

本发明的目的可以通过采用如下的技术措施来实现，提供一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法，其步骤包括：就薄壁箱梁的顶板、腹板和底板，分别获取刚度矩阵，叠加后得到薄壁箱梁的总刚度矩阵；获取薄壁箱梁顶板、腹板和底板的阻尼矩阵和质量矩阵，叠加后得到薄壁箱梁的总阻尼矩阵和总质量矩阵；将薄壁箱梁的总刚度矩阵、总阻尼矩阵和总质量矩阵代入有限单元动力平衡方程，计算得到薄壁箱梁的自振频率；通过薄壁箱梁的自振频率，与预设的梁箱自振频率阈值进行比较，以对薄壁箱梁结构动力特性进行分析。

其中，获取薄壁箱梁顶板刚度矩阵的步骤包括：

在顶板上选取4个节点，设定其位移参数为基本位移参数，

δ＝[δ1δ2δ3δ4]^t(1)

其中，δi＝[uiviωiθxiθyi](i＝1，2，3，4)；ui、vi、ωi、θxi、θyi分别为每个节点位移的轴向分量、横向分量、竖向分量、绕x轴的转动分量和绕y轴的转动分量；

分别设置节点位移的各个分量的差值函数；

其中，取一个符合实测结果的3次多项式作为顶板节点的轴向位移插值函数，其表达式如下：

其中，为局部坐标系下的节点横向坐标，方向同y轴方向，k为顶板单元宽度；

顶板节点竖向和转动位移插值函数取拉格朗日插值函数，其表达式如下：

n＝[ni，j]5×20(3)

顶板的刚度矩阵由平面应力状态下顶板的刚度矩阵和薄板小挠度弯曲状态下顶板的刚度矩阵叠加而成，即：

其中，t、d分别为顶板单元厚度、顶板单元长度；bp、dp、bb分别为平面应力状态下单元的应变矩阵、弹性矩阵以及薄板小挠度弯曲状态下单元的应变矩阵。

其中，获取薄壁箱梁腹板和底板刚度矩阵的步骤包括：

取一次多项式和三次多项式分别作为腹板节点的轴向位移插值函数与竖向位移插值函数，如下所示：

m＝[1-ζζ](6)

n＝[1-3ζ²+2ζ³(ζ-2ζ²+ζ²)d3ζ²-2ζ³(3ζ³-ζ²)d](7)

其中，ζ＝x/d，为局部坐标系下的节点轴向坐标(方向同x轴方向)；

腹板节点位移参数根据腹板单元与顶板单元的变形协调关系导出，根据公式(1)-(3)，选取的腹板节点(7)和(8)的轴向位移参数表示为：

腹板节点(7)和(8)转角位移参数可表示为：

腹板节点(7)和(8)竖向位移参数表示为：

其中，s为翼板宽度；a为薄壁箱梁上口宽度。

设定薄壁箱梁的其中一端为m端，m端左、右腹板绕y轴的转角位移参数与节点(7)、(8)绕y轴的转角相同，可分别表示为：

m端底板质心处的转角位移参数由底板两位移节点(5)和(6)绕y轴的转角位移参数线性内插得到，即：

腹板的左、右板件质心处的轴向位移参数表示为：

底板位移节点(5)和(6)的轴向位移参数表示为：

u5＝u8-hθ1cm(19)

u6＝u7-hθrcm(20)

其中，h为腹板高度。

底板质心处的轴向位移参数可根据节点(5)和(6)的轴向位移参数线性内插得到，即：

m端腹板左、右板件质心处的竖向位移与顶板节点(7)、(8)的竖向位移参数数值相同，即：

ω1cm＝ω8(22)

ωrcm＝ω7(23)

同理，m端底板质心处的竖向位移由节点(5)、(6)的竖向位移线性内插得到，由于节点(5)、(6)的竖向位移与腹板左、右板件质心处的竖向位移相同，因此，底板子单元m端质心处的竖向位移可表示为：

对于腹板左板件，假定分别为板件质心处的轴向位移列阵和竖向位移列阵：

轴向位移参数与竖向位移参数表示为：

其中，a为和δ的转换矩阵，由公式(1)、(9)、(14)、(17)、(25)联立求解得到；b为和δ的转换矩阵，由公式(1)、(11)、(14)、(22)、(26)求解得到；下面给出矩阵a和b各元素：

其余元素均为0；

对于腹板右板件，假定为板件质心处的轴向位移列阵和竖向位移列阵：

则板件轴向位移参数与竖向位移参数可表示为：

其中，c为和δ的转换矩阵，可由公式(1)，(8)，(15)，(18)，(29)联立求解得到；d为和δ的转换矩阵，可由公式(1)，(10)，(15)，(23)，(30)联立求解得到；

假定分别为底板质心处的轴向位移列阵和竖向位移列阵：

底板质心处的轴向位移和竖向位移分别表示为：

其中，e为和δ的转换矩阵，由公式(1)，(8)，(9)，(14)，(15)，(19)-(21)，(33)求解得到；f为和δ的转换矩阵，由公式(1)，(12)-(16)，(22)-(24)，(34)求解得到；

薄壁箱梁总应变能方程表示为：

∏＝∏r+∏w+∏f-f^etδ(37)

其中，∏r为顶板应变能、∏w、∏f分别为腹板、底板的应变能，可以分别表示为以下形式：

其中，ea1、ear、eab分别为腹板左、右板件，底板的抗压刚度；eiy1、eiyr、eiyb分别为腹板左、右板件，底板的抗弯刚度；分别为腹板、底板的刚度矩阵；f^et为外荷载列阵；

根据最小势能原理，对公式(39)变分并令结果为0，可得到腹板、底板的刚度矩阵分别如下式所示：

其中，计算得到薄壁箱梁的自振频率的步骤包括：

振动系统的有限单元动力平衡方程表示为如下形式：

其中，m^e，c^e，k^e分别为质量矩阵、阻尼矩阵、刚度矩阵；δ、p^e分别为相应的加速度列阵、速度列阵、位移列阵；

采用一致质量矩阵分析结构的振动特性，假定材料的质量密度为ρ，则顶板子单元的一致质量矩阵可表示为：

其中，au为顶板面积；

腹板一致质量矩阵表示为：

底板一致质量矩阵表示为：

其中，a1，ar，af分别为腹板左、右板件，底板子单元面积；

总质量矩阵可通过叠加获得；总刚度矩阵可通过叠加获得；

作用于顶板节点上的外荷载以集中荷载和均布荷载的形式作用于结构上，则由此两种荷载产生的荷载列阵可以表示为：

有限单元动力平衡方程中基本参数均已推导完成，代入求解即可得到薄壁箱梁的自振频率。

其中，根据计算得到的薄壁箱梁的自振频率与预设的城市人行天桥自振频率阈值进行比较，自振频率阈值为3hz，当计算得到的自振频率大于3hz，薄壁箱梁不会与行人载荷产生共振，若计算得到的自振频率小于3hz时，薄壁箱梁可能与行人载荷产生共振。

区别于现有技术，本发明的考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法在传统有限单元法的基础上，通过减少单元基本位移参数和节点自由度提高了结构动力分析效率，并且引入一个更加符合实际变形规律的轴向位移分布函数以提高结构动力分析精度。该方法能有效解决目前对于考虑剪力滞后效应的薄壁箱梁结构动力分析过程中存在的建模工作量大、计算效率低、计算精度差等问题。

附图说明

图1是本发明提供的一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法的流程示意图。

图2是本发明提供的一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法所涉及的薄壁箱梁的结构示意图。

图3是本发明提供的一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法中具体实施例所涉及的薄壁箱梁的结构示意图。

图4是本发明提供的一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法中薄壁箱梁的剖面结构示意图。

具体实施方式

下面结合具体实施方式对本发明的技术方案作进一步更详细的描述。显然，所描述的实施例仅仅是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

参阅图1，图1是本发明提供的一种考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法的流程示意图。该方法的步骤包括：

s110：就薄壁箱梁的顶板、腹板和底板，分别获取刚度矩阵，叠加后得到薄壁箱梁的总刚度矩阵。薄壁箱梁的结构如图2所示。

获取薄壁箱梁顶板刚度矩阵的步骤包括：

在顶板上选取4个节点，设定其位移参数为基本位移参数，四个节点为图2中的节点1-4。

δ＝[δ1δ2δ3δ4]^t(1)

其中，δi＝[uiviωiθxiθyi](i＝1，2，3，4)；ui、vi、ωi、θxi、θyi分别为每个节点位移的轴向分量、横向分量、竖向分量、绕x轴的转动分量和绕y轴的转动分量；

分别设置节点位移的各个分量的差值函数；

其中，取一个符合实测结果的3次多项式作为顶板节点的轴向位移插值函数，其表达式如下：

其中，为局部坐标系下的节点横向坐标，方向同y轴方向，k为顶板单元宽度；

顶板节点竖向和转动位移插值函数取拉格朗日插值函数，其表达式如下：

n＝[ni，j]5×20(3)

顶板的刚度矩阵由平面应力状态下顶板的刚度矩阵和薄板小挠度弯曲状态下顶板的刚度矩阵叠加而成，即：

其中，t、d分别为顶板单元厚度、顶板单元长度；bp、dp、bb分别为平面应力状态下单元的应变矩阵、弹性矩阵以及薄板小挠度弯曲状态下单元的应变矩阵。

获取薄壁箱梁腹板和底板刚度矩阵的步骤包括：

取一次多项式和三次多项式分别作为腹板节点的轴向位移插值函数与竖向位移插值函数，如下所示：

m＝[1-ζζ](6)

n＝[1-3ζ²+2ζ³(ζ-2ζ²+ζ²)d3ζ²-2ζ³(3ζ³-ζ²)d](7)

其中，ζ＝x/d，为局部坐标系下的节点轴向坐标(方向同x轴方向)；

腹板节点位移参数根据腹板单元与顶板单元的变形协调关系导出，根据公式(1)-(3)，选取的腹板节点(7)和(8)为图2中腹板与顶板连接边缘位置的点，轴向位移参数表示为：

腹板节点(7)和(8)转角位移参数可表示为：

腹板节点(7)和(8)竖向位移参数表示为：

其中，s为翼板宽度；a为薄壁箱梁上口宽度。

m端左、右腹板绕y轴的转角位移参数与节点(7)、(8)绕y轴的转角相同，可分别表示为：

m端底板质心处的转角位移参数由底板两位移节点(5)和(6)绕y轴的转角位移参数线性内插得到，即：

腹板的左、右板件质心处的轴向位移参数表示为：

底板位移节点(5)和(6)的轴向位移参数表示为：

u5＝u8-hθ1cm(19)

u6＝u7-hθrcm(20)

其中，h为腹板高度。

底板质心处的轴向位移参数可根据节点(5)和(6)的轴向位移参数线性内插得到，即：

m端腹板左、右板件质心处的竖向位移与顶板节点(7)、(8)的竖向位移参数数值相同，即：

ω1cm＝ω8(22)

ωrcm＝ω7(23)

同理，m端底板质心处的竖向位移由节点(5)、(6)的竖向位移线性内插得到，由于节点(5)、(6)的竖向位移与腹板左、右板件质心处的竖向位移相同，因此，底板子单元m端质心处的竖向位移可表示为：

对于腹板左板件，假定分别为板件质心处的轴向位移列阵和竖向位移列阵：

轴向位移参数与竖向位移参数表示为：

其中，a为和δ的转换矩阵，由公式(1)、(9)、(14)、(17)、(25)联立求解得到；b为和δ的转换矩阵，由公式(1)、(11)、(14)、(22)、(26)求解得到；下面给出矩阵a和b各元素：

其余元素均为0；

对于腹板右板件，假定为板件质心处的轴向位移列阵和竖向位移列阵：

则板件轴向位移参数与竖向位移参数可表示为：

其中，c为和δ的转换矩阵，可由公式(1)，(8)，(15)，(18)，(29)联立求解得到；d为和δ的转换矩阵，可由公式(1)，(10)，(15)，(23)，(30)联立求解得到；

假定分别为底板质心处的轴向位移列阵和竖向位移列阵：

底板质心处的轴向位移和竖向位移分别表示为：

其中，e为和δ的转换矩阵，由公式(1)，(8)，(9)，(14)，

其中，e为和δ的转换矩阵，由公式(1)，(8)，(9)，(14)，(15)，(19)-(21)，(33)求解得到；f为和δ的转换矩阵，由公式(1)，(12)-(16)，(22)-(24)，(34)求解得到；

薄壁箱梁总应变能方程表示为：

∏＝∏r+∏w+∏f-f^etδ(37)

其中，∏r为顶板应变能、∏w、∏f分别为腹板、底板的应变能，分别表示为以下形式：

其中，ea1、ear、eab分别为腹板左、右板件，底板的抗压刚度；eiy1、eiyr、eiyb分别为腹板左、右板件，底板的抗弯刚度；分别为腹板、底板的刚度矩阵；f^et为外荷载列阵；

根据最小势能原理，对公式(39)变分并令结果为0，可得到腹板、底板的刚度矩阵分别如下式所示：

s120：获取薄壁箱梁顶板、腹板和底板的阻尼矩阵和质量矩阵，叠加后得到薄壁箱梁的总阻尼矩阵和总质量矩阵。

采用一致质量矩阵分析结构的振动特性，假定材料的质量密度为ρ，则顶板子单元的一致质量矩阵可表示为：

其中，au为顶板面积；

腹板一致质量矩阵表示为：

底板一致质量矩阵表示为：

其中，a1，ar，af分别为腹板左、右板件，底板子单元面积；

总质量矩阵可通过叠加获得；总刚度矩阵可通过叠加获得；

作用于顶板节点上的外荷载以集中荷载和均布荷载的形式作用于结构上，则由此两种荷载产生的荷载列阵可以表示为：

s130：将薄壁箱梁的总刚度矩阵、总阻尼矩阵和总质量矩阵代入有限单元动力平衡方程，计算得到薄壁箱梁的自振频率。

对于一个振动系统而言，有限单元动力平衡方程可表示为如下形式：

其中，m^e，c^e，k^e分别为质量矩阵、阻尼矩阵、刚度矩阵；δ，p^e分别为相应的加速度列阵，速度列阵、位移列阵。

通过前两步骤，有限单元动力平衡方程中基本参数均已推导完成，代入求解即可得到薄壁箱梁的自振频率。

s140：通过薄壁箱梁的自振频率，与预设的梁箱自振频率阈值进行比较，以对薄壁箱梁结构动力特性进行分析。

以人行天桥为例进行说明，人行天桥的主要活载荷为行人，行人在行走时有其步行频率，步行频率不论男女老幼，差别不大，一般在2hz左右。为避免主桥的固有自振频率与步行频率较接近而引起主桥振动及挠度过大，引起行人不适，甚至危及人形天桥安全，因此，《城市人行天桥与人行地道技术规范》第254条规定：“为避免共振，减少行人不安全感，天桥上部结构竖向自振频率不应小于3hz。”因此本发明以3hz作为经验阈值。

根据计算得到的薄壁箱梁的自振频率与预设的城市人行天桥自振频率阈值进行比较，自振频率阈值为3hz，当计算得到的自振频率大于3hz，薄壁箱梁不会与行人载荷产生共振，若计算得到的自振频率小于3hz时，薄壁箱梁可能与行人载荷产生共振。

下述为本发明的一个实施例：

一座主跨跨径为30米的简支混凝土薄壁箱梁桥，主梁截面形式为单箱单室，顶板翼缘板长度为3.55米，厚度为0.25米。底板宽度和厚度分别为7.1米和0.25米。腹板高度和厚度分别为2米和0.4米。实例立面图和横断面图如图3所示，剖面图如图4所示。结构几何尺寸和材料参数见表1。

表1：结构几何尺寸和材料参数

运用本发明的方法计算得到该工程实例考虑剪力滞后效应的前6阶自振频率，为了对比验证本方法的计算精度和计算效率，采用大型通用有限元软件ansys建立本工程实例的壳单元有限元模型，结构沿顺桥向划分为30个单元，单元采用shell181单元。将采用两种方法计算得到的结果进行对比分析，如表2所示。

表2实例前6阶自振频率对比

通过上表对比分析结果可以得出：运用本方法计算薄壁箱梁结构的自振频率，其结果与采用壳单元有限元模型计算得到的结果基本吻合，说明本文方法具有较高的计算精度；同时从自由度个数可以看出，运用本文方法建立模型，节点自由度总个数为310，采用通用有限元软件aanys建立模型，节点自由度总个数为2232，因此，本方法具有更高的计算效率。

区别于现有技术，本发明的考虑剪力滞后作用的薄壁箱梁结构动力特性分析方法在传统有限单元法的基础上，通过减少单元基本位移参数和节点自由度提高了结构动力分析效率，并且引入一个更加符合实际变形规律的轴向位移分布函数以提高结构动力分析精度。该方法能有效解决目前对于考虑剪力滞后效应的薄壁箱梁结构动力分析过程中存在的建模工作量大、计算效率低、计算精度差等问题。

以上仅为本发明的实施方式，并非因此限制本发明的专利范围，凡是利用本发明说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本发明的专利保护范围内。