汽车工程学报

模态数	仿真时间/h	梁长/m	仿真时间/h
5	0.33	30	5.71
10	1.85	60	8.28
20	5.25	80	10.19
30	7.82	120	15.75

模态数	仿真时间/h	梁长/m	仿真时间/h
5	0.33	30	5.71
10	1.85	60	8.28
20	5.25	80	10.19
30	7.82	120	15.75

基于无限长梁模型的车路振动耦合系统半解析解

PDF下载

赵妍 ¹ , 张步云 ¹ , 刘志强 ³ , 胡正林 ⁴

汽车工程学报 | 系统动力学专栏 2024,14(2): 217-225

收起

汽车工程学报 | 系统动力学专栏 2024, 14(2): 217-225

基于无限长梁模型的车路振动耦合系统半解析解

全屏

赵妍¹, 张步云¹, 刘志强³, 胡正林⁴

作者信息

¹ 江苏大学汽车工程研究院镇江 212013

² 江苏大学汽车与交通工程学院镇江 212013

³ 江苏恒力制动器制造有限公司靖江 214500

赵妍(1996-),女,河南周口人,硕士研究生,主要研究方向为车路耦合振动分析、车辆舒适性和道路友好性优化与控制。Tel: 18252586072 E-mail：2391850117@qq.com

通讯作者:

张步云(1987-),男,江苏淮安人,博士,副教授,主要研究方向为振动测试与控制、车辆随机动力学分析。Tel: 13914560131 E-mail：zhangby@ujs.edu.cn

Semi-Analytical Solution for Vehicle-Road Vibrational Coupling System Based on an Infinite-Length Beam Model

Yan ZHAO¹, Buyun ZHANG¹, Zhiqiang LIU³, Zhenglin HU⁴

Affiliations

¹ Automotive Engineering Research Institute Jiangsu University Zhenjiang 212013 China

² Department of Mechanical Engineering Wayne State University Detroit 48202 USA

³ School of Automotive and Traffic Engineering Jiangsu University Zhenjiang 212013 China

doi: 10.3969/j.issn.2095–1469.2024.02.06

文章导航

摘要

收起

针对无限长道路与车辆耦合系统响应计算复杂难题,考虑地基的弹性特性与道路不平度,建立基于无限长欧拉伯努利梁模型的车路振动耦合系统。进而以车辆为参考点建立移动坐标系,提出通过积分变换推导耦合系统振动响应解析解的方法,并应用留数定理对其进行数值计算,获得车辆垂向位移、加速度、路面振动响应等系统响应的半解析解。与传统应用模态叠加法的有限长道路与车辆耦合响应相比,具有更高的计算效率与精度,系统参数化研究也证明了该半解析解的有效性。

关键词

车路振动耦合 / 无限长梁模型 / 积分变换 / 留数定理 / 半解析解

Abstract

收起

To solve the complex challenge of response calculation for the coupling system between the infinitelength road and vehicle, the elastic characteristics of foundation and road roughness are considered in the analysis, and a vehicleroad vibration coupling system is established based on an infinite length EulerBernoulli beam model. Then, the moving coordinate system was set up using the vehicle as the reference point. The analytical solution of vibration response of the coupled system was derived by integral transformation. The numerical calculations were carried out by applying the residue theorem, and the semianalytical solutions for the vehicle's vertical displacement, acceleration and road vibration response were obtained. Compared to the traditional modal superposition method used for the coupling response of finite

Key words

vehicle-road dynamic coupling / infinite beam model / integral transformation / residue theorem / semi-analytic solution

引用本文

赵妍, 张步云, 刘志强, 胡正林. 基于无限长梁模型的车路振动耦合系统半解析解. 汽车工程学报, 2024 , 14 (2) : 217 -225 . DOI: 10.3969/j.issn.2095–1469.2024.02.06

Yan ZHAO, Buyun ZHANG, Zhiqiang LIU, Zhenglin HU. Semi-Analytical Solution for Vehicle-Road Vibrational Coupling System Based on an Infinite-Length Beam Model[J]. Chinese Journal of Automotive Engineering, 2024 , 14 (2) : 217 -225 . DOI: 10.3969/j.issn.2095–1469.2024.02.06

正文

收起

车辆对道路的作用影响分析可归为移动载荷问题 (Moving Load Problems, MLP) 研究, 其目的在于通过建立车辆与道路的运动微分方程, 进而求解系统的振动响应。将车辆作用在道路上的力简化为移动集中力 ^{[ 1 ]} 、移动线载荷 ^{[ 2 ]} 、移动分布载荷 ^{[ 3 ]} 比较常见, 但其并不能准确地描述车辆与道路互相的耦合作用。因此,建立四分之一车辆 ^{[ 4 ]} 、半车或整车模型 ^{[ 5 ]} 就显得尤为必要,考虑车辆与道路的振动耦合对精确获得车辆振动响应 ^{[ 6 ]} 、悬架系统优化设计 ^{[ 7 ]} 以及评估对道路的损伤 ^{[ 8 ]} 具有更重要的意义,受到了国内外研究人员的广泛关注 ^{[ 5 , 9 - 11 ]} , 其中建立合理准确的车路振动耦合模型是关键内容。

车辆通常采用离散自由度模型 ^{[ 12 ]} 进行合理简化, 将路面不平度以及车路之间的耦合振动视为垂向的主要激励。道路包含路面和地基两部分, 建模较复杂。路面通常建模为一维梁模型、二维板模型或三维空间模型, 地基有 Winkler 模型、Pasternak 模型、变刚度地基或弹性半空间地基模型等 ^{[ 13 ]} 。其中梁模型应用较广泛, 分为有边界的有限长梁模型与无限长梁模型。对于有限长梁与车辆耦合的振动研究,已有较多文献论述。SUN Lu 等 ^{[ 2 ]} 分析了弹性 Winkler 地基上受移动线荷载作用的梁响应, 通过傅里叶变换将梁的控制方程转化为代数方程, 并利用复变函数理论将解转化为时域解, 得到了移动线荷载的临界速度。MALLIK 等 ^{[ 14 ]} 研究了弹性地基上的梁在匀速运动集中荷载作用下的稳态响应,得出了梁的最大位移和弯矩。DING Hu 等 ^{[ 15 ]} 研究了非线性弹性地基上梁在集中荷载作用下的动力响应 Galerkin 方法的收敛性, 分析不同边界条件和跨度对收敛性和动力响应影响。KRISHNANUNNI 等 ^{[ 16 ]} 结合 Galerkin 方法、有限差分格式和 Newmark 方法, 提出了车辆-路面系统动力响应的迭代解耦方法。而由于有限长梁模型的局限性, 越来越多的学者考虑利用无限长梁进行建模。

