天然氣工程課程設(shè)計報告

1、天然氣工程課程設(shè)計報告姓 名： 王 巖 學 號： 2008* 班 級： 054*-* 專 業(yè)： 勘查技術(shù)與工程（油氣井）學 院： 工程學院 指導老師： 寧伏龍老師 目錄一、 程序運行圖示二、 程序代碼分析三、 該井節(jié)點分析四、 未來天然氣工業(yè)的發(fā)展趨勢和我國面臨的一些瓶頸問題一、程序運行圖示1.歡迎界面以及基本參數(shù)的輸入2.酸性氣體偏差系數(shù)計算(hy模型/dpk模型/dpr模型與wa校正法)3.天然氣相平衡計算(露點模型/泡點模型/等溫閃蒸模型)4天然氣兩相管流壓力計算(mukherjeebrill模型)5.天然氣水合物生成預測及抑制劑用量計算二、程序代碼分析1.氣體偏差系數(shù)計算（hy模型，d

2、pk模型，dpr模型）1.1hy模型+wa校正法1.1.1hy偏差系數(shù)計算模型 該法以 starling-carnahan 狀態(tài)方程為基礎(chǔ) ,通過對 standing-katz 圖版進行擬合 ,得到以下關(guān)系式:z = 0. 06125* ppr / (*tpr ) *exp - 1. 2 *(1 - 1/ tpr )2r 為擬對比密度 ,可用牛頓迭代法從如下方程求得:0 = (r +r2 +r3 r4)/(1 - r )3 -(14. 76 / tpr -9. 76 / tpr2+ 4. 58/tpr3)*r2+(90. 7/tpr-242. 2/tpr2+ 42. 4/tpr3)* r (2

3、. 18+2. 82/ tpr)-0. 01652( ppr / tpr )* exp - 1.2*(1 - 1/ tpr )2 該法應用范圍是:1. 2 tpr 3 ;0. 1 ppr 24. 0。部分代碼如下：private sub command1_click()dim n(11) as double -定義各組分dim i as integerfor i = 0 to 11n(i) = val(歡迎界面.摩爾分數(shù)(i) / 100next idim p as double, t as double-定義壓力、溫度p = val(絕對工作壓力)t = val(絕對工作溫度)dim tc(

4、11) as double, pc(11) as double-定義臨界溫度、臨界壓力tc(0) = 373.5833 pc(11) = 295.5315dim tpc as double, ppc as double-計算擬臨界壓力、擬臨界溫度tpc = 0ppc = 0for i = 0 to 11 tpc = tpc + n(i) * tc(i) ppc = ppc + n(i) * pc(i)next idim tpr as double, ppr as double-計算擬對比壓力、擬對比溫度tpr = t / tpcppr = p / ppcdim r as double-計算對比

5、密度(牛頓迭代法)dim x as double, x1 as double, f as double, f1 as doubleabsolution = 1x = 0do while absolution 0.00001f = (-0.06152) * (ppr / tpr) * exp(-1.2) * (1 - (1 / tpr) 2) + (x + x 2 + x 3 - x 4) / (1 - x) 3) - (14.76 / tpr) - (9.76 / (tpr 2) + (4.58 / (tpr 3) * (x 2) + (90.7 / tpr) - (242.2 / (tpr

6、2) + (42.4 / (tpr 3) * (x (2.18 + (2.82 / tpr)f1 = (1 + 4 * x + 4 * (x 2) - 4 * (x 3) + x 4) / (1 - x) 4) - (29.52 / tpr) - (19.52 / (tpr 2) + (9.16 / (tpr 3) * x + (2.18 + (2.82 / tpr) * (90.7 / tpr) - (242.2 / (tpr 2) + (42.4 / (tpr 3) * (x (1.18 + (2.82 / tpr)x1 = x - f / f1absolution = abs(x1 -

7、x)x = x1loopr = xdim z as double-用hy模型計算偏差系數(shù)z = (1 + r + r 2 - r 3) / (1 - r) 3) - (14.76 / tpr - 9.76 / (tpr 2) + 4.58 / (tpr 3) * r + (90.7 / tpr - 242.2 / (tpr 2) + 42.4 / (tpr 3) * (r (1.18 + (2.82 / tpr)hy模型z = format(z, 0.0000)-hy模型求得未校正的z1.1.2wa校正法wa校正法引入?yún)?shù),主要考慮了一些常見的極性分子( h2 s、 co2 )的影響 ,希望用

8、此參數(shù)來彌補常用計算方法的缺陷。參數(shù)的關(guān)系式如下:= 15 (m - m2) + 4. 167 ( n0. 5- n2)式中:m 為氣體混合物中 h2 s與 co2 的摩爾分數(shù)之和; n 為氣體混合物中 h2 s的摩爾分數(shù)。 根據(jù) wichert2aziz 的觀點 ,每個組分的臨界溫度和臨界壓力都應與參數(shù)有關(guān) ,臨界參數(shù)的校正關(guān)系式如下所示:t ci= tci - p ci= pci * t ci/ tci 式中: tci為i 組分的臨界溫度 ,k; pci為i 組分的臨界壓力 ,kpa ; tci為 i 組分的校正臨界溫度 , k; pci為 i組分的校正臨界壓力 ,kpa。部分代碼如下：d

9、im nhc as double, nh as double, as double, tc1(11) as double, pc1(11) as double-wa校正nhc = n(0) + n(2)nh = n(0) = 15 * (nhc - nhc 2) + 4.167 * (nh 0.5 - nh 2)for i = 0 to 11 tc1(i) = tc(i) - pc1(i) = pc(i) * tc1(i) / tc(i)next idim tpc1 as double, ppc1 as double, z1 as doubletpc1 = 0ppc1 = 0for i = 0

10、 to 11 tpc1 = tpc1 + n(i) * tc1(i) ppc1 = ppc1 + n(i) * pc1(i)next idim tpr1 as double, ppr1 as doubleif t = 17.24 then dim t1 as double t1 = t + 1.94 * (p / 2760 - (2.1 * (10 (-8) * (p 2) tpr1 = t1 / tpc1 ppr1 = p / ppc1 z1 = (1 + r + r 2 - r 3) / (1 - r) 3) - (14.76 / tpr1 - 9.76 / (tpr1 2) + 4.58

11、 / (tpr1 3) * r + (90.7 / tpr1 - 242.2 / (tpr1 2) + 42.4 / (tpr1 3) * (r (1.18 + (2.82 / tpr1) hy模型waz = format(z1, 0.0000)else tpr1 = t / tpc1 ppr1 = p / ppc1 z1 = (1 + r + r 2 - r 3) / (1 - r) 3) - (14.76 / tpr1 - 9.76 / (tpr1 2) + 4.58 / (tpr1 3) * r + (90.7 / tpr1 - 242.2 / (tpr1 2) + 42.4 / (tp

12、r1 3) * (r (1.18 + (2.82 / tpr1) hy模型waz = format(z1, 0.0000)-校正后的偏差系數(shù)zend ifend sub1.2dak模型+wa校正法1.2.1dak偏差系數(shù)計算模型該模型與 dranchuk-purvis-robinsion 計算法相同 ,但相對密度應采用牛頓迭代法從下式求出:0=1 + (a1 + a2/tpr+ a3/tpr3+ a4/tpr4+a5/tpr5)* r +(a6+a7/tpr+ a8/tpr2)*r2- a9*(a7/tpr+ a8/tpr3) r5 +a10/tpr3*r 2*(1 + a11*r2 )* e

13、xp ( - a11*r2 ) - 0. 27*ppr/(r *tpr) 式中: ai為給定系數(shù)。 該法應用范圍是:1. 0 tpr 3 ,0. 2 ppr 30 ;或者0. 7 tpr 1. 0 , ppr 0.00001f = (-0.27) * ppr / tpr + x + (a1 + a2 / tpr + a3 / (tpr 3) + a4 / (tpr 4) + a5 / (tpr 5) * (x 2) + (a6 + a7 / tpr + a8 / (tpr 2) * (x 3) - a9 * (a7 / tpr + a8 / (tpr 2) * (x 6) + a10 * (1

14、 + a11 * (x 2) * (x 3) / (tpr 3) * exp(-a11 * (x 2)f11 = 1 + 2 * (a1 + a2 / tpr + a3 / (tpr 3) + a4 / (tpr 4) + a5 / (tpr 5) * x + 3 * (a6 + a7 / tpr + a8 / (tpr 2) * (x 2) - 6 * a9 * (a7 / tpr + a8 / (tpr 2) * (x 5)f12 = (a10 / (tpr 3) * (3 * (x 2) + a11 * (3 * (x 4) - 2 * a11 * (x 6) * exp(-a11 *

15、(x 2)f1 = f11 + f12x1 = x - f / f1absolution = abs(x1 - x)x = x1loopr = xdim z as double-用dak模型計算偏差系數(shù)z = 1 + (a1 + a2 / tpr + a3 / (tpr 3) + a4 / (tpr 4) + a5 / (tpr 5) * r + (a6 + a7 / tpr + a8 / (tpr 2) * (r 2) - a9 * (a7 / tpr + a8 / (tpr 2) * (r 5) + a10 * (1 + a11 * (x 2) * (x 2) / (tpr 3) * ex

16、p(-a11) * (x 2)dak模型z = format(z, 0.0000)-dpk模型求得未校正的z1.2.2wa校正法同上，略。1.3dpr模型+wa校正法1.3.1dpr偏差系數(shù)計算模型dranchuk、purvis 和 robinsion 根據(jù) benedict-webb-rubin 狀態(tài)方程,將偏差系數(shù)轉(zhuǎn)換為對比壓力和對比溫度的函數(shù) ,于 1974 年推導出了帶 8 個常數(shù)的經(jīng)驗公式 ,其形式為:z = 1 +( a1 + a2/ tpr+ a3/tpr3)* r+ (a4 + a5/tpr)* r2+(a5 * a6/tpr)* r 5 - a7/tpr3* r2*(1 +

17、a8* r2 )* exp ( - a8* r2)式中: ai 為給定系數(shù); ppr為擬對比壓力 ,無因次; tpr為擬對比溫度 ,無因次。 dpr法使用 newton-raphson迭代法解非線性問題可得到偏差系數(shù)的值。這種方法的使用范圍是:1. 05 tpr 3 ;0. 2 ppr 30。部分代碼如下：-計算對比密度dim a1 as double, a2 as double, a3 as double, a4 as double, a5 as double, a6 as double, a7 as double, a8 as doublea1 = 0.31506237a2 = -1.04

18、67099a3 = -0.57832729a4 = 0.53530771a5 = -0.61232032a6 = -0.10488813a7 = 0.68157001a8 = 0.68446549dim x, x1 as double, f as double, f1 as doubleabsolution = 1x = 0do while absolution 0.00001f = -0.27 * ppr / tpr + x + (a1 + a2 / tpr + a3 / (tpr 3) * (x 2) + (a4 + a5 / tpr) * (x 3) + a5 * a6 / tpr *

19、(x 6) + a7 / (tpr 3) * (1 + a8 * (x 2) * (x 3) * exp(-a8 * (x 2)f1 = 1 + 2 * (a1 + a2 / tpr + a3 / (tpr 3) * x + 3 * (a4 + a5 / tpr) * (x 2) + 6 * a5 * a6 / tpr * (x 5) + a7 / (tpr 3) * (3 * (x 2) + a8 * (3 * (x 4) - 2 * a8 * (x 6) * exp(-a8 * (x 2)x1 = x - f / f1absolution = abs(x1 - x)x = x1loopr

20、= xdim z as double-用dpr模型計算偏差系數(shù)z = 1 + (a1 + a2 / tpr + a3 / (tpr 3) * x + (a4 + a5 / tpr) * (x 2) + a5 * a6 / tpr * (x 5) + a7 / (tpr 3) * (1 + a8 * (x 2) * (x 2) * exp(-a8) * (x 2)dpr模型z = format(z, 0.0000)-dpr型求得未校正的z1.3.2wa校正法同上，略。2.天然氣相平衡計算（露點，泡點，等溫閃蒸模型）2.1等溫閃蒸模型流體相平衡模型主要由三部分構(gòu)成：描述平衡氣液相組成、物質(zhì)的來那個

21、（摩爾數(shù)）及平衡常數(shù)與溫度、壓力關(guān)系的物料平衡方程組；描述平衡氣液相組成、物質(zhì)的量、平衡常數(shù)與逸度關(guān)系的熱力學平衡條件方程組及用于相平衡計算的狀態(tài)方程組機構(gòu)體系。三大方程組聯(lián)立即可分析油氣烴類體系相平衡模型。等溫閃蒸模型即油氣烴類體系的相態(tài)變化處于部分汽化和部分液化的狀態(tài)?？衫脀ilson平衡常數(shù)公式以及pr狀態(tài)方程求解分析。部分代碼如下：private sub command3_click()dim n(11) as double-獲得各組分摩爾分數(shù)dim i as integerfor i = 0 to 11 n(i) = val(歡迎界面.摩爾分數(shù)(i) / 100next idim

22、tc(11) as double, pc(11) as double-獲得臨界參數(shù)dim w(11) as double-偏心因子w(0) = 0.1dim kij(11, 11) as double-獲得二元交互作用參數(shù)kij(0, 0) = 0dim t as double, p as double-獲得一組t、p值t = val(溫度)p = val(壓力)dim k(11) as double, tr(11) as double, pr(11) as double-用wilson公式求解各組分平衡常數(shù)for i = 0 to 11 pr(i) = p / pc(i) tr(i) = t

23、 / tc(i)next ifor i = 0 to 11 k(i) = (exp(5.37 * (1 + w(i) * (1 - 1 / tr(i) / pr(i)next iresult = false-do while result = false-迭代相平衡常數(shù)kdim f as double, f1 as double, v as double, v1 as double-用牛頓迭代求vdim newtonv as doublev = 0.7newtonv = 1do while newtonv 0.001 for i = 0 to 11 f = f + n(i) * (k(i) -

24、 1) / (1 + (k(i) - 1) * v) f1 = f1 - n(i) * (k(i) - 1) 2) / (1 + (k(i) - 1) * v) 2) next iv1 = v - f / f1newtonv = abs(v1 - v) / v1)v = v1loopdim xx(11) as double, yy(11) as double-分別求各組分的液相摩爾數(shù)x和氣相摩爾數(shù)yfor i = 0 to 11 xx(i) = n(i) / (1 + (k(i) - 1) * v) yy(i) = n(i) * k(i) / (1 + (k(i) - 1) * v)next

25、idim x0 as double, x(11) as double-液相歸一化處理x0 = 0for i = 0 to 11 x0 = x0 + xx(i)next ifor i = 0 to 11 x(i) = xx(i) / x0next idim y0 as double, y(11) as double- 氣相歸一化處理y0 = 0for i = 0 to 11 y0 = y0 + yy(i)next ifor i = 0 to 11 y(i) = yy(i) / y0next idim r as double-計算pr狀態(tài)方程系數(shù)ai,bi,am,bm,am,bmr = 82.06

26、dim a(11) as double, b(11) as doublefor i = 0 to 11 a(i) = 0.45724 * (r 2) * (tc(i) 2) / pc(i) b(i) = 0.0778 * r * tc(i) / pc(i)next idim (11) as double, m(11) as doublefor i = 0 to 11 m(i) = 0.37464 + 1.54226 * w(i) - 0.26992 * (w(i) 2) (i) = (1 + m(i) * (1 - tr(i) 0.5) 2next idim aml as double, bm

27、l as double-液相混合物的平均引力系數(shù)am、平均斥力系數(shù)bmaml = 0bml = 0for i = 0 to 11 for j = 0 to 11 aml = aml + x(i) * x(j) * (a(i) * a(j) * (i) * (j) 0.5 * (1 - kij(i, j) next jnext ifor i = 0 to 11 bml = bml + x(i) * b(i)next idim amg as double, bmg as double-氣相混合物的平均引力系數(shù)am、平均斥力系數(shù)bm amg = 0bmg = 0for i = 0 to 11 for

28、 j = 0 to 11 amg = amg + y(i) * y(j) * (a(i) * a(j) * (i) * (j) 0.5 * (1 - kij(i, j) next jnext ifor i = 0 to 11 bmg = bmg + y(i) * b(i)next idim ag as double, bg as double-氣相zv-計算氣、液相的偏差系數(shù)zv、zlag = amg * p / (r 2 * t 2)bg = bmg * p / (r * t)dim zg0 as double, zg1 as doublezg0 = 0anther = 1do while

29、anther 0.0001 g = zg0 3 - (1 - bg) * zg0 2 + (ag - 2 * bg - 3 * bg 2) * zg0 - (ag * bg - bg 2 - bg 3) g1 = 3 * zg0 2 - 2 * (1 - bg) * zg0 + (ag - 2 * bg - 3 * bg 2) zg1 = zg0 - g / g1 anther = abs(zg1 - zg0) / zg1) zg0 = zg1loopdim zv as doublezv = zg0dim al as double, bl as double-液相zl -計算氣、液相的偏差系數(shù)

30、zv、zlal = aml * p / (r 2 * t 2)bl = bml * p / (r * t)dim zl0 as double, zl1 as doublezl0 = 0other = 1do while other 0.0001 h = zl0 3 - (1 - bl) * zl0 2 + (al - 2 * bl - 3 * bl 2) * zl0 - (al * bl - bl 2 - bl 3) h1 = 3 * zl0 2 - 2 * (1 - bl) * zl0 + (al - 2 * bl - 3 * bl 2) zl1 = zl0 - h / h1 other =

31、 abs(zl1 - zl0) / zl1) zl0 = zl1loopdim zl as doublezl = zl0dim (11) as double-氣相fig-氣、液相中各組分逸度fig、filfor i = 0 to 11 (i) = 0 for j = 0 to 11 (i) = (i) + y(j) * (a(i) * a(j) * (i) * (j) 0.5) * (1 - kij(i, j) next jnext idim fg(11) as doublefor i = 0 to 11 fg(i) = y(i) * p * exp(b(i) / bmg * (zv - 1)

32、 - log(zv - bg) - ag / (2 * 2 0.5 * bg) * (2 * (i) / amg - b(i) / bmg) * log(zv + 2.414 * bg) / (zv - 0.414 * bg) / (2 * bg * 2 0.5)next idim (11) as double-液相fil -氣、液相中各組分逸度fig、filfor i = 0 to 11 (i) = 0 for j = 0 to 11 (i) = (i) + x(j) * (a(i) * a(j) * (i) * (j) 0.5) * (1 - kij(i, j) next jnext id

33、im fl(11) as doublefor i = 0 to 11 fl(i) = x(i) * p * exp(b(i) / bml * (zl - 1) - log(zl - bl) - al / (2 * 2 0.5 * bl) * (2 * (i) / aml - b(i) / bml) * log(zl + 2.414 * bl) / (zl - 0.414 * bl) / (2 * bl * 2 0.5)next idim k1(11) as double-牛頓迭代調(diào)整平衡常數(shù)for i = 0 to 11 k1(i) = k(i) * fl(i) / fg(i)next i f

34、or i = 0 to 11 k(i) = k1(i)next i-判斷逸度系數(shù)是否滿足條件result = (abs(fl(0) - fg(0) = 0.001) and (abs(fl(1) - fg(1) = 0.0001) and (abs(fl(2) - fg(2) = 0.01) and (abs(fl(3) - fg(3) = 0.001) and (abs(fl(4) - fg(4) = 0.001) and (abs(fl(5) - fg(5) = 0.001) and (abs(fl(6) - fg(6) = 0.001) and (abs(fl(7) - fg(7) =

35、0.001) and (abs(fl(8) - fg(8) = 0.001) and (abs(fl(9) - fg(9) = 0.001) and (abs(fl(10) - fg(10) = 0.001) and (abs(fl(11) - fg(11) 0.001 for i = 0 to 11 f = f + n(i) * (k(i) - 1) / (1 + (k(i) - 1) * v) f1 = f1 - n(i) * (k(i) - 1) 2) / (1 + (k(i) - 1) * v) 2) next iv1 = v - f / f1newtonv = abs(f1 - f)

36、v = v1loopend sub2.3泡點計算 泡點壓力計算步驟同露點壓力類似，只需把x(i) = n(i) / k(i)；y(i) = n(i)換成y(i) = n(i) * k(i)；x(i) = n(i)即可，此處不再贅述。3.天然氣兩相管流壓力計算（mukherjeebrill模型）1985年mukherjee和brill在beggsbrill（1973）究工作的基礎(chǔ)上，改進實驗條件，對傾斜管兩相流進行了深入的研究，提出了適用于傾斜管及垂直管兩相流應用方便的持液率及摩阻系數(shù)經(jīng)驗公式。這是目前不僅適用于垂直管，同時也適用于水平傾斜管的兩相管流方面較全面的研究成果。m-b模型的壓力梯度方

37、程為： mukherjeebrill持液率公式共有三個：一個用于水平流；另外兩個用分別用于下降的分層流和其它流型。mukherjeebrill持液率形式簡單，只是控制流型的三個無因次變量的函數(shù)。 式中管斜角（與水平方向的夾角）；無因次液相粘度: 無因次液相速度: 無因次氣相速度: 用mukherjeebrill確定摩阻系數(shù)時，需判別其流態(tài)，對于定向井，流體是向上或水平流動的，流態(tài)主要為泡流段塞流和霧狀流，其流態(tài)判別式為： 若，則為泡流段塞流，兩相摩阻系數(shù)為無滑脫摩阻系數(shù)； 無滑脫雷諾數(shù):無滑脫混合物密度: 無滑脫混合物粘度:無滑脫持液率:若，則為霧狀流，兩相摩阻系數(shù)為相對持液率和無滑脫摩阻系數(shù)

38、的函數(shù)，即： 程序框圖如下：部分代碼如下：dim c1 as double, c2 as double, c3 as double, c4 as double, c5 as double, c6 as doubleprivate sub option1_click()-m-b持液率公式回歸系數(shù) c1 = -0.380113 c2 = 0.129875 c3 = -0.119788 c4 = 2.343227 c5 = 0.475686 c6 = 0.288657end subprivate sub option2_click() c1 = -1.330282 c2 = 4.808139 c3

39、= 4.171584 c4 = 56.262268 c5 = 0.079951 c6 = 0.504887end subprivate sub option3_click() c1 = -0.516644 c2 = 0.789805 c3 = 0.551627 c4 = 15.519214 c5 = 0.371771 c6 = 0.393952end subprivate sub command1_click()-定義以及輸入需要的基本參數(shù)const as double = 0.06const g as double = 9.81const pi as double = 3.141592653const e as double = 0.000016 管壁絕對粗糙度dim p0 as double, t0 as double, gt as double, d as double, qsc as double, ql as double, ug as