数值模拟1各种方法总结1.docx

上传人:b****6 文档编号:4301512 上传时间:2022-11-29 格式:DOCX 页数:13 大小:47.44KB
下载 相关 举报
数值模拟1各种方法总结1.docx_第1页
第1页 / 共13页
数值模拟1各种方法总结1.docx_第2页
第2页 / 共13页
数值模拟1各种方法总结1.docx_第3页
第3页 / 共13页
数值模拟1各种方法总结1.docx_第4页
第4页 / 共13页
数值模拟1各种方法总结1.docx_第5页
第5页 / 共13页
点击查看更多>>
下载资源
资源描述

数值模拟1各种方法总结1.docx

《数值模拟1各种方法总结1.docx》由会员分享,可在线阅读,更多相关《数值模拟1各种方法总结1.docx(13页珍藏版)》请在冰豆网上搜索。

数值模拟1各种方法总结1.docx

数值模拟1各种方法总结1

油藏数值模拟基本过程

一、数值模拟发展概况

30年代人们开始研究地下流体渗流规律并将理论用于石油开发;

50年代在模似计算的方法方面,取得较大进展;

60年代起步,人们开始用计算机解决油田开发上的一些较为简单间题,由于当时计算机的速度只有每秒几万到几十万次,实际上只能做些简单的科学运算;

70年后主要体现于计算机的快速升级带动了油藏数模的迅猛发展,大型标量机计算速度达到100--500万次,内存也高增主约16兆字节。

在理论上黑油模型计算方法更趋成熟,D.W.Peaceman的<油藏数值模似基础>以及K.Aziz和A.Settari的<油藏模似>等主要著作都是在这个阶段出版的,但仍受到计算机速度和内存的限制,使用的方法一般仅限于IMPES及半隐式等,只能解决中小型油藏的模拟应用问题;

80年代则是油藏数值模似技术飞跃发展的年代,解决不同类型油藏的数模计算方法及软件相应问世,同时超级向量机的诞生,使计算机速度达到亿次,甚至几十亿次,内存高达10—20亿字节。

90年代特别是后期,油藏模似软件各模块功能也有了惊人的发展,主要体现为向一体化方面发展;即集地震、测井、油藏工程(数模)、工艺及地面集输、经济评价等为一体的大型软件方面发展。

目前油藏数值模似软件基本上形成了一套能处理各种类型油气藏和各种不同开采方式的软件系列。

?

黑油模型已被广泛用于各种常规油气藏的模拟;

?

裂缝模型可用来解决除砂岩以外的灰岩、花岗岩、凝灰岩和变质岩的裂缝性油气藏开发问题;

?

组分模型用于凝析气藏、轻质油、挥发油藏的开发设计和混相驱的研究;

?

热采模型用于稠(重)油油藏蒸气吞吐、蒸汽驱和就地燃烧的设计;

?

化学驱模型用于在注入水中添加聚合物、表面活性剂、碱等各种化学剂进行三次采油提高采收率的计算和设计。

油藏数值模拟方法的新突破

随着计算机运算速度的提高,向量算法的出现和应用是软件设计上一个划时代的发展。

预处理共轭梯度法更快速、有效地解各种更为复杂和困难的大型稀疏线性方程组。

网格化方面不局限于静态和动态的局部网格加密技术,不规则网格、PEBI网格的出现更好的解决了在边界、断层插值计算以及面与面垂直正交的新型数模计算方法,更快速收敛。

此外多重网格法、混合有限元法、流线法等都在逐步完善和发展。

同时,并行处理技术给大中型油田数模工作带来了生机。

二、数值模拟基本原理

以渗流力学为基础建立数值模型,即通过一组方程组,在一定的假设条件下,考虑油藏构造形态、断层位置、砂体分布、储层孔渗饱等参数的变化;流体高压物性变化;不同岩石类型;不同渗流驱替特征曲线(相渗);井筒垂直管流等描述油藏真实的物理过程。

主要包括:

运动方程、状态方程和连续方程。

1、运动方程

引入单相(或多相)流的达西定律方程:

单相流达西定律:

多相流达西定律:

其中:

K:

空气渗透率;达西

A:

截面积;厘米2

L:

长度;厘米

△P:

压差;大气压

Kro、Krw:

油、水相对渗透率

μo:

