内外流耦合作用下的开采立管涡激振动特性分析方法



1.本发明涉及油气开发领域,具体涉及外部海流与内部多相流耦合作用下的开采立管涡激振动特性分析方法。


背景技术:



2.深水油气立管是油气开发系统中重要的部件,是连接海底和开采平台的重要结构,立管在来流的作用下产生交替泄放的尾涡,从而导致涡激振动现象的发生,若长时间振动,立管强度会降低,易发生失效事故,因此海洋立管的涡激振动响应得到学者的持续关注,对海洋立管在外部环境载荷作用下产生的涡激振动响应已有大量的研究,特别是针对大长细比的柔性立管。这些研究主要分析了涡激振动的模态、频率以及驻波和行波响应等,为深入理解涡激振动机理和响应规律奠定了基础。
3.与传统含纯液流立管相比,管内气液两相流会引起立管的固有频率降低,振动幅度和振动频率增加;管内排量、两相混合密度的增加,会引起立管的固有频率和振动响应频率同时降低,其模态响应保持不变;而随着进气比的增加,同样会引起立管固有频率降低,但振动频率增大,进而诱发立管产生更高阶模态的振动。
4.目前大多数立管的涡激振动响应研究只考虑外部海洋环境条件或管内流体单一方面的影响。考虑内外流耦合作用下的立管涡激振动的研究甚少,立管在油气开采过程中内部多相流动态相变与外部海洋环境载荷参数变化耦合作用下立管涡激振动响应特性与机理认识还不清晰。


技术实现要素:



5.本发明的目的在于提供内外流耦合作用下的开采立管涡激振动特性分析方法,为多相流开采立管的安全性提供有效的理论指导。
6.为实现上述目的,内外流耦合作用下的开采立管涡激振动特性分析方法,提出了考虑内外流耦合作用下的开采立管数值分析模型和管内多相流模型分析方法,包括开采立管力学模型分析方法、流体力模型分析方法、尾流振子模型分析方法和管内两相流模型分析方法,采用有限元法和newmark-β法对模型进行求解,具体求解步骤如下所示:
7.s1:根据外部海流速度大小、立管长度、立管内外径、内流速度等参数对深海含两相流立管进行模拟,模拟开采立管涡激振动过程;
8.s2:分析开采立管单元和管内流体单元段受力,建立考虑内外流耦合作用下的开采立管力学模型,此基础上,建立考虑管内流型对立管动力响应影响的两相流模型,计算管内流体质量和速度分布;
9.s3:将计算出的管内流体质量和速度分布情况与开采立管模型结合,基于模型的边界条件采用newmark-β法,求解考虑内外流耦合作用下的开采立管力学模型分析方法;
10.s4:根据计算出的考虑内外流耦合作用下的立管动力响应,模拟分析立管的涡激振动强弱以及振动频率的大小,总结了内外流耦合作用下的开采立管涡激振动响应规律。
11.进一步的,s1中的两相流立管管内的两相流能够对开采立管产生动态激励作用,在开采立管运移两相流的时候,进一步将开采立管简化为两端简支的内输多相流立管;立管置于横向流动的海流中产生大变形弯曲,同时由于海水绕过圆柱管道会因旋涡脱落而产生涡激力作用于开采立管上,进而使得立管产生涡激振动;
12.进一步的,s2中在分析相流动基本方程推导过程中,将两相分别按单相流处理并计入相间的相互作用,再将各相的方程加以合并,具体步骤为:
13.根据流体微元单元段受力分析图,可得管内多相流流体微元段x方向上的力平衡方程:
[0014][0015]
开采立管微元段x方向的力平衡方程:
[0016][0017]
剪切力q与弯矩m可以表示为:
[0018][0019][0020]
联立方程(1)、(2)、(3)和(4)可以得到两相流的深水立管顺流向(x方向)运动微分方程为:
[0021][0022]
同理可得到开采立管横流方向(y方向)上的运动微分方程:
[0023][0024]
式中:ei—弯曲刚度(n
·
m2);x—顺流向位移(m);y—横流向位移(m);z—立管轴向位置(m);t—截面张力(n);c—结构阻尼;f
x
(z,t)—顺流向海洋环境载荷(n);
[0025]fy
(z,t)—横流向海洋环境载荷(n);v
l
、vg—气、液体的速度(m/s);mr、m
l
与mg—单位长度立管质量、单位长度立管内液体与气体质量(kg);
[0026]
在轴向方向上,开采立管轴向上需要考虑立管自身重力与管内流体的重力,因此立管截面有效张力沿程是变化的。开采立管除了受到顶部张力与自身浮力外,其上部立管还要承受下部立管的重力,因此截面上的张力t=t(z),表示为:
[0027][0028]
式中:t
top
—顶部张力(n);ρw—海水密度(kg/m3);d0—立管外径(m)。
[0029]
进一步的,s3中的边界条件为使用newmark-β法计算开采立管在全局坐标系下振动控制矩阵方程
[0030]
其中[m]、[c]、[k]分别为质量矩阵,阻尼矩阵与刚度矩阵;分别为立管加速度、速度与位移矢量;{f}为外界的流体力载荷。
[0031]
进一步的,s4中在开采立管振动求解过程中,保持运算时间步一致,将尾流振子模型与立管控制方程耦合迭代计算,可以分别得到开采立管顺流向与横流向的位移时程响应,通过对位移时程曲线进行快速傅里变换,可以得到开采立管顺流向与横流向的频率响应。
[0032]
进一步的,s4中使用微分控制方程计算悬挂立管的无阻尼自由振动,得出模型的固有频率,对气液两相流和纯液流下的涡激振动特性对比分析、不同排量下的涡激振动特性对比分析、不同两相混合密度下的涡激振动特性对比分析和不同进气比下的涡激振动特性对比分析进行总结。
[0033]
基于上述分析方法,过建立内外流耦合作用下的立管系统动力学分析模型,分析含气液两相流的立管的涡激振动特性。研究了排量、混合流体密度、进气比等因素对涡激振动的影响,本发明具有以下结论:
[0034]
(1)本发明提供的内外流耦合作用下的开采立管涡激振动特性分析方法,可以得出,液两相流会引起立管的固有频率降低,在相同的外部激励下,导致其振动幅值更大,振动频率也更大,且表现出的模态响应也更高阶。因此,在实际工程问题中,应重点关注管内气液两相流对立管涡激振动的影响。
[0035]
(2)本发明提供的内外流耦合作用下的开采立管涡激振动特性分析方法,可以得出,随着管内流体排量和流体混合密度的增加,立管轴向张力减小,固有频率降低,振动幅值增大。由于此时发生频率锁定现象,振动频率呈现与固有频率同步降低的趋势,模态振型保持不变。
[0036]
(3)本发明提供的内外流耦合作用下的开采立管涡激振动特性分析方法,可以得出,进气比的增大导致立管刚度减小,振动幅值增加。随着管内含气量的增加,管内两相流会诱发立管更高一阶模态的振动,振动频率和模态振型均增大。
附图说明
[0037]
图1是本发明模型求解流程示意图;
[0038]
图2是本发明物理简化模型示意图示意图;
[0039]
图3是本发明立管微元段受力图;
[0040]
图4是本发明流体微元段受力图;
[0041]
图5是本发明海流下开采立管某截面某时刻的运动示意图;
[0042]
图6是本发明开采立管外部的流体力示意图;
[0043]
图7是本发明流体单元质量守恒的物理模型示意图;
[0044]
图8是本发明气液两相流和纯液流下立管中点处运动轨迹对比图;
[0045]
图9是本发明不同排量下立管中点处运动轨迹对比图;
[0046]
图10是本发明不同两相混合密度下立管中点处运动轨迹对比图;
[0047]
图11是本发明不同进气比下立管中点处运动轨迹对比图;
具体实施方式
[0048]
为了更好地了解本发明的目的、结构及功能,下面结合附图,对本发明的内外流耦合作用下的开采立管涡激振动特性分析方法做进一步详细的描述。
[0049]
在本发明的描述中,需要说明的是,术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系,仅是为了便于描述本发明和简化描述,而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作,因此不能理解为对本发明的限制。此外,术语“第一”、“第二”、“第三”仅用于描述目的,而不能理解为指示或暗示相对重要性。
[0050]
如图1-4所示,采立管在海洋环境中除受自重外,开采立管还受到外部海流的和内部多相流的作用。外部海流会在立管表面交替泄放漩涡从而产生涡激振动;内部多相流的质量、密度等参数随时间发生变化,管内两相流能够对开采立管产生动态激励作用,有开采立管力学模型分析方法、流体力模型分析方法、尾流振子模型分析方法和管内两相流模型分析方法,采用有限元法和newmark-β法对模型进行求解,具体求解步骤如下所示:
[0051]
s1:根据外部海流速度大小、立管长度、立管内外径、内流速度等参数对深海含两相流立管进行模拟,模拟开采立管涡激振动过程;
[0052]
s2:分析开采立管单元和管内流体单元段受力,建立考虑内外流耦合作用下的开采立管力学模型,此基础上,建立考虑管内流型对立管动力响应影响的两相流模型,计算管内流体质量和速度分布;
[0053]
s3:将计算出的管内流体质量和速度分布情况与开采立管模型结合,基于模型的边界条件采用newmark-β法,求解考虑内外流耦合作用下的开采立管力学模型分析方法;
[0054]
s4:根据计算出的考虑内外流耦合作用下的立管动力响应,模拟分析立管的涡激振动强弱以及振动频率的大小,总结了内外流耦合作用下的开采立管涡激振动响应规律。
[0055]
s1中的两相流立管管内的两相流能够对开采立管产生动态激励作用,在开采立管运移两相流的时候,进一步将开采立管简化为两端简支的内输多相流立管;立管置于横向流动的海流中产生大变形弯曲,同时由于海水绕过圆柱管道会因旋涡脱落而产生涡激力作用于开采立管上,进而使得立管产生涡激振动;
[0056]
s2中在分析相流动基本方程推导过程中,将两相分别按单相流处理并计入相间的相互作用,再将各相的方程加以合并,具体步骤为:
[0057]
根据流体微元单元段受力分析图,可得管内多相流流体微元段x方向上的力平衡
方程:
[0058][0059]
开采立管微元段x方向的力平衡方程:
[0060][0061]
剪切力q与弯矩m可以表示为:
[0062][0063][0064]
联立方程(1)、(2)、(3)和(4)可以得到两相流的深水立管顺流向(x方向)运动微分方程为:
[0065][0066]
同理可得到开采立管横流方向(y方向)上的运动微分方程:
[0067][0068]
式中:ei—弯曲刚度(n
·
m2);x—顺流向位移(m);y—横流向位移(m);z—立管轴向位置(m);t—截面张力(n);c—结构阻尼;f
x
(z,t)—顺流向海洋环境载荷(n);
[0069]fy
(z,t)—横流向海洋环境载荷(n);v
l
、vg—气、液体的速度(m/s);mr、m
l
与mg—单位长度立管质量、单位长度立管内液体与气体质量(kg);
[0070]
在轴向方向上,开采立管轴向上需要考虑立管自身重力与管内流体的重力,因此立管截面有效张力沿程是变化的。开采立管除了受到顶部张力与自身浮力外,其上部立管还要承受下部立管的重力,因此截面上的张力t=t(z),表示为:
[0071][0072]
式中:t
top
—顶部张力(n);ρw—海水密度(kg/m3);d0—立管外径(m)。
[0073]
进一步的,s3中的边界条件为使用newmark-β法计算开采立管在全局坐标系下振动控制矩阵方程
[0074]
其中[m]、[c]、[k]分别为质量矩阵,阻尼矩阵与刚度矩阵;分别为立管加速度、速度与位移矢量;{f}为外界的流体力载荷。
[0075]
采用经典的morison方程来讨论作用于开采立管的流体力。假设开采立管某深度位置的外部海流流速为uc,某时刻开采立管在海流下的截面示意图如图5所示。
[0076]
假设海流与开采立管之间的相对速度为vr,根据流体在立管顺流向与横流向上的相对速度,可以得到海流相对于立管的相对速度,基于相对流速可以得到作用于立管上的稳态拖曳力和升力。拖曳力方向沿着流体相对速度方向,升力方向垂直与拖曳力方向,作用于开采立管的外部流体力示意图如图6所示。
[0077]
根据图6可得:
[0078][0079][0080]
根据开采立管受力示意图,结合莫里森方程,可得到作用于立管上的稳态拖曳力以及升力:
[0081]
作用于开采立管的稳态拖曳力为:
[0082][0083]
作用于开采立管的稳态升力为:
[0084][0085]
根据受力分析,可得到作用于立管x与y方向上的流体力分量为:
[0086][0087]
作用于开采立管上的脉动拖曳力fd'和脉动升力f
l
'为:
[0088][0089]
其中,cd与c
l
分别为脉动拖曳力系数与脉动升力系数。
[0090]
因此,结合式(12)与(13)并联合式(9)可得作用于开采立管上的外部流体力为:
[0091]
[0092]
由于开采立管外截面是圆形截面,因此稳态升力系数通常取值为0。稳态拖曳力系数取值为1.2,将其带入式子(14)中。同时,假设忽略高阶项的影响,整理可得流体力表达式为:
[0093][0094]
式中,脉动拖曳力系数cd与脉动升力系数c
l
是通过尾流振子变量求得。
[0095]
立管运动微分方程中,mg、m
l
与νg、ν
l
是求解考虑管内气、液体分布的立管运动方程最为重要的物理参数,其取值大小取决于管内的多相流分布。在开采立管工作期间,管内流体沿开采立管上返到海面,随着温度升高且管内压力逐渐降低,气体逐渐膨胀、破裂,含气率增加,流型发生变化。针对气液两相流动问题,本文建立了计算流动参数的数值模型。
[0096]
气相质量守恒物理模型如图7所示。根据质量守恒定律,控制单元的质量变化是控制单元的输入质量减去输出质量。对于气相,控制单元的输入质量为:
[0097]
ρgν
geg
adt+qgdtdz(21)控制单元的输出质量为:
[0098][0099]
气体空隙率变化引起的内部质量变化为:
[0100][0101]
因此,管内气相连续性方程表示为:
[0102][0103]
开采立管内液相的连续性方程为:
[0104][0105]
式中:ρ—气/液相的密度(kg/m3);
[0106]
ν—气/液相的速度(m/s);
[0107]eg
—气体孔隙率(/);
[0108]el
—持液率;(/)
[0109]
a—管内截面积(m2);
[0110]bm
—管内液体局部体积系数;
[0111]rms
—局部溶液气体-液体比率;
[0112]
ρ
gs
—标准条件下的气体密度(kg/m3);
[0113]
qg—气体的产量(/)。
[0114]
根据动量守恒,物体动量随时间的变化率等于施加在物体上的外力之和。因此,可以得到如下动量方程:
[0115][0116]
对于气体,动量方程可以表示为:
[0117][0118]
对于液体,动量方程可以表示为:
[0119][0120]
在控制单元中,
[0121]eg
+e
l
=1(29)
[0122][0123]
因此,总动量方程可以写成:
[0124][0125]
在开采期间,管内的压力随着水深、管内流体质量、速度的变化而改变,管内的压降可以表示为:
[0126][0127]
管内气液相的质量可分别表示为:
[0128][0129]
式中:f
x
—为流动摩擦系数(/);
[0130]
d0—开采立管外径(m);
[0131]di
—开采立管内径(m);
[0132]
气液两相流的主要模式是气泡、段塞、搅动和环状流模式[36],表达形式如下:
[0133]
泡状流:
[0134]
ν
sg
≤0.429
×
ν
sl
+0.357
×
ν
oo
(34)
[0135][0136]
段塞流:
[0137]
ν
sg
>0.429
×
ν
sl
+0.357
×
ν
oo
(36)
[0138]
搅动流:
[0139][0140]
环状流:
[0141][0142]
式中:ν
sg
—气体的表面速度(m/s);
[0143]
ν
sl
—液体的表面速度(m/s);
[0144]
ν
oo
—气泡的极限上升速度(m/s);
[0145]
σ—表面张力(n/m)。
[0146]
s4中在开采立管振动求解过程中,保持运算时间步一致,将尾流振子模型与立管控制方程耦合迭代计算,可以分别得到开采立管顺流向与横流向的位移时程响应,通过对位移时程曲线进行快速傅里变换,可以得到开采立管顺流向与横流向的频率响应。
[0147]
s4中使用微分控制方程计算悬挂立管的无阻尼自由振动,得出模型的固有频率,对气液两相流和纯液流下的涡激振动特性对比分析、不同排量下的涡激振动特性对比分析、不同两相混合密度下的涡激振动特性对比分析和不同进气比下的涡激振动特性对比分析进行总结。
[0148]
如图8-11所示,对气液两相流和纯液流下的涡激振动特性对比分析,可知,与纯液流相比,在相同的外部激励下,由于管内气液两相流对立管的激励作用更强,导致其振动幅值更大,振动频率也更大,且表现出的模态响应也更高阶。因此,在实际工程问题中,应重点关注管内气液两相流对立管涡激振动的影响。
[0149]
对不同排量下的涡激振动特性对比分析,可知,在相同的外部激励下,随着排量的增加,管内气液两相混合速度的增大对立管刚度影响显著,从而其导致有效截面张力减小,振动幅值更大,振动频率减小,但表现出的模态响应保持不变。在实际工程中,在保证产量的情况下,可以适当降低管内流体排量,有利于减小立管涡激振动。
[0150]
对不同两相混合密度下的涡激振动特性对比分析,可知,在相同的外部激励下,随着两相混合密度的增加,立管的振动幅度增大,振动响应频率降低,模态振型保持不变。
[0151]
对不同进气比下的涡激振动特性对比分析,可知,在相同外部激励下,随着进气比的增加,立管的固有频率降低,振动响应频率和振动幅值均增大,且进气比的增加会激发立管更高阶的模态振型。
[0152]
本发明可以得出以下结论:
[0153]
通过建立内外流耦合作用下的立管系统动力学分析模型,分析含气液两相流的立管的涡激振动特性,研究了排量、混合流体密度、进气比等因素对涡激振动的影响,可以知道,气液两相流会引起立管的固有频率降低,在相同的外部激励下,导致其振动幅值更大,振动频率也更大,且表现出的模态响应也更高阶,因此,在实际工程问题中,应重点关注管
内气液两相流对立管涡激振动的影响。
[0154]
随着管内流体排量和流体混合密度的增加,立管轴向张力减小,固有频率降低,振动幅值增大。由于此时发生频率锁定现象,振动频率呈现与固有频率同步降低的趋势,模态振型保持不变。
[0155]
进气比的增大导致立管刚度减小,振动幅值增加。随着管内含气量的增加,管内两相流会诱发立管更高一阶模态的振动,振动频率和模态振型均增大。
[0156]
可以理解,本发明使通过一些实施例进行描述的,本领域技术人员知悉的,在不脱离本发明的精神和范围情况下,可以对这些特征和实施例进行各种改变或等效替换。另外,在本发明的教导下,可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此,本发明不受此处所公开的具体实施例的限制,所有落入本技术的权利要求范围内的实施例都属于本发明所保护的范围内。

技术特征:


1.内外流耦合作用下的开采立管涡激振动特性分析方法,其特征在于,提出了考虑内外流耦合作用下的开采立管数值分析模型和管内多相流模型分析方法,包括开采立管力学模型分析方法、流体力模型分析方法、尾流振子模型分析方法和管内两相流模型分析方法,采用有限元法和newmark-β法对模型进行求解,具体求解步骤如下所示:s1:根据外部海流速度大小、立管长度、立管内外径、内流速度等参数对深海含两相流立管进行模拟,模拟开采立管涡激振动过程;s2:分析开采立管单元和管内流体单元段受力,建立考虑内外流耦合作用下的开采立管力学模型,此基础上,建立考虑管内流型对立管动力响应影响的两相流模型,计算管内流体质量和速度分布;s3:将计算出的管内流体质量和速度分布情况与开采立管模型结合,基于模型的边界条件采用newmark-β法,求解考虑内外流耦合作用下的开采立管力学模型分析方法;s4:根据计算出的考虑内外流耦合作用下的立管动力响应,模拟分析立管的涡激振动强弱以及振动频率的大小,得出内外流耦合作用下的开采立管涡激振动响应规律。2.根据权利要求1所述的内外流耦合作用下的开采立管涡激振动特性分析方法,其特征在于,所述s1中的两相流立管管内的两相流能够对开采立管产生动态激励作用,在开采立管运移两相流的时候,进一步将开采立管简化为两端简支的内输多相流立管;立管置于横向流动的海流中产生大变形弯曲,同时由于海水绕过圆柱管道会因旋涡脱落而产生涡激力作用于开采立管上,进而使得立管产生涡激振动。3.根据权利要求1所述的内外流耦合作用下的开采立管涡激振动特性分析方法,其特征在于,所述s2中在分析相流动基本方程推导过程中,将两相分别按单相流处理并计入相间的相互作用,再将各相的方程加以合并,具体步骤为:根据流体微元单元段受力分析图,得到管内多相流流体微元段x方向上的力平衡方程为:开采立管微元段x方向的力平衡方程位:剪切力q与弯矩m为:剪切力q与弯矩m为:联立上述方程得到两相流的深水立管顺流向(x方向)运动微分方程为:得到开采立管横流方向(y方向)上的运动微分方程为:
式中:ei—弯曲刚度(n
·
m2);x—顺流向位移(m);y—横流向位移(m);
z
—立管轴向位置(m);t—截面张力(n);c—结构阻尼;f
x
(z,t)—顺流向海洋环境载荷(n);f
y
(z,t)—横流向海洋环境载荷(n);v
l
、v
g
—气、液体的速度(m/s);m
r
、m
l
与m
g
—单位长度立管质量、单位长度立管内液体与气体质量(kg);在轴向方向上,开采立管轴向上根据立管自身重力与管内流体的重力,立管截面有效张力沿程是变化,开采立管除了受到顶部张力与自身浮力外,其上部立管还要承受下部立管的重力,因此截面上的张力t=t(z),表示为:式中:t
top
—顶部张力(n);ρ
w
—海水密度(kg/m3);d0—立管外径(m)。4.根据权利要求1所述的内外流耦合作用下的开采立管涡激振动特性分析方法,其特征在于,所述s3中的边界条件为使用newmark-β法计算开采立管在全局坐标系下振动控制矩阵方程其中[m]、[c]、[k]分别为质量矩阵,阻尼矩阵与刚度矩阵;{u}分别为立管加速度、速度与位移矢量;{f}为外界的流体力载荷。5.根据权利要求1所述的内外流耦合作用下的开采立管涡激振动特性分析方法,其特征在于,所述s4中在开采立管振动求解过程中,保持运算时间步一致,将尾流振子模型与立管控制方程耦合迭代计算,分别得到开采立管顺流向与横流向的位移时程响应,通过对位移时程曲线进行快速傅里变换,得到开采立管顺流向与横流向的频率响应。6.根据权利要求1所述的内外流耦合作用下的开采立管涡激振动特性分析方法,其特征在于,所述s4中使用微分控制方程计算悬挂立管的无阻尼自由振动,得出模型的固有频率,对气液两相流和纯液流下的涡激振动特性对比分析、不同排量下的涡激振动特性对比分析、不同两相混合密度下的涡激振动特性对比分析和不同进气比下的涡激振动特性对比分析进行总结。

技术总结


本发明公开了内外流耦合作用下的开采立管涡激振动特性分析方法,本文建立了考虑内外流耦合作用下的开采立管数值分析模型,以确定含两相流立管的涡激振动特性。采用有限元法和Newmark-β法对模型进行分析求解,总结了内外流耦合作用下的开采立管涡激振动响应规律。本发明通过外部海流与内部多相流耦合作用下的开采立管涡激振动特性分析方法为多相流开采立管的安全性提供有效的理论指导。立管的安全性提供有效的理论指导。立管的安全性提供有效的理论指导。


技术研发人员:

毛良杰 严巨熙 黄鑫 付强

受保护的技术使用者:

西南石油大学

技术研发日:

2022.10.21

技术公布日:

2022/12/30

本文发布于:2024-09-21 21:49:08,感谢您对本站的认可!

本文链接:https://www.17tex.com/tex/2/49777.html

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系,我们将在24小时内删除。

标签:管内   流体   两相   模型
留言与评论(共有 0 条评论)
   
验证码:
Copyright ©2019-2024 Comsenz Inc.Powered by © 易纺专利技术学习网 豫ICP备2022007602号 豫公网安备41160202000603 站长QQ:729038198 关于我们 投诉建议