天然气工程课程设计报告.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 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-迭代相平衡常数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) - 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和气相摩尔数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 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,bmr = 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、平均斥力系数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-气相混合物的平均引力系数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、zlag = 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、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 = 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) - 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 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 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)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）究工作的基础上，改进实验条件，对倾斜管两相流进行了深入的研究，提出了适用于倾斜管及垂直管两相流应用方便的持液率及摩阻系数经验公式。这是目前不仅适用于垂直管，同时也适用于水平倾斜管的两相管流方面较全面的研究成果。m-b模型的压力梯度方程为： mukherjeebrill持液率公式共有三个：一个用于水平流；另外两个用分别用于下降的分层流和其它流型。mukherjeebrill持液率形式简单，只是控制流型的三个无因次变量的函数。 式中管斜角（与水平方向的夹角）；无因次液相粘度: 无因次液相速度: 无因次气相速度: 用mukherjeebrill确定摩阻系数时，需判别其流态，对于定向井，流体是向上或水平流动的，流态主要为泡流段塞流和雾状流，其流态判别式为： 若，则为泡流段塞流，两相摩阻系数为无滑脱摩阻系数； 无滑脱雷诺数:无滑脱混合物密度: 无滑脱混合物粘度:无滑脱持液率:若，则为雾状流，两相摩阻系数为相对持液率和无滑脱摩阻系数的函数，即： 程序框图如下：部分代码如下：dim c1 as double, c2 as double, c3 as double, c4 as double, c5 as double, c6 as doubleprivate sub option1_click()-m-b持液率公式回归系数 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 = 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()-定义以及输入需要的基本参数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