针对格子玻尔兹曼方法数值模拟中曲面边界的处理方法

1.本发明涉及流体力学技术领域，具体为一种计算流体力学中相关数值模拟针对曲面物理边界的处理方法，用于翼型、圆柱绕流及发动机内流等的数值模拟。


背景技术：

2.格子玻尔兹曼方法是一种介观的计算流体力学算法，基于流体分子在网格格点的碰撞迁移理论而求解分布函数，从而获得网格格点上的宏观物理量，比如速度，密度，压力等。在实际数值模拟当中，对于翼型、圆柱绕流及发动机内流等复杂外形通常具备的曲面边界格点(如图一所示)，没有办法参与整个流场的演化过程(碰撞迁移)。为此，通常使用插值格式对这些曲面边界点人工赋值(曲面边界处理)。现有的较为流行的处理格式有三种，分别是mei格式
1.,chen格式
2.以及guo格式
3.。其中mei格式和guo格式把分布函数分为平衡态和非平衡态两部分，然后通过创造虚拟点来插值求解物面点的分布函数。其缺点，第一不符合物理实际，因为虚拟点的位置实际上是物体内部，不存在流体分子的迁移碰撞；第二虚拟点的物理量(速度和密度)插值格式比较简单，没有包含足够信息。chen格式通过对流体点的线性内插值直接求解物面点的分布函数，插值格式也比较简单没有包含物面点附近足够多的信息且没有将分布函数按非平衡态和平衡态来处理，在物理层面不够精确。以上三种曲面边界格式在处理一般数值模拟问题上表现优异，但是对于追求高精度高效率数值模拟方面，尤其是在高性能翼型设计、新一代发动机设计等方面仍有巨大优化或发展空间。为了满足越来越高要求的翼型绕流数值模拟以及发动机内腔内流数值模拟等实际应用需求，提高针对具体工程应用背景的数值模拟精度和效率，需要在数值模拟过程中，开发高精度高效率的计算方法。


技术实现要素：

