湍流模型-大涡模拟

1 大涡模拟

目前计算机的计算能力仍对数值模拟紊流时所采用的网格尺度提出了严格的限制条件。人们可以获得尺度大于网格尺度的紊流结构，但却无法模拟小于该网格尺度的紊动结构。大涡模拟的思路是：直接数值模拟大尺度紊流运动，而利用次网格尺度模型模拟小尺度紊流运动对大尺度紊流运动的影响[2]。大涡模拟较直接数值模拟占计算机的内存小，模拟需要的时间也短，并且能够得到较雷诺平均模型更多的信息。所以随着计算机的发展，大涡模拟越来越收到国内外研究者的关注，并且认为大涡模拟将是最有前景的湍流模型。

使用大涡模拟的时候，要注意以下4个问题[3]:

1) 用于N-S 方程进行过滤的函数。

2) 彻底经过经验封闭的模型(包括传统亚格子模型和其它封闭方法)。

3) 足够多的边界条件和初始条件。

4) 使控制方程在空间和时间上离散的合适数值方法。

不可压缩常粘性系数的紊流运动控制方程为N-S 方程[4]：

(1-1) 式中：S 拉伸率张量，表达式为：2/)//(i j j i ij x u x u S ∂∂+∂∂=；γ分子粘性系数；ρ流体密度。根据LES 基本思想，必须采用一种平均方法以区分可求解的大尺度涡和待模化的小尺度涡，即将方程(1-1)中变量u 变成大尺度可求解变量u 。