基于无限长梁模型的车路耦合振动研究, 重点在于车辆模型的简化与响应计算方法的研究两方面。如 KARGARNOVIN 等 ^{[ 17 ]} 研究了非线性粘弹性地基上无限长梁在简谐移动荷载作用下的响应, 提出了一种适用于频域的求解方法, 涉及格林函数和留数定理用于积分的计算。丁虎等 ^{[ 18 ]} 研究了三次非线性粘弹性地基上无限长 Timoshenko 梁在运动集中力作用下的响应, 采用阿多米亚分解法和摄动法结合复傅里叶变换求解微分方程, 推导出近似的积分形式的封闭解。JANG 等 ^{[ 19 ]} 基于几何非线性 Karman 近似, 建立了非线性积分方程组, 提出了一种分析局部荷载作用下弹性地基上无限长非线性梁中等大挠度的新方法。YU Haitao 等 ^{[ 20 ]} 建立了粘弹性地基上无限长梁在任意动荷载作用下的动力响应解析解,并验证了所提一般解的有效性。由此看出, 学者们提出了多种作用力简化模型与响应计算方法, 且多在频域中开展分析, 对于时域解的计算仍存在较大困难, 其关键在于精确的解析解求解困难较大。ZHEN Bin 等 ^{[ 21 ]} 研究了无限长欧拉伯努利梁在非线性移动荷载作用下的稳态响应, 得到了不包含复杂无穷积分的封闭形式解。ZHANG Yin 等 ^{[ 3 ]} 讨论了粘弹性地基上具有移动分布质量荷载的无限长梁的稳态响应, 得到稳态解的半解析解。 FROIO 等 ^{[ 22 ]} 分析了 Pasternak 弹性地基上的均匀无限长欧拉-伯努利梁在匀速运动的集中载荷作用下的稳态响应, 通过傅里叶变换导出了一个通用的封闭形式解析解。DIMITROVOVÁ ^{[ 4 ]} 用积分变换和轮廓积分法, 推导并验证了由双参数粘弹性地基支撑下的二自由度系统在无限梁上运动引起的横向振动的半解析解。这表明对于无限长梁模型在移动载荷作用下而言, 半解析解是一种更加合适的分析方法。但对于无限长梁与车辆振动耦合, 尤其是考虑路面不平度与振子等车辆模型的耦合振动研究, 目前仍为鲜见。

综上所述, 基于有限长梁的车路耦合研究在数值求解过程中除了耗费较大的计算成本, 还需分析结果的收敛性。而传统基于无限长梁的研究虽然考虑了耦合的作用, 但车辆大多数被简化为已知移动力或移动载荷, 且未考虑道路表面不平度的影响, 难以表示车辆与道路真实的耦合作用情况。因此, 本文针对无限长道路与车辆耦合系统响应计算复杂难题, 考虑地基的弹性特性与道路不平度, 建立基于无限欧拉-伯努利梁模型的车路振动耦合系统。进而以车辆为参考点建立移动坐标系, 提出通过积分变换推导耦合系统振动响应解析解的方法, 并应用留数定理对其进行数值分析处理, 获得车辆垂向位移、加速度、路面振动响应等系统响应的半解析解。

1 耦合系统模型及运动微分方程

收起

本文建立的车路耦合振动系统模型, 遵循以下假设:

1)采用无限长的欧拉-伯努利梁；

2)将车辆悬架简化为垂向振子模型;

3)梁的截面恒定,材料均匀,单位长度质量恒定;

4)梁满足小挠度理论假设；

5)车辆匀速行驶,由左向右,如图 1 所示。

图 1 中: $m$ 为车身质量; $z\left( t\right)$ 为车身垂直运动位移; ${k}_{t}$ 和 ${c}_{t}$ 分别为车辆悬架的弹簧刚度和阻尼; $\xi \left( x\right)$ 为路面不平度; ${k}_{f}$ 为弹性地基的刚度; $\bar{w}\left( {x, t}\right)$ 为道路横向变形,车辆行驶速度为 $v$ 。根据牛顿第二定律, 系统运动微分方程为:

(1)

$ m\ddot{z}\left( t\right) = - {f}_{a}\left( t\right) \text{。} $

(2)

$ {\rho A}{\bar{w}}_{,{tt}}\left( {x, t}\right) + \mathrm{{EI}}{\bar{w}}_{,{xxxx}}\left( {x, t}\right) + {k}_{f}\bar{w}\left( {x, t}\right) \\ = {f}_{a}\left( t\right) \delta \left( {x - {x}_{c}}\right) \text{。} \\ {f}_{a}\left( t\right) = {c}_{t}\left\lbrack {\dot{z}\left( t\right) - \dot{\xi }\left( {x}_{c}\right) - {\bar{w}}_{, t}\left( {{x}_{c}, t}\right) }\right\rbrack + {k}_{t}\lbrack z\left( t\right) - \\ \left. {\xi \left( {x}_{c}\right) - \bar{w}\left( {{x}_{c}, t}\right) }\right\rbrack \text{。} $

0197d13f-5613-7ed1-96d5-af28cf465cce_2_943_231_580_461_0.jpg

图 1 车辆与有限长梁及无限长梁振动耦合系统模型

式中: EI 为梁的弯曲刚度; ${\rho A}$ 为单位长度梁的质量; ${x}_{c}$ 为车辆行驶的距离, ${x}_{c} = {vt}$ ; 变量上的点 ( ) 为关于时间的导数算子,其下标(ot)和 (, ${xxxx}$ ) 分别为偏微分方程中关于时间 $t$ 和距离 $x$ 的偏导数算子; ${f}_{a}\left( t\right)$ 为车辆悬架受到来自于路面的作用力; ${f}_{a}\left( t\right) \delta \left( {x - {x}_{c}}\right)$ 表示在车辆所在处对于道路的作用力。

