首页> 中国专利> 超深井稠油掺稀比例确定的优化方法及其掺稀混配器

超深井稠油掺稀比例确定的优化方法及其掺稀混配器

摘要

本发明是超深井稠油掺稀比例确定优化的方法及其掺稀混配器,方法包括以下步骤:1)建立多相流综合压降计算模型和修正模型;2)推导出掺稀井井筒中流体温度梯度的计算方法;3)进行超深井弯曲井段的多相流计算,并计算出井眼形状中参数的分布规律;4)建立稠油掺稀生产动态分析模型;5)对稠油掺稀比例及生产参数进行优化。本发明实现对超深井稠油掺稀生产动态的模拟和分析。定量分析、模拟超深井稠油掺稀生产动态,得出压力、密度、粘度等参数沿井筒的分布规律。并可根据油井的产能,对稠油掺稀比例及生产参数进行优化。优选出最佳的掺稀比例及掺入压力。具有极强的工业实用性。

著录项

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2013-06-26

    授权

    授权

  • 2010-05-12

    实质审查的生效 IPC(主分类):E21B43/16 申请日:20080928

    实质审查的生效

  • 2010-03-31

    公开

    公开

说明书

【技术领域】

本发明涉及一种稠油掺混方法及其及工具,特别是一种超深井稠油掺稀比例确定的优化方法和使用掺稀混配器。

【背景技术】

众所周知,采取降粘的措施是开采稠油的有效手段。这些降粘措施包括火烧油层、蒸汽吞吐、化学降粘、电加热等。特别是,超深井的稠油具有其特殊性。由于油层埋藏深度在5000m以上,地层温度高达130℃,以上,因此,原油在储层中具有较好的流动性。但是,原油进入井筒后再向地面流动的过程中随着井筒温度的降低,稠油的粘度随之增大,将逐渐失去流动性。因此,对超深井的稠油开采而言,行之有效的降粘手段应该在井筒举升过程中完成。研究表明,大部分超深井稠油的开采需要采取井筒降粘的辅助措施,而掺稀工艺是各种降粘工艺中应用效果最好、适用性最广的降粘工艺。随着油田开发进程的不断深入,掺稀井含水的上升,且掺稀井数不断增多,如何确定合理的掺稀比例,以最小的掺稀量,获得最大的增油量是掺稀工艺亟待解决的问题。

【发明内容】

因此,本发明的找到一种超深井稠油掺稀比例确定及优化方法,利用该方法能够达到如下目的:

1、能够实现对超深井稠油掺稀生产动态的模拟和分析。

2、可定量分析、模拟超深井稠油掺稀生产动态,得出压力、密度、粘度等参数沿井筒的分布规律。

3、并可根据油井的产能,对稠油掺稀比例及生产参数进行优化。

4、优选最佳的掺稀比例及掺入压力。

5、研制出一种新型井下掺混器;

6、能够指导现场采取合理的工作制度,划分掺稀的经济界限,取得更高的经济效益。

具体技术方案:超深井稠油掺稀比例确定优化的方法,包括以下步骤:

1)、建立多相流综合压降计算模型和超深井稠油原油物性模拟的修正模型;

2)、推导出超深井稠油掺稀井的井筒中流体温度梯度的计算方法;

3)、进行超深井弯曲井段的多相流计算,并计算出井筒状中的温度、压力、密度、粘度参数的分布规律;

4)、建立稠油掺稀生产动态分析模型;

5)、对稠油掺稀比例及生产参数进行数据优化。

所述的步骤1)中多相流综合压降计算模型的建立进一步包括:

步骤1.1-1)、根据质量守恒、动量守恒和能量守恒定理建立气液两相流的基本控制方程模型;

步骤1.1-2)、根据管道中气液两相流动的描述简化步骤1.1-1)的模型,首先,若管道的长度远远大于管道直径,就可认为管道横截面的流动参数的平均值不变,即流动参数仅在轴线方向变化;其次,在考虑流动结构后,可将流场分为若干区域,在各个区域内,将流动处理为单相的或多相混合的;

步骤1.1-3)、根据压力梯度表达式,确定倾斜管多相流压降的计算方法。

步骤1.1-1)根据质量守恒、动量守恒和能量守恒定理建立气液两相流的基本控制方程模型具体包括:

A)、由质量守恒原理,得到两相流动的连续方程可表示为

ddtV1(t)ρ1dV1+ddtV2(t)ρ2dV2=-A1(t)ρ1U1·n1dA1-A2(t)ρ2U2·n2dA2---(1-1)

式中

ρk——为k相(k=1,2)的密度,kg/m3

Uk——为k相(k=1,2)的速度矢量,m/s;

nk——相应表面的单位外法向矢量;

在多相流动的流场中考虑一个固结于惯性坐标系被界面AI(t)切割的确定的控制体V,它由相界面AI(t)分为两个分别由两相占据的子空间域V1(t)和V2(t),V1(t)由表面A1(t)和AI(t)围成,V2(t)由表面A2(t)和AI(t)围成,由A1(t)和A2(t)组成控制面;

B)、由动量守恒原理知,

控制体V中动量对时间的变化率等于:

(1)单位时间通过控制面A进入的流体动量通量;

(2)作用在控制体中和控制面上所有外力的合力,这些外力由体积力(即重力)和表面力(即应力张量)组成:

ddtV1(t)ρ1U1dV1+ddtV2(t)ρ2U2dV2=-A1(t)ρ1U1(U1·n1)dA1

-A2(t)ρ2U2(U2·n2)dA2+V1(t)ρ1fdV1+V2(t)ρ2fdV2---(1-2)

+A1(t)n1·T1dA1+A2(t)n2·T2dA2

式中

f——是作用在流体上的单位质量力,N;

Tk——为第k相流体在控制面上的应力张量,它与应变率之间的关系由流体的本构方程给出;

通常,应力张量可表示为:T=-pI+τ

式中

p——为压力,Pa;

I——为单位张量;

τ——为粘性应力张量;

对牛顿流体:

τik=μ(vixk+vkxi-23δikvlxl)+κδikvlxl---(1-3)

式中

μ——为剪切粘性系数,Pa·s;

k——称为第二粘性系数,或整体粘性系数(bulk coefficient ofviscosity);

C)、由能量守恒原理知,

控制体中单位时间多相流体总能量(动能和内能之和)的增量(即总能量对时间的变化率)等于:

(1)单位时间通过控制面流入流体的能量通量;

(2)作用在控制体和控制面上的所有外力(体力和表面力)在单位时间对流体所作的功;

(3)单位时间通过控制面传入的热通量三部分之和;表示为:

ddtV1(t)ρ1(12U12+e1)dV1+ddtV2(t)ρ1(12U22+e2)dV2=-A1(t)ρ1(12U12+e1)U1·n1dA1

-A2(t)ρ2(12U22+e2)U2·n2dA+V1(t)ρ1f·U1dV1

+V2(t)ρ2f·U2dV2+A1(t)U1·(n1·T1)dA1---(1-4)

+A2(t)U2·(n2·T2)dA2-A1(t)q1·n1dA1

-A2(t)q2·n2dA2

式中

ek——为第k相单位质量流体的内能,kJ/kg;

qk——为第k相热流矢量,即单位时间通过单位面积的热通量,kJ/(m2.s)。

所述的步骤1.1-2)根据管道中气液两相流动的描述简化步骤1.1-1)中的计算方法的步骤进一步包括:

步骤1.1-2.1)、当采用均相流模型时,

如果多相混合物能看成均匀介质,流动参数可取为各相介质的平均值,进而可将混合物平衡方程简化为均相流平衡方程,

这时要求:①气相和液相的速度相等,即vg=vl=v

②多相介质已达到热力学平衡状态,压力、密度等互为单值函数;

这条假设在等温流动中是成立的,在不等温稳定流动中是基本成立的,在不稳定流动中则是近似的;

对于等截面管道中的稳定流动此方程式可简化为:

-dpdz=τ0χA+ddz(ρhv2)+hsinθ---(1-5)

这是计算压降时最常用的方程式;

步骤1.1-2.2)、当采用分相流模型时,

分相流模型是把多相流动看成为各相分开的流动,介质参数分别取各自独立的物性参数;为此,需要分别建立各相的流体动力特性方程式;这就需要预先确定各相占有过流段面的份额(即真实含气率)以及介质与管壁的摩擦阻力和各相介质之间的摩擦阻力;这些数据利用实验找出经验关系式;

分流模型的基本假设为:

①各相介质分别有各自的平均流速,其数值根据所占的断面积计算;

②尽管各相之间可能有质量交换,但相间处于热力学平衡状态,压力、密度均为单值函数;

对于等截面圆形管道中的稳定流动,可推得压力梯度方程式为

-dpdz=τ0χA+m·2ddz[(1-x)2ρl(1-φ)+x2ρgφ]+TPsinθ---(1-6)

对分相流模型,定常气液两相流的基本方程为:

z(Σk=12ρkvkA)=0

-dpdz=τ0χA+(Σk=12ρkvk)2ddzΣk=12xk2ρkφk+gΣk=12ρkφksinθ

ddz[(Σk=12ρkvk)Ahm]+ddz{(Σk=12ρkvk)Avm22}+g(Σk=12ρkvk)Asinθ=q·χ

定解条件为:给定管道入口的压力、温度、气相流量、液相流量、管道直径、管道长度、出口压力;

式中:管道各截面的气相和液相速度vg、vl、混合物压力p、气相密度ρg、气相和液相内能eg、el,以及,气相和液相温度Tg、Tl,干度x,空隙率φ,沿管壁的热通量壁面摩擦力τ0,反映汽化潜热的焓值h1g、内热发生率qv

所述的步骤1.1-3)进一步包括以下步骤:

步骤1.1-3.1)、计算压力梯度表达式:

A)、均流模型的压力梯度表达式

根据均流模型等截面管道中的稳定流动压力梯度表达式,总压力梯度是摩阻压力梯度、重位压力梯度和加速压力梯度之和。或写成如下形式:

-dpdz=(pz)1+(pz)2+(pz)3---(1-7)

在本发明中,采用Beggs Brill(贝格斯布里尔)方法计算总压力梯度,其中包括流态判别、持液率、沿程阻力系数及压降的计算。

B)、分流模型的压力梯度表达式

如果流动为分层流,采用层流模型来计算压力梯度;

-dpdz=2f0VlD(MA)2φ02+ρTPgsinθ+(MA)2ddz[x2Vgφ+(1-x)2Vl1-φ]---(1-8)

式中

φ02——无因次参数;用Xiao(1990)等人的计算方法来确定的;

由(1-8)式可见,在求得式中φ02值及空隙率φ的相关规律后,就可求得倾斜气液两相分层流的压力梯度;

步骤1.1-3.2)、确定综合压降的计算方法:

根据上述压力梯度方程,就可以求出倾斜管多相流动的压降,其计算流程如下所示:

(A)、从入口出发,并给定一个Δp;

(B)、计算出在p、T下的物性参数,进行流动型态判别,当为

(C1)、分离流,则:利用分流模型压差公式(1-5)计算压差,

(C2)、间歇流,则:利用均流模型压差公式(1-6)计算压差,

(C3)、分散流,则:利用均流模型压差公式(1-6)计算压差,

(D)、算出ΔL,重复上述步骤一直到∑ΔL等于两个测压点之间的距离。

所述的步骤1)中超深井稠油原油物性模拟的修正模型的建立进一步包括:

步骤1.2-1)、溶解油气比采用Vazquez和Beggs的研究结果,即

Ss=0.1781C1δgs(0.1450p)C2exp{C3[δoAPI1.8t+492]}---(1-9)

式中Ss为溶解油气比,m3/m3;p为压力(绝对),kPa;t为温度,℃;δgs为689.5kPa表压下天然气的相对密度;δoAPI为标准状况下,原油的API密度;Ci为系数;

采用以下方法进行修正

Cik=Ki·Ci (i=1,2,3)     (1-10)

式中Ki为修正系数;Cik为修正后的系数值;

步骤1.2-2)、原油体积系数Bo的典型计算方法是Standing公式和Vazquez-Beggs公式:

当压力p≤pb时,这两个公式所得到的模拟曲线具有相似的变化规律,但近似相差某一常数;

当压力p>pb时,由于Vazquez-Beggs公式考虑了饱和压力pb的影响,在pb处曲线出现转折,而Standing曲线不存在转折点;

考虑饱和压力的影响,油田原油体积系数Bo的模拟计算公式表示为

Bo=BoS+BoV2+[0.0892(1-pb12)-0.0035]---(1-11)

