湍流模型-大涡模拟

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