原油粘度;厘泊

转变为运动方程:

2、连续性方程

研究流入流出单元体中质量的变化方程。

根据物质平衡原理,流入单元体中的流体流量减去流出单元体流体流量等于单元体流体质量变化。

单元体中流体质量发生变化的速率

将达西定律代入连续性方程得出油、水的流动方程:

油水

3、状态方程

为求解上述两个方程中的两个未知数P、Sw,引入状态方程:

(1)

(2)h、、为常数

(3)

(4)、为的函数So=1-Sw

由连续性方程推导出油、水的流动方程,从而引入压力方程:

饱和度方程引入含水百分数:

最后把两维两相流动方程写成如下形式:

解上述方程通过把微分方程离散化变代数形式来求解。

三、数值模拟的目的

(一)、为什么开展油藏数值模拟工作

研究和开发一个油田是一个复杂的综合性的科技问题,高精度的地震资料的处理解释提供研究区域的构造、断层、边界及其走向,但地震纵向分辨率受到限制,不能很好的反映一个同相轴(地震道)中沉积砂体的物性变化特征;测井可较好的反映到小于1米以下沉积砂体的物性特征,提供可靠的地层对比结果。

但作为新老油田开发方案的研究及剩余油分布的研究,是地震、地质、测井理论方法都无法做到的。

地质上仅定性或半定量分析,测井用于生产监测不能以点带面。

惟独油藏数值模拟工作可再现生产历史,定量分析剩余油潜力;并做到室内研究投入少、时间短,还可进行开发方案优选及经济评价工作。

所以总公司强调开发方案的部署一定要开展数值模拟工作。

值得强调的是油藏数值模拟工作提倡一体化,注重前期的地震解释和测井解释即油藏描述工作。

(二)、油藏数值模拟的目的

在进行油藏数值模拟工作前,首先应根据油田开发过程中存在难以解决的实际问题,提出开展此项工作的目的及意义,即最终所要达到解决问题的目标是什么?

一般通过油藏数值模拟可进行以下研究工作:

1.初期开发方案的模拟

1).评价开发方式;如:

枯竭开采、注水开发等。

2).选择合理井网、开发层系、确定井位;

3).选择合理的注采方式、注采比;

4).对油藏和流体敏感性研究。

2.对已开发油田历史模拟

1).核实地质储量,确定基本的驱替机理(如:

是天然驱,还是注水开发。

);

2).确定产液量和生产周期;

3).确定油藏和流体特性;

4).提出问题、潜力所在区域。

3.动态预测

1).开发指标预测及经济评价

2).评价提高采收率的方法(如:

一次采油、注水、注气、化学驱等)

3).剩余油饱和度分布规律的研究,再现生产历史动态

诸如:

研究剩余油饱和度分布范围和类型;

单井调整:

改变液流方向、注采井别、注水层位;

扩大水驱油效率和波及系数;

4).潜力评价和提高采收率的方向

诸如:

?

确定井位、加密井的位置;

?

确定油田开发最大产液量、产量对采收率的影响;

?

确定地面和井的设备。

5).专题和机理问题的研究

诸如:

?

对比注水、注气和天然枯竭开采动态;

?

研究各种注水方式的效果;

?

研究井距、井网对油藏动态的影响;

?

研究不同开发层系对油藏动态的影响;

?

研究注水速度对产油量和采收率的影响;

?

研究油藏平面性质和层间非均质性对油藏动态的影响;

?

验证油藏的面积和地质储量;

?

校验油藏数据;

?

为谈判和开发提供必要的数据。

注意

无论是对油藏进行初期开发方案、已开发油田历史模拟,还是动态预测的数值模拟工作,都要求油藏工程师要有针对性的拟定出能解决油田开发实际问题的数值模拟工作详细计划,及其开展此项工作的目的和应达到的目标是什么。

四、选择适当的数值模型及相类

对一般油藏而言,可具有两相(油+水、油+气、气+水)和三相(油、气、水)的模型,而油藏维数可具有一维平面、垂直模型、两维平面(剖面)、两维径相(锥进)模型和三维平面等模型。

但实际油藏其地下流体渗流机理、岩石及流体性质等地质特征不同;生产过程中的开采方式、机理不一等复杂问题,促使非常规油藏模型孕育而生。

如:

注意

