基于时程阵列交互迭代的车桥耦合振动分析方法及系统与流程

文档序号:32213495发布日期:2022-11-16 06:48阅读:31来源:国知局
基于时程阵列交互迭代的车桥耦合振动分析方法及系统与流程

1.本发明涉及车桥耦合振动分析技术领域,具体涉及基于时程阵列交互迭代 的车桥耦合振动分析方法及系统。


背景技术:

2.汽车或列车等车辆通过桥梁时,对桥梁产生冲击,引起桥梁的振动,桥梁 的振动又反过来对运行车辆的平稳性和安全性产生影响,车辆和桥梁组成一个 相互影响,相互关联的耦合振动统一系统,单独对车辆和桥梁进行研究,都不 能正确反映车辆和桥梁的振动响应。
3.目前,在车辆-桥梁耦合振动研究中主要是建立车辆桥梁整体动力方程, 同时求解车辆和桥梁的动力响应,这种方法很好地解决了车辆桥梁的耦合问题, 但需要在每一积分步内对整体动力方程质量、刚度和阻尼矩阵进行修正,对于 列车还要首先预测列车运行的“蛇行波”。另一种方法将车辆和桥梁作为两个 子系统,通过两个系统之间的协调关系采用迭代方法分别求得车辆和桥梁的动 力响应。这种方法针对汽车,建立车辆模型阻尼矩阵时需预先知道竖向力的大 小,而竖向力是通过迭代来求解,因此,在每个时间步,不仅要进行迭代,同 时每个迭代步还要进行阻尼矩阵的修正。针对列车求解轮轨力需要预先知道轮 对和钢轨的位移速度,而轮对和钢轨的位移速度本身是未知量,因此需要通过 先假定位移和速度值,再不断迭代,逐步收敛于真实解。
4.以上分析方法通过编制自己的逐步数值积分程序来具体实现。与任何数值 积分方法一样,逐步积分法的精度依赖于所取时间增量的长度,为针对时间步 的积分迭代。且以上分析方法为了可靠地进行分析,时间增量必须足够地短, 导致运算量都较大。


技术实现要素:

