基于改进型FVM-LBFS方法的流场计算方法与流程

基于改进型fvm
‑
lbfs方法的流场计算方法
技术领域
1.本发明的一种基于改进型fvm
‑
lbfs方法的流场计算新方法，属飞行器气动计算领域。


背景技术：

2.20世纪80年代后期，随着计算机技术的飞速发展，计算机运行速度的越来越快，从介观层次开展流体力学的数值模拟方法开始逐步兴起。其中典型代表之一的格子boltzman方法(lattice boltzmann method，lbm)自1988年提出以来得到研究人员的广泛关注，在近20余年的时间内得到了飞速的发展。与基于连续性假设的n
‑
s方程不同的是，格子boltzman方法是一种介于流体微观分子动力学计算方法和宏观连续方程计算方法之间的介观尺度流体计算力学方法，因此与传统cfd方法相比，格子boltzman方法不受连续介质假设的限制，可实现跨尺度流动模拟。
3.格子boltzman方程是对boltzman方程在速度空间、物理空间和时间上进行离散得到的，其描述了粒子以一定速度在下一个时间步内运动到相邻点的过程。格子boltzmann方法的一个优势在于：它将传统方法求解非线性偏微分方程组(如euler方程组、n
‑
s方程组等)的困难，在模型建立时一次性完成，因而在数值模拟中只需要处理格子boltzmann方程简单线性方程或方程组，与宏观方法相比，方程个数和形式得到大大简化。标准的格子boltzmann方法的求解过程只包括简单的线性迁移过程和碰撞过程两个部分，仅需代数运算，无需求解压力泊松方程，因此程序实现起来也非常简单。鉴于格子boltzmann方法具有许多独特的优势，其已经发展成为一种新的模拟复杂流动和物理现象的计算工具。因此，在许多研究领域，例如可压缩/不可压缩流动、高超声速流动、多相流动、多孔介质流动、化学反应流动、稀薄气体流动等，基于格子boltzman模型的数值方法都得到了很好的应用，而且也揭示了许多复杂流动问题的物理作用机理。
4.格子boltzmann方法的实现主要基于格子boltzmann模型(lattice boltzmann model)和相应的数值方法两个部分，而格子boltzmann模型中的离散速度模型和相应的平衡分布函数形式是关键，当前国内外研究学者在此方面做出了大量的研究工作。传统的格子boltzmann模型大多数是针对不可压缩流动建立的，例如1992年qian等提出的ddqm系列基本模型(d表示空间维度，m表示离散速度个数)。为了拓展lbm方法在可压缩流动中的应用，可压缩格子boltzmann模型也不断被提出，例如多速度模型、比热容比可调模型、耦合双分布模型等。基于这些模型，lbm方法实现了对ma>0.3流动的数值模拟，但很多格子boltzmann模型目前也只适用于低马赫数可压缩流动(由于热力过程而导致的气动密度明显变化的流动)。这是因为这些模型大都基于maxwellian分布函数的泰勒级数展开建立的，截断误差与马赫数直接相关，导致不适用于高马赫数流动。为了克服这一问题，qu等人利用圆函数代替麦克斯韦分布函数，建立了一系列一维和二维模型，如d1q5l2模型、d2q13l2模型，并将这些模型成功地应用于激波管、双马赫数反射等强激波超声速流场的数值模拟。dellar也推导出了一些二维的可压缩格子boltzmann模型，如不分裂七速度模型、4 3分裂
速度模型。然而这些模型中含有许多自由参数，格子速度需要进行人为给定，这样导致在复杂高马赫数流动数值模拟中应用的通用性较差。鉴于此，shu与其研究团队通过理论推导，给出了一个基于矩守恒形式的非自由参数格子boltzmann模型推导平台，通过该平台得到了一系列求解可压缩无粘流动的非自由参数格子boltzmann模型，如d1q3模型、d1q4模型、d2q8模型、d2q12模型等。利用这些模型可以局部精确地恢复宏观方程，在当地单元界面建立一维可压缩无粘流动riemann通量求解器。
5.由于格子均匀性的要求，传统的lbm方法对计算网格的均匀性要求非常高，因此不适用于复杂外形的计算。而且对于高雷诺数流动，由于迁移时间步长的限制，计算网格会非常密集导致计算效率大幅度降低。因此，当采用格子boltzmann方法对可压缩流动进行数值模拟时，大多数采用基于离散速度boltzmann方程(discrete velocity boltzmann equation，dvbe)的方法计算来及时更新粒子分布函数。qu，q.li，ubertini以及stiebler等很多学者在dvbe方法上做了大量的研究工作。但需要注意的是，在dvbe方法中，粒子分布函数的个数远多于守恒变量的个数，且分布函数形式非常复杂，这样导致dvbe方法的直接求解过程繁琐，相比传统cfd方法的计算效率差很多。为了开发一种更有效的具有lbm性质的可压缩流动求解器，shu和他的研究团队提出了一种基于有限体积方法的格子boltzmann通量求解器(lbfs)。众所周知，有限体积方法(fvm)在cfd数值模拟中有其独特优势：fvm从积分形式的守恒方程出发，推导过程物理概念清晰，在每个控制体内保持物理量的守恒特性；且对于计算网格(结构/非结构网格)都具有很好的灵活性，适用于复杂几何外形的计算。因此，在lbfs方法中，避免直接求解dvbe方程的复杂过程，利用传统的fvm的优势对宏观控制方程进行离散，采用一维可压缩格子boltzmann模型的局部解重构单元界面的无粘通量。数值计算结果表明，lbfs方法在强激波可压缩流动模拟中具有良好的正值性和高效性。在现有的lbfs方法中，为了捕捉强激波等强间断物理流动，利用单元界面分布函数的非平衡部分来引入数值耗散来保证计算的稳定性。然而，这样直接引入的数值耗散量太大，严重干扰边界层等光滑区域的真实解，特别是对于精确预测高超声速流场的热流密度而言是不可接受的。因此，如何保证lbfs方法在高超声速流动中同时稳定捕捉强激波流动和准确预测气动热，还需要进一步完善，这也是本发明研究的主要内容。
6.本发明针对现有基于格子boltzmann模型数值方法模拟高马赫流动的不足之处，将提出了一种改进型fvm
‑
lbfs方法用于高超声速流动的数值模拟并给出详细推导过程。该方法采用有限体积法对宏观n
‑
s控制方程进行空间离散，利用可压缩lb模型来构造当地一维riemann通量求解器，计算过程简单高效。改进型fvm
‑
lbfs方法将采用新提出的压力和温度同时控制的改进型开关控制函数实现了对现有lbfs方法中无粘通量数值粘性的精确控制，并引入网格单元长细比修正系数拓展非均匀网格适用范围。最后，选取一系列高速数值算例对所提出计算方法的可靠性与准确性进行验证与分析。


