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

1、作業(yè)一：計算井筒多相流的Beggs-Brill方法基本參數(shù)：套管尺寸：5 1/2, 油管：2 1/2, GOR=24m3/ m3, o=0.9659, o=100cp, 井底溫度63，井口溫度25，井口壓力1MPa，產(chǎn)液137方/天，油層深度1490m，井底流壓8MPa，含水率29，泵的單級揚程4.6m。設計步驟：截圖：Sub BeggsBrill(cal As Single, rw As Single, ro As Single, rg As Single, yco2 As Single, yh2s As Single, yn2 As Single, ynacl As Single, _pp

2、c As Single, tpc As Single, d2 As Single, d1 As Single, eras As Single, lon As Single, alfa0 As Single, T1 As Single, pcal1 As Single, T2 As Single, qw As Single, qo As Single, Qsc As Single, daltal As Single, PCAL2 As Single)Do If lsum + daltal > lon Then daltal = lon - lsum End If lmid = lsum +

3、 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 Loop Until lsum >= lon PCAL2 = p2 End Sub psc = 0.101325 tsc = 293 densc = 1.204 * rg deno = 1000 * ro denw = 1000 * r

4、w qow = qo + qw If qo > 0 Then gor = Qsc / qo Call PVTO(ro, rg, gor, T, p, pbo, rso, Bo, uo, sto) End If If qw > 0 Then gwr = Qsc / qw Call PVTW(T, p, ynacl, gwr, rsw, bw, uw, stw) End If If qow > 0 Then GLR = Qsc / qow RS = (rso * qo + rsw * qw) / qow If RS >= GLR Then RS = GLR End If q

5、gsc = qow * (GLR - RS) Else qgsc = Qsc End If If qgsc > 0 Then ppr = p / ppc tpr = T / tpc Call dprzf(ppr, tpr, dr, z, cpr) bg = psc / tsc * z * T / p Call lgeug(rg, yn2, yco2, yh2s, p, T, z, ug) qg = qgsc * bg wg = densc * qgsc densg = wg / qg Else qg = 0 densg = 0 bg = 0 End If Ql = qo * Bo + q

6、w * bw wl = deno * qo + denw * qw + densc * RS * qow If Ql > 0 Then densl = wl / Ql fo = qo * Bo / Ql End If ul = uo * fo + uw * (1 - fo) stal = sto * fo + stw * (1 - fo)End Sub gor = gorin * 5.615 T = tin * 1.8 - 460 p = 145 * pin api = 141.5 / ro - 131.5 If api <= 30 Then a1 = 0.0362 b1 = 1.

7、0937 c1 = 25.724 Else a1 = 0.0178 b1 = 1.187 c1 = 23.931 End If tp = T + 460 pbo = (gor / (a1 * rg * Exp(c1 * api / tp) (1 / b1) rso = a1 * rg * p b1 * Exp(c1 * api / tp) 'ft3/bbl If p >= pbo Or rso >= gor Then rso = gor End If If api <= 30 Then a1 = 1.761 * 10 (-5) b1 = 0.0004677 c1 =

8、-1.1811 * 10 (-8) Else a1 = 0.000011 b1 = 0.000467 c1 = 1.377 * 10 (-9) End If 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 If p >= pbo Or rso >= gor Then Bo = Bo * Exp(co * (pbo - p) zz =

9、 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 ' If pbo > 0 And (p >= pbo Or rso >= gor) Then uo = uo * (p / pbo) (2.6 * p 1.187 * Exp(-0.0000898 * p - 11.513) End

10、If sto = (42.4 - 0.047 * T - 0.267 * api) * Exp(-0.0007 * p) / 1000 pbo = pbo / 145 rso = rso / 5.615End Sub 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

11、) * T 2 scnacl = 1 - (0.0753 - 0.000173 * T) * ynacl rws = (a1 + b1 * p + c1 * p 2) * scnacl / 5.615 If rws > gwr Then rws = gwr End If If rws <= 0 Then 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) -

12、 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 End If tp = T - 60 tl = 5.1 * 10 (-8) * p + (5.47 * 10

13、 (-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) *

14、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.2046End Sub'

15、;程序功能：計算當?shù)貤l件下的有關參數(shù)的過程 ga = 9.807don = 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) ' If Vsl > 0 Then nv = Vsl * (densl / (ga * stal) 0.25 ' End If denn = na * de

