天然气工程课程设计报告.doc
天然气工程课程设计报告姓 名： 王 岩 学 号： 2008* 班 级： 054*-* 专 业： 勘查技术与工程（油气井）学 院： 工程学院 指导老师： 宁伏龙老师 目录一、 程序运行图示二、 程序代码分析三、 该井节点分析四、 未来天然气工业的发展趋势和我国面临的一些瓶颈问题一、程序运行图示1.欢迎界面以及基本参数的输入2.酸性气体偏差系数计算(HY模型/DPK模型/DPR模型与WA校正法)3.天然气相平衡计算(露点模型/泡点模型/等温闪蒸模型)4天然气两相管流压力计算(MukherjeeBrill模型)5.天然气水合物生成预测及抑制剂用量计算二、程序代码分析1.气体偏差系数计算（HY模型，DPK模型，DPR模型）1.1HY模型+WA校正法1.1.1HY偏差系数计算模型 该法以 Starling-Carnahan 状态方程为基础 ,通过对 Standing-Katz 图版进行拟合 ,得到以下关系式: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. 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(欢迎界面.摩尔分数(i) / 100Next iDim p As Double, T As Double-'定义压力、温度'p = Val(绝对工作压力)T = Val(绝对工作温度)Dim Tc(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-'计算对比密度(牛顿迭代法)'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 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 - x)x = x1Loopr = xDim z As Double-'用HY模型计算偏差系数'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校正法引入参数,主要考虑了一些常见的极性分子( H2 S、 CO2 )的影响 ,希望用此参数来弥补常用计算方法的缺陷。参数的关系式如下:= 15 (M - M2) + 4. 167 ( N0. 5- N2)式中:M 为气体混合物中 H2 S与 CO2 的摩尔分数之和; N 为气体混合物中 H2 S的摩尔分数。 根据 Wichert2Aziz 的观点 ,每个组分的临界温度和临界压力都应与参数有关 ,临界参数的校正关系式如下所示:T ci= Tci - p ci= pci * T ci/ Tci 式中: Tci为i 组分的临界温度 ,K; pci为i 组分的临界压力 ,kPa ; Tci为 i 组分的校正临界温度 , K; pci为 i组分的校正临界压力 ,kPa。部分代码如下：Dim 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 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 / (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 / (Tpr1 3) * (r (1.18 + (2.82 / Tpr1) HY模型WAZ = Format(z1, "0.0000")-校正后的偏差系数ZEnd IfEnd Sub1.2DAK模型+WA校正法1.2.1DAK偏差系数计算模型该模型与 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 )* exp ( - A11*r2 ) - 0. 27*ppr/(r *Tpr) 式中: Ai为给定系数。 该法应用范围是:1. 0 Tpr 3 ,0. 2 ppr 30 ;或者0. 7 Tpr 1. 0 , ppr < 1. 0。部分代码如下：-'计算对比密度（牛顿迭代法）'Dim A1 As Double, A2 As Double, A3 As Double, A4 As Double, A5 As Double, A6 As Double, A7 As Double, A8 As Double, A9 As Double, A10 As Double, A11 As DoubleA1 = 0.3265A2 = -1.07A3 = -0.5339A4 = 0.01569A5 = -0.05165A6 = 0.5475A7 = -0.7361A8 = 0.1844A9 = 0.1056A10 = 0.6134A11 = 0.721Dim r As DoubleDim f11 As DoubleDim f12 As DoubleDim x As Double, 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) + 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 + 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 * (x 2)F1 = f11 + f12x1 = x - F / F1absolution = Abs(x1 - x)x = x1Loopr = xDim z As Double-'用DAK模型计算偏差系数'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) * Exp(-A11) * (x 2)DAK模型Z = Format(z, "0.0000")-DPK模型求得未校正的Z1.2.2WA校正法同上，略。1.3DPR模型+WA校正法1.3.1DPR偏差系数计算模型Dranchuk、Purvis 和 Robinsion 根据 Benedict-Webb-Rubin 状态方程,将偏差系数转换为对比压力和对比温度的函数 ,于 1974 年推导出了带 8 个常数的经验公式 ,其形式为:Z = 1 +( A1 + A2/ Tpr+ A3/TPr3)* r+ (A4 + A5/Tpr)* r2+(A5 * A6/Tpr)* r 5 - A7/TPr3* r2*(1 + A8* r2 )* exp ( - A8* r2)式中: Ai 为给定系数; ppr为拟对比压力 ,无因次; Tpr为拟对比温度 ,无因次。 DPR法使用 Newton-Raphson迭代法解非线性问题可得到偏差系数的值。这种方法的使用范围是: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.0467099A3 = -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 * (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 = xDim z As Double-'用DPR模型计算偏差系数'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等温闪蒸模型流体相平衡模型主要由三部分构成：描述平衡气液相组成、物质的来那个（摩尔数）及平衡常数与温度、压力关系的物料平衡方程组；描述平衡气液相组成、物质的量、平衡常数与逸度关系的热力学平衡条件方程组及用于相平衡计算的状态方程组机构体系。三大方程组联立即可分析油气烃类体系相平衡模型。等温闪蒸模型即油气烃类体系的相态变化处于部分汽化和部分液化的状态。可利用Wilson平衡常数公式以及PR状态方程求解分析。部分代码如下：Private Sub Command3_Click()Dim n(11) As Double-'获得各组分摩尔分数'Dim i As IntegerFor i = 0 To 11 n(i) = Val(欢迎界面.摩尔分数(i) / 100Next iDim Tc(11) As Double, Pc(11) As Double-'获得临界参数'Dim W(11) As Double-'偏心因子'W(0) = 0.1Dim kij(11, 11) As Double-'获得二元交互作用参数'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公式求解各组分平衡常数'For i = 0 To 11 Pr(i) = p / Pc(i) Tr(i) = T / 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-'迭代相平衡常数K'Dim F As Double, F1 As Double, V As Double, V1 As Double-'用牛顿迭代求V'Dim newtonV As DoubleV = 0.7newtonV = 1Do While newtonV > 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(V1 - V) / V1)V = V1LoopDim xx(11) As Double, yy(11) As Double-'分别求各组分的液相摩尔数x和气相摩尔数y'For 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 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状态方程系数ai,bi,am,bm,Am,Bm'R = 82.06Dim 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, bml As Double-'液相混合物的平均引力系数am、平均斥力系数bm'aml = 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-'气相混合物的平均引力系数am、平均斥力系数bm 'amg = 0bmg = 0For i = 0 To 11 For 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'-'计算气、液相的偏差系数ZV、ZL'Ag = amg * p / (R 2 * T 2)Bg = bmg * p / (R * T)Dim Zg0 As Double, Zg1 As DoubleZg0 = 0anther = 1Do While 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 '-'计算气、液相的偏差系数ZV、ZL'Al = 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 = Abs(Zl1 - Zl0) / Zl1) Zl0 = Zl1LoopDim Zl As DoubleZl = Zl0Dim (11) As Double-'气相fig'-'气、液相中各组分逸度fig、fil'For 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) - 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、fil'For 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 iDim 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-牛顿迭代调整平衡常数For i = 0 To 11 K1(i) = K(i) * fl(i) / fg(i)Next i For i = 0 To 11 K(i) = K1(i)Next i-判断逸度系数是否满足条件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) <= 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) Loop-'输出'For i = 0 To 11 等温闪蒸液相i组分摩尔分数(i) = Format(x(i), "0.00000") 等温闪蒸气相i组分摩尔分数(i) = Format(y(i), "0.00000")Next i等温闪蒸气相总摩尔数 = Format(V, "0.00000")Dim L As Double-液相总摩尔数LL = 1 V等温闪蒸液相总摩尔数 = Format(L, "0.00000")等温闪蒸气相偏差系数 = Format(ZV, "0.00000")等温闪蒸液相偏差系数 = Format(Zl, "0.00000")End Sub2.2露点模型露点模型与等温闪蒸模型原理一样，
