CN110569541A - 管道水锤分析方法 - Google Patents

Abstract

本发明涉及流体力学技术领域，为提出一种分析管道水锤的新方法，求解出移动坐标***下的水锤方程，实现对管道中水锤现象的模拟。同时在满足数值精度的前提下，能够更方便地模拟管道水锤问题。为此，本发明采取的技术方案是，管道水锤分析方法，步骤如下：步骤一、初始化：初始化相关的变量、粒子、上下游虚粒子信息；步骤二、列出求解方程并迭代计算：步骤三，输出结果。本发明主要应用于流体力学工程设计场合。

Description

管道水锤分析方法

技术领域

本发明涉及流体力学技术领域，具体是涉及一种分析管道水锤的新方法。

背景技术

管道水锤问题在不考虑水的弹性时可以用刚性柱理论、特征线法或一些网格方法来进行求解，但是，在大多情况下，为了更精确的求解问题，水体的弹性有时必须考虑，而特征线法或网格方法通常要追踪界面的移动，即而会产生复杂的计算问题，并且计算精度也随着插值误差的累积而降低。用密度求和形式的无网格粒子方法解决管道水锤问题时，既能够充分的考虑水的弹性，而且计算过程也更加的方便，在实际应用中能够更好的模拟管道水锤问题。

发明内容

为克服现有技术的不足，本发明旨在提出一种分析管道水锤的新方法，求解出移动坐标***下的水锤方程，实现对管道中水锤现象的模拟。同时在满足数值精度的前提下，能够更方便地模拟管道水锤问题。为此，本发明采取的技术方案是，管道水锤分析方法，步骤如下：

步骤一、初始化：初始化相关的变量、粒子、上下游虚粒子信息；

步骤二、列出求解方程并迭代计算：

求解水锤方程为：

其中P是压力，V是流体的速度，ρ是流体的密度，c是声波速度，g是重力加速度，θ是管道倾角，D是管道内直径，λ是摩擦系数t和x表示时间和沿着管道的距离分布；

当用密度求和法求解水锤中的压力时，控制方程为：

ρ_i＝∑m_jW_ij (4)

P_i-P₀＝c²(ρ_i-ρ₀) (5)

其中ρ₀是给定参考压力P₀下的流体密度，ρ_i表示i粒子的密度，m_j表示j粒子的质量，P_i表示i粒子处的压力，W指光滑函数；

上游边界条件为

其中P(0，t)表示在任意时刻t管道上位置为0处的压力，P_R表示水库压力，V(0，t)表示在任意时刻t管道上位置为0处的速度；

下游边界条件为

V(L,t)＝0 (8)

上式表示在任意时刻t，管道位置L处的速度为0；

在光滑粒子流体动力学的方法中，任意函数f(x)的积分形式表示为：

其中δ(x-x')为狄拉克函数，Ω为包含x的积分体积，具体表达式如下：

若用光滑函数W(x-x',h)代替狄拉克函数，则f(x)的积分表示式为：

其中h是光滑函数的光滑长度，在SPH(光滑粒子流体动力学)方法的推导中，核近似算子一般用角括弧标记，即函数的积分表示为：

<f(x)>＝∫_Ωf(x')W(x-x',h)dx' (12)

为了得到函数导数的近似，我们用f'(x)取代f(x)，即函数导数近似为:

又因为：

[f'(x')]W(x-x',h)＝[f(x')W(x-x',h)]'-f(x')W'(x-x',h) (14)

所以：

将积分形式转化为求和形式，进行粒子近似得到：

其中N为光滑函数紧支域中的粒子个数，为粒子对应的长度，导数W'(x-x_j,h)与粒子j相关；又因为：

所以：

又因为：

所以，在粒子i处的函数的粒子近似式写为：

利用光滑粒子流体动力学方法对控制方程(3)、(4)、(5)离散得到：

ρ_i＝∑m_jW_ij (21)

P_i-P₀＝c²(ρ_i-ρ₀) (23)

其中，V_i表示i粒子处的速度，W_ij,x表示光滑函数对x的一阶导数，Π_ij为人工黏性项，Monaghan型黏性项具体表达式为：

其中，

上述公式中，c和V分别代表人工声速和粒子的速度矢量，α和β为常数，Monaghan型黏性项是一种在光滑粒子流体动力学中最常用的黏性项；

根据控制方程，计算不同时刻每个粒子的压力和速度信息；

步骤三，输出结果：

1)结束粒子循环，保存并输出中间结果；

2)结束时间循环，并输出最终结果。

利用黎曼Riemann解构造的新的人工粘性来取代人工粘性项，具体实现如下：

将方程(22)中的压力项P_i和P_j用Riemann解替代，得到用Riemann解构造的新的人工粘性来取代传统SPH中的人工粘性项的控制方程：

其中，

根据离散形式的控制方程，计算不同时刻每个粒子的压力和速度信息：

根据控制方程，计算不同时刻每个粒子的压力和速度信息，具体计算过程为：

1)对时间变量循环；

2)对粒子循环；