梁的长度为无限长, 其边界条件和初始条件为:

(4)

$ \mathop{\lim }\limits_{{x \rightarrow \pm \infty }}\bar{w}\left( {x, t}\right) = 0,\;\mathop{\lim }\limits_{{x \rightarrow \pm \infty }}\frac{{\partial }^{k}\bar{w}\left( {x, t}\right) }{\partial {x}^{k}} = 0\left( {k = 1,2,3}\right) 。 $

(5)

$ {\left. \bar{w}\left( x, t\right) \right| }_{t = 0} = 0,{\left. {\bar{w}}_{, t}\left( x, t\right) \right| }_{t = 0} = {0}_{ \circ } $

当车辆到达 ${x}_{c}$ 处时, $\xi \left( {x}_{c}\right) = {P}_{0}\sin \left( {{2\pi }{x}_{c}/\lambda }\right)$ , 其中, $\lambda$ 是波长, ${P}_{0}$ 是路面不平度幅值。将 ${x}_{c} = {vt}$ 代入 $\xi \left( {x}_{c}\right)$ 可得时域上的路面激励:

(6)

$ \bar{\xi }\left( t\right) = {P}_{0}\sin \left( {{2\pi vt}/\lambda }\right) \text{。} $

可知,激励频率为 ${\omega }_{\xi } = {2\pi v}/\lambda$ 。车路耦合系统的振动特性是由路面不平度激发的, 但与单纯恒移动载荷问题不同, 当考虑路面的柔性特点时, 耦合作用力可视为车辆与无限长道路系统的内力。为更加清晰地描述系统的特性, 以车辆为参考点建立移动的坐标系,令 $r = x - {x}_{c}$ 。继而可将梁的横向变形 $\bar{w}\left( {x, t}\right)$ 替换为:

(7)

$ w\left( {r, t}\right) = w\left( {x - {x}_{c}, t}\right) = w\left( {x - {vt}, t}\right) = \bar{w}\left( {x, t}\right) 。 $

式(2)和式(3)中相关表达式变为:

(8)

$ {\bar{w}}_{, t}\left( {x, t}\right) = - v{w}_{, r}\left( {r, t}\right) + {w}_{, t}\left( {r, t}\right) 。 $

(9)

$ {\bar{w}}_{, u}\left( {x, t}\right) = {v}^{2}{w}_{,{rr}}\left( {r, t}\right) - {2v}{w}_{,{rt}}\left( {r, t}\right) + {w}_{, u}\left( {r, t}\right) 。 $

(10)

$ {\bar{w}}_{, x}\left( {x, t}\right) = {w}_{, r}\left( {r, t}\right) ,{\bar{w}}_{,{xxxx}}\left( {x, t}\right) = {w}_{,{rrrr}}\left( {r, t}\right) 。 $

(11)

$ \bar{w}\left( {{x}_{c}, t}\right) = w\left( {0, t}\right) ,{\bar{w}}_{, t}\left( {{x}_{c}, t}\right) = {w}_{, t}\left( {0, t}\right) 。 $

将式 (6)、式 (7) $\sim$ (11) 分别代入式 (2) ~ (3), 得到梁在移动坐标系下的偏微分运动方程:

(12)

$ \mathrm{{EI}}{w}_{,{rrrr}}\left( {r, t}\right) + {\rho A}{v}^{2}{w}_{,{rr}}\left( {r, t}\right) - {2\rho Av}{w}_{,{rt}}\left( {r, t}\right) + \\ {k}_{f}w\left( {r, t}\right) + {\rho A}{w}_{, u}\left( {r, t}\right) = {f}_{s}\left( {r, t}\right) \text{。} \\ {f}_{s}\left( {r, t}\right) = \left\{ {{c}_{t}\left\lbrack {\dot{z}\left( t\right) - {\omega }_{\xi }{P}_{0}\cos \left( {{\omega }_{\xi }t}\right) - {w}_{, t}\left( {0, t}\right) }\right\rbrack + }\right. \\ \left. {{k}_{t}\left\lbrack {z\left( t\right) - {P}_{0}\sin \left( {{\omega }_{\xi }t}\right) - w\left( {0, t}\right) }\right\rbrack }\right\} \delta \left( r\right) 。 $

由于在移动坐标系中车辆为坐标原点, 所以对于车辆悬架而言 $r = 0$ 。因此,车辆垂直运动位移变为:

(14)

$ m\ddot{z}\left( t\right) + {c}_{t}\left\lbrack {\dot{z}\left( t\right) - {\omega }_{\xi }{P}_{0}\cos \left( {{\omega }_{\xi }t}\right) - {w}_{, t}\left( {0, t}\right) }\right\rbrack + \\ {k}_{t}\left\lbrack {z\left( t\right) - {P}_{0}\sin \left( {{\omega }_{\xi }t}\right) - w\left( {0, t}\right) }\right\rbrack = 0\text{ 。 } $

相应地, 初始条件和边界条件可写为:

(15)

$ \mathop{\lim }\limits_{{r \rightarrow \pm \infty }}w\left( {r, t}\right) = 0,\;\mathop{\lim }\limits_{{r \rightarrow \pm \infty }}\frac{{\partial }^{i}w\left( {r, t}\right) }{\partial {r}^{i}} = 0\left( {i = 1,2,3}\right) 。 $

(16)

$ {\left. w\left( r, t\right) \right| }_{t = 0} = 0,{\left. {w}_{, t}\left( r, t\right) \right| }_{t = 0} = {0}_{ \circ } $

2 微分方程的积分变换

收起

为求解微分方程, 简化运算, 方便对谐波信号进行分析计算, 运用积分变换将微积分公式变换成代数公式, 可用于分析系统的性质, 判定系统的稳定性, 实现时域和频域的变换。

2.1 时间域上的拉普拉斯变换

