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

摘要

一种弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法，包括：步骤S1、一计算模块对于给定的低雷诺数不可压缩流，判断式（1）是否成立，

说明书

技术领域

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

背景技术

应用光滑粒子流体动力学(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）是否成立，

$ϵ=\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、若步骤S1中判断公式（1）成立，计算模块在模拟流场中生成 或继承粒子（这是因为在刚开始模拟时是生成粒子，而后面的模拟中 是继承粒子），并采用多边界切线技术处理固壁边界，根据该技术的 需要随之生成虚拟粒子；

步骤S3、计算模块根据公式（2）计算每个粒子在任意位置x点处的体 积力F_A2，并且入口处的体积力F_A1 由公式（3）确定；

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

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

其中，dA表示面积元，ρ为流体密度，Δp为压力差；

步骤S4、输出模块输出步骤S3中任意位置x点处的体积力F_A2。

优选的，在所述步骤S2和S3之间还包括，

步骤S21、计算模块搜索相邻的粒子，并计算光滑函数；

步骤S22、计算模块确定是否用求和方法计算密度值，如果是则执行步 骤S23，否则执行步骤S24；

步骤S23、由式（6）或式（7）得出密度值；

${\rho }_a=\underset{b}{\Sigma }{m}_{b}W\left(\right|r_a-r_b|,h)---\left(6\right)$

${\rho }_a=\frac{\underset{b}{\Sigma }{m}_b{W}_{\mathrm{ab}}}{\underset{b}{\Sigma }\left(\frac{m_b}{{\rho }_b}\right)W_{\mathrm{ab}}}---\left(7\right)$

其中，m_b为粒子b的质量，W_ab为粒子b对粒子a产生影响的光滑函数。

步骤S24、由式（8）计算出密度变化率；

$\frac{d{\rho }_a}{\mathrm{dt}}=\underset{b}{\Sigma }{m}_b(u_a-u_b)·{▿}_a{W}_{\mathrm{ab}}---\left(8\right).$

优选的，在步骤S3和S4之间还包括，

步骤S31、由式（5）、式（9）和式（10）分别计算出动态压力、动态 压力梯度和粘度；

p=c²ρ              (5)

${(\frac{1}{\rho }▿p)}_a=\underset{b}{\Sigma }{m}_b(\frac{p_a}{{\rho }_a^2}+\frac{p_b}{{\rho }_b^2}){▿}_a{W}_{\mathrm{ab}}---\left(9\right)$

${\left(\frac{\mu }{\rho }{▿}^{2}u\right)}_a=\underset{b}{\Sigma }{m}_b\frac{({\mu }_a+{\mu }_b)}{{\rho }_a{\rho }_b}\frac{r_{\mathrm{ab}}·{▿}_a{W}_{\mathrm{ab}}}{{\left|r_{\mathrm{ab}}\right|}^2+{\eta }^2}(u_a-u_b)---\left(10\right)$

其中，c为人工声速，m_b为粒子a的相邻粒子b的质量，μ_a、μ_b分别表 示粒子a，b的动态粘度，，▽_aW_ab为粒子b对粒子a产生影响的光滑函数 的梯度，|r_ab|为粒子a，b之间的距离，η=0.1h，h为光滑长度；

步骤S32、再判断式（1）是否有效，如果有效则执行步骤S33，否则输 出模块报错；

步骤S33、计算流体动量和能量的增加值；

步骤S34、更新位置、速度、能量、密度的数值；

步骤S35、判断雷诺数R_e是否小于1，如果小于1则输出模块输出任意位 置x点处的体积力F_A2，否则输出模块报错。

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

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

附图说明

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

图1是本发明所述弯曲边界上低雷诺数不可压缩流中压力差的仿真方法 流程图。

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

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

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

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

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

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

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

具体实施方式

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

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

参阅图2所示，为了方便对本发明的描述，本发明考虑一个轴对称的、 非平直但变化缓慢的曲线边界的管道内流动的低雷诺数不可压缩流， 并且：假设管道壁为固壁边界，其轴向剖面如图1所示，采用柱坐标系 ，以x和r分别表示轴向和径向坐标；并给定管道的入口出口两端压差 △p= p1－p2，其中，p1为入口处的压力、p2为出口处的压力。

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

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

$ϵ=\frac{\delta }{x_0}<<1,$或R₀<< x₀，L>> R₀，R_e<1            （1）；

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

同时，对流体必须满足的连续性和运动方程进行量级估计，并进行数 学推导后可得：

其中， dA表示面积元，F_A1 、F_A2分别表示A1、A2截面处的体积力。

则根据上式（2），如果F_A1已知，即可容易地求出流场中任意x点处的 体积力F_A2。

然而，在一般的流场的已知条件中很少直接给出F_A1，更多的情形是给 出入口与出口之间的压力差△p。因此，在这些情形中，必须要建立F _A1与压力差△p之间的关系，其应该满足公式（3）。

其中，ρ为流体密度。 

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

下面再结合附图1对本发明的仿真方法流程作进一步介绍：

第一步、计算模块对于上述给定的一个低雷诺数不可压缩流， 判断 式（1）是否成立，

第二步、若第一步中判断式（1）成立，计算模块在流场生成或继承粒 子（这是因为在刚开始模拟时是生成粒子，而后面的模拟中是继承粒 子），并采用多边界切线技术处理固壁边界，根据该技术的需要随之 生成虚拟粒子；

第三步、计算模块搜索相邻的粒子，并计算光滑函数；

第四步、计算模块确定是否用求和方法计算密度值，如果是则执行第 五步，否则执行第六步；

第五步、由式（6）或式（7）得出密度值；

其中，m_b为粒子b的质量，W_ab为粒子b对粒子a产生影响的光滑函数。

第六步、由式（8）计算出密度变化率；

第七步、计算模块根据式（2）计算每个粒子在任意位置x点处的体积 力F_A2，并且入口处的体积力F_A1 由（3）式确定；

第八步、由式（5）、式（9）和式（10）分别计算出动态压力、动态 压力梯度和粘度；

p=c²ρ            (5)

其中，c为人工声速，m_b为粒子a的相邻粒子b的质量，μ_a为，μ_b分别 表示粒子a，b的动态粘度，▽_aW_ab为粒子b对粒子a产生影响的光滑函数 的梯度，|r_ab|为粒子a，b之间的距离，η=0.1h，h为光滑长度。

第九步、再判断式（1）是否有效，如果有效则执行第十步，否则输出 模块报错；

第十步、计算流体动量和能量的增加值；

第十一步、更新位置、速度、能量、密度的数值；

第十二步、判断雷诺数R_e是否小于1，如果小于1则执行第十三步，否 则输出模块报错。

第十三步、输出模块输出步骤S3中任意位置x点处的体积力F_A2。

为了检验本发明的方法，下面通过它去模拟三个低雷诺数不可压缩流 ，包括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^-4 N/m²，通过式（1）、式（2）和式（3）模拟的结果如图3所示 ，并且从图3中可以看出，本发明的方法与用常数体积力模拟的结果相 一致，且都与理论解吻合良好。

实施例二

本实施例采用局部膨胀管道流，如图4所示，可以考虑一段轴对称的血 管，其压差驱动流体在管道内流动，压差为△p= p₁－ p₂= 1.939 006287×10^-3 N/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轴速度分别如图5和图6所示。并且从图5和图6 中可以看出，采用本发明的方法的模拟结果与理论解吻合良好，但采 用常数体积力的模拟结果却不能与理论解吻合。

实施例三

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

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

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