天然气管道运行模拟及仿真技术研究1011202045 蔡永军 科学计算选讲结课论文为了预测天然气管道运行状态，制定合理的管输计划，更好的配置设备开机，天然气管道输送过程中需要进行工况模拟及仿真。

实际工作中需要建立压缩机、阀门等设备的模型，确定管段的控制方程、气体的状态方程，针对给出的初始条件和边界条件，筛选确定天然气管网数学模型的离散方法与非线性方程组的求解算法寻找合理的非线性方程的求解算法，得到合理的数值解。

1天然气管道仿真数学模型 1.1管段的控制方程对于管道中的任意管段，经过适当的简化可以用下列公式来描述： 连续性方程：()0A At xρρμ∂∂+=∂∂ （1） 运动方程：2()(.)sin()2A P A A A g A t x x Dρμρμμρμρθλ∂∂∂+=---∂∂∂ （2）能量方程：221(())(.())22sin()()W Ph A h PAAA g Dk T T txxμμρρμρρμθπ∂-+∂+∂+=----∂∂∂ （3）式中：A ——管道的横截面积，m 2；ρ——流体密度，kg/m 3； t ——时间，s ； x ——坐标，m ； u ——速度，m/s ； P ——压力，Pa ； θ——管道倾角，rad ； λ——水力摩阻系数； D ——管道内径，m ； T ——流体温度，k ；k 1——流体至管壁的换热系数；h ——比焓；T w ——管壁的温度，k 。

1.2 阀门控制方程阀门控制方程如下：120dw up g M M MC h h ρ-=-== （4） 式中： M up ——阀门入口质量流量，kg/s ；M dw ——阀门入口质量流量，kg/s C g ——阀门系数；P up ——阀的入口压力，Pa ；P dw ——阀的出口压力，Pa 。

1.3压缩机控制方程简化后的压缩机控制方程如下222111001()()dw up fuel m mdw up n n a b Q c Q n n M M M T T εε-=++-== （5）式中：ε——压缩机压比；m ——多变压缩指数；n ——压缩机的实际转速，rpm ； n 0——压缩机的额定转速，rpm ； a 1, b 1, c 1——系数；Q ——给定状态下的体积流量，m 3/s ； 1.4 理想调节阀阀控制方程理想调节阀控制方程如下：12dw up dw M M P c h h -=== （6） 2气体的状态方程采用BWRS 气体状态方程，如下：230000023436222()()()(1)exp()C D E d P RT B RT A bRT a T T T T d c a T Tρρρραργργρ=+--+-+--++++- （7）式中：P ——系统压力，KPa ； T ——系统温度，K ；ρ——混合气体密度，Kmol/m3； R ——气体常数，8.3143KJ/( Kmol .K)。

A 0、B 0、C 0、D 0、E 0、a 、b 、c 、d 、α、γ为方程的是一个参数，根据（8）确定。

1/21/2000110011/21/23000111/21/24000111/21/250001131/311/31(1)(1)(1)(1)n ni j i i ij i j ni ii nni j i i iji j nni j i i ij i j nni j i i ij i j n i i i n i i i A x x A A k B x B C x x C C k D x x D D k E x x E E k a x a b x b ============-==-=-=-⎡⎤=⎢⎥⎣⎦⎡⎤=⎢⎣⎦∑∑∑∑∑∑∑∑∑∑∑331/3131/3131/3131/31n i i i n i i i n i i i n i i i c x c d x d x x γγαα====⎥⎡⎤=⎢⎥⎣⎦⎡⎤=⎢⎥⎣⎦⎡⎤=⎢⎥⎣⎦⎡⎤=⎢⎥⎣⎦∑∑∑∑ （8）式中：x i 、x j ——混合气体中i 和j 组分的摩尔分数；k ij ——为i 、j 组分间的交互作用系数。

3气体的焓方程气体的焓方程如下：0000002342522422456(2)1417(23)(6)25(3)exp()]2C D E h h B RT A T T Td d bRT a a a T T c T ρρργργργργ=+--+++--+++--- （9） 4 管道周边的热力模型管道的有效土壤厚度采用等效圆筒法，传热半径由下式计算：20.521122((()1)1)H H R R R D D-=+-- （10）式中：R 2-R 1——土壤厚度，m ；R 1——从管道中心至土壤层的半径，m ； H ——至管道中心的实际埋深，m ； D ——管道直径，m 。

管道和周围环境的瞬态热力模型计算式如下：()/r r p t k rT r C T ρ= （11）式中：k ——周围环境导热系数；r ——传热半径； T r ——r 处的气体温度； C p ——气体定压比热； T t ——t 时刻的气体温度。

单位管长热流量由下式表示。

通过该公式计算管壁在任意节点的温度。

2012112()()ln(()/)w w k T T k D T T R R R πφπ-==-- （12）式中：k 2——管壁至土壤换热系数；K 1——流体至管壁换热系数； T w ——管壁温度； T 0——R 2处的温度； T ——气体温度。