3)对初始化后的粒子根据离散形式的控制方程得到流体粒子的压力信息，并更新流体粒子的压力信息；

4)根据方程(6)(7)计算并更新上游虚粒子的压力信息；

5)根据离散形式的控制方程计算并更新流体粒子的速度信息；

6)根据上下游边界流体粒子速度分别更新上下游虚粒子的速度；

7)根据粒子的速度更新粒子位置、粒子的压力、上下游虚粒子速度和压力信息。

本发明的特点及有益效果是：

本发明提供了一种分析管道水锤的新方法，该方法通过密度求和法求解水锤压力，用SPH方法求解水锤速度，进而求解出移动坐标***下的水锤方程。同时针对SPH方法中的人工粘性项本发明利用Riemann解构造的新的人工粘性对水锤方程进行求解，新的粘性无需再人工调节粘性的系数。能够充分考虑水的弱可压缩性带来的影响，在满足数值精度的前提下，能够更方便地模拟管道水锤问题。

附图说明：

图1是管道水锤的物理模型图；

图2是求解管道水锤问题的计算流程图；

图3是用求解管道水锤新方法得出的压力结果图。

图4是用新的人工黏性项求解水锤的压力结果图。

具体实施方式

本发明解决的技术问题是提供一种分析管道水锤的新方法，该方法通过密度求和法求解水锤压力，用SPH(光滑粒子流体动力学)方法求解水锤速度，进而求解出移动坐标***下的水锤方程。实现对管道中水锤现象的模拟。同时针对SPH方法中的人工粘性项本发明利用Riemann解构造的新的人工粘性对水锤方程进行求解。能够充分考虑水的弱可压缩性带来的影响，在满足数值精度的前提下，能够更方便地模拟管道水锤问题。为了解决上述技术问题，本发明的技术方案是：

一种分析管道水锤的新方法，包括以下步骤：

步骤一、初始化:初始化相关的变量、粒子、上下游虚粒子信息；

步骤二、列出求解方程并迭代计算：

求解水锤方程为：

其中P是压力，V是流体的速度，ρ是流体的密度，c是声波速度，g是重力加速度，θ是管道倾角，D是管道内直径，λ是摩擦系数t和x表示时间和沿着管道的距离分布；

当用密度求和法求解水锤中的压力时，控制方程为：

ρ_i＝∑m_jW_ij (4)

P_i-P₀＝c²(ρ_i-ρ₀) (5)

其中ρ₀是给定参考压力P₀下的流体密度，ρ_i表示i粒子的密度，m_j表示j粒子的质量，P_i表示i粒子处的压力，W指光滑函数；

上游边界条件为

其中P(0，t)表示在任意时刻t管道上位置为0处的压力，P_R表示水库的压力，V(0，t)表示在任意时刻t管道上位置为0处的速度；

下游边界条件为

V(L,t)＝0 (8)

上式表示在任意时刻t，管道位置L处的速度为0；

在光滑粒子流体动力学的方法中，任意函数f(x)的积分形式表示为：

