计算井筒多相流的BeggsBrill方法.docx

计算井筒多相流的BeggsBrill方法.docx

《计算井筒多相流的BeggsBrill方法.docx》由会员分享，可在线阅读，更多相关《计算井筒多相流的BeggsBrill方法.docx（30页珍藏版）》请在冰豆网上搜索。

计算井筒多相流的BeggsBrill方法

作业一：

计算井筒多相流的Beggs-Brill方法

基本参数：

套管尺寸：

51/2’’,油管：

21/2’’,GOR=24

/

=0.9659,

=100cp,井底温度63℃，井口温度25℃，井口压力1MPa，产液137方/天，油层深度1490m，井底流压8MPa，含水率29％，泵的单级扬程4.6m。

设计步骤：

截图：

SubBeggsBrill（calAsSingle,rwAsSingle,roAsSingle,rgAsSingle,yco2AsSingle,yh2sAsSingle,yn2AsSingle,ynaclAsSingle,_ppcAsSingle,tpcAsSingle,d2AsSingle,d1AsSingle,erasAsSingle,lonAsSingle,alfa0AsSingle,T1AsSingle,pcal1AsSingle,T2AsSingle,qwAsSingle,qoAsSingle,QscAsSingle,daltalAsSingle,PCAL2AsSingle）

Do

Iflsum+daltal>lonThen

daltal=lon-lsum

EndIf

lmid=lsum+daltal/2

T=T1+（T2-T1）*lmid/lon

Do

daltap0=daltap

p=p1+Sgn（cal）*daltap/2'p利用假设的daltap计算，带入以后的p

lsum=lsum+daltal

p2=p1+Sgn（cal）*daltap

p1=p2

LoopUntillsum>=lon

PCAL2=p2

EndSub

psc=0.101325

tsc=293

densc=1.204*rg

deno=1000*ro

denw=1000*rw

qow=qo+qw

Ifqo>0Then

gor=Qsc/qo

CallPVTO（ro,rg,gor,T,p,pbo,rso,Bo,uo,sto）

EndIf

Ifqw>0Then

gwr=Qsc/qw

CallPVTW（T,p,ynacl,gwr,rsw,bw,uw,stw）

EndIf

Ifqow>0Then

GLR=Qsc/qow

RS=（rso*qo+rsw*qw）/qow

IfRS>=GLRThen

RS=GLR

EndIf

qgsc=qow*（GLR-RS）

Else

qgsc=Qsc

EndIf

Ifqgsc>0Then

ppr=p/ppc

tpr=T/tpc

Calldprzf（ppr,tpr,dr,z,cpr）

bg=psc/tsc*z*T/p

Calllgeug（rg,yn2,yco2,yh2s,p,T,z,ug）

qg=qgsc*bg

wg=densc*qgsc

densg=wg/qg

Else

qg=0

densg=0

bg=0

EndIf

Ql=qo*Bo+qw*bw

wl=deno*qo+denw*qw+densc*RS*qow

IfQl>0Then

densl=wl/Ql

fo=qo*Bo/Ql

EndIf

ul=uo*fo+uw*（1-fo）

stal=sto*fo+stw*（1-fo）

EndSub

gor=gorin*5.615

T=tin*1.8-460

p=145*pin

api=141.5/ro-131.5

Ifapi<=30Then

a1=0.0362

b1=1.0937

c1=25.724

Else

a1=0.0178

b1=1.187

c1=23.931

EndIf

tp=T+460

pbo=（gor/（a1*rg*Exp（c1*api/tp）））^（1/b1）

rso=a1*rg*p^b1*Exp（c1*api/tp）'ft^3/bbl

Ifp>=pboOrrso>=gorThen

rso=gor

EndIf

Ifapi<=30Then

a1=1.761*10^（-5）

b1=0.0004677

c1=-1.1811*10^（-8）

Else

a1=0.000011

b1=0.000467

