一种载流管道动力特性的计算方法与流程

文档序号:12364614阅读:385来源:国知局
一种载流管道动力特性的计算方法与流程

本发明涉及工程计算领域,尤其涉及一种载流管道动力特性的计算方法。



背景技术:

管道作为物料输送的一种特种设备在现代化工业生产和人民生活中起着很重要的作用,管道事故时有发生,严重影响着人民的生命和财产安全。管道振动对安全生产造成很大的威胁。强烈的管道振动会使管路附件,特别是管道的连接部位和管道附件的连接部位等处发生松动和破裂,造成严重事故。

现阶段,对管道内的流体流动会影响管道的横向弯曲动力特性——固有频率进行计算时,得出的计算结果显示出,随着流速的增大管道的横向振动固有频率降低,但是与实际的情况不相符,在试验中,一根弯曲的软管会随着水流流速的增大而伸直,即管道的固有频率随流速的增大而提高。



技术实现要素:

本发明要解决的技术问题是如何克服现有技术的不足,提供一种载流管道动力特性的计算方法。

本发明为实现上述目的采用的技术方案是:一种载流管道动力特性的计算方法,包括如下计算步骤:

(1)对管道微元段进行动力学分析,得出管道微元段的动力学方程,即公式1

<mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>t</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <mi>E</mi> <mi>I</mi> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>4</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>4</mn> </msup> </mrow> </mfrac> <mo>-</mo> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>-</mo> <mn>2</mn> <mi>&rho;</mi> <mi>A</mi> <mi>v</mi> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <mi>x</mi> <mo>&part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mn>0</mn> </mrow>

式中:

m——管道单位长度的质量;

y——管道的横向变形(挠度);

x——管道长度方向的坐标;

ρ——管道内流体密度;

A——管道内壁横截面积;

v——管道内流流速;

EI——管道横截面抗弯模量;

t——时间;

(2)对公式1进行求解,由于公式1是建立在小变形假定基础上的,截面转动的角速度是高阶小量,因此,忽略科氏加速度后得公知振型方程,即公式2:

φ(x)-αφ″(x)-βφ(x)=0

式中:

φ(x)——振型函数,

<mrow> <mi>&alpha;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> </mrow> <mrow> <mi>E</mi> <mi>I</mi> </mrow> </mfrac> <mo>,</mo> <mi>&beta;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>&omega;</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>E</mi> <mi>I</mi> </mrow> </mfrac> <mo>,</mo> </mrow>

ω——载流管道横向弯曲振动固有频率;

(3)对公式2进行求解,设φ(x)=epx,代入公式2得关于p的4次方程:

p4-αp2-β=0

令p2=s,则上式为s的一元二次方程

s2-αs-β=0

由上式求得

<mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mi>&alpha;</mi> <mo>+</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> <mo>,</mo> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mi>&alpha;</mi> <mo>-</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> </mrow>

则公式2的4个解为:

<mrow> <msub> <mi>p</mi> <mrow> <mn>1</mn> <mo>,</mo> <mn>2</mn> </mrow> </msub> <mo>=</mo> <mo>&PlusMinus;</mo> <msqrt> <mfrac> <mrow> <mi>&alpha;</mi> <mo>+</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> </msqrt> <mo>,</mo> <msub> <mi>p</mi> <mrow> <mn>3</mn> <mo>,</mo> <mn>4</mn> </mrow> </msub> <mo>=</mo> <mo>&PlusMinus;</mo> <mi>i</mi> <msqrt> <mfrac> <mrow> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> <mo>-</mo> <mi>&alpha;</mi> </mrow> <mn>2</mn> </mfrac> </msqrt> </mrow>

<mrow> <mi>&delta;</mi> <mo>=</mo> <msqrt> <mfrac> <mrow> <mi>&alpha;</mi> <mo>+</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> </msqrt> <mo>,</mo> <mi>&epsiv;</mi> <mo>=</mo> <msqrt> <mfrac> <mrow> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> <mo>-</mo> <mi>&alpha;</mi> </mrow> <mn>2</mn> </mfrac> </msqrt> </mrow>

由此得

φ(x)=C1eδx+C2e-δx+C3eiεx+C4e-iεx

由欧拉公式

<mrow> <mi>cosh</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mi>x</mi> </msup> <mo>+</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>x</mi> </mrow> </msup> </mrow> <mn>2</mn> </mfrac> <mo>,</mo> <mi>sinh</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mi>x</mi> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>x</mi> </mrow> </msup> </mrow> <mn>2</mn> </mfrac> </mrow>

<mrow> <mi>cos</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mrow> <mi>i</mi> <mi>x</mi> </mrow> </msup> <mo>+</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mi>x</mi> </mrow> </msup> </mrow> <mn>2</mn> </mfrac> <mo>,</mo> <mi>sin</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mrow> <mi>i</mi> <mi>x</mi> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mi>x</mi> </mrow> </msup> </mrow> <mrow> <mn>2</mn> <mi>i</mi> </mrow> </mfrac> </mrow>

得到振型函数,即公式3:

φ(x)=D1chδx+D2shδx+D3cosεx+D4sinεx

式中:

D1、D2、D3、D4——待定系数;

δ、ε——管道两端约束条件有关的参数,

管道两端的约束条件包括:两端简支、两端固定、一端固定、一端固定一端简支和一端固定一端自由。

例如:两端简支条件为挠度为零和弯矩为零,即φ(0)=φ(l)=0,φ″(0)=φ″(l)=0,将它们代入上式即可求出δ,ε,而D1、D2、D3、D4这些均为结构动力学理论中的公知技术在上述求解过程中同时得到;

(4)根据管道两端实际约束条件,计算得出ε参数值;

(5)根据管道两端的约束条件有公式3求得关于ω的公式,即公式4:

<mrow> <mi>&omega;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>&epsiv;</mi> <mn>2</mn> </msup> <mi>E</mi> <mi>I</mi> </mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> </mrow> </mfrac> <msqrt> <mrow> <mn>1</mn> <mo>+</mo> <mfrac> <mrow> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> </mrow> <mrow> <msup> <mi>&epsiv;</mi> <mn>2</mn> </msup> <mi>E</mi> <mi>I</mi> </mrow> </mfrac> </mrow> </msqrt> </mrow>

(6)针对公式4中载流管道的EI、m、A进行必要的测算,并且根据实际情况测算出参数ρ、v的值,代入公式4中得出载流管道横向弯曲振动固有频率,即ω值。

本发明正确的描述了载流管道内流流速对管道横向弯曲振动固有频率的影响,提出了正确的计算方法,降低了管道在使用中存在的安全风险。

附图说明

图1本发明试验方法的试验结果图示。

图2本发明管道微元段的受力和运动示意图,

其中,T——张力,N、dN——剪力和剪力增量,M、dM——弯矩和弯矩增量,

dθ——管道微元段截面转角,——管道微元段出口和出口处的流体动量,——管道微元段横向加速度,dx——管道微元段长度

具体实施方式

建立正确的载流管道动力特性计算方法,能够正确的反应载流管道内的流体流动对管道横向弯曲动力特性(固有频率)的影响。

目前的载流管道动力特性是基于下列方程计算的:

式中:

m——管道单位长度的质量;

y——管道的横向变形(挠度);

x——管道长度方向的坐标;

ρ——管道内流体密度;

A——管道内壁横截面积;

v——管道内流流速;

由于管道的振动属于小变形,其截面转动的角速度非常小,因此,科氏加速度项此时,方程1简化为:

<mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>t</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <mi>E</mi> <mi>I</mi> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>4</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>4</mn> </msup> </mrow> </mfrac> <mo>+</mo> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>=</mo> <mn>0</mn> </mrow>

并计算得出载流管道的固有频率计算公式:

<mrow> <mi>&omega;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>&epsiv;</mi> <mn>2</mn> </msup> <mi>E</mi> <mi>I</mi> </mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> </mrow> </mfrac> <msqrt> <mrow> <mn>1</mn> <mo>-</mo> <mfrac> <mrow> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> </mrow> <mrow> <msup> <mi>&epsiv;</mi> <mn>2</mn> </msup> <mi>E</mi> <mi>I</mi> </mrow> </mfrac> </mrow> </msqrt> </mrow>

式中:ε——与管道两端约束条件有关的参数,管道两端的约束条件为

简支边界条件时,

从上述的固有频率计算公式可以看出,随着流速的增大管道的横向振动固有频率降低,但这与自然现象不符,也与试验现象不符;一根弯曲的软管会随着水流流速增大而伸直,如消防水龙带。图1中是本发明试验方法的试验结果,其表征了载流管道的固有频率随流速的增大而提高。

综上,目前的载流管道横向弯曲振动固有频率的计算方法是不正确的,原因在“方程1”不能正确的描述事实现象。

载流管道发生弯曲振动时,由于管道的弯曲变形,使得管内流体的动量发生变化,对于直管,尽管流速变化引起动量的变化,但由于动量的方向与管道轴线平行,因此,不会使管道产生横向作用力。但对于弯曲的管道或管道发生横向变形时,即使流速为常数(定长流速),由于流速方向的变化也将引起流体动量发生变化,而这部分变化的动量其方向垂直于流速(流速方向变化引起的动量增量垂直于流速),即垂直于管道,从而产生的横向作用力。

由于动量的增量等于力的冲量,因此,这部分动量的增量是管道的张力冲量在动量增量方向的投影。图2中示出了管道微元段的受力和运动示意图,对其进行分析,由此可得到管道微元段的动力学方程,即公式1

<mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>t</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>+</mo> <mi>E</mi> <mi>I</mi> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>4</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>4</mn> </msup> </mrow> </mfrac> <mo>-</mo> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <msup> <mi>x</mi> <mn>2</mn> </msup> </mrow> </mfrac> <mo>-</mo> <mn>2</mn> <mi>&rho;</mi> <mi>A</mi> <mi>v</mi> <mfrac> <mrow> <msup> <mo>&part;</mo> <mn>2</mn> </msup> <mi>y</mi> </mrow> <mrow> <mo>&part;</mo> <mi>x</mi> <mo>&part;</mo> <mi>t</mi> </mrow> </mfrac> <mo>=</mo> <mn>0</mn> </mrow>

式中:

m——管道单位长度的质量;

y——管道的横向变形(挠度);

x——管道长度方向的坐标;

ρ——管道内流体密度;

A——管道内壁横截面积;

v——管道内流流速;

EI——管道横截面抗弯模量;

t——时间;

(2)对公式1进行求解,由于公式1是建立在小变形假定基础上的,截面转动的角速度是高阶小量,因此,忽略科氏加速度后得公知振型方程,即公式2:

φ(x)-αφ″(x)-βφ(x)=0

式中:

φ(x)——振型函数,

<mrow> <mi>&alpha;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> </mrow> <mrow> <mi>E</mi> <mi>I</mi> </mrow> </mfrac> <mo>,</mo> <mi>&beta;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>&omega;</mi> <mn>2</mn> </msup> <mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> </mrow> </mrow> <mrow> <mi>E</mi> <mi>I</mi> </mrow> </mfrac> <mo>,</mo> </mrow>

ω——载流管道横向弯曲振动固有频率;

(3)对公式2进行求解,设φ(x)=epx,代入公式2得关于p的4次方程:

p4-αp2-β=0

令p2=s,则上式为s的一元二次方程

s2-αs-β=0

由上式求得

<mrow> <msub> <mi>s</mi> <mn>1</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mi>&alpha;</mi> <mo>+</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> <mo>,</mo> <msub> <mi>s</mi> <mn>2</mn> </msub> <mo>=</mo> <mfrac> <mrow> <mi>&alpha;</mi> <mo>-</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> </mrow>

则公式2的4个解为:

<mrow> <msub> <mi>p</mi> <mrow> <mn>1</mn> <mo>,</mo> <mn>2</mn> </mrow> </msub> <mo>=</mo> <mo>&PlusMinus;</mo> <msqrt> <mfrac> <mrow> <mi>&alpha;</mi> <mo>+</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> </msqrt> <mo>,</mo> <msub> <mi>p</mi> <mrow> <mn>3</mn> <mo>,</mo> <mn>4</mn> </mrow> </msub> <mo>=</mo> <mo>&PlusMinus;</mo> <mi>i</mi> <msqrt> <mfrac> <mrow> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> <mo>-</mo> <mi>&alpha;</mi> </mrow> <mn>2</mn> </mfrac> </msqrt> </mrow>

<mrow> <mi>&delta;</mi> <mo>=</mo> <msqrt> <mfrac> <mrow> <mi>&alpha;</mi> <mo>+</mo> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> </mrow> <mn>2</mn> </mfrac> </msqrt> <mo>,</mo> <mi>&epsiv;</mi> <mo>=</mo> <msqrt> <mfrac> <mrow> <msqrt> <mrow> <msup> <mi>&alpha;</mi> <mn>2</mn> </msup> <mo>+</mo> <mn>4</mn> <mi>&beta;</mi> </mrow> </msqrt> <mo>-</mo> <mi>&alpha;</mi> </mrow> <mn>2</mn> </mfrac> </msqrt> </mrow>

由此得

φ(x)=C1eδx+C2e-δx+C3eiεx+C4e-iεx

由欧拉公式

<mrow> <mi>cosh</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mi>x</mi> </msup> <mo>+</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>x</mi> </mrow> </msup> </mrow> <mn>2</mn> </mfrac> <mo>,</mo> <mi>sinh</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mi>x</mi> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>x</mi> </mrow> </msup> </mrow> <mn>2</mn> </mfrac> </mrow>

<mrow> <mi>cos</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mrow> <mi>i</mi> <mi>x</mi> </mrow> </msup> <mo>+</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mi>x</mi> </mrow> </msup> </mrow> <mn>2</mn> </mfrac> <mo>,</mo> <mi>sin</mi> <mi> </mi> <mi>x</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>e</mi> <mrow> <mi>i</mi> <mi>x</mi> </mrow> </msup> <mo>-</mo> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mi>x</mi> </mrow> </msup> </mrow> <mrow> <mn>2</mn> <mi>i</mi> </mrow> </mfrac> </mrow>

得到振型函数,即公式3:

φ(x)=D1chδx+D2shδx+D3cosεx+D4sinεx

式中:

D1、D2、D3、D4——待定系数;

δ、ε——管道两端约束条件有关的参数,

管道两端的约束条件包括:两端简支、两端固定、一端固定、一端固定一端简支和一端固定一端自由。

其中一种约束条件:两端简支条件为挠度为零和弯矩为零,即φ(0)=φ(l)=0、φ″(0)=φ″(l)=0,将它们代入上式即可求出δ,ε,而D1、D2、D3、D4这些均为结构动力学理论中的公知技术在上述求解过程中同时得到;

(4)根据管道两端实际约束条件,计算得出ε参数值;

(5)根据管道两端的约束条件有公式3求得关于ω的公式,即公式4:

<mrow> <mi>&omega;</mi> <mo>=</mo> <mfrac> <mrow> <msup> <mi>&epsiv;</mi> <mn>2</mn> </msup> <mi>E</mi> <mi>I</mi> </mrow> <mrow> <mo>(</mo> <mi>m</mi> <mo>+</mo> <mi>&rho;</mi> <mi>A</mi> <mo>)</mo> </mrow> </mfrac> <msqrt> <mrow> <mn>1</mn> <mo>+</mo> <mfrac> <mrow> <msup> <mi>&rho;Av</mi> <mn>2</mn> </msup> </mrow> <mrow> <msup> <mi>&epsiv;</mi> <mn>2</mn> </msup> <mi>E</mi> <mi>I</mi> </mrow> </mfrac> </mrow> </msqrt> </mrow>

(6)针对公式4中载流管道的EI、m、A进行必要的测算,并且根据实际情况测算出参数ρ、v的值,代入公式4中得出载流管道横向弯曲振动固有频率,即ω值。

上述实施例只是为了说明本发明的技术构思及特点,其目的是在于让本领域内的普通技术人员能够了解本发明的内容并据以实施,并不能以此限制本发明的保护范围。凡是根据本发明内容的实质所作出的等效的变化或修饰,都应涵盖在本发明的保护范围内。

当前第1页1 2 3 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1