16、nsl + (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 'End Sub If na > 0 Then ls1 = 316 * na 0.302 ls2 = 0.0009252 * na (-2.4684) ls3 = 0.1 * na (-1.4516) ls4 = 0.5 * na (-6.738)

17、 If na < 0.00001 Then ident = 5 ElseIf na > 0.99999 Then ' ident = 6 ElseIf na < 0.01 And nr < ls1 Or na >= 0.01 And nr < ls2 Then ' ident = 1 ElseIf (na >= 0.01 And na < 0.4) And (nr > ls3 And nr <= ls1) Or na >= 0.4 And (nr > ls3 And nr <= ls4) Then i

18、dent = 2 ElseIf na < 0.4 And nr >= ls1 Or na >= 0.4 And nr >= ls4 Then ident = 3 ElseIf na >= 0.01 And (nr >= ls2 And nr <= ls3) Then ident = 4 End If Else ident = 5 End IfEnd Sub If alfa >= 0 Then Select Case ident Case 1 dd = 0.011 ee = -3.768 ff = 3.539 gg = -1.614 Case 2

19、dd = 2.96 ee = 0.305 ff = -0.4473 gg = 0.0978 Case 3 dd = 1 ee = 0 ff = 0 gg = 0 End Select Else dd = 4.7 ee = -0.3692 ff = 0.1244 gg = -0.5056 End If cc = (1 - na) * Log(dd * na ee * nv ff * nr gg) ' If cc < 0 Then cc = 0 End If tm = Sin(1.8 * alfa) hui = 1 + cc * (tm - tm 3 / 3) ' If hu

20、i < 0 Then hui = 0 End If Select Case ident Case 1 dd = 0.98 ee = 0.4846 cc = 0.0868 Case 2 dd = 0.845 ee = 0.5351 cc = 0.0173 Case 3 dd = 1.065 ee = 0.5824 cc = 0.0609 End Select hl0 = dd * na ee * nr (-cc) ' If hl0 < na Then hl0 = na End If hlsta = 0.918 * hl0 * hui If hlsta < 0 Then

21、hlsta = 0 End If If hlsta > 1 Then hlsta = 1 End IfEnd Sub ga = 9.807 Select Case ident Case 4 ls2 = 0.0009252 * na (-2.4684) ls3 = 0.1 * na (-1.4516) ident = 1 Call HOLDUP(alfa, ident, na, nv, nr, hlsta) hlsta1 = hlsta ident = 2 Call HOLDUP(alfa, ident, na, nv, nr, hlsta) hlsta2 = hlsta tm = (ls

22、3 - nr) / (ls3 - ls2) hlsta = tm * hlsta1 + (1 - tm) * hlsta2 Case 5 hlsta = 0 Case 6 hlsta = 1 Case Else Call HOLDUP(alfa, ident, na, nv, nr, hlsta) End Select If hlsta < na Then hlsta = na End If denm = densl * hlsta + densg * (1 - hlsta) tm = denm * ga * Sin(alfa) If hlsta = 0 Or hlsta = 1 The

23、n sb = 0 Else yb = na / hlsta 2 yr = Log(yb) If yb > 1 And yb < 1.2 Then sb = Log(2.2 * yb - 1.2) Else sb = yr / (-0.0523 + 3.182 * yr - 0.8725 * yr 2 + 0.01853 * yr 4) '=16 End If End If fm = fmn * Exp(sb) Tf = fm * denn * Vm 2 / (2 * don) coef = 1 - denm * Vm * vsg / (1000000 * p) T1 = (

24、tm + Tf) / (1000000 * coef) daltap = daltal * T1 End Sub 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

25、 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 If dr = 0 Or dr = 1 Then Exit Sub End If Loop Until Abs(f / fp) < 0.00001 z

