关于生成大涡模拟非定常进口条件的研究一直以来都是一个难题.在很多计算流体力学的数值模拟中,例如使用大涡模拟对叶轮机进行的数值模拟,计算结果在很大程度上受进口条件影响^{[1, 2]}.大涡模拟进口的流场需要符合湍流的统计特性,生成大涡模拟进口条件的方法要尽可能地容易操作,这样针对不同的进口情况能快速有效地生成满足要求的流场.近些年来,已经有一些方法可以生成大涡模拟的进口条件.其中一类是直接在平均流场上叠加一个在时间和空间上无相关关系的脉动信号,但这类方法得到的流场无法满足湍流的统计特性,与真实的湍流相差甚远.另一类方法是在进出口采用循环条件,初始速度场是在平均速度场上添加一个微小的扰动,然后代入到Navier-Stokes方程中计算,在流向计算长度足够长的情况下,能够得到充分发展的湍流,但由于流向计算长度太长,这种方法的计算成本较高.还有一类通过傅里叶函数合成方法,但这类方法受限于进口的几何形状,应用起来不够方便^[3].另外一类是涡方法,通过在进口处安放漩涡场,利用大小不同而且能够移动的漩涡场来生成脉动速度场^{[4, 5, 6]}.由于涡方法容易操作并且在流向上不需要太长的长度就能生成比较符合要求的湍流,所以本文针对此类方法进行了研究.采用由里昂流体与声学实验室的Sergent提出的涡方法^[4]对槽道湍流进行探索性数值计算实验.在改进涡方法中,主要通过Langevin方程来控制漩涡的运动以模仿实验方法中的蜂窝器,期望通过较短的流向长度来生成充分发展的湍流.

1 涡方法 1.1 涡方法的基本原理

非定常的脉动速度场需要至少满足以下条件:平均脉动速度为零;湍动能 为脉动速度);耗散率 为瞬时变形率张量,ν为运动黏度)与RANS(Reynolds-Averaged Navier-Stokes)计算结果匹配.

涡方法的主要目的是通过漩涡场生成湍流脉动速度场,基本原理构图如图 1所示.

图 1 涡方法基本原理构图 Fig. 1 Configuration of vortex method

图选项

二维情况下,涡方法是基于涡量的动力学方程:

理想中的漩涡在空间中的分布方式为狄拉克分布,那么式(3)的解析解可以通过对涡量进行卷积运算得到:

将式(4)代入到式(2),就可以得到速度场和漩涡场的关系表达式^[4](毕奥-萨伐尔定律):

在实际的三维情况下,漩涡在空间中的分布不可能像理想中的狄拉克分布那样,漩涡在空间中的分布由形状方程ξ_σ和环量值Γ来定义,其中σ是漩涡的作用半径.在任意一个位置x由漩涡i得到的旋量:

由式(6)，得到相应的脉动速度场为

全场的速度和漩涡可以用以下的关系式给出:

1.2 运用涡方法的基本流程

这个方法最早是由Sergent提出的,并使用商用计算流体力学软件Fluent进行了验证^[7].涡方法的使用流程具体如下.

在进口界面上随机分布一个漩涡场.

针对所添加的每一个漩涡,确定其环量值以及作用半径(σ和Γ).

在数值模拟计算的每一个网格上,计算出由漩涡场生成的脉动速度场,然后将其叠加到平均速度的剖面上.

控制漩涡的运动.例如,在进口界面上平移或者使其绕着某一轴旋转.

每隔一段时间(例如k/ε),可以调转漩涡的方向.

通过漩涡场计算出相应的脉动速度场.

2 改进后的涡方法的参数设定 2.1 漩涡的形状函数

漩涡的形状函数采用式(11)给出的形式,这种形状函数是改进的高斯函数型,与原始的涡方法中对形状函数的设定保持一致.

2.2 漩涡的作用半径

原始的涡方法中,Sergent用式(12)来给出漩涡的作用半径σ,这个表达式基于k和ε:

在改进后的涡方法中,对漩涡半径的参数设定是基于均匀各向同性湍流的假设.通过匹配由速度场表达的湍动能以及由漩涡场来表达的耗散率,可以得到近似于泰勒尺度的漩涡半径:

2.3 漩涡的环量和密度分布

基于均匀各向同性湍流的假设所得到的漩涡的环量与密度的关系表达式为

为了更为精确地确定在不同位置添加漩涡的数量,当地的漩涡数量应该和当地的湍动能和耗散率相关联.由此提出漩涡场的密度分布方式:

2.4 对漩涡场的运动的控制

原始的涡方法中,漩涡的移动是完全随机的.这样产生的一个副作用是,它会打乱脉动速度场在时间和空间的相关性.因此改进后的涡方法通过Langevin方程^{[8, 9]}来控制漩涡的运动:

实质上,这个过程是模仿实验方法中的蜂窝器来生成湍流^[10].蜂窝器通过网格来扰动平均速度场,进而生成脉动速度,这与添加的漩涡场的作用是一样的.蜂窝器的网格所具备的分形结构与通过控制方程来控制漩涡场的运动具有相似性.以Langevin方程控制漩涡的运动,来重构流场在时间和空间的相关关系,以便尽快地生成湍流,这与蜂窝器的网格的分形结构所起的作用是一样的.

3 数值模拟实验