对于任一函数 $g\left( {r, t}\right)$ ,定义 $\widetilde{G}\left( {r, s}\right)$ 表示它对时间域 $t$ 的拉普拉斯变换:

(17)

$ \widetilde{G}\left( {r, s}\right) = {\int }_{0}^{\infty }g\left( {r, t}\right) {e}^{-{st}}\mathrm{\;d}t $

将其应用于运动微分方程式 (12) $\sim \left( {14}\right)$ , 代入初始条件, 可得:

(18)

$ \mathrm{{EI}}{\widetilde{W}}_{,{rrrr}}\left( {r, s}\right) + {\rho A}{v}^{2}{\widetilde{W}}_{,{rr}}\left( {r, s}\right) - {2\rho Avs}{\widetilde{W}}_{, r}\left( {r, s}\right) + \\ \left( {{k}_{f} + {\rho A}{s}^{2}}\right) \widetilde{W}\left( {r, s}\right) = {\widetilde{F}}_{s}\left( {r, s}\right) 。 \\ {\widetilde{F}}_{s}\left( {r, s}\right) = \left\lbrack {\left( {{c}_{t}s + {k}_{t}}\right) \widetilde{Z}\left( s\right) - {P}_{0}{\omega }_{\xi }\frac{{c}_{t}s + {k}_{t}}{{s}^{2} + {\omega }_{\xi }^{2}} - \left( {{c}_{t}s + }\right. }\right. \\ \left. {\left. {k}_{t}\right) \widetilde{W}\left( {0, s}\right) }\right\rbrack \delta \left( r\right) \text{。} \\ \left( {m{s}^{2} + {c}_{t}s + {k}_{t}}\right) \widetilde{Z}\left( s\right) - {P}_{0}{\omega }_{\xi }\frac{{c}_{t}s + {k}_{t}}{{s}^{2} + {\omega }_{\xi }^{2}}\left( {{c}_{t}s + }\right. \\ \left. {k}_{t}\right) \widetilde{W}\left( {0, s}\right) = 0\text{。} $

### 2.2 空间域上的傅里叶变换

对于函数 $\widetilde{G}\left( {r, s}\right)$ ,定义 $\overline{\widetilde{G}}\left( {q, s}\right)$ 表示其对空间域 $r$ 的傅里叶变换:

(21)

$ \overline{\widetilde{F}}\left( {q, s}\right) = {\int }_{-\infty }^{\infty }\widetilde{F}\left( {r, s}\right) {e}^{-{iqr}}\mathrm{\;d}r。 $

将其应用到式(18) $\sim \left( {20}\right)$ 中,代入初始边界条件, 得:

(22)

$ \mathrm{{EI}}{q}^{4}\overline{\widetilde{W}}\left( {q, s}\right) - {\rho A}{v}^{2}{q}^{2}\overline{\widetilde{W}}\left( {q, s}\right) - {2\rho Avsiq}\overline{\widetilde{W}}\left( {q, s}\right) + \\ \left( {{k}_{f} + {\rho A}{s}^{2}}\right) \overline{\widetilde{W}}\left( {q, s}\right) = {\overline{\widetilde{F}}}_{s}\left( {0, s}\right) 。 \\ {\overline{\widetilde{F}}}_{s}\left( {0, s}\right) = \left( {{c}_{t}s + {k}_{t}}\right) \widetilde{Z}\left( s\right) - {P}_{0}{\omega }_{\xi }\frac{{c}_{t}s + {k}_{t}}{{s}^{2} + {\omega }_{\xi }^{2}}\left( {{c}_{t}s + }\right. \\ \left. {k}_{t}\right) \widetilde{W}\left( {0, s}\right) \text{。} $

由上式可得道路横向运动的双重积分变换为:

(24)

$ \overline{\widetilde{W}}\left( {q, s}\right) = \frac{{\overline{\widetilde{F}}}_{s}\left( {0, s}\right) }{D\left( {q, s}\right) }。 $

其中

(25)

$ D\left( {q, s}\right) = \mathrm{{EI}}{q}^{4} - {\rho A}{v}^{2}{q}^{2} - {2\rho Avsiq} + {k}_{f} + {\rho A}{s}^{2}\text{。} $

2.3 空间域上的逆傅里叶变换

式 (21) 的傅里叶逆变换为:

(26)

$ \widetilde{F}\left( {r, s}\right) = \frac{1}{2\pi }{\int }_{-\infty }^{\infty }\overline{\widetilde{F}}\left( {q, s}\right) {e}^{iqr}\mathrm{\;d}q $

因此, 对式 (24) 应用傅里叶逆变换得:

(27)

$ \widetilde{W}\left( {r, s}\right) = \frac{1}{2\pi }{\int }_{-\infty }^{\infty }\overline{\widetilde{W}}\left( {q, s}\right) {e}^{iqr}\mathrm{\;d}q = \frac{1}{2\pi }{\int }_{-\infty }^{\infty }\left\lbrack \left( {{c}_{t}s + }\right. \right. \\ \left. {k}_{t}\right) \widetilde{Z}\left( s\right) - {P}_{0}{\omega }_{\xi }\frac{{c}_{t}s + {k}_{t}}{{s}^{2} + {\omega }_{\xi }^{2}} - \left( {{c}_{t}s + }\right. \\ \left. {\left. {k}_{t}\right) \widetilde{W}\left( {0, s}\right) }\right\rbrack \frac{{e}^{iqr}\mathrm{\;d}q}{D\left( {q, s}\right) } \circ $

此式为拉普拉斯域的解。为消去 $\widetilde{W}\left( {0, s}\right)$ 项, 引入 $r = 0$ ,得:

(28)

$ \widetilde{W}\left( {0, s}\right) = \frac{1}{2\pi }\left\lbrack {\left( {{c}_{t}s + {k}_{t}}\right) \widetilde{Z}\left( s\right) - {P}_{0}{\omega }_{\xi }\frac{{c}_{t}s + {k}_{t}}{{s}^{2} + {\omega }_{\xi }^{2}} - }\right. \\ \left. {\left( {{c}_{t}s + {k}_{t}}\right) \widetilde{W}\left( {0, s}\right) }\right\rbrack {\int }_{-\infty }^{\infty }\frac{\mathrm{d}q}{D\left( {q, s}\right) }。 $