26、 = 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)End Sub mg = 28.97 * rg R = 0.008314 g1 = yh2s * (5.7 * rg - 1.7) * 0.00001 g2 = yco2 * (5 * rg + 1.

27、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)End Sub作業(yè)二 抽油桿柱設計基本參數(shù)

28、配產(chǎn)：30m3/d,泵深：2300，動液面：950m，沖程：3.6，沖次6次/min，泵徑：56mm，油管內(nèi)徑：73mm，地層水粘度：1,8mPa.s，地層水密度：1050kg/m3，抽油桿密度：7850kg/m3，最小桿徑：19mm，桿使用系數(shù)：0.95設計步驟：截圖：Dim hc As Single '沉沒度，mDim SF As Single '桿使用系數(shù)Dim s As Single '沖程，mDim ql As Single '產(chǎn)液量，m3/dDim dp As Single '泵徑，mmDim de As Single '柱塞與泵筒間

29、隙平均值，mDim uw As Single '地層水粘度，mPa.sDim n As Single '沖次，次/minDim rouw As Single '地層水密度，kg/m3Dim rous As Single '鋼的密度，kg/m3Dim t As Single '抽油桿最小抗張強度，MPaDim Pl As Single, PLi As Single, deltaPL As Single '輸入應力范圍比、計算應力范圍比、兩級抽油桿頂端面應力范圍比的最大允許差值Dim dt As Single '油管內(nèi)徑，mmDim fp A

30、s Single, fr As Single, fo As Single '柱塞截面積，抽油桿截面積，閥孔截面積，m2Dim nrod As Integer '桿級數(shù)Dim L0 As Single, L1 As Single '計算所得桿柱長度，m 已設計出的桿柱總長度，mDim dfo As Single '閥孔徑，mDim qr() As Single '每米抽油桿質(zhì)量，kg/mDim Li() As Single '每級桿長度，mDim dr() As Single '抽油桿直徑,mmDim b As Single '失重

31、系數(shù)Dim m As Single '油管內(nèi)徑與抽油桿直徑之比Dim vmax As Single '抽油桿柱最大下行速度，m/sDim vp As Single, vf As Single '柱塞最大運動速度、閥孔液流速度，m/sDim NRe As Single, ksi As Single '雷諾數(shù)，閥流量系數(shù)Dim Pmax0 As Single, Pmin0 As Single, PL0 As Single '抽油桿上下沖程底端載荷，N，底端應力范圍比Dim () As Single, dL() As Single, d() As Single

32、, d() As Single '井斜角，井斜深、井斜角、方位角變化量Dim Pmaxi1() As Single, Pmini1() As Single, PLi1() As Single '第一級桿節(jié)點最大、最小載荷，N, 最大應力范圍比 定義動態(tài)數(shù)組Dim Pmaxi2() As Single, Pmini2() As Single, PLi2() As Single '第二級桿節(jié)點最大、最小載荷，N，最大應力范圍比Dim Pmaxi3() As Single, Pmini3() As Single, PLi3() As Single '第三級桿節(jié)點最大、

33、最小載荷，N，最大應力范圍比Dim Fli As Single '變徑處產(chǎn)生的浮力，NDim max As Single, min As Single, all As Single '節(jié)點最大應力、最小應力、許用應力，MPaDim f As Single '摩擦系數(shù)Dim Fcp As Single '泵柱塞與泵筒間摩擦力，NDim Fvl As Single '液流通過一個游動閥產(chǎn)生的柱塞下行阻力，NDim i As Single, j As Single If Me.cmbMaterial.Text = "K" Then t =

34、570 Else If Me.cmbMaterial.Text = "C" Then t = 610 Else If Me.cmbMaterial.Text = "D" Then t = 794 Else If Me.cmbMaterial.Text = "H" Then t = 980 Else If Me.cmbMaterial.Text = "" Then Exit Sub End If Pl = Val(txtPL.Text) ' deltaPL = Val(txtDeltaPL.Text) 

