油气储运工程外文翻译应用特征方法的天然气管道破裂问题的计算课件

1 中文译文中文译文 应用特征方法的天然气管道破裂问题的计算 J.A. Olorunmaiye ，机械工程系，Ilorin 大学，P.M.B. 1525, Ilorin(Nigeria) N.E. Imide，机械工程系，Abubakar Tafawa Balewa 大学，Bauchi (Nigeria) （完稿于 1992 年 7 月 4 日，校正后定稿于 1992 年 12 月 2 日） 摘 要 将在长距离、高压、伴随有突然断裂的天然气管道内的流动模拟为不稳定的一维 等温流动。集合的双曲线偏微分方程已由数值特征方法所解答。应用具备二次插值的 线性特性时，其数值计算的准确性是满足要求的。我们也发现，等温流动中特征曲率 的作用并不像 Flatt 的绝热流的报道中所说的那样明显。为了评估天然气管道破裂的危 险，有必要了解断点处气体的外泄速率。预测的断裂口处泄漏气体的质量流量要比应 用绝热流理论预算的小 18%，然而，这与早些时候利用等温流动理论的工作者得出的 结果却很好的吻合，其理论计算方法是基于加权残差的。 1. 说 明 天然气管道的瞬态压力波动通常是由阀门或者压缩机的突然开启与关闭引起的。 管道系统内存在压力波的传输与反射，从而会因为压力过大导致管道破裂。天然气工 业对长输管道在事故破裂后的气体不稳定流动是相当感兴趣的，因为管道内含有大量 的可燃气体，它的破裂也就意味着相当大的危害。为计算气体扩散范围和空气天然 气混合气的地面浓度，有必要了解气体的泄漏速率。破裂后，由于在低压段的流动条 件下，膨胀波进入管道后会造成流动逆转，低压段管道内的气体流动将比高压段更复 杂。在图 1 中，假定了破裂相当于管道截面大小，并立即向左方延伸（其高压管段） ， 阀门在断裂的瞬间完全关闭。 Fannelsp and Ryhming 1研究了高压下发生破裂时天然气管道内的不稳定流动。他 们将流动划分为三个时间状态。紧接着破裂发生的“早期”状态受波动进程影响，并且 开口端的压力接近周围环境压力。“中期”状态伴随着一个内部压力峰值的出现，大约 2 对应于流动逆转速度为 0 的位置。当压力峰顶到达低压尽端（截断处）时，“后期”状 态便出现了。在这一时段内，压力从关闭端到开口端单调递减。 图 1 破裂后膨胀波在管道内传播的物理环境 在假定等温流动的条件下，他们用积分法分析了管道内的不稳定流动。并且忽略 惯性条件得到了一组抛物线方程。不过，他们的方法无法解决“早期”状态的问题同时 在内部流动问题上也毫无建树，因为压力和压力流速积在 x=0 至 x=lo 上的空间分布 是被假定的，其目的是获得解决方法。 Lang and Fannelsp 2 使用了一些加权残差族中近似的程序在前人的成果上有所提 高。他们的流动没有以 Fannelsp and Ryhming 的工作为先验，但是却成为了解决方法 的一部分。他们能够获得更准确的结果，但是该方法在“早期”状态的预测中会产生因 变量的震荡。 Ryhming 3 利用了匹配渐进展开的方法来研究长输管道突然破裂后气体的不稳定 等温流动。他证明了流动处于临界状态下，管壁的摩擦使得流速顺着管道在破裂处变 的异常。他的复合第一阶解给出了不同时段的流速分布图。因为他的分析集中于预算 流速分布，并未考虑压力或者密度的变化，所以此方法不能预算“早期”状态的外泄速 率。 Flatt4研发了一种特征方法来估算长输管道破裂后天然气的不稳定绝热流动。因 为使用了全套的双曲线偏微分阶级方程，他的方法能够处理三个时态的流动问题。 然而，像管道破裂产生的高压变化，流体与管道周围媒介之间的热传递是相当可 观的，除非管道加有防冻保温层。事实上，流动过程会更接近等温而不是绝热状态， 3 因为密度和临界流速的补偿误差发生在前者。目前的努力都是针对于开发一种数值计 算方法来解决双曲线偏微分阶级方程以描述破裂发生时的不稳定等温流动。因此，这 就要向 Ryhming, Fannelep 和 Lang 提供了“早期”状态的分析中所缺少的数据表示祝贺。 它对三个时态的流动分析也同样有利。 Flatt4并没有将绝热模型应用到预测直到“后期”状态时的流动，因为这需要很长 的计算时间。由于等温流动方程（2 个）要比绝热方程（3 个）简单，这样用等温方程 解决问题会花费较少的时间，从而使其在计算所有时态内的流动时更经济。 2. 数学模型 对于等径管内的一维流动，其质量和动量守恒方程为： 式中： 将 u2表示为 u | u|以保证摩擦力方向总是与运动方向相反的。 理想气体等温流动的状态转移方程是： 联立式（1）和（2）,消掉密度 ，带入式（4）中，则连续性方程和动量方程变为： 和 其中，等温条件下的声速表示为 尽管等温条件下的声速是恒定的，膨胀波仍然会分散传播，如图 1。因为波头在流 体粒子中的传播具有正速度，然而波尾以声速 c 的传播相对于流体粒子来说却是负速 度。 以管道长度，等温声速，环境压力和温度为参考常量，方程（5）和（6）的无因 4 次形式为： 3. 数值方法 式（8）和（9）构成了一套拟线性的双曲线偏微分阶级方程，可以较好地利用特 征方法来解决。早些时候，Hartree 的方法也被 Olorunmaiye and Kentfield 5应用于这 项工作上。在该方法中，矩形网格加于集成领域并且方程沿着特征方向而整合。这里 的非独立变量有 P 和 U。方程（8）和（9）式的特征曲线由下式给出： 式（10）给出了压力波相对于流体以等温声速扩散时，其在 X-Z 坐标系中的传播方向。 特征协调方程有相对的斜率（U+1）和（U-1） ，分别为： 不同流速下到达格点的特征如图 2 所示。斜率为（U+1）和（U-1）的特征值分别 被标为 OV 和 TV。式（11）和（12）的有限分差近似值为： 重排式（13） 、 （14）得 变量的双脚注意味着在变量被两种脚注标示时采用其末点所表达的意义（见图 2） 。 5 图 2 连接到内部网格点的特征 4. 边界条件 流速|U|1 时，在右端只有特征值 OV 到达格点而且在左端只有 TV 到达格点。 4.1 右端 假设引起不稳定流的变化出现在左端，如图 1，在波动到达 Z* 之前，右端的压力 和流速是不受左端影响的。因此， 其中，UR和 PR是初始条件下右端的流量变量。 模拟在最糟糕的情况下，当压力波到达这个位置时右端的阀门被立即关闭。 该方程和关于 OV 的兼容性方程被用来获得右端的 PV值。 4.2 左端 在左端，要考虑两种截然不同的时间序列。一旦压力足够大于外部压力，流动便 会堵塞，而且 因为等温是一个粗略的近似条件，所以应该在堵塞发生的时候处理泄漏问题。为 了避免急速扩张的天然气温度下降，从周围环境至天然气有非常高的热传导速率将会 是很必要的。然而，由于持续时间短暂，泄漏堵塞时所损失的气量和三个时态内总的 天然气损失量相比是很小的。因此，假设的误差不影响最终的计算。 当压力降低到外部压力水平时，其状态为： 4.3 其他边界条件 该模型可能用到的其他边界条件有： 6 ()将静态压力描述为时间的函数。 ()将速度描述为时间的函数，这可以被用来模拟阀门的关闭。 其他的因变量可以通过利用恰当的通用方程得到。 7 5. 初始条件 先于不稳定流动开始之前，初始条件是管道内的稳定状态。长输管道中，与 4 fx/d 相比，忽略 ln( PIN/P)2（见参考6中的表（6.42）6），最初的恒定等温流的压力和速 度近似由下式给出： 6. 稳定性和准确性标准 时间步长的选择需要与 Courant-Friedrichs-Lewy 稳定性标准吻合，它要求： 最大允许时间步长的 80%用于该项计算的。 Flatt 4所定义的基于质量守恒定律的精度标准被用于这项工作。该精度标准被定 义为： 理想的值是 1。因此，误差为： Flatt 4介绍了一个系数 Kf，从而使流动处于声速时，对于边界格点的通用方程的 摩擦条件有所增加。Kf = 1 使关系保持恒定，然而 Kf = 0 则表示粘度条件被完全忽略。 他发现，在 Ryhming 的论文中所讨论的临界流动下破裂端的异常状态以 Kf= 1 计 算时引起的误差要比 Kf=0 计算时大。 7. 数值方案 针对特性曲线和配伍关系的方程已经由标准程序导出。然而，实际的特征方法的 数值执行方案不是标准的。在“后期”状态时，线性和二次插值可以用于获得在 OV 或 TV 特征底部的因变量的值。特征值也可以假定为线性的，或者曲率的影响可以像 Flatt 7所建议的那样。 8 采用的三个数值方案有： （）线性特性和线性插值方案； （）线性特性和二次插值方案； （）弯曲特性和二次插值方案。 更多关于计算的详细描述可以从其他处得知8。 8. 结果和讨论 为以上三种方案所写的 FORTRAN 语言程序运行于 IBM 兼容个人电脑的广域网程 序上，该电脑配置有使用 87 编译器的数学协处理器。 应用该程序的物理条件是长为 145km 内径为 0.87m 的水下管道，进口压力 133 个 大气压，出口压力 55 个大气压，气体温度 To=281K，流量约 650kg/s。外部压力 Po=6atm（相当于 50 米水深）。管道内充满了比热 y=1.3，气体常数 R=428.2 J/kg.K 的 天然气。这是在文献 1、2、4 中的案例。管道中的流动最初是稳态等温的，直到时间 t=0 即管道在高压端被突然破坏时。破坏发生前，恒定的摩擦系数 f=0.0018 造成了 78 巴的压降，此摩擦系数被用于非稳定流的计算。 8.1 三种方案的比较 表 1 显示了当非稳定流动状态持续 100s 时不同数量格点的误差和计算时间。Kf的 巨大影响显而易见。当 Kf =1 时，得到的误差更大而且其计算时间比 Kf=0 时更长。因 此当流动处于临界状态时，在随后的计算中像 Flatt 建议的那样，对于左端格点边界条 件的通用方程中忽略了粘度条件（Kf =0）。 表 1 在 t=0 到 t=100s 的时间区间内，对不同数量网点的三种数值方案的计算时间 和误差的比较 误差，t=100s 时的 e 值计算时间（s）网点 数量 Kf 值 LCLILCQICCQILCLILCQICCQI 00.3570.4290.43269.9275.7483.04 201 10.9671.0330.8889.4292.71119.91 00.2280.2310.233251.4280.56309.73 401 10.8370.6710.669282.6310.39342.52 00.1460.1240.124966.91087.741202.59 801 10.4370.3340.331148.081278.041296.9 00.0880.0660.0663927.394289.464745.56 1601 10.2310.170.1694431.554564.585223.91 9 续表 1 LCLI=线性特征和线性插值方案 LCQI=线性特征和二次插值方案 CCQI=弯曲特性和二次插值方案 如期望的那样，格点的数量越大，各种种方案的误差越小。在更高的格点数量时， 两种利用二次插值的方案产生的误差比另外一种利用线性插值的方案要少。 列表 1 中的误差值可能不会对应用于此项工作的数值方法的准确性有一个恰当的 评估。整体性的误差标准可能会比式（24）中提供的局部误差标准更好。即使式 （24）的数值评估本身是对准确性敏感的。 图 3 显示了有 401 个格点的三种数值方案预算的破裂端的质量流量的比较。因为 应用二次插值的两种方案预算的结果非常接近，而且计算时间是比较短的，所以最好 是采用二次插值方案的线性特征。 在此方案中的特征曲率的影响并不如 Flatt 的工作中那么重要。因为等温流动中的 声速 c 保持不变，那么特征 OV 或 TV（Uc）的斜率变化只是由速度量的小变化引起 的，特别是对于一个细网格。这不同于绝热流动的情况相对的特征斜率为（Ua） 而且流速和等熵声速（a）在 T 和 V、O 和 V 之间变化，从而使斜率变化比较细微。 其次，等温流动有两个个特征族，而绝热流动则有三个特征族。 图 3 使用了 401 个格点的三种方案所预算破裂端质量流量的比较。Flatt 应用绝热理论 得出的结果也一并给出。 10 8.2 试验案例 考虑巨大的长度/直径比，这种类型的问题不能执行于实验测试上，所以没有实验 结果可用来测试该模型。故而寻找其他方法来评估该模型是非常必要的。 图 4 伴随管道突然破裂产生的无摩擦不稳定等温流动模型预算的 Z=0.25 时管内流速分 布的比较，结果由封闭形式方法给出。 Ryhming3给出了“早期”状态得到的无摩擦不稳定等温流动的封闭形式解决方法， 处于最初管道内没有流动的时候： 考虑其几何结构。 该解决方案适用于 XZ，速度 U=0，因为促使左侧流向破裂端的膨胀 波还没有到达该点。图 4 和 5 显示了管道最初在 22.17 的无量纲压力下充满天然气时， 封闭形式方案的模型得出结果的比较。用这种方法得出结果的一致性是非常好的。 在第二个测试案例中，管道压力被加载到 22.17（参照与管道埋深处 6 个大气压的 环境压力），并且利用隔膜将左端从压力为 9.05 的低压储液槽中分离。左侧边界的隔 膜在 t=0 时被突然击穿。当瞬态结束后，流动应该与式（21）和（22）中的稳定等温 流动相一致。该模型预算的不同时期管内压力和流速的分布见图 6 和图 7.考虑到计算 中利用了粗网格的计算，可知在 t=10000s 时压力和流速的分布与稳定状态的解决方案 吻合的相当好。 11 图 5 伴随管道突然破裂产生的无摩擦不稳定等温流动模型所预算的在时间 X=0.25 时管 内流速变化的比较，结果由封闭形式方法给出。 图 6 连接两个储液槽中不同时期的流速分布（计算中用到了 101 个格点） 12 图 7 连接两个储液槽中不同时期的压力分布（计算中用到了 101 个格点） 图 8 管道破裂后的压力分布，用到了 401 个格点计算参数: t=O, 100, 200, 300, 400, 500, 700, 900, 1000, 2000 s. 13 8.3 关于北海道天然气管道的更多研究成果 破裂后的不同时期管线上的压力分布如图 8 所示。有趣的是，破裂发生后很长一 段时间（至少 1000s），离裂端的无量纲长度为 0.75 的地方，其压力并不改变。在 Flatt 的研究成果里，这样的一个点可以在他的绝热流动分析中距离破裂端 0.68 处找到， 条件是阀门处于低压端并且在 t=0 时关闭。如果管道压力在此点被监测，那么很长一 段时间里在这种标明的压力下将不会得到可预见的改变，即使在压力波到达该点之后 也是如此。 左右端压力随时间的变化可见图 9 和 10。一旦压力波到达，右端的压力便开始上 升，阀门将随之关闭，并且在压力下降之前阀门压力会达到峰值 13.7。 图 3 所表示的是 t=100s 和 t=300s 之间破裂端外泄的质量流量，由 Flatt 应用不稳定 绝热理论得出。目前利用等温流动理论得出的结果要比 Flatt 的结果低 18%。部分原因 是等熵声速高于等温声速这个事实。 图 11 显示了用不同数量的格点所得出的无量纲质量流量。共花费了 7 小时 57 分 钟 17 秒来计算实际流动时间为 2000s 的流量，用到了 401 个格点。使用 201 和 101 个 格点花费了 5 小时 54 分和 3 小时 45 分 61 秒，分别来计算实际流动时间为 7000s 和 20000s 的流量。这些结果与 Lang 和 Fannelsp 利用频谱配置方法和勒让德多项式得出 的结果很好的吻合，可见于图 11。 从图 3 中可知，最初的泄漏质量流量为 7,550 kg/s。到 100s 时流量减至 1630 kg/s。5 小时 30 分钟后，流量就小于 15 kg/s 了，见图 11. 14 图 9 破裂点随时间变化的预测压力变化 图 10 管道端部右侧预测压力 图 11 破裂处相对流量与时间的关系 15 9. 结论 本文提出了一种基于不稳定等温流动理论的数学模型，并且已通过特征方法所解 决。该模型所预测的结果与其他工作者的预测相一致。 虽然所预测的管道破裂端的流量比利用绝热流动得到的值低 18%，但它与 Lang 和 Fannelep 利用加权残差来求解不稳定等温流动方程所得出的结果很好的吻合。 距离管线破裂端 0.75 处，压力在很长时间内没有改变。 等温流动中的特征曲率并不像绝热流动中那么明显。因此，在不稳定等温流动的 计算中，没必要考虑特征曲率的影响。 该模型适用于与管道运行有关的其他不稳定流动中，例如先于关闭和维修而控制 外泄，和管道任一端突然的压力变化。这些操作中的压力波动不可能比破裂引起的波 动大。因此，该模型能够很容易处理这些情况。唯一有必要的是采用合适的边界条件。 鸣谢 感谢 Ilorin 大学计算机中心和工学部在项目初期和后期为该项目提供计算机设备。 第一作者非常感激在 Ilorin 大学内进行该项目研究时所获得的编号 No. 8: 15 : 303 参议 员研究经费。 16 符号 这里给出了适用的无因次形式。 罗马字符 c 等温声速.参考速度 d D=d/lo 管道直径 e 误差 f 摩擦系数 F F = Flo /c2 摩擦力每单位质量 I, M, N ,V 网格点 lo 管道长度.参考长度 O.T 特征根 P P= P/Po 管道横截面的平均压力 Po 环境压力.参考压力 R 气体常数 t Z=tc/lo 时间 To 环境温度.参考温度 u U= u/c 流速 x X=x/ lo 轴向距离 Z* 波从左端传播到右端需要的无因次时间 希腊字符 t Z 时间步长大小 x X 空间网格大小 精度标准 D=/o 密度 o 参考压力和温度下的密度 上标 IN 管道左端的初始条件 O 特征 OV 的底点 OV 特征 OV R 右端的因变量 T 特征 TV 的底点 TV 特征 TV V 网格点 V 下标 关于适当参考值的归一化 运算符 D/DZ 实质衍生 +/Z 相互斜率为（U+1）的特征区分 -/Z 相互斜率为（U+1）的特征区分 17 参考文献 1 T.K. Fannelrap and I.L. Ryhming, Massive release of gas from long pipelines, J. Energy,6 (1982) 132-140. 2 E. Lang and T.K. Fannelap, Efficient computation of the pipeline break problem. 3rd Symp. on Fuid-Transients in Fluid Structure Int