5.针对现有技术存在的不足,本发明提出基于时程阵列交互迭代的车桥耦合 振动分析方法及系统,提高了运算效率。
6.第一方面,本发明提供一种基于时程阵列交互迭代的车桥耦合振动分析方 法。
7.在第一种可实现方式中,基于时程阵列交互迭代的车桥耦合振动分析方法, 包括:基于ansys平台对桥梁子系统进行分析,获取桥梁初始响应时程阵列; ansys平台将桥梁初始响应时程阵列分享给matlab,matlab平台根据 桥梁初始响应时程阵列对车辆子系统进行分析,获取等效路面不平度时程阵列; 对桥梁子系统和车辆子系统进行耦合分析,根据等效路面不平度时程阵列获取 桥梁和车辆的振动响应。
8.结合第一种可实现方式,在第二种可实现方式中,获取桥梁初始响应时程 阵列,包括:在ansys平台基础上,根据预设的桥梁结构信息进行有限元建 模;确定车对桥的初始作用时程矩阵;在桥梁结构模型中施加车对桥的初始作 用时程矩阵的激励,获取桥梁初始响应时程阵列。
9.结合第二种可实现方式,在第三种可实现方式中,确定车对桥的初始作用 时程矩
阵,包括:根据车辆运动方程确定桥对车的初始作用时程矩阵;根据桥 对车的初始作用时程矩阵确定车对桥的初始作用时程矩阵。
10.结合第二种可实现方式,在第四种可实现方式中,在桥梁结构模型中施加 车对桥的初始作用时程矩阵的激励,获取桥梁初始响应时程阵列,包括:采用 插值法依次提取车辆前后两轴与桥梁主梁接触点的响应、以及四个桥梁主梁相 关节点的响应;将车辆前后两轴与桥梁主梁接触点的响应以及四个桥梁主梁相 关节点的响应组成桥梁初始响应时程阵列。
11.结合第一种可实现方式,在第五种可实现方式中,根据ansys平台分享 的桥梁初始响应时程阵列获取等效路面不平度时程阵列,包括:在matlab 平台上,根据桥梁初始响应时程阵列和预设的模拟路面不平度时程,获得等效 路面不平度时程阵列。
12.结合第一种可实现方式,在第六种可实现方式中,根据等效路面不平度时 程阵列获取桥梁和车辆的振动响应,包括:根据等效路面不平度时程阵列获取 桥梁振动响应时程阵列;判断桥梁振动响应时程阵列是否收敛;在桥梁振动响 应时程阵列不收敛的情况下,根据桥梁振动响应时程阵列获取更新后的等效路 面不平度时程阵列,将更新后的等效路面不平度时程阵列作为等效路面不平度 时程阵列重复上述过程;在桥梁振动响应时程阵列收敛的情况下,获取车辆振 动响应时程阵列。
13.结合第六种可实现方式,在第七种可实现方式中,根据等效路面不平度时 程阵列获取桥梁振动响应时程阵列,包括:根据等效路面不平度时程阵列获取 车桥接触点处的桥对车的作用时程矩阵;根据桥对车的作用时程阵列确定车对 桥的作用时程阵列;根据车对桥的作用时程阵列获取桥梁振动响应时程阵列。
14.结合第六种可实现方式,在第八种可实现方式中,根据等效路面不平度时 程阵列获取车桥接触点处的桥对车的作用时程矩阵,包括:
15.通过计算获得车桥接触点处的桥对车 的作用时程矩阵;
16.其中,f
vbzi
为车轮承受的竖向作用时程矩阵,f
vbhi
为车轮承受的横向作用时 程矩阵,c
lzi
为轮胎与路面接触点处的竖向阻尼系数,k
lzi
为轮胎与路面接触点 处的刚度系数;为车轮与路面接触点处路面的竖向速度,z
ci
为车轮与路面 接触点处路面的竖向位移;为车轮与路面接触点处车轮的竖向速度,z
si
为车 轮与路面接触点处车轮的竖向位移;c
hi
为轮胎与路面接触点处的横向阻尼系数, k
hi
为轮胎与路面接触点处的弹性系数;为车轮与路面接触点处路面的横向速 度,y
ci
为车轮与路面接触点处路面的横向位移;为车轮与路面接触点处车轮 的横向速度,y
si
为车轮与路面接触点处车轮的横向位移。
17.结合第七种可实现方式,在第九种可实现方式中,在桥梁振动响应时程阵 列收敛的情况下,获取车辆振动响应时程阵列,包括:
18.在桥梁振动响应时程阵列收敛的情况下,根据桥对车的作用时程阵列获取 车辆振动响应时程阵列。
19.第二方面,本发明提供一种基于时程阵列交互迭代的车桥耦合振动分析系 统。
20.在第十种可实现方式中,基于时程阵列交互迭代的车桥耦合振动分析系统, 包括:桥梁子系统分析模块,被配置为基于ansys平台对桥梁子系统进行分 析,获取桥梁初始响应阵列;车辆子系统分析模块,被配置为将ansys平台 与matlab平台进行数据分享,matlab平台获得桥梁初始响应阵列,并对 车辆子系统进行分析,根据桥梁初始响应阵列获取等效路面不平度时程阵列; 耦合分析模块,被配置为对桥梁子系统和车辆子系统进行耦合分析,根据等效 路面不平度时程阵列获取桥梁和车辆的振动响应。
21.由上述技术方案可知,本发明的有益技术效果如下:
22.1.结合并充分利用matlab平台高效的矩阵计算功能和ansys结构分 析软件专业的建模、有限元计算、前后处理功能,将车、桥两者作为一个相互 作用而又统一的系统,利用matlab和ansys之间的数据共享,分别对车 辆子系统、桥梁子系统进行分析,然后对车辆子系统、桥梁子系统进行耦合分 析,最终得到桥梁和车辆的振动响应。这样,避免了对时间步的积分迭代,减 少了运算量,提高运算效率。
23.2.利用matlab强大的矩阵计算功能,以时程阵列的方式进行数据交换 更新,只需要在一开始根据车辆过桥的总时间设置好时间步,避免了逐步积分 法对时间增量的稳定敏感性,运算效率大大提高;同样地,时程阵列方式对于 复杂车辆、多车辆无非是矩阵列、阵列数的增加,建模可轻易实现,具有普遍 适用性和拓展性。
附图说明
24.为了更清楚地说明本发明具体实施方式或现有技术中的技术方案,下面将 对具体实施方式或现有技术描述中所需要使用的附图作简单介绍。在所有附图 中,类似的元件或部分一般由类似的附图标记标识。附图中,各元件或部分并 不一定按照实际的比例绘制。
25.图1为本发明提供的一种基于时程阵列交互迭代的车桥耦合振动分析方 法的示意图;
26.图2为本发明提供的一种获取桥梁和车辆的振动响应的流程图;
27.图3为本发明提供的一种基于时程阵列交互迭代的车桥耦合振动分析系 统的结构示意图;
28.图4为本发明提供的一种匀速移动常力通过简支梁模型示意图;
29.图5-a为本发明提供的一种应用krylov的理论解计算的跨中位移的曲线图;
30.图5-b为本发明提供的一种应用文献解计算的跨中位移的曲线图;
31.图5-c为本发明提供的一种采用本技术的分析方法计算的跨中位移曲线图。
具体实施方式
32.下面将结合附图对本发明技术方案的实施例进行详细的描述。以下实施例 仅用于更加清楚地说明本发明的技术方案,因此只作为示例,而不能以此来限 制本发明的保护范围。
33.需要注意的是,除非另有说明,本技术使用的技术术语或者科学术语应当 为本发明所属领域技术人员所理解的通常意义。
34.在一些实施例中,预先设置车辆信息和桥梁结构相关参数,并形成相应数 据文件
以便平台调用,预设设置公路等级信息,根据公路等级信息随时模拟生 成路面不平度,并形成数据文件以便平台调用。
35.结合图1所示,本实施例提供了一种基于时程阵列交互迭代的车桥耦合振 动分析方法,包括:
36.步骤s01、基于ansys平台对桥梁子系统进行分析,获取桥梁初始响应 阵列;
37.步骤s02、ansys平台与matlab平台进行数据分享后,matlab平 台对车辆子系统进行分析,根据ansys平台分享的桥梁初始响应时程阵列获 取等效路面不平度时程阵列;
38.步骤s03、对桥梁子系统和车辆子系统进行耦合分析,根据等效路面不平 度时程阵列获取桥梁和车辆的振动响应。
39.结合并充分利用matlab(matrix&laboratory,矩阵工厂)平台高效的矩 阵计算功能和ansys(analysis of systems,系统分析)结构分析软件专业的 建模、有限元计算、前后处理功能,将车、桥两者作为一个相互作用而又统一 的系统,利用matlab和ansys之间的数据共享,分别对车辆子系统、桥 梁子系统进行分析,然后对车辆子系统、桥梁子系统进行耦合分析,最终得到 桥梁和车辆的振动响应。避免对时间步的积分迭代,提高了运算效率。
40.可选地,获取桥梁初始响应阵列,包括:在ansys平台基础上,根据预 设的桥梁结构信息进行有限元建模;车对桥的初始作用时程矩阵;在桥梁结构 模型中施加车对桥的初始作用时程矩阵的激励,获取桥梁初始响应时程阵列。
41.在一些实施例中,在ansys平台上,对预设的桥梁结构相关参数对应的 数据文件进行调用,进行有限云建模,获得桥梁结构模型。
42.可选地,确定车对桥的初始作用时程矩阵,包括:根据车辆运动方程确定 桥对车的初始作用时程矩阵;根据桥对车的初始作用时程矩阵确定车对桥的初 始作用时程矩阵。
43.可选地,车辆运动方程通过以下公式表示:
[0044][0045]
在上式中,mv为车辆的总体质量,cv为阻尼,kv为刚度矩阵,为车 辆的加速度,为车辆的速度,uv为车辆的位移,f
vb
为桥对车的作用力时程 矩阵。
[0046]
在对车辆和桥梁进行统一分析的时候,车辆运动方程的质量、刚度矩阵、 阻尼都是随着汽车在桥上的运行而变化。采用统一的耦合运动方程,则分析中 每一时步都需对系数矩阵进行重新生成与分解;同时,随着上桥车辆数量的增 加,耦联的自由度也越来越多,使得计算工作量增大,编程也比较困难。本申 请以车轮与桥面接触点为界,将车与桥梁作为两个子系统,分别对车辆子系统 和桥梁子系统进行分析,再根据两个子系统间的耦合关系进行平衡迭代。这样, 车辆及桥梁各自运动方程的系数矩阵在各时步保持不变,大大降低了运算难度。 而且,利用matlab强大的矩阵计算功能,以时程阵列的方式进行数据交换 更新,只需要在一开始根据车辆过桥的总时间设置好时间步,避免了逐步积分 法对时间增量的稳定敏感性,运算效率进一步提高;同样地,时程阵列方式对 于复杂车辆、多车辆无非是矩阵列、阵列数的增加,建模可轻易实现,具有普 遍适用性和拓展性。
[0047]
可选地,根据相互作用原理,桥对车的作用力与车对桥的作用力大小相同, 方向相反,根据桥对车的初始作用时程矩阵f
vb
即可获得车对桥的初始作用时 程矩阵。
[0048]
在一些实施例中,车辆的初始加速度初始速度和初始位移均 为0,
将初始加速度、初始速度、初始位移代入车辆运动方程,获得桥对车的 初始作用时程矩阵为0,根据相互作用原理,则车对桥的初始作用时程矩阵f
bv
0(t) 为0。
[0049]
可选地,在桥梁结构模型中施加车对桥的初始作用时程矩阵的激励,获取 桥梁初始响应阵列,包括:采用插值法依次提取车辆前后两轴与桥梁主梁接触 点的响应、以及四个桥梁主梁相关节点的响应;将车辆前后两轴与桥梁主梁接 触点的响应以及四个桥梁主梁相关节点的响应组成桥梁初始响应阵列。这样, 可以考虑到多车辆以不同车速、间距等因素上、下桥行进作用的整个过程,实 现车桥接触点的几何、力学关系传递,提高了分析的准确性以及适用性。
[0050]
在一些实施例中,获得车对桥的初始作用时程矩阵后,运用ansys瞬态 分析计算获得桥梁初始响应阵列,在结果提取时,采用插值法依次提取车辆前 后两轴与桥梁主梁接触点的响应、以及四个桥梁主梁相关节点的响应,将车辆 前后两轴与桥梁主梁接触点的响应以及四个桥梁主梁相关节点的响应组成桥 梁初始响应阵列。
[0051]
可选地,获取等效路面不平度时程阵列,包括:在matlab平台上,根 据桥梁初始响应阵列和模拟路面不平度时程,获得等效路面不平度时程阵列。
[0052]
在一些实施例中,在matlab平台上,将桥梁初始响应阵列组合模拟路 面不平度时程后,形成等效路面不平度时程阵列。等效路面不平度时程阵列表示为其中,为桥梁竖向变形和路 面不平度的综合响应,为车轮与路面接触点处路面的竖向加速度,为车 轮与路面接触点处路面的竖向速度,z
ci
为车轮与路面接触点处路面的竖向位移;为是桥梁横向变形和路面不平度的综合响应,为车轮与路面接 触点处路面的横向加速度,为车轮与路面接触点处路面的横向速度,y
ci
为车 轮与路面接触点处路面的横向位移。
[0053]
可选地,根据等效路面不平度时程阵列获取桥梁和车辆的振动响应,包括: 根据等效路面不平度时程阵列获取桥梁振动响应时程阵列;判断桥梁振动响应 时程阵列是否收敛;在桥梁振动响应时程阵列不收敛的情况下,根据桥梁振动 响应时程阵列获取更新后的等效路面不平度时程阵列,将更新后的等效路面不 平度时程阵列作为新的等效路面不平度时程阵列重复上述过程;在桥梁振动响 应时程阵列收敛的情况下,获取车辆振动响应时程阵列。
[0054]
可选地,根据等效路面不平度时程阵列获取桥梁振动响应时程阵列,包括: 根据等效路面不平度时程阵列获取车桥接触点处的桥对车的作用时程矩阵;根 据桥对车的作用时程阵列确定车对桥的作用时程阵列;根据车对桥的作用时程 阵列获取桥梁振动响应时程阵列。
[0055]
可选地,通过以下公式获得车桥接触点处的桥对车的作用时程矩阵:
[0056][0057][0058]
上式中,f
vbzi
为车轮承受的竖向作用时程矩阵,f
vbhi
为车轮承受的横向作用 时程矩阵,c
lzi
为轮胎与路面接触点处的竖向阻尼系数,k
lzi
为轮胎与路面接触 点处的刚度系数;为车轮与路面接触点处路面的竖向速度,z
ci
为车轮与路 面接触点处路面的竖向位
移;为车轮与路面接触点处车轮的竖向速度,z
si
为 车轮与路面接触点处车轮的竖向位移;c
hi
为轮胎与路面接触点处的横向阻尼系 数,k
hi
为轮胎与路面接触点处的弹性系数;为车轮与路面接触点处路面的横 向速度,y
ci
为车轮与路面接触点处路面的横向位移;为车轮与路面接触点处 车轮的横向速度,y
si
为车轮与路面接触点处车轮的横向位移。
[0059]
在一些实施例中,根据相互作用原理,通过桥对车的作用时程阵列f
vb
确 定出车对桥的作用时程阵列f
bv
。在ansys平台,在桥梁结构模型中施加车对 桥的作用时程阵列f
bv
的激励,获取桥梁振动响应阵列其中,为第n次更新的桥梁结构的加速度、为第n次更新的桥梁结构的速度、为第n次更新的桥梁结构的位移。ansys平台将桥梁振动响应阵列分享给 matlab平台。matlab平台根据桥梁振动响应时程阵列对车辆子系统进行 分析,获取更新后的等效路面不平度时程阵列。
[0060]
可选地,在桥梁振动响应时程阵列收敛的情况下,获取车辆振动响应时程 阵列,包括:在桥梁振动响应时程阵列收敛的情况下,根据桥对车的作用时程 阵列获取车辆振动响应时程阵列。
[0061]
在一些实施例中,已知桥对车的作用时程阵列f
vb
,运用matlab平台的 runge-kutta法求解车辆运动方程式获得车辆振动响应时程 阵列其中,为第n次更新的车辆的加速度、为第n 次更新的车辆的速度、为第n次更新的车辆的位移。
[0062]
结合图2所示,在一些实施例中,获取桥梁和车辆的振动响应的流程包括:
[0063]
步骤s11、在ansys平台基础上,根据预设的桥梁结构信息进行有限元 建模,确定车对桥的初始作用时程矩阵,在桥梁结构模型中施加车对桥的初始 作用时程矩阵的激励,获取桥梁初始响应时程阵列;
[0064]
步骤s12、在matlab平台上,根据桥梁初始响应时程阵列和预设的模 拟路面不平度时程,获得等效路面不平度时程阵列;
[0065]
步骤s13、根据等效路面不平度时程阵列获取车桥接触点处的桥对车的作 用时程矩阵;
[0066]
步骤s14、根据桥对车的作用时程阵列确定车对桥的作用时程阵列;
[0067]
步骤s15、根据车对桥的作用时程阵列获取桥梁振动响应时程阵列;
[0068]
步骤s16、判断桥梁振动响应时程阵列是否收敛;在桥梁振动响应时程阵 列不收敛的情况下,执行步骤s17;在桥梁振动响应时程阵列收敛的情况下, 执行步骤s18;
[0069]
步骤s17、根据桥梁振动响应时程阵列获取更新后的等效路面不平度时程 阵列,将更新后的等效路面不平度时程阵列作为等效路面不平度时程阵列返回 步骤s13;
[0070]
步骤s18、根据桥对车的作用时程阵列获取车辆振动响应时程阵列。
[0071]
这样,以时程阵列的方式进行数据交换更新,只需要在一开始根据车辆过 桥的总时间设置好时间步,一般设为1/500即可满足精度要求,避免了逐步积 分法对时间增量的稳定敏感性,运算效率大大提高;而且,时程阵列方式对于 复杂车辆、多车辆无非是矩阵列、阵列数的增加,建模可轻易实现,具有普遍 适用性和拓展性。
[0072]
结合图3所示,本实施例提供一种基于时程阵列交互迭代的车桥耦合振动 分析系统,包括桥梁子系统分析模块101、车辆子系统分析模块102和耦合分 析模块103。桥梁子系
统分析模块101被配置为基于ansys平台对桥梁子系 统进行分析,获取桥梁初始响应阵列;车辆子系统分析模块102被配置为 ansys平台将桥梁初始响应阵列分享给matlab平台,matlab平台根据 桥梁初始响应阵列对车辆子系统进行分析,获取等效路面不平度时程阵列;耦 合分析模块103被配置为对桥梁子系统和车辆子系统进行耦合分析,根据等效 路面不平度时程阵列获取桥梁和车辆的振动响应。
[0073]
利用matlab平台高效的矩阵计算功能和ansys结构分析软件专业的 建模、有限元计算、前后处理功能,将车、桥两者作为一个相互作用而又统一 的系统,分别建立车、桥子系统的运动方程,通过车桥之间的力学平衡与几何 相容条件实现数据互通及迭代,对各子系统交互求解,最终得到协调一致的车 桥耦合振动分析结果。这样,不需要进行全域分析,从而避免了对庞大的系统 矩阵的运算,进而提高了运算效率,通过时程阵列的方式实现数据交换更新, 进一步提高了运算效率。避免了针对时间步的积分迭代,不需要反复调整收敛 参数,进而提高了收敛的速度。同时,还可以模拟更多复杂的车体,适用性更 强。
[0074]
充分利用并结合现有软件平台的技术优势,在其基础上开发自己所需的算 法功能,具有广泛的应用前景,相比而言,自编程序较复杂,无法或很难建立 复杂的桥梁模型,应用受到限制。ansys前处理模块提供了一个强大的实体 建模及网格划分工具,其丰富的单元库可支持大多数结构根据需要构建精确的 有限元模型。其次,ansys还专门设有与众所熟知的一些大型通用分析软件 的数据接口,方便数据的交换。此外,ansys后处理器可对分析结果进行各 种便利操作。以上优点使其非常适合作为自主程序的运行平台。
[0075]
在一些实施例中,对本技术的基于时程阵列交互迭代的车桥耦合振动分析 方法的验证过程如下:
[0076]
桥梁分析模型为如图4所示的跨度l=1.1938m的简支梁,截面积 a=0.51e-2m2(等截面),惯性矩i=0.9448e-8m4,材料弹性模量e=10.4804e10 n/m2,材料ρ=2.9602e3 kg/m3。桥梁有限元模型由100个beam4单元模拟,不 考虑桥梁阻尼的影响。
[0077]
可选地,等截面简支梁的自振频率理论计算式为:
[0078][0079]
将桥梁定义参数代入上式即可得到桥梁分析模型的理论自振频率。
[0080]
表1为简支梁前10阶振动频率(hz)
[0081][0082]
在一些实施例中,采用本技术的分析方法计算桥梁前100阶模态,表1 列出了其中前10阶模态的频率。此外,(杨建荣,2008)用其自编程序所求得的 解也列于表1中。如表1所示,序号1对应的理论解为8.93,对应的文献解为 8.93,采用本分析方法计算的本解为8.93;序号5对应的理论解为223.16,对 应的文献解为223.12,采用本分析方法计算的本解为223.12;序号10对应的 理论解为892.63,对应的文献解为892.04,采用本分析方法计算
的本解为 892.06。对比桥梁模型的前10阶振动频率,采用本分析方法的计算结果与理 论频率几乎完全相同。
[0083]
可选地,定义速度参数α为α=t/τ,其中,t为桥梁的第1阶振动周期,τ 表示车辆过桥的耗时。由此,速度参数也可表示为α=tv/l,其中l为桥梁跨 径,v为车速。当移动常力f=52.9442n以不同的车速过桥时,采用本分析方法 计算桥梁的动力响应。
[0084]
图5-a为应用krylov的理论解计算跨中位移曲线图,图5-b为应用文献解 (杨建荣,2008)计算跨中位移曲线图,图5-c为采用本分析方法计算跨中位移曲 线图,三种方法计算结果几乎完全一致。以上计算结果表明,采用本技术的分 析方法在计算移动常力作用下的桥梁动力响应结果正确。
[0085]
通过算例分析表明,本分析方法仅经一次交互迭代即可获得较高精度及可 靠的数值结果,大大提高了运算速度。
[0086]
最后应说明的是:以上实施例仅用以说明本发明的技术方案,而非对其限 制;尽管参照前述各实施例对本发明进行了详细的说明,本领域的普通技术人 员应当理解:其依然可以对前述各实施例所记载的技术方案进行修改,或者对 其中部分或者全部技术特征进行等同替换;而这些修改或者替换,并不使相应 技术方案的本质脱离本发明各实施例技术方案的范围,其均应涵盖在本发明的 权利要求和说明书的范围当中。
当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1