CN110826188A - 一种基于gpu加速的天然气管网水力参数仿真方法 - Google Patents

Abstract

本发明公开了一种基于GPU加速的天然气管网水力参数仿真方法，首先对天然气管网的水力数学模型进行线性化和离散处理；将离散处理后每根管道的离散方程、非管元件的模型方程和水力参数边界条件联立后得到一个大型稀疏代数方程组；将大型稀疏代数方程组基于DIMENS算法进行分解，划分为若干个能并行计算的小型块三对角矩阵和1个小型稀疏矩阵；在图形处理器GPU上实现管道预求解方程组求解过程的并行；得到各管道离散方程的通解，将通解代入边界方程组使之能封闭求解；在完成所述边界方程组的求解后进行回代求解，获得该管道所有差分节点的水力参数。该方法具有计算速度快、效率高的优点，能够满足实时仿真的需求。

Description

一种基于GPU加速的天然气管网水力参数仿真方法

技术领域

本发明涉及天然气管网技术领域，尤其涉及一种基于GPU加速的天然气管网水力参数仿真方法。

背景技术

近年来，天然气消费市场的不断膨胀促使我国天然气工业进入了快速发展的黄金时期，而天然气管道作为天然气最主要的输送方式，也得到了大规模的建设。目前，我国新建天然气干线管道2.7万公里，使得我国天然气管道总长达到7.2万公里，已经初步建成一个四通八达的全国性管网，与此同时城镇燃气管道长度达到47.2万公里，而随着城市燃煤锅炉改造及大气污染治理等要求，城市燃气管网的规模将进一步扩大。天然气管网里程日益增加，各管道相互联通，管网拓扑结构不断趋于复杂化，任何管道或设备运行状态变化均会波及整个管网，这对管网的设计、运行、管理等均提出了新的挑战。若是能提前预测天然气管网内流动参数，则可以根据这些流动参数的变化趋势来制定相应运行方案，有可效应对那些新的挑战。在天然气流动参数中，压力、流量等水力参数非常重要，主要体现在以下几个方面：(1)管道中天然气流量测算的准确性决定了工程中天然气计量的准确性，而天然气计量的准确性直接关系到天然气企业的实际效益；(2) 天然气在管道中的流动压力影响着天然气的流动状态，它是驱动天然气流动、克服管道摩阻的重要因素，天然气的压力同时也关系到管网中各项设备的运行安全性，当天然气流动压力超过了设备的允许值时往往会造成设备损坏甚至爆炸等事故的产生。可见，水力参数的正确预测，可以为合理规划管网、优化调度运营、保障天然气供应等提供重要的参考依据，对天然气管网输送业务有重要的实际意义。

天然气管网中水力参数预测的方案是在天然气温度参数确定的情况下，通过求解天然气在管道内部流动的连续性方程和动量方程以及非管道元件(压缩机或阀门等)的流量、压力平衡方程，以得到管道内部的水力参数，这个过程称为天然气管网的水力仿真。天然气在管道内部流动的连续性方程和动量方程为非线性偏微分方程，无法直接给出解析解，因此工程上常采用数值方法来求解，典型的天然气管网水力数值解法的流程如下：首先建立起描述天然气管网气体流动的水力数学模型；由于该数学模型是非线性的，直接针对该模型进行求解的稳定性一般不太理想，因此第二步一般是将管网的水力模型进行线性化；接下来运用特征线法、隐式差分等方法对线性化后的水力模型进行离散得到差分方程组，再联立边界及节点平衡条件形成一个大型的稀疏线性方程组；最后结合一些常用的一些矩阵求解技术如Krylov子空间迭代法等进行整体求解。

上述现有技术的方法虽然能够实现天然气管网的水力仿真，但是对于大型的天然气管网而言，所形成的方程组规模通常会非常大，例如对于一个有1万个离散节点的天然气管网而言，方程组的规模将达到2万阶。在求解大型稀疏线性方程组时，若采用串行计算，常用的矩阵求解技术的求解效率仍然有限，尤其将其应用在天然气管网实时仿真中时，若管网流动状态变化剧烈导致选取的时间步长较小，往往求解计算的时间就超过了选取的时间步长，从而无法满足实时仿真的需求。

发明内容

本发明的目的是提供一种基于GPU加速的天然气管网水力参数仿真方法，该方法具有计算速度快、效率高的优点，能够满足实时仿真的需求，对天然气流动状态波动持续时间较长或波动较大的工况均可高效地进行水力参数预测。