式中BoS为Standing公式计算出的原油体积系数,m3/m3;BoV为Vazquez-Beggs公式计算出的原油体积系数,m3/m3

步骤1.2-3)、在计算稠油的粘度时,首先要根据实验室的粘温测试数据,根据Andrde粘温关系式得出稠油粘度随温度的变化关系:

μt=aeb/T    (1-12)

式中μt为Andrde公式计算出的原油粘度,mPa.s;T为测试的绝对温度,K;a,b为常数;根据室内实验测定两个温度下的粘度值代入上式就可求得常数a,b,实测值较多时,采用最小二乘法求得a和b值;

步骤1.2-4)、考虑压力对粘度的影响时常用下式表示为

μp=μt×100.15p[0.0239+0.01638(μtγ0)0.278]---(1-13)

式中p为单元段中平均的原油压力,Pa;μp为脱气原油在压力p下的粘度,mPa.s;

步骤1.2-5)、考虑溶解气对粘度的影响时,溶解气的影响表征为:

μo=AμpB    (1-14)

其中A=10.715(5.615Ss+100)-0.515,B=5.44(5.615Ss+150)-0.338

所述步骤1.2-5)中进一步包括:在实际应用中通过油田上实测的PVT数据进行修正:

将μp、B及A括号中的数据看作常数,即令C=μpB,则通过实测数据,用插值法求得C1、C2,无模拟数据时,C1、C2取理论值,即C1=10.715,C2=-0.515。

步骤2)推导出超深井稠油掺稀井井筒中流体温度梯度的计算方法具体的包括:在常规生产时井筒流体温度分布计算:通常,自喷与常规有杆泵抽油是指井筒不加热抽油,温度沿井筒的分布,可用下列公式计算:

θ=[(θ0-ted+mH)-mWKi]eKiW(H-1)+(ted-mh)+mWKi---(1-15)

式中,θ为深度h处产液的温度,℃;θ0为井口的温度,℃;m为地温梯度,℃/m;H为井深,m;W为产出液水当量(流量与比热之积);k1为油气混合物与地层间传热系数,W/(m·℃);ted为井底地层温度,℃;h为从井底向上起算的距离,m。

所述的步骤2)推导出超深井稠油掺稀井井筒中流体温度梯度的计算方法具体的包括:在油套环空掺流体工艺井筒流体温度分布计算:计算井筒温度分布时将井身结构分为两段:

第一段为热载体循环段,从井口至泵的入口处,在该段内,井筒温度分布可根据正、反循环情况进行计算;第二段为泵的入口至油层中部,由于该段流体流动近似于常规井中流体流动,因此,可按常规井进行分析计算井筒温度分布;

热流体在流向地层深处的过程中,由于传导、对流作用,一方面通过套管、水泥环,向地层传递热量,另一方面又要经过与地层产出液体混合,加热产出液体,根据传热学知识与能量平衡原理,可得出能量平衡方程式:

t=C1er1h+C2er2h+to+mh+(W-W2kh3-W2kh1)mθ=(1+W2r1kh1)C1er1h+(1+W2r2kh1)C2er2h+to+mh+W-W2kh3m---(1-16)

式中C1与C2的值分别为:

C1={[ti-to-(W-W2kh3-W2kh1)m]er2hf+to+mhf+(W-W2kt3-W2kt1)m-tf}/(er2hf-er1hf)C2={[ti-to-(W-W2kh3-W2kh1)m]er1hf+to+mhf+(W-W2kt3-W2kt1)m-tf}/(er1hf-er2hf)

式中r1与r2的值分别为:

r1=12[(kh1+kh3W-kh1W2)+(kh1+kh3W-kt1W2)2+4kh1kh3WW2]r2=12[(kh1+kh3W-kh1W2)-(kh1+kh3W-kh1W2)2+4kh1kh3WW2].

所述的步骤2)推导出超深井稠油掺稀井井筒中流体温度梯度的计算方法具体的包括:油管掺流体工艺井筒流体温度分布计算

t=(1-Wr1kh1)C1er1h+(1-Wr2kh1)C2er2h+to+mh+W-W2kh3mθ=C1er1h+C2er2h+to+mh+(W-W2kh3+Wkh1)m---(1-17)

式中C1与C2的值分别为:

C1={[ti-to-W-W2kh3m]er2hf+to+mhf+W-W2kh3m-hf}/{(1-Wr1kh1)(er2hf-er1hf)}C2={[ti-to-W-W2kh3m]er1hf+to+mhf+W-W2kh3m-tf}/((1-Wr2kh1)(er1hf-er2hf))r1,2=12[(kh1W-kh1+kh3W2)±(kh1W-kh1+kh3W2)2+4kh1kh3WW2]

式中W为地面产出混合液体的水当量,W/℃;W1为油层产出液体水当量,W/℃;W2为井筒注入液体水当量,W/℃;t为沿井深任一点处注入液体的温度,℃;θ为沿井深任一点处混合液体的温度,℃;h为由井口算起沿井筒的深度,m;hf为注入液体在井筒中掺入点深度,m;kh1为油管内流体与环空中流体之间的传热系数,W/(m℃);kh3为环空内流体与地层之间的传热系数,W/(m℃);to为地表年平均温度(即恒温层温度),℃;ti为注入液在地面时的温度,℃;tf为掺入深度处注入流体的温度,℃;θ′f为在掺入深度处油层产出流体的温度,℃;θf为在掺入深度处混合流体的温度,℃;m为地温梯度,℃/m。

所述的步骤3)具体包括建立井眼轨迹计算模型:

弯曲井筒内的压降计算是以压力梯度方程式为基础,在该式中,涉及到了管道与水平方向的夹角θ,在计算压降时,可以把弯曲井段分成若干个斜直段,而每段的倾斜角是不同的,要得到各段的平均井斜角,就需要计算出该井段两端点处的井斜角;此外,由于流体的温度是与垂深相联系的,所以还需要有计算井眼轨迹垂深的方法;

如果将井眼轨迹假设为空间螺旋线或自然曲线,则有[6]

α=αA+αB-αALB-LA(L-LA)---(1-18)

ΔH=LB-LAαB-αA(sinα-sinαA)---(1-19)

式中L为计算点的井深,m;α为计算点的井斜角,°;ΔH为计算点与上测点间的垂深增量,m。

