CN105653860B - 可压缩气体与不可压缩液体多介质界面追踪数值方法 - Google Patents
可压缩气体与不可压缩液体多介质界面追踪数值方法 Download PDFInfo
- Publication number
- CN105653860B CN105653860B CN201511022502.5A CN201511022502A CN105653860B CN 105653860 B CN105653860 B CN 105653860B CN 201511022502 A CN201511022502 A CN 201511022502A CN 105653860 B CN105653860 B CN 105653860B
- Authority
- CN
- China
- Prior art keywords
- interface
- incompressible
- multimedium
- compressible
- liquid
- 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
- 238000000034 method Methods 0.000 title claims abstract description 91
- 239000007788 liquid Substances 0.000 title claims abstract description 59
- 230000035939 shock Effects 0.000 claims abstract description 24
- 239000000463 material Substances 0.000 claims abstract description 19
- 230000002123 temporal effect Effects 0.000 claims abstract description 10
- 238000004088 simulation Methods 0.000 claims abstract description 5
- 238000013213 extrapolation Methods 0.000 claims abstract description 3
- 238000012795 verification Methods 0.000 abstract 1
- 239000007789 gas Substances 0.000 description 49
- 238000004364 calculation method Methods 0.000 description 13
- 238000007596 consolidation process Methods 0.000 description 8
- 239000012530 fluid Substances 0.000 description 6
- 238000010276 construction Methods 0.000 description 5
- 230000000694 effects Effects 0.000 description 5
- 230000008901 benefit Effects 0.000 description 3
- 230000010355 oscillation Effects 0.000 description 3
- 230000000644 propagated effect Effects 0.000 description 3
- 230000004907 flux Effects 0.000 description 2
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 2
- 239000010931 gold Substances 0.000 description 2
- 229910052737 gold Inorganic materials 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 230000007480 spreading Effects 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 230000003044 adaptive effect Effects 0.000 description 1
- 230000005540 biological transmission Effects 0.000 description 1
- 238000000354 decomposition reaction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 239000003292 glue Substances 0.000 description 1
- 230000003993 interaction Effects 0.000 description 1
- 230000009022 nonlinear effect Effects 0.000 description 1
- 238000005192 partition Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16Z—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
- G16Z99/00—Subject matter not provided for in other main groups of this subclass
Landscapes
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开一种可压缩气体与不可压缩液体多介质界面追踪数值方法,单一的气体和液体介质区域均采用间断伽辽金方法进行计算,同时由于多介质流场的时间步长一般是非均匀变化的,本发明在采用间断伽辽金方法求解不可压缩Navier‑Stokes方程时,推导出基于非均匀时间步长的时间离散格式。多介质界面边界条件采用修正虚拟流体方法进行定义。在物质界面法向构造一种新的可压缩与不可压缩黎曼问题来预测界面状态,预测的界面状态不仅用于推进多介质界面,而且可以直接更新界面附近可压缩气体一侧的密度并外推得到对应的虚拟流体状态。通过多个经典数值算例的验证,发现对于气体中出现激波的多介质问题本发明能得到更为精确的数值模拟结果。
Description
技术领域:
本发明涉及一种可压缩气体与不可压缩液体多介质界面追踪数值方法,其属于计算流体力学,多相流数值模拟领域。
背景技术:
由于可压缩气体与不可压缩液体的物质属性存在较大差异,数值模拟相关气液多介质问题存在较大困难。在早期的研究中,为降低计算难度气体和液体均作为可压缩介质并且液体的状态方程为Tait方程。然而,由于界面两侧的物质属性差异较大,在物质界面处易产生非物理振荡,而且由于液体中的声速较大会导致计算时间步长过小从而降低计算效率。当气体和液体均作为不可压缩介质时,对于含有激波的高速流动问题难以求解,因而这一方法只适应于求解低速多介质流动问题。而当气体作为可压缩介质而液体作为不可压缩介质时,可以较好地考虑到物质属性的差异以及能够处理气体中含有激波的高速流动问题。当然这一方法的前提是正确地处理多介质界面。
近二十年来,间断伽辽金(DG)方法在求解单介质流场上逐渐成为研究热点。在对DG方法进行的研究中,较著名的是由Cockburn和Shu等人提出的龙格库塔间断伽辽金(RKDG)方法:在时间离散上采用总变差减小的(TVD)高阶龙格库塔(RK)方法,在空间离散上通过求解精确或近似黎曼问题获得网格界面处的数值通量,应用总变差有界(TVB)限制器控制数值解的虚假振荡。相比之下,DG方法在不可压缩流场的求解上起步较晚。Liu等人最初以涡量函数的形式研究高阶不可压缩DG方法。后来Cockburn等人提出了求解Stokes,Oseen,Navier-Stokes方程的局部间断伽辽金(LDG)方法。Ferrer等人提出了一种更简单高效的求解非定常不可压Navier-Stokes方程的DG方法,该方法采用Interior Penalty(IP)空间离散方法求解椭圆形方程。在多介质界面处理上,早期基于γ模型、质量分数模型、level set函数的方法在多介质界面处往往会产生较大的数值耗散。然而,对于界面追踪方法,界面是显式追踪的而且界面位置在计算过程中能够始终保持非常清晰。原始的虚拟流体方法(GFM)主要用于定义界面边界条件从而使得每种介质可以单独求解。但是当压力或者速度在界面附近存在较大的梯度时,该方法的表现效果较差。后来,在处理水气多介质界面问题时,Fedkiw等提出新的虚拟流体方法(new GFM)。在该方法中,虚拟流体速度取决于液体而虚拟流体压力由气体决定。New GFM的优点是简单且易于推广到高维情况。然而界面边界条件应当非线性地取决于界面两侧的物质属性,所以直接把真实流体状态作为虚拟流体状态并不是一种严格的方法。而且,对于强激波打击物质界面问题,应当首先预测界面状态来考虑界面处波系的相互作用。其后又有许多学者提出了一些改进版本的虚拟流体方法,例如修正虚拟流体方法(MGFM)。该方法通过求解界面黎曼问题来定义虚拟流体状态。MGFM在处理可压缩多介质问题上是非常有效的。
当可压缩气体与不可压缩液体相互作用时,由于物质界面两侧控制方程和物质属性存在较大的差异,需要给出更为严格的物质界面边界条件。New GFM给出的界面边界条件过于简单,没有充分考虑到界面两侧物质属性的非线性影响。MGFM通过在物质界面处构造黎曼问题来预测界面状态,在求解可压缩多介质问题上取得了较大的成功。然而当液体作为不可压缩介质时,需要在界面处构造可压缩与不可压缩黎曼问题。所以本发明在界面的法向构造出新的可压缩与不可压缩黎曼问题来预测界面状态,并且预测的界面状态可以直接用于推进界面以及定义虚拟流体状态。在定义了物质界面边界条件后,可压缩气体和不可压缩液体分别采用DG方法进行求解。同时由于现有的不可压缩DG方法的时间步长为均匀不变的,这往往难以适应多介质问题变时间步长的要求,所以本发明推导出基于非均匀时间步长的不可压缩时间离散格式。采用可压缩与不可压缩DG方法的目的是由于其具有较好的空间紧致性,这使得所需定义的虚拟流体状态数目大为减少,在处理复杂多介质界面问题上具有较大优势。
发明内容:
为了克服现有技术中的不足,本发明提供了一种可压缩气体与不可压缩液体多介质界面追踪数值方法。
本发明采用如下技术方案:一种可压缩气体与不可压缩液体多介质界面追踪数值方法,其包括如下步骤:
步骤1:利用一系列标志点离散多介质界面,在标志点法向构造可压缩与不可压缩黎曼问题预测界面状态;
步骤2:直接利用标志点上预测的界面状态来更新界面附近可压缩气体一侧的密度,并分别外推得到对应的虚拟流体状态;
步骤3:采用龙格库塔间断伽辽金方法求解可压缩气体,对于不可压缩Navier-Stokes方程,推导基于非均匀时间步长的时间离散格式,采用不可压缩间断伽辽金方法对液体介质进行计算;
步骤4:由标志点上预测的界面状态推进多介质界面,重构界面,得到新的界面位置和标志点并根据新的界面确定整个多介质流场的数值解;
步骤5:判断是否满足计算终止条件,如不满足此条件,返回步骤1继续进行下一时刻的数值计算。
进一步地,在可压缩与不可压缩黎曼问题的构造上,把不可压缩液体当作声速趋于无穷大的可压缩液体,当所有介质为可压缩时,状态方程可以统一表示为:
p=(γ-1)ρe-γB
其中B是物质参数,对于理想气体,B为零,考虑黎曼问题的初始条件:
其中QL=[ρL,uL,pL,γL,BL]T,QR=[ρR,uR,pR,γR,BR]T,不失一般性,假设左侧为气体,右侧为液体,首先采用双激波近似方法求解该黎曼问题:
其中:
而预测的界面速度定义为:
其中:
由于液体的声速:
假设为无穷大,但是ρR为有限值,所以:
kR=γR(pR+BR)→+∞
将其代入界面速度的定义式以及基于双激波近似方法得到的黎曼解中,可以得到:
fR(pI,QR)=0
对于左侧气体有:
uI=uL-fL(pI,QL)
可以得到:
uI=uR
进一步地,采用不可压缩间断伽辽金方法求解液体介质时,推导出一种基于非均匀时间步长的时间离散格式,对Navier-Stokes方程的时间导数项采用二阶后退差分离散:
其中γ0=α0+α1,下标“n-1”,“n”和“n+1”分别表示tn-1,tn和tn+1时刻,参数γ0,α0,α1为:
其中Δtn-1=tn-tn-1,Δtn=tn+1-tn,时间导数项可以进一步表示为:
其中f是Navier-Stokes方程的非线性项、粘性项和压力项,非线性项采用显式处理,粘性项和压力项采用隐式处理,可以得到不可压缩Navier-Stokes方程的时间离散格式为:
其中:
通过设置Δtn≠Δtn-1,即可得到基于非均匀时间步长的不可压缩Navier-Stokes方程的时间离散格式,满足多介质数值模拟中非均匀时间步长的要求。
本发明具有如下有益效果:
(1)通过构造新的可压缩与不可压缩黎曼问题来预测物质界面状态,可以在液体为不可压缩时应用MGFM定义更为精确的界面边界条件;
(2)推导出基于非均匀时间步长的不可压缩时间离散格式,使之能够应用于多介质中不可压缩流体的数值计算;
(3)采用空间紧致性较好的可压缩与不可压缩DG方法进行计算使得在物质界面附近所需虚拟流体状态信息大为减少,提高了流场计算的精确性;
(4)预测的物质界面状态不仅用于推进界面,而且可直接用于定义虚拟流体状态,编程简单,计算效率较高,且易于向高维推广。
附图说明:
图1:标志点处黎曼问题的构造。
图2:可压缩流体密度的更新。
图3:数值结果(算例8.1):压力(左)和速度(右)。
图4:数值结果(算例8.2):压力(左)和速度(右)。
图5:数值结果(算例8.3):压力(左)和速度(右)。
图6:数值结果(算例8.4):压力(左)和速度(右)。
图7:数值结果(算例8.5):压力(左)和速度(右)。
图8:不可压密度ρ=1000kg/m3的数值结果(算例8.6):压力(左)和x方向速度(右)。
图9:不可压密度ρ=1000kg/m3的压力云图(算例8.6)。
图10:不可压密度ρ=10kg/m3的数值结果(算例8.6):压力(左)和x方向速度(右)。
图11:不可压密度ρ=10kg/m3的界面形状(算例8.6)。
图12:不可压密度ρ=1000kg/m3的数值结果(算例8.7):压力(左)和x方向速度(右)。
图13:不可压密度ρ=1000kg/m3的压力云图(算例8.7)。
图14:不可压密度ρ=10kg/m3的数值结果(算例8.7):压力(左)和x方向速度(右)。
图15:不可压密度ρ=10kg/m3的界面形状(算例8.7)。
具体实施方式:
下面结合附图和算例对发明内容作进一步说明。
(1)控制方程
可压缩气体的控制方程为Euler方程:
其中:U=[ρ,ρu,ρv,E]T,F(U)=[F1(U),F2(U)],F1(U)=[ρu,ρu2+p,ρuv,(E+p)u]T,F2(U)=[ρv,ρuv,ρv2+p,(E+p)v]T。ρ是密度,u是x方向速度,v是y方向速度,p是压力,E是单位体积的总能量,其定义为:
E=ρe+ρ(u2+v2)/2 (2)
其中e是单位质量的内能。状态方程为:
p=(γ-1)ρe (3)
对于理想气体γ是比热比。
不可压缩液体的控制方程为Navier-Stokes方程:
其中,V=(u,v)是速度矢量,μ是粘性系数。ρ和μ在不可压缩液体中假设为常数。Navier-Stokes方程的解记为:
(2)自适应时间步长
多介质全流场的时间步长取为可压缩气体与不可压缩液体时间步长的最小值:
Δt=min(Δtc,Δti) (6)
对于可压缩气体,时间步长定义为:
其中Cf是CFL系数,a是局部声速,Δx和Δy是网格步长。对于不可压缩液体,时间
步长定义为:
其中是特征长度(网格尺寸),是特征速度,ki是不可压缩DG方法的多项式阶数。
(3)求解Euler方程的RKDG方法
对于矩形网格Ii,j=[xi-1/2,xi+1/2]×[yj-1/2,yj+1/2],网格中心记为(xi,yj),网格步长为Δx=xi+1/2-xi-1/2,Δy=yj+1/2-yj-1/2。数值解Uh(x,y,t)和基函数由给定(其中是网格Ii,j上多项式次数小于或等于kc的多项式集合)。网格Ii,j上的当地正交基函数定义为:
数值解Uh(x,y,t)可以表示为:
其中定义为:
其中把式(10)代入式(1),对式(1)乘以基函数并在网格Ii,j上分部积分,可以得到:
其中是网格Ii,j上边e的单位法向量。采用Lax-Friedrichs通量近似。右端积分项采用高斯求积法求解。式(12)的半离散形式表示为:
Ut=L(U) (13)
采用三阶TVD RK方法进行时间离散:
当流场出现间断时应用斜率限制器克服可能出现的虚假振荡问题。对于标量方程,记:
其中:
对进行修正:
其中是修正minmod函数:
其中是TVB限制器常量。对于和上式中Δx应替换为Δy。对于方程组,需要对限制器进行当地特征分解。
(4)求解Navier-Stokes方程的DG方法
不可压缩区域记为Ω,其边界可以为Dirichlet类型(ΓD)或者Neumann类型(ΓN)等。对不可压缩区域进行矩形划分Ωh={K},外边界记为内边界记为Γh。数值解Wh(x,y,t)和基函数由给定(其中是网格K上多项式次数小于或等于ki的多项式集合)。记网格K=Ii,j,数值解Wh(x,y,t)可以表示为:
其中是多项式系数,基函数的定义与式(9)一致。在现有的不可压缩DG方法中,求解单介质流动问题的时间步长是均匀不变的。然而对于多介质问题,总的时间步长(式(6))应由可压缩介质与不可压缩介质共同决定,所以有必要推导出基于非均匀时间步长的不可压缩时间离散格式。对Navier-Stokes方程的时间导数项采用二阶后退差分离散:
其中γ0=α0+α1。下标“n-1”,“n”和“n+1”分别表示tn-1,tn和tn+1时刻。参数γ0,α0,α1为:
其中Δtn-1=tn-tn-1,Δtn=tn+1-tn。式(19)进一步表示为:
其中f是式(4)的右端项。非线性项采用显式处理,粘性项和压力项采用隐式处理。这样可以得到不可压缩Navier-Stokes方程的时间离散格式:
其中:
可以发现,如果令Δtn≠Δtn-1,式(24)即为非均匀时间步长的时间离散格式。引入中间变量和式(24)被分裂成:
对式(28)求散度并考虑得到:
式(29)重新表示为:
由式(30)求解更新中间变量并根据式(31)求解Vn+1。采用对称IP伽辽金(SIPG)方法求解式(30)和(31):
其中,q∈H1(Ω)为标量函数(可以直接推广到矢量函数),α是波数,n是上外法向矢量,g∈L2(Ω),LD∈H1/2(ΓD)和LN∈L2(ΓN)表示边界条件。式(32)的弱解可以通过下式求得:
其中,ne是边e的外法向矢量,[[·]]和{{·}}表示相邻网格K1和K2之间的跳跃值和平均值:
在边界上有:
参数λ定义为:
其中ηK为网格K的面积与周长之比。
(5)界面追踪方法
如图1所示,在tn时刻可压缩气体(介质1)和不可压缩液体(介质2)被界面Π分开。标记点是界面和网格线的交点。以标记点E为例,点A(xA,yA)和点B(xB,yB)在不同介质中,与标记点E相距为Δn。这里:
xA=xE-Δn·NEx,yA=yE-Δn·NEy
(40)
xB=xE+Δn·NEx,yB=yE+Δn·NEy
其中是标记点E的单位法向矢量。点A和点B的状态由式(10)和式(18)得到:
其中KA和KB分别是包含点A和点B的网格。计算点A和点B的密度、法向速度和压力,分别表示为和沿着标志点E的法线方向构造可压缩与不可压缩黎曼问题,其初始条件为:
求解该黎曼问题并记其解为其中上标I表示界面,下标L和R表示界面左右两侧。由于不可压缩介质的刚性较强,界面的切向速度应由液体决定。此处标志点E的切向速度定义为:其中是点B的切向速度。这样标志点E以速度vf=(uI,vI)被推进到下一时刻。其他标志点的推进方式是类似的。
(6)可压缩与不可压缩黎曼问题
在修正虚拟流体方法中,所有的介质为可压缩的,通过构造可压缩黎曼问题来预测界面状态。然而当液体作为不可压缩介质时,应当在界面处构造可压缩与不可压缩黎曼问题。考虑到当马赫数趋于零时,可压缩流体的解趋于不可压缩流体的解。本发明把不可压缩液体当作声速趋于无穷大的可压缩液体来构造可压缩与不可压缩黎曼问题。当所有介质为可压缩时,状态方程可以统一表示为:
p=(γ-1)ρe-γB (43)
其中B是物质参数。对于理想气体,B为零,上式退化为式(3)。考虑黎曼问题的初始条件:
其中QL=[ρL,uL,pL,γL,BL]T,QR=[ρR,uR,pR,γR,BR]T。不失一般性,假设左侧为气体,右侧为液体。首先采用双激波近似方法求解该黎曼问题:
其中:
而预测的界面速度定义为:
其中:
由于液体的声速:
假设为无穷大,但是ρR为有限值,所以:
kR=γR(pR+BR)→+∞ (53)
将其代入式(46),(48),(51),(47)和(49)得到:
fR(pI,QR)=0 (56)
对于左侧气体有:
uI=uL-fL(pI,QL) (59)
可以得到:
uI=uR (60)
对于式(57)通过迭代方法求解pI。可由式(45)得到。在上述可压缩与不可压缩黎曼问题的解中,有这一结果满足不可压缩介质密度保持不变的特点。uI=uR表明了多介质界面的速度应当由不可压缩介质决定。
(7)修正虚拟流体方法
在修正虚拟流体方法中,在界面附近通过构造黎曼问题来预测界面状态。由于在界面追踪部分已经构造标志点处的黎曼问题,预测的界面状态可以直接用于定义虚拟流体状态。如图2所示,可压缩气体(介质1)和不可压缩液体(介质2)被界面Π分开。是标志点的法向矢量,是界面附近可压缩气体中网格P的法向矢量。网格P的密度可由与其法向夹角最小的标志点上的界面状态进行更新。如果标志点E是与网格P法向夹角最小的标志点,那么对于标记点E通过等熵修正更新网格P密度的平均值。界面附近其它可压缩气体的密度更新方法类似。虚拟流体状态通过求解如下方程得到:
其中表示密度,法向速度,切向速度和压力。如图2所示,在求解介质1的边界条件时采用“+”号,否则采用“-”号。采用迭代方法求解该方程,空间上采用一阶迎风格式离散,时间上采用三阶TVD RK方法进行离散。由于DG方法空间紧致性很好,此处仅采用2~3个网格作为虚拟流体区域。通过把虚拟流体状态投影到基函数空间,得到虚拟网格状态的平均值,这样即可确定界面边界条件。
(8)新方法验证
数值算例采用均匀结构网格。可压缩与不可压缩DG方法的多项式阶数均取为2。式(7)中的CFL数Cf取为0.18。RKDG方法中的TVB限制器常量取为0.1。比热比γ取为1.4,粘性系数μ=0.001137kg/(m s)。数值算例中同时给出了基于new GFM的计算结果作为比较。为了进行更为细致的比较,在压力和速度结果下方给出了计算结果的局部放大图。
一维算例
所有一维算例的计算域为[0m,1m],网格数为200。如果气体和液体均作为可压缩介质,算例8.1~8.3实际上为可压缩黎曼问题,可以得到其相应的解析解以作进一步比较。
算例8.1.气液界面位于x=0.4m,初始条件为:气体:ρ=1kg/m3,u=1000m/s,p=1×105Pa,液体:ρ=1009kg/m3,u=0m/s,p=2×107Pa。图3给出了t=2.8×10-4s时压力和速度的计算结果。观察可压缩黎曼问题的解析解,可以发现在液体中形成了稀疏波(但在速度结果中不是很明显)。然而,当液体为不可压缩介质时,稀疏波立刻传出计算域。基于MGFM得到的激波位置与可压缩解析解基本一致,而基于new GFM得到的激波在可压缩气体中传播较慢。这是由于在MGFM中,通过求解可压缩与不可压缩黎曼问题预测界面状态。对于在可压缩气体中产生激波的问题,界面状态事先被预测出来。这导致在计算的最初几个时间步中产生强度较大的激波,所以基于MGFM得到的激波传播较快。
算例8.2.气液界面位于x=0.4m,初始条件为:气体:ρ=1kg/m3,u=800m/s,p=1×105Pa,液体:ρ=1000kg/m3,u=0m/s,p=1×105Pa。图4给出了t=2.8×10-4s时压力和速度的计算结果。可以发现在可压缩黎曼问题的解析解中,液体中形成了激波(但在速度结果中不是很明显)。但是当液体为不可压缩介质时,该激波立刻以无穷大的速度传出计算域。基于MGFM得到的激波位置与可压缩解析解基本一致,而基于new GFM得到的激波在可压缩气体中传播较慢。
算例8.3.气液界面位于x=0.5m,初始条件为:气体:ρ=10kg/m3,u=-1000m/s,p=1×107Pa,液体:ρ=1002kg/m3,u=0m/s,p=5×106Pa。图5给出了t=2×10-4s时压力和速度的计算结果。在可压缩黎曼问题的解析解中,稀疏波仍然在液体中传播(但在速度结果中不是很明显)。可以发现,当在可压缩气体中出现稀疏波时,这些方法得到的稀疏波位置基本一致。
算例8.4.液滴位于0.4<x<0.6的区域中。初始条件为:气体:ρ=1.226kg/m3,u=0m/s,p=1×105Pa,液滴:ρ=1000kg/m3,u=100m/s,p=1×105Pa。图6给出了t=7.5×10-4s时压力和速度的计算结果。由于液滴向右运动,在右侧的气体中会形成激波,而在左侧的气体中会形成稀疏波。基于MGFM和new GFM得到的稀疏波位置基本一致,而基于MGFM得到的激波位置几乎比new GFM领先一个网格步长。
算例8.5.气体中位于x=0.1m的右行激波打击位于0.4<x<0.6的液滴。初始条件为:波前气体:ρ=1.58317kg/m3,u=0m/s,p=98066.5Pa,波后气体:ρ=2.124kg/m3,u=89.981m/s,p=148407.3Pa,液滴:ρ=1000kg/m3,u=0m/s,p=98066.5Pa。图7给出了t=1.75×10-3s时压力和速度的计算结果。由于液滴的刚性较强,在左侧的气体中会反射激波,而在液滴中形成速度很快、强度较弱的传递激波。可以发现基于MGFM得到的激波位置几乎仍然比new GFM领先一个网格步长。
二维算例
二维多介质问题的计算域为[0m,1m]×[0m,1m],网格数为100×100。半径为0.2m的圆形液滴位于计算域的中央。由于流场上下对称,在二维算例中仅仅模拟下面半个流场的流动问题。
算例8.6.初始条件为:气体:ρ=1.226kg/m3,u=0m/s,v=0m/s,p=1×105Pa,液滴:ρ=1000kg/m3,u=100m/s,v=0m/s,p=1×105Pa。图8给出了t=5×10-4s时沿着y=0.5m压力和x方向速度的计算结果。由于液滴向右运动,在右侧的气体中会形成激波,而在左侧的气体中会形成稀疏波。基于MGFM得到的激波位置几乎比new GFM领先一个网格步长。图9给出了相同时刻的压力云图(图中虚线表示界面位置)。可以看出由于液滴的刚性较大,界面的形状基本保持不变。为了进一步考察液滴的密度对于计算结果的影响,这里将液滴密度替换为ρ=10kg/m3,压力和x方向速度的计算结果如图10所示。基于new GFM得到的激波位置几乎落后一个网格步长。图11给出在t=2.5×10-3s、网格尺寸为0.02m,0.1m,0.005m时的界面形状(图中虚线表示初始时刻界面形状)。可以看出较轻的液滴更易发生变形且速度降低较快。密网格上已经出现了Kelvin-Helmholtz不稳定现象。这主要是因为在界面处没有考虑粘性力的影响从而不能保证界面两侧切向速度连续。然而粗网格的数值粘性在一定程度上弱化了这一现象。对上述不同网格尺寸下的液滴损失率进行研究,计算结果为:0.0065,0.0032,0.0011。可以看出随着网格的不断细化,水滴的损失率逐渐变小。
算例8.7.气体中位于x=0.1m的右行激波打击液滴。初始条件为:波前气体:ρ=1.58317kg/m3,u=0m/s,v=0m/s,p=98066.5Pa,波后气体:ρ=2.124kg/m3,u=89.981m/s,v=0m/s,p=148407.3Pa,液滴:ρ=1000kg/m3,u=0m/s,v=0m/s,p=98066.5Pa。图12给出了t=1.25×10-3s时沿着y=0.5m压力和x方向速度的计算结果。基于new GFM得到的激波位置几乎落后一个网格步长。相同时刻的压力云图如图13所示(图中虚线表示界面位置),可以较清晰地看出入射激波和反射激波。由于液滴的刚性较强,液滴几乎保持静止。和算例8.6相似,这里把液滴的密度替换为ρ=10kg/m3,压力和速度的结果如图14所示。可以发现基于MGFM得到的激波位置领先一个网格步长。图15给出了在t=2.5×10-3s、网格尺寸为0.02m,0.1m,0.005m时的界面形状(图中虚线表示初始时刻界面形状)。通过观察发现,Kelvin-Helmholtz不稳定现象再次出现在细网格上。计算上述不同网格尺寸下的液滴损失率:0.0175,0.0041,0.003,可以看出损失率随着网格的细化而变小。
以上所述仅是本发明的优选实施方式,应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理的前提下还可以作出若干改进,这些改进也应视为本发明的保护范围。
Claims (3)
1.一种可压缩气体与不可压缩液体多介质界面追踪数值方法,其特征在于:包括如下步骤
步骤1:利用一系列标志点离散多介质界面,在标志点法向构造可压缩与不可压缩黎曼问题预测界面状态;
步骤2:直接利用标志点上预测的界面状态来更新界面附近可压缩气体一侧的密度,并分别外推得到对应的虚拟流体状态;
步骤3:采用龙格库塔间断伽辽金方法求解可压缩气体,对于不可压缩Navier-Stokes方程,推导基于非均匀时间步长的时间离散格式,采用不可压缩间断伽辽金方法对液体介质进行计算;
步骤4:由标志点上预测的界面状态推进多介质界面,重构界面,得到新的界面位置和标志点并根据新的界面确定整个多介质流场的数值解;
步骤5:判断是否满足计算终止条件,如不满足此条件,返回步骤1继续进行下一时刻的数值计算。
2.如权利要求1所述的可压缩气体与不可压缩液体多介质界面追踪数值方法,其特征在于:在可压缩与不可压缩黎曼问题的构造上,把不可压缩液体当作声速趋于无穷大的可压缩液体,当所有介质为可压缩时,状态方程可以统一表示为:
p=(γ-1)ρe-γB
其中p是压力,ρ是密度,e是单位质量的内能,γ和B是物质参数,对于理想气体,γ是比热比,B为零,考虑黎曼问题的初始条件:
其中Q0表示黎曼问题的初始条件,QL表示黎曼问题界面左侧的初始条件,QR表示黎曼问题界面右侧的初始条件,QL=[ρL,uL,pL,γL,BL]T,QR=[ρR,uR,pR,γR,BR]T,其中ρL、uL和pL分别表示黎曼问题界面左侧初始条件的密度、速度和压力,γL和BL是黎曼问题界面左侧初始条件的物质参数,ρR、uR和pR分别表示黎曼问题界面右侧初始条件的密度、速度和压力,γR和BR是黎曼问题界面右侧初始条件的物质参数,不失一般性,假设左侧为气体,右侧为液体,首先采用双激波近似方法求解该黎曼问题:
其中和分别表示预测的界面左侧的密度和预测的界面右侧的密度,pI表示预测的界面压力,WL和WR定义为:
而预测的界面速度定义为:
其中:
由于液体的声速:
假设为无穷大,但是ρR为有限值,所以:
kR=γR(pR+BR)→+∞
将其代入界面速度的定义式以及基于双激波近似方法得到的黎曼解中,可以得到:
fR(pI,QR)=0
对于左侧气体有:
uI=uL-fL(pI,QL)
可以得到:
uI=uR
采用迭代方法求解式可以得到预测的界面压力pI,代入式得到同时有uI=uR,这样即求得可压缩与不可压缩黎曼问题的解。
3.如权利要求1所述的可压缩气体与不可压缩液体多介质界面追踪数值方法,其特征在于:采用不可压缩间断伽辽金方法求解液体介质时,推导出一种基于非均匀时间步长的时间离散格式,对Navier-Stokes方程的时间导数项采用二阶后退差分离散:
其中V是速度矢量,γ0=α0+α1,上标“n-1”,“n”和“n+1”分别表示tn-1,tn和tn+1时刻,参数γ0,α0,α1为:
其中Δtn-1=tn-tn-1,Δtn=tn+1-tn,时间导数项可以进一步表示为:
其中f是Navier-Stokes方程的非线性项、粘性项和压力项,非线性项采用显式处理,粘性项和压力项采用隐式处理,可以得到不可压缩Navier-Stokes方程的时间离散格式为:
其中表示速度矢量梯度,表示压力梯度和密度之比,β0和β1分别为:
通过设置Δtn≠Δtn-1,即可得到基于非均匀时间步长的不可压缩Navier-Stokes方程的时间离散格式,满足多介质数值模拟中非均匀时间步长的要求。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201511022502.5A CN105653860B (zh) | 2015-12-30 | 2015-12-30 | 可压缩气体与不可压缩液体多介质界面追踪数值方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201511022502.5A CN105653860B (zh) | 2015-12-30 | 2015-12-30 | 可压缩气体与不可压缩液体多介质界面追踪数值方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105653860A CN105653860A (zh) | 2016-06-08 |
CN105653860B true CN105653860B (zh) | 2018-08-28 |
Family
ID=56490719
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201511022502.5A Expired - Fee Related CN105653860B (zh) | 2015-12-30 | 2015-12-30 | 可压缩气体与不可压缩液体多介质界面追踪数值方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105653860B (zh) |
Families Citing this family (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108153984B (zh) * | 2017-12-27 | 2021-04-13 | 中国空气动力研究与发展中心计算空气动力研究所 | 一种基于流场密度阶跃的高精度间断迦辽金人工粘性激波捕捉方法 |
CN108021533A (zh) * | 2017-12-29 | 2018-05-11 | 电子科技大学 | 一种基于广义坐标系求解任意色散材料电磁特性的方法 |
CN109002624B (zh) * | 2018-07-26 | 2022-11-08 | 上海交通大学 | 超声速刚性燃烧流动双自适应解耦优化模拟方法及系统 |
CN110309541B (zh) * | 2019-05-31 | 2023-04-07 | 中国航天空气动力技术研究院 | 一种变比热气体的不同介质多组份流场界面条件构建方法 |
CN112861263B (zh) * | 2021-02-22 | 2024-02-13 | 西北工业大学 | 一种适用于可压缩两相流的计算模拟方法 |
-
2015
- 2015-12-30 CN CN201511022502.5A patent/CN105653860B/zh not_active Expired - Fee Related
Non-Patent Citations (8)
Title |
---|
A front-tracking method with projected interface conditions for compressible multi-fluid flows;Hiroshi Terashima等;《Computers & Fluids》;20101231;第39卷(第10期);第1804-1814页 * |
MGFM 方法在运动固壁边界问题中的应用;徐娜 等;《科技风》;20110131;第203-204页 * |
一维多介质可压缩流动的守恒型界面追踪方法;王东红 等;《计算数学》;20090531;第31卷(第2期);第118-126页 * |
基于带限制器的RKDG法的可压缩流数值模拟;杨磊 等;《扬州大学学报(自然科学版)》;20140831;第17卷(第3期);第24-28页 * |
基于黎曼问题的界面追踪方法;刘希;《中国优秀硕士论文数据库 基础科学辑》;20150215;第2015年卷(第2期);A004-84 * |
复杂界面的界面追踪法;周袁媛;《中国优秀硕士论文数据库 基础科学辑》;20040915;第2004年卷(第3期);A004-31 * |
带TVB限制器的RKDG方法与浸入边界方法在Euler方程中的应用;倪建;《中国优秀硕士论文数据库 基础科学辑》;20150115;第2015年卷(第1期);A002-124 * |
运动界面的界面追踪法及其并行实现;周春晨;《中国优秀硕士论文数据库 基础科学辑》;20071215;第2007年卷(第6期);A004-36 * |
Also Published As
Publication number | Publication date |
---|---|
CN105653860A (zh) | 2016-06-08 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN105653860B (zh) | 可压缩气体与不可压缩液体多介质界面追踪数值方法 | |
Shetabivash et al. | Numerical analysis of droplet impact onto liquid film | |
De Rosis et al. | Aeroelastic study of flexible flapping wings by a coupled lattice Boltzmann-finite element approach with immersed boundary method | |
Karapetsas et al. | Efficient modelling of droplet dynamics on complex surfaces | |
Liu et al. | The dynamic behaviors of a bubble in a confined domain | |
Lee et al. | Fluid–structure interaction analysis on a flexible plate normal to a free stream at low Reynolds numbers | |
Niu et al. | Simulations of the shock waves and cavitation bubbles during a three-dimensional high-speed droplet impingement based on a two-fluid model | |
Xu et al. | Spread and recoil of liquid droplets impacting on solid surfaces with various wetting properties | |
Vincent et al. | Penalty and Eulerian–Lagrangian VOF methods for impact and solidification of metal droplets plasma spray process | |
Hu et al. | Modeling of keyhole dynamics and analysis of energy absorption efficiency based on Fresnel law during deep-penetration laser spot welding | |
Karakus et al. | An adaptive fully discontinuous Galerkin level set method for incompressible multiphase flows | |
Wang et al. | Asymmetric gravity–capillary solitary waves on deep water | |
Li et al. | Hydrodynamic analysis of the energy dissipation of droplets on vibrating superhydrophobic surfaces | |
Gao et al. | A high-order fluid–structure interaction framework with application to shock-wave/turbulent boundary-layer interaction over an elastic panel | |
Xie et al. | A numerical simulation of VIV on a flexible circular cylinder | |
Xu et al. | A geometry-adaptive immersed boundary–lattice Boltzmann method for modelling fluid–structure interaction problems | |
Fontes et al. | Numerical simulation of a water droplet splash: Effects of density interpolation schemes | |
Ma et al. | Modelling cavitating flows using an eulerian-lagrangian approach and a nucleation model | |
CN105183965B (zh) | 用于预测雾化过程的大涡模拟方法 | |
Parashar et al. | Lagrangian investigations of velocity gradients in compressible turbulence: lifetime of flow-field topologies | |
Ma et al. | Efficient solution of bimaterial Riemann problems for compressible multi-material flow simulations | |
Wang et al. | Modeling of droplet-based processing for the production of high-performance particulate materials using the level set method | |
Dong et al. | An experimental study about drag crisis phenomenon on teardrop model | |
Sotomayor-Zakharov et al. | A dynamic collection efficiency formulation implemented in the icing code DICEPS | |
Pournaderi et al. | The effect of the surface inclination on the hydrodynamics and thermodynamics of leidenfrost droplets |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
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: 20180828 |