本发明的目的是通过以下技术方案实现的：

一种基于GPU加速的天然气管网水力参数仿真方法，所述方法包括：

步骤1、对天然气管网的水力数学模型进行线性化和离散处理；

步骤2、将离散处理后每根管道的离散方程、非管元件的模型方程和水力参数边界条件联立后得到一个大型稀疏代数方程组；

步骤3、将所述大型稀疏代数方程组基于DIMENS算法进行分解，划分为若干个能并行计算的小型块三对角矩阵和1个小型稀疏矩阵；其中，每个小型块三对角矩阵为对应管道的预求解方程组，方程形式为A_iU_i-1+B_iU_i+C_iU_i+1＝D_i，U为天然气水力参数，A,B,C,D 为预求解方程组的系数矩阵，下标i代表该管道的第i个差分节点；小型稀疏矩阵为内部及外部边界条件共同构成的边界方程组；

步骤4、在图形处理器GPU上建立粗粒度线程映射模型，将每个管道的预求解方程组的求解任务映射至一个对应的GPU线程块Block上，接着建立细粒度线程映射模型，实现管道预求解方程组求解过程的并行；

步骤5、在管道预求解方程组完成计算后得到各管道离散方程的通解，将所述通解代入所述边界方程组使之能封闭求解；

步骤6、在完成所述边界方程组的求解后获得管网内所有管道的两端点水力参数U_0,k和 U_N+1,k，再将U_0,k和U_N+1,k回代至该管道的预求解方程组中，获得该管道所有差分节点的水力参数。

由上述本发明提供的技术方案可以看出，上述方法具有计算速度快、效率高的优点，能够满足实时仿真的需求，对天然气流动状态波动持续时间较长或波动较大的工况均可高效地进行水力参数预测，对天然气管网的运行管理有着重要的实际意义。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例提供的基于GPU加速的天然气管网水力参数仿真方法流程示意图；

图2为本发明实施例所述管道的离散示意图；

图3为本发明实施例所述粗粒度层的线程模型示意图；

图4为本发明实施例所述管道预求解方程组的示意图；

图5为本发明实施例所述PCR算法在GPU上对应的线程映射模型示意图；

图6为本发明实施例所述边界方程组的系数矩阵示意图；

图7为本发明所举实例的管网拓扑示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

下面将结合附图对本发明实施例作进一步地详细描述，如图1所示为本发明实施例提供的基于GPU加速的天然气管网水力参数仿真方法流程示意图，所述方法包括：

步骤1、对天然气管网的水力数学模型进行线性化和离散处理；

在该步骤中，对天然气管网的水力数学模型进行线性化的过程为：

对天然气管网的水力数学模型进行线性化是针对天然气在管道中流动的连续性方程和动量方程来实现，其中，连续性方程表示为：

动量方程表示为：

其中，m为质量流量，kg/s；p为压力，Pa；ρ为天然气密度，kg/m³；A为管道横截面面积，m²；d为管道内径，m；s为管道高程，m；θ为管道倾斜角，rad；λ为摩擦阻力系数；g为重力加速度，m/s²；T为气体温度，K；

和

分别表示了定温条件下压力随密度的变化率以及定密度条件下压力随温度的变化率，该变化率可通过天然气状态方程获得；

经过线性化处理后的连续性方程和动量方程表示为：

上述对天然气管网的水力数学模型进行离散处理的过程为：

首先将天然气管网中的每根管道划分成很多个小管段，每个小管段称为一个离散段；

将天然气管网中的非管道元件作为一个离散段，其离散方程表示为：

两个离散段间的连接点称为离散节点，将天然气管网用很多个具有编号的离散节点代替；其中，在完成管道的时间和空间离散后，将每根管道的线性化连续性方程和动量方程分别基于隐式差分法进行离散，离散后的管道离散方程表示为：

I为2×2的单位矩阵。

具体实现中，对离散节点的编号说明如下：

天然气管网中所有离散节点的编号形式为(i,j)，其中j为元件在天然气管网中的编号，i为离散节点在元件中的编号，则(i,j)即为元件j中的离散节点i；

天然气管网中元件的编号从1到M，天然气管网中编号为j的管道划分的离散段数目为 N_j，按管道j的起点到终点方向，将离散节点分别编号1到N_j+1；其中，i,j均为自然数；M为大于1的自然数；N_j为大于等于1的自然数。

