技术领域
[0001] 本
发明涉及
流体力学技术领域，具体是涉及一种基于
天然气管道传输问题的无网格粒子方法。
背景技术
[0002] 针对天然气管道气体的输运中，由于气体的可压缩性强，所以针对气体的计算和模拟和其他流体有明显的不同，对于这一问题，常用的方法主要有基于特征线法的界面追踪技术和有限体积法等。为了降低气态垂直界面假设引起的计算误差，通过假定气体状态在单个网格距离内是线性变化的，提出了改进的界面追踪技术，但这影响了计算的效率和
稳定性，也可能因插值误差的积累而降低计算
精度，针对以上问题，本发明提出了一种基于拉格朗粒子模型的无网格粒子方法。
发明内容
[0003] 为克服
现有技术的不足，本发明旨在提出提供一种求解天然气管道传输问题的方法。为此，本发明采取的技术方案是，预测管道气锤的方法，考虑气体的强可压缩性带来的影响，采用光滑粒子流体动力学方法求解移动
坐标系统下的激波管物理方程，模拟气体对管道的冲击。
[0004] 具体步骤如下：
[0005] 步骤一，初始化:初始化系统的相关变量和粒子信息；
[0006] 步骤二，生成粒子信息；
[0008] 根据天然气管道问题的原理，得出物理模型的数学模型，即拉格朗日形式下的控制方程为：
[0009]
[0010]
[0011]
[0012] 其中P是气体压力,ν是气体的速度，ρ是气体的
密度，e是气体的
能量，对于理想气体，其
状态方程为：
[0013] p＝(γ-1)ρe (4)
[0014] 将问题简化，相应的状态方程为：
[0015]
[0016] 此时，能量项将从控制方程中解耦出来；
[0017] 为了求解上述方程(5)，从而预测天然气管道中气锤问题，采取步骤如下：
[0018] 光滑粒子流体动力学方法中，函数f(x)的积分表示式为：
[0019] f(x)＝∫Ωf(x′)δ(x-x′)dx′ (6)
[0020] 其中δ(x-x′)为狄拉克函数，Ω为包含x的积分体积，若用光滑函数W(x-x′,h)取代狄拉克函数，则f(x)的积分表示式为：
[0021] f(x)≈∫Ωf(x′)W(x-x′,h)dx′ (7)
[0022] 函数导数积分表示为：
[0023]
[0024] 又因为
[0025]
[0026] 所以
[0027]
[0028] 然后进行粒子近似得到：
[0029]
[0030] 因为
[0031]
[0032] 所以
[0033]
[0034] 又有
[0035]
[0036] 所以，在粒子i处的函数的粒子近似式写为：
[0037]
[0038] 利用光滑粒子流体动力学方法将方程(1)～(3)离散可得(16)～(18)式[0039]
[0040]
[0041]
[0042] p＝(γ-1)ρe (19)
[0043]
[0044] 其中∏ij为Monaghon型人工
粘度，Hij为人工热量项，γ、ξ为常系数；
[0045] 根据方程(16)、(17)和(18)计算不同时刻每个粒子的密度、速度和能量信息，进而根据(19)式求解出粒子的压强；
[0046] 步骤四，输出结果：
[0047] 1)每一个时间步结束，保存中间结果，并输出中间结果；
[0048] 2)结束时间循环，并输出最终结果。
[0049] 根据方程(16)、(17)和(18)计算不同时刻每个粒子的密度、速度和能量信息，进而根据(19)式求解出粒子的压强，具体计算过程为：
[0050] 1)循环每个时间步；
[0051] 2)初始化计算域内所有粒子
属性信息后，搜索目标粒子的邻近粒子，获得初始值的速度和
质量，然后通过方程(16)计算该粒子的密度的导数，最后通过欧拉时间积分更新粒子的密度信息；
[0052] 3)通过上一步得到的更新的密度信息、初始设置的压力值和Monaghon型人工粘度方程计算得到速度导数，然后通过欧拉时间积分更新粒子的速度信息，最后通过速度信息更新粒子的
位置；
[0053] 4)根据上面两步得到的密度、速度粒子信息和Monaghon型人工粘度方程及人工热量项方程计算得到能量导数，然后通过欧拉时间积分更新粒子的能量信息；
[0054] 5)通过上面得到的密度和能量信息结合方程(19)，更新粒子的压力；
[0055] 6)然后通过(20)式更新粒子搜索的光滑长度信息。
[0056] 本发明的特点及有益效果是：
[0057] 本发明采用光滑粒子流体动力学方法求解移动坐标系统下的激波管物理方程，充分考虑了气体的强可压缩性带来的影响，在满足数值精度的前提下能够更方便地模拟气体对管道的冲击问题。
附图说明：
[0059] 图2天然气管道左边冲击右边扩张问题物理模型。
[0060] 图3 CASE1中天然气管道在0.6s时刻密度、压强、速度和能量随位移的分布图。
[0061] 图4一维天然气管道气锤问题物理模型。
[0062] 图5 CASE2中天然气管道中间点压力和密度随时间分布。
具体实施方式
[0063] 本发明解决的技术问题是提供一种求解天然气管道传输问题的拉格朗日粒子方法，方法采用光滑粒子流体动力学方法求解移动坐标系统下的激波管物理方程，充分考虑了气体的强可压缩性带来的影响，在满足数值精度的前提下能够更方便地模拟气体对管道的冲击问题。
[0064] 步骤一，初始化:初始化系统的相关变量和粒子信息；
[0065] 步骤二，生成粒子信息；
[0066] 步骤三，列出求解方程并迭代计算：
[0067] 根据天然气管道问题的原理，可以得出物理模型的数学模型，即拉格朗日形式下的控制方程为：
[0068]
[0069]
[0070]
[0071] 其中P是气体压力,ν是气体的速度，ρ是气体的密度，e是气体的能量。(1)～(3)中的拉格朗日控制方程是未闭合的，需要添加状态方程来保证系统完
[0072] 整性。对于理想气体，其状态方程为：
[0073] p＝(γ-1)ρe (4)
[0074] 对于实际管道中的气锤问题，如果
温度变化不大，可以将问题简化，相应的状态方程为：
[0075]
[0076] 此时，能量项将从控制方程中解耦出来。
[0077] 为了求解上述方程，从而预测天然气管道中气锤问题，本发明提出了一种拉格朗日粒子方法。具体内容如下：
[0078] 光滑粒子流体动力学方法中，函数f(x)的积分表示式为：
[0079] f(x)＝∫Ωf(x′)δ(x-x′)dx′ (6)
[0080] 其中δ(x-x′)为狄拉克函数，Ω为包含x的积分体积。若用光滑函数W(x-x′,h)取代狄拉克函数，则f(x)的积分表示式为：
[0081] f(x)≈∫Ωf(x′)W(x-x′,h)dx′ (7)
[0082] 函数导数积分表示为：
[0083]
[0084] 又因为
[0085]
[0086] 所以
[0087]
[0088] 然后进行粒子近似得到：
[0089]
[0090] 因为
[0091]
[0092] 所以
[0093]
[0094] 又有
[0095]
[0096] 所以，在粒子i处的函数的粒子近似式可写为：
[0097]
[0098] 利用光滑粒子流体动力学方法将方程(1)～(3)离散可得(16)～(18)式[0099]
[0100]
[0101]
[0102] p＝(γ-1)ρe (19)
[0103]
[0104] 其中Πij为Monaghon型人工粘度，Hij为人工热量项，γ、ξ为常系数。
[0105] 根据方程(16)、(17)和(18)计算不同时刻每个粒子的密度、速度和能量信息，进而根据(19)式求解出粒子的压强。具体计算过程如下所示：
[0106] 具体计算过程为：
[0107] 6)循环每个时间步；
[0108] 7)初始化计算域内所有粒子属性信息后，搜索目标粒子的邻近粒子，获得初始值的速度和质量，然后通过方程(16)计算该粒子的密度的导数，最后通过欧拉时间积分更新粒子的密度信息；
[0109] 8)通过上一步得到的更新的密度信息、初始设置的压力值和Monaghon型人工粘度方程计算得到速度导数，然后通过欧拉时间积分更新粒子的速度信息，最后通过速度信息更新粒子的位置；
[0110] 9)根据上面两步得到的密度、速度粒子信息和Monaghon型(由Monaghon在1992年提出解决数值震荡的方法)人工粘度方程及人工热量项方程计算得到能量导数，然后通过欧拉时间积分更新粒子的能量信息；
[0111] 10)通过上面得到的密度和能量信息结合方程(19)，更新粒子的压力；
[0112] 6)然后通过(20)式更新粒子搜索的光滑长度信息。
[0113] 步骤四，输出结果：
[0114] 1)每一个时间步结束，保存中间结果，并输出中间结果；
[0115] 2)结束时间循环，并输出最终结果。
[0116] 进一步地，在上述方案中，所述初始化变量信息和运行参数具体设置如下：
[0117] CASE 1：天然气管道左边冲击右边扩张问题
[0118] 本实验模拟问题的物理模型见图2，计算域是长度为30m的一维空间，管道的两侧均是打开状态，左端设置为入口，右端设置的是出口，管道的中间位置设置了一层隔板来模拟天然气管道中的
阀门。在初始时刻，隔板两侧存在两种不同状态的气体，通过打开中间隔板来模拟实验问题。本实验中我们在出入口都添加了虚粒子，用来弥补边界粒子处的缺失问题。
[0119] 在CASE1问题中，我们初始条件设置为：PL＝7,VL＝0,ρL＝1,PR＝10,ρR＝1,[0120] VR＝0。在本实验中，我们在管道左端布置1000个质量相同粒子，管道的右端布置1000个质量相同的粒子。实验模拟的时长为1s，时间积分步长为10e-5s。Monaghon型人工粘度的参数α＝1，β＝1；人工热量项的系数为g1＝0.2，g2＝0.4。
[0121] 基于管道方程的SPH离散格式(15)～(17)，模拟一维天然气管道的拉格朗日粒子方法为：
[0122] 步骤一，初始化。初始化系统的相关变量和运行参数，具体包括：
[0123] 如图1所示，天然气管道长度L为30m，管道左端0～15m处的压力、速度、密度分别为：7Pa、0m/s、1kg/m3,右端15～30m处依次分别为：10Pa、0m/s、1kg/m3。粒子间距Δx0＝0.01m，计算时间步长为10e-5s，计算时间为1s。在本实验中，采用可变光滑长度进行计算，初始光滑长度h0＝1.5Δx。取三次样条函数作为核函数。
[0124] 步骤二，生成粒子信息，具体包括：
[0125] 初始化粒子步骤中，总共生成气体粒子2000个(不包括虚粒子)，在计算域的左右两侧初始边界布置虚粒子，保证距左右两测初始处2h范围内的气体粒子能够进行正常运算，虚粒子的密度、速度、压力和能量与气体粒子相同。
[0126] 步骤三，列出求解方程并迭代计算：
[0127]
[0128]
[0129]
[0130] 其中P是气体压力,ν是气体的速度，ρ是气体的密度，e是气体的能量。(1)～(3)中的拉格朗日控制方程是未闭合的，需要添加状态方程来保证系统的完整性。对于理想气体，其状态方程为：
[0131] p＝(γ-1)ρe (4)
[0132] 对于实际管道中的气锤问题，如果温度变化不大，可以将问题简化，相应的状态方程为:
[0133]
[0134] 此时，能量项将从控制方程中解耦出来。
[0135] 光滑粒子流体动力学方法中，函数f(x)的积分表示式为：
[0136] f(x)＝∫Ωf(x′)δ(x-x′)dx′ (6)
[0137] 其中δ(x-x′)为狄拉克函数，Ω为包含x的积分体积。若用光滑函数W(x-x′,h)取代狄拉克函数，则f(x)的积分表示式为：
[0138] f(x)≈∫Ωf(x′)W(x-x′,h)dx′ (7)
[0139] 函数导数积分表示为：
[0140]
[0141] 又因为
[0142] 所以
[0143]
[0144]
[0145] 然后进行粒子近似得到：
[0146]
[0147] 因为
[0148]
[0149] 所以
[0150]
[0151] 又有
[0152]
[0153] 所以，在粒子i处的函数的粒子近似式可写为：
[0154]
[0155] 利用光滑粒子流体动力学方法将方程(1)～(3)离散可得(16)～(18)式[0156]
[0157]
[0158]
[0159] p＝(γ-1)ρe (19)
[0160]
[0161] 其中Πij为Monaghon型人工粘度，Hij为人工热量项，γ、ξ为常系数。
[0162] 根据方程(16)、(17)和(18)计算不同时刻每个粒子的密度、速度和能量信息，进而根据(19)式求解出粒子的压强。具体计算过程如下所示：
[0163] 具体计算过程为：
[0164] 1)循环每个时间步；
[0165] 2)初始化计算域内所有粒子属性信息后，搜索目标粒子的邻近粒子，获得初始值的速度和质量，然后通过方程(16)计算该粒子的密度的导数，最后通过欧拉时间积分更新粒子的密度信息；
[0166] 3)通过上一步得到的更新的密度信息、初始设置的压力值和Monaghon型人工粘度方程计算得到速度导数，然后通过欧拉时间积分更新粒子的速度信息，最后通过速度信息更新粒子的位置；
[0167] 4)根据上面两步得到的密度、速度粒子信息和Monaghon型人工粘度方程及人工热量项方程计算得到能量导数，然后通过欧拉时间积分更新粒子的能量信息；
[0168] 5)通过上面得到的密度和能量信息结合方程(19)，更新粒子的压力；
[0169] 6)然后通过(20)式更新粒子搜索的光滑长度信息。
[0170] 步骤四，输出结果：
[0171] 1)每一个时间步结束，保存中间结果，并输出中间结果；
[0172] 2)结束时间循环，并输出最终结果。
[0173] CASE 2：管道气锤问题
[0174] 本实验模拟问题的物理模型见图4，天然气管道长度L为37.5m，气罐内压力PR为250Kpa，管道内气体的初始速度为10m/s，管道直径D为0.0221m，气体的密度ρ为1.2kg/m3，
摩擦系数f为0.02。天然气传输的初始条件为V(x,0)＝V0和 0
在数值模拟过程中，初始均匀分布了376个粒子(包括上下游的虚粒子)，计算时间步长为
0.0001s，计算总时间为46.524s。
[0175] 在一个实例中，具体步骤如下：
[0176] 步骤一，初始化：初始化相关的变量和粒子(包括添加的虚粒子)信息。具体包括：
[0177] 1)初始化与问题相关的变量信息：初始化变量信息管道直径D为0.0221m，天然气3
管道长度L为37.5m，气罐内压力PR为250Kpa，气体的密度ρ为1.2kg/m ，管道内气体的初始速度为10m/s，计算时间步长为0.0001s，计算总时长为46.524s等；
[0178] 2)初始化气体粒子信息，在流体域均匀分布粒子，并添加初始的信息：初始化流体粒子，管道均匀分布共372个气体粒子，气体粒子信息为，V(x,0)＝V0和0
[0179] 3)初始化虚粒子信息，在流体上下游边界分别布上两层虚粒子并根据边界条件添加初始信息：上下游边界各两个虚粒子，上游虚粒子压力为PR，下游虚粒子压力为0，上下游虚粒子初始速度均为V0。
[0180] 步骤二，列出求解方程并迭代计算,步骤同CASE1中的步骤三，将状态方程替换成方程(5)。
[0181] 步骤三同CASE1中步骤四。
[0182] 尽管上面结合图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多
变形，这些均属于本发明的保护之内。