所述的步骤4)具体包括处理掺稀动态参数生产参数:

步骤4.1)、产液量的计算:

在掺稀生产过程中,油井的真实产量应该等于井口实测产量减掉掺入稀油产量,即

Qc0=Q实测-Qx0-Qw    (1-20)

式中Qc0为油井的真实日产量,m3/d;Q实测为油井的实测日产液量,m3/d;Qx0为油井日掺入的稀油量,m3/d;Qw为日产水量,m3/d。

步骤4.2)、含水百分比的计算:

在掺稀生产过程中,根据含水百分比的定义可知,

式中fw为体积含水率;

步骤4.3)、掺稀比例的计算:

fx0=Qx0Qx0+Qc0---(1-22)

式中fx0为掺入稀油的体积分数。

步骤4.4)、掺稀井井筒中分段压力分布规律的确定:

在掺稀采油工艺中,如果将过流断面取在油井井底入口、井口和掺稀油入口,对所取控制体进行能量分析,根据能量平衡原理总可得到如下的方程,即

(井底产油入口能量-掺入前产出液的能量损失)+(掺入液入口能量-掺入前掺入液能量损失)-掺混后混合液的能量损失+泵给流体的能量=井口混合液的剩余能量                     (1-23)

由式(1-23)可推导出掺稀井井筒压力的分布规律如方程(1-24),即

pf=pk+Δp4+Δp3-Δp2+Δp1        (1-24)

式中pf为油井的流压,MPa;pk为井口油压,MPa;Δp4为掺混后混合液的压降,MPa;Δp3掺入点以下泵排出口以上油气水的压降,MPa,Δp2泵的增压,MPa,如果是自喷井,则Δp2=0;Δp1地层到泵吸入口油气水的压降,MPa。

所述的步骤5)对稠油掺稀比例及生产参数进行数据优化具体的为掺稀井的动态生产参数的优化:在油井举升工况分析过程中,当油压为已知时,可以以井底为求解点;首先根据油井的产能,做出油井的流入曲线A,即IPR曲线;若已知油压,则设定一组产量,即改变掺稀比例,通过井筒多相流计算可得一组流压,此压力与产量的关系曲线便为的油井的流出曲线B;把这两条曲线绘制在同一坐标系中,其交点b便为该系统在所给条件下在井底得到的解,即在所给条件下可获得的油井最大产量及相应的井底流压;与其相对应的掺稀比例就是最优掺稀比例。

本发明对于因原油入井后流动困难,需要采取井筒辅助降粘措施的油井均适用。井筒辅助降粘措施包括井筒掺稀降粘、掺活性水降粘及掺水降粘。

本发明中还设计了一种超深井稠油掺稀混配器,

所述混配器包括内管,外管,上、下封堵,混配器上、下接头和一组掺混出口;所述内管、外管通过上封堵和下封堵同心焊接固定形成喷嘴座;所述一组掺混出口为斜孔,开设在所述喷嘴座的内外管壁上;

所述稀油通过掺混出口喷出后,与井筒内流动流体充分混合;稀油以自由射流的方式射入,一边射入稀油,一边混合,与抽油管内稠油进行动量交换。

在具体的设计中,所述一组掺混出口为2-8个;且所述每个掺混出口包括内管斜孔,喷嘴架和喷嘴;所述内管斜孔设置在所述内管壁上,所述喷嘴设置在外管壁上,且喷嘴与井筒流体流动方向成角度;所述喷嘴与井筒流体流动方向的角度为30度~60度。所述喷嘴架为设置在所述内、外管之间的通路,用于连接所述内管斜孔和喷嘴。

在实际的应用中,所述的掺混出口为4个,所述的喷嘴架11为圆锥体,且所述喷嘴为圆锥体。

在具体的应用中,所述混配器通过所述上、下接头与油管螺纹连接,且在混配器下方固定有封隔器;

所述混配器设置在井筒内部的方式有反掺稀和正掺稀;反掺稀方式为所述混配器内管内径与油管内径相同,以保证采油管不变径;正掺稀方式为所述混配器外管外径与油管外径相同。

基于该方法及使用本申请所述的掺混器,在现场进行了五十余井次的实验,结果表明,利用该方进行的掺稀降粘举升工艺,较之以往掺稀比例确定方法相比,在其它生产条件都不改变的条件下,所需掺稀比例平均降低25.8%,可见掺稀效益明显提高。

本发明对于因原油入井后流动困难,需要采取井筒辅助降粘措施的油井均适用。

井筒辅助降粘措施包括井筒掺稀降粘、掺活性水降粘及掺水降粘。

随着油田开发进程的不断深入及快速发展的经济对原油的需求,稠油的开采份额越来越大,国内外需采取井筒辅助降粘措施的油井会越来越多,因此该发明在国内外均具有广阔的应用前景。

为了能更进一步了解本发明的特征以及技术内容,以下结合本发明的具体实施例及附图进行说明。但所举附图及实施例并非用来对本发明加以限制。

【附图说明】

图1是本发明方法中的控制体与控制面的示意图。

图2是本发明方法中的压降计算流程图。

图3是本发明方法的掺稀井井筒内压力分布图。

图4是本发明方法的流压与产量之间的关系曲线。

图5是本发明方法的掺稀井产量与压力之间的关系图。

图6是本发明方法的最优掺稀比例的确定示意图。

图7是本发明掺稀混配器的正视图的结构示意图。

图8是本发明掺稀混配器的左视图的结构示意图。

图9是本发明掺稀混配器的主视图的结构示意图。

图10是本发明与掺稀混配器连接的掺稀井筒结构示意图的反掺方式示意图。

图11是本发明与掺稀混配器连接的掺稀井筒结构示意图的正掺方式示意图。

【具体实施方式】

本发明提出了超深井变流变模式的多相流综合压降计算模型和超深井稠油原油物性模拟的修正模型,其结果与PVT曲线具有良好的相关性。根据掺稀方式给出了超深井稠油掺稀井井筒中流体温度梯度的计算方法。通过引入井眼轨迹的描述与计算技术,成功地解决了超深井弯曲井段的多相流计算问题,从而可以计算出任意井眼形状中的参数分布规律。在此基础上,建立了稠油掺稀生产动态分析模型,实现了对超深井稠油掺稀生产动态的模拟和分析,并据此可对稠油掺稀比例及生产参数进行优化,为使加入的稀油能够更加均匀地掺混到稠油中去,发明出一种新型井下掺混器,现场应用收到了良好的效果。

