以分离涡模拟^[1](DES,detached eddy method)为代表的整体式(Global)RANS/LES混合方法自1997年被提出以来得到了长足的发展。通过对传统湍流模型的简单改造,它们能够根据湍流长度尺度和网格疏密分布自动划分RANS区和LES区:壁面附近由RANS模化,尽可能降低计算量;远场分离区RANS模型方程转化为类亚格子模型的形式,显著降低了湍流黏性,起到了类似LES隐式滤波的效果。

值得注意的是,这类以“LES在RANS之上”(二者交界面与流向平行)为特点的整体式RANS／LES混合方法的一个固有缺陷就是RANS区和LES区交界面不可控。这里的“不可控”,除了体现在交界面位置很难事先确定进而难以精确规划LES区域外,更重要的是在交界面并未对往来RANS区和LES区之间的交互信息进行修正,特别是RANS区较大的湍流黏性直接进入LES区,严重抑制了后者对湍流的解析能力,延迟了湍流的发展,即出现了所谓的灰区(grey area)。针对这一情况,一种直观的解决思路就是在二者交界面附加某种“合适”的湍流脉动,这样流动从RANS区一进入LES区瞬间即由模化量转变为解析量,大大缩短了从RANS到LES的过度区间(即为灰区的长度)。

鉴于此,嵌入式(embedded)RANS／LES混合方法^{[2, 3, 4]}得到了越来越多关注。与整体式方法不同,嵌入式混合方法需要事先人为将若干LES区嵌入全场的RANS区——这里“嵌入”的目的是让最终的RANS和LES分区交界面尽可能与流向垂直,即“LES在RANS之后”,从而保证了流动方向与从RANS到LES信息传递方向的一致性,这样在二者交界面能够更为合理更为容易地引入额外的湍流信息来促进LES区的湍流进一步向下游发展这种方法特别适用于一些几何外形相对简单的局部精细流动模拟。

本文以低雷诺数槽道流动为算例,为便于分析,直接将其流场简单分割为上游RANS区和下游LES区,以SST-DES为基本模型,重点观测采用合成涡方法(SEM,synthetic eddy method)生成的交界面湍流脉动对数值模拟结果的影响。

点击浏览公式

点击浏览公式

点击浏览公式

对于嵌入式混合方法而言,在上游RANS区

点击浏览公式

即为标准SST两方程模型。在下游LES区

点击浏览公式

即为Yoshizawa亚格子模式。另外,借鉴文献的思路,由于嵌入式混合方法中RANS区和LES区沿流向依次分布,上游RANS产生的过大湍动能直接冲击下游LES区,其不利影响甚至超过了整体式混合方法。因此有必要对交界面处的湍动能进行限制:

点击浏览公式

比耗散率则由Smagorinsky亚格子模式表征:

点击浏览公式

限于篇幅,SST-DES方法中的其他参数选择不再详述,可参考文献[5]。

正如上文所述,在绝大多数流向运动占主导地位的数值模拟中,这种“LES在RANS之后”的嵌入式方法让下游LES区直接承接上游RANS的模化结果,必须辅之以额外的湍流脉动信息抵消其不利影响。即在二者的交界面,有

点击浏览公式

式中,a为小于1的常数——这是由于添加的湍流脉动速度分量u′,在交界面附近产生了额外的速度梯度,相当于改变了动量方程的源项,导致最终得到的脉动强度往往比预期的要大。这样,选用一种正确的、合理的方法生成交界面的湍流脉动速度成为运用嵌入式混合方法的核心所在,除了生成的脉动场平均值为零外,它还应当满足下面2个最基本的条件:

1) 脉动速度场要与目标流场物理契合。也就是说,它的若干物理特性,如能谱、长度尺度和各向异性程度等要尽可能接近目标流场;

2) 脉动速度场必须满足时间关联性。在数值模拟中,每一个时间步长都有新的脉动速度添加到交界面,而这些脉动速度不可能是时间独立的“白噪声”。

合成涡方法是生成湍流脉动的一种简洁高效的方法^{[4, 6]}。该方法假设目标脉动速度是由空间的若干个随机分布的点涡共同诱导产生的,可写做

点击浏览公式

点击浏览公式

N和X^k分别为为空间点涡的总数和三维坐标,U为流动的时均速度。经过时间步长dt后,点涡坐标运动了Udt(交界面坐标不变),在新的位置重新计算在交界面的诱导速度,保证了每个时间步的二维脉动场与上个时间步结果的时间关联。f为点涡诱导速度分布的形函数

点击浏览公式

点击浏览公式

式中,σ为交界面每个单元的长度尺度(也是计算该单元脉动速度时点涡产生的最大诱导距离),可直接采用RANS的计算结果。另外在(9)式中还有一个表征脉动强度的重要随机变量ε

点击浏览公式

即以相等概率(50％)取值为-1或1。这样最终脉动速度场的期望(平均脉动速度值)与方差(雷诺应力张量的分量)为

点击浏览公式