令

(29)

$ K\left( s\right) = {\int }_{-\infty }^{\infty }\frac{\mathrm{d}q}{D\left( {q, s}\right) }, K\left( {r, s}\right) = {\int }_{-\infty }^{\infty }\frac{{e}^{iqr}}{D\left( {q, s}\right) }\mathrm{d}q。 $

将未知变量 $\widetilde{W}\left( {0, s}\right)$ 和 $\widetilde{Z}\left( s\right)$ 代入式 (28),得:

(30)

$ \left\lbrack {{2\pi } + K\left( s\right) \left( {{c}_{t}s + {k}_{t}}\right) }\right\rbrack \widetilde{W}\left( {0, s}\right) - K\left( s\right) \left( {{c}_{t}s + }\right. \\ \left. {k}_{t}\right) \widetilde{Z}\left( s\right) = - K\left( s\right) {P}_{0}{\omega }_{\xi }\frac{{c}_{t}s + {k}_{t}}{{s}^{2} + {\omega }_{\xi }^{2}} \circ $

继而将其与式(20)联立,可得关于 $\widetilde{W}\left( {0, s}\right)$ 和 $\widetilde{Z}\left( s\right)$ 的方程:

(31)

$ {\Delta }_{wz} = \left\lbrack {{2\pi } + K\left( s\right) \left( {{c}_{t}s + {k}_{t}}\right) }\right\rbrack \left( {m{s}^{2} + {c}_{t}s + {k}_{t}}\right) - \\ K\left( s\right) {\left( {c}_{t}s + {k}_{t}\right) }^{2} = K\left( s\right) \left( {{c}_{t}s + {k}_{t}}\right) m{s}^{2} + {2\pi }\left( {m{s}^{2} + }\right. \\ \left. {{c}_{t}s + {k}_{t}}\right) \text{。} \\ \widetilde{W}\left( {0, s}\right) = \frac{-K\left( s\right) {P}_{0}{\omega }_{\xi }\left( {{c}_{t}s + {k}_{t}}\right) m{s}^{2}}{{\Delta }_{wz}\left( {{s}^{2} + {\omega }_{\xi }^{2}}\right) }。 \\ \widetilde{Z}\left( s\right) = \frac{{2\pi }{P}_{0}{\omega }_{\xi }\left( {{c}_{t}s + {k}_{t}}\right) }{{\Delta }_{wz}\left( {{s}^{2} + {\omega }_{\xi }^{2}}\right) } \circ $

将式(32)和式(33)代入方程式(27),得 $\widetilde{W}\left( {r, s}\right)$ 的表达式为:

(34)

$ \widetilde{W}\left( {r, s}\right) = \frac{-K\left( {r, s}\right) }{{\Delta }_{wz}\left( {{s}^{2} + {\omega }_{\xi }^{2}}\right) }m{P}_{0}{\omega }_{\xi }\left( {{c}_{t}s + {k}_{t}}\right) {s}^{2}\text{。} $

2.4 时间域上的逆拉普拉斯变换

式 (17) 的拉普拉斯逆变换为:

(35)

$ f\left( t\right) = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{a - {iT}}^{a + {iT}}\widetilde{F}\left( s\right) {e}^{st}\mathrm{\;d}s $

因此, 对式 (33) 和 (34) 应用拉普拉斯逆变换可分别得到车辆和无限长梁在时域内的响应。

(36)

$ w\left( {r, t}\right) = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{a - {iT}}^{a + {iT}}\widetilde{W}\left( {r, s}\right) {e}^{st}\mathrm{\;d}s \\ = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{a - {iT}}^{a + {iT}}\frac{{\omega }_{\xi }K\left( {r, s}\right) }{{\Delta }_{wz}\left( {{s}^{2} + {\omega }_{\xi }^{2}}\right) }m{P}_{0}\left( {{c}_{t}s + {k}_{t}}\right) {s}^{2}{e}^{st}\mathrm{\;d}s \circ $

(37)

$ z\left( t\right) = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{a - {iT}}^{a + {iT}}\frac{{2\pi }{P}_{0}{\omega }_{\xi }\left( {{c}_{t}s + {k}_{t}}\right) }{{\Delta }_{wz}\left( {{s}^{2} + {\omega }_{\xi }^{2}}\right) }{e}^{st}\mathrm{\;d}s \circ $

式中: $s = {i\omega };t > 0;a$ 为大于 0 的实数。

当 $a$ 是正实数且必须大于所有奇点的实部时, 被积函数收敛。将 $s = {i\omega }$ 代入式(25)、(29)和 (36),则 $D\left( {q, s}\right)$ 转化为 $D\left( {q,\omega }\right) , K\left( s\right)$ 转化为 $K\left( \omega \right) , K\left( {r, s}\right)$ 转化为 $K\left( {r,\omega }\right)$ ,涉及的相关公式转化为:

(38)

$ D\left( {q,\omega }\right) = \mathrm{{EI}}{q}^{4} - {\rho A}{v}^{2}{q}^{2} + {2\rho Av\omega q} - {\rho A}{\omega }^{2} + {k}_{f} \circ $

(39)

$ K\left( \omega \right) = {\int }_{-\infty }^{\infty }\frac{1}{D\left( {q,\omega }\right) }\mathrm{d}q, K\left( {r,\omega }\right) {\int }_{-\infty }^{\infty }\frac{{e}^{iqr}}{D\left( {q,\omega }\right) }\mathrm{d}q。 $

(40)

$ {\Delta }_{wz}\left( \omega \right) = - K\left( \omega \right) \left( {{c}_{t}{i\omega } + {k}_{t}}\right) m{\omega }^{2} + {2\pi }\left( {-m{\omega }^{2} + }\right. \\ \left. {{c}_{t}{i\omega } + {k}_{t}}\right) \text{。} $

综上所述, 得车路耦合系统振动响应的解析解:

(41)

$ w\left( {r, t}\right) = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{a - {iT}}^{a + {iT}}\widetilde{W}\left( {r, s}\right) {e}^{st}\mathrm{\;d}\left( s\right) 。 $