其中δ(x-x')为狄拉克函数，Ω为包含x的积分体积，具体表达式如下：

若用光滑函数W(x-x',h)代替狄拉克函数，则f(x)的积分表示式为：

其中h是光滑函数的光滑长度，在SPH(光滑粒子流体动力学)方法的推导中，核近似算子一般用角括弧标记，即函数的积分表示为：

<f(x)>＝∫_Ωf(x')W(x-x',h)dx' (12)

为了得到函数导数的近似，我们用f'(x)取代f(x)，即函数导数近似为:

又因为：

[f'(x')]W(x-x',h)＝[f(x')W(x-x',h)]'-f(x')W'(x-x',h) (14)

所以：

将积分形式转化为求和形式，进行粒子近似得到：

其中N为光滑函数紧支域中的粒子个数，为粒子对应的长度(二维情况下为粒子对应的面积，三维情况下为粒子对应的体积)，导数W'(x-x_j,h)与粒子j相关。

又因为：

所以：

又因为：

所以，在粒子i处的函数的粒子近似式可写为：

利用光滑粒子流体动力学方法对控制方程(3)、(4)、(5)离散可得：

ρ_i＝∑m_jW_ij (21)

P_i-P₀＝c²(ρ_i-ρ₀) (23)

其中，V_i表示i粒子处的速度，W_ij,x表示光滑函数对x的一阶导数，Π_ij为人工黏性项，Monaghan型黏性项具体表达式为：

其中，

上述公式中，c和V分别代表人工声速和粒子的速度矢量，α和β为常数，通常取1.0左右，但是在不同流体状态下取值应适当调整，Monaghan型黏性项是一种在光滑粒子流体动力学中最常用的黏性项，它可以将动能转换为热能，提供冲击波产生的耗散，还可以解决粒子相互接近时的非物理穿透。

在用以上人工粘性求解问题时需人工调整人工粘性的系数，给实验带来诸多不便。针对此问题本发明利用Riemann解构造的新的人工粘性来取代传统SPH中的人工粘性项。具体实现如下：

将方程(22)中的压力项P_i和P_j用Riemann解替代，得到用Riemann解构造的新的人工粘性来取代传统SPH中的人工粘性项的控制方程：

其中，

根据离散形式的控制方程，计算不同时刻每个粒子的压力和速度信息：

具体计算过程为：

1)对时间变量循环；

2)对粒子循环；

3)对初始化后的粒子根据离散形式的控制方程得到流体粒子的压力信息，并更新流体粒子的压力信息；

4)根据方程(6)(7)计算并更新上游虚粒子的压力信息；

5)根据离散形式的控制方程计算并更新流体粒子的速度信息；

6)根据上下游边界流体粒子速度分别更新上下游虚粒子的速度；

7)根据粒子的速度更新粒子位置、粒子的压力、上下游虚粒子速度和压力信息；

步骤三，输出结果：

1)结束粒子循环，保存并输出中间结果；

2)结束时间循环，并输出最终结果。

进一步地，在上述方案中，步骤一所述初始化相关的变量、粒子、上下游虚粒子信息。具体包括：

1)初始化与问题相关的变量信息；

2)初始化流体粒子信息，在流体域均匀分布粒子，并添加初始的信息；

3)初始化虚粒子信息，在流体上下游边界分别布上两层虚粒子并根据边界条件添加初始信息。

进一步地，在以上方案中，初始化变量信息、粒子信息、上下游虚粒子信息具体参数设置如下：管道直径D为0.797m，管道长度L为20m，水库初始时的流量为0.5m3/s,水库压力P_R为1MPa，密度ρ为1000kg/m3，重力加速度g为9.8m/s2，声波波速c为1250.7m/s，初始速度为1.002m/s,粒子总质量为100。水柱的初始条件为V(x,0)＝0和P(x,0)＝P_R；在数值模拟过程中，初始均匀分布了205个粒子，这其中包括上下游两端的虚粒子。计算时间步长为0.0001s，计算总时长为5s。

进一步地，在以上方案中，初始化流体粒子信息具体参数设置如下：管道均匀分布共201个流体粒子，流体粒子信息为，V(x,0)＝0和P(x,0)＝P_R。

进一步地，在以上方案中，所述初始化虚粒子信息参数具体设置如下：上下游边界各两个虚粒子，上游虚粒子压力为P_R，下游虚粒子压力等于最靠近下游端的的流体粒子的压力，上游虚粒子初始速度均为1.002m/s,下游虚粒子采用镜像粒子法，其初始速度为-1.002m/s。

下面结合附图对本发明做进一步详细地描述。

如图1所示，本发明的试验用管道直径D为0.797m，管道长度L为20m，水库初始时的流量为0.5m3/s,水库压力P_R为1MPa，密度ρ为1000kg/m3，重力加速度g为9.8m/s2，声波波速c为1250.7m/s，初始速度为1.002m/s,粒子总质量为100。水柱的初始条件为V(x,0)＝0和P(x,0)＝P_R；在数值模拟过程中，初始均匀分布了205个粒子，这其中包括上下游两端的虚粒子。计算时间步长为0.0001s，计算总时长为5s。

分析管道水锤的新方法为：

步骤一，初始化:初始化相关的变量、粒子、上下游虚粒子信息；具体包括：

1)初始化与问题相关的变量信息：管道直径D为0.797m，管道长度L为20m，水库初始时的流量为0.5m3/s,水库压力P_R为1MPa，密度ρ为1000kg/m3，重力加速度g为9.8m/s2，声波波速c为1250.7m/s，初始速度为1.002m/s,粒子总质量为100。水柱的初始条件为V(x,0)＝0和P(x,0)＝P_R；在数值模拟过程中，初始均匀分布了205个粒子(包括每端的虚粒子)，计算时间步长为0.0001s，计算总时长为5s。