也就是主分量均为1(实际应用中可根据情况将其放大或缩小)的各项同性速度场。对于一些已有DNS或实验结果的经典算例来说,还可采用如下步骤进一步优化脉动速度场^[2]:

1) 从DNS或实验数据库中提取出各项异性程度最大位置的雷诺应力张量。对于展向流动不明显的算例,可直接选取〈u′v′〉绝对值最大处。

2) 计算该参考应力张量的特征主方向。

3) 将点涡和交界面单元的空间坐标均旋转到该主方向,运用上文的合成涡方法得到新坐标系下的脉动速度场(仍为各项同性)。

4) 将脉动速度旋转回原始坐标系下,此时得到与步骤1)中提取应力张量相符的各项异性场,但其主分量在每个坐标轴方向上都是不变的,并未考虑到实际算例中壁面的影响。

5) 进一步用DNS或实验得到的脉动速度型曲线(图 1)中波浪形实线即为〈u′v′〉的分布)在壁面方向规约主分量,得到最终的湍流脉动速度场。

Re=950下的槽道流动是一个测试湍流模型表现的常规算例,其DNS的结果数据库可参考文献^[7],本文选取用于合成涡方法的参考雷诺应力为

点击浏览公式

槽道流计算域在流向(x)、法向(y)和展向(z)的长度分别为无量纲长度6.0,2.0和2.0,与之对应的网格单元数为128×80×64,距离壁面最近一层网格y⁺约为0.8,流向和展向网格均匀分布——该网格事实上近似达到了全场LES的分辨率。RANS和LES分区与图 1相同,即交界面为x=1.5处的yz平面,RANS区和LES区分别采用如(4)式和(5)式的SST-DES控制方程。为减少计算量,本算例中流向长度相对较短,RANS区的入口速度型直接写作:

点击浏览公式

图 1 槽道流动计算域及分区示意图

图选项

在(10)式中取其平均值U=17.5。算例共迭代了20 000步,从12 000步开始统计平均,生成的脉动速度场相当于一个有20 000×80×64个单元的长方体区域,每个时间步长有一个包含脉动速度信息的yz平面“推进”到流场交界面并附加在通过RANS得到的模化速度上。

为便于描述,下文用solver1和solver2分别指采用附加脉动速度和无附加脉动速度的嵌入式混合方法的模拟结果。图 2和图 3分别给出了湍流黏性系数和雷诺应力主分量〈u′v′〉的最大值沿流向变化曲线。当附加脉动速度时,LES区湍流黏性几乎立即下降到0附近,而solver2中湍流黏性仍停滞在高位,严重抑制了LES的解析能力。与之印证的是,solver1中的〈u′v′〉从RANS区的近似为0在通过交界面后快速上升到10左右,与(16)式中给定的7.992 5基本吻合。而solver2中的〈u′v′〉则在LES区始终处于缓慢爬升的过程,这一过度区域的长度甚至远超过上游RANS的长度。

图 2 湍流黏性系数分布

图选项

图 3 雷诺应力主分量〈u′v′〉的最大值沿流向分布

图选项

图 4和图 5显示时均速度在x=5.0截面和壁面摩擦阻力因数在流向的变化趋势。与上文的分析类似,相比与solver2,solver1不但在LES区出口附近呈现出非常贴合RANS入口的速度型(即为(17)式),而且通过交界面附加脉动的作用,摩擦阻力值在较短的调节过度后即恢复到RANS区的水平。而solver2受上游过大湍流粘性影响,LES区网格分布明明足够细密却依然无法得到合理的结果,解析湍流发展缓慢,摩擦阻力因数严重偏低。

图 4 时均速度分布(x=5.0)

图选项

图 5 摩擦阻力因数分布

图选项

进一步,图 5依次给出了3个关键雷诺应力分量〈u′u′〉,〈v′v′〉和〈u′v′〉在x=5.0界面的法向分布。总得来看,solver1得到的脉动速度曲线与DNS的结果非常接近,再一次证明了基于合成脉动的嵌入式方法的正确性,保全了LES的解析湍流能力。而solver2最终在LES区得到的湍流脉动仍旧非常微弱,如果想改善其结果只能继续加长计算域流向长度来让LES区越过灰区,避免来自RANS区的数值污染。

图 6 雷诺应力分量在x=5.0截面的法向分布

图选项

本文论述了一种基于SST-DES模型,采用合成涡方法生成RANS和LES交界面脉动信息的嵌入式混合方法,并以低雷诺数槽道流动为算例验证了其正确性。相比传统整体式方法,嵌入式方法通过人为事先划分RANS和LES区域,在二者交界面附加额外的湍流脉动信息以实现降低RANS对LES的影响,促进LES发展的目的。显然,这个额外湍流脉动信息的合理与否对嵌入式方法的最终模拟结果有着根本性的影响。就本文内容而言,合成涡方法一方面基本达到了上述目标,另一方面依旧不可避免地引入了DNS信息作为辅助,如何将其进一步改进,提高其先验性仍有待进一步工作的研究。