c1=1.377*10^（-9）

EndIf

tp=T-60

Bo=1+a1*tp*（api/rg）+（b1+c1*tp*（api/rg））*rso

co=（-1433+5*rso+17.2*T-1180*rg+12.61*api）/1000001/p

Ifp>=pboOrrso>=gorThenBo=Bo*Exp（co*（pbo-p））

zz=3.0324-0.0203*api

yy=10^zz

xx=yy*T^（-1.163）

uod=10^xx-1

aa=10.715*（rso+100）^（-0.515）

bb=5.44*（rso+150）^（-0.338）

'uo=aa*uod^bb

uo=100'

Ifpbo>0And（p>=pboOrrso>=gor）Then

uo=uo*（p/pbo）^（2.6*p^1.187*Exp（-0.0000898*p-11.513））

EndIf

sto=（42.4-0.047*T-0.267*api）*Exp（-0.0007*p）/1000

pbo=pbo/145

rso=rso/5.615

EndSub

T=tin*1.8-460'F

p=145*pin

a1=2.12+0.00345*T-0.0000359*T^2

b1=0.0107-0.0000526*T+1.48*10^（-7）*T^2

c1=-8.75*10^（-7）+3.9*10^（-9）*T-1.02*10^（-11）*T^2

scnacl=1-（0.0753-0.000173*T）*ynacl

rws=（a1+b1*p+c1*p^2）*scnacl/5.615

Ifrws>gwrThen

rws=gwr

EndIf

Ifrws<=0Then

a1=0.9947+0.0000058*T+1.02*10^（-6）*T^2

b1=-4.228*10^（-6）+1.8376*10^（-8）*T-6.77*10^（-11）*T^2

c1=1.3*10^（-10）-1.3855*10^（-12）*T+4.285*10^（-15）*T^2

Else

a1=0.9911+6.350001*10^（-5）*T+8.5*10^（-7）*T^2

b1=-1.093*10^（-6）-3.497*10^（-9）*T+4.75*10^（-12）*T^2

c1=-5*10^（-11）+6.429*10^（-13）*T-1.43*10^（-15）*T^2

EndIf

tp=T-60

tl=5.1*10^（-8）*p+（5.47*10^（-6）-1.95*10^（-10）*p）*tp

scfm=（tl+（-3.23*10^（-8）+8.5*10^（-13）*p）*tp^2）*ynacl+1

bw=（a1+b1*p+c1*p^2）*scfm

tl=（T^0.5-0.0135*T）*（0.00276*ynacl-0.000344*ynacl^1.5）

scvis=1-0.00187*ynacl^0.5+0.000218*ynacl^2.5+tl

scpr=1+3.5*10^（-12）*p^2*（T-40）

uw=0.02414*scvis*scpr*10^（446/（T+208））

ax=1-6.94829*10^（-4）*p+8.878355*10^（-7）*p^2

ax=ax+2.36093*10^（-13）*p^4-6.42967*10^（-10）*p^3

ax=ax-4.19405*10^（-17）*p^5+2.86787*10^（-21）*p^6

stw=ax*（0.1746-0.00021*T）/2.2046

EndSub

'程序功能：

计算当地条件下的有关参数的过程

ga=9.807

don=d2-d1'

aer=3.14159*（d2^2-d1^2）/4

Vsl=Ql/86400/aer'm/s’

vsg=qg/86400/aer'

Vm=Vsl+vsg'

na=Vsl/Vm'

nr=Vm^2/（ga*don）'

IfVsl>0Then

nv=Vsl*（densl/（ga*stal））^0.25'

EndIf

denn=na*densl+（1-na）*densg'

un=（na*ul+（1-na）*ug）/1000

ren=denn*Vm*don/un'

fmn=（1/（1.14-2*Log（eras/don+21.25/ren^0.9）/Log（10）））^2'

EndSub

Ifna>0Then

ls1=316*na^0.302

ls2=0.0009252*na^（-2.4684）

ls3=0.1*na^（-1.4516）

ls4=0.5*na^（-6.738）

Ifna<0.00001Then

ident=5

ElseIfna>0.99999Then'

ident=6

ElseIfna<0.01Andnr=0.01Andnr

ident=1

ElseIf（na>=0.01Andna<0.4）And（nr>ls3Andnr<=ls1）Orna>=0.4And（nr>ls3Andnr<=ls4）Then

ident=2

ElseIfna<0.4Andnr>=ls1Orna>=0.4Andnr>=ls4Then

ident=3

ElseIfna>=0.01And（nr>=ls2Andnr<=ls3）Then

ident=4

EndIf

Else

ident=5

EndIf

EndSub

Ifalfa>=0Then

SelectCaseident

Case1

dd=0.011

ee=-3.768

ff=3.539

gg=-1.614

Case2

dd=2.96

ee=0.305

ff=-0.4473

gg=0.0978

Case3

dd=1

ee=0

ff=0

gg=0

EndSelect

Else

dd=4.7

ee=-0.3692

ff=0.1244

gg=-0.5056

EndIf

cc=（1-na）*Log（dd*na^ee*nv^ff*nr^gg）'

Ifcc<0Then

cc=0

EndIf

tm=Sin（1.8*alfa）

hui=1+cc*（tm-tm^3/3）'

Ifhui<0Then

hui=0

EndIf

SelectCaseident

Case1

dd=0.98

ee=0.4846

cc=0.0868

Case2

dd=0.845

ee=0.5351

cc=0.0173

Case3

dd=1.065

ee=0.5824

cc=0.0609

EndSelect

hl0=dd*na^ee*nr^（-cc）'

Ifhl0

hl0=na

EndIf

hlsta=0.918*hl0*hui

Ifhlsta<0Then

hlsta=0

EndIf

Ifhlsta>1Then

hlsta=1

EndIf

EndSub

ga=9.807

SelectCaseident

Case4

ls2=0.0009252*na^（-2.4684）

ls3=0.1*na^（-1.4516）

ident=1

CallHOLDUP（alfa,ident,na,nv,nr,hlsta）

hlsta1=hlsta

ident=2

CallHOLDUP（alfa,ident,na,nv,nr,hlsta）

hlsta2=hlsta

tm=（ls3-nr）/（ls3-ls2）

hlsta=tm*hlsta1+（1-tm）*hlsta2

Case5

hlsta=0

Case6

hlsta=1

CaseElse

CallHOLDUP（alfa,ident,na,nv,nr,hlsta）

EndSelect

Ifhlsta

hlsta=na

EndIf

denm=densl*hlsta+densg*（1-hlsta）

tm=denm*ga*Sin（alfa）

Ifhlsta=0Orhlsta=1Then

sb=0

Else

yb=na/hlsta^2

yr=Log（yb）

Ifyb>1Andyb<1.2Then

sb=Log（2.2*yb-1.2）

Else

sb=yr/（-0.0523+3.182*yr-0.8725*yr^2+0.01853*yr^4）'=======16

EndIf

EndIf

fm=fmn*Exp（sb）

Tf=fm*denn*Vm^2/（2*don）

coef=1-denm*Vm*vsg/（1000000*p）

T1=（tm+Tf）/（1000000*coef）

daltap=daltal*T1

EndSub

a1=0.31506237

a2=-1.0467099

a3=-0.57832729

a4=0.53530771

a5=-0.61232032

a6=-0.10488813

a7=0.681570001

a8=0.68446549

a=a1+a2/tr+a3/tr^3

B=a4+a5/tr

C=a5*a6*tr

D=a7/tr^3

e=a8

G=-0.27*pr/tr

dr=-G

Do

ed=e*dr^2

exa=Exp（-ed）

f=G/dr+1+a*dr+B*dr^2+C*dr^5

f=f+D*dr^2*（1+ed）*exa

fp=-G/dr^2+a+2*B*dr+5*C*dr^4

fp=fp+2*D*dr*（1+ed-ed^2）*exa

dr=dr-f/fp

Ifdr=0Ordr=1Then

ExitSub

EndIf

LoopUntilAbs（f/fp）<0.00001

z=0.27*pr/（dr*tr）

ed=e*dr^2

exa=Exp（-ed）

dzdr=a+2*B*dr+5*C*dr^4

dzdr=dzdr+2*D*dr*（1+ed-ed^2）*exa

cr=1/pr-0.27/（z^2*tr）*dzdr/（1+dr*dzdr/z）

EndSub

mg=28.97*rg

R=0.008314

g1=yh2s*（5.7*rg-1.7）*0.00001

g2=yco2*（5*rg+1.7）*0.00001

g3=yn2*（5*rg+4.7）*0.00001

tr=1.8*T

G=（9.399999+0.02*mg）*tr^1.5

G=0.0001*G/（209+19*mg+tr）+g1+g2+g3

j=1.7-197.2/tr-0.002*mg

h=（3.5+986/tr+0.01*mg）*0.001^j

den=mg*p/（z*R*T）

ug=G*Exp（h*den^j）

EndSub

作业二抽油杆柱设计

基本参数配产：

30m3/d,泵深：

2300，动液面：

950m，冲程：

3.6，冲次6次/min，泵径：

56mm，油管内径：

73mm，地层水粘度：

1,8mPa.s，地层水密度：

1050kg/m3，抽油杆密度：

7850kg/m3，最小杆径：

19mm，杆使用系数：

0.95

设计步骤：

截图：

DimhcAsSingle'沉没度，m

DimSFAsSingle'杆使用系数

DimsAsSingle'冲程，m

DimqlAsSingle'产液量，m^3/d

DimdpAsSingle'泵径，mm

DimdeAsSingle'柱塞与泵筒间隙平均值，m

DimuwAsSingle'地层水粘度，mPa.s

DimnAsSingle'冲次，次/min

DimrouwAsSingle'地层水密度，kg/m^3

DimrousAsSingle'钢的密度，kg/m^3

DimtAsSingle'抽油杆最小抗张强度，MPa

DimPlAsSingle,PLiAsSingle,deltaPLAsSingle'输入应力范围比、计算应力范围比、两级抽油杆顶端面应力范围比的最大允许差值

DimdtAsSingle'油管内径，mm

DimfpAsSingle,frAsSingle,foAsSingle'柱塞截面积，抽油杆截面积，阀孔截面积，m^2

DimnrodAsInteger'杆级数

DimL0AsSingle,L1AsSingle'计算所得杆柱长度，m已设计出的杆柱总长度，m

DimdfoAsSingle'阀孔径，m

Dimqr（）AsSingle'每米抽油杆质量，kg/m

DimLi（）AsSingle'每级杆长度，m

Dimdr（）AsSingle'抽油杆直径,mm

DimbAsSingle'失重系数

DimmAsSingle'油管内径与抽油杆直径之比

DimvmaxAsSingle'抽油杆柱最大下行速度，m/s

DimvpAsSingle,vfAsSingle'柱塞最大运动速度、阀孔液流速度，m/s

DimNReAsSingle,ksiAsSingle'雷诺数，阀流量系数

DimPmax0AsSingle,Pmin0AsSingle,PL0AsSingle'抽油杆上下冲程底端载荷，N，底端应力范围比

Dimα（）AsSingle,dL（）AsSingle,dα（）AsSingle,dφ（）AsSingle'井斜角，井斜深、井斜角、方位角变化量

DimPmaxi1（）AsSingle,Pmini1（）AsSingle,PLi1（）AsSingle'第一级杆节点最大、最小载荷，N,最大应力范围比定义动态数组

DimPmaxi2（）AsSingle,Pmini2（）AsSingle,PL