建立多相流综合压降计算模型和超深井稠油原油物性模拟的修正模型:

首先,建立多相流综合压降的计算模型。

1、气液两相倾斜管流的基本控制方程的建立,主要介绍气液两相流动的基本理论。

1.1、气液两相流的基本控制方程

为得到两相流控制方程,我们首先考虑包括两个相的一个确定的控制体的积分平衡定律,再由该体积分导出对每一相有效的偏微分方程。

如图1所示,在多相流动的流场中考虑一个固结于惯性坐标系被界面AI(t)切割的确定的控制体V,它由相界面AI(t)分为两个分别由两相占据的子空间域V1(t)和V2(t),V1(t)由表面A1(t)和AI(t)围成,V2(t)由表面A2(t)和AI(t)围成,由A1(t)和A2(t)组成控制面。假设控制体内气液两相具有明显的相界面,各种交换只能通过相界面进行,且控制体内气液两相符合以下平衡定律。为简便起见,这里忽略表面张力的作用。

1).质量守恒

假定所取的控制体内没有源和汇。由质量守恒原理,得到两相流动的连续方程可表示为

ddtV1(t)ρ1dV1+ddtV2(t)ρ2dV2=-A1(t)ρ1U1·n1dA1-A2(t)ρ2U2·n2dA2---(1-1)

式中

ρk——为k相(k=1,2)的密度,kg/m3

Uk——为k相(k=1,2)的速度矢量,m/s;

nk——相应表面的单位外法向矢量。

2).动量守恒

控制体V中动量对时间的变化率等于:(1)单位时间通过控制面A进入的流体动量通量;(2)作用在控制体中和控制面上所有外力的合力,这些外力由体积力(即重力)和表面力(即应力张量)组成:

ddtV1(t)ρ1U1dV1+ddtV2(t)ρ1U2dV2=-A1(t)ρ1U1(U1·n1)dA1

-A2(t)ρ2U2(U2·n2)dA2+V1(t)ρ1fdV1+V2(t)ρ2fdV2(1-2)

+A1(t)n1·T1dA1+A2(t)n2·T2dA2

式中

f——是作用在流体上的单位质量力,N;

Tk——为第k相流体在控制面上的应力张量,它与应变率之间的关系由流体的本构方程给出。通常,应力张量可表示为:

T=-pI+τ

式中

p——为压力,Pa;

I——为单位张量;

τ——为粘性应力张量。

对牛顿流体:

τik=μ(vixk+vkxi-23δikvlxl)+κδikvlxl---(1-3)

式中

μ——为剪切粘性系数,Pa.s;

k——称为第二粘性系数,或整体粘性系数(bulk coefficient ofviscosity)。

3).能量守恒

由能量守恒原理知,控制体中单位时间多相流体总能量(动能和内能之和,这里未计及势能)的增量(即总能量对时间的变化率)等于:(1)单位时间通过控制面流入流体的能量通量;(2)作用在控制体和控制面上的所有外力(体力和表面力)在单位时间对流体所作的功;(3)单位时间通过控制面传入的热通量三部分之和。表示为:

ddtV1(t)ρ1(12U12+e1)dV1+ddtV2(t)ρ2(12U22+e2)dV2=-A1(t)ρ1(12U12+e1)U1·n1dA1

-A2(t)ρ2(12U22+e2)U2·n2dA+V1(t)ρ1f·U1dV1

+V2(t)ρ2f·U2dV2+A1(t)U1·(n1·T1)dA1---(1-4)

+A2(t)U2·(n2·T2)dA2-A1(t)q1·n1dA1

-A2(t)q2·n2dA2

式中

ek——为第k相单位质量流体的内能,kJ/kg;

qk——为第k相热流矢量,即单位时间通过单位面积的热通量,kJ/(m2.s)。

我们一共得到三组(气相、液相、界面)九个方程(六个平衡方程、三个界面方程)来描述气液两相流动。这三组方程描述的是气相、液相和界面在流场中的瞬时状态,就整体流场而言,它们是不连续的,不能直接通过常用的数学方法来求解。为了达到求解的目的,必须对它们进行连续化使其成为类似单相流的方程组。所谓连续化,就是用平均的方法,把方程改造成在一定尺度内既保留原流场的起伏变化特性,又表现出连续性。这里的平均,可以按时间间隔进行,称为时间平均;也可以按空间(长度、面积或体积)间隔进行,称为空间平均。不管按时间平均还是按空间平均,都必须满足前面提出的原则,即在最大尺度内需保证平均后仍保留原流场的波动特性,在最小尺度内达到必要的连续特性。

1.2、管道中气液两相流动的描述

如果仅仅研究管道中气液两相流动,上述描述多相流的方程组可以得到很大的简化。首先,若管道的长度远远大于管道直径,就可认为管道横截面的流动参数的平均值不变,即流动参数仅在轴线方向变化。其次,在考虑流动结构后,可将流场分为若干区域,在各个区域内,将流动处理为单相的或多相混合的。例如当流型已知且各相分离较好时(如分层流或环状流),采用这种方法就比较简便。这样就可以得到不同的描述流动的模型。

为了方便,实际使用上述流动模型时,还往往根据流动性质的不同,将其划分为均相流模型和分相流模型。下面将给出具体的表达式。

1).均相流模型

如果多相混合物能看成均匀介质,流动参数可取为各相介质的平均值,进而可将混合物平衡方程简化为均相流平衡方程。

这时要求:①气相和液相的速度相等,即vg=vl=v。

②多相介质已达到热力学平衡状态,压力、密度等互为单值函数。

这条假设在等温流动中是成立的,在不等温稳定流动中是基本成立的,在不稳定流动中则是近似的。

对于等截面管道中的稳定流动此方程式可简化为:

-dpdz=τ0χA+ddz(ρhv2)+hsinθ---(1-5)

这是计算压降时最常用的方程式。

2).分相流模型

分相流模型是把多相流动看成为各相分开的流动,介质参数分别取各自独立的物性参数。为此,需要分别建立各相的流体动力特性方程式。这就需要预先确定各相占有过流段面的份额(即真实含气率)以及介质与管壁的摩擦阻力和各相介质之间的摩擦阻力。目前,这些数据主要是利用实验找出经验关系式。

基本假设:

分流模型的基本假设为:

①各相介质分别有各自的平均流速,其数值根据所占的断面积计算;

②尽管各相之间可能有质量交换,但相间处于热力学平衡状态,压力、密度等均为单值函数。

以环状流动为例对分流模型进行了说明;但是,分流模型的基本方程式与流动所具有的特殊流动结构无关。

对于等截面圆形管道中的稳定流动,可推得压力梯度方程式为

-dpdz=τ0χA+m·2ddz[(1-x)2ρl(1-φ)+x2ρgφ]+TPsinθ---(1-6)

对分相流模型,定常气液两相流的基本方程为:

z(Σk=12ρkvkA)=0

-dpdz=τ0χA+(Σk=12ρkvk)2ddzΣk=12xk2ρkφk+gΣk=12ρkφksinθ

ddz[(Σk=12ρkvk)Ahm]+ddz{(Σk=12ρkvk)Avm22}+g(Σk=12ρkvk)Asinθ=q·χ

定解条件为:给定管道入口的压力、温度、气相流量、液相流量、管道直径、管道长度、出口压力。

需要求解的未知量有:管道各截面的气相和液相速度vg、vl、混合物压力p、气相密度ρg、气相和液相内能eg、el,相关未知量有:气相和液相温度Tg、Tl,干度x,空隙率φ,沿管壁的热通量壁面摩擦力τ0,反映汽化潜热的焓值h1g、内热发生率qv

一共14个未知量,已经有六个方程,需要补充与均相流模型类似的8个方程,且为本领域技术人员能知悉的,故,这里不再赘述。

这样,气液两相流分相流模型也是一个封闭的定解问题。

1.3、倾斜管多相流压降的计算方法

在石油工程中,倾斜管多相流的工艺计算主要包括流型判别、持液率计算和压降计算。流型判别和持液率计算的目的是得到尽可能准确的压降分布。在采油过程中,准确的压降数据对于优化设计采油方案、监测油井的生产动态和改善油田的开发效果具有重要意义。在稠油注蒸汽热采过程中,压降计算是进行注汽设计的基础,对油井注汽工况的分析具有重要意义。不仅可以进行注汽参数敏感性分析,优选注汽参数、合理设计汽量,而且可为井下分层配汽器的设计提供必要的理论基础。本文在压降计算中,给出一种综合压降计算模型,即针对不同流动型态的特点,分别采用以下两种压降计算方法。

1.3.1、压力梯度表达式

1).均流模型的压力梯度表达式

根据均流模型等截面管道中的稳定流动压力梯度表达式,总压力梯度是摩阻压力梯度、重位压力梯度和加速压力梯度之和。或写成如下形式:

-dpdz=(pz)1+(pz)2+(pz)3---(1-7)

在本发明中,采用Beggs Brill(贝格斯布里尔)方法计算总压力梯度,其中包括流态判别、持液率、沿程阻力系数及压降的计算。

B)、分流模型的压力梯度表达式

如果流动为分层流,采用层流模型来计算压力梯度;

-dpdz=2f0VlD(MA)2φ02+ρTPgsinθ+(MA)2ddz[x2Vgφ+(1-x)2Vl1-φ]---(1-8)

式中

φ02——无因次参数;用Xiao(1990)等人的计算方法来确定的;

由(1-8)式可见,在求得式中φ02值及空隙率φ的相关规律后,就可求得倾斜气液两相分层流的压力梯度;

步骤1.1-3.2)、确定综合压降的计算方法:

根据上述压力梯度方程,就可以求出倾斜管多相流动的压降,其计算流程如下所示:

(A)、从入口出发,并给定一个Δp;

(B)、计算出在p、T下的物性参数,进行流动型态判别,当为

(C1)、分离流,则:利用分流模型压差公式(1-5)计算压差,

(C2)、间歇流,则:利用均流模型压差公式(1-6)计算压差,

(C3)、分散流,则:利用均流模型压差公式(1-6)计算压差,

(D)、算出ΔL,重复上述步骤一直到∑ΔL等于两个测压点之间的距离。

步骤1)中超深井稠油原油物性模拟的修正模型的建立进一步包括:

步骤1.2-1)、溶解油气比采用Vazquez和Beggs的研究结果,即

Ss=0.1781C1δgs(0.1450p)C2exp{C3[δoAPI1.8t+492]}---(1-9)

式中Ss为溶解油气比,m3/m3;p为压力(绝对),kPa;t为温度,℃;δgs为689.5kPa表压下天然气的相对密度;δoAPI为标准状况下,原油的API密度;Ci为系数;采用以下方法进行修正

Cik=Ki·Ci  (i=1,2,3)       (1-10)

式中Ki为修正系数;Cik为修正后的系数值;

步骤1.2-2)、原油体积系数Bo的典型计算方法是Standing公式和Vazquez-Beggs公式:

当压力p≤pb时,这两个公式所得到的模拟曲线具有相似的变化规律,但近似相差某一常数;

当压力p>pb时,由于Vazquez-Beggs公式考虑了饱和压力pb的影响,在pb处曲线出现转折,而Standing曲线不存在转折点;

考虑饱和压力的影响,油田原油体积系数Bo的模拟计算公式表示为

Bo=BoS+BoV2+[0.0892(1-pb12)-0.0035]---(1-11)

式中BoS为Standing公式计算出的原油体积系数,m3/m3;BoV为Vazquez-Beggs公式计算出的原油体积系数,m3/m3

步骤1.2-3)、在计算稠油的粘度时,首先要根据实验室的粘温测试数据,根据Andrde粘温关系式得出稠油粘度随温度的变化关系:

μt=aeb/T    (1-12)

式中μt为Andrde公式计算出的原油粘度,mPa.s;T为测试的绝对温度,K;a,b为常数;根据室内实验测定两个温度下的粘度值代入上式就可求得常数a,b,实测值较多时,采用最小二乘法求得a和b值;

步骤1.2-4)、考虑压力对粘度的影响时常用下式表示为

μp=μt×100.15p[0.0239+0.01638(μtγ0)0.278]---(1-13)

式中p为单元段中平均的原油压力,Pa;μp为脱气原油在压力p下的粘度,mPa.s;

步骤1.2-5)、考虑溶解气对粘度的影响时,溶解气的影响表征为:

μo=AμpB    (1-14)

其中A=10.715(5.615Ss+100)-0.515,B=5.44(5.615Ss+150)-0.338