与雷j ij i j j i i x S x P x u u t u ∂∙∂+∂∂-=∂∂+∂∂)2(1γρ

诺时间平均不同的是LES 采用空间平均方法。设将变量i u 分解为方程(1-1)中i u 和

次网格变量(模化变量)'i u ,即'+=i i i u u u ,i u 可以采用leonard 提出的算式表示

为:

(1-2)

式中)(x x G '-称为过滤函数，显然G(x)满足

常用的过滤函数有帽型函数(top —hat)、高斯函数等。帽型函数因为形式简单而被广泛使用

(1-3) 这里∆为网格平均尺度，三维情况下，3/1321)(∆∆∆=∆，1∆，2∆，3∆分别为x 1，x 2，x 3 方向的网格尺度。当0→∆时，LES 即转变为DNS 。

将过滤函数作用与N-S 方程的各项，得到过滤后的紊流控制方程组：

(1-4) 由于无法同时求解出变量i u 和j i u u ，所以将j i u u 分解成ij j i j i u u u u τ+∙=，ij τ即称为次网格剪切应力张量(亦称为亚格子应力)。

由此动量方程又可写成：

(1-5)

式中ij τ代表了小窝对大涡的影响。

上述叙述的过滤器属于非均匀过滤器，实际应用中还有均匀过滤器，例如盒式过滤器、高斯过滤器、谱空间低通过滤器等等。

为了能够对ij τ进行模化，学者们提出了亚格子模型。

2 亚格子模型

大涡模拟的基本思想就是对可解尺度湍流(或者讲大尺度湍流)直接数值求解，但对不可解尺度湍流对可解湍流的影响由亚格子模型进行模化。亚格子模型一般x d x u x x G x u i i '''-=

⎰+∞∞-)()()(⎰+∞∞-=1)(dx x G ⎪⎩⎪⎨⎧∆>'-∆≤'-∆='-2/02//1)(x x x x x x G j ij i j j i i x S x P x u u t u ∂∙∂+∂∂-=∂∂+∂∂)

2(1)(γρj ij j ij i j j i i x x S x P x u u t u ∂∂-∂∂+∂∂-=∂∙∂+∂∂τγρ)2(1)(

有以下集中类型[5]:唯象论的亚格子涡粘和涡扩散模型及其改进模型、结构性亚格子模式、理性亚格子模式和其它亚格子模式。目前，在大涡模拟中经常广泛采用的亚格子模型有标准的Smagorinsky 模型、动态涡粘性模型、动态混合模型、尺度相似模型、梯度模型、选择函数模型等[6]。其中Smagorinsky 模型被广泛应用。

2.1 亚格子涡粘和涡扩散模型[1]

不可压缩湍流的亚格子涡粘和涡扩散模型采用分子粘性和分子热扩散形式，即

kk ij ij t ij S τδντ3

12+= (2-1) i

t i x T ∂∂=θκ (2-2) 以上公式中t ν和t κ分别称为亚格子涡粘系数和亚格子涡扩散系数；)]/()/[()2/1(i j j i ij x u x u S ∂∂+∂∂∙=是可接尺度的变形率张量。式(2-1)第2项是为了满足不可压缩的连续方程，当ij S 收缩是(ij S =0)等式两边可以相等。涡粘和涡扩散模型的最大优点是计算方便，只要增加一个涡粘系数和涡扩散系数的模块，就可以利用N-S 方程的数值计算方法和程序。此外，整体上亚格子湍动能耗散或亚格子标量能量耗散总是正值，因此涡粘和涡扩散模型的计算稳定性和鲁棒性也较好。

将亚格子应力的涡粘模型公式(2-1)代入到(1-5)式中，变形得

)])([()3(i

j i i t i kk i i i j i x u x u x p x x u u t u ∂∂+∂∂+∂∂++∂∂-=∂∂+∂∂νντρ (2-3) 0=∂∂i

x u i (2-4) 2.2 Smagorinsky 模型

Smagorinsky 模型是由Smagorinsky 于1963年提出来的，该模型是第一个亚格子模型。文献[7]中是这样介绍Smagorinsky 模型的：

广泛用于大涡模拟中的涡粘模型认为亚格子应力的表达式如下：

(2-6)

ij T kk ij ij S ντδτ231-=-

式中)]/()/[()2/1(i j j i ij x u x u S ∂∂+∂∂∙=是可接尺度的变形率张量，T ν是涡粘系数。1963年Smagorinsky 定义了涡粘系数：

(2-7)

式中2/1)2(ij ij S S S =是变形率张量的大小，∆是过滤尺度，C S 无量纲参数，称为Smagorinsky 系数。

需要指出的是(2-7)式是根据各向同性湍流的能量输运推到的公式，在实际应用中会发现Smagorinsky 模型的一个致命的缺陷就是耗散过大。故文献[8]描述的动态Smagorinsky 模型可以弥补一些Smagorinsky 模型的缺点。

动态Smagorinsky 模型是基于为了减小Smagorinsky 模型过大耗散的Germano 等式而得来的，1991年Lilly 进行了改进。文献[1]对动态Smagorinsky 模型进行了详细的阐述。

为了表示简便，以1∆过滤的可解速度用上标“—”表示，以2∆过滤的可解速度用上标“~”表示，一次过滤的Smagorinsky 模型的亚格子偏应力公式为

(2-8)

式中C D 是取代Smagorinsky 系数的动态系数，1∆是一次过滤的过滤长度。

假设过滤尺度1∆、2∆都是在惯性子区范围内，则以2∆尺度过滤的亚格子应力的系数应当和1∆果过滤的系数相等，即

(2-9) 这里需要Garmano 等式，故附上式(2-10)。 (2-10)

由于g ij fg ij )()(ττ=，将式(2-8)和(2-9)带入(2-10)，有

(2-11)

可令

(2-12)

S C S T 2)(∆=ν()ij D ij f ij ij f kk f ij S S C S 2

12)(2)(31∆==-τδττ()ij D ij g ij ij g kk g ij S S C S 2

22)(2)(3

1∆==-τδττg f ij fg ij ij L ])[()(ττ-=]~~[231212

2ij ij D ij kk ij S S S S C L L ∆-∆=-δij D ij kk ij M C L L =-δ31

由(2-12)不能直接结算系数C D ，因为它是超定方程。有几种方法解决超定问题。

1) 变形率张量收缩法

将式(2-7)两边同乘以可解尺度的变形率张量，关于C D 的方程就确定了。实际计算表明，这种方法计算的模式系数很不规则，计算的稳定性较差。

2) 最小误差法

令式(2-12)两边的平方差最小，即

0)3

1(2=--∂∂ij D ij ij ij D M C L L C δ (2-13) 由上式可得

ij ij ij

ij D M M L M C = (2-14)

最小误差法较之变形率张量收缩法有很大的该进，但是还是有缺陷，故用下式改进。

ij ij ij

ij D M M L M C = (2-15)

以上介绍的物理模型都是在物理空间的大涡模拟，文献[1]还介绍了谱空间的涡粘模式，这里就不多做介绍。

3 定解条件

虽然有了亚格子模型，大涡模拟方程已经封闭了，但是还需要合适、并且足够多的定解条件才可以完全求解出结果。定解条件包括初始条件和边界条件。初始条件要考虑湍流是均匀的还是切变的。边界条件包括：固体壁面、周期条件、渐近条件、进口条件、出口条件和可压缩湍流的附加边界条件等等。在文献[1]和相关的论著中有更详细的讲解这里就不介绍了。

至于使控制方程在空间和时间上离散的合适数值方法并不属于本论文的研究内容。在文献[9]中，D. Fauconnier 等人介绍了一种动态最优化的有限差分格式，该动态格式是由泰勒级数展开式推导而来的，经过论证这种方法可以减小离散误差，并且得到的结果和理论上的预测值很接近。

一、 湍流数值模拟方法简介

目前的湍流数值模拟可以分为直接数值模拟方法（DNS ）和非直接数值模拟方法。所谓直接数值模拟方法是指直接求解瞬时的湍流控制方程（N-S 方程），无需对湍流流动作任何简化或近似；而非直接数值模拟方法是不直接计算湍流的脉动特性，而是设法对湍流作某

种程度的近似和简化处理，依赖所采用的近似和简化方法不同，非直接数值模拟可以分为大涡模拟和Reynolds平均法。

直接数值模拟（DNS）理论上可以得到相对准确的计算结果，但是由于湍流是多尺度的不规则流动，要获得所有尺度的流动信息，需要很高的空间和时间分辨率，也就是需要巨大的计算机内存和耗时很大的计算量，目前还无法应用于真正意义上的工程计算。

工程上广泛应用的是Reynolds平均法，这种方法的核心是不直接求解瞬时的N-S方程，而是想办法求解时均化的Reynolds方程。这种方法的优点是不计算各种尺度的湍流脉动，只计算平均运动，因此它的空间分辨率低，计算工作量小。缺点是模型没有普适性。