R13/R26矩方法的介观尺度壁面边界条件
doi: 10.11887/j.cn.202403010
杨伟奇1 , 杨惠2
1. 国防科技大学 空天科学学院, 湖南 长沙 410073
2. 国防科技大学 计算机学院, 湖南 长沙 410073
基金项目: 国家自然科学基金资助项目(U1730247,12302382) ; 湖南省自然科学基金资助项目(2022JJ40542)
Wall boundary condition for the R13/R26 moment method at mesoscopic level
YANG Weiqi1 , YANG Hui2
1. College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073 , China
2. College of Computer Science and Technology, National University of Defense Technology, Changsha 410073 , China
摘要
宏观尺度NSF(Navier-Stokes-Fourier)、R13/R26矩方程边界条件在中、大克努森数Kn来流条件下计算精度大幅度降低,也极易发散。针对这一难题,提出R13/R26矩方程的介观尺度边界条件,在靠近壁面处重构速度分布函数,并输入介观尺度Boltzmann模型方程;基于离散速度法求解宏观参数,所得到的宏观参数作为R13/R26矩方程的壁面边界条件。仿真结果表明:基于介观尺度边界条件的R13/R26矩方法相较原边界条件计算精度最大提高59.84%,同时,所提出的边界条件将矩方法对Kn的适用范围拓展到1.0。
Abstract
Wall boundary conditions for the macroscopic equations, i.e. the NSF(Navier-Stokes-Fourier) equations, R13/R26 moment equations, lose their accuracy dramatically and are easy to diverge, especially in the middle and high Knudsen number regimes.To overcome these difficulties, a wall boundary condition for the R13/R26 moment method was proposed at the mesoscopic level. The velocity distribution function was reconstructed and feedback into the Boltzmann model equation in the near-wall region, and the wall boundary condition for the R13/R26 moment method was calculated on the basis of solving the Boltzmann equation with the discrete velocity method. Results indicate that: the proposed wall boundary condition is able to increase the computational accuracy up to 59.84% compared with the classical approach. Meanwhile, it is able to get the steady-state solution for the Knudsen number up to 1.0.
描述气体非平衡程度与稀薄程度的克努森数Kn定义为分子平均自由程λ与特征长度L0的比值[1]。通常情况下,稀薄气体动力学输运特性可以通过介观尺度Boltzmann方程来进行描述,但由于其高维度积分-微分特性以及复杂碰撞项结构,过大的计算量导致其难以应用在实际工程中[2]
同时,宏观方程与拓展流体力学方程均可以从介观尺度Boltzmann方程推导得到[2-3],两种主要的推导方法分别为:基于CE(Chapman-Enskog)的分布函数展开方法[4-5]与Grad提出的基于高阶矩的分布函数展开方法[6-7]
CE展开法是将分子速度分布函数展开为Kn的级数,进而构建本构关系,其零阶、一阶展开形式分别对应Euler、Navier-Stokes方程[3-4]。当来流Kn<0.1时,NSF(Navier-Stokes-Fourier)方程组与速度滑移、温度跳跃边界条件可以近似描述非平衡稀薄流的输运特性。随着来流Kn的增加,NSF方程组的计算精度可以通过引入CE的高阶展开来提高,例如:CE二阶、三阶展开形式分别对应Burnett与超Burnett方程组[5]。然而,Grad指出基于CE展开的速度分布函数,无论展开的级数有多高,其得到的宏观方程仅能描述近连续流的流动特性[8]。同时,Karlin等认为Burnett方程违反了Boltzmann方程中隐含的物理规律,而超Burnett方程没有改变这一物理特性[9]
因此,Grad提出基于Hermite多项式,将速度分布函数展开为高阶矩的级数,推导得到高阶矩的控制方程,形成Grad矩方法[7]。Reitebuch等利用Grad矩方法研究了平板间Couette流动[10]。Struchtrup等对Grad矩方法进行规整化处理,得到R13矩方法[311]。Gu等讨论了R13矩方法在微尺度稀薄流动中的Kn适用范围并指出:R13矩方法能够捕获Kn<0.25时气体的非平衡流动输运特性[12]。进一步,Gu等将规整化矩方法推广到R26矩方法[13],并成功应用至管道流动[1213-14]、热驱流动[15-16]、多孔介质流动[16-17]以及稀薄条件下的圆柱绕流模拟中[18]
一直以来,Maxwell壁面边界条件广泛应用在介观尺度Boltzmann方程中[19],但反观宏观尺度拓展流体力学方程,由于其控制方程与高阶矩较多且高阶矩没有明确的物理意义等,其壁面边界条件推导难度极大。Liu等使用极小值原理研究了基于最小熵增的矩方程边界条件[20]。Gu等基于介观尺度Maxwell壁面边界条件,开拓性地提出了R13矩方程组中宏观量的壁面边界条件,为矩方程的工程化应用铺平了道路[21]。然而,Torrilhon等发现Gu等对R13矩方程组进行了过度约束,导致靠近壁面处的切应力与温度剖面出现不规则异常,故其对R13矩方程边界条件进行了修正[22]。进一步,Gu等利用同样的方法为R26矩方法推导了壁面边界条件[13]
近期,研究人员在宏观与介观两个层面对比了NSF、R13、R26矩方程组与Boltzmann方程的计算精度与计算效率[1214-1523],深入对比了一阶、三阶、五阶Hermite多项式重构分布函数的重构误差,并分析了各宏观方程壁面处边界条件的计算误差[1423]。结果表明:对于微管道流动,矩方法在壁面边界处过度计算了18.7%的滑移速度与39.1%的热流,壁面边界条件的计算精度直接限制了矩方法在全流场的数值模拟精度[1]。经典矩方法壁面边界条件误差的来源在于:分布函数被截断为Hermite多项式的三阶或五阶形式,高阶矩信息丢失。基于此,本文提出R13/R26矩方法的介观边界条件,旨在引入更为精确的壁面边界处分布函数,以改进矩方法在全流场的计算精度。
1 矩方法与壁面边界条件
1.1 矩方法
经典流体力学中密度ρ、速度ui、温度T对应分子速度分布函数f的前五阶矩。流体力学控制方程可以从Boltzmann方程中推导得到[2-3],即质量、动量、能量守恒方程,如式(1)~(3)所示。
ρt+ρuixi=0
(1)
ρuit+ρuiujxj+σijxj=-pxi+ρai
(2)
ρTt+ρuiTxi+23Rqixi=-23Rpuixi+σijujxi
(3)
其中: txi分别为时间与空间坐标;下标ij为通用张量表示方法; ai为加速度驱动力;压力p服从理想气体状态方程p=ρRTR为理想气体常数。在矩方法中,应力张量σij与热流qi作为自由变量引入拓展流体力学中,此时R13矩方法自由变量为{ρuiTσijqi},其控制方程可由Boltzmann方程推导得到,如式(4)~(5)所示。
σijt+ukσijxk+mijkxk=-pμσij-2puixj+Σij
(4)
qit+ujqixj+12Rijxj=-23pμqi-52pRTxi+Qi
(5)
其中,μ为黏性,
Σij=-4qi5xj-2σkiujxk
(6)
Qi=-7σikR2Txk+σijρpxj+σjkxk-2572qkuixk+qkukxi+qiukxk-16Δxi-mijkujxk-RTσikxk
(7)
mijkΔRij为高阶矩张量,无明确物理意义。为封闭矩方程组,Grad在G13矩方法中令mijk=Δ=Rij=0;Struchtrup等在R13矩方法中为其推导了高阶矩形式[11]。借鉴R13矩封闭思想,Gu等将拓展流体力学方程推广到R26矩方法中,并基于Boltzmann方程推导了mijkRijΔ的控制方程,如式(8)~(10)所示。此时R26矩方法自由变量拓展为{ρuiTσijqimijkRijΔ}。
mijkt+ulmijkxl+ϕijklxl=-32pμmijk+Mijk-3pσij/ρxk
(8)
Rijt+ukRijxk+ψijkxk=-76pμRij+Rij-285pqi/ρxj
(9)
Δt+Δuixi+Ωixi=-32pμΔ-8pqi/ρxi+
(10)
其中:MijkRij为R26矩方程的源项,其表达式参考文献[13]φijklψijkΩi为高阶矩张量,无明确物理意义,Gu等在R26矩方法中为其推导了高阶矩表达式[13],在此不作赘述;下标ijkl、〈、〉均为张量符号[6]
1.2 矩方法壁面边界条件
在矩方法壁面边界条件的推导中,首先假设所有壁面均为漫反射,所有与壁面碰撞的分子均不穿过壁面。在R26矩方程边界条件推导过程中,速度分布函数被截断为Hermite多项式的五阶展开形式,并将垂直于壁面、指向气体的方向向量记作ni,平行于壁面的方向向量记作τi。此时,气体的滑移速度uτ与温度跳跃边界条件可以写为:
(11)
(12)
其中,
(13)
σnnσqτmnnτRnnψnnτΩτφnnnn分别是σijqimijkRijψijkΩiφijkl沿壁面切线与法线的分量。除去下划线的部分之后,矩方法边界条件退化为NSF滑移速度边界与温度跳跃边界条件。下划线部分是高阶矩对边界条件信息的补充,这些信息使得R26矩方程边界条件更加准确。在此基础上,Gu与Emerson首先给出了应力张量σij、热流张量qi以及高阶矩mijkRijψijkΩi的边界条件,其表达式如文献[13]所示。
本文采用的无量纲格式为:v0=2RT0u¯=uv0x¯=xL0t-=v0tL0a¯=L0av02C-i=Civ0T-=TT0f-=v03fρ0ρ-=ρρ0p-=pp0μ-=μμ0。其中v0为参考速度,Ci为分子离散速度的标量形式,μ-x¯等分别代表对应参数的无量纲量。
高阶矩的无量纲格式为:σ-ij=σijρ0v02q-i=qiρ0v03m-ijk=mijkρ0v03ϕ-ijkl=ϕijklρ0v04R-ij=Rijρ0v04ψ-ijk=ψijkρ0v05Ω-i=Ωiρ0v05Δ-=Δρ0v04
2 壁面边界条件的介观表示
矩方法的介观尺度边界条件中,由于分布函数在壁面处被截断为三阶或五阶形式,现有宏观边界条件对中、大Kn情况计算精度较低,且极容易发生发散现象。因此,本文提出在靠近壁面处基于介观Boltzmann模型方程,结合Maxwell壁面边界条件计算宏观参数,并代入矩方程中作为矩方程的边界条件,如图1所示。
1矩方法的介观尺度壁面边界条件迭代方法
Fig.1Iteration scheme of mesoscopic level wall boundary condition for the moment method
图1中,实心圆点区域为矩方法计算区域,基于宏观R13/R26矩方法控制方程求解流场输运特性。空心圆圈层为矩方法/介观方法的耦合重构区,在本层进行宏观参数与介观分布函数数据之间的交互与重构。空心方块区域为介观方法计算区域,靠近壁面处基于介观Maxwell镜面漫反射边界条件求解。
空心圆圈层的数据交换格式包含:
1)边界处介观尺度分布函数f=ftxC)到宏观参数ρTuiσijqimijkRijΔ的数据交换。当得到分布函数f时,流体力学宏观参数可以写为:
ρ=fdCui=1ρCifdCT=23ρc2dCσij=cicjdC-δijρT2qi=12c2cifdCmijk=cicjckfdCRij=cicjc2fdC-7RTσijΔ=c4fdC-15pRTϕijkl=cicjckclfdCψijk=cicjckc2fdC-9RTmijkΩi=c4cifdC-28RTq
(14)
其中,ci=Ci-uic2=Ci-ui2δij为狄拉克函数。
2)宏观参数到介观尺度分布函数的数据交换基于Hermite多项式,如式(15)~(16)所示。
f(3)=fM1+ciqipRTc25RT-1+mijkcicjck6p(RT)2+σijcicj2pRT+Rijcicj4p(RT)2c27RT-1+Δ8pRTc415(RT)2-2c23RT+1
(15)
f(5)=f(3)+fMϕijklcicjckcl24p(RT)3+ψijkcicjck12p(RT)3c29RT-1+ciΩi40p(RT)2c47(RT)2-2c2RT+5
(16)
其中,fM为Maxwell平衡态分布函数[19]
本文提出的方法与经典矩方法边界条件的区别在于:经典矩方法壁面边界条件基于截断的分布函数进行推导,其分布函数损失了计算精度,在推导宏观参数时又进行了二次近似,进一步增大了计算误差;而本文提出的方法直接采用介观尺度Maxwell壁面边界条件,分布函数与壁面处宏观量更为精确,仅在宏观/介观边界重构处损失一次计算精度。同时,壁面处的非平衡效应对流场的影响远大于远场参数的影响,矩方法计算精度也对边界条件更为敏感,因此,本文提出的方法有利于提高矩方法的计算精度。
3 算例与讨论
针对本文提出的介观尺度边界条件,选取无限长平板间Couette流、非平衡热驱流动典型算例进行验证。为保持宏观、介观计算条件一致,本文所有算例共用同一套网格,所有气体介质均选为氩气,Maxwell壁面边界条件采用全漫反射边界条件。算例在物理空间、速度空间网格进行了网格无关性测试,认为介观尺度Boltzmann模型方程的解是标准解。黏性服从Sutherland法则:
μ=μrefTTref1.5Tref+ST+S
(17)
式中,Sutherland法则中的参考温度与参考黏性分别为Tref=273 K、μref=21.25×10-6 Pa·s,Sutherland常数取S=144 K。
3.1 无限长平板间Couette流动
R13/R26矩方法的介观尺度边界条件首先在无限长平板间Couette流动中进行验证。上下平板分别位于y=L0/2与y=-L0/2,参考长度与参考温度分别为L0=10-5 m、T0=273 K,上下平板移动速度分别为uupper=1 m/s、ulower=-1 m/s。在物理空间上沿y方向分布101个非均匀网格点,并在靠近壁面处进行加密。介观Boltzmann模型方程基于离散速度法求解,无限长速度空间被截断为-62RT062RT03,并基于Newton-Cotes积分方法将其离散为64×64×24个非均匀网格点。
Kn取0.01、0.2、0.5、1.0时,分别采用宏观NSF、R13/R26矩方法、改进R13/R26矩方法与介观尺度Boltzmann方法进行迭代求解。考虑Couette流动的对称性,仅给出y∈[0,L0/2]区域内速度分布曲线,对比计算结果如图2所示。
图2(a)可以看出:在滑移流区域(Kn=0.01),所有宏观方程与介观Boltzmann模型方程均能得到一致计算结果。随着Kn增大到过渡流时,如图2(b)、(c)所示,宏观方程在壁面处计算精度均大幅降低,在宏观方法计算精度对比中:R26矩方法>R13矩方法>NSF方法。所提出的改进R13与改进R26矩方法相较原矩方法最大分别提高了20.66%与12.51%的计算精度。当Kn=1.0时,R26矩方法无法获得稳态解,所提出的改进R13矩方法相较于原矩方法最大提高了29.23%的计算精度,并且将R26矩方法Kn的适用范围拓展到1.0,如图2(d)所示。
2无限长平板间Couette流动中速度ux沿y轴分布
Fig.2ux-velocity along the y-axis in the Couette flow between two infinite plates
3.2 非平衡热驱流动
R13/R26矩方法介观尺度边界条件在非平衡热驱流动中进行进一步验证。方腔上下壁面分别位于y=L0/2与y=-L0/2处,左右壁面位于x=-2.5L0x=2.5L0处。左右壁面温度设置为Tw=233 K、Te=313 K,上下壁面温度按照线性增长从233 K增加到313 K,如图3所示。
3非平衡热驱流动示意图
Fig.3Schematic of thermally induced non-equilibrium flows
在本算例中参考长度与参考温度分别为L0=10-5 m、T0=273 K,初始压力p0通过Kn反算得到。物理空间网格划分为101×51个非均匀网格点,无限长速度空间被截断为-62RT062RT03,并基于Newton-Cotes积分方法将其离散为64×64×24个非均匀网格点。非平衡热驱流动中x轴中心线处速度ux分布、温度场与流线分布如图4图5所示。
4Kn=0.1时,非平衡热驱流动中x轴中心线处速度ux沿y轴分布
Fig.4ux-velocity along y axis of the x axis centerline for thermally induced non-equilibrium flows at Kn=0.1
5Kn=0.1时,非平衡热驱流动中温度场与流线分布
Fig.5Temperature field and streamline distribution for thermally induced non-equilibrium flows at Kn=0.1
图4图5可以看出:非平衡热驱流动中稀薄效应对宏观方程计算精度影响极大,所有宏观方程均无法捕获准确的边界处滑移速度。改进R13与改进R26矩方法相较原方法最大分别提高了59.84%与50.18%的计算精度,尤其改进R26矩方法可以准确预测壁面滑移速度,捕获速度分布独特的双峰现象,其计算结果可以基本恢复介观尺度Boltzmann方程结果。
4 结论
本文针对宏观尺度R13/R26矩方法提出了介观尺度边界条件。仿真结果表明:所提出的边界条件能够有效提高矩方法的计算精度,尤其当0.5<Kn<1.0时,原矩方法边界条件计算精度大幅度降低,此时R26矩方法极难收敛。基于介观尺度边界条件的R13/R26矩方法最大提高了59.84%的计算精度,并且拓展了R26矩方法对Kn的适用范围,改进的R26矩方法能够准确地捕获非平衡热驱流动中速度独特的双峰分布现象。
致谢
英国达斯伯里实验室(Daresbury Laboratory)顾晓军教授在关于R13/R26矩方法、耦合方法的研究中提供了帮助和指导,谨致谢意!
1矩方法的介观尺度壁面边界条件迭代方法
Fig.1Iteration scheme of mesoscopic level wall boundary condition for the moment method
2无限长平板间Couette流动中速度ux沿y轴分布
Fig.2ux-velocity along the y-axis in the Couette flow between two infinite plates
3非平衡热驱流动示意图
Fig.3Schematic of thermally induced non-equilibrium flows
4Kn=0.1时,非平衡热驱流动中x轴中心线处速度ux沿y轴分布
Fig.4ux-velocity along y axis of the x axis centerline for thermally induced non-equilibrium flows at Kn=0.1
5Kn=0.1时,非平衡热驱流动中温度场与流线分布
Fig.5Temperature field and streamline distribution for thermally induced non-equilibrium flows at Kn=0.1
YANG W Q, GU X J, WU L,et al. A hybrid approach to couple the discrete velocity method and method of moments for rarefied gas flows[J]. Journal of Computational Physics,2020,410:109397.
LIU W, ZHANG Y B, ZENG J N,et al. Further acceleration of multiscale simulation of rarefied gas flow via a generalized boundary treatment[J]. Journal of Computational Physics,2024,503:112830.
STRUCHTRUP H, FREZZOTTI A. Twenty-six moment equations for the Enskog-Vlasov equation[J]. Journal of Fluid Mechanics,2022,940: A40.
CHAPMAN S, COWLING T G, PARK D. The mathematical theory of non-uniform gases[J]. The Mathematical Gazette,1954,38(323):63-64.
BOBYLEV A V. On some properties of generalized Burnett equations of hydrodynamics[J]. Journal of Statistical Physics,2022,188:6.
GRAD H. Note on N-dimensional Hermite polynomials[J]. Communications on Pure and Applied Mathematics,1949,2(4):325-330.
STRUCHTRUP H, FREZZOTTI A. Grad′s 13 moments approximation for Enskog-Vlasov equation[C]//Proceedings of the 31st International Symposium on Rarefied Gas Dynamics,2019.
GRAD H. Asymptotic theory of the Boltzmann equation[J]. Physics of Fluids,1963,6(2):147-181.
KARLIN I V, GORBAN A N. Hydrodynamics from Grad′s equations:what can we learn from exact solutions?[J]. Annalen Der Physik,2002,514(10/11):783-833.
REITEBUCH D, WEISS W. Application of high moment theory to the plane Couette flow[J]. Continuum Mechanics and Thermodynamics,1999,11:217-225.
STRUCHTRUP H,ÖTTINGER H C. Thermodynamically admissible 13-moment equations[J]. Physics of Fluids,2022,34(1):017105.
GU X J, BARBER R W, EMERSON D R. How far can 13 moments go in modeling microscale gas phenomena?[J]. Nanoscale and Microscale Thermophysical Engineering,2007,11(1/2):85-97.
GU X J, EMERSON D R. A high-order moment approach for capturing non-equilibrium phenomena in the transition regime[J]. Journal of Fluid Mechanics,2009,636:177-216.
YANG W Q, TANG S, YANG H. Analysis of the moment method and the discrete velocity method in modeling non-equilibrium rarefied gas flows:a comparative study[J]. Applied Sciences,2019,9(13):2733.
YANG W Q, GU X J, EMERSON D R,et al. Modelling thermally induced non-equilibrium gas flows by coupling kinetic and extended thermodynamic methods[J]. Entropy,2019,21(8):816.
SHENG Q, TANG G H, GU X J,et al. Simulation of thermal transpiration flow using a high-order moment method[J]. International Journal of Modern Physics C,2014,25(11):1450061.
WU L, HO M T, GERMANOU L,et al. On the apparent permeability of porous media in rarefied gas flows[J]. Journal of Fluid Mechanics,2017,822:398-417.
GU X J, BARBER R W, JOHN B,et al. Non-equilibrium effects on flow past a circular cylinder in the slip and early transition regime[J]. Journal of Fluid Mechanics,2019,860:654-681.
MAXWELL J C.Ⅶ. On stresses in rarified gases arising from inequalities of temperature[J]. Philosophical Transactions of the Royal Society of London,1879,170(170):231-256.
LIU I S, RINCON M A. A boundary value problem in extended thermodynamics:one-dimensional steady flows with heat conduction[J]. Continuum Mechanics and Thermodynamics,2004,16(1):109-124.
GU X J, EMERSON D R. A computational strategy for the regularized 13 moment equations with enhanced wall-boundary conditions[J]. Journal of Computational Physics,2007,225(1):263-283.
TORRILHON M, STRUCHTRUP H. Boundary conditions for regularized 13-moment-equations for micro-channel-flows[J]. Journal of Computational Physics,2008,227(3):1982-2011.
YANG W Q, GU X J, EMERSON D,et al. Comparative study of the discrete velocity and the moment method for rarefied gas flows[C]//Proceedings of the 31st International Symposium on Rarefied Gas Dynamics,2018:23-27.