步骤1.2-5)中进一步包括:在实际应用中通过油田上实测的PVT数据进行修正:

将μp、B及A括号中的数据看作常数,即令C=μpB,则通过实测数据,用插值法求得C1、C2,无模拟数据时,C1、C2取理论值,即C1=10.715,C2=-0.515。

步骤2)推导出超深井稠油掺稀井井筒中流体温度梯度的计算方法具体的包括:在常规生产时井筒流体温度分布计算:通常,自喷与常规有杆泵抽油是指井筒不加热抽油,温度沿井筒的分布,可用下列公式计算:

θ=[(θ0-ted+mH)-mWKt]eKtW(H-1)+(ted-mh)+mWKt---(1-15)

式中,θ为深度h处产液的温度,℃;θ0为井口的温度,℃;m为地温梯度,℃/m;H为井深,m;W为产出液水当量(流量与比热之积);k1为油气混合物与地层间传热系数,W/(m·℃);ted为井底地层温度,℃;h为从井底向上起算的距离,m。

步骤2)推导出超深井稠油掺稀井井筒中流体温度梯度的计算方法具体的包括:在油套环空掺流体工艺井筒流体温度分布计算:计算井筒温度分布时将井身结构分为两段:

第一段为热载体循环段,从井口至泵的入口处,在该段内,井筒温度分布可根据正、反循环情况进行计算;第二段为泵的入口至油层中部,由于该段流体流动近似于常规井中流体流动,因此,可按常规井进行分析计算井筒温度分布;

热流体在流向地层深处的过程中,由于传导、对流作用,一方面通过套管、水泥环,向地层传递热量,另一方面又要经过与地层产出液体混合,加热产出液体,根据传热学知识与能量平衡原理,可得出能量平衡方程式:

t=C1er1h+C2er2h+to+mh+(W-W2kh3-W2kh1)mθ=(1+W2r1kh1)C1er1h+(1+W2r2kh1)C2er2h+to+mh+W-W2kh3m---(1-16)

式中C1与C2的值分别为:

C1={[ti-to-(W-W2kh3-W2kh1)m]er2hf+to+mhf+(W-W2kt3-W2kt1)m-tf}/(er2hf-er1hf)C2={[ti-to-(W-W2kh3-W2kh1)m]er1hf+to+mhf+(W-W2kt3-W2kt1)m-tf}/(er1hf-er2hf)

式中r1与r2的值分别为:

r1=12[(kh1+kh3W-kh1W2)+(kh1+kh3W-kt1W2)2+4kh1kh3WW2]r2=12[(kh1+kh3W-kh1W2)-(kh1+kh3W-kh1W2)2+4kh1kh3WW2].

步骤2)推导出超深井稠油掺稀井井筒中流体温度梯度的计算方法具体的包括:油管掺流体工艺井筒流体温度分布计算

t=(1-Wr1kh1)C1er1h+(1-Wr2kh1)C2er2h+to+mh+W-W2kh3mθ=C1er1h+C2er2h+to+mh+(W-W2kh3+Wkh1)m---(1-17)

式中C1与C2的值分别为:

C1={[ti-to-W-W2kh3m]er2hf+to+mhf+W-W2kh3m-hf}/{(1-Wr1kh1)(er2hf-er1hf)}C2={[ti-to-W-W2kh3m]er1hf+to+mhf+W-W2kh3m-tf}/((1-Wr2kh1)(er1hf-er2hf))r1,2=12[(kh1W-kh1+kh3W2)±(kh1W-kh1+kh3W2)2+4kh1kh3WW2]

式中W为地面产出混合液体的水当量,W/℃;W1为油层产出液体水当量,W/℃;W2为井筒注入液体水当量,W/℃;t为沿井深任一点处注入液体的温度,℃;θ为沿井深任一点处混合液体的温度,℃;h为由井口算起沿井筒的深度,m;hf为注入液体在井筒中掺入点深度,m;kh1为油管内流体与环空中流体之间的传热系数,W/(m℃);kh3为环空内流体与地层之间的传热系数,W/(m℃);to为地表年平均温度(即恒温层温度),℃;ti为注入液在地面时的温度,℃;tf为掺入深度处注入流体的温度,℃;θ′f为在掺入深度处油层产出流体的温度,℃;θf为在掺入深度处混合流体的温度,℃;m为地温梯度,℃/m。

步骤3)具体包括建立井眼轨迹计算模型:

弯曲井筒内的压降计算是以压力梯度方程式为基础,在该式中,涉及到了管道与水平方向的夹角θ,在计算压降时,可以把弯曲井段分成若干个斜直段,而每段的倾斜角是不同的,要得到各段的平均井斜角,就需要计算出该井段两端点处的井斜角;此外,由于流体的温度是与垂深相联系的,所以还需要有计算井眼轨迹垂深的方法;如果将井眼轨迹假设为空间螺旋线或自然曲线,则有

α=αA+αB-αALB-LA(L-LA)---(1-18)

ΔH=LB-LAαB-αA(sinα-sinαA)---(1-19)

式中L为计算点的井深,m;α为计算点的井斜角,°;ΔH为计算点与上测点间的垂深增量,m。

步骤4)具体包括处理掺稀动态参数生产参数:

步骤4.1)、产液量的计算:

在掺稀生产过程中,油井的真实产量应该等于井口实测产量减掉掺入稀油产量,即

Qc0=Q实测-Qx0-Qw    (1-20)

式中Qc0为油井的真实日产量,m3/d;Q实测为油井的实测日产液量,m3/d;Qx0为油井日掺入的稀油量,m3/d;Qw为日产水量,m3/d。

步骤4.2)、含水百分比的计算:

在掺稀生产过程中,根据含水百分比的定义可知,

式中fw为体积含水率;

步骤4.3)、掺稀比例的计算:

fx0=Qx0Qx0+Qc0---(1-22)

式中fx0为掺入稀油的体积分数。

步骤4.4)、掺稀井井筒中分段压力分布规律的确定:

在掺稀采油工艺中,如果将过流断面取在油井井底入口、井口和掺稀油入口,对所取控制体进行能量分析,根据能量平衡原理总可得到如下的方程,即

(井底产油入口能量-掺入前产出液的能量损失)+(掺入液入口能量-掺入前掺入液能量损失)-掺混后混合液的能量损失+泵给流体的能量=井口混合液的剩余能量                                                              (1-23)