3.为解决现有技术存在的问题，满足高精度高效率的数值模拟需求同时符合物理实际情况，本发明基于实际数值模拟计算需求综合以上三种曲面处理格式，提出一种新的适用于lbm数值模拟的曲面边界优化处理方法，能够用于开发满足复杂流场高精度高效率数值模拟的计算平台。基于该方法，进而提出了相应的翼型绕流数值模拟以及发动机内腔内流数值模拟方法，进行高性能翼型设计、新一代发动机设计等实际应用领域，提高了翼型绕流数值模拟以及发动机内腔内流数值模拟的精度和效率，推动本学科的发展和创新。
4.本发明的技术方案为：
5.一种针对格子玻尔兹曼方法数值模拟中曲面边界的处理方法，对物面边界中的点f通过以下过程求解点f沿α方向的分布函数m
α
：
6.步骤1：将点f沿α方向的分布函数分为平衡态分布函数和非平衡态分布函数
7.8.其中平衡态分布函数通过平衡态分布函数function(ρ,u)计算公式得到；
9.步骤2：采用以下步骤求解点f的非平衡态分布函数
10.步骤2.1：使用线性内插值来综合点ff和w的信息
[0011][0012]
其中为嵌入深度；
[0013]
步骤2.2：利用线性外插值来综合点fff和ff的信息
[0014][0015]
步骤2.3：平均线性内外插值的结果得到：
[0016][0017]
其中：采用以下过程计算w的信息
[0018]
点w的分布函数分为非平衡态和平衡态两部分
[0019][0020]
而由物面无滑移条件限制确定得到其中w点的平衡态分布函数通过平衡态分布函数function(ρ,u)计算公式得到；
[0021]
同时使用线性内外插值求解点w的分布函数
[0022][0023][0024]
从而得到其中b点沿-α方向的分布函数在经过碰撞迁移之后由f点沿-α方向的分布函数迁移而来，f点沿-α方向的分布函数在经过碰撞迁移之后由ff点沿-α方向的分布函数迁移而来，ff点沿-α方向的分布函数在经过碰撞迁移之后由fff点沿-α方向的分布函数迁移而来；
[0025]
从而得到：
[0026][0027]
进而得到
[0028][0029]
进一步的，步骤1中，计算平衡态分布函数maeq时，点f的密度和速度用点f周围的点来插值得到：
[0030][0031]
进一步的，步骤3中，运用点ff和f的信息来插值点w的速度和密度，进而通过平衡态分布函数计算得到w点的平衡态分布函数
[0032]
有益效果
[0033]
本发明提出的针对格子玻尔兹曼方法数值模拟中曲面边界的处理方法能够满足lbm复杂流场，复杂外形数值模拟中遇到的曲面边界处理问题，同时比以上三种现有曲面边界处理格式更为精确，由于每一步计算都更为精确，所以收敛速度也较快，实现了数值模拟的高精度和高效率，在翼型设计、发动机设计等实际应用中，能够取得较好应用效果。
[0034]
本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。
附图说明
[0035]
本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：
[0036]
图1展示了直角网格计算域中的曲面物面边界；
[0037]
图2表述了针对格子玻尔兹曼方法的曲面物面边界处理格式的基本构造；
[0038]
图3描述了圆柱绕流的定常状态(re＝40)流线图；
[0039]
图4描述了圆柱绕流的非定常周期性状态(re＝3900)流线图；
[0040]
图5描述了naca0012翼型的数值模拟结果(re＝3000,α＝0
°
)；(1)上表面压力系数，(2)流线图和压力云图；
[0041]
图6描述了naca0012翼型的数值模拟结果(re＝100000,α＝0
°
)；(1)升力系数，(2)阻力系数，(3)流线图，(4)压力云图；
[0042]
图7描述了sd7003翼型(re＝60000,α＝14
°
)和u型腔(re＝1000)的数值模拟结果；(1)sd7003翼型，re＝60000,α＝14
°
流线图，(2)u型方腔内流流线图，re＝1000。
具体实施方式
[0043]
本发明针对基于lattice boltzmann method(lbm)算法的复杂外形数值模拟中遇到的曲面物面边界处理问题，提出了一种新的曲面边界处理方法，在实际工程计算中不仅可以处理翼型绕流等外流流场，还适用于拥有曲面物理边界的腔体内流流场数值模拟。
[0044]
下面分别以圆柱绕流数值模拟、naca0012翼型扰流流场数值模拟、sd7003翼型扰流流场数值模拟、u型腔体内流流场数值模拟为例，说明新的曲面边界处理方法。
[0045]
实施例1，圆柱绕流数值模拟：
[0046]
在进行圆柱绕流数值模拟时，建立圆柱三维模型、网格划分等cfd常规操作后，对圆柱物面边界中的点通过以下过程求解分布函数：
[0047]
如图二所示，以点f沿α方向的分布函数m
α
求解为例，点fff,ff和f为流体点，每一次演化过后，点fff和ff的信息都得了更新。而点f沿α方向的分布函数m
α
未知，因为它是从
虚拟点b迁移而来，而b点为虚拟点，没有分布函数。因此我们只能通过点fff,ff以及w的信息来构造f点沿α方向的分布函数m
α
。
[0048]
首先把点f沿α方向的分布函数分为平衡态分布函数和非平衡态分布函数来处理。
[0049][0050]
其中平衡态分布函数可以通过平衡态分布函数function(ρ,u)计算公式得到，所以我们只需要知道点f的密度和速度就可以，在本方法中点f的密度和速度可以用点f周围的点来插值
[0051][0052]
这样我们就可得到点f的平衡态分布函数接下来我们来求解点f的非平衡态分布函数
[0053]
首先我们使用线性内插值综合点ff和w的信息
[0054][0055]
其中为嵌入深度，然后再利用线性外插值综合点fff和ff的信息
[0056][0057]
然后平均线性内外插值的结果得到
[0058][0059]
每一次演化之后公式(4)中除了这一项其他都是知道的，因此接下来我们需要求解
[0060]
首先我们知道，点w的分布函数可以分为非平衡态和平衡态两部分
[0061][0062]
或
[0063][0064]
同时我们知道由于物面无滑移条件限制
[0065][0066]
因此公式(6)变为
[0067][0068]
由于w点的平衡态分布函数可以通过平衡态分布函数function(ρ,u)计算公式得到，我们只需运用点ff和f的信息来插值点w的速度和密度，如果考虑静止物面则物面速度为0，如果物面不静止，只需要对应具体情况就可以确定点w的速度。然后，我们再一次同时
使用线性内外插值求解点w的分布函数
[0069][0070][0071]
综合公式(9)和(10)
[0072][0073]
其中b点沿-α方向的分布函数在经过碰撞迁移之后由f点沿-α方向的分布函数迁移而来,同理f点沿-α方向的分布函数在经过碰撞迁移之后由ff点沿-α方向的分布函数迁移而来,ff点沿-α方向的分布函数在经过碰撞迁移之后由fff点沿-α方向的分布函数迁移而来。将公式(11)带入公式(8)
[0074][0075]
把公式(12)代入公式(4)我们得到
[0076][0077]
至此点f沿α方向的分布函数m
α
求解完毕，综合使用了点fff,ff,f以及w的信息构造该点非平衡态分布函数，同时使用红圈标点的综合信息构造平衡态分布函数。
[0078]
相比较chen格式，本发明在构造分布函数的时候区别对待非平衡态和平衡态分布函数。并且新格式同时采用了线性内插值和外插值综合考虑周边点的信息和影响。相比较mei格式和guo格式，在构造平衡态分布选取了更为合理的周围众多点作为目的点速度和密度的插值信息来源，没有强行构造虚拟点(与物理实际不符)的速度和密度。本格式中b点沿-α方向的分布函数只是一个用来做数据交换(赋值)的概念。同时计算非平衡态分布函数时mei格式和guo格式只使用了点f的信息作为插值信息来源，而本格式选用内外线性插值显得更为合理且同时考虑了周边多点的综合信息。
[0079]
[0080][0081]
表1.使用本发明的圆柱绕流数值计算结果与其他数据的对比(re＝40)
[0082][0083][0084]
表1给出了使用本发明提出的新型曲面边界处理格式后的圆柱绕流在re＝40时定常流动的流动数据与其他结果的对比。由表可见，新型曲面边界处理格式的计算精度优于
其他三种流行的曲面边界处理格式。由于新型格式区别于其他格式，分为前期预处理和后期构造两部分，所以相较于其他格式每一个迭代步所需时间较长，但是由于新格式本身具有的高精度性在每一个迭代步之后更接近理论解所以可以加速收敛，从而整体上所耗费的时间与其他三种格式基本相同。
[0085]
图3展示了圆柱绕流re＝40时定常流动的流线图，可以明显观察到两个稳定的附着涡。
[0086]
表2.使用新格式的圆柱绕流hopf流动分岔点与其他结果的对比
[0087][0088]
表2给出了不同文献及三种流行格式的对圆柱绕流hopf流动分岔点的预测与新格式计算结果的对比。可以看到，新格式的计算结果与其他文献的计算结果非常吻合，说明新格式的计算精度较高。
[0089]
表3.使用新格式的圆柱绕流计算结果与其他文献的对比(re＝3900).
[0090][0091][0092]
图4和表3展示了使用新格式的圆柱绕流的非定常周期性流动状态的计算结果。图4(a)描述了升力系数随时间的变化曲线，t0,t1,t2,t3,t4表示了一个完整周期，对应每个时间节点图4(b)、(c)、(d)、(e)、(f)依次展示了圆柱绕流re＝3900时的瞬时流线图。
[0093]
实施例2，naca0012翼型的扰流流场数值模拟，包括re＝3000,攻角为0
°
以及re＝100000,攻角为0
°
两个状态；
[0094]
实施例3，sd7003翼型的扰流流场数值模拟，状态为re＝60000,攻角为14
°
；
[0095]
实施例4，u型腔体内流流场数值模拟，状态为re＝1000；
[0096]
上述实施例2,3,4中，对naca0012翼型、sd7003翼型以及u型腔体物面边界中的点采用与实施例1中相同的方法求解分布函数，从而得到的数值模拟结果如图5～图7所示。
[0097]
图5展示了naca0012翼型的扰流流场数值模拟，re＝3000,攻角为0
°
时的上表面压力系数、流线图和压力云图。图6描述了naca0012翼型的扰流流场数值模拟，re＝100000,攻角为0
°
时的升阻力系数、流线图和压力云图。图7展示了sd7003翼型的扰流流场数值模拟，re＝60000,攻角为14
°
时的流线图和同时也给出了u型腔体内流流场的流线图(re＝1000)。可以看出，本发明提出的曲面边界格式适用于翼型绕流外流场和u型腔体内流流场的数值模拟。
[0098]
尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。
[0099]
本发明中引用的参考文献为：
[0100]
[1]mei,r.w.,yu,d.z.and shyy,w.an accurate curved boundary treatment in the lattice boltzmann method.journal of computational physics,155(2),307-330.(1999)doi:10.1006/jcph.1999.6334
[0101]
[2]chen,y.,xia,z.h.and cai,q.d.lattice boltzmann method with tree-structured mesh and treatment of curved boundaries.chinese journal of computational physics,27(1):23-30.(2010)
[0102]
[3]guo,z.l.,zheng,c.g.and shi,b.c.an extrapolation method for boundary conditions in lattice boltzmann method.physics of fluids,14(6),2007-2010.(2002)doi:10.1063/1.1471914
[0103]
[4]f.bouchon,t.dubois,and n.james,“a second-order cut-cell method for the numerical simulation of 2d flow past obstacles,”computers&fluids,65,80-91.(2012)
[0104]
[5]k.taira,and t.colonius,“the immersed boundary method:a projection approach,”journal of computational physics,225(2),2118-2137.(2007)
[0105]
[6]m.n.linnick,and h.f.fasel,“a high-order immersed interface method for simulating unsteady incompressible flows on irregular domains,”journal of computational physics,204(1),157-192.(2005)
[0106]
[7]z.l.guo,b.c.shi,and n.c.wang,“a nonuniform lattice boltzmann method based on domain decomposition,”chinese journal of computational physics,18(2),181-184.(2001)
[0107]
[8]s.behara,and s.mittal,“flow past a circular cylinder at low reynolds number:oblique vortex shedding,”physics of fluids,22(5),article number:054102.(2010)(critical value re＝47,re＝100,st,cl and cd)
[0108]
[9]t.k.prasanth,and s.mittal,“vortex-induced vibrations of a circular cylinder at low reynolds numbers,”journal of fluid mechanics,594,463-491.(2008)(critical value re＝47,re＝100,st,cl and cd)
[0109]
[10]m.morzynski,and f.thiele,“solution of the eigenvalue problems resulting from global non-parallel flow stability analysis,”journal of fluid mechanics,169(1-2),161-176.(1999)(critical value re＝47.0)
[0110]
[11]c.norberg,“an experimental investigation of the flow around a circular cylinder:influence of aspect ratio,”journal of fluid mechanics,258,287-316.(1994)(critical value re＝47.4)
[0111]
[12]j.h.williamson,“oblique and parallel modes of vortex shedding in the wake of a circular at low reynolds numbers,”journal of fluid mechanics,206,579-627.(1989)(critical value re＝47.9)
[0112]
[13]kravchenko,g.and moin,p.numerical studies of flow over a circular cylinder at red＝3900.phys.fluids,12(2),403.(2000)
[0113]
[14]you,d.and moin,p.a dynamic global-coefficient subgrid-scale eddy-viscosity model for large eddy simulation in complex geometries.physics of fluids,19(6),065110.(2007)
[0114]
[15]rajani,b.n.,kandasamy,a.and majumdar,s.les of flow past circular cylinder at re＝3900.journal of applied fluid mechanics,9(3),1421-1435.(2016)