通过一个雷诺数为395的槽道湍流算例来检验改进后的涡方法.数值模拟使用的程序为法国里昂流体与声学实验室(LMFA)开发的Turb’Flow.槽道的大小为2πh×2h×πh,h是槽道半高(h=0.01m).3个方向的网格数为49×89×41.流向和展向网格均匀分布,y方向的网格分布如下:

图 2 槽道的网格分布示意图 Fig. 2 Meshes of channel

图选项

为了加速计算,整个槽道并行分配为8个网格块.槽道y方向的两个面是无滑移壁面,x方向是流向,进口处的边界条件通过涡方法添加,出口处为压力边界.z方向采用周期性边界条件.

数值模拟计算所需要的平均场由RANS计算得到,所采用的RANS模型是Wilcox k-ω模型^[11].大涡模拟中使用的亚格子模型为WALE模型^{[12, 13]}.

为了和原始的涡方法进行对比,第1种模拟计算采用了在进口界面随机分布的800个漩涡来生成漩涡场,所有漩涡在进口界面的平移也是随机的,按照式(11)、式(13)和式(14)来设定漩涡的形状函数、作用半径和环量值,这与原始的涡方法是一致的.第2种计算采用了改进后的涡方法,按照式(11)、式(13)和式(14)设定漩涡的形状函数、作用半径和环量值,但漩涡场采用式(15)的密度分布方式,总的漩涡数目依然为800个,并且运用Langevin方程(16)来控制漩涡场的运动.试验结果的参考数据来自于DNS(Direct Numerical Simulation)的计算结果^{[14, 15, 16]}.

4 计算结果及分析 4.1 平均速度场

实验所得到的平均速度场的结果由图 3和图 4给出.在进口处,由于平均速度场是根据RANS的计算结果给定,所以无论是在平均速度的线性区还是对数区,它都比较符合直接数值模拟的结果.在出口处,平均速度也能和直接数值模拟的数据相吻合.因为湍流沿着流向不断发展过程中,脉动速度越来越符合湍流的脉动统计特性,对平均速度场统计的影响越来越小.

图 3 进口平均速度场的统计结果 Fig. 3 Mean velocity of inlet

图选项

图 4 出口平均速度场的统计结果 Fig. 4 Mean velocity of outlet

图选项

脉动涡量ω′的结果由图 5和图 6给出.很容易观察到,在进口处,传统的涡方法不能够得到准确的统计结果.而采用密度分布方式并且由Langevin方程来控制漩涡运动的改进涡方法能够生成与DNS统计数据相吻合的脉动涡量场,这是由于采用密度分布漩涡的方式能准确地根据不同位置的需求量来分布漩涡的数目,从而得到更接近真实湍流的脉动速度场.

图 5 进口处的涡量统计结果 Fig. 5 Vorticity of inlet

图选项

图 6 出口处的涡量统计结果 Fig. 6 Vorticity of outlet

图选项

流体向下游流动的过程中,湍动能在不同尺度的脉动间传递,直至传递到Kolmogorov尺度并耗散掉,这个过程中能量会衰减.并且由于进口处添加的漩涡场产生的脉动速度场在统计上不能完全符合湍流的统计特性,这种初始场在通过Navier-Stokes方程求解的过程中,会摧毁掉初始场所包含的流场信息.但从曲线的走势来看,出口处,曲线更加平滑,更加趋近于DNS的统计曲线.对比DNS数据,在大小上约有一个量级的误差.

4.3 雷诺应力

图 7和图 8给出了进出口处的雷诺应力的统计结果.进口处,改进后的涡方法中,漩涡场的分布跟当地的湍动能和耗散率相关,所以得到的脉动速度的统计更接近DNS结果.在出口处,由于脉动场得到了一定程度的发展,所以雷诺应力的统计曲线更为平滑,改进后的涡方法得到的〈uu〉和〈uv〉的统计曲线更为接近DNS结果.同时从统计结果可以看出,流体在向下游流动的过程中,湍动能有所衰减,初始场所包含的信息在通过Navier-Stokes方程求解的过程中在慢慢丢失,流体更加符合湍流真实的统计特性.

图 7 进口处的雷诺应力统计结果 Fig. 7 Reynolds stress of inlet

图选项

图 8 出口处的雷诺应力统计结果 Fig. 8 Reynolds stress of outlet

图选项

实验证明,涡方法能够通过RANS的结果来生成脉动场.改进后的涡方法采用密度分布方式来分布进口漩涡场并通过Langevin方程控制漩涡的运动.

1) 相比于原始的涡方法,对比涡量和雷诺应力的统计结果,不论是在槽道的进口处,还是出口处,改进后的涡方法生成的脉动速度场更加接近真实的湍流.

2) 在出口处,湍流得到一定程度的发展,从统计结果来看,更加符合真实湍流的统计特性,这个过程中,湍流所需要的流向发展长度约为10h,对比传统的使用循环条件^[2]来生成大涡模拟进口条件的方法,涡方法显得更加有效.

3) 另一方面,从进出口处雷诺应力的统计数据来看(〈w′w′〉和〈v′v′〉),改进后的涡方法所生成的脉动速度场还不能完全和DNS数据契合,这方面还需要进一步的探索和改进.

致谢 感谢北京航空航天大学的方乐老师以及法国里昂流体与声学实验室的高峰博士对我的帮助和指导.