2)初始化流体粒子信息，在流体域均匀分布粒子，并添加初始的信息：初始化流体粒子，管道均匀分布共201个流体粒子，流体粒子信息为，V(x,0)＝0和P(x,0)＝P_R。

3)初始化虚粒子信息，在流体上下游边界分别布上两层虚粒子并根据边界条件添加初始信息：上下游边界各两个虚粒子，上游虚粒子压力为P_R，下游虚粒子压力等于最靠近下游端的的流体粒子的压力，上游虚粒子初始速度均为1.002m/s,下游虚粒子采用镜像粒子法，其初始速度为-1.002m/s。

步骤二，列出求解方程并迭代计算：

求解水锤方程为：

其中P是压力，V是流体的速度，ρ是流体的密度，c是声波速度，g是重力加速度，θ是管道倾角，D是管道内直径，λ是摩擦系数t和x表示时间和沿着管道的距离分布；

当用密度求和法求解水锤中的压力时，控制方程为：

ρ_i＝∑m_jW_ij (4)

P_i-P₀＝c²(ρ_i-ρ₀) (5)

其中ρ₀是给定参考压力P₀下的流体密度，ρ_i表示i粒子的密度，m_j表示j粒子的质量，P_i表示i粒子处的压力，W指光滑函数；

上游边界条件为

其中P(0，t)表示在任意时刻t管道上位置为0处的压力，P_R表示水库的压力，V(0，t)表示在任意时刻t管道上位置为0处的速度；

下游边界条件为

V(L,t)＝0 (8)

上式表示在任意时刻t，管道位置L处的速度为0；

在光滑粒子流体动力学的方法中，任意函数f(x)的积分形式表示为：

其中δ(x-x')为狄拉克函数，Ω为包含x的积分体积，具体表达式如下：

若用光滑函数W(x-x',h)代替狄拉克函数，则f(x)的积分表示式为：

其中h是光滑函数的光滑长度，在SPH(光滑粒子流体动力学)方法的推导中，核近似算子一般用角括弧标记，即函数的积分表示为：

<f(x)>＝∫_Ωf(x')W(x-x',h)dx' (12)

为了得到函数导数的近似，我们用f'(x)取代f(x)，即函数导数近似为:

又因为：

[f'(x')]W(x-x',h)＝[f(x')W(x-x',h)]'-f(x')W'(x-x',h) (14)

所以：

将积分形式转化为求和形式，进行粒子近似得到：

其中N为光滑函数紧支域中的粒子个数，为粒子对应的长度(二维情况下为粒子对应的面积，三维情况下为粒子对应的体积)，导数W'(x-x_j,h)与粒子j相关。

又因为：

所以：

又因为：

所以，在粒子i处的函数的粒子近似式可写为：

利用光滑粒子流体动力学方法对控制方程(3)、(4)、(5)离散可得：

ρ_i＝∑m_jW_ij (21)

P_i-P₀＝c²(ρ_i-ρ₀) (23)

其中，V_i表示i粒子处的速度，W_ij,x表示光滑函数对x的一阶导数，Π_ij为人工黏性项，Monaghan型黏性项具体表达式为：

其中，

上述公式中，c和V分别代表人工声速和粒子的速度矢量，α和β为常数，通常取1.0左右，但是在不同流体状态下取值应适当调整，Monaghan型黏性项是一种在光滑粒子流体动力学中最常用的黏性项，它可以将动能转换为热能，提供冲击波产生的耗散，还可以解决粒子相互接近时的非物理穿透。

在用以上人工粘性求解问题时需人工调整人工粘性的系数，给实验带来诸多不便。针对此问题本发明利用Riemann解构造的新的人工粘性来取代传统SPH中的人工粘性项。具体实现如下：

将方程(22)中的压力项P_i和P_j用Riemann解替代，得到用Riemann解构造的新的人工粘性来取代传统SPH中的人工粘性项的控制方程：

其中，

根据离散形式的控制方程，计算不同时刻每个粒子的压力和速度信息：

具体计算过程为：

3)对时间变量循环；

4)对粒子循环；

3)对初始化后的粒子根据离散形式的控制方程得到流体粒子的压力信息，并更新流体粒子的压力信息；

4)根据方程(6)(7)计算并更新上游虚粒子的压力信息；

