PARAMETER FREEZING PRECISE EXPONENTIAL INTEGRATOR AND ITS APPLICATION IN NONLINEAR VEHICLE-BRIDGE COUPLED VIBRATION
-
摘要: 描述车桥耦合作用的基本问题是一个时变系统问题, 且很多工况下需考虑非线性特性, 使得该问题难以得到解析解, 甚至数值解也可能很复杂. 针对该问题的求解, 提出了一种参数冻结精细指数积分法, 将其应用于车桥耦合动力学模型的数值分析中. 该方法结合了精细积分和指数积分特点, 并将时变系数矩阵在每一积分步参数冻结, 用于获得系统振动响应的数值解. 考虑汽车轮胎与桥面的力和位移耦合关系、桥面沥青铺装层、桥梁材料黏弹性和几何非线性特性, 建立了车桥耦合动力学模型, 并应用参数冻结精细指数积分法对该模型进行了求解. 通过与近似解析解、辛Runge-Kutta算法以及经典的Newmark-β数值积分法计算结果进行对比, 验证了所提出方法计算结果的有效性和准确性. 在此基础上, 制作了缩尺车桥耦合系统模型, 测试了跨中挠度响应, 进一步验证了理论建模和所提算法的有效性和实用性. 通过数值计算分析了所提算法的数值特性, 结果表明: 提出的参数冻结精细指数积分法不仅可以处理时变、非线性问题, 且具有良好的数值计算精度和长时间数值稳定性; 由于精细积分的特点, 参数冻结精细指数积分法的计算时间步长可以取的较大, 可有效提高计算效率. 因此, 所提出的参数冻结精细指数积分法预期可成为求解车桥耦合动力学问题的一种新的高效算法.
-
关键词:
- 精细指数积分法 /
- 车桥耦合振动 /
- 参数冻结 /
- 非线性时变系统 /
- 辛Runge-Kutta算法
Abstract: The basic problem describing vehicle-bridge interaction is a time-varying system problem, which the nonlinear characteristics need to be considered in many working conditions. It is difficult to get the analytical solution to this problem, and even the numerical solution may be complex. In this paper, focus on the solution to this problem, a parameter freezing precise exponential integrator is proposed and its application to vehicle-bridge coupled problem is analyzed. This algorithm combines the characteristics of precise integration and exponential integration, and freezes the parameters of time-varying coefficient matrix at each integration step to obtain the numerical solutions of vibration responses of the system. Considering the coupling relationship of the force and displacement between the vehicle tire and bridge deck, the asphalt pavement of the bridge deck, the viscoelastic and geometric nonlinear characteristics of the bridge, the vehicle-bridge coupled dynamic model is established. The parameter freezing precise exponential integrator is used to solve this dynamic model. By comparing with the approximate analytical solution, the solution calculated by the symplectic Runge-Kutta algorithm and the classical Newmark-β numerical integration method, the validity and accuracy of the proposed method are verified. On this basis, a scaled vehicle-bridge coupling system model is made. The mid-span deflection response is tested, which further verifies the validity and practicability of the theoretical model and proposed algorithm. The numerical characteristics of the proposed algorithm are analyzed by numerical calculation. Numerical results show that the proposed parameter freezing precise exponential integrator can not only deal with time-varying and nonlinear problems, but also has good numerical calculation accuracy and long-time numerical stability. Due to the characteristics of the precise integration, the time step can be set to be larger for the parameter freezing precise exponential integrator, which can effectively improve the calculation efficiency. Therefore, the proposed parameter freezing precise exponential integrator is expected to be a new and efficient algorithm for solving vehicle-bridge coupling dynamics problem. -
引 言
在桥梁系统众多力学问题中, 车桥相互作用问题是一个长期的热点课题[1-6]. 其核心在于如何准确地预测车桥耦合振动行为. 当车辆在桥梁上运行时, 车辆和桥梁之间是相互作用和互相影响的, 这称之为车桥耦合振动. 车桥耦合振动行为的研究可为桥梁结构的设计、运营、维护和管理提供必要的理论基础、分析方法以及评估手段, 具有非常重要的工程应用价值[3]. 关于车桥耦合系统动力学问题的研究, 早期研究将移动的车辆对桥梁施加的力建模为移动载荷[7-8]. 在这种模型中, 移动的车辆和桥梁之间不存在动态耦合关系, 车辆载荷只是移动在桥面上的外部激励. 然而, 当车辆模型中考虑惯性力时, 桥梁和车辆在接触点处的响应是耦合的, 即构成车桥耦合系统[9]. 且由于车辆位置随时间的不断改变, 描述车桥耦合系统的动力学方程组具有时变特性[10-11]. 通常情况下, 即使是线性时变系统也很难得到精确解, 因此各种近似方法被应用于获得车桥耦合系统的动力学响应.
车桥耦合系统通常采用有限元法[12-16]或模态叠加法[17-19]进行建模和分析. 通过相应方法得到时变常微分动力学方程以后, 进一步采用直接积分法进行求解, 主要有线性加速度法、Wilson法和Newmark法等[20]. Zhu等[21]采用精细积分法计算了移动载荷通过连续梁桥时的动力学响应. Zhai[22]基于Newmark法发展了两种高效的显式数值积分方法, 提高了大规模计算时的效率和精度. Stoura等[23]提出了一种求解车桥耦合动力学模型的动态划分方法, 该方法如结合适当的数值格式可显著降低计算量. 处理各种时变系统动力学问题时, Newmark法[10]和Runge-Kutta法[24]是应用较多的数值方法. 此外, 控制问题中经典的时间参数冻结技术也被应用于慢时变系统动力学问题[25]. 最近, Ge等[26]针对时变动力学系统的计算问题, 提出了求解车桥耦合系统受迫振动的时间参数冻结方法, 该方法可较好地处理线性时变车桥耦合问题.
在车桥耦合问题的计算中, 为了数值计算的准确性, 在某些工况下需考虑非线性因素, 如桥梁跨度大, 当载荷作用下产生大变形时的几何非线性[27]. 当考虑非线性因素时, 描述车桥相互作用问题的动力学方程组具有典型的时变、非线性特性, 对其高效、精确的求解具有重要的意义.
1994年, 钟万勰[28]针对线性常系数微分方程问题的计算, 提出了精细积分法, 该方法在求解指数矩阵时用矩阵加法代替乘法, 巧妙地避免了舍入误差. 由于其优异的数值特性, 近年来精细积分法得到了广泛的研究和应用[29-32], 并发展出了增维精细积分法[33]、广义精细积分法[34]、辛精细积分法[35-36]和快速精细积分法[37]等方法. 针对半线性微分方程问题, 近年来发展了多种指数积分方法及其相应的计算格式[38], 在处理高振荡、刚性问题中得到了良好的应用[39-40]. 邓子辰等[41]结合精细积分法和指数积分法, 提出了精细指数积分法, 并将其应用于卫星编队飞行动力学问题中弱非线性方程的求解, 但并未考虑时变问题.
车桥耦合动力学问题解析求解难度大, 通常采用数值积分方法进行求解, 如经典的Newmark-β法、Runge-Kutta法. 但是这两种方法受稳定性影响较大, 数值积分步长需选取较小, 这会增加计算时间, 降低计算效率. 且这些算法本身具有一定的数值耗散性, 在长时间数值积分后可能会导致计算结果的失效. 因此, 本文针对上述问题, 将精细指数积分法与时间参数冻结技术相结合, 提出一种新的计算方法——参数冻结精细指数积分法, 并将其拓展应用于车桥耦合系统动力学问题的求解中. 充分发挥不同方法的特点和优势, 使其能够高效、精确地求解具有时变、非线性特性的车桥耦合动力学方程组, 预期可为车桥耦合问题的动力学仿真计算提供一种新的方法.
1. 考虑桥面铺装层的车桥耦合模型
考虑汽车通过跨度为L的桥梁, 建模时桥梁采用包含铺装层和桥体的层合梁模型, 汽车采用车桥耦合振动问题建模中常用的两自由度1/4汽车悬架模型[10,24,42]. 车桥耦合模型示意图如图1所示.
考虑汽车自桥梁左端向右以匀速v行驶, 1/4车模型有两个描述垂向振动的动力学方程[10]
$$ {m_2}{\ddot z_2} + {k_2}\left( {{z_2} - {z_1}} \right) + {c_2}\left( {{{\dot z}_2} - {{\dot z}_1}} \right) = 0 $$ (1) $$ \begin{split} & {m_1}{{\ddot z}_1} + {k_2}\left( {{z_1} - {z_2}} \right) + {c_2}\left( {{{\dot z}_1} - {{\dot z}_2}} \right) + \\ &\qquad {k_1}\left( {{z_1} - {w_1}} \right) + {c_1}\left( {{{\dot z}_1} - {{\dot w}_1}} \right) = 0\end{split} $$ (2) 其中$ {m_1} $和$ {m_2} $分别为非簧载质量和簧载质量, $ {z_1} $和$ {z_2} $分别为非簧载质量位移和簧载质量位移, $ {k_1} $和$ {k_2} $分别为轮胎和悬架刚度系数, $ {c_1} $和$ {c_2} $分别为轮胎和悬架阻尼系数, $ {w_1} $为车桥耦合振动引起的车轮与桥面接触点处桥面的二次位移激励值.
车轮与桥面接触的动态轮胎力可表示为
$$ {F_1} = {k_1}\left( {{z_1} - {w_1}} \right) + {c_1}\left( {{{\dot z}_1} - {{\dot w}_1}} \right) - \left( {{m_1} + {m_2}} \right)g $$ (3) 基于达朗贝尔原理, 考虑图1中的具体层合梁模型, 梁的振动方程可给出如下
$$ \rho A\ddot w - M'' - F\left( {x,t} \right) = 0 $$ (4) $$ \rho = \frac{{{\rho _1}\left( {{h_t} - {h_m}} \right) + {\rho _2}\left( {{h_m} - {h_b}} \right)}}{{{h_t} - {h_b}}} $$ (5) $$ M = \int_A {{\sigma _x}z} {\text{d}}A $$ (6) 其中$ \rho $为梁材料等效密度, $ {\rho _1} $和$ {\rho _2} $分别为桥面铺装层和桥体的密度, $ {h_t},{h_m}和{h_b} $分别表示铺装层上表面、铺装层与桥体连接面和桥体下表面的位置坐标, $ A $为梁横截面积, $ w $表示梁的横向位移场, M为弯矩, $ {\sigma _x} $为梁轴向应力, $ F\left( {x,t} \right) $为梁横向外载荷, 符号(.)表示对时间t求导, 符号(')表示对空间x求导.
当式(6)做积分运算时, z值0点在层合梁中性轴处, 即图1中x轴位置. 因此, 需根据桥面铺装层和桥体具体参数确定$ {h_t},{h_m}和{h_b} $的坐标. 假设铺装层厚度为$ {h_1} $, 桥体厚度为$ {h_2} $, 桥面铺装层和桥体的弹性模量比值为
$$ K = {E_1}/{E_2} $$ (7) 其中$ {E_1} $和$ {E_2} $分别表示桥面铺装层和桥体的弹性模量. 中性轴距离上边缘距离$ {h_{pt}} $和下边缘的距离$ {h_{pb}} $以及$ {h_t},{h_m},{h_b} $的坐标分别为
$$ {h_{dt}} = \left[ {K{h_1}/2 + {h_2}\left( {{h_1} + {h_2}/2} \right)} \right]/\left( {K{h_1} + {h_2}} \right) $$ (8) $$ {h_{db}} = {h_1} + {h_2} - {h_{dt}} $$ (9) $$ {h_t} = {h_{dt}},{h_m} = {h_{dt}} - {h_1},{h_b} = - {h_{db}} $$ (10) 桥梁结构自身具有一定的对外界能量耗散和减振能力, 其对振动能量的耗散能力可用黏性阻尼系数来体现[43]. 本文桥梁建模考虑了Kelvin黏弹性理论以及几何非线性特性, 梁轴向应力−应变关系和应变−位移关系可给出如下
$$\qquad\qquad {\sigma _x} = E\left( {{\varepsilon _x} + \eta {{\dot \varepsilon }_x}} \right) $$ (11) $$ \qquad\qquad {\varepsilon _x} = - zw'' + \frac{1}{2}{\left( {w'} \right)^2} $$ (12) 其中$ {\varepsilon _x} $为梁轴向应变, $ \eta $为黏性阻尼系数.
将式(3)、式(6)、式(11)和式(12)代入式(4)可以得到
$$ \begin{split} & \rho A\ddot w + {\xi _1}w{''''} - 2{\xi _2}\left[ {{{\left( {w{''}} \right)}^2} + w{'}w{'''}} \right]- \\ &\qquad {\xi _3}\left( {w{'''}\dot w{'} + 2w{''}\dot w{''} + w{'}\dot w{'''}} \right) + {\xi _4}\dot w{''''} = \\ &\qquad {F_1}\delta \left( {x - vt} \right) \end{split} $$ (13) $$\left. \begin{aligned} & {\xi _1} = \frac{b}{3}\left[ {{E_2}\left( { - {h_b^3} + {h_m^3}} \right) + {E_1}\left( { - {h_m^3} + {h_t^3}} \right)} \right] \\ & {\xi _2} = \frac{b}{4}\left[ {{E_2}\left( { - {h_b^2} + {h_m^2}} \right) + {E_1}\left( { - {h_m^2} + {h_t^2}} \right)} \right] \\ & {\xi _3} = \frac{b}{2}\left[ {{E_2}{\eta _2}\left( { - {h_b^2} + {h_m^2}} \right) + {E_1}{\eta _1}\left( { - {h_m^2} + {h_t^2}} \right)} \right] \\ & {\xi _4} = \frac{b}{3}\left[ {{E_2}{\eta _2}\left( { - {h_b^3} + {h_m^3}} \right) + {E_1}{\eta _1}\left( { - {h_m^3} + {h_t^3}} \right)} \right] \end{aligned}\right\} $$ (14) 其中$ \delta $为狄拉克函数, 满足关系
$$ \delta (x - {x_0}) = \left\{ \begin{split} & {1,\quad \quad x = {x_0}} \\ & {0,\quad \quad x \ne {x_0}} \end{split} \right. $$ (15) b为梁宽度, $ {\eta }_{1}, {\eta }_{2} $分别表示铺装层和桥体的黏性阻尼系数. 式(13)中如不考虑梁几何非线性, 则$ {\xi _2} $和$ {\xi _3} $等于0, 如不考虑梁的黏弹性, 则$ {\xi _4} $等于0.
采用伽辽金法可将偏微分方程(13)进行简化, 本文考虑梁两端为简支边界条件, 则梁横向位移函数可表示为
$$ w(x,t) = \sum\limits_{i = 1}^N {\sin \left( {i\text{π} x/L} \right) \cdot } {q_i}(t) $$ (16) 其中$ \sin (i\text{π} x/L) $为梁简支边界条件下的第i阶模态函数, $ {q_i}(t) $为梁振动的第i阶广义位移. 对于受移动载荷作用的简支梁桥, 其一阶模态函数对横向振动起决定性作用[8]. 因此, 本文选取了桥梁一阶模态进行研究, 即N = 1. 应用伽辽金法, 将式(16)代入式(13), 乘以第1阶模态函数并沿梁长积分, 则可给出描述梁横向振动的非线性常微分方程如下
$$ \begin{split} & {\psi _1}{{\ddot q}_1} + {\psi _2}{q_1} - {\psi _3}q_1^2 - {\psi _4}{q_1}{{\dot q}_1} + {\psi _5}{{\dot q}_1} = \\ &\qquad - {\phi _1}\left[ {{k_1}\left( {{\phi _1}{q_1} - {z_1}} \right) + {c_1}\left( {{\phi _1}{{\dot q}_1} - {{\dot z}_1}} \right) + \left( {{m_1} + {m_2}} \right)g} \right] \end{split} $$ (17) 其中
$$\left. \begin{split} & {\psi _1} = \frac{{\rho AL}}{2},{\psi _2} = \frac{{{\text{π} ^4}{\xi _1}}}{{2{L^3}}},{\psi _3} = \frac{{4{\text{π} ^3}{\xi _2}}}{{3{L^3}}} \\ & {\psi _4} = \frac{{4{\text{π} ^3}{\xi _3}}}{{3{L^3}}},{\psi _5} = \frac{{{\text{π} ^4}{\xi _4}}}{{2{L^3}}},{\phi _1} = \sin \frac{{\text{π} vt}}{L} \end{split} \right\}$$ (18) 整理式(1)、式(2)和式(17), 即为本文中描述车桥耦合振动的动力学方程组, 其中广义位移向量为$ {\boldsymbol{X}} = {\left[ {{z_1},{z_2},{q_1}} \right]^{\text{T}}} $
$$\left.\begin{split} & {m_2}{{\ddot z}_2} + {k_2}\left( {{z_2} - {z_1}} \right) + {c_2}\left( {{{\dot z}_2} - {{\dot z}_1}} \right) = 0 \\ & {m_1}{{\ddot z}_1} + {k_2}\left( {{z_1} - {z_2}} \right) + {c_2}\left( {{{\dot z}_1} - {{\dot z}_2}} \right) + \\ &\qquad {k_1}\left( {{z_1} - {\phi _1}{q_1}} \right) + {c_1}\left( {{{\dot z}_1} - {\phi _1}{{\dot q}_1}} \right) = 0 \\ & {\psi _1}{{\ddot q}_1} + {\psi _2}{q_1} - {\psi _3}q_1^2 - {\psi _4}{q_1}{{\dot q}_1} + {\psi _5}{{\dot q}_1} = \\ &\qquad - {\phi _1}\left[ {{k_1}\left( {{\phi _1}{q_1} - {z_1}} \right) + {c_1}\left( {{\phi _1}{{\dot q}_1} - {{\dot z}_1}} \right) + \left( {{m_1} + {m_2}} \right)g} \right] \end{split}\right\} $$ (19) 通过观察式(19)可以发现, 由于系数中显含时间t, 且由于$ {\psi _3} $和$ {\psi _4} $的存在, 该动力学方程组为时变和非线性的, 很难得到解析解, 通常采用数值方法进行求解. 式(19)可整理成矩阵的形式, 如下
$$ {\boldsymbol{M}} \cdot \ddot {\boldsymbol{X}} + {\boldsymbol{C}} \cdot \dot {\boldsymbol{X}} + {\boldsymbol{K}} \cdot {\boldsymbol{X}} = {\boldsymbol{F}} $$ (20) 其中
$$\left. \begin{split} & {\boldsymbol{M}} = \left[ {\begin{array}{*{20}{c}} 0&{{m_2}}&0 \\ {{m_1}}&0&0 \\ 0&0&{{\psi _1}} \end{array}} \right] \\ & {\boldsymbol{C}} = \left[ {\begin{array}{*{20}{c}} { - {c_2}}&{{c_2}}&0 \\ {{c_1} + {c_2}}&{ - {c_2}}&{ - {c_1}{\phi _1}} \\ { - {c_1}{\phi _1}}&0&{{\psi _5} + {c_1}{\phi _1}^2} \end{array}} \right] \\ & {\boldsymbol{K}} = \left[ {\begin{array}{*{20}{c}} { - {k_2}}&{{k_2}}&0 \\ {{k_1} + {k_2}}&{ - {k_2}}&{ - {k_1}{\phi _1}} \\ { - {k_1}{\phi _1}}&0&{{\psi _2} + {k_1}{\phi _1}^2} \end{array}} \right] \\ & {\boldsymbol{F}} = \left[ {\begin{array}{*{20}{c}} 0 \\ 0 \\ {{\psi _3}q_1^2 + {\psi _4}{q_1}{{\dot q}_1} - {\phi _1}\left( {{m_1} + {m_2}} \right)g} \end{array}} \right] \end{split} \right\}$$ (21) 系数矩阵$ {\boldsymbol{M}},\;{\boldsymbol{C}},\;{\boldsymbol{K}},\;{\boldsymbol{F}} $分别为质量矩阵、阻尼矩阵、刚度矩阵和载荷列阵. 为便于编程计算, 式(19)中的非线性项写入式(20)载荷列阵$ {\boldsymbol{F}} $中, 在每一数值积分步, 可采用上一积分步的计算值近似求解非线性项. 式(20)即可采用常规数值积分方法进行求解, 如经典的Newmark-β法.
2. 参数冻结精细指数积分方法的应用
将方程组(19)降阶为一阶动力学方程组, 可表示为如下形式
$$ {\boldsymbol{\dot u}} = {\boldsymbol{L}}\left( t \right){\boldsymbol{u}} + {\boldsymbol{N}}\left( {{\boldsymbol{u}},t} \right) $$ (22) 其中, $ {\boldsymbol{u}} = {\left[ {{z_1},{z_2},{q_1},{{\dot z}_1},{{\dot z}_2},{{\dot q}_1}} \right]^{\text{T}}} $为状态向量, $ {\boldsymbol{L}}\left( t \right) $为时变系数矩阵, $ {\boldsymbol{N}}\left( {{\boldsymbol{u}},t} \right) $为方程组的非线性部分, 分别为
$$\left.\begin{aligned} & {\boldsymbol{L}} = \left[ {\begin{array}{*{20}{c}} { {{{\boldsymbol{l}}_{11}}} }&{ {{{\boldsymbol{l}}_{12}}}} \\ { {{{\boldsymbol{l}}_{21}}} }&{ {{{\boldsymbol{l}}_{22}}} } \end{array}} \right] \\ & {{\boldsymbol{l}}_{11}} = \left[ {\begin{array}{*{20}{c}} 0&0&0 \\ 0&0&0 \\ 0&0&0 \end{array}} \right],\quad {{\boldsymbol{l}}_{12}} = \left[ {\begin{array}{*{20}{c}} 1&0&0 \\ 0&1&0 \\ 0&0&1 \end{array}} \right] \\ & {{\boldsymbol{l}}_{21}} = \left[ {\begin{array}{*{20}{c}} { - \dfrac{{{k_1} + {k_2}}}{{{m_1}}}}&{\dfrac{{{k_2}}}{{{m_1}}}}&{\dfrac{{{k_1}}}{{{m_1}}}{\phi _1}} \\ {\dfrac{{{k_2}}}{{{m_2}}}}&{ - \dfrac{{{k_2}}}{{{m_2}}}}&0 \\ {\dfrac{{{k_1}}}{{{\psi _1}}}{\phi _1}}&0&{ - \dfrac{{{\psi _2}}}{{{\psi _1}}} - \dfrac{{{k_1}}}{{{\psi _1}}}\phi _1^2} \end{array}} \right] \\ & {{\boldsymbol{l}}_{22}} = \left[ {\begin{array}{*{20}{c}} { - \dfrac{{{c_1} + {c_2}}}{{{m_1}}}}&{\dfrac{{{c_2}}}{{{m_1}}}}&{\dfrac{{{c_1}}}{{{m_1}}}{\phi _1}} \\ {\dfrac{{{c_2}}}{{{m_2}}}}&{ - \dfrac{{{c_2}}}{{{m_2}}}}&0 \\ {\dfrac{{{c_1}}}{{{\psi _1}}}{\phi _1}}&0&{ - \dfrac{{{\psi _5}}}{{{\psi _1}}} - \dfrac{{{c_1}}}{{{\psi _1}}}\phi _1^2} \end{array}} \right] \end{aligned}\right\} $$ (23) $$ {\boldsymbol{N}}\left( {{\boldsymbol{u}},t} \right) = \left\{ {\begin{array}{*{20}{c}} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ {\dfrac{{{\psi _3}}}{{{\psi _1}}}q_1^2 + \dfrac{{{\psi _4}}}{{{\psi _1}}}{q_1}{{\dot q}_1} - \dfrac{{\left( {{m_1} + {m_2}} \right)g}}{{{\psi _1}}}{\phi _1}} \end{array}} \right\} $$ (24) 对于式(22)中线性部分$ {\boldsymbol{\dot u}} = {\boldsymbol{L}}\left( t \right){\boldsymbol{u}} $的处理, 可以采用精细积分法的思想[28], 其能高效地处理齐次线性定常微分方程. 将时间以等步长$ \tau $离散, 即$ {t_k} = k\tau $, $ {{\boldsymbol{u}}_k} = {\boldsymbol{u}}\left( {{t_k}} \right) $, $ k = 0,1,2, \cdot \cdot \cdot $, 如当$ {\boldsymbol{L}} $为常系数矩阵时, $ {\boldsymbol{L}} = {{\boldsymbol{L}}_C} $, 则有
$$ {{\boldsymbol{u}}_{k + 1}} = {\boldsymbol{A}}{{\boldsymbol{u}}_k},\;{\boldsymbol{A}} = \exp \left( {{{\boldsymbol{L}}_C}\tau } \right) $$ (25) 如能精确求解矩阵$ {\boldsymbol{A}} $, 即可由式(25)计算出每一离散时间步的数值解$ {{\boldsymbol{u}}_k} $. 对于矩阵$ {\boldsymbol{A}} $的求解, 可采用精细积分法[28], 其思路是利用指数函数的加法原理, 将时间步长分为$ {2^n} $等分, 如令$ \Delta t = \tau /{2^n} $, 则有[41]
$$ {\boldsymbol{A}} = \exp \left( {{{\boldsymbol{L}}_C}\tau } \right) = {\left[ {\exp \left( {{{\boldsymbol{L}}_C}\Delta t} \right)} \right]^{{2^n}}} $$ (26) 由于$ \Delta t = \tau /{2^n} $是一个极小的时间段, 如对$ \exp \left( {{{\boldsymbol{L}}_C}\Delta t} \right) $进行泰勒展开, 取其前5项就能达到足够的精度, 即
$$\left. \begin{aligned} & \exp \left( {{{\boldsymbol{L}}_C}\Delta t} \right) \approx {\boldsymbol{I}} + {{\boldsymbol{A}}_1} \\ & {{\boldsymbol{A}}_1} = {{\boldsymbol{L}}_C}\Delta t + {\left( {{{\boldsymbol{L}}_C}\Delta t} \right)^2}/2!+ \\ &\qquad {\left( {{{\boldsymbol{L}}_C}\Delta t} \right)^3}/3! + {\left( {{{\boldsymbol{L}}_C}\Delta t} \right)^4}/4! \end{aligned} \right\}$$ (27) 其中, $ {\boldsymbol{I}} $为单位矩阵, $ {{\boldsymbol{A}}_1} $是一个很小的矩阵. 为避免计算过程中的舍入误差, 由式(26)和式(27)有
$$ {\boldsymbol{A}} = {\left( {{\boldsymbol{I}} + {{\boldsymbol{A}}_1}} \right)^{{2^n}}} = {\left( {{\boldsymbol{I}} + 2{{\boldsymbol{A}}_1} + {\boldsymbol{A}}_1^2} \right)^{{2^{n - 1}}}} $$ (28) 如定义矩阵序列
$$ {{\boldsymbol{A}}_{i + 1}} = 2{{\boldsymbol{A}}_i} + {\boldsymbol{A}}_i^2,\;i = 1,2,3, \cdot \cdot \cdot $$ (29) 则由式(28)可得递推关系如下
$$\begin{split} & {\boldsymbol{A}} = {\left( {{\boldsymbol{I}} + {{\boldsymbol{A}}_1}} \right)^{{2^n}}} = {\left( {{\boldsymbol{I}} + {{\boldsymbol{A}}_2}} \right)^{{2^{n - 1}}}}= \\ &\qquad {\left( {{\boldsymbol{I}} + {{\boldsymbol{A}}_3}} \right)^{{2^{n - 2}}}} = \cdot \cdot \cdot = {\left( {{\boldsymbol{I}} + {{\boldsymbol{A}}_{n + 1}}} \right)^{{2^0}}} \end{split} $$ (30) 由式(29)求得矩阵序列$ {{\boldsymbol{A}}_{n + 1}} $的值, 代入式(30), 与$ {\boldsymbol{I}} $阵相加即可求得矩阵$ {\boldsymbol{A}} $的值. 此过程中, 利用了指数函数的性质, 巧妙地运用加法代替乘法, 避免了严重的舍入误差, 使算法具有很高的精度, 通常取$ n = 20 $即可保证良好的数值性能. 且由于在计算中的每一时间步$ \tau $内插入了$ {2^n} $个点, 即使步长$ \tau $值取的较大, 也不影响计算的精确性[44].
精细积分法在计算一般的常系数齐次线性微分方程时, 能够达到极高的精度. 但车桥耦合问题的动力学方程组系数矩阵具有时变性, 即式(22)中$ {\boldsymbol{L}}\left( t \right) $为时变系数矩阵. 此时可以采用时间参数冻结的思想, 基于时间参数冻结技术的基本概念, 在每一个离散步长计算间隔内, 假设$ {\boldsymbol{L}}\left( t \right) $为常数矩阵, 则系统定常, 即可应用精细积分法求得每一时间区间内的响应. 如式(25), 在计算$ {{\boldsymbol{u}}_{k + 1}} $的数值结果时, 依据平均的思想, 可假设
$$ {{\boldsymbol{L}}_C} = {{\boldsymbol{L}}_{k + 1}} = \frac{{{\boldsymbol{L}}\left( {{t_{k + 1}}} \right) + {\boldsymbol{L}}\left( {{t_k}} \right)}}{2} $$ (31) 式(31)中由于在每一离散步内引入了近似, 因此, 当时间间隔$ \tau $越小时, 近似解越接近于精确解.
进一步考虑非线性部分的处理, 式(22)中对时间离散后, 解的形式可表示为
$$\left. \begin{split} & {{\boldsymbol{u}}_{k + 1}} = {\boldsymbol{A}}{{\boldsymbol{u}}_k} + \int_{{t_k}}^{{t_{k + 1}}} {\exp \left[ {{{\boldsymbol{L}}_C}\left( {{t_{k + 1}} - s} \right)} \right]} {\boldsymbol{N}}\left( {{\boldsymbol{u}}\left( s \right),s} \right)\mathrm{d}s \\ & {\boldsymbol{A}} = \exp \left( {{{\boldsymbol{L}}_C}\tau } \right) \end{split}\right\}$$ (32) 其中, 右端第1项为线性部分, 即采用上述参数冻结后的精细积分法进行求解, 右端第2部分积分项为杜哈梅积分. 对式(32)采用不同的方法求解, 可得到不同的指数积分法, 应用较多的是指数Runge-Kutta法, 本文即采用指数Runge-Kutta法进行计算. s级的指数Runge-Kutta法可表示为[41]
$$\left. \begin{aligned} & {{\boldsymbol{u}}_{k + 1}} = \exp \left( {{{\boldsymbol{L}}_C}\tau } \right){{\boldsymbol{u}}_k} + \tau \sum\limits_{i = 1}^s {{{\boldsymbol{b}}_i}\left( {{{\boldsymbol{L}}_C}\tau } \right)} {{\boldsymbol{K}}_i} \\ & {{\boldsymbol{K}}_i} = {\boldsymbol{N}}\left[ {\exp \left( {{c_i}{{\boldsymbol{L}}_C}\tau } \right){{\boldsymbol{u}}_k} + \tau \sum\limits_{j = 1}^s {{{\boldsymbol{a}}_{ij}}\left( {{{\boldsymbol{L}}_C}\tau } \right){{\boldsymbol{K}}_i},{t_k} + {c_i}\tau } } \right] \end{aligned} \right\}$$ (33) 其中, $ {{\boldsymbol{a}}_{ij}}\left( {{{\boldsymbol{L}}_C}\tau } \right) $和$ {{\boldsymbol{b}}_i}\left( {{{\boldsymbol{L}}_C}\tau } \right) $为待定矩阵, $ {c_i} $为待定常数. $ {{\boldsymbol{a}}_{ij}}\left( {{{\boldsymbol{L}}_C}\tau } \right) $和$ {{\boldsymbol{b}}_i}\left( {{{\boldsymbol{L}}_C}\tau } \right) $可以根据积分因子法按照如下的规律给出[41]
$$\left. \begin{aligned} & {{\boldsymbol{a}}_{ij}}\left( {{{\boldsymbol{L}}_C}\tau } \right) = {a_{ij}}\exp \left[ {\left( {{c_i} - {c_j}} \right){{\boldsymbol{L}}_C}\tau } \right] \\ & {{\boldsymbol{b}}_i}\left( {{{\boldsymbol{L}}_C}\tau } \right) = {b_i}\exp \left[ {\left( {1 - {c_i}} \right){{\boldsymbol{L}}_C}\tau } \right] \end{aligned} \right\}$$ (34) 其中, 系数$ {a_{ij}},{b_i},{c_i} $与经典Runge-Kutta法中系数相同, 即每一种经典Runge-Kutta法都对应有一种指数积分法. 式(33)中关于指数矩阵的求解均可由精细积分法进行精确计算, 这种将精细积分法和指数积分法结合得到的方法称为精细指数积分法[41]. 本文基于精细指数积分法与时间参数冻结技术, 创新性地将二者结合到一起, 提出参数冻结精细指数积分法, 并将其应用于车桥耦合动力学问题的求解中. 参数冻结精细指数积分法的计算重点在于每一迭代步中对时变矩阵$ {\boldsymbol{L}}\left( t \right) $的近似处理和对指数矩阵$ {\boldsymbol{A}} $的精细积分求解, 其数值计算格式容易构造, 且可以方便地处理时变、非线性问题.
3. 算例分析
如前文所述, 每一种Runge-Kutta法都可构造相对应的指数积分法, 本文基于简单的四级四阶显式Runge-Kutta法(RK4-4)[41]构造了相应的参数冻结精细指数积分格式(FPEI4-4).
首先, 为验证本文车桥耦合模型及理论计算的有效性, 基于车桥耦合实验平台, 进行了实验测试. 用于测试的实验装置参数为: 跨度$ 3{\text{ m}} $的两端简支实验桥, 采用空心铝板制作, 其密度为$ 2700{\text{ kg}}/{{\text{m}}^3} $, 横截面积为$ 9.24 \times {10^{ - 4}}{{\text{ m}}^2} $, 弹性模量为$ 70{\text{ GPa}} $, 惯性矩为$ 2.67 \times {10^{ - 7}}{{\text{ m}}^4} $; 实验小车质量为$ 2{\text{ kg}} $, 轮间距为$ 9.5 \times {10^{ - 2}}{\text{ m}} $, 车速为$ 0.2{\text{ m/s}} $. 实验时在桥跨中布置位移传感器, 利用ARTDAQ信号采集显示软件进行数据采集, 采样频率为1000 Hz. 实验测试模型及测试现场如图2和图3所示. 图2中, 1为西门子V90伺服电机, 2为传动齿轮, 3为传动链条, 4为挡板, 5为小车, 6为实验桥支座, 7为桥板, A和C为辅助梁, B为实验梁.
实验测试结果(test result, TR)与基于1/4车−桥耦合模型的FPEI4-4格式数值计算结果进行了对比, FPEI4-4格式选取了时间步长为$ \tau=1.0 \times 10^{-3}\text{ s} $, 结果如图4所示.
从图4中可以看出, 实验测试和数值计算结果的演化趋势基本吻合, 挠度峰值误差很小. 图4中的结果对比能够说明本文理论建模和数值计算的有效性, 其具有实际应用价值. 此外, 图4中数值计算与实验测试结果误差主要来源于车辆简化模型的误差, 以及实验装置的制作误差、外界环境因素等. 如采用更加精细化的车辆模型, 可减小理论计算和实验结果的误差, 但即使是选用简单的1/4车辆模型进行理论分析, 其依然能够较为精确地预测挠度峰值结果.
理论结果计算的精确性主要取决于所建立的车桥耦合动力学模型的精细化程度, 本文所提算法主要关注于对已建立的动力学微分方程的高效求解. 后文应用所构造的FPEI4-4格式与经典的Newmark-β算法进行了数值计算, 对比了2种算法在求解车桥耦合动力学方程组(19)时的数值特性.
数值模拟中采用了重型载货汽车模型, 具体车辆参数见表1[45], 表2给出了桥的参数, 包括桥面沥青铺装层参数和桥体参数.
表 1 汽车参数Table 1. Vehicle parametersSymbol Value m2/kg 18000 m1/kg 2000 k2/(N·m−1) 9000000 c2/(N·s·m−1) 80000 k1/(N·m−1) 36000000 c1/(N·s·m−1) 700 表 2 桥面沥青铺装层和桥体参数Table 2. Pavement and bridge parametersSymbol Value L/m 40 b/m 8 h1/m 0.1 h2/m 2 E1/GPa 8 E2/GPa 34.5 $ {\rho _1} $/(kg·m−3) 2000 $ {\rho _2} $/(kg·m−3) 2500 $ {\eta _1} $ 0.1 $ {\eta _2} $ 0.01 由于文中车桥耦合动力学方程组具有时变和非线性的特性, 难以得到解析解, 因此选用了小步长下高精度的二级4阶隐式辛Runge-Kutta算法(SRK2-4)的计算结果作为参考.
s级隐式Runge-Kutta法的一般格式为[46]
$$\left.\begin{aligned} & {{\boldsymbol{u}}_{k + 1}} = {{\boldsymbol{u}}_k} + \tau \sum\limits_{i = 1}^s {{c_i}} {{\boldsymbol{K}}_i} \\ & {{\boldsymbol{K}}_i} = f\left( {{{\boldsymbol{u}}_k} + \tau \sum\limits_{j = 1}^s {{b_{ij}}{{\boldsymbol{K}}_j},{t_n} + {a_i}\tau } } \right), {r = 1,2, \cdots ,s} \end{aligned}\right\} $$ (35) 其中, $ i,j = 1,2, \cdot \cdot \cdot ,s $; $ {c_j} \geqslant 0 $; $ \tau $为数值积分步长
$$ \sum\limits_{i = 1}^s {{c_i}} = 1,\quad \sum\limits_{j = 1}^s {{{{b}}_{ij}} = {a_i}} ,\quad \sum\limits_{j = 1}^s {{a_j}} = 1 $$ (36) 式(35)中当其系数满足以下条件时为辛Runge-Kutta法
$$ {c_i}{c_j} - {b_{ij}}c{}_i - {b_{ji}}{c_j} = 0,\;\; {i,j = 1,2, \cdots, s}$$ (37) 当系数$ {b_{ij}}和\;c{}_i $取不同值时, 可得到不同的辛Runge-Kutta格式. 比较常用的SRK2-4积分格式, 其系数表示如下
$$\left. \begin{split} & \left[ {\begin{array}{*{20}{c}} {{b_{11}}}&{{b_{12}}} \\ {{b_{21}}}&{{b_{22}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\dfrac{1}{4}}&{\dfrac{1}{4} - \dfrac{{\sqrt 3 }}{6}} \\ {\dfrac{1}{4} + \dfrac{{\sqrt 3 }}{6}}&{\dfrac{1}{4}} \end{array}} \right] \\ & \left[ {\begin{array}{*{20}{c}} {{c_1}}&{{c_2}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\dfrac{1}{2}}&{\dfrac{1}{2}} \end{array}} \right] \end{split} \right\}$$ (38) 同时, 为了能够说明针对动力学微分方程模型数值计算结果的有效性, 构造了本文模型中桥梁位移的近似解析解. 基于文献[47]中的思路, 式(19)中第3个方程可简化为
$$ \begin{split} & {\psi _1}{{\ddot q}_1} + {\psi _2}{q_1} - {\psi _3}q_1^2 - {\psi _4}{q_1}{{\dot q}_1} + {\psi _5}{{\dot q}_1} = \\ &\qquad - {\phi _1}\left( {{m_1} + {m_2}} \right)g \end{split} $$ (39) 由于本文工况中广义位移$ {q_1} $值和广义速度$ {\dot q_1} $值较小, 非线性项$ - {\psi _3}q_1^2 - {\psi _4}{q_1}{\dot q_1} $对振动的影响较弱, 若略去两项非线性项, 并引入参数
$$\left.\begin{split} & P = - \left( {{m_1} + {m_2}} \right)g,\;\omega = \frac{{\text{π} v}}{L},\;{\omega _N} = \sqrt {\frac{{{\psi _2}}}{{{\psi _1}}}} ,\;s = \frac{\omega }{{{\omega _N}}} \\ & \zeta = \frac{{{\psi _5}}}{{2\sqrt {{\psi _1}{\psi _2}} }},\;{\omega _D} = {\omega _N}\sqrt {1 - {\zeta ^2}} ,\;{U_0} = \frac{P}{{{\psi _2}}}\end{split} \right\}$$ (40) 则可得到桥梁广义位移$ {q_1} $的近似解析解(approximate analytical solution, AAS)为
$$\left.\begin{aligned} & {q_1}\left( t \right) = {{\mathrm{e}}^{ - \zeta {\omega _N}t}}\left( {{C_1}\cos {\omega _D}t + {C_2}\sin {\omega _D}t} \right)+ \\[-1pt] & \qquad {C_3}\cos \omega t + {C_4}\sin \omega t \\[-1pt] & {C_1} = {q_1}\left( 0 \right) + \frac{{2{U_0}\zeta s}}{{{{\left( {1 - {s^2}} \right)}^2} + {{\left( {2\zeta s} \right)}^2}}} \\[-1pt] & {C_2} = \frac{{{{\dot q}_1}\left( 0 \right) + \zeta {\omega _N}{q_1}\left( 0 \right)}}{{{\omega _D}}} - \frac{{{U_0}s\left( {1 - {s^2} - 2{\zeta ^2}} \right)}}{{\sqrt {1 - {\zeta ^2}} \left[ {{{\left( {1 - {s^2}} \right)}^2} + {{\left( {2\zeta s} \right)}^2}} \right]}} \\[-1pt] & {C_3} = \frac{{ - 2{U_0}\zeta s}}{{{{\left( {1 - {s^2}} \right)}^2} + {{\left( {2\zeta s} \right)}^2}}} \\[-1pt] & {C_4} = \frac{{{U_0}\left( {1 - {s^2}} \right)}}{{{{\left( {1 - {s^2}} \right)}^2} + {{\left( {2\zeta s} \right)}^2}}} \end{aligned} \right\}$$ (41) 图5给出了车速分别为10 和50 m/s时, 近似解析解、FPEI4-4格式和SRK2-4格式计算的桥梁跨中挠度的结果, 其中计算时FPEI4-4格式选取了大时间步长$ \tau = 1.0 \times {10^{ - 2}}{\text{ s}} $、SRK2-4格式选取了小时间步长$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $.
![]()
图 5 近似解析解、时间步长为$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $的FPEI4-4格式和时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $的SRK2-4格式计算跨中挠度结果Figure 5. The results of mid-span deflection calculated by analytical approximate solution, FPEI4-4 scheme with $ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $ and SRK2-4 scheme with $ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $从图5中可以看出, 每种工况下3条曲线的演化趋势相同, 车速较小时(v = 10 m/s)曲线吻合得更好, 这充分体现了数值计算结果的有效性. 此外, 从局部放大图中可以看出, FPEI4-4格式和SRK2-4格式的计算结果始终几乎重合, 而与近似解析解有一定的偏差. 这是因为近似解析解忽略了耦合因素和非线性因素, 而数值格式是直接基于原微分方程(19)计算的原因.
从图5中的对比也可以体现出, 当车速较大时, 车桥耦合和非线性因素对桥梁的振动特性影响更为明显. 由于所构造的近似解析解忽略了部分因素, 不能严格地反应振动响应的演化趋势, 因此后文为更好地体现所提出方法的数值特性, 选用了SRK2-4格式的数值计算结果进行参照对比.
假设车速分别为10 和20 m/s, 利用SRK2-4格式计算了文中车桥耦合动力学方程组. 图6中给出了当选取时间步长分别为$ \tau =1.0 \times {10^{ - 4}} $和$ 1.0 \times {10^{ - 5}}$时, SRK2-4格式计算的桥梁跨中挠度结果之差.
![]()
图 6 时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $和$ \tau =1.0 \times {10^{ - 5}}{\text{ s}} $时SRK2-4格式计算跨中挠度结果之差Figure 6. The difference between the results of mid-span deflection calculated by SRK2-4 scheme under time step $ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $and $ \tau =1.0 \times {10^{ - 5}}{\text{ s}} $在数值积分运算中, 当选取数值积分步长越小时, 数值计算结果越收敛于解析解. 从图6中可以看出, 选取较小时间步长为$ \tau =1.0 \times {10^{ - 4}}$和$1.0 \times {10^{ - 5}} $时, 两种车速工况下的数值结果差值已经非常小, 达到了10−17 m量级, 此时已非常接近解析解. 因此, 可将选取小步长下的SRK2-4格式计算结果作为参照解. 由于SRK2-4格式为隐式格式, 为提高计算效率, 本文选用了积分时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $时的SRK2-4格式进行数值计算. 后文中的误差均表示不同时间步长下相应算法的计算结果与时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $时SRK2-4格式计算结果的差值.
3.1 不同数值积分时间步长时的计算结果对比
对比了选取不同时间步长时, 应用FPEI4-4格式和Newmark-β算法计算文中车桥耦合模型的精度. 选取车速为相对较快的20 m/s过桥, 以时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $时的SRK2-4格式计算结果作为参照解, 应用本文构造的FPEI4-4格式和经典的Newmark-β算法得到了桥梁跨中挠度计算结果与参照解的误差, 如图7所示.
![]()
图 7 FPEI4-4格式和Newmark-β算法求解跨中挠度计算误差Figure 7. The numerical errors of mid-span deflection calculated by FPEI4-4 scheme and Newmark-β algorithm从图7可以看出, 在选取时间步长为$ \tau =1.0 \;\times {10^{ - 3}}{\text{ s}} $时, FPEI4-4格式与参照解的误差在10−11 m量级, 在选取时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $时, 与参照解的误差达到了10−13 m量级. 作为一种显式计算格式, 相比于相同数值积分步长时的经典Newmark-β算法, 本文所构造的FPEI4-4格式具有明显更小的误差、更高的计算精度.
当进一步增大时间步长, 取时间步长为$ \tau = 5.0 \times {10^{ - 3}}{\text{ s}} $时, 如图8所示, Newmark-β算法计算结果已出现明显的数值发散现象, 导致计算失败. 如继续增大时间步长, Newmark-β算法计算结果会直接发散. 而从图8中还可以看出, 即使将时间步长增加到很大的$ \tau =1.0 \times {10^{ - 1}}{\text{ s}} $时, FPEI4-4格式同样能给出较好的计算结果, 这体现了该计算格式优异的数值特性. 由于精细积分法自身的特性, FPEI4-4格式即使选取大时间步长, 也不会出现数值发散现象, 而大步长下的数值误差主要来源于对系数矩阵$ {\boldsymbol{L}}\left( t \right) $参数冻结时引入的误差.
![]()
图 8 FPEI4-4格式与Newmark-β算法得到的桥梁跨中挠度结果Figure 8. The results of mid-span deflection calculated by FPEI4-4 scheme and Newmark-β algorithm表3进一步给出了2种算法在选取不同时间步长时的计算时间和最大误差. 从表3中可以看出, Newmark-β算法作为一种经典的动力学方程数值求解方法, 在相同的小时间步长下, 其比FPEI4-4格式的计算速度更快. 而FPEI4-4格式则具有明显的精度优势. 此外可以看出FPEI4-4格式具有更好的数值稳定性, 能够选取较大的积分时间步长来加快计算速度, 且在大步长下仍能保持较好的计算精度, 而在大步长情况下Newmark-β算法将直接发散使得计算失败. 如选取$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $的FPEI4-4格式与$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $的Newmark-β算法对比, FPEI4-4格式在计算时间(0.031098 s)远小于Newmark-β算法(0.869262 s)的情况下, 其最大误差仅约为Newmark-β算法的1/50, 具有更好的效率和精度优势.
表 3 Newmark-β算法和FPEI4-4格式计算跨中挠度对比Table 3. Comparison of mid-span deflection calculated by Newmark-β algorithm and FPEI4-4 schemeTime step/s Computation time/s Maximum error/m Newmark-β algorithm $ {10^{ - 4}} $ 0.869262 $ - 2.222\;306 \times {10^{ - 7}} $ $ {10^{ - 3}} $ 0.078846 $ 2.269\;753 \times {10^{ - 6}} $ $ {10^{ - 2}} $ — — FPEI4-4 scheme $ {10^{ - 4}} $ 2.849518 $ 4.373\;681 \times {10^{ - 13}} $ $ {10^{ - 3}} $ 0.319489 $ 4.373\;730 \times {10^{ - 11}} $ $ {10^{ - 2}} $ 0.031098 $ 4.409\;290 \times {10^{ - 9}} $ $ {10^{ - 1}} $ 0.003465 $ - 1.741\;726 \times {10^{ - 6}} $ 后文中为了说明本文FPEI4-4格式的优势, 直接选取了大时间步长进行计算, 即取$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $, 在该时间步长下传统的Newmark-β算法将直接发散.
3.2 数值积分的长时间稳定性
为了更好地体现本文所构造格式的长时间数值特性, 进行了长时间的数值积分运算结果对比.
首先关注长时间振动解的稳定性质. 因车在桥上时间较短, 当车下桥以后, 由式(19)描述的车桥耦合作用将消失, 其后桥将做有初位移和初速度的自由振动, 其初位移和初速度由车下桥时桥振动位移和速度给出, 其值均不大, 靠近相图(0, 0)点.
关注于桥自由振动阶段, 此时动力学方程为
$$ {\psi _1}{\ddot q_1} + {\psi _2}{q_1} - {\psi _3}q_1^2 - {\psi _4}{q_1}{\dot q_1} + {\psi _5}{\dot q_1} = 0 $$ (42) 如令$ x = {q_1},\;y = {\dot q_1} $, 可得一阶微分方程组为
$$ \left. \begin{split} & \frac{{{\text{d}}x}}{{{\text{d}}t}} = y \\ & \frac{{{\text{d}}y}}{{{\text{d}}t}} = - \frac{{{\psi _5}}}{{{\psi _1}}}y + \frac{{{\psi _4}}}{{{\psi _1}}}xy + \frac{{{\psi _3}}}{{{\psi _1}}}{x^2} - \frac{{{\psi _2}}}{{{\psi _1}}}x \end{split} \right\} $$ (43) 原点(0, 0)为满足式(43)的奇点, 是式(43)的零解.
如考虑桥梁黏性阻尼, 式(43)的线性近似方程组特征方程为
$$ \left| {\begin{array}{*{20}{c}} {0 - \lambda }&1 \\ { - \dfrac{{{\psi _2}}}{{{\psi _1}}}}&{ - \dfrac{{{\psi _5}}}{{{\psi _1}}} - \lambda } \end{array}} \right| = 0 $$ (44) 或
$$ {\psi _1}{\lambda ^2} + {\psi _5}\lambda + {\psi _2} = 0 $$ (45) 可得
$$ {\lambda _{1,2}} = \frac{{ - {\psi _5} \pm \sqrt {\psi _5^2 - 4{\psi _1}{\psi _2}} }}{{2{\psi _1}}} $$ (46) 代入数值可得特征方程两个根均具有负实部, 可知零解$ x = y = 0 $为渐近稳定的, 此时奇点(0, 0)为稳定焦点, 相图中的轨线表现为盘旋趋近原点(0, 0).
如忽略桥梁黏性阻尼, 则线性近似方程组特征方程为
$$ \left| {\begin{array}{*{20}{c}} {0 - \lambda }&1 \\ { - \dfrac{{{\psi _2}}}{{{\psi _1}}}}&{0 - \lambda } \end{array}} \right| = 0 $$ (47) 或
$$ {\psi _1}{\lambda ^2} + {\psi _2} = 0 $$ (48) 可得
$$ {\lambda _{1,2}} = \pm \frac{{\sqrt {\psi _5^2 - 4{\psi _1}{\psi _2}} }}{{2{\psi _1}}} $$ (49) 代入数值可得特征方程的2个根均为纯虚根, 此时奇点(0, 0)为中心型奇点, 零解为稳定但非渐近稳定的, 其相图中的轨线表现为以(0, 0)为中心的圆.
综上所述, 由于初始条件趋近(0, 0)点, 且(0, 0)点为稳定奇点, 因此长时间积分计算的解是稳定的. 这里忽略了桥梁黏弹性特性, 则汽车在过桥以后, 桥梁将做幅值相等的无耗散往复自由振动.
首先给出了Newmark-β算法在选取数值积分步长$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $时计算的跨中挠度时程曲线, 以及跨中挠度−速度相图, 如图9所示. 算例中选取车速为20 m/s, 汽车下桥后继续计算了100 s的桥梁自由振动, 计算耗时为35.948393 s.
从图9中可以看出, 即使选择了较小的时间步长($ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $), Newmark-β算法在长时间积分时, 仍具有明显的数值耗散性. 图9(a)中显示, Newmark-β算法计算的汽车下桥后桥梁自由振动数值解具有发散的特性, 长时间数值积分后会导致计算结果的失效. 从图9(b)中也可以看出, 在自由振动的周期阶段, 跨中挠度−速度相图表现出了明显的发散现象.
选取相同的工况, 图10进一步给出了显式FPEI4-4格式计算的桥梁跨中挠度以及挠度−速度相图的数值结果. 数值积分计算时, 选用了大时间步长$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $.
从图10中可以看出, 即使选择了较大的时间步长, 本文所提出的显式FPEI4-4格式仍具有良好的数值特性, 数值计算结果保持了长时间的数值稳定性, 无发散迹象. 由于此时零解为稳定但非渐近稳定的, 其相图中的轨线表现为以(0, 0)为中心的圆, 如图10(b)中右侧所示, 图10(b)左侧为车桥耦合初始阶段的轨线. 由于精细积分法自身的特点, 可以选取大时间步长进行数值积分运算, 本算例中步长取为$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $, 计算耗时仅为1.370017 s. 文中虽然对时变矩阵$ {\boldsymbol{L}}\left( t \right) $在每一离散步进行了参数冻结近似, 但结果显示在选取较大步长时依然具有良好的数值精度. 通过与时间步长$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $时的隐式SRK2-4格式计算结果对比, 自由振动阶段的最大误差稳定维持在$ 1.38 \times {10^{ - 6}}{\text{ m}} $量级. 本文选取的计算模型矩阵维数较小, 如当处理高维时变车桥耦合动力学问题时, 选取大步长所带来的计算效率优势将更加明显.
3.3 不同时变程度的计算结果对比
由于汽车为移动载荷, 导致车桥耦合动力学方程组具有时变性, 由式(22)可以看出, 当汽车和桥梁参数确定的情况下, 系数矩阵时变的快慢和车速v直接相关. 本文数值格式构造时, 在每一离散步长时间间隔内对时变矩阵$ {\boldsymbol{L}}\left( t \right) $进行了参数冻结处理, 将时变系数矩阵近似为每一离散步内的常数矩阵. 这种近似带来的误差和矩阵时变的快慢以及步长的选取有关, 因此进一步分析了系数矩阵具有不同时变程度时(不同车速)所构造的FPEI4-4格式的数值特性.
图11给出了选取大时间步长为$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $时, 不同车速下FPEI4-4格式求解的桥梁跨中挠度与参照解的计算误差. 从图11中可以明显看出, 车速越快, FPEI4-4格式计算结果与参照解的计算结果差别越大. 当车速为10 m/s时, 最大差值为$ 7.754\;9 \times {10^{ - 10}} $ m; 车速为20 m/s时, 最大差值为$ 4.409\;3 \times {10^{ - 9}} $ m; 车速为50 m/s时, 最大差值为$ 1.432\;6 \times {10^{ - 8}} $ m. 由于车速与时变矩阵变化快慢直接相关, 车速越大, 参数冻结引起的误差就会越大, 这符合预期. 但即使在车速为50 m/s时(时速为180 km/h), 误差仍可保持在10−8 m量级, 而此车速下跨中挠度的最大值为$ 1.743\;2 \times {10^{ - 3}} $ m, 跨中计算结果仍是可接受的. 因此本文所构造的FPEI4-4格式同样适于计算快车速、快时变问题, 具有实际应用性.
![]()
图 11 $ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $, 不同车速时FPEI4-4格式求解跨中挠度计算误差Figure 11. The numerical errors of mid-span deflection calculated by FPEI4-4 scheme under different vehicle speeds with time step $ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $3.4 时变系数矩阵$ {\boldsymbol{L}}\left( t \right) $冻结方式的对比
本文应用的数值格式在处理时变问题时, 采用了参数冻结技术, 对每一积分步的时变系数矩阵进行了近似处理, 本节对比分析选取不同的近似方式时对数值计算结果的影响. 下式列举了每一离散步对时变矩阵$ {\boldsymbol{L}}\left( t \right) $进行参数冻结的3种方式, 本文基于平均的思想, 应用的近似矩阵$ {{\boldsymbol{L}}_{C1}} $. 在此, 对比分析分别选取$ {{\boldsymbol{L}}}_{C1},\;{{\boldsymbol{L}}}_{C2}和{{\boldsymbol{L}}}_{C3} $时的数值特性
$$\left.\begin{split} & {{\boldsymbol{L}}_{C1}} = \left( {{\boldsymbol{L}}\left( {{t_{k + 1}}} \right) + {\boldsymbol{L}}\left( {{t_k}} \right)} \right)/2 \\ & {{\boldsymbol{L}}_{C2}} = {\boldsymbol{L}}\left( {{t_{k + 1}}} \right) \\ & {{\boldsymbol{L}}_{C3}} = {\boldsymbol{L}}\left( {{t_k}} \right) \end{split}\right\} $$ (50) 图12给出了选取不同近似矩阵时, 应用显式FPEI4-4格式计算的桥梁跨中挠度和挠度−速度相图的数值结果, 计算时选取车速为20 m/s, 积分步长为$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $, 考虑了桥梁的黏弹性特性.
![]()
图 12 $ \tau = 1.0 \times{10^{ - 2}}{\text{ s}}$, FPEI4-4格式选取不同时变矩阵冻结方式的计算结果Figure 12. The numerical results calculated by FPEI4-4 scheme under different parameter freezing forms with time step $\tau = 1.0 \times {10^{ - 2}}{\text{ s}}$从图12中可以看出, 由于考虑了桥梁的黏弹性, 在汽车下桥以后, 跨中位移振动迅速衰减为0, 在相图中曲线逐渐趋于(0, 0)点, 这和实际情况相符. 由于此时零解为渐近稳定的, 奇点(0, 0)为稳定焦点, 如图12(b)中右侧所示, 相图中的轨线表现为盘旋趋近原点(0, 0), 图12(b)左侧为车桥耦合初始阶段的轨线. 同时从图中可以看出, 应用本文构造的显式FPEI4-4格式进行数值积分计算时, 选取3种不同参数冻结方式时的图线几乎吻合, 均具有良好的数值计算效果. 虽然在对时变矩阵$ {\boldsymbol{L}}\left( t \right) $进行参数冻结时引入了误差, 但数值结果显示, 即使选择了较大的步长, 依然可以得到较为满意的结果, 这体现了本文所提出算法在计算车桥耦合振动问题时的优势.
进一步分析了采用不同参数冻结形式时的计算结果和积分步长$ \tau =1.0 \times {10^{ - 4}}{\text{ s}}$时SRK2-4格式计算参照解的误差, 结果如图13所示. 从图13中可以看出, 采用$ {{\boldsymbol{L}}_{C2}} $和$ {{\boldsymbol{L}}_{C3}} $参数冻结形式时产生的误差峰值相近, 而采用$ {{\boldsymbol{L}}_{C1}} $的时变参数冻结方式产生的误差要明显小于$ {{\boldsymbol{L}}_{C2}} $和$ {{\boldsymbol{L}}_{C3}} $, 具有更好的数值计算效果. 表4进一步给出了FPEI4-4格式选取不同积分时间步长时, 不同参数冻结形式与参照解的最大误差. 从表4中同样可以看出, 选取$ {{\boldsymbol{L}}_{C1}} $的时变参数冻结方式具有更好的数值精度, 而选用$ {{\boldsymbol{L}}_{C2}} $和$ {{\boldsymbol{L}}_{C3}} $参数冻结形式时精度稍差且峰值误差接近, 而此车速下跨中挠度的最大值为$ 1.395\;9 \times {10^{ - 3}} $m, 不同冻结方式的跨中计算结果均是可接受的. 基于平均化的思想, 在处理不同的车桥耦合动力学问题时, 建议采用$ {{\boldsymbol{L}}_{C1}} $的时变参数冻结方式.
![]()
图 13 $ \tau = 1.0 \times {10^{ - 2}}{\text{ s}}$, FPEI4-4格式不同时变矩阵冻结方式求解跨中挠度的计算误差Figure 13. The numerical errors of mid-span deflection calculated by FPEI4-4 scheme under different parameter freezing forms with time step $\tau = 1.0 \times {10^{ - 2}}{\text{ s}}$表 4 不同冰冻近似方式的对比Table 4. Comparison of different parameter freezing formsTime step/s Form Maximum error/m
$1.0 \times {10^{ - 2} }$LC1 $ 4.409\;290 \times {10^{ - 9}} $ LC2 $ 8.127\;175 \times {10^{ - 8}} $ LC3 $ 8.418\;344 \times {10^{ - 8}} $
$1.0 \times {10^{ - 3} }$LC1 $ 4.373\;730 \times {10^{ - 11}} $ LC2 $ 8.131\;892 \times {10^{ - 9}} $ LC3 $ 8.139\;901 \times {10^{ - 9}} $ 4. 结 论
本文基于车桥耦合振动系统中时变、非线性动力学方程组的求解问题, 提出了一种新的数值计算方法, 即时间参数冻结的精细指数积分方法. 基于简单的四级四阶显式Runge-Kutta格式(RK4-4)构造了相应的时间参数冻结精细指数积分格式(FPEI4-4格式). 以考虑了桥面铺装层、桥梁结构几何非线性和黏弹性特性的车桥耦合问题为例, 推导了描述该问题的时变、非线性动力学方程组, 应用本文提出的FPEI4-4格式进行了求解. 首先, 通过与实验测试结果进行对比, 体现了本文理论建模和数值计算的有效性和实用价值. 进一步, 通过与近似解析解、二级4阶隐式辛Runge-Kutta格式(SRK2-4)计算的参照解以及经典的Newmark-β数值积分格式计算结果进行对比, 体现了本文所提算法的良好数值特性. 文中主要结论如下:
(1) 所提出的时间参数冻结精细指数积分方法可方便处理时变、非线性动力学问题, 即使选取较大的积分步长也能得到满意的结果, 有效地提高了计算效率;
(2) 本文针对车桥耦合动力学问题提出的参数冻结精细指数积分法计算速度快, 具有长时间的数值积分稳定性, 能够有效地保持动力学问题的本质特性;
(3) 针对快时变问题, 对时变系数矩阵参数冻结引起的误差会增大, 但本文格式数值积分结果的误差仍可保持在较小量级, 计算结果仍是可接受的;
(4) 在每一离散步, 对于时变系数矩阵的参数冻结方式, 建议选用上一积分时刻和本积分时刻矩阵的均值, 但即使选择上一时刻矩阵值或当前时刻矩阵值进行近似, 也能取得较满意的数值积分结果.
本文所提出的参数冻结精细指数积分法预期可为时变、非线性的车桥耦合动力学问题研究提供一种新的计算途径.
-
![]()
图 5 近似解析解、时间步长为$ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $的FPEI4-4格式和时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $的SRK2-4格式计算跨中挠度结果
Figure 5. The results of mid-span deflection calculated by analytical approximate solution, FPEI4-4 scheme with $ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $ and SRK2-4 scheme with $ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $
![]()
图 6 时间步长为$ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $和$ \tau =1.0 \times {10^{ - 5}}{\text{ s}} $时SRK2-4格式计算跨中挠度结果之差
Figure 6. The difference between the results of mid-span deflection calculated by SRK2-4 scheme under time step $ \tau =1.0 \times {10^{ - 4}}{\text{ s}} $and $ \tau =1.0 \times {10^{ - 5}}{\text{ s}} $
![]()
图 7 FPEI4-4格式和Newmark-β算法求解跨中挠度计算误差
Figure 7. The numerical errors of mid-span deflection calculated by FPEI4-4 scheme and Newmark-β algorithm
![]()
图 8 FPEI4-4格式与Newmark-β算法得到的桥梁跨中挠度结果
Figure 8. The results of mid-span deflection calculated by FPEI4-4 scheme and Newmark-β algorithm
![]()
图 11 $ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $, 不同车速时FPEI4-4格式求解跨中挠度计算误差
Figure 11. The numerical errors of mid-span deflection calculated by FPEI4-4 scheme under different vehicle speeds with time step $ \tau =1.0 \times {10^{ - 2}}{\text{ s}} $
![]()
图 12 $ \tau = 1.0 \times{10^{ - 2}}{\text{ s}}$, FPEI4-4格式选取不同时变矩阵冻结方式的计算结果
Figure 12. The numerical results calculated by FPEI4-4 scheme under different parameter freezing forms with time step $\tau = 1.0 \times {10^{ - 2}}{\text{ s}}$
![]()
图 13 $ \tau = 1.0 \times {10^{ - 2}}{\text{ s}}$, FPEI4-4格式不同时变矩阵冻结方式求解跨中挠度的计算误差
Figure 13. The numerical errors of mid-span deflection calculated by FPEI4-4 scheme under different parameter freezing forms with time step $\tau = 1.0 \times {10^{ - 2}}{\text{ s}}$
表 1 汽车参数
Table 1 Vehicle parameters
Symbol Value m2/kg 18000 m1/kg 2000 k2/(N·m−1) 9000000 c2/(N·s·m−1) 80000 k1/(N·m−1) 36000000 c1/(N·s·m−1) 700 
下载: 导出CSV
表 2 桥面沥青铺装层和桥体参数
Table 2 Pavement and bridge parameters
Symbol Value L/m 40 b/m 8 h1/m 0.1 h2/m 2 E1/GPa 8 E2/GPa 34.5 $ {\rho _1} $/(kg·m−3) 2000 $ {\rho _2} $/(kg·m−3) 2500 $ {\eta _1} $ 0.1 $ {\eta _2} $ 0.01
下载: 导出CSV
表 3 Newmark-β算法和FPEI4-4格式计算跨中挠度对比
Table 3 Comparison of mid-span deflection calculated by Newmark-β algorithm and FPEI4-4 scheme
Time step/s Computation time/s Maximum error/m Newmark-β algorithm $ {10^{ - 4}} $ 0.869262 $ - 2.222\;306 \times {10^{ - 7}} $ $ {10^{ - 3}} $ 0.078846 $ 2.269\;753 \times {10^{ - 6}} $ $ {10^{ - 2}} $ — — FPEI4-4 scheme $ {10^{ - 4}} $ 2.849518 $ 4.373\;681 \times {10^{ - 13}} $ $ {10^{ - 3}} $ 0.319489 $ 4.373\;730 \times {10^{ - 11}} $ $ {10^{ - 2}} $ 0.031098 $ 4.409\;290 \times {10^{ - 9}} $ $ {10^{ - 1}} $ 0.003465 $ - 1.741\;726 \times {10^{ - 6}} $
下载: 导出CSV
表 4 不同冰冻近似方式的对比
Table 4 Comparison of different parameter freezing forms
Time step/s Form Maximum error/m
$1.0 \times {10^{ - 2} }$LC1 $ 4.409\;290 \times {10^{ - 9}} $ LC2 $ 8.127\;175 \times {10^{ - 8}} $ LC3 $ 8.418\;344 \times {10^{ - 8}} $
$1.0 \times {10^{ - 3} }$LC1 $ 4.373\;730 \times {10^{ - 11}} $ LC2 $ 8.131\;892 \times {10^{ - 9}} $ LC3 $ 8.139\;901 \times {10^{ - 9}} $
下载: 导出CSV
-
[1] Zhu XQ, Law SS. Recent developments in inverse problems of vehicle-bridge interaction dynamics. Journal of Civil Structural Health Monitoring, 2016, 6(1): 107-128 doi: 10.1007/s13349-016-0155-x
[2] Yang YB, Yang JP. State-of-the-art review on modal identification and damage detection of bridges by moving test vehicles. International Journal of Structural Stability and Dynamics, 2018, 18(2): 1850025 doi: 10.1142/S0219455418500256
[3] 李小珍, 辛莉峰, 王铭等. 车−桥耦合振动2019年度研究进展. 土木与环境工程学报(中英文), 2020, 42(5): 126-138 (Li Xiaozhen, Xin Lifeng, Wang Ming, et al. State-of-the-art review of vehicle-bridge interactions in 2019. Journal of Civil and Environmental Engineering, 2020, 42(5): 126-138 (in Chinese) Li Xiaozhen, Xin Lifeng, Wang Ming, et al. State-of-the-art review of vehicle-bridge interactions in 2019. Journal of Civil and Environmental Engineering, 2020, 42(5): 126-138 (in Chinese)
[4] 周勇军, 薛宇欣, 李冉冉等. 桥梁冲击系数理论研究和应用进展. 中国公路学报, 2021, 35(4): 31-50 (Zhou Yongjun, Xue Yuxin, Li Ranran et al. State-of-the-art of theory and applications of bridge dynamic load allowance. China Journal of Highway and Transport, 2021, 35(4): 31-50 (in Chinese) Zhou Yongjun, Xue Yuxin, Li Ranran et al . State-of-the-art of theory and applications of bridge dynamic load allowance. China Journal of Highway and Transport,2021 ,35 (4 ):31 -50 (in Chinese)[5] 杨永斌, 王志鲁, 史康等. 基于车辆响应的桥梁间接测量与检测研究综述. 中国公路学报, 2021, 34(4): 1-12 (Yang Yongbin, Wang Zhilu, Shi Kang, et al. Research progress on bridge indirect measurement and monitoring from moving vehicle response. China Journal of Highway and Transport, 2021, 34(4): 1-12 (in Chinese) Yang Yongbin, Wang Zhilu, Shi Kang, et al . Research progress on bridge indirect measurement and monitoring from moving vehicle response. China Journal of Highway and Transport,2021 ,34 (4 ):1 -12 (in Chinese)[6] Han WS, Liu XD, Guo XL, et al. Research status and prospect of wind-vehicle-bridge coupling vibration system. Journal of Traffic and Transportation Engineering (English Edition), 2022, 9(3): 319-338 doi: 10.1016/j.jtte.2021.05.002
[7] Fryba L. Vibration of Solids and Structures under Moving Loads. London: Thomas Telford, 1999
[8] Yang YB, Yau JD, Hsu LC. Vibration of simple beams due to trains moving at high speeds. Engineering Structures, 1997, 19(11): 936-944 doi: 10.1016/S0141-0296(97)00001-1
[9] 夏禾, 张楠, 郭薇薇等. 车桥耦合振动工程. 北京: 科学出版社, 2014 (Xia He, Zhang Nan, Guo Weiwei, et al. Coupling Vibrations of Train-bridge System. Beijing: China Science Publishing and Media Ltd., 2014 (in Chinese) Xia He, Zhang Nan, Guo Weiwei, et al. Coupling Vibrations of Train-bridge System. Beijing: China Science Publishing and Media Ltd., 2014 (in Chinese)
[10] 张宇, 王嘉伟, 李韶华等. 基于智能驾驶电动汽车的多车−桥梁耦合系统动力学特性研究. 力学学报, 2022, 54(9): 1-13 (Zhang Yu, Wang Jiawei, Li Shaohua, et al. Dynamic characteristics of multi-vehicle-bridge coupling system based on intelligent driving electric vehicle. China Journal of Theoretical and Applied Mechanics, 2022, 54(9): 1-13 (in Chinese) Zhang Yu, Wang Jiawei, Li Shaohua, et al . Dynamic characteristics of multi-vehicle-bridge coupling system based on intelligent driving electric vehicle. China Journal of Theoretical and Applied Mechanics,2022 ,54 (9 ):1 -13 (in Chinese)[11] Ren JY, Chen YH, Sun ZQ, et al. A vehicle-bridge interaction vibration model considering bridge deck pavement. Journal of Low Frequency Noise, Vibration and Active Control, 2022, in press. doi: 10.1177/14613484221122736
[12] Camara A, Kavrakov I, Nguyen K, et al. Complete framework of wind-vehicle-bridge interaction with random road surfaces. Journal of Sound and Vibration, 2019, 458: 197-217 doi: 10.1016/j.jsv.2019.06.020
[13] 殷新锋, 晏万里, 仁厚乾等. 车载作用下公路桥梁耦合振动精细化建模及验证分析. 湖南大学学报, 2021, 48(9): 129-137 (Yin Xinfeng, Yan Wanli, Ren Houqian, et al. Fine modeling of coupled vibration of highway bridge under vehicle loading and verification analysis. Journal of Hunan University, 2021, 48(9): 129-137 (in Chinese) Yin Xinfeng, Yan Wanli, Ren Houqian, et al . Fine modeling of coupled vibration of highway bridge under vehicle loading and verification analysis. Journal of Hunan University,2021 ,48 (9 ):129 -137 (in Chinese)[14] 冯东明, 黎剑安, 吴刚等. 考虑路面不平度的车桥耦合系统快速建模与动力分析方法. 东南大学学报, 2022, 52(6): 1088-1094 (Feng Dongming, Li Jianan, Wu Gang, et al. Fast modeling and dynamic analysis method for vehicle-bridge interaction system considering road roughness. Journal of Southeast University, 2022, 52(6): 1088-1094 (in Chinese) Feng Dongming, Li Jianan, Wu Gang, et al . Fast modeling and dynamic analysis method for vehicle-bridge interaction system considering road roughness. Journal of Southeast University,2022 ,52 (6 ):1088 -1094 (in Chinese)[15] 王艳, 刘哲, 周瑞娇等. 基于车桥耦合振动分析的大跨径曲弦桁梁桥加固方案评价研究. 振动与冲击, 2022, 41(14): 199-209 (Wang Yan, Liu Zhe, Zhou Ruijiao, et al. Evaluation of the reinforcement scheme for a long-span curved truss bridge based on vehicle-bridge coupled vibration analysis. Journal of Vibration and Shock, 2022, 41(14): 199-209 (in Chinese) doi: 10.13465/j.cnki.jvs.2022.14.027 Wang Yan, Liu Zhe, Zhou Ruijiao, et al . Evaluation of the reinforcement scheme for a long-span curved truss bridge based on vehicle-bridge coupled vibration analysis. Journal of Vibration and Shock,2022 ,41 (14 ):199 -209 (in Chinese) doi: 10.13465/j.cnki.jvs.2022.14.027[16] 崔圣爱, 张猛, 曾光等. 基于车桥耦合振动的双层六线铁路桥梁疲劳性能评估. 东南大学学报, 2023, 53(1): 67-75 (Cui Shengai, Zhang Meng, Zeng Guang, et al. Fatigue performance assessment of double-deck six-line railway bridge based on vehicle-bridge coupling vibration. Journal of Southeast University, 2023, 53(1): 67-75 (in Chinese) Cui Shengai, Zhang Meng, Zeng Guang, et al . Fatigue performance assessment of double-deck six-line railway bridge based on vehicle-bridge coupling vibration. Journal of Southeast University,2023 ,53 (1 ):67 -75 (in Chinese)[17] Gao H, Yang BB. Parametric vibration of a flexible structure excited by periodic passage of moving oscillators. Journal of Applied Mechanics, 2020, 87: 071001 doi: 10.1115/1.4046781
[18] Konig P, Salcher P, Adam C. An efficient model for the dynamic vehicle-track-bridge-soil interaction system. Engineering Structures, 2022, 253: 113769 doi: 10.1016/j.engstruct.2021.113769
[19] 阳洋, 许文明, 卢会城等. 基于车桥耦合理论的梁式桥阻尼比识别研究. 力学学报, 2022, 54(5): 1387-1402 (Yang Yang, Xu Wenming, Lu Huicheng, et al. Research on damping ratio identification of beam bridge based on vehicle bridge coupling theory. China Journal of Theoretical and Applied Mechanics, 2022, 54(5): 1387-1402 (in Chinese) doi: 10.6052/0459-1879-21-691 Yang Yang, Xu Wenming, Lu Huicheng, et al . Research on damping ratio identification of beam bridge based on vehicle bridge coupling theory. China Journal of Theoretical and Applied Mechanics,2022 ,54 (5 ):1387 -1402 (in Chinese) doi: 10.6052/0459-1879-21-691[20] Zhang Q, Vrouwenvelder A, Wardenier J. Numerical simulation of train-bridge interactive dynamics. Computers and Structures, 2001, 79(10): 1059-1075 doi: 10.1016/S0045-7949(00)00181-4
[21] Zhu XQ, Law SS. Precise time-step integration for the dynamic response of a continuous beam under moving loads. Journal of Sound and Vibration, 2001, 240: 962-970 doi: 10.1006/jsvi.2000.3184
[22] Zhai WM. Two simple fast integration methods for large-scale dynamic problems in engineering. International Journal for Numerical Methods in Engineering, 1996, 39(24): 4199-4214 doi: 10.1002/(SICI)1097-0207(19961230)39:24<4199::AID-NME39>3.0.CO;2-Y
[23] Stoura CD, Paraskevopoulos E, Dimitrakopoulos EG, et al. A dynamic partitioning method to solve the vehicle-bridge interaction problem. Computers and Structures, 2021, 251: 106547 doi: 10.1016/j.compstruc.2021.106547
[24] 李韶华, 冯桂珍, 丁虎. 考虑胎路多点接触的电动汽车−路面耦合系统振动分析. 力学学报, 2021, 53(9): 2554-2568 (Li Shaohua, Feng Guizhen, Ding Hu. Vibration analysis of electric vehicle-road coupling system considering tire-road multi-point contact. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(9): 2554-2568 (in Chinese) Li Shaohua, Feng Guizhen, Ding Hu . Vibration analysis of electric vehicle-road coupling system considering tire-road multi-point contact. Chinese Journal of Theoretical and Applied Mechanics,2021 ,53 (9 ):2554 -2568 (in Chinese)[25] Liu K, Kujath MR. Response of slowly time-varying systems to harmonic excitation. Journal of Sound and Vibration, 1994, 177: 423-432 doi: 10.1006/jsvi.1994.1444
[26] Ge CQ, Tan CA, Lin K, et al. On the efficacy and accuracy of freezing technique for the response analysis of time-varying vehicle-bridge interaction systems. Journal of Sound and Vibration, 2021, 514: 116453 doi: 10.1016/j.jsv.2021.116453
[27] Sheng GG, Wang X. The geometrically nonlinear dynamic responses of simply supported beams under moving loads. Applied Mathematical Modelling, 2017, 48: 183-195 doi: 10.1016/j.apm.2017.03.064
[28] 钟万勰. 结构动力方程的精细时程积分法. 大连理工大学学报, 1994, 34(2): 131-136 (Zhong Wanxie. On precise time-integration method for structural dynamics. Journal of Dalian University of Technology, 1994, 34(2): 131-136 (in Chinese) Zhong Wanxie . On precise time-integration method for structural dynamics. Journal of Dalian University of Technology,1994 ,34 (2 ):131 -136 (in Chinese)[29] 高强, 谭述君, 钟万勰. 精细积分方法研究综述. 中国科学:技术科学, 2016, 46(12): 1207-1218 (Gao Qiang, Tan Shujun, Zhong Wanxie. A survey of the precise integration method. Scientia Sinica Technologica, 2016, 46(12): 1207-1218 (in Chinese) doi: 10.1360/N092016-00205 Gao Qiang, Tan Shujun, Zhong Wanxie . A survey of the precise integration method. Scientia Sinica Technologica,2016 ,46 (12 ):1207 -1218 (in Chinese) doi: 10.1360/N092016-00205[30] Tian W, Yang ZC, Gu YS. Dynamic analysis of an aeroelastic airfoil with freeplay nonlinearity by precise integration method based on Pade approximation. Nonlinear Dynamics, 2017, 89(3): 2173-2194 doi: 10.1007/s11071-017-3577-z
[31] Zhou XK, Duan ML, Granados JJ. Precise integration method for natural frequencies and mode shapes of ocean risers with elastic boundary conditions. Applied Mathematical Modelling, 2018, 61: 709-725 doi: 10.1016/j.apm.2018.05.017
[32] Zhu L, Ye KQ, Huang DW, et al. An adaptively filtered precise integration method considering perturbation for stochastic dynamics problems. Acta Mechanica Solida Sinica, 2023, 36: 317-326 doi: 10.1007/s10338-023-00381-4
[33] 顾元宪, 陈飚松, 张洪武. 结构动力方程的增维精细积分法. 力学学报, 2000, 32(4): 447-456 (Gu Yuanxian, Chen Biaosong, Zhang Hongwu. Precise time-integration with dimension expanding method. China Journal of Theoretical and Applied Mechanics, 2000, 32(4): 447-456 (in Chinese) Gu Yuanxian, Chen Biaosong, Zhang Hongwu . Precise time-integration with dimension expanding method. China Journal of Theoretical and Applied Mechanics,2000 ,32 (4 ):447 -456 (in Chinese)[34] 富明慧, 刘祚秋, 林敬华. 一种广义精细积分法. 力学学报, 2007, 39(5): 672-677 (Fu Minghui, Liu Zuoqiu, Lin Jinghua. A generalized precise time step integration method. China Journal of Theoretical and Applied Mechanics, 2007, 39(5): 672-677 (in Chinese) Fu Minghui, Liu Zuoqiu, Lin Jinghua . A generalized precise time step integration method. China Journal of Theoretical and Applied Mechanics,2007 ,39 (5 ):672 -677 (in Chinese)[35] Huang YA, Deng ZC, Yao LX. An improved symplectic precise integration method for analysis of the rotating rigid-flexible coupled system. Journal of Sound and Vibration, 2007, 299(1-2): 229-246 doi: 10.1016/j.jsv.2006.07.009
[36] Hu WP, Huai YL, Xu MB, et al. Mechanoelectrical flexible hub-beam model of ionic-type solvent-free nanofluids. Mechanical Systems and Signal Processing, 2021, 159: 107833 doi: 10.1016/j.ymssp.2021.107833
[37] Wu F, Gao Q, Zhong WX. Fast precise integration method for hyperbolic heat conduction problems. Applied Mathematics and Mechanics-English Edition, 2013, 34(7): 791-800 doi: 10.1007/s10483-013-1707-6
[38] Hochbruck M, Ostermann A. Exponential integrators. Acta Numerica, 2010, 19(1): 209-286
[39] Cox SM, Matthews PC. Exponential time differencing for stiff systems. Journal of Computational Physics, 2002, 176(2): 430-455 doi: 10.1006/jcph.2002.6995
[40] Gu XL, Jiang CL, Wang YS, et al. Efficient energy-preserving exponential integrators for multi-component Hamiltonian systems. Journal of Scientific Computing, 2022, 92(1): 26 doi: 10.1007/s10915-022-01874-z
[41] 邓子辰, 李庆军. 精细指数积分法在卫星编队飞行动力学中的应用. 北京大学学报, 2016, 52(4): 669-675 (Deng Zichen, Li Qingjun. Precise exponential integrator and its application in dynamics of spacecraft formation flying. Acta Scientiarum Naturalium Universitatis Pekinensis, 2016, 52(4): 669-675 (in Chinese) Deng Zichen, Li Qingjun . Precise exponential integrator and its application in dynamics of spacecraft formation flying. Acta Scientiarum Naturalium Universitatis Pekinensis,2016 ,52 (4 ):669 -675 (in Chinese)[42] Zhang JN, Yang SP, Li SH, et al. Study on crack propagation path of asphalt pavement under vehicle-road coupled vibration. Applied Mathematical Modelling, 2022, 101: 481-502 doi: 10.1016/j.apm.2021.09.004
[43] 李梦琪, 张锋, 冯德成等. 车桥耦合下钢桥面沥青铺装层动力响应研究. 工程力学, 2019, 36(12): 177-187 (Li Mengqi, Zhang Feng, Feng Decheng et al. Dynamic response of steel deck asphalt pavement considering vehicle-bridge coupling effect. Engineering Mechanics, 2019, 36(12): 177-187 (in Chinese) Li Mengqi, Zhang Feng, Feng Decheng et al . Dynamic response of steel deck asphalt pavement considering vehicle-bridge coupling effect. Engineering Mechanics,2019 ,36 (12 ):177 -187 (in Chinese)[44] Zhong WX. On precise integration method. Journal of Computational and Applied Mathematics, 2004, 163(1): 59-78 doi: 10.1016/j.cam.2003.08.053
[45] Lin K, Tan CA, Ge CQ, et al. Time-varying transmissibility analysis of vehicle-bridge interaction systems with application to bridge-friendly vehicles. International Journal of Structural Stability and Dynamics, 2022, 22(5): 2250022
[46] 邓子辰, 曹珊珊, 李庆军等. 基于辛Runge-Kutta方法的太阳帆塔动力学特性研究. 中国科学: 技术科学, 2016, 46(12): 1242-1253 (Deng Zichen, Cao Shanshan, Li Qingjun et al. Dynamic behavior of sail tower SPS based on the symplectic Runge-Kutta method. Scientia Sinica Technologica, 2016, 46(12): 1242-1253 (in Chinese) doi: 10.1360/N092016-00163 Deng Zichen, Cao Shanshan, Li Qingjun et al . Dynamic behavior of sail tower SPS based on the symplectic Runge-Kutta method. Scientia Sinica Technologica,2016 ,46 (12 ):1242 -1253 (in Chinese) doi: 10.1360/N092016-00163[47] Yau JD, Yang YB, Kuo SR. Impact response of high speed rail bridges and riding comfort of rail cars. Engineering Structures, 1999, 21(9): 836-844 doi: 10.1016/S0141-0296(98)00037-6
-
期刊类型引用(1)
1. Yu ZHANG，Chao ZHANG，Shaohua LI，Shaopu YANG. Research on the application of the parameter freezing precise exponential integrator in vehicle-road coupling vibration. Applied Mathematics and Mechanics(English Edition). 2025(02): 373-390 . 
必应学术
其他类型引用(2)
计量
- 文章访问数: 391
- HTML全文浏览量: 112
- PDF下载量: 68
- 被引次数: 3