35、9; hp = Val(txtHp.Text) ' hd = Val(txtHd.Text) ' hc = hp - hd ' s = Val(txtS.Text) n = Val(txtN.Text) vmax = pi * s * n / 3600 ' ql = Val(txtQl.Text) dp = Val(txtDp.Text) / 1000 de = 0.056 / 1000 dt = Val(txtDt.Text) / 1000 uw = Val(txtUw.Text) / 1000 rouw = Val(txtRouw.Text) rous =

36、Val(txtRows.Text) SF = Val(txtSF.Text) '桿使用系數(shù) b = (rous - rouw) / rous '失重系數(shù) f = 0.015 '摩擦系數(shù) fp = dp 2 * pi / 4 '柱塞面積 fr = (Val(txtDrMin.Text) / 1000) 2 * pi / 4 '油管截面積 Fcp = 0.94 * dp / de - 140 '計算泵柱塞與泵筒間摩擦力，N '*計算液流通過一個游動閥產(chǎn)生的柱塞下行阻力，N P110* Call selectDfo(dp, dfo) '

37、根據(jù)泵徑確定閥孔徑 fo = dfo 2 * pi / 4 '閥孔截面積，m vp = pi * s * n / 3600 ' vf = fp * vp / fo '閥孔液流速度 NRe = rouw * dfo * vf / uw '雷諾數(shù) Call NRe_u(NRe, ksi) '根據(jù)雷諾數(shù)計算閥孔流量系數(shù) P110 Fvl = rouw * fp * (fp * s * n / fo) 2 / (729 * ksi 2) '液流通過一個游動閥產(chǎn)生的柱塞下行阻力，N P110 For i = 0 To hp + 1 ReDim Preser

38、ve (i) ReDim Preserve dL(i) ReDim Preserve d(i) ReDim Preserve d(i) (i) = 0 dL(i) = 1 '每一節(jié)點長度設置為1 d(i) = 0 d(i) = 0 Next iPmax0 = Fcp + (fp - fr) * hp * rouw * g - (hc * rouw * g * fp - Fvl) Pmin0 = -Fcp - hp * rouw * g * fr max = Abs(Pmax0) / fr / 1000000 min = Abs(Pmin0) / fr / 1000000 all = (

39、t / 4 + 0.5625 * min) * SF PL0 = (max - min) / (all - min) '底端應力范圍比 If PL0 > Pl Then '底端應力范圍比 > 輸入應力范圍比 MsgBox "所選材料和最小桿徑的抽油桿強度不夠，請調(diào)整！", vbOKOnly + vbExclamation, "提示" Exit Sub End If L0 = 0 '計算所得桿柱長度，m L1 = 0 '已設計出的桿柱總長度，m Do nrod = 0 dr(nrod) = Val(txtDrMi

40、n.Text) / 1000 Call selectqr(dr(nrod) * 1000, qr(nrod) '根據(jù)桿徑選擇每米長桿重 Li(nrod) = hp '令泵深 = 第一級桿柱長度'計算L1，Pmax1，Pmin1，PL1，一級桿 （ 一級桿長，最大載荷，最小載荷，應力范圍比） Call RodDesign(dr(nrod), qr(nrod), dt, b, fp, s, n, uw, vmax, hp, f, t, SF, Pl, , dL, d, d, 0, 0, Pmax0, Pmin0, Pmaxi1, Pmini1, PLi1, L0) '

41、;每一個rod design出來一個L0 If L0 < hp Then '這里的L0是上一個sub算出來的 Li(nrod) = L0 '將L0賦給Li（0） nrod = nrod + 1 '* ReDim Preserve dr(nrod) As Single ReDim Preserve qr(nrod) As Single ReDim Preserve Li(nrod) As Single '* dr(nrod) = dr(nrod - 1) + 0.003 '下一級桿柱直徑 = 上一級 + 3mm Call selectqr(dr(nrod) * 1000, qr(nrod) '根據(jù)桿徑選擇每米長桿重 Li(nrod) = hp - Li(nrod - 1) '下一級桿柱長度 = 泵深 - 上一級 Fli = rouw * g * (dr(nro