由式(1-23)可推导出掺稀井井筒压力的分布规律如方程(1-24),即

pf=pk+Δp4+Δp3-Δp2+Δp1       (1-24)

式中pf为油井的流压,MPa;pk为井口油压,MPa;Δp4为掺混后混合液的压降,MPa;Δp3掺入点以下泵排出口以上油气水的压降,MPa,Δp2泵的增压,MPa,如果是自喷井,则Δp2=0;Δp1地层到泵吸入口油气水的压降,MPa。

根据上述的原油物性模拟以及井筒内的温度分布规律、密度分布规律、粘度分布规律等理论,可以求解压力梯度方程式,直到油层为止,从而得到各参数沿井筒的分布规律。但是,上述计算过程是在已知井筒内温度梯度的前提下进行的,而温度梯度又与掺稀有关,所以需要使用迭代法求解。

如果井口油压和井底流压都是已知的,则可以通过计算油井的物性参数,对掺稀井的实际工作特性进行分析。

数据优化:对稠油掺稀比例及生产参数进行数据优化。

根据能量平衡原理并结合本文所给掺稀井井筒中分段压降的计算方法,可对掺稀比例、掺稀后的降粘率、掺入稀油的比重及产液的含水、掺入深度等进行参数敏感性分析,从而搞清各参数对掺稀井动态的影响,为稠油掺稀降粘方案的制定提供理论依据。另外,可通过对掺稀井动态生产的模拟结果,对油井的掺稀比例进行优化。

在油井举升工况分析过程中,当油压为已知时,可以以井底为求解点。首先根据油井的产能,做出油井的流入曲线A,即IPR曲线。若已知油压,则设定一组产量,即改变掺稀比例,通过井筒多相流计算可得一组流压,此压力与产量的关系曲线便为的油井的流出曲线B。把这两条曲线绘制在同一坐标系中,其交点b便为该系统在所给条件下在井底得到的解,即在所给条件下可获得的油井最大产量及相应的井底流压(如图4所示)。与其相对应的掺稀比例就是最优掺稀比例。

基于上述发明内容,可研制开发相关稠油掺稀生产动态分析与优化设计模拟系统,并在现场进行了50余井次的实验和应用,收到了良好的效果。现以某油田的一口环空掺稀自喷井XB-1井为例,其基础数据如表1所示,计算结果如图5,6所示。该井的最佳掺稀比为1.15,与现场试验结果吻合较好。

应用本发明的方法,得出:

(1)流动型态是影响多相流动能量损失的重要因素。不同的流态不但影响两相流动的力学关系,而且还影响其传热和传质性质。当产气量较低时,井筒中的流态多为段塞流和泡状流。

(2)原油粘度与溶解油气比具有很好的相关性,而压力、温度等因素对原油粘度的影响可认为已隐含在溶解油气比的计算中,本方法给出了超深井稠油原油物性模拟的修正模型,其结果与PVT曲线具有良好的相关性。

(3)本方法把井眼轨迹的描述与计算技术溶入多相管流的计算中,成功地解决了超深井弯曲井段的多相流计算问题,从而可以计算出任意井眼形状中的参数分布规律。

(4)水对井筒中油、气流动的总影响还不甚清楚,目前多按油、水混合物的体积含水率来处理水的影响。但是,各相及混相密度的计算必须考虑随压力、温度等的变化规律。

(5)根据能量平衡原理及多相流体力学原理建立了稠油掺稀生产动态分析模型,实现了对超深井稠油掺稀生产动态的模拟和分析,并据此可对稠油掺稀比例及生产参数进行优化,现场应用收到了良好的效果。为深层稠油掺稀井的掺稀生产方案设计提供了科学的方法和依据。

使用本发明的方法采用一种掺稀混配器,如图7、8、9所示,该掺稀混配器包括:外管1、上封堵2、内管3、混配器上接头4、第一掺混孔5、第二掺混孔6、第四掺混孔7、第三掺混孔8、混配器下接头9、下封堵10,喷嘴架11。

掺混器结构为:内管3、外管1通过上封堵2、下封堵10同心焊接在一起,统称喷嘴座;喷嘴座内外壁上开有4个均匀分布的斜孔(第一掺混孔5、第二掺混孔6、第四掺混孔7、第三掺混孔8、),圆锥形的喷嘴架11通过内外壁的小孔安装在喷嘴座上,圆锥形的喷嘴与井筒流体流动方向成30~60°。喷嘴通过焊接安装在喷嘴架11上,喷嘴形状为圆锥形。稀油通过喷嘴喷出后,与井筒内流动流体充分混合。边界是抽油管,稀油以自由射流的方式射入,一边射入稀油,一边混合,与抽油管内稠油进行动量交换。

使用掺稀混配器时,井下掺稀管柱的下部需有封隔器,其特点是掺入压力不直接作用于油层,不会产生倒灌现象。掺稀方法分为套管掺稀,也叫反掺;和油管掺稀,也叫正掺两种。对于稠油自喷井,采用混配器掺稀,有正掺、反掺两种工艺形式,其井筒结构示意图如图10、11所示。其中10表示:掺入稀油,20表示:地层产液,30表示:封隔器,40表示:掺混器,50表示:混合产液,60表示:油管外径。掺混器通过螺纹与油管连接,在反掺时掺混器内管内径与油管内径相同,以保证采油管不变径;在正掺时,掺混器外管外径与油管外径相同,以保证油套环形空间不变径,最终目的减少举升井筒内的局部压力损失。

表1 XB-1井基础数据

  油层中深[m]  5522.67  原油粘度[mPa.s]  19000  气相对密度  0.81  油层中温[m]  128.03  原油粘温系数A  4.3075  掺入点深度[m]  5449.65  套管尺寸[mm]  177.8  原油粘温系数B  0.513  含水率[%]  0.15  油层静压[MPa]  61.81  掺液粘温系数A  3.561  掺液相对密度  0.85  饱和压力[MPa]  16  掺液粘温系数B  0.4274  掺液粘度[mPa.s]  50  原油相对密度  0.991  生产气油比  22  掺液温度[C]  25

综上所述,对于本领域的普通技术人员来说,可以根据本发明的技术方案和技术构思作出其他各种相应的改变和变形,而所有这些改变和变形都应属于本发明所述的权利要求的保护范围。

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号