(42)

$ = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{-{ia} - T}^{-{ia} + T}i\widetilde{W}\left( {r,\omega }\right) {e}^{i\omega t}\mathrm{\;d}\omega \\ = \mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{-{ia} - T}^{-{ia} + T}\frac{-i{\omega }_{\xi }m{P}_{0}\left( {{c}_{t}{i\omega } + {k}_{t}}\right) {\omega }^{2}K\left( {r,\omega }\right) {e}^{i\omega t}\mathrm{\;d}\omega }{{2\pi i}{\Delta }_{wz}\left( {{\omega }_{\xi }^{2} - {\omega }^{2}}\right) } \\ z\left( t\right) = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{a - {iT}}^{a + {iT}}\widetilde{Z}\left( s\right) {e}^{st}\mathrm{\;d}\left( s\right) \\ = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{-{ia} - T}^{-{ia} + T}i\widetilde{Z}\left( \omega \right) {e}^{i\omega t}\mathrm{\;d}\omega \\ = \frac{1}{2\pi i}\mathop{\lim }\limits_{{T \rightarrow \infty }}{\int }_{-{ia} - T}^{-{ia} + T}i\frac{{2\pi }{P}_{0}{\omega }_{\xi }\left( {i{c}_{t}\omega + {k}_{t}}\right) }{{\Delta }_{wz}\left( {{\omega }_{\xi }^{2} - {\omega }^{2}}\right) }{e}^{i\omega t}\mathrm{\;d}\omega \text{。} $

## 3 留数定理在半解析解中的应用

为计算积分形式的解析解, 本文应用留数定理,需要求出被积函数 $i\widetilde{W}\left( {r,\omega }\right)$ 和 $i\widetilde{Z}\left( \omega \right)$ 的极点, 这两个被积函数具有相同的分母 ${\Delta }_{wz}\left( {{\omega }_{\xi }^{2} - {\omega }^{2}}\right) ,{\Delta }_{wz}$ 的定义如式(40)所示。显然, 分母为 0 时, 其中两个根为 $\pm {\omega }_{\xi }$ ,其他根定义为 ${\omega }_{i}$ 。因此,如果 ${\omega }_{i1}$ 是 ${\Delta }_{wz} = 0$ 的根,那么 ${\omega }_{i2} = {\left( -{\omega }_{i1}\right) }^{ * }$ 也是 ${\Delta }_{wz} = 0$ 根, 其中 $*$ 表示共轭,可知频率 ${\omega }_{i}$ 的数量是偶数个,这也表明一对频率 ${\omega }_{i}$ 总是有相同的虚部和相反的实部。解析解中 ${e}^{i\omega t}$ 表明,频率 ${\omega }_{i}$ 的实部将形成振动的谐波项,而虚部将指示频率 ${\omega }_{i}$ 引起的振动的振幅是否会随时间逐渐停止, 变得不稳定或保持不变, 这将分别发生在正、负和零虚部的情况下。所以在不考虑瞬态解的情况下, 系统的时域解可以写成对应于激励频率 ${\omega }_{\xi }$ 的谐波项与对应于频率 ${\omega }_{i}$ 的谐波项之和。

根据以上分析, 在不考虑瞬态解的情况下, 无限长梁的挠度变化和车辆的振动情况可由 $w\left( {r, t}\right)$ 和 $z\left( t\right)$ 给出:

(43)

$ w\left( {r, t}\right) = \sum \operatorname{res}\left( {i\widetilde{W}\left( {r,\omega }\right) ,\omega }\right) 。 $

(44)

$ z\left( t\right) = \sum \operatorname{res}\left( {i\widetilde{Z}\left( \omega \right) ,\omega }\right) 。 $

上述求和涉及被积函数的所有极点。由路面不平度产生的项是在 $\omega = \pm {\omega }_{\xi }$ 处,只考虑非负 $\omega$ ,则由路面不平度产生的振动谐波项为:

(45)

$ {w}_{1}\left( {r, t}\right) = - \frac{i{\omega }_{\xi }^{2}m{P}_{0}\left( {{c}_{t}i{\omega }_{\xi } + {k}_{t}}\right) K\left( {r,{\omega }_{\xi }}\right) }{2{\Delta }_{wz}\left( {\omega }_{\xi }\right) }{e}^{i{\omega }_{\xi }t} \circ $

(46)

$ {z}_{1}\left( t\right) = \frac{{i\pi }{P}_{0}\left( {i{c}_{t}{\omega }_{\xi } + {k}_{t}}\right) }{{\Delta }_{wz}\left( {\omega }_{\xi }\right) }{e}^{i{\omega }_{\xi }t}。 $

其中

(47)

$ {\Delta }_{wz}\left( {\omega }_{\xi }\right) = - K\left( {\omega }_{\xi }\right) \left( {{c}_{t}i{\omega }_{\xi } + {k}_{t}}\right) m{\omega }_{\xi }^{2} + {2\pi }\left( {{c}_{t}i{\omega }_{\xi } + {k}_{t} - }\right. \\ \left. {m{\omega }_{\xi }{}^{2}}\right) \text{。} $

振动的第二谐波项与计算的每个感应频率 ${\omega }_{i}$ 有关, ${\Delta }_{wz}\left( {\omega }_{i}\right) = 0$ ,无限长梁和车辆的响应可以表示为:

(48)

$ {w}_{2}\left( {r, t}\right) = \nabla \left\lbrack \left( {\omega - }\right. \right. \\ \left. {{\left. {\omega }_{i}\right) }^{n}\frac{i{\omega }_{\xi }{\omega }^{2}m{P}_{0}\left( {{c}_{t}{i\omega } + {k}_{t}}\right) K\left( {r,\omega }\right) }{{\Delta }_{wz}\left( \omega \right) \left( {{\omega }_{\xi }^{2} - {\omega }^{2}}\right) }{e}^{i\omega t}}\right\rbrack \circ \\ {z}_{2}\left( t\right) = \nabla \left\lbrack {{\left( \omega - {\omega }_{i}\right) }^{n}\frac{{i2\pi }{P}_{0}{\omega }_{\xi }\left( {i{c}_{t}\omega + {k}_{t}}\right) }{{\Delta }_{wz}\left( \omega \right) \left( {{\omega }_{\xi }^{2} - {\omega }^{2}}\right) }{e}^{i\omega t}}\right\rbrack \circ $