5)根据离散形式的控制方程计算并更新流体粒子的速度信息；

6)根据上下游边界流体粒子速度分别更新上下游虚粒子的速度；

7)根据粒子的速度更新粒子位置、粒子的压力、上下游虚粒子速度和压力信息；

步骤三，输出结果：

1)结束粒子循环，保存并输出中间结果；

2)结束时间循环，并输出最终结果。

尽管结合图对本发明进行了描述，但本发明并不仅仅局限于上述的具体实施方式，上述的所述的具体实施方式只是一种而已(是示意性的)，它不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。

Claims (2)

1.一种管道水锤分析方法，其特征是，步骤如下：

步骤一、初始化：初始化相关的变量、粒子、上下游虚粒子信息；

步骤二、列出求解方程并迭代计算：

求解水锤方程为：

其中P是压力，V是流体的速度，ρ是流体的密度，c是声波速度，g是重力加速度，θ是管道倾角，D是管道内直径，λ是摩擦系数t和x表示时间和沿着管道的距离分布；

当用密度求和法求解水锤中的压力时，控制方程为：

ρ_i＝∑m_jW_ij (4)

P_i-P₀＝c²(ρ_i-ρ₀) (5)

其中ρ₀是给定参考压力P₀下的流体密度，ρ_i表示i粒子的密度，m_j表示j粒子的质量，P_i表示i粒子处的压力，W指光滑函数；

上游边界条件为

其中P(0，t)表示在任意时刻t管道上位置为0处的压力，P_R表示水库压力，V(0，t)表示在任意时刻t管道上位置为0处的速度；

下游边界条件为

V(L,t)＝0 (8)

上式表示在任意时刻t，管道位置L处的速度为0；

在光滑粒子流体动力学的方法中，任意函数f(x)的积分形式表示为：

其中δ(x-x')为狄拉克函数，Ω为包含x的积分体积，具体表达式如下：

若用光滑函数W(x-x',h)代替狄拉克函数，则f(x)的积分表示式为：

其中h是光滑函数的光滑长度，在光滑粒子流体动力学SPH方法的推导中，核近似算子一般用角括弧标记，即函数的积分表示为：

<f(x)>＝∫_Ωf(x')W(x-x',h)dx' (12)

为了得到函数导数的近似，我们用f'(x)取代f(x)，即函数导数近似为:

又因为：

[f'(x')]W(x-x',h)＝[f(x')W(x-x',h)]'-f(x')W'(x-x',h) (14)

所以：

将积分形式转化为求和形式，进行粒子近似得到：

其中N为光滑函数紧支域中的粒子个数，为粒子对应的长度，导数W'(x-x_j,h)与粒子j相关；又因为：

所以：

又因为：

所以，在粒子i处的函数的粒子近似式写为：

利用光滑粒子流体动力学方法对控制方程(3)、(4)、(5)离散得到：

ρ_i＝∑m_jW_ij (21)

P_i-P₀＝c²(ρ_i-ρ₀) (23)

其中，V_i表示i粒子处的速度，W_ij,x表示光滑函数对x的一阶导数，Π_ij为人工黏性项，Monaghan型黏性项具体表达式为：

其中，

上述公式中，c和V分别代表人工声速和粒子的速度矢量，α和β为常数，Monaghan型黏性项是一种在光滑粒子流体动力学中最常用的黏性项；

根据控制方程，计算不同时刻每个粒子的压力和速度信息；

步骤三，输出结果：

1)结束粒子循环，保存并输出中间结果；

2)结束时间循环，并输出最终结果。

2.如权利要求1所述的管道水锤分析方法，其特征是，利用黎曼Riemann解构造的新的人工粘性来取代人工粘性项，具体实现如下：

将方程(22)中的压力项P_i和P_j用Riemann解替代，得到用Riemann解构造的新的人工粘性来取代传统SPH中的人工粘性项的控制方程：

其中，

根据离散形式的控制方程，计算不同时刻每个粒子的压力和速度信息：

根据控制方程，计算不同时刻每个粒子的压力和速度信息，具体计算过程为：

1)对时间变量循环；

2)对粒子循环；

3)对初始化后的粒子根据离散形式的控制方程得到流体粒子的压力信息，并更新流体粒子的压力信息；

4)根据方程(6)(7)计算并更新上游虚粒子的压力信息；

5)根据离散形式的控制方程计算并更新流体粒子的速度信息；

6)根据上下游边界流体粒子速度分别更新上下游虚粒子的速度；

7)根据粒子的速度更新粒子位置、粒子的压力、上下游虚粒子速度和压力信息。