如何确定多相问题?

需了解油藏流体性质、生产特征,以及采油机理。

如:

?

气藏是否带油环和下伏水区;

?

有无明显相态变化的混相驱油藏;

?

有无表现出凝析现象的相态特征;

?

气顶及水锥效应是否占优势的油藏等;

由此,在模拟区已确定的同时,需要认真考虑选择模型及相类的问题。

五、数据录取的准备工作

在数据录取工作中,建立精细的油藏地质模型是模拟工作成败的关键;而对生产数据的正确分析、合理取舍、根据生产情况合理分层和划分时间阶段步长,即录入生产数据卡工作也极为重要。

油藏模拟工作要研究油田开发中的问题很多,不可能一律对待,也不可能在一个模型里都予以满足,需根据研究任务和客观许可条件建立相应的地质模型。

(一)建立油藏地质模型

1.在油藏描述的基础上建立油藏地质模型

开展油藏描述工作,对油藏的地质、油层非均质特征,沉积相的详细描述和研究,根据油藏沉积相研究建立该油藏特征的沉积模式。

油藏描述分析的目的是综合所有的测井、岩心和生产测试等资料来得出一个与全油田一致的储集层模型。

对各种未知的基本参数例如:

对顶面深度、砂厚、孔、渗、饱等空间分布的评价中最大限度地发挥现有测井资料的作用,同时将这些参数结合所需储集层的几何特性参量进行计算,并结合地质沉积相分析提供出更为精细、完善的油田地质模型。

确定一个油藏地质模型所需的许多参数,在油藏勘探试采阶段或初期刚投入开发阶段用有限几口井的资料进行计算、解释及建模其精度是不高的,这些参数初始误差越大,则通过历史拟合达到令人满意的油藏特征描述所需的时间越长。

所以,无论在对老区或新区进行数值模拟时,应对所选区块选用所有的井(特别是“关键井”)开展油藏描述工作,进行全面的分析研究。

2.运用油藏描述技术对地质沉积模型进行地层网格化

根据实际油田数值模拟所建立网格的需要,进行油藏描述工作时,即可对沉积地质模型进行相应的网格化,这样网格化的模型既能代表油藏沉积相的地质特征又有利于计算机数模作。

3.运用油藏描述技术提供相应沉积模式的地质参数

根据建立的油层沉积模式,提供砂岩、孔、渗、饱等参数的分布模式,并核实地质储量。

(二)网格选择

通过以上工作,可建立更为精细、完善的地质模型。

而在目前的实际工作中,大多数的老区块并没有开展过此项工作。

仅能做得是将各井点参数按网格划分的层系进行厚度加权平均,其结果是厚层贡献值大,薄层贡献值低,实际模拟中很有可能掩盖这样或那样的一些矛盾。

如:

现有测试资料反映,注水开发油田,注入水沿高渗薄层(厚度为0.5m)突进,造成水淹。

而要降低这样的矛盾,只能把焦点倾注于网格层系的划分上,这便是以下所要谈到的问题。

在一个区块的模拟研究中,合理有效的选择网格必须考虑以下几个问题。

1.网格的定向

?

网格的界限要与天然的非流动边界相符合,包括整个系统的矩形网格应最大可能的重迭在油藏上;

网格应包含有所有的井位(包括即将完善的新井、扩边井);

网格方向要与流体流动的主要方向(沿主河道方向,即平行渗透率的主轴)和油藏内天然势能梯度吻合;

网格的定向尽量减少死接点数目。

2.网格的尺寸

网格越多,每个时间步长中所需计算的数学问题越多,机时费用越多;

当时间步长由最大饱和度所控制时,较小的网格通常使最大可允许时间步长减小;

一般邻井之间至少要有2~3个空网格或更多,使其能反映油藏结构和参数在空间的连续变化,同时足够的网格能控制和跟踪流体界面的运动;

如果模拟前考虑井网加密方案,应确定适当的加密井井位和网格尺寸;

3.纵向网格的划分

?

除考虑本身按沉积韵律划分的层系外,还要考虑生产过程中的整体改造工作,如:

补层、压裂酸化、堵水等等;同时,对生产特征(如:

底水锥进、气顶等)都应合理考虑,对一个层系中的细分小层问题更是如此;

4.不规则网格的选择