式中: $\nabla = \frac{1}{\left( {n - 1}\right) !}\mathop{\lim }\limits_{{\omega \rightarrow {\omega }_{i}}}\frac{{d}^{n - 1}}{{d}_{{\omega }^{n - 1}}};n$ 为 ${\omega }_{i}$ 的重根数。

4 仿真分析

收起

为了验证前文推导的半解析解, 以车辆与无限长路面振动耦合系统稳态响应求解为例, 将其与模态叠加法相比较研究。其中道路为沥青混凝土路面,车辆模型取自重型车辆。

耦合系统参数设置如下。

道路参数: 长度分别设 $L = {500}\text{、}{800}\text{、}{1000}\mathrm{\;m}$ 以及无限长,弹性模量 $E = {2.8} \times {10}^{10}\mathrm{\;N}/{\mathrm{m}}^{2}$ ,梁高 $h = {0.3}\mathrm{\;m}$ ,梁宽 $b = 6\mathrm{\;m}$ ,惯性矩 $I = {0.0135}{\mathrm{\;m}}^{4}$ , 沥青混凝土密度 $\rho = {2.5} \times {10}^{3}\mathrm{\;{kg}}/{\mathrm{m}}^{3}$ ,地基刚度 ${k}_{f} =$ ${4.8} \times {10}^{7}\mathrm{\;N}/\mathrm{m}$ 。

车辆参数: 车辆质量 $m = 1 \times {10}^{4}\mathrm{\;{kg}}$ ,刚度 ${k}_{t} =$ ${2.1} \times {10}^{6}\mathrm{\;N}/\mathrm{m}$ ,阻尼系数 ${c}_{t} = {700}\mathrm{\;N} \cdot \mathrm{s}/\mathrm{m}$ ,车速 $v =$ ${20}\mathrm{\;m}/\mathrm{s}$ 。

4.1 无限长路面与有限长路面车路耦合振动系统响应对比

将车辆与无限长路面振动耦合系统的仿真结果与车辆与有限长路面振动耦合系统的仿真结果进行对比分析,其中有限长路面选取 ${500}\mathrm{\;m}$ ,采用的模态数为 30 阶。考虑来自路面的激励为主要激励, 此时车辆与路面接触点处的路面位移、车辆垂向位移、车辆加速度的时域对比如图 2 所示。

图 2 中有限长路面的仿真结果是基于模态叠加法得到的, 无限长路面的仿真结果是基于上述提出的半解析解得到的。由图 2 可知, 本文半解析解与应用模态叠加法求得的一定路面长度和模态数下的振动耦合系统的数值解具有高度的一致性, 两者在峰值处存在偏差。经计算, 两种求解方法的平均偏差约为 ${3.6}\% \sim {4.8}\%$ ,偏差与模态数的大小和有限长路面的选取长度有关。

路面不平度通常表现为随机激励, 可用谐波叠加法进行复现, 即将其看成一系列具有不同频率和幅值的正弦波叠加。本文所涉及的模型为线性系统, 满足线性叠加原理, 因此, 上述提出的半解析解也可用于求解路面随机激励响应。此外, 半解析解法的求解效率明显高于模态叠加法, 其主要体现在仿真所耗时间上。分别考虑一定路面长度下不同模态数对求解耦合系统响应所耗时间的影响, 以及一定模态数下不同路面长度对求解耦合系统响应所耗时间的影响, 选取几组数据进行对比, 对比情况见表 1 。

由表 1 可知, 用模态叠加法计算所需时间随模态数和路面长度的增加而增加。而半解析解具有快速获得高精度结果的优点, 通过半解析解获取车路耦合振动系统响应所需时间约为 ${0.37}\mathrm{\;h}$ ,并且不需要测试结果的收敛性, 用于振动分析时其计算效率大大提高。

4.2 参数研究

为进一步研究上述偏差与模态数和路面长度之间的关系,分别取模态数 $N = {30}\text{、}{50}\text{、}{100}, L =$ ${500}\text{、}{800}\text{、}{1000}\mathrm{\;m}$ 进行仿真,仍考虑车辆与路面接触点处的路面位移、车辆垂向位移、车辆加速度, 得到时域对比, 如图 3 所示。

由图 3 a 可知,随着模态数的增加,应用模态叠加法求得的有限长路面与车辆振动耦合系统的路面位移响应随着模态数的增加而增大, 并逐渐接近于无限长路面与车辆振动耦合系统的振动响应, 两种求解方法的一致性增大。由图 3 b 和 c 可知,随着模态数的增加, 有限长路面与车辆振动耦合系统的车辆位移和加速度响应变化相对较小。

进一步研究道路长度的选取对结果的影响。由

5 结论

收起

1)本文考虑柔性路面特性的车路耦合振动系统, 并针对无限长道路与车辆耦合系统响应计算复杂难题, 提出通过积分变换推导耦合系统振动响应解析解的方法, 并应用留数定理对其进行数值计算, 获得车辆垂向位移、加速度、路面振动响应等半解析解。

2)将无限长路面与车辆耦合振动系统的半解析解与传统应用模态叠加法的有限长道路与车辆耦合振动系统响应相比,两者拥有高度一致性,证明了半解析解的有效性。

3)通过分析模态数和路面长度对有限长路面和车辆耦合系统振动响应的影响, 证明建立足够长的道路模型并采用合适的模态数才能使其计算更贴近无限长路面的响应。半解析解具有更高的计算效率与精度, 研究结果对于提高车辆的舒适性以及道路友好性优化与控制具有重要的理论指导意义。

基金

收起

国家自然科学基金项目(51705205)
中国博士后科学基金面上项目(2022M711384)
江苏省高等学校基础学科重大项目(21KJA460002)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

CIFUENTES

A O

. Dynamic Response of a Beam Excited by a Moving Mass[J]. Finite Elements in Analysis and Design, 1989,5:237-246.