上述天然气管网中的编号为j的压缩机和阀门等短小的非管道元件划分的离散段数目为N_j＝1，因此，将该类元件实际工作时天然气流进和流出位置所对应的离散节点分别编号为1和2，如图2所示为本发明实施例所述管道的离散示意图，参考图2：n为计算过程中时层的编号；n时层的天然气水力参数已知或已求解，n+1时层的天然气水里参数为待预测。

步骤2、将离散处理后每根管道的离散方程、非管元件的模型方程和水力参数边界条件联立后得到一个大型稀疏代数方程组；

在该步骤中，除了管道模型以外，天然气管网中还有阀门、压缩机等元件，前者为管道中天然气流动提供动力，后者控制气体流向，可统称为非管元件。阀门和压缩机这两类元件的种类繁多，结构和性能千差万别，而且不同种类或型号都有其特定的工作特性曲线，所以天然气经过非管道元件时水力参数变化的方程的具体数学表达式也不相同，该非管元件的模型方程表示为：

f(U_in,U_out)＝0

其中，f为数学公式的通用表达式，代表某种数学计算过程；U_in为非管道元件进口的水力参数组；U_out为非管道元件出口的水力参数组；

针对非管道元件模型进行如下线性化处理为：

进一步的，非管道元件线性化后的代数方程表示为：

式中，(1,j)和(2,j)为编号为j的非管道元件的进口和出口位置；U为天然气水力参数，包含压力p和质量流量m；上标n为计算过程中时层的编号，n时层的天然气水力参数已知或已求解，n+1时层的天然气水力参数为待预测。

除上述方程外，天然气管网中还包括水力参数边界条件，所述水力参数边界条件包括两类，第一类是已知天然气在供气源和分输点的气体流量或者压力值；第二类是内部连接点处的流量平衡条件和压力平衡条件，其数学表达式分别为：

第一类边界条件：p＝p(t),m＝m(t)

第二类边界条件：∑m_in＝∑m_out,p_i,1＝p_i,2＝…＝p_i,sum

其中，p为天然气压力，单位为Pa；m为天然气质量流量，单位为kg/s；p(t)和m(t)为t时刻供气源和分输点处的天然气压力和质量流量；下标in代表流入连接点，下标out 代表流出连接点；下标(i,x)中i代表第i个连接点，x代表与连接点相连的第x个元件， sum为与该连接点相连的元件总数。

通过求解上述大型稀疏代数方程组即能获得管网所有元件的水力参数。

步骤3、将所述大型稀疏代数方程组基于DIMENS算法进行分解，划分为若干个能并行计算的小型块三对角矩阵和1个小型稀疏矩阵；

其中，每个小型块三对角矩阵为对应管道的预求解方程组，方程形式为 A_iU_i-1+B_iU_i+C_iU_i+1＝D_i，U为天然气水力参数，A,B,C,D为预求解方程组的系数矩阵，下标i代表该管道的第i个差分节点；小型稀疏矩阵为内部及外部边界条件共同构成的边界方程组。

步骤4、在图形处理器GPU(Graphic Processing Unit)上建立粗粒度线程映射模型，将每个管道的预求解方程组的求解任务映射至一个对应的GPU线程块Block上，接着建立细粒度线程映射模型，实现管道预求解方程组求解过程的并行；

在该步骤中，CUDA(Compute Unified Device Architecture)是目前广泛使用的GPU 编程框架，它的核心是一种两层线程模型(线程块Block+线程Thread)，上述并行策略在GPU上执行时需要首先建立起相应的线程映射模型，即建立起GPU线程与计算任务之间的映射关系，如图3所示为本发明实施例所述粗粒度层的线程模型示意图，图3中每个管道的预求解方程组的求解任务映射至一个对应的线程块上，在CUDA中，线程块Block之间的计算是并发执行的，根据上述原理，在GPU上进行粗粒度并行策略的过程为：

首先对一个有M_p根管道、M_c个压缩机、M_v个阀门、M_e个气源(包括分输点)的管网，一共可以列出M_p个管道预求解方程组、2M_c个压缩机方程、2M_v个阀门方程、M_e个气源边界条件方程，2M_p+2M_c+2M_v+M_e个内部连接点边界条件方程；

向GPU申请M_p个线程块Block，并申请相应大小的共享内存(Shared Memory)用于存储对应编号管道的预求解方程组的系数矩阵A,B,C,D；

每个Block针对对应的管道预求解方程组分别执行细粒度层的并行计算，待计算完毕后释放计算资源并等待其它Block完成计算(即进行同步)；

将细粒度层并行计算得到的M_p个方程通解与其它方程联立组成一个规模为 4M_p+4M_c+4M_v+2M_e的边界方程组，然后根据实际数据规模对其进行GPU并行求解或者CPU串行求解；

再次向GPU申请M_p个线程块Block，将边界方程组求解得到的每根管道两端的水力参数存储在对应编号Block的共享内存上，并进行管道预求解方程组在细粒度层上的并行回代；

待所有Block完成计算后对计算结果进行处理后退出。

另外，如图4所示为本发明实施例所述管道预求解方程组的示意图，该方程组呈块三对角型，本实施例中使用PCR(并行循环规约)算法实现其求解过程的并行，该方法具有线程利用率和活跃度高的特点，与GPU的并行特点相匹配。如图5所示为本发明实施例所述PCR算法在GPU上对应的线程映射模型示意图，方程组中每个变量的计算都由1个线程完成，而线程间的计算是并发执行的，仅需在每次规约时进行线程同步，由PCR算法的原理可知，对规模为2N的管道预求解方程组，仅需次规约计算就能够获得所有变量的值，符号

代表向上取整，根据上述原理，在GPU上进行细粒度层并行计算的实施步骤如下：

对每个线程块Block，假设与之对应的管道共有N+1个差分节点，则该Block需要申请 N+1个线程，线程的编号从1～N+1，每个线程的寄存器中保存对应节点的方程系数 A_i,B_i,C_i,D_i；

基于并行循环规约PCR算法对该管道预求解方程进行第一次规约，规约过程中编号为 i的线程负责计算第i个差分节点的新方程系数A_i’,B_i’,C_i’,D_i’；规约的方法参考PCR算法原理，线程i计算完成后需与其它线程进行同步；

所有线程同步后即完成了一次规约，规约后形成两个方程组，分别由原方程中的0,2,4…N号差分节点的水力参数变量和1,3,4…N+1号差分节点的水力参数变量组成，两个方程组的规模相对于原方程组减半，接下来对这两个方程重复上述规约操作直至新形成的方程组只剩下一个变量；

对仅剩一个变量的方程B_iX_i＝D_i，由对应的i号线程直接进行计算从而获得X_i的值；

在得到所有X_i的值后释放计算资源，退出计算。

步骤5、在管道预求解方程组完成计算后得到各管道离散方程的通解，将所述通解代入所述边界方程组使之能封闭求解；

如图6所示为本发明实施例所述边界方程组的系数矩阵示意图，与管道预求解方程组不同的是，边界方程组的系数矩阵并不规则，求解此类方程组一般采用Krylov子空间迭代法，该方法中存在大量的矩阵向量乘的操作，这类操作具有较好的并发性，因此同样可以利用GPU进行加速。

步骤6、在完成所述边界方程组的求解后获得管网内所有管道的两端点水力参数U_0,k和 U_N+1,k，将U_0,k和U_N+1,k回代至该管道的预求解方程组中，获得该管道所有差分节点的水力参数。

其中，U代表水力参数组

下标k代表第k根管道，下标的第二个参数代表第i根管道的差分节点编号。

具体实现中，上述并行回代也是细粒度层的并行，同样可以利用GPU来完成，其求解实施步骤如下：

①对每个线程块Block，假设与之对应的管道共有N+1个差分节点，则该Block需要申请N+1个线程，线程的编号从0～N+1，每个线程的寄存器中保存对应节点的方程系数α_i,β_i,γ_i,该系数是通过前面PCR算法计算获得的；

②编号为i的线程负责计算

计算完毕后等待其它线程完成计算并进行同步；

③所有线程计算完毕后结束该Block，并释放计算资源。

另外，在获得管道所有差分节点的水力参数后，所述方法还包括：

由所得到的水力参数画出曲线图来描述和分析管道内部的天然气温度。

下面通过两个方面来说明上述方案的原理：

第一，上述方案可利用GPU实现并行计算的原理：原始DIMENS算法实现了管道与管道计算间的解耦，在边界方程组求解之前各管道可以独立计算获得管道预求解方程组的通解,相互之间互不影响，因此是可以并行求解的，这一层级的并行属于粗粒度层的并行，而在每根管道的预求解方程组的计算中，由于方程系数矩阵呈块三对角形，因此可以利用PCR算法实现块三对角形方程组的并行计算，这一层级的并行属于细粒度层的并行，此外，对边界方程组的求解同样可以使用CUDA官方提供的并行矩阵运算库结合Krylov子空间迭代法进行计算，在完成边界方程组的计算后可以将计算结果并行回代至各管道预求解方程组的通解表达式中，该部分的并行也属于细粒度层的并行。

第二，上述方案具有非常好的加速效果的原理：GPU具有海量线程的特点，因此要求算法的并行粒度要尽可能得小从而才能最大化地利用GPU的计算能力，本发明通过两层级的并行任务划分使得每个线程实际仅需要进行多次对应方程系数的计算和约简(对边界方程组的并行来说也是数次乘法以及加和运算)，因此并行计算的粒度非常小，从而能够获得十分显著的加速效果。

下面以具体的实例对上述方法的过程进行详细说明，如图7所示为本发明所举实例的管网拓扑示意图，参考图7，具体实施过程为：

1、天然气管网水力数学模型的线性化和离散，天然气管网水力数学模型的线性化主要是针对天然气在管道中流动的连续性方程和动量方程。其中连续性方程的形式如下：

动量方程的形式如下：

其中，m为质量流量，kg/s；p为压力，Pa；ρ为天然气密度，kg/m³；A为管道横截面面积，m²；d为管道内径，m；s为管道高程，m；θ为管道倾斜角，rad；λ为摩擦阻力系数；g为重力加速度，m/s²；T为气体温度，K；

和

分别表示了定温条件下压力随密度的变化率以及定密度条件下压力随温度的变化率，该变化率可通过天然气状态方程获得；

经过线性化后的连续性方程和动量方程可以统一写作如下的形式：

对所述图7中管网的管道按照1～5进行编号，编号顺序对应图中的l₁,l₂,l₃,l₄,l₅，接下来按照上述方法实施例对每根管道的上述数学模型进行离散，各管道编号和所划分的离散段数目如下表1所示：

表1管道编号及划分段数

管道编号	1	2	3	4	5
						对应图中标记	l<sub>1</sub>	l<sub>2</sub>	l<sub>3</sub>	l<sub>4</sub>	l<sub>5</sub>
划分段数	N<sub>1</sub>	N<sub>2</sub>	N<sub>3</sub>	N<sub>4</sub>	N<sub>5</sub>

在完成管道的时间和空间离散后，将每根管道的线性化连续性方程和动量方程分别在上述网格的基础上基于隐式差分法进行离散，离散后的方程形式如下：

I为2×2的单位矩阵。

除了管道模型以外，该天然气管网中还有阀门、压缩机等元件，前者为管道中天然气流动提供动力，后者控制气体流向，可统称为非管元件。按照前述说明将天然气在非管道元件(压缩机或阀门等)的水力参数变化方程写成以下的通式f(U_in,U_out)＝0，表示天然气在非管道元件进口处的水力参数与出口处的水力参数之间的联系；其中，f为数学公式的通用表达式，代表某种数学计算过程。U_in为非管道元件进口的水力参数组； U_out为非管道元件出口的水力参数组。

对所述管网中的非管元件按照6～7进行编号，编号顺序对应图中的c₁,v₁。针对非管道元件j，进行如下线性化处理，

非管道元件线性化后的代数方程可整理为：

式中，f为数学公式的通用表达式，代表某种数学计算过程；(1,j)和(2,j)为编号为j的非管道元件的进口和出口位置；U为天然气水力参数，包含压力p和质量流量 m；上标“n”为计算过程中时层的编号；n时层的天然气水力参数已知或已求解，n+1时层的天然气水力参数为待预测。

除了上述方程以外，天然气管网中还包括了水力参数边界条件，这些边界条件分为两类，第一类是已知天然气在供气源和分输点的气体流量或者压力值，第二类是内部连接点处的流量平衡条件和压力平衡条件。对所述管网中的气源(分输点)按照8～11进行编号，编号顺序对应图中的e₁,e₂,e₃,e₄，其中气源1为压力控制，其它气源为流量控制，则第一类边界条件包括：

气源1的压力边界条件：p₈＝p₈(t)

气源2的流量边界条件：m₉＝m₉(t)

气源3的流量边界条件：m₁₀＝m₁₀(t)

分输点4的流量边界条件：m₁₁＝m₁₁(t)

第二类边界条件包括：

管道1起点处的流量平衡方程和压力平衡方程：m₈＝m_(1,1),p₈＝p_(1,1)

管道1终点处的流量平衡方程和压力平衡方程：

管道2起点处的流量平衡方程和压力平衡方程：m_(2,5)＝m_(1,2),p_(2,5)＝p_(1,2)

管道2终点处(管道4终点处、管道5起点处)的流量平衡方程和压力平衡方程：

管道3起点处的流量平衡方程和压力平衡方程：m₉＝m_(1,3),p₉＝p_(1,3)

管道3终点处的流量平衡方程和压力平衡方程：

管道4起点处的流量平衡方程和压力平衡方程：m_(2,6)＝m_(1,4),p_(2,6)＝p_(1,4)

管道5终点处的流量平衡方程和压力平衡方程：

其中，p为天然气压力，单位为Pa；m为天然气质量流量，单位为kg/s；p(t)和m(t)为t时刻供气源和分输点处的天然气压力和质量流量；下标(i,x)代表第x个元件的第i 个节点。

将每根管道的离散方程、非管元件的模型方程和管网边界条件联立后可以得到一个大型的稀疏代数方程组，方程组规模为2(N₁+N₂+N₃+N₄)+26，求解该方程组即能获得管网所有元件的水力参数。

2、天然气管道计算任务的解耦即粗粒度层并行：对上述大型稀疏代数方程组，基于DIMENS算法进行任务分解，原始的方程组被拆分为5个块三对角矩阵和1个小型稀疏矩阵，每个块三对角矩阵即为对应管道的预求解方程组，方程形式为 A_iU_i-1+B_iU_i+C_iU_i+1＝D_i,下标i代表该管道的第i个差分节点，5个块三对角矩阵的规模分别为2(N₁+1),2(N₂+1),2(N₃+1),2(N₄+1),2(N₅+1)；小型稀疏矩阵则为内部及外部边界条件共同构成的边界方程组，方程组规模为36。其中管道预求解方程组计算是可以并行执行的，这种按照管道为单位进行计算任务的划分和并行的过程即为天然气管网水力仿真粗粒度层的并行策略。

粗粒度层并行策略的实施步骤如下：

①对上述管网一共可以列出5个管道预求解方程组,2个压缩机方程，2个阀门方程，4个气源边界条件方程，18个内部连接点边界条件方程。

②向GPU申请5个线程块Block，编号分别为Block₁,Block₂,Block₃,Block₄,Block₅,并申请相应大小的共享内存(Shared Memory)用于存储对应编号管道的预求解方程组的系数矩阵A,B,C,D。

③每个Block针对对应的管道预求解方程组分别执行细粒度层的并行计算，待计算完毕后释放计算资源并等待其它Block完成计算(即进行同步)。

④将细粒度层并行计算得到的5个方程通解与其它方程联立组成一个规模为36的边界方程组，然后根据实际数据规模对其进行GPU并行求解或者CPU串行求解。

⑤再次向GPU申请5个线程块Block，编号分别为Block₁,Block₂,Block₃,Block₄,Block₅,将边界方程组求解得到的每根管道两端的水力参数存储在对应编号Block的共享内存上，并进行管道预求解方程组在细粒度层上的并行回代。

⑥待所有Block完成计算后对计算结果进行处理后退出。

3、管道预求解方程组求解过程的并行即细粒度层并行：管道预求解方程组的形式呈块三对角型，可以用方程A_iU_i-1+B_iU_i+C_iU_i+1＝D_i进行表示。本发明使用PCR(并行循环规约)算法实现其求解过程的并行，方程组中每个变量的计算都由1个线程完成，而线程间的计算是并发执行的，仅需在每次规约时进行线程同步。由PCR算法的原理可知，对规模为2N的管道预求解方程组，仅需

次规约计算就能够获得所有变量的值，根据上述原理，细粒度层并行策略的实施步骤如下：

①对每个线程块Block_k，下标k代表线程块的编号，与之对应的管道共有N_k+1个差分节点，则Block_k需要申请N_k+1个线程，线程的编号从1～N_k+1，每个线程的寄存器中保存对应节点的方程系数A_i,B_i,C_i,D_i。

②对该管道预求解方程进行第一次规约，规约过程中编号为i的线程主要负责计算第i个差分节点的新方程系数A_i’,B_i’,C_i’,D_i’，规约的方法参考PCR算法原理，线程i计算完成后需与其它线程进行同步。

③所有线程同步后即完成了一次规约，规约后形成两个方程组，分别由原方程中的 0,2,4…N_k号差分节点的水力参数变量和1,3,4…N_k+1号差分节点的水力参数变量组成，两个方程组的规模相对于原方程组减半，接下来对这两个方程重复步骤②中的规约操作直至新形成的方程组只剩下一个变量。

④对仅剩一个变量的方程B_iX_i＝D_i,可由对应的i号线程直接进行计算从而获得X_i的值。

⑤在得到所有X_i的值后释放计算资源，退出计算。

4、管网边界方程组的求解：管道预求解方程组完成计算后将得到各管道离散方程的通解，将这些通解代入管网边界方程组使之能够封闭求解。

5、管道预求解方程的并行回代求解：在完成边界方程组的求解后可以获得管网内所有管道的两端点水力参数U_0,k和U_N+1,k，U代表水力参数组

下标k代表第k根管道，下标的第一个参数代表第k根管道的差分节点编号。接下来将U_0,k和U_N+1,k回代至该管道的预求解方程组A_iU_i-1+B_iU_i+C_iU_i+1＝D_i的通解表达式中获得该管道的所有差分节点的水力参数，该回代过程同样是可以在GPU上进行并行的，其求解实施步骤如下：

①对线程块Block_k，假设与之对应的管道共有N_k+1个差分节点，则Block_k需要申请N_k+1个线程，线程的编号从0～N_k+1，每个线程的寄存器中保存对应节点的方程系数α_i,β_i,γ_i,该系数是通过前面PCR算法计算获得的。

②编号为i的线程负责计算计算完毕后等待其它线程完成计算并进行同步。

③所有线程计算完毕后结束该Block，并释放计算资源。

6、结果展示：由所得到结果画出曲线图来描述和分析管道内部的天然气温度。

值得注意的是，本发明实施例中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (10)

1.一种基于GPU加速的天然气管网水力参数仿真方法，其特征在于，所述方法包括：

步骤1、对天然气管网的水力数学模型进行线性化和离散处理；

步骤2、将离散处理后每根管道的离散方程、非管元件的模型方程和水力参数边界条件联立后得到一个大型稀疏代数方程组；

步骤3、将所述大型稀疏代数方程组基于DIMENS算法进行分解，划分为若干个能并行计算的小型块三对角矩阵和1个小型稀疏矩阵；其中，每个小型块三对角矩阵为对应管道的预求解方程组，方程形式为A_iU_i-1+B_iU_i+C_iU_i+1＝D_i，U为天然气水力参数，A,B,C,D为预求解方程组的系数矩阵，下标i代表该管道的第i个差分节点；小型稀疏矩阵为内部及外部边界条件共同构成的边界方程组；

步骤4、在图形处理器GPU上建立粗粒度线程映射模型，将每个管道的预求解方程组的求解任务映射至一个对应的GPU线程块Block上，接着建立细粒度线程映射模型，实现管道预求解方程组求解过程的并行；

步骤5、在管道预求解方程组完成计算后得到各管道离散方程的通解，将所述通解代入所述边界方程组使之能封闭求解；

步骤6、在完成所述边界方程组的求解后获得管网内所有管道的两端点水力参数U_0,k和U_N+1,k，再将U_0,k和U_N+1,k回代至该管道的预求解方程组中，获得该管道所有差分节点的水力参数。

2.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在步骤1中，所述对天然气管网的水力数学模型进行线性化的过程为：

对天然气管网的水力数学模型进行线性化是针对天然气在管道中流动的连续性方程和动量方程来实现，其中，连续性方程表示为：

动量方程表示为：

其中，m为质量流量；p为压力；ρ为天然气密度；A为管道横截面面积；d为管道内径；s为管道高程；θ为管道倾斜角；λ为摩擦阻力系数；g为重力加速度；T为气体温度；

和

分别表示定温条件下压力随密度的变化率以及定密度条件下压力随温度的变化率；

经过线性化处理后的连续性方程和动量方程表示为：

3.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在步骤1中，所述对天然气管网的水力数学模型进行离散处理的过程为：

首先将天然气管网中的每根管道划分成很多个小管段，每个小管段称为一个离散段；

将天然气管网中的非管道元件作为一个离散段，其离散方程表示为：

两个离散段间的连接点称为离散节点，将天然气管网用很多个具有编号的离散节点代替；其中，管道的离散方程表示为：

I为2×2的单位矩阵。

4.根据权利要求3所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，所述离散节点的编号具体为：

天然气管网中所有离散节点的编号形式为(i,j)，其中j为元件在天然气管网中的编号，i为离散节点在元件中的编号，则(i,j)即为元件j中的离散节点i；

天然气管网中元件的编号从1到M，天然气管网中编号为j的管道划分的离散段数目为N_j，按管道j的起点到终点方向，将离散节点分别编号1到N_j+1；

其中，i,j均为自然数；M为大于1的自然数；N_j为大于等于1的自然数。

5.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在步骤2中，所述非管元件包括压缩机和阀门，该非管元件的模型方程表示为：

f(U_in,U_out)＝0

其中，f为数学公式的通用表达式，代表某种数学计算过程；U_in为非管道元件进口的水力参数组；U_out为非管道元件出口的水力参数组；

针对非管道元件模型进行如下线性化处理为：

进一步的，非管道元件线性化后的代数方程表示为：

式中，(1,j)和(2,j)为编号为j的非管道元件的进口和出口位置；U为天然气水力参数，包含压力p和质量流量m；上标n为计算过程中时层的编号，n时层的天然气水力参数已知或已求解，n+1时层的天然气水力参数为待预测。

6.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在步骤2中，所述水力参数边界条件包括两类，第一类是已知天然气在供气源和分输点的气体流量或者压力值；第二类是内部连接点处的流量平衡条件和压力平衡条件，其数学表达式分别为：

第一类边界条件：p＝p(t),m＝m(t)

第二类边界条件：∑m_in＝∑m_out,p_i,1＝p_i,2＝…＝p_i,sum

其中，p为天然气压力；m为天然气质量流量；p(t)和m(t)为t时刻供气源和分输点处的天然气压力和质量流量；下标in代表流入连接点，下标out代表流出连接点；下标(i,x)中i代表第i个连接点，x代表与连接点相连的第x个元件；sum为与该连接点相连的元件总数。

7.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在步骤4中，在GPU上建立粗粒度线程映射模型的过程为：

对一个有M_p根管道、M_c个压缩机、M_v个阀门、M_e个气源的管网，一共可以列出M_p个管道预求解方程组、2M_c个压缩机方程、2M_v个阀门方程、M_e个气源边界条件方程，2M_p+2M_c+2M_v+M_e个内部连接点边界条件方程；

向GPU申请M_p个线程块Block，并申请相应大小的共享内存用于存储对应编号管道的预求解方程组的系数矩阵A,B,C,D；

每个Block针对对应的管道预求解方程组分别执行细粒度层的并行计算，待计算完毕后释放计算资源并等待其它Block完成计算；

将细粒度层并行计算得到的M_p个方程通解与其它方程联立组成一个规模为4M_p+4M_c+4M_v+2M_e的边界方程组，然后根据实际数据规模对其进行GPU并行求解或者CPU串行求解；

再次向GPU申请M_p个线程块Block，将边界方程组求解得到的每根管道两端的水力参数存储在对应编号Block的共享内存上，并进行管道预求解方程组在细粒度层上的并行回代；

待所有Block完成计算后对计算结果进行处理后退出。

8.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在步骤4中，在GPU上进行细粒度层并行计算的实施步骤如下：

对每个线程块Block，假设与之对应的管道共有N+1个差分节点，则该Block需要申请N+1个线程，线程的编号从1～N+1，每个线程的寄存器中保存对应节点的方程系数A_i,B_i,C_i,D_i；

基于并行循环规约PCR算法对该管道预求解方程进行第一次规约，规约过程中编号为i的线程负责计算第i个差分节点的新方程系数A_i’,B_i’,C_i’,D_i’；

所有线程同步后即完成了一次规约，规约后形成两个方程组，分别由原方程中的0,2,4…N号差分节点的水力参数变量和1,3,4…N+1号差分节点的水力参数变量组成，两个方程组的规模相对于原方程组减半，接下来对这两个方程重复上述规约操作直至新形成的方程组只剩下一个变量；

对仅剩一个变量的方程B_iX_i＝D_i，由对应的i号线程直接进行计算从而获得X_i的值；

在得到所有X_i的值后释放计算资源，退出计算。

9.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在步骤5中，在对边界方程组进行求解过程中，具体设定一个阈值δ，当所述边界方程组规模小于δ阶时，所述边界方程组的求解在CPU上执行；当所述边界方程组规模大于δ阶时，所述边界方程组的求解在GPU上执行。

10.根据权利要求1所述基于GPU加速的天然气管网水力参数仿真方法，其特征在于，在获得管道所有差分节点的水力参数后，所述方法还包括：

由所得到的水力参数画出曲线图来描述和分析管道内部的天然气温度。

Images