天然气管道运行模拟及仿真技术研究

天然气管道运行模拟及仿真技术研究1011202045蔡永军科学计算选讲结课论文为了预测天然气管道运行状态，制定合理的管输计划，更好的配置设备开机，天然气管道输送过程中需要进行工况模拟及仿真。实际工作中需要建立压缩机、阀门等设备的模型，确定管段的控制方程、气体的状态方程，针对给出的初始条件和边界条件，筛选确定天然气管网数学模型的离散方法与非线性方程组的求解算法寻找合理的非线性方程的求解算法，得到合理的数值解。1天然气管道仿真数学模型1.1管段的控制方程对于管道中的任意管段，经过适当的简化可以用下列公式来描述：连续性方程：()0AAtx（1）运动方程：2()(.)sin()2APAAAgAtxxD（2）能量方程：221(())(.())22sin()()WPhAhPAAAgDkTTtxx（3）式中：A——管道的横截面积，m2；ρ——流体密度，kg/m3；t——时间，s；x——坐标，m；u——速度，m/s；P——压力，Pa；θ——管道倾角，rad；λ——水力摩阻系数；D——管道内径，m；T——流体温度，k；k1——流体至管壁的换热系数；h——比焓；Tw——管壁的温度，k。1.2阀门控制方程阀门控制方程如下：120()0dwupgupdwdwMMMCPPPhh（4）式中：Mup——阀门入口质量流量，kg/s；Mdw——阀门入口质量流量，kg/sCg——阀门系数；Pup——阀的入口压力，Pa；Pdw——阀的出口压力，Pa。1.3压缩机控制方程简化后的压缩机控制方程如下222111001()()dwupfuelmmdwupnnabQcQnnMMMTT（5）式中：——压缩机压比；m——多变压缩指数；n——压缩机的实际转速，rpm；n0——压缩机的额定转速，rpm；a1,b1,c1——系数；Q——给定状态下的体积流量，m3/s；1.4理想调节阀阀控制方程理想调节阀控制方程如下：120dwupdwMMPchh（6）2气体的状态方程采用BWRS气体状态方程，如下：230000023436222()()()(1)exp()CDEdPRTBRTAbRTaTTTTdcaTT（7）式中：P——系统压力，KPa；T——系统温度，K；ρ——混合气体密度，Kmol/m3；R——气体常数，8.3143KJ/(Kmol.K)。A0、B0、C0、D0、E0、a、b、c、d、α、γ为方程的是一个参数，根据（8）确定。1/21/2000110011/21/23000111/21/24000111/21/250001131/311/31(1)(1)(1)(1)nnijiiijijniiinnijiiijijnnijiiijijnnijiiijijniiiniiiAxxAAkBxBCxxCCkDxxDDkExxEEkaxabxb331/3131/3131/3131/31niiiniiiniiiniiicxcdxdxx（8）式中：xi、xj——混合气体中i和j组分的摩尔分数；kij——为i、j组分间的交互作用系数。3气体的焓方程气体的焓方程如下：0000002342522422456(2)1417(23)(6)25(3)exp()]2CDEhhBRTATTTddbRTaaaTTcT（9）4管道周边的热力模型管道的有效土壤厚度采用等效圆筒法，传热半径由下式计算：20.521122((()1)1)HHRRRDD（10）式中：R2-R1——土壤厚度，m；R1——从管道中心至土壤层的半径，m；H——至管道中心的实际埋深，m；D——管道直径，m。管道和周围环境的瞬态热力模型计算式如下：()/rrptkrTrCT（11）式中：k——周围环境导热系数；r——传热半径；Tr——r处的气体温度；Cp——气体定压比热；Tt——t时刻的气体温度。单位管长热流量由下式表示。通过该公式计算管壁在任意节点的温度。2012112()()ln(()/)wwkTTkDTTRRR（12）式中：k2——管壁至土壤换热系数；K1——流体至管壁换热系数；Tw——管壁温度；T0——R2处的温度；T——气体温度。5水力摩阻系数计算式管段控制方程涉及的水力摩阻系数λ采用F.Colebrook－White公式计算，该公式表达如下：101/1.73852log(2/18.574/(*))eeDR（13）式中，e/D——管道粗糙度和内径的无因次比；Re——雷诺数。6控制方程的离散化由管道控制方程与气体状态方程组成的非线性偏微分方程组，一般不能得出管流气体基本变量的解析解，因此有必要应用计算数学的方法求解偏微分方程组的数值解。本专题中选用中心隐式差分法对控制方程进行离散化。确定采用的基本变量为气体的密度（ρ）、速度（u）和温度（T）。6.1离散形式引进变量φ，φ代表三个流动基本中的任意一个。在时间步长为Δx,空间步长为Δt的情况下，以空间i和时间网格点t采用中心隐式差分格式，则有以下离散形式：对于基本流动变量：11114kkkkiiii（14）基本流动变量对时间的一阶偏导数：11112kkkkiiiitt（15）基本流动变量对空间的一阶偏导数：11112kkkkiiiixx（16）基本流动变量对时间的二阶偏导数：121212211122222(2)2(2)(2)16kkkkkkkkkiiiiiiiiitt（17）基本流动变量对空间的二阶偏导数：111222221212122(2)2(2)(2)16kkkkkkkkkiiiiiiiiixx（18）基本流动变量对空间及时间的二阶偏导数：2222216kkkkiiiixtxt（19）6.2离散后的控制方程离散后的控制方程如下：离散后的连续性方程：111111111111022kkkkkkkkkkkkiiiiiiiiiiiiuuuutx（20）离散后的运动方程：11111111111111221122111111112111122()()()()2()()0244kkkkkkkkkkkkiiiiiiiiiiiikkkkkkkkiiiiiiiikkkkkkkkiiiiiiiiuuuuPPPPtxuuuuxuuuuD（21）离散后的能量方程：11111221111311111111111122111111111111-()-(-())22211()((2211-()-(-())222kkkkkkkkkkiiiiiiiiiikkkkkkkiiiiiiikkkkkkkkkkiiiiiiiiiihPuhPutuhuuhuhPuhPut3111113311111))211()(())222(4)0kikkkkkkkkiiiiiiiikkkkiiiiwxuhuuhuxkTTTTTD（22）6.3初始条件与边界条件初始条件指系统开始运行时的初始压力、流量或温度的分布状态。边界条件指某一管段起始节点和终止节点上的约束条件。主要包括：(1)管段端点上的输油泵、压缩机或阀门等的出入口压力、流量、温度、转速、压比或开度设定值；(2)气源对应节点的压力、流量或温度设定值；(3)分输点对应节点的压力、流量或温度设定值；(4)节点处压力、流量或温度的一致性；(5)节点处压力、流量或温度的范围控制值；(6)管道物理元件周围的温度场状况。7非线性方程组的求解算法离散后的控制方程配合边界条件和初始条件才能封闭，封闭后形成了非线性方程组，对于该非线性方程组选取牛顿迭代法进行求解。若采用()Cxxb的矩阵形式（其中C(x)为非线性方程组的系数矩阵），则123(,,,...,)Tnxxxxx为需要求解的向量，123(,,,...,)Tnbbbbb为等式右边的向量。（1）牛顿拉普森迭代法设迭代函数列123(,,,...,)TnFFFFF迭代变量123(,,,...,)Tnxxxxx迭代增量123(,,,...,)Tnxxxxx迭代函数123(,,,...,)iinFFxxxx牛顿拉普森迭代公式如下：1kkkxxx（5.7-1）对于迭代函数F，将求解非线性方程组问题转化成为寻根问题，也即要求下式成立：123(,,,...,)0TnFFFFF（5.7-2）对任意点x0和它的相邻点/邻域（x0+△x），通过泰勒展开式我们有：2001()()()1,2,...,niiijjFFxxFxxxinxj（5.7-3）若采用矩阵形式，则有：200()()()FxxFxJxx（5.7-4）其中J为n×n的雅可比矩阵且iijjFJx。如果略去其中的高次项2()x，并要求0()0Fxx。我们得出：1xJF（5.7-5）至此，可以按照牛顿拉普森迭代法的求解步骤进行计算。（2）牛顿＋最速搜索迭代法引入目标函数0.5fFF，通过简单的数学运算，可以得到牛顿迭代沿着此目标函数的梯度方向f，始终可以发现一个数值α能使得目标函数的值下降，即：1()0fxFJJFFF。所以本专题研究采用如下迭代公式：1kkkxxx（5.7-6）该方法较牛顿拉普森方法具有收敛速度快，且全局收敛的特点。8仿真运行8.1仿真运算的基本过程仿真运行的基本过程如图1仿真运行所示。运行系统状态数据和不同视角构成单文档－多视的关系。图中的兰色带箭头线条表示通过不同的视角和核心功能层接口，可以监视或编辑仿真系统的状态数据，并下达计算命令；粉色线条表示运算中仿真模拟器和数据模块进行数据交互；绿色线条表示初始化过程加载数据。简单人机界面调度模块编译模块数据模块视角视角表格视角视角视角命令行视角编连文件运行系统状态数据加载系统及其初始参数监视或编辑状态数据仿真模拟器计算状态数据核心功能层接口图1仿真运行过程仿真运行基本过程如下：(1)通过各个接口或视角，下达加载系统命令；(2)调度模块命令编译模块加载编连文件；(3)如果需要，通过各个接口或视角，对系统参数进行进一步初始化；(4)通过各个接口或视角，下达一轮计算命令；(5)仿真模拟器开始一轮计算，并输出结果到数据模块；调度模块通知相应接口和视角计算结束；(6)相应接口或视角获取关心的数据；(7)如果需要进行新一轮计算，相应接口和视角可以对部分参数进行重新设定，并下达新一轮计算命令，系统将回到第5步。8.2仿真运算的实现结构仿真运算由仿真模拟器作为核心模块来实现，两个直接的辅助模块是数据模块和编译模块。这3个模块相互协同进行仿真运算，基本过程如图2所示。编译模块仿真模拟器编连文件模型文件元件库控制方程组建立此系统的联立组1输出：方程组＋元件信息数据模块2加载数据控制方程组2加载方程组输入参数基本输出参数其它输出参数3数值计算45元件参数关系计算函数库其它输出参数数值化计算6476输入参数预处理46元件参数关系计算函数库外部边界方程组外部边界方程组