需按实际情况酌情考虑,如井网密度大的井区,相邻的油、水井之间用一个空白网格分开,或处在相邻的网格中也是允许的,必要情况下可考虑井点网格加密。

(三)数据录入准备

1.表格数据

主要指岩石物性、流体性质。

油气PVT数据表(高压物性分析);

水及岩石PVT性质(高压物性分析);

油水相对渗透曲线;

毛管压力曲线(岩心压汞实验)等;

井筒流动数据。

2.网格数据

油藏顶面海拔深度;

砂层厚度(有效);

孔、渗、饱参数、岩石类型等。

3.动态数据

完井数据:

射孔、补孔、压裂、堵水、解堵日期、层位、井指数等;

生产数据:

平均日产油、日产水、日产气、平均油气比和含水比等;

压力数据:

井底流压、网格压力等。

动态监测资料(分层测试、吸水、产液剖面等)

4.其它数据

主要包括算法选择、输入输出控制、油水井约束界限、油井定压定产等参数。

为对今后数模工作数据资料录入规范化,特制定以下几类表格形式(参见附表)。

注意

对于黑油油藏,PVT数据极为重要,其数据由地层体积系数、溶解油气比和粘度作为压力的函数表所组成;

表格数据要求变量与自变量之间的关系要光滑,如:

不光滑的油水两相渗透率曲线,将导致拟合的含水曲线差枝不齐,导致迭代不收敛等问题。

六、历史拟合方法及技巧

数值模拟过程(特别是历史拟合)是一项复杂的、消耗人力和机时的繁琐工作,如不遵循一定步骤,掌握一定技巧,可能陷入难以解脱的矛盾之中。

一般认为,同时拟合全区和单井的压力、含水和油气比难以办到,必须将历史拟合过程分解为相对比较容易的步骤进行。

历史拟合一般采取以下几个步骤:

1确定模型参数的可调范围;

2对模型参数全面检查;

3历史拟合;

1).全区和单井压力拟合;

2).全区和单井含水拟合;

3).单井生产指数拟合。

(一)确定模型参数的可调范围

确定模型参数的可调范围是一项重要而细致的工作,需收集和分析一切可以利用的资料。

首先分清哪些参数是确定的,哪些参数是可调的。

资料及专家介绍:

孔隙度允许修改范围±30%;

渗透率视为不定参数,可修改范围±3倍或更多;

有效厚度,由于源于测井资料,与取心资料对比偏高30%左右,主要是钙质层和泥质夹层没有完全挑出来,视为不定参数,可调范围-30%左右;

流体压缩系数源于实验室测定,变化范围小,视为确定参数;

岩石压缩系数源于实验测定,但受岩石内饱和流体和应力状态的影响,有一定变化范围;同时砂岩中与有效厚度相连的非有效部分,也有一定孔隙和流体在内,在油气运移中起一定弹性作用。

因而,允许岩石压缩系数可以扩大一倍;

相对渗透率曲线视为不定参数,允许作适当修改;

油、气的PVT性质,视为确定参数;

油水界面,在资料不多的情况下,允许在一定范围内修改。

(二)对模型参数全面检查工资

油藏数值模拟的数据很多,出现错误的可能性很大。

为此,在进行历史拟合之前,对模型数据进行全面检查是十分必要的。

数据检查包括模拟器自动检查和人工检查两方面,缺一不可。

模拟器自动检查包括:

1、各项参数上下界的检查

对各项参数上下界的检查,发现某一参数超过界限,打出错误信息。

1).检查原始地质储量并与容积法计算进行比较;

N=7758?

A×h×Φ×Soi/Boi

2).检查所有原始油藏性质图和输入数据。

2、平衡检查

在全部模型井的产率(注入率)都指定为零的情况下,进行一次模拟计算,其结果应是油藏状态参数(压力场和饱和度场)应该与油藏初始状态参数一致,无任何明显变化,流体应该是处于平衡状态。

否则表明参数有了问题,需重新检查模拟卡中的相关参数。

人工检查包括:

1)、不同来源的资料相互对扣;

2)、日产(注)量、月产(注)量和累积量相互对扣;

3)、物质平衡检查,即分析全区压力变化与累积净注入量(或亏空)的关系是否一致;

4)、对串槽井的产水量进行修正。

(三)历史拟合分两步

