CN103473418B - 弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法 - Google Patents
弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法 Download PDFInfo
- Publication number
- CN103473418B CN103473418B CN201310428248.3A CN201310428248A CN103473418B CN 103473418 B CN103473418 B CN 103473418B CN 201310428248 A CN201310428248 A CN 201310428248A CN 103473418 B CN103473418 B CN 103473418B
- Authority
- CN
- China
- Prior art keywords
- particle
- formula
- reynolds number
- rho
- flow
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
一种弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法,包括:步骤S1、一计算模块对于给定的低雷诺数不可压缩流,判断式(1)是否成立,或R0<< x0,L>> R0,Re<1(1);步骤S2、若S1中的式(1)成立,计算模块在模拟流场中生成或继承粒子,采用多边界切线技术处理固壁边界并根据该技术的需要生成虚拟粒子;步骤S3、计算模块根据式(2)计算每个粒子在任意位置x点处的体积力FA2,且入口处的体积力FA1 由式(3)确定; 步骤S4、一输出模块输出S3中的x点处的体积力FA2。本发明可方便的将弯曲边界上的低雷诺数不可压缩流中压力差转化为任意一处的体积力。
Description
技术领域
本发明涉及计算机仿真技术领域,更具体地说,特别涉及一种弯曲边界上低雷诺数不可压缩流中压力差的仿真方法。
背景技术
应用光滑粒子流体动力学(Smoothed Particle Hydrodynamics ,简称SPH)模拟低雷诺数不可压缩流时,求解驱动流体运动的压力梯度是很重要的,因为压力在Navier-Stokes方程中只是表现为梯度。在弱可压缩SPH(Weakly Compressible SPH,简称WCSPH)算法中,总压力通常被分解为动态压力和静水压力,因此总压力的梯度也就可以通过这两个压力的梯度来获得。
对于WCSPH方法来说,模拟动态压力梯度是简单而又直接的,而静水压力梯度通常被看作是一个体积力。Morris在1997年用WCSPH研究了低雷诺数不可压缩流,他的测试算例是Poiseuille流和绕柱流,所得的结果与有限差分法的结果吻合得很好。刘谋斌和他的同事在2005年用有限粒子法也模拟了Poiseuille流,结果也相当不错。他们都把静水压力梯度(或者静水压力差)转化为体积力。
对于边界平直的低雷诺数不可压缩流,这个转化是简单的,因为在这些情形中,流场中的静水压力梯度是一个常数,相应的体积力可以简单地由入口与出口的压力之差除以流场的长度来得到。然而,对于边界弯曲的低雷诺数不可压缩流来说,静水压力梯度是不均匀的,各处的静水压力梯度并不是常数,怎样计算各处相应的体积力就成了一个问题。因此,需要研究一种弯曲边界上低雷诺数不可压缩流中压力差的仿真方法。
发明内容
本发明的目的在于针对现有技术存在不能处理弯曲边界上的低雷诺数不可压缩流中的压力差的技术问题,提供一种弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法。
为了达到上述目的,本发明采用的技术方案如下:
弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法,包括计算模块和输出模块,在该方法中,计算模块给定一低雷诺数不可压缩流,该低雷诺数不可压缩流在具有曲线边界的管道内流动,所述的具有曲线边界的管道为轴对称的、非平直的并且管道壁为固壁边界的管道,且该管道的出入口两端压差为△p,其中△p = p1-p2,p1为入口处的压力、p2为出口处的压力,并采用以x和r分别表示轴向坐标和径向坐标的柱坐标系,,该方法具体包括以下步骤,
步骤S1、计算模块对于上述给定的低雷诺数不可压缩流,判断公式(1)是否成立,
其中,δ为壁厚,x0为任意一点的轴向坐标,R0为管道平直处的半径,L为管道的长度,Re为雷诺数;
步骤S2、若步骤S1中判断公式(1)成立,计算模块在模拟流场中生成或继承粒子(这是因为在刚开始模拟时是生成粒子,而后面的模拟中是继承粒子),并采用多边界切线技术处理固壁边界,根据该技术的需要随之生成虚拟粒子;
步骤S3、计算模块根据公式(2)计算每个粒子在任意位置x点处的体积力FA2,并且入口处的体积力FA1 由公式(3)确定;
其中,dA表示面积元,ρ为流体密度,Δp为压力差;
步骤S4、输出模块输出步骤S3中任意位置x点处的体积力FA2。
优选的,在所述步骤S2和S3之间还包括,
步骤S21、计算模块搜索相邻的粒子,并计算光滑函数;
步骤S22、计算模块确定是否用求和方法计算密度值,如果是则执行步骤S23,否则执行步骤S24;
步骤S23、由式(6)或式(7)得出密度值;
其中,mb为粒子b的质量,Wab为粒子b对粒子a产生影响的光滑函数。
步骤S24、由式(8)计算出密度变化率;
优选的,在步骤S3和S4之间还包括,
步骤S31、由式(5)、式(9)和式(10)分别计算出动态压力、动态压力梯度和粘度;
p=c2ρ (5)
其中,c为人工声速,mb为粒子a的相邻粒子b的质量,μa、μb分别表示粒子a,b的动态粘度,,▽aWab为粒子b对粒子a产生影响的光滑函数的梯度,|rab|为粒子a,b之间的距离,η=0.1h,h为光滑长度;
步骤S32、再判断式(1)是否有效,如果有效则执行步骤S33,否则输出模块报错;
步骤S33、计算流体动量和能量的增加值;
步骤S34、更新位置、速度、能量、密度的数值;
步骤S35、判断雷诺数Re是否小于1,如果小于1则输出模块输出任意位置x点处的体积力FA2,否则输出模块报错。
优选的,所述公式(2)是对管道内的流体满足连续性和运动方程进行量级计算得出。
与现有技术相比,本发明的优点在于:本发明可以方便的将弯曲边界上的低雷诺数不可压缩流中压力差转化为任意一处的体积力。
附图说明
下面结合附图和实施例对本发明作进一步说明。
图1是本发明所述弯曲边界上低雷诺数不可压缩流中压力差的仿真方法流程图。
图2是本发明的轴对称非平直边界的管段结构示意图。
图3是本发明的实施例一中的Poiseuille流的模拟结果图。
图4是本发明的实施例二中的局部膨胀管流的结构示意图。
图5是本发明的实施例二中的局部膨胀管流在X轴上的模拟结果图。
图6是本发明的实施例二中的局部膨胀管流在Y轴上的模拟结果图。
图7是本发明的实施例三中的倾斜平板流的结构示意图。
图8是本发明的实施例三中的倾斜平板流在X轴上的模拟结果图。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步的详细说明。
提出本发明的目的在于:对于边界平直的低雷诺数流,静水压力梯度在整个流场中是一个常数,相应的体积力可以简单地由入口与出口的压力之差除以流场的长度来得到。然而,对于边界弯曲的流场来说,各处的静水压力梯度并不是常数,各处的体积力不能简单地由入口与出口的压力之差除以流场的长度来得到。因此,必须提出一个方法计算各处相应的体积力。
参阅图2所示,为了方便对本发明的描述,本发明考虑一个轴对称的、非平直但变化缓慢的曲线边界的管道内流动的低雷诺数不可压缩流,并且:假设管道壁为固壁边界,其轴向剖面如图1所示,采用柱坐标系,以x和r分别表示轴向和径向坐标;并给定管道的入口出口两端压差△p= p1-p2,其中,p1为入口处的压力、p2为出口处的压力。
在这一流场中沿X轴方向不同位置的静水压力梯度不是处处相等的,而怎样把流场两端的压差转化为粒子的体积力是对这一流场进行SPH模拟需要解决的问题。
设定管道的直径很小,并水平放置,这样管道内的流体在入口与出口之间的压力梯度驱动下流动,并且其重力的影响可以忽略。此时应当满足,
其中,δ为壁厚,x0为任意一点的轴向坐标,R0为管道平直处的半径,L为管道的长度,Re为雷诺数;
同时,对流体必须满足的连续性和运动方程进行量级估计,并进行数学推导后可得:
其中, dA表示面积元,FA1 、FA2分别表示A1、A2截面处的体积力。
则根据上式(2),如果FA1已知,即可容易地求出流场中任意x点处的体积力FA2。
然而,在一般的流场的已知条件中很少直接给出FA1,更多的情形是给出入口与出口之间的压力差△p。因此,在这些情形中,必须要建立FA1与压力差△p之间的关系,其应该满足公式(3)。
其中,ρ为流体密度。
通过上述的式(2)和式(3),并将式(3)代入式(2)中得出管道的任意一截面处的体积力FA2;即可将流场两端的压差△p转化为任意一处粒子的体积力,使得用SPH模拟弯曲边界的低雷诺数不可压缩流可以顺利进行。
下面再结合附图1对本发明的仿真方法流程作进一步介绍:
第一步、计算模块对于上述给定的一个低雷诺数不可压缩流, 判断式(1)是否成立,
第二步、若第一步中判断式(1)成立,计算模块在流场生成或继承粒子(这是因为在刚开始模拟时是生成粒子,而后面的模拟中是继承粒子),并采用多边界切线技术处理固壁边界,根据该技术的需要随之生成虚拟粒子;
第三步、计算模块搜索相邻的粒子,并计算光滑函数;
第四步、计算模块确定是否用求和方法计算密度值,如果是则执行第五步,否则执行第六步;
第五步、由式(6)或式(7)得出密度值;
其中,mb为粒子b的质量,Wab为粒子b对粒子a产生影响的光滑函数。
第六步、由式(8)计算出密度变化率;
第七步、计算模块根据式(2)计算每个粒子在任意位置x点处的体积力FA2,并且入口处的体积力FA1 由(3)式确定;
第八步、由式(5)、式(9)和式(10)分别计算出动态压力、动态压力梯度和粘度;
p=c2ρ (5)
其中,c为人工声速,mb为粒子a的相邻粒子b的质量,μa为,μb分别表示粒子a,b的动态粘度,▽aWab为粒子b对粒子a产生影响的光滑函数的梯度,|rab|为粒子a,b之间的距离,η=0.1h,h为光滑长度。
第九步、再判断式(1)是否有效,如果有效则执行第十步,否则输出模块报错;
第十步、计算流体动量和能量的增加值;
第十一步、更新位置、速度、能量、密度的数值;
第十二步、判断雷诺数Re是否小于1,如果小于1则执行第十三步,否则输出模块报错。
第十三步、输出模块输出步骤S3中任意位置x点处的体积力FA2。
为了检验本发明的方法,下面通过它去模拟三个低雷诺数不可压缩流,包括Poiseuille流、局部膨胀管道流和倾斜平板流。为了比较,这些例子同样用常数体积力的WCSPH方法进行模拟,并且模拟的流体都是水,具有相同的初始密度ρ0=103kg/m3;同时,三个实施例都给出了理论解以供比较。
实施例一
本实施例采用Poiseuille流,即将两块无限大的平板分别放在坐标y=-0.5YL 和y=0.5YL处,其中,YL是两块平板之间的距离。其初始静止的流体由平行于X轴的体积力F(相应于入口与出口的压差)驱动,最后将达到一个稳定状态。
在本实施例中,YL=10-3m,流场长度为d=5×10-4m,静水压力差为△p =10-4 N/m2,通过式(1)、式(2)和式(3)模拟的结果如图3所示,并且从图3中可以看出,本发明的方法与用常数体积力模拟的结果相一致,且都与理论解吻合良好。
实施例二
本实施例采用局部膨胀管道流,如图4所示,可以考虑一段轴对称的血管,其压差驱动流体在管道内流动,压差为△p= p1- p2= 1.939006287×10-3 N/m2,管道的半径是关于位置x的函数,可以表示为(4)式:
其中,R0=0.5×10-3m,δ=0.2R0, L=6×10-3m, X0=2×10-3m,ε=δ/X0=0.05。
模拟结果,X轴速度和Y轴速度分别如图5和图6所示。并且从图5和图6中可以看出,采用本发明的方法的模拟结果与理论解吻合良好,但采用常数体积力的模拟结果却不能与理论解吻合。
实施例三
本实施例采用倾斜平板流,如图7所示,在本实施例中,dBC=4mm, 2l1=0.5mm, α=3.503°,且B与C之间的静水压力差为△p =1.21665968×10-3 N/m2。
采用本发明的方法,沿X轴上的流体粒子的水平速度的模拟结果如图8所示;从图8中可以看出,用本发明的方法的模拟结果再一次与理论解相符,但常数体积力情况下的模拟却不能相符。
虽然结合附图描述了本发明的实施方式,但是专利所有者可以在所附权利要求的范围之内做出各种变形或修改,只要不超过本发明的权利要求所描述的保护范围,都应当在本发明的保护范围之内。
Claims (4)
1.一种弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法,该方法是基于计算模块和输出模块,其特征在于:在该方法中,计算模块给定一低雷诺数不可压缩流,该低雷诺数不可压缩流在具有曲线边界的管道内流动,所述的具有曲线边界的管道为轴对称的、非平直的并且管道壁为固壁边界的管道,且该管道的出入口两端压差为△p,其中△p=p1-p2,p1为入口处的压力、p2为出口处的压力,并采用以x和r分别表示轴向坐标和径向坐标的柱坐标系,该方法具体包括以下步骤,
步骤S1、计算模块对于上述给定的低雷诺数不可压缩流,判断公式(1)是否成立,
或R0<<x0,L>>R0,Re<1 (1);
其中,δ为壁厚,x0为任意一点的轴向坐标,R0为管道平直处的半径,L为管道的长度,Re为雷诺数;
步骤S2、若步骤S1中判断公式(1)成立,计算模块在模拟流场中生成或继承粒子,并采用多边界切线技术处理固壁边界,根据该技术的需要随之生成虚拟粒子;
步骤S3、计算模块根据公式(2)计算每个粒子在任意位置x点处的体积力FA2,并且入口处的体积力FA1由公式(3)确定;
f(r,x)=[r2-R2(x)];
其中,dA表示面积元,ρ为流体密度,Δp为压力差,A1表示流场入口处垂直于X轴的截面,A2表示流场中任意坐标处垂直于X轴的截面,R(x)为管道曲线边界处的半径;
步骤S4、输出模块输出步骤S3中任意位置x点处的体积力FA2。
2.根据权利要求1所述的弯曲边界上低雷诺数不可压缩流中压力差的仿真方法,其特征在于:在所述步骤S2和S3之间还包括,
步骤S21、计算模块搜索相邻的粒子,并计算光滑函数;
步骤S22、计算模块确定是否用求和方法计算密度值,如果是则执行步骤S23,否则执行步骤S24;
步骤S23、由式(6)或式(7)得出密度值;
其中,mb为粒子b的质量,Wab为粒子b对粒子a产生影响的光滑函数,ra代表粒子a的位置向量,rb代表粒子b的位置向量,h表示光滑长度,ρb表示b的密度,h表示光滑长度;步骤S24、由式(8)计算出密度变化率;为粒子b对粒子a产生影响的光滑函数的梯度;
其中:ua,ub分别表示粒子a,b的速度矢量。
3.根据权利要求2所述的弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法,其特征在于:在步骤S3和S4之间还包括,
步骤S31、由式(5)、式(9)和式(10)分别计算出动态压力、动态压力梯度和粘度;
p=c2ρ (5)
其中,c为人工声速,mb为粒子a的相邻粒子b的质量,μa、μb分别表示粒子a,b的动态粘度,为粒子b对粒子a产生影响的光滑函数的梯度,|rab|为粒子a,b之间的距离,η=0.1h,h为光滑长度,μ表示层流的动态粘度,代表做拉普拉斯运算,表示压力的梯度;
步骤S32、再判断式(1)是否有效,如果有效则执行步骤S33,否则输出模块报错;
步骤S33、计算流体动量和能量的增加值;
步骤S34、更新位置、速度、能量、密度的数值;
步骤S35、判断雷诺数Re是否小于1,如果小于1则输出模块输出任意位置x点处的体积力FA2,否则输出模块报错。
4.根据权利要求3所述的弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法,其特征在于:所述公式(2)是对管道内的流体满足连续性和运动方程进行量级计算得出。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310428248.3A CN103473418B (zh) | 2013-09-18 | 2013-09-18 | 弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310428248.3A CN103473418B (zh) | 2013-09-18 | 2013-09-18 | 弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103473418A CN103473418A (zh) | 2013-12-25 |
CN103473418B true CN103473418B (zh) | 2017-02-08 |
Family
ID=49798266
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310428248.3A Expired - Fee Related CN103473418B (zh) | 2013-09-18 | 2013-09-18 | 弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103473418B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105587959B (zh) * | 2016-01-18 | 2018-08-28 | 北京航空航天大学 | 一种控制低雷诺数u型弯管流动分离的变截面方法 |
CN108446419B (zh) * | 2018-01-17 | 2022-02-22 | 天津大学 | 一种超声速边界层特征厚度估算方法 |
CN110569541A (zh) * | 2019-08-01 | 2019-12-13 | 天津大学 | 管道水锤分析方法 |
CN110705185A (zh) * | 2019-09-29 | 2020-01-17 | 天津大学 | 预测管道气锤的方法 |
Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101484922A (zh) * | 2006-04-05 | 2009-07-15 | 财团法人Seoul大学校产学协力财团 | 使用导数质点来模拟流体详细运动的方法 |
-
2013
- 2013-09-18 CN CN201310428248.3A patent/CN103473418B/zh not_active Expired - Fee Related
Patent Citations (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101484922A (zh) * | 2006-04-05 | 2009-07-15 | 财团法人Seoul大学校产学协力财团 | 使用导数质点来模拟流体详细运动的方法 |
Also Published As
Publication number | Publication date |
---|---|
CN103473418A (zh) | 2013-12-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Porter et al. | Multicomponent interparticle-potential lattice Boltzmann model for fluids with large viscosity ratios | |
CN107590336B (zh) | 燃气管道泄漏对内部流场影响的数值模拟方法 | |
Du et al. | Multi-relaxation-time lattice Boltzmann model for incompressible flow | |
Xie et al. | A mesoscale 3D CFD analysis of the liquid flow in a rotating packed bed | |
Song et al. | Experimental study of a separating, reattaching, and redeveloping flow over a smoothly contoured ramp | |
CN103473418B (zh) | 弯曲边界上的低雷诺数不可压缩流中的压力差的仿真方法 | |
Gungor et al. | Scaling and statistics of large-defect adverse pressure gradient turbulent boundary layers | |
Ndinisa et al. | Computational fluid dynamics simulations of Taylor bubbles in tubular membranes: model validation and application to laminar flow systems | |
Tijsseling et al. | Improved one-dimensional models for rapid emptying and filling of pipelines | |
Epely-Chauvin et al. | Numerical modelling of plunge pool scour evolution in non-cohesive sediments | |
Bazilevs et al. | Computation of the flow over a sphere at Re= 3700: A comparison of uniform and turbulent inflow conditions | |
Amiri-Jaghargh et al. | DSMC simulation of low knudsen micro/nanoflows using small number of particles per cells | |
Rashidi et al. | Numerical and experimental study of a ventilated supercavitating vehicle | |
Araújo et al. | Simulation of slug flow systems under laminar regime: Hydrodynamics with individual and a pair of consecutive Taylor bubbles | |
Zeng et al. | A hybrid RANS-LES model for combining flows in open-channel T-junctions | |
Zhou | A lattice Boltzmann method for solute transport | |
Shams et al. | Prediction of small-scale cavitation in a high speed flow over an open cavity using large-eddy simulation | |
Rumsey et al. | Time-accurate computations of isolated circular synthetic jets in crossflow | |
Dias et al. | Application of URANS-VOF models in hydrodynamics study of oscillating water column | |
Coussirat et al. | Study of available turbulence and cavitation models to reproduce flow patterns in confined flows | |
Ilyushin et al. | Large eddy simulation of the turbulent round jet dynamics | |
Ebrahimi et al. | A computational study of turbulent airflow and tracer gas diffusion in a generic aircraft cabin model | |
Isoz | CFD study of gas flow through structured separation columns packings Mellapak 250. X and Mellapak 250. Y | |
CN105973319A (zh) | 一种控制棒驱动机构排污系统水力特性计算方法 | |
Smith | BASIC hydraulics |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170208 Termination date: 20210918 |