[2]

SUN

. An Explicit Representation of Steady State Response of a Beam on an Elastic Foundation to Moving Harmonic Line Loads[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2003,27(1):69-84.

[3]

ZHANG

Yin

. Steady State Response of an Infinite Beam on a Viscoelastic Foundation with Moving Distributed Mass and Load[J]. Science China Physics, Mechanics & Astronomy, 2020,63(8):71-85.

[4]

DIMITROVOVÁ

. Semi-Analytical Solution for a Problem of a Uniformly Moving Oscillator on an Infinite Beam on a Two-Parameter Visco-Elastic Foundation[J]. Journal of Sound and Vibration, 2019,438:257-290.

[5]

ZHANG

Songhan

, FAN

Wei

, YANG

Caijin

. Semi-Analytical Solution to the Steady-State Periodic Dynamic Response of an Infinite Beam Carrying a Moving Vehicle[J]. International Journal of Mechanical Sciences, 2022,226:107409.1-107409.16.

[6]

FRÝBA

. Vibration of Solids and Structures Under Moving Loads[M]. Groningen: Noordhoff International Publishing, 2013.

[7]

RYS

, BURNOS

. Study on the Accuracy of Axle Load Spectra Used for Pavement Design[J]. International Journal of Pavement Engineering, 2021,23(11):3706-3715.

[8]

MARTÍNEZ

CASAS J

, MAZZOLA

, BAEZA

, et al. Numerical Estimation of Stresses in Railway Axles Using a Train-Track Interaction Model[J]. International Journal of Fatigue, 2013,47:18-30.

[9]

LIU

, MCCULLOUGH

B F

, OEY

H S

. Response of Rigid Pavements Due to Vehicle-Road Interaction[J]. Journal of Transportation Engineering, 2000,126(3):237-242.

[10]

KWASNIEWSKI

, LI

Hongyi

, WEKEZER

, et al. Finite Element Analysis of Vehicle-Bridge Interaction[J]. Finite Elements in Analysis and Design, 2006,42(11):950-959.

[11]

ZHAIT

, CAI

. Dynamic Interaction Between a Lumped Mass Vehicle and a Discretely Supported Continuous Rall Track[J]. Computers & Structures, 1997,63(5):987-997.

[12]

Y S

, YANG

Y B

. Steady-State Response and Riding Comfort of Trains Moving over a Series of Simply Supported Bridges[J]. Engineering Structures, 2003,25(2):251-265.

[13]

杨绍普, 陈立群, 李韶华. 车辆道路耦合系统动力学研究[M]. 北京: 科学出版社, 2012.

YANG

Shaopu

, CHEN

Liqun

, LI

Shaohua

. Dynamics of Vehicle-Road Coupled System[M]. Beijing: Science Press, 2012. (in Chinese)

[14]

MALLIK

A K

, CHANDRA

, SINGH

A B

. Steady-State Response of an Elastically Supported Infinite Beam to a Moving Load[J]. Journal of Sound and Vibration, 2006,291(3-5):1148-1169.

[15]

DING

, CHEN

Liqun

, YANG

Shaopu

. Convergence of Galerkin Truncation for Dynamic Response of Finite Beams on Nonlinear Foundations Under a Moving Load[J]. Journal of Sound and Vibration, 2012,331(10):2426-2442.

[16]

KRISHNANUNNI

, RAO

B N

. Decoupled Technique for Dynamic Response of Vehicle-Pavement Systems[J]. Engineering Structures, 2019,191:264-279.

[17]

KARGARNOVIN

M H

, YOUNESIAN

, THOMPSON

, et al. Response of Beams on Nonlinear Viscoelastic Foundations to Harmonic Moving Loads[J]. Computers & Structures, 2005,83(23-24):1865-1877.

[18]

DING

, YANG

Yan

, CHEN

Liqun

, et al. Vibration of Vehicle-Pavement Coupled System Based on a Timoshenko Beam on a Nonlinear Foundation[J]. Journal of Sound and Vibration, 2014,333(24):6623-6636.

[19]

JANG

T S

, SUNG

H G

. A New Semi-Analytical Method for the Non-Linear Static Analysis of an Infinite Beam on a Non-Linear Elastic Foundation: A General Approach to a Variable Beam Cross-Section[J]. International Journal of Non-Linear Mechanics, 2012,47(4):132-139.

[20]

Haitao

, YUAN

Yong

. Analytical Solution for an Infinite Euler-Bernoulli Beam on a Viscoelastic Foundation Subjected to Arbitrary Dynamic Loads[J]. Journal of Engineering Mechanics, 2014,140(3):542-551.

[21]

ZHEN

Bin

, XU

Jian

, SUN

Jianqiao

. Analytical Solutions for Steady State Responses of an Infinite Euler-Bernoulli Beam on a Nonlinear Viscoelastic Foundation Subjected to a Harmonic Moving Load[J]. Journal of Sound and Vibration, 2020,476:115271.1-115271.21.

[22]

FROIO

, RIZZI

, SIMÕES

F M F

, et al. Universal Analytical Solution of the Steady-State Response of an Infinite Beam on a Pasternak Elastic Foundation Under Moving Load[J]. International Journal of Solids and Structures, 2018,132-133:245-263.

2024年第14卷第2期

PDF下载

298

123

引用本文

BibTeX

文章信息

doi: 10.3969/j.issn.2095–1469.2024.02.06

接收时间：2023-02-05
首发时间：2025-07-21

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-02-05
修回日期：2023-03-14

基金

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

中国博士后科学基金面上项目(2022M711384)

江苏省高等学校基础学科重大项目(21KJA460002)

作者信息

¹ 江苏大学汽车工程研究院镇江 212013

² 江苏大学汽车与交通工程学院镇江 212013

³ 江苏恒力制动器制造有限公司靖江 214500

通讯作者:

张步云(1987-),男,江苏淮安人,博士,副教授,主要研究方向为振动测试与控制、车辆随机动力学分析。Tel: 13914560131 E-mail：zhangby@ujs.edu.cn

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/qcgcxb/CN/10.3969/j.issn.2095–1469.2024.02.06

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT