大涡模拟

1大涡模拟

目前计算机的计算能力仍对数值模拟紊流时所采用的网格尺度提出了严格

的限制条件。人们可以获得尺度大于网格尺度的紊流结构，但却无法模拟小于该

网格尺度的紊动结构。大涡模拟的思路是：直接数值模拟大尺度紊流运动，而利

用次网格尺度模型模拟小尺度紊流运动对大尺度紊流运动的影响

[2]

。大涡模拟较

直接数值模拟占计算机的内存小，模拟需要的时间也短，并且能够得到较雷诺平

均模型更多的信息。所以随着计算机的发展，大涡模拟越来越收到国内外研究者

的关注，并且认为大涡模拟将是最有前景的湍流模型。

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

[3]

:

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

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

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

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

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

[4]

：

u

i

u

i

u

j

1P

(



2S

ij)



tx

j



x

i

xj

(1-1)

式中：S拉伸率张量，表达式为：

S

ij

(u

i

/x

j

u

j

/x

i

)/2

；

分子粘性系数；



流体密度。根据LES基本思想，必须采用一种平均方法以区分可求解的大尺度

涡和待模化的小尺度涡，即将方程(1-1)中变量u变成大尺度可求解变量

u

。与雷

诺时间平均不同的是LES采用空间平均方法。设将变量

u

i

分解为方程(1-1)中

u

i和



次网格变量(模化变量)

u

i



,即

u

i

u

i

u

i

,

u

i

可以采用leonard提出的算式表示

为:

u

i

(x)

单而被广泛使用







G(xx



)u

i

(x



)dx

(1-2)

式中

G(xx



)

称为过滤函数，显然G(x)满足







G(x)dx1

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





1/



G

(x



x

)





0





xx



/2

xx

英语口语翻译

/2

(1-3)

这里



为网格平均尺度，三维情况下，

(

1



2



3

)

1/3

，



1

，



2

，



3

分别

为x

1

，x

2

，x

3

方向的网格尺度。当

0

时，LES即转变为DNS。

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

u

i



(u

i

u

j

)



1P



(



2S

ij)

tx

j



x

i

xj

(1-4)

由于无法同时求解出变量

u

i

和

u

i

u

j

，所以将

u

i

u

j

分解成

u

i

u

j

u

i

u

j





ij，



ij

即称为次网格剪切应力张量(亦称为亚格子应力)。

由此动量方程又可写成：

(2S

ij

)



ij

u

i

(u

i

u

j

)

1P







tx

j



x

i

x

j

xj

(1-5)

式中



ij

代表了小窝对大涡的影响。

上述叙述的过滤器属于非均匀过滤器，实际应用中还有均匀过滤器，例如盒

式过滤器、高斯过滤器、谱空间低通过滤器等等。

为了能够对



ij

进行模化，学者们提出了亚格子模型。

2亚格子模型

大涡模拟的基本思想就是对可解尺度湍流(或者讲大尺度湍流)直接数值求解，

但对不可解尺度湍流对可解湍流的影响由亚格子模型进行模化。亚格子模型一般

有以下集中类型

[5]

:唯象论的亚格子涡粘和涡扩散模型及其改进模型、结构性亚格

子模式、理性亚格子模式和其它亚格子模式。目前，在大涡模拟中经常广泛采用

的亚格子模型有标准的Smagorinsky模型、动态涡粘性模型、动态混合模型、尺

度相似模型、梯度模型、选择函数模型等

[6]

。其中Smagorinsky模型被广泛应用。

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

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

即

1



ij

2



t

S

ij





ij



kk

(2-1)

3

T

i



t



(2-2)

xi

以上公式中



t

和



t

分别称为亚格子涡粘系数和亚格子涡扩散系数；

S

ij

(1/2)[(u

i

/x

j

)(u

j

/x

i

)]

是可接尺度的变形率张量。式(2-1)第2项是为

了满足不可压缩的连蛇跟什么属相配 续方程，当

S

ij

收缩是(

S脑神经能恢复吗

ij

=0)等式两边可以相等。涡粘和涡

扩散模型的最大优点是计算方便，只要增加一个涡粘系数和涡扩散系数的模块，

就可以利用N-S方程的数值计算方法和程序。此外，整体上亚格子湍动能耗散或

亚格子标量能量耗散总是正值，因此涡粘和涡扩散模型的计算稳定性和鲁棒性也

较好。

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

u

i

u

i

u

i

u

j

p



kk

u

j

()[(







t

)()]

(2-3)

tx

i

x

i



3x

i

x

i

xi

ui

0

(2-4)

xi

2.2Smagorinsky模型

Smagorins梦到捉鱼 ky模型是由Smagorinsky于1963年提出来的，该模型是第一个亚

格子模型。文献[7]中是这样介绍Smagorinsky模型的：

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

1



ij





ij



kk

2



T

Sij

3

(2-6)

式中

S

ij

(1/2)[(u

i

/x

j

)(u

j

/x

i

)]

是可接尺度的变形率张量，



T

是涡粘系数。

1963年Smagorinsky定义了涡粘系数：



T

(C

S

)

2S

(2-7)

式中

S(2S

ij

S

ij

)

1/2

是变形率张量的大小，



是过滤尺度，C

S

无量纲参数，称为

Smagorinsky系数。

需要指出的是(2-7)式是根据各向同性湍流的能量输运推到的公式，在实际应

用中会发现Smagorinsky模型的一个致命的缺陷就是耗散过大。故文献[8]描述的

动态Smagorinsky模型可以弥补一些Smagorinsky模型的缺点。

动态Smagorinsky模型是基于为了减小Smagorinsky模型过大耗散的

Germano等式而得来的，1991年Lilly进行了改进。文献[1]对动态Smagorinsky

模型进行了详细的阐述。

为了表示简便，以



1

过滤的可解速度用上标“—”表示，以



2

过滤的可解

速度用上标“~”表示，一次过滤的Smagorinsky模型的亚格子偏应力公式为





ijf

1

(



kk

)

f



ij

2(



ij

)

f

S

ij

2C

D



2

1

SS

ij3

(2-8)

式中C

D

是取代Smagorinsky系数的动态系数，



1

是一次过滤的过滤长度。

假设过滤尺度



1

、



2

都是在惯性子区范围内，则以



2

尺度过滤的亚格子应

力的系数应当和



1

果过滤的系数相等，即





ijg

1

(



kk

)

g



ij

2(



ij

)

g

S

ij

2C

D



2

2

SS

ij3

L

(2-9)

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

ij

(



ij

)

fg

[(



ij

)

f

]g

(2-10)

由于

(



ij

)

fg

(



ij

)

g

，将式(2-8)和(2-9)带入(2-10)，有

L



L

ij

1

3

~~2



2C[SS

k威尔克斯 kijDij1

SS

ij]

2

2

(2-11)

可令

1

L

ij

L

kk



ij

C

D

Mij

3

(2-12)

由(2-12)不能直接结算系数C

D

，因为它是超定方程。有几种方法解决超定问题。

1)变形率张量收缩法

将式(2-7)两边同乘以可解尺度的变形率张量，关于C

D

的方程就确定了。实

际计算表明，这种方法计算的模式系数很不规则，计算的稳定性较差。

2)最小误差法

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

1

(L

ij

L

ij



ij

C

D

M

ij

)

2

0

(2-13)

C

D3

由上式可得

C

D

M

ij

L

ij

M

ij

Mij

(2-14)

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

改进。

C

D

M

ij

L

ij

M

ij

Mij

(2-15)

以上介绍的物理模型都是在物理空间的大涡模拟，文献[1]还介绍了谱空间

的涡粘模式，这里就不多做介绍。

3定解条件

虽然有了亚格子模型，大涡模拟方程已经封闭了，但是还需要合适、并且足

够多的定解条件才可以完全求解出结果。定解条件包括初始条件和边界条件。初

始条件要考虑湍流是均匀的还是切变的。边界条件包括：固体壁面、周期条件、

渐近条件、进口条件、出口条件和可压缩湍流的附加边界条件等等。在文献[1]

和相关的论著中有更详细的讲解这里就不介绍了。

至于使控制方程在空间和时间上离散的合适数值方法并不属于本论文的研

究内容。在文献[9]中，nier等人介绍了一种动态最优化的有限差分格式，

该动态格式是由泰勒级数展开式推导而来的，经过论证这种方法可以减小离散误

差，并且得到的结果和理论上的预测值很接近。

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

目前的湍流数值模拟可以分为直接数值模拟方法（DNS）镁和稀硫酸反应的化学方程式 和非直接数值模拟方法。所谓

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

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

大涡模拟和Reynolds平均法。

直接数值模拟（DNS）理论上可以得到相对准确的计算结果，但是由于湍流是多尺度的

不规则流动，要获得所有尺度的流动信息，需要很高的空间和时间分辨率，也就是需要巨

大的计算机内存和耗时很大的计算量，目前还无法应用于真正意义上的工程计算。

工程上广泛应用的是Reynolds平均法，这种方法的核心是不直接求解瞬时的N-S方程，

而是想办法求解时均化的Reynolds方程。这种方法的优点是不计算各种尺度的湍流脉动，

只计算平均运动，因此它的空间分辨率低，计算工作量小。缺点是模型没有普适性。