,首先是拟合全区压力到拟合单井压力,然后是饱和度(全区和单井含水)拟合。

1.压力拟合

油藏中产生压力分布是由于流体场产、注流动后重新分布的结果,由达西定律所描述:

(x=0,g,w)

公式中包含出现的不定参数K和Kr,改变其中之一将有效地改变压力的重新分布,由此当网格单元井的压力形式出现异常上升或下降时,表明不是真实地层条件下的K和Kr值。

而地层平均压力水平的主要影响因素是由油藏总孔隙体积和总压缩率来确定,两者有如下关:

而:

由此可以确定造成不正确压力大小的不定参数是Cr和Φ,通常改变Cr可以观测到压力的显著变化,而±30%的Φ变化对计算出来的油藏动态参数(包括压力)可能差别不大,是次要因素。

A:

拟合全区或单井油层压力

方法

1.修改Cr、Φ、h、S(饱和度)参数

1)一般如果一个区块有室内岩心Cr数据,建议不改变其值大小,如按拟合情形需改变时,也不能超过一倍为好;

2)如果拟合油藏压力水平过高,则往往表示油藏地质储量过高。

此时需减少Φ、h以及S值,以达到使地层压力水平降低的目的;

2.修改渗透率值改变流体流动方向

以达到改变油层压力的目的;

增加低压带的渗透率,以达到提高低压带的压力,反之亦然。

检查原始地层压力梯度、原油体积系数、脱气油密度,以校正地层压力水平;

1).原始地层压力与深度关系将直接影响到整个地层的压力水平。

如果油藏压力水平过高或偏低,首先应检查输入的基准面深度和相应的油相压力,是否符合压力梯度关系,如果不符合则需修改这个相参数。

2).从压力梯度数据中效核地下原油密度;

压力梯度:

dP/dD=(P2-P1)/(D2-D1)

通过单位换算为地下原油密度,与实际输入值比较,如(脱气原油密度)有误也需要修改。

3).检查原油体积系数

检查地面脱气原油密度,如果正确,则按:

Bo=ρ地面/ρ地下

求出原油体积系数,进行修改。

总之,先拟合全区压力,然后拟合单井压力形式。

而单井压力拟合主要是靠修改井局部地区的渗透率或方向渗透率。

注意在进行全区压力拟合时,要兼顾单井点的情况进行修改。

而且,同时要照顾到单井点的含水拟合。

这样,压力拟合阶段对方向渗透率的修改就有利于以后含水的拟合,节省机时费用。

2.全区和单井含水拟合(饱和度的拟合)

压力拟合达到满意的效果后,将进行全区和单井含水拟合。

油藏中流体饱和度的分布,影响井的注入量和采出量,即影响油水比和油气比。

瞬时油水比(WOR)和油气比(GOR)由一下公式计算:

B:

全区和单井含水的拟合

其方法步骤如下:

1.整相渗曲线

首先,要检查相渗曲线是否平滑,才不至于导致计算含水值过高或过低。

其次,对初期拟合含水偏低的情形,可适当左移水相渗透率曲线(即抬高水相渗透率曲线),对高含水期拟合含水偏低的情形,可适当右移油相渗透率曲线(即抬高油相渗透率曲线)。

反之亦然。

(如图所示)

再次,尽可能根据分采层的含水上升率曲线反推几组相渗曲线,用以代表不同类型产层的渗流机理。

由于我国油藏大多属于陆相湖盆沉积,物源近、以及多物源方向供给碎屑物质,造成沉积相带窄,非均质性严重,砂体类型也多。

加上沉积受多级旋回的控制,形成多层系含油的特点,无论从岩性还是从岩相上变化都很悬殊。

纵向上各层间渗透率差别很大,平面上连通性差、砂岩体往往在短距离内就尖灭、交叉或迭加。

而东部多为断陷含油气盆地,断层发育,构造复杂,致使各断块间油气水分布关系难以摸清。

而且,更为严峻的问题是开发过程中对达到一个油藏或小到一个断块的五项渗流特征参数的岩心实验资料录取少,一般应用一组相渗曲线很难表征地下油藏各层系、小层内流体的真实渗流机理,更别说对无资料的地区要借用相同或类似地区的资料来使用。

2.调整毛管压力曲线