5 水力摩阻系数计算式管段控制方程涉及的水力摩阻系数λ采用F.Colebrook －White 公式计算，该公式表达如下：101/ 1.73852log (2/))e e D R =-+ （13）式中，e/D ——管道粗糙度和内径的无因次比；R e ——雷诺数。

6控制方程的离散化由管道控制方程与气体状态方程组成的非线性偏微分方程组，一般不能得出管流气体基本变量的解析解，因此有必要应用计算数学的方法求解偏微分方程组的数值解。

本专题中选用中心隐式差分法对控制方程进行离散化。

确定采用的基本变量为气体的密度（ρ）、速度（u ）和温度（T ）。

6.1离散形式引进变量φ，φ代表三个流动基本中的任意一个。

在时间步长为Δx , 空间步长为Δt 的情况下，以空间i 和时间网格点t 采用中心隐式差分格式，则有以下离散形式：对于基本流动变量：11114k k k k i i i i φφφφφ+++++++=（14）基本流动变量对时间的一阶偏导数：11112k k k ki i i it tφφφφφ++++-+-∂=∂∆ （15）基本流动变量对空间的一阶偏导数：11112k k k ki i i ix xφφφφφ++++-+-∂=∂∆ （16）基本流动变量对时间的二阶偏导数：121212211122222(2)2(2)(2)16k k k k k k k k k i i i i i i i i i t tφφφφφφφφφφ++++++++++++-++-++-+∂=∂∆ （17） 基本流动变量对空间的二阶偏导数：111222221212122(2)2(2)(2)16k k k k k k k k k i i i i i i i i ix xφφφφφφφφφφ++++++++++++-++-++-+∂=∂∆ （18） 基本流动变量对空间及时间的二阶偏导数：2222216k k k k i i ii x t x tφφφφφ+++++--∂=∂∂∆∆ （19） 6.2 离散后的控制方程离散后的控制方程如下： 离散后的连续性方程：111111*********k k k k k k k k k k k k i i i i i i i i i i i iu u u u t xρρρρρρρρ++++++++++++-+--+-+=⨯∆⨯∆ （20）离散后的运动方程：11111111111111221122111111112111122()()()()2()()0244k k k k k k k k k k k k i i i i i i i i i i iik k k k k k k k i i i i i i i i k k k k k k k k i i i i i i i i u u u u P P P P t xu u u u xu u u u D ρρρρρρρρρρρρλ++++++++++++++++++++++++++++-+--+-+⨯∆⨯∆-+-+⨯∆+++++++⨯= （21） 离散后的能量方程：11111221111311111111111122111111111111-()-(-())22211()((2211-()-(-())222k k k k k k k k k k i i i i i i i i i i k k k k k k ki i i i i i i k k k k k k k k k k i i i i i i i i i i h P u h P u tu h u u h u h P u h P u tρρρρρρρρρρ++++++++++++++++++++++++++++++++++⨯∆+-++++⨯∆3111113311111))211()(())222(4)0k i k k k k k k kk i i i i i i i i k k k k i i i i w x u h u u h u xk T T T T T Dρρ+++++++++⨯∆+-++⨯∆++++-= （22）6.3 初始条件与边界条件初始条件指系统开始运行时的初始压力、流量或温度的分布状态。

边界条件指某一管段起始节点和终止节点上的约束条件。

主要包括：(1) 管段端点上的输油泵、压缩机或阀门等的出入口压力、流量、温度、转速、压比或开度设定值；(2) 气源对应节点的压力、流量或温度设定值； (3) 分输点对应节点的压力、流量或温度设定值； (4) 节点处压力、流量或温度的一致性； (5) 节点处压力、流量或温度的范围控制值；(6) 管道物理元件周围的温度场状况。

7非线性方程组的求解算法离散后的控制方程配合边界条件和初始条件才能封闭，封闭后形成了非线性方程组，对于该非线性方程组选取牛顿迭代法进行求解。

若采用()C x x b =的矩阵形式（其中C(x)为非线性方程组的系数矩阵），则123(,,,...,)T n x x x x x =为需要求解的向量，123(,,,...,)T n b b b b b =为等式右边的向量。

（1） 牛顿拉普森迭代法设迭代函数列123(,,,...,)Tn F F F F F = 迭代变量123(,,,...,)T n x x x x x = 迭代增量123(,,,...,)Tn x x x x x ∂=∂∂∂∂ 迭代函数123(,,,...,)i i n F F x x x x = 牛顿拉普森迭代公式如下：1k k k x x x δ+=- （5.7-1）对于迭代函数F ， 将求解非线性方程组问题转化成为寻根问题，也即要求下式成立：123(,,,...,)0T n F F F F F == （5.7-2）对任意点x 0和它的相邻点/邻域（x 0+△x ）， 通过泰勒展开式我们有：2001()()() 1,2,...,nii i j j F F x x F x x x i n xjδδοδ=∂+=++=∂∑ （5.7-3）若采用矩阵形式，则有：200()()() F x x F x J x x δδοδ+=++ （5.7-4）其中 J 为n ×n 的雅可比矩阵且iij jF J x ∂=∂。