技术实现要素：

7.本发明主要采用数值计算手段，围绕高超声速气动加热问题，发展了一种基于改进型fvm
‑
lbfs方法的流场计算新方法。
8.针对现有基于格子boltzmann模型数值方法模拟高雷诺数流动难的问题，提出了一种改进型fvm
‑
lbfs方法用于高超声速流动的数值模拟。该方法将传统有限体积方法与
lbm方法相结合，采用有限体积法对宏观n
‑
s控制方程进行空间离散，利用一维可压缩格子boltzmann模型的局部解重构单元界面的无粘通量，计算过程简单高效。改进型fvm
‑
lbfs方法将采用新提出的压力和温度同时控制的改进型开关控制函数实现了对现有lbfs方法中无粘通量数值粘性的精确控制，并引入网格单元长细比修正系数拓展非均匀网格适用范围。通过一系列的数值模拟算例验证与分析，改进型fvm
‑
lbfs方法在高超声速流动数值模拟中可以同时稳定捕捉复杂强激波流动和准确预测气动热参数，很好地拓展了基于格子boltzmann模型数值方法在高超声速流动数值计算领域的应用范围。
9.本发明的技术方案：
10.基于改进型fvm
‑
lbfs方法的流场计算方法，步骤如下：
11.采用平均雷诺n
‑
s(rans)方程作为控制方程，其积分形式(无源项)如下，
[0012][0013]
式中：w为守恒量；f
c
为无粘通量；f
v
为粘性通量；w、f
c
、f
v
的定义为：
[0014][0015][0016]
其中ρ、p和t分别表示密度、压力、温度，且满足理想气体状态方程p＝ρrt，r是理想气体常数。u、v、w分别为控制体内流体三个方向的速度，n
x
、n
y
、n
z
分别为控制体单元面ds的外法向单位向量分量，v＝n
x
u n
y
v n
z
w为面法向速度，τ
ij
为摩擦应力张量分量，k为热传导系数，e＝p/ρ(γ
‑
1)为控制体内流体的能量，h为总焓，γ为比热比；h＝e p/ρ。
[0017]
从上可知，从方程(1)可以看出对其求解的关键步骤在于对无粘通量f
c
和粘性通量f
v
的计算。本发明将采用基于非自由参数d1q4模型的lbfs方法计算无粘通量f
c
，考虑到粘性通量的物理性质，将采用中心差分方法计算粘性通量f
v
。下面将详细介绍和推导lbfs方法求解无粘通量f
c
的计算过程。
[0018]
lbfs方法是基于格子boltzmann模型建立的，为了准确模拟高超声速流动，合理选取可压缩lb模型对lbfs方法是非常重要的。本发明将选取非自由参数的d1q4模型用于本发明的lbfs方法，下面简要阐述d1q4模型的推导过程：
[0019]
以使用maxwell分布函数来获得还原euler方程或n
‑
s方程所需的各阶矩关系，于恢复一维euler方程，可以将动量守恒关系式简化成如下求和形式，
[0020][0021]
[0022][0023][0024]
式中，为第i方向的粒子平衡分布函数，ξ
i
为相应方向的粒子速度，ρ和u分别为一维宏观流动的密度和速度，式中表示运动粒子的特定速度。在进行格子boltzmann模型的推导时，首先需要确定lb模型的粒子速度分布形态。可以看出，可压缩d1q4模型由两个粒子速度(d1，d2)和四个平衡分布函数)和四个平衡分布函数组成，因此总共存在6个待求的未知参数。由于方程(1)只含有4个独立方程，想要求解这6个未知量，还需要再引入如下两个独立方程
[0025][0026][0027]
方程(3)和(4)为两个更高阶的动量守恒关系式，这两个附加的高阶矩关系方程可以对可压缩d1q4模型进行进一步的约束，使得其能更加真实地反映物理流动属性，特别是对于高超声速流动中的强激波问题的计算具有更好的正值性。
[0028]
由此，联立方程(2)
‑
(4)，可以得到如下方程组形式，
[0029][0030][0031][0032][0033][0034][0035]
求解方程组(5)可以得出可压缩d1q4模型中4个平衡分布函数和2个离散粒子速度表达式分别为
[0036][0037][0038][0039][0040]
[0041][0042]
fvm
‑
lbfs构建n
‑
s方程与非自由参数d1q4模型之间的联系，根据恢复euler方程的物理守恒关系式(2)，关于法向方向速度的守恒变量和无粘通量f
c
可以通过如下式子计算得到，
[0043][0044][0045]
式中，ξ
i
表示第i方向的离散粒子相速度，在一维非自由参数d1q4模型中，有ξ1＝d1，ξ2＝
‑
d1，ξ3＝d2，ξ4＝
‑
d2。v
n
为界面法向速度，e为单位质量内能；则表示矩关系向量，其表达式为
[0046][0047]
另外，e
p
为粒子势能，f
i
(r,t)表示t时刻，在空间位置r处第i方向的粒子分布函数，通过联立以上方程可以分别得到整个守恒变量w、无粘通量f
c
与关于法向方向速度的守恒变量和无粘通量之间的联系，如下
[0048][0049]
其中，分别为式子(8)中对应的第1，2，3项，u
τ
为切向方向速度，u
τx
、u
τy
、u
τz
分别为u
τ
在三个方向的分量，、分别为(9)中对应的第1、2、3项；
[0050]
在本发明改进型fvm
‑
lbfs方法中，核心部分在于无粘通量的推导与求解，本发明从chapman
‑
enskog展开分析出发，在lbfs方法中同时保留粒子分布函数在单元界面上的平衡部分和非平衡部分，本发明lbfs方法仅用于计算单元界面的无粘通量，分布函数在单元界面上的非平衡部分对通量的贡献可以看作是人工数值粘性，数值粘性的控制需要根据流场的物理特征进行准确控制，下面介绍无粘通量计算的推导过程:
[0051]
在本发明lbfs中，格子boltzmann模型用于在单元界面处建立局部riemann通量求解器来计算无粘通量。可以得到单元界面处关于法向速度的无粘通量表达式如下：
[0052][0053]
式中，符号“*”表示数值解，下标“j 1/2”表示单元界面位置，上角标i表示单元界面处平衡分布函数对无粘通量产生的贡献，上角标ii表示单元界面周围点的平衡分布函数
对无粘通量产生的贡献。进一步地，将切向方向速度对无粘通量的贡献合并到式(12)中，单元界面处无粘通量的完整表达式可以写成如下形式
[0054][0055]
δt为粒子迁移时间步长，本发明将分布函数非平衡部分的贡献视为数值粘性，τ0可作为数值粘性的权重系数；
[0056]
单元界面处的密度、法向动量和法向速度对对能量产生的贡献，如下
[0057][0058]
根据相容性条件可知，守恒变量的计算不依赖于粒子分布函数的非平衡部分。因此，方程(14)中分布函数非平衡部分可以写成，
[0059][0060]
可以得到，
[0061][0062]
从方程(16)可以看出，单元界面上的守恒变量可以通过平衡分布函数g
i
(
‑
ξ
i
δt,t
‑
δt)计算得出。进一步地，，我们可以得到单元界面上的守恒变量表达式，
[0063][0064]
式中，分别表示单元界面左右的平衡分布函数；
[0065]
由于在方程(17)中应用了一维非自由参数lb模型，因此仅能获得单元界面上与法向速度相关的守恒变量。为了计算单元界面上切向速度的贡献，本发明采用了一种可行的近似方法，如下所示，
[0066][0067]
式中，分别为单元界面左侧、单元界面右侧和单元界面上的切向速度。其中，和可以分别由重构后的单元界面左右两侧的总速度与法向速度的差值计算得到。通过联立方程，我们可以直接计算得到单元界面上的守恒变量ρ
*
，和p
*
。进一步地，将计算得到的守恒变量ρ
*
，和p
*
代入方程(6)
‑
(7)就可以计算得到单元界面上的平衡分布函数g
i
(0,t)。随着平衡分布函数g
i
(
‑
ξ
i
δt,t
‑
δt)和g
i
(0,t)被确定，单元界面上相应的无粘通量就可以显式地计算出来。
[0068]
为了准确控制整个计算域内单元界面分布函数的非平衡部分产生的数值粘性的大小，需要构造一个全局开关函数来计算数值粘性权重系数τ0的值大小。
[0069]
由于高超声速流动中的热流计算的准确性依赖于激波和温度边界层的精确捕捉，
而激波和温度边界层的精确捕捉又对计算中的数值粘性非常敏感。根据分子动力学理论，无量纲松弛时间τ0与压力和温度都相关，因此为了更准确控制计算域内数值粘性权重系数τ0的值，本发明将同时引入单元界面左右两侧的控制体单元的温度值和压力值对权重系数τ0的值进行控制。基于以上理论分析，提出了一种新的改进型开关控制函数来计算τ0的值，定义如下，
[0070][0071][0072][0073][0074][0075]
式中，tanh(x)是双曲函数中的双曲正切函数，tanh(x)函数值可以在[0,1]之间光滑过度。p
l
和p
r
分别为单元界面处左、右两侧的压力值，c为松弛因子常数，其取值范围为[1,100]，p
l
，p
r
和t
l
，t
r
分别为单元界面处左右两侧的压力值和温度值，δ为网格单元长细比修正系数，ω
l
、ω
r
分别为单元界面左右两侧控制体单元的体积，则为左右控制体单元的平均体积，l
l,r
为左右控制单元中心的距离。n
fl
和n
fr
分别表示单元界面左侧网格单元和右侧网格单元的面数，分别为第i个面上的权重系数。
[0076]
其具体求解流程如下：
[0077]
(1)计算各网格单元守恒变量的导数，通过线性重构得到单元界面两侧的守恒变量；
[0078]
(2)将单元界面两侧的守恒变量代入方程(6)
‑
(7)，计算出平衡分布函数和
[0079]
(3)通过方程(17)
‑
(18)计算出单元界面上的守恒变量，进一步地，采用方程(17)计算原始变量ρ
*
，和p
*
；
[0080]
(4)采用方程组(19)计算数值粘性权重系数τ0，然后通过方程(13)计算出单元界面处总的无粘通量的
[0081]
(5)计算出单元界面处的粘性通量的f
v
；
[0082]
(6)通过时间推进方法求解n
‑
s方程，这一步可以计算出下一时间步网格单元中心新的守恒变量值；
[0083]
(7)重复步骤(1)
‑‑
(7)，直到得到满足条件的收敛解。
[0084]
本发明具有如下有益效果：
[0085]
1.发明了可压缩格子boltzmann模型用于高超声速气动热计算
[0086]
采用一种非自由参数一维d1q4模型。该模型的分布函数构造简单，计算效率高，而
且粒子速度通过高阶动量矩关系式来确定，不需人工给定自由参数，更加合理。
[0087]
2.发明了一种新型基于网格效应的数值粘性控制函数
[0088]
本发明将同时引入单元界面左右两侧的控制体单元的温度值和压力值对权重系数τ0的值进行控制。基于以上理论分析，提出了一种新的改进型开关控制函数来计算τ0的值，
[0089]
3.发明了基于改进型fvm
‑
lbfs方法的流场计算新方法
[0090]
本发明为了拓展基于格子boltzmann模型数值方法在高超声速流域的应用范围，本发明提出了一种改进型fvm
‑
lbfs方法用于高超声速流动的数值模拟。该方法采用有限体积法对宏观n
‑
s控制方程进行空间离散，采用基于非自由参数d1q4模型的lbfs方法计算无粘通量，采用二阶中心差分方法计算粘性通量。通过新提出的压力和温度同时控制的改进型开关控制函数实现了对现有lbfs方法中无粘通量数值粘性的精确控制，并引入网格单元长细比修正系数拓展非均匀网格适用范围。
附图说明
[0091]
图1为计算网格。
[0092]
图2(a)为流场压力云图。
[0093]
图2(b)为流场温度云图。
[0094]
图2(c)为流场数值粘性控制函数云图。
[0095]
图3(a)圆柱表面压强分布。
[0096]
图3(b)圆柱表面热流分布。
[0097]
图4为界面迁移示意图。
具体实施方式
[0098]
下面结合应用实例对本发明进行说明。
[0099]
1.技术参数
[0100]
wieting.a.r(allen r w 1987nasa tm
‑
100484)于1987年在美国nasa兰利研究中心的8英尺的高焓风洞完成了不锈钢圆管前缘气动加热试验(wieting a r,holden m s 1987aiaa 22nd thermophysics conference honolulu hawaii june 8
‑
10,1987)，该试验已经被多次用于验证流场结构传热耦合计算的准确性。本发明同样选取条件完全相同的无限长不锈钢圆管模型，圆管尺寸为外径r
out
＝0.0381m,来流条件计算参数在表1中给出。
[0101]
表1来流条件参数
[0102][0103]
2.模型网格
[0104]
针对本算例，选取320
×
320的非均匀结构网格，且近壁面第一层网格高度的选取根据当地网格雷诺数确定，本算例来流雷诺数为1.835
×
105，为了计算方便取近壁面法向第一层网格雷诺数re
cell
＝1.835，符合文献中热边界层准确捕捉对网格雷诺数的要求。图1
为ma
∞
＝8.03高超声速圆柱绕流计算网格示意图。
[0105]
3.数值模拟结果
[0106]
图中分别给出了采用本发明所改进的lbfs方法计算得到的ma＝8.03高超声速圆柱流动的流场压力、温度和开关控制函数云图。从图2(a)和图2(b)我们可以看到，改进的lbfs方法可以稳定和清晰地捕捉到圆柱前的强激波区域，且未发出现任何明显数值振荡。另外，从图2(c)可以发现，本发明所改进的开关控制函数在高超声速流动中能起到了很好的数值粘性控制作用，在光滑流域内数值粘性权重系数τ0值接近于零，而在强激波区域周围τ0值接近于1，因此，τ0值的分布和预期理论保持一致，达到了控制函数的设计目标。
[0107]
图中分别给出了本发明计算得到的沿圆柱表面的压力(上)和热流(下)分布(以驻点值进行归一化)。另外，图中还分别给出了实验数据分布、采用ausm 格式的计算结果以及xu等人采用气体动力学bgk格式的计算结果，可以观察到采用本发明改进的lbfs方法预测的压力和热流分布与这些参考值吻合较好。此外，表2给出了采用lbfs方法计算得到的驻点处压力和热通量与参考文献xu.et al和yang.et al计算结果的比较。从表中可以看出，改进型的lbfs计算的结果与参考值吻合很好。
[0108]
通过以上数据对比表明，对于高超声速流动问题的数值模拟，本发明所改进型的lbfs方法可以同时稳定地捕捉强激波和准确地计算薄热边界层流动信息。
[0109]
表1驻点压力与热流值对比
[0110]