拟合含水还需检查毛管压力曲线,以改变束缚水饱和度和初始含水饱和度分布。

3.局部井点含水的拟合

通过如下修改:

?

改变含水区地质储量,如调整这些地区的孔隙度Φ、渗透率K或流体S值的大小,以达到含水饱和度的拟合;

?

减少与水区连通部位的渗透率值,以控制含水上升的目的;

?

在局部地区含水拟合差别较大时,可调整X、Y方向渗透率,即AKX、AKY、BKX、BKY、AKZ(纵向),以达到在不增加地质储量的条件下,增加或减少流体沿某一方向的流动性,实现含水的拟合。

油藏数值模拟流程图

附件1:

关于实测压力的皮斯曼校正

由于模拟计算结果输出的网格压力通常是折算成模拟区域内的基准面处的压力,不同于通常测试或计算的地层压力,不能进行比较。

因此压力校正第一步是把代表不同深度的每口井实测地层压力折算成与模型相同的基准面深度处的压力,目的是消除深度的影响。

Pwsd------折算到D处的压力,MPa;

Pws----关井后稳定压力,Mpa;

Dd----基准面深度,m;

Dms----油层中部海拔深度,m;

Pd----油的压力梯度,Mpa/m。

ρ0=r0●ρw

ρg=rg●ρair

ρ0、ρg、ρw---分别为油、气、水的密度g/cm3

r0----空气相对密度g/cm3

rg----气相对密度g/cm3

Rs----气油比m3/m3

B0----原宥体积系数。

拟合完成后,能与实际压力比较的是井网格压力,但它又不同于压力恢复时测得的短期关井井底压力,也不同于通常计算的地层压力,不能进行直接比较。

如何进行校正?

1、如果油井测得压力恢复数据,并出现直线段后,用下式校正实测压力后才可与网格压力对比:

式中:

k---网格节点平均渗透率

Δt----关井时间

Δx、Δy-------x、y方向的网格步长

Ct----总压缩系数

μ----粘度

φ----孔隙度

2、如果油井测得压力数据是一个压力点,该点位于压力恢复曲线直线段上,

则:

首先将质量总流率(产量)Q换算成体积总流率q,

又:

q=Q[fw(1?

1/r0)+1/r0],然后用下式计算出校正的实际地层压力。

式中:

P0----校正的实际地层压力

Pwsd----折算到Dd处的压力

q-----关井前稳定总产率

----总流度

H----射开厚度

-----孔隙度

------关井时间

、------网格x、y方向的步长

附件2:

关于烃类有效孔隙体积的计算

在渗流过程中,地层的孔隙体积实际上是发生变化的,它与地层压力和压实作用有关:

式中:

----压缩因子

----节点当前平均压力

----原始孔隙体积

----岩石压缩系数

----节点当前压力

----基准面压力

----变化后的孔隙体积

当时,说明孔隙体积没有变化,在简化三维三相黑油模型中,也常常设,即岩石是不可压缩的,但这对模拟结果可能带来某种影响,在空间某一微小的单元里,虽然孔隙体积变化很微弱,但它可能影响到烃类有效孔隙体积的大小,也就是说,对剩余油的分布有一定影响。

因此,需进行烃类有效孔隙体积的计算。

为什么?

通常情况下,数值模拟提供各时间阶段的剩余油分布来,并不意味着剩余油饱和度大的地方,就是挖潜或调整的重点,因为剩余油饱和度与油层厚度没有直接关系,因此,其分布状况往往掩盖了厚油层与薄油层的区别,为了使剩余油分布从定性向定量转化,可用烃类(包括油和水)有效孔隙体积PV来表示:

PV=DX?

DY?

H?

?

?

(Sor?

Bo?

r0/r0r+Sgr?

Bg?

rg/rgr)

其中:

DX、DY–---为网格节点X、Y方向步长,

H---油层有效厚度,

Sor、Sgr---为剩余油、气饱和度,

Bo、Bg—为油、气体积系数,

r0、r0r和rg、rgr---为油、气地面和地下密度。

 

双孔双渗,就是模型中有基质和裂缝两种孔隙体积,基质孔隙是主要的储油空间,裂缝是主要的流动通道,基质和裂缝都有孔隙体积和渗透率,所以叫双孔双渗。

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 初中教育 > 理化生

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1