弯曲边界上的低雷诺数不可压缩流中的压力差的处理方法

技术领域

本发明涉及流体力学技术领域，更具体地说，特别涉及一种弯曲边界上的低雷诺数不可 压缩流中的压力差的处理方法。

背景技术

应用光滑粒子流体动力学(Smoothed Particle Hydrodynamics，简称SPH)模拟低雷诺数不 可压缩流时，求解驱动流体运动的压力梯度是很重要的，因为压力在Navier-Stokes方程中只 是表现为梯度。在弱可压缩SPH（Weakly Compressible SPH，简称WCSPH）算法中，总压力通 常被分解为动态压力和静水压力，因此总压力的梯度也就可以通过这两个压力的梯度来获得。

对于WCSPH方法来说，模拟动态压力梯度是简单而又直接的，而静水压力梯度通常被看 作是一个体积力。Morris在1997年用WCSPH研究了低雷诺数不可压缩流，他的测试算例是 Poiseuille流和绕柱流，所得的结果与有限差分法的结果吻合得很好。刘谋斌和他的同事在2005 年用有限粒子法也模拟了Poiseuille流，结果也相当不错。他们都把静水压力梯度（或者静水 压力差）转化为体积力。

对于边界平直的低雷诺数不可压缩流，这个转化是简单的，因为在这些情形中，流场中 的静水压力梯度是一个常数，相应的体积力可以简单地由入口与出口的压力之差除以流场的 长度来得到。然而，对于边界弯曲的低雷诺数不可压缩流来说，静水压力梯度是不均匀的， 各处的静水压力梯度并不是常数，怎样计算各处相应的体积力就成了一个问题。因此，需要 研究一种弯曲边界上低雷诺数不可压缩流中压力差的处理方法。

发明内容

本发明的目的在于针对现有技术存在不能处理弯曲边界上的低雷诺数不可压缩流中 的压力差的技术问题，提供一种弯曲边界上的低雷诺数不可压缩流中的压力差的处理方 法。

为了达到上述目的，本发明采用的技术方案如下：

弯曲边界上的低雷诺数不可压缩流中的压力差的处理方法，在该方法中，低雷诺数不可 压缩流在具有曲线边界的管道内流动，所述的具有曲线边界的管道为轴对称的、非平直的并 且管壁为固壁边界的管道，且该管道的出入口两端压差为△p，其中△p=p₁－p₂，p₁为入口 处的压力、p₂为出口处的压力，并采用以x和r分别表示轴向坐标和径向坐标的柱坐标系，该 方法具体包括以下步骤，

步骤S1，确认是否满足下列公式（1），若满足，则执行步骤S2，

$ϵ=\frac{\delta }{x_0}<<1,\frac{R_0}{x_0}=O$或R₀＜＜x₀，L＞＞R₀，R_e＜1    (1)；

其中，δ为壁厚，x₀为任意一点的轴向坐标，R₀为管道平直处的半径，L为管的长度，R_e为雷诺数；

步骤S2，计算管道的入口截面处的体积力F_A1，

$F_{A1}=\frac{\mathrm{\Delta p}}{\rho {\int }_{-0.5L}^{0.5L}\frac{\underset{A1}{\int }f(r,x)\mathrm{dA}}{\underset{A2}{\int }f(r,x)\mathrm{dA}}\mathrm{dx}}---\left(2\right);$

其中，ρ为流体密度；

步骤S3，将步骤S2中的公式（2）代入下面的公式（3），得出管道的任意一截面处的体 积力F_A2，

$\frac{F_{A2}}{F_{A1}}\approx \frac{\underset{A1}{\int }f(r,x)\mathrm{dA}}{\underset{A2}{\int }f(r,x)\mathrm{dA}}---\left(3\right);$

其中，dA表示面积元。

优选的，所述公式（3）是对管道内的流体满足连续性和运动方程进行量级计算得出。

与现有技术相比，本发明的优点在于：本发明能够方便的将弯曲边界上的低雷诺数不可 压缩流中压力差转化为任意一处的体积力。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的轴对称非平直边界的管段结构示意图。

图2是本发明的实施例一中的Poiseuille流的模拟结果图。

图3是本发明的实施例二中的局部膨胀管流的结构示意图。

图4是本发明的实施例二中的局部膨胀管流在X轴上的模拟结果图。

图5是本发明的实施例二中的局部膨胀管流在Y轴上的模拟结果图。

图6是本发明的实施例三中的倾斜平板流的结构示意图。

图7是本发明的实施例三中的倾斜平板流在X轴上的模拟结果图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步的详细说明。

提出本发明的目的在于：对于边界平直的低雷诺数流，静水压力梯度在整个流场中是一 个常数，相应的体积力可以简单地由入口与出口的压力之差除以流场的长度来得到。然而， 对于边界弯曲的流场来说，各处的静水压力梯度并不是常数，各处的体积力不能简单地由入 口与出口的压力之差除以流场的长度来得到。因此，必须提出一个方法计算各处相应的体积 力。

参阅图1所示，为了方便对本发明的描述，本发明中提供的雷诺数不可压缩流在具有曲 线边界的管道内流动，并且具有曲线边界的管道为轴对称的、非平直的并且管壁为固壁边界 的管道，给定该管道的出入口两端压差为△p，其中△p=p₁－p₂，p₁为入口处的压力、p₂为 出口处的压力，并采用以x和r分别表示轴向坐标和径向坐标的柱坐标系。

在这一流场中沿X轴方向不同位置的静水压力梯度不是处处相等的，而怎样把流场两端 的压差转化为粒子的体积力是对这一流场进行SPH模拟需要解决的问题。

设定管道的直径很小，并水平放置，这样管道内的流体在入口与出口之间的压力梯度驱 动下流动，并且其重力的影响可以忽略。

此时，本发明的处理方法具体为：

首先，确认管道内是否满足下列公式（1），当满足公式（1）时，可进行后续处理：

$ϵ=\frac{\delta }{x_0}<<1,\frac{R_0}{x_0}=O$或R₀＜＜x₀，L＞＞R₀，R_e＜1    (1)；

其中，δ为壁厚，x₀为任意一点的轴向坐标，R₀为管平直处的半径，L为管的长度，R_e为雷诺数；

其次，对管道内流体必须满足的连续性和运动方程进行量级估计，并进行数学推导（该 数学推导过程一种计算推导方法，此处不再赘述）后可得管道的入口截面处的体积力F_A1和管 道的任意一截面处的体积力F_A2满足公式（2）：

$\frac{F_{A2}}{F_{A1}}\approx \frac{\underset{A1}{\int }f(r,x)\mathrm{dA}}{\underset{A2}{\int }f(r,x)\mathrm{dA}}---\left(2\right);$

其中，dA表示面积元；

这样，从式（2）中可以看出，如果F_A1已知，则可以很容易地求出流场中任意x点处的体 积力F_A2；然而，一般的流场的已知条件中很少直接给出F_A1，更多的情形是给出入口与出口之 间的压力差。

最后，建立F_A1与压力差△p之间的关系，即得出管道的入口截面处的体积力F_A1满足式（3）；

$F_{A1}=\frac{\mathrm{\Delta p}}{\rho {\int }_{-0.5L}^{0.5L}\frac{\underset{A1}{\int }f(r,x)\mathrm{dA}}{\underset{A2}{\int }f(r,x)\mathrm{dA}}\mathrm{dx}}---\left(3\right);$

其中，ρ为流体密度。

通过上述的式（2）和式（3），并将式（3）代入式（2）中得出管道的任意一截面处的体 积力F_A2；即可将流场两端的压差△p转化为任意一处粒子的体积力，使得用SPH模拟弯曲边界 的低雷诺数不可压缩流可以顺利进行。

为了检验本发明的方法，下面通过它去模拟三个低雷诺数不可压缩流，包括Poiseuille流、 局部膨胀管流和倾斜平板流。为了比较，这些例子同样用常数体积力的WCSPH方法进行模拟， 并且模拟的流体都是水，具有相同的初始密度ρ0=10³kg/m^3；同时，三个实施例都给出了理论 解以供比较。

实施例一

本实施例采用Poiseuille流，即将两块无限大的平板分别放在坐标y=-0.5YL和y=0.5YL 处，其中，YL是两块平板之间的距离。其初始静止的流体由平行于X轴的体积力F（相应于 入口与出口的压差）驱动，最后将达到一个稳定状态。

在本实施例中，YL=10^-3m，流场长度为d=5×10^-4m，静水压力差为△p=10^-4N/m²，通过 式（1）、式（2）和式（3）模拟的结果如图2所示，并且从图2中可以看出，本发明的方法 与用常数体积力模拟的结果相一致，且都与理论解吻合良好。

实施例二

本实施例采用局部膨胀管流，如图3所示，可以考虑一段轴对称的血管，其压差驱动流 体在管道内流动，压差为△p=p₁－p₂=1.939006287×10^-3N/m²，管道的半径是关于位置x的 函数，可以表示为（4）式：

$\frac{R}{R_0}=\left(\begin{array}{cc}1& \left|x\right|>X_0\\ 1+\frac{\delta }{{2R}_0}(1+\cos \frac{\mathrm{\pi x}}{X_0})& \left|x\right|\le X_0\end{array}\right),---\left(4\right);$

其中，R₀=0.5×10^-3m，δ=0.2R₀，L=6×10^-3m，X₀=2×10^-3m，ε=δ/X₀=0.05。

模拟结果，X轴速度和Y轴速度分别如图4和图5所示。并且从图4和图5中可以看出， 采用本发明的方法的模拟结果与理论解吻合良好，但采用常数体积力的模拟结果却不能与理 论解吻合。

实施例三

本实施例采用倾斜平板流，如图6所示，在本实施例中，d_BC=4mm,2l₁=0.5mm,α＝ 3.503°，且B与C之间的静水压力差为△p=1.21665968×10^-3N/m²。

采用本发明的方法，沿X轴上的流体粒子的水平速度的模拟结果如图7所示；从图7中 可以看出，用本发明的方法的模拟结果再一次与理论解相符，但常数体积力情况下的模拟却 不能相符。

虽然结合附图描述了本发明的实施方式，但是专利所有者可以在所附权利要求的范围之 内做出各种变形或修改，只要不超过本发明的权利要求所描述的保护范围，都应当在本发明 的保护范围之内。