壁函数对大气边界层数值模拟结果的影响

壁函数对大气边界层数值模拟结果的影响
由于湍流边界层在靠近壁面区域存在动态的涡旋结构，所以在高雷诺数情况下大涡模拟(LES)要有非常高的分辨率，但是解决这些涡旋结构使得大涡模拟的计算量和直接数值模拟一样大。

为了避免这种情况，我们可以用壁函数来模拟近壁的网格，它在固体边界上给大涡模拟提供了近似边界条件，使大涡模拟在高雷诺数条件下可以增大网格分辨率。

本文的目的就是为了模拟不同壁函数对湍流边界层数值模拟结果的影响，模拟流动的控制方程为纳维斯托克斯方程。

模拟采用全隐式的解耦方法,在LU分解和近似分解的基础上，速度项和压力项被解耦，同时保留了时间二阶精度。

在未加壁函数时，文章中使用了两种不同的亚格子模型分析了不同亚格子模型对湍流边界层数值模拟结果的影响，最后分析了不同壁函数、不同网格对结果的影响，并且将得到的结果与Lee和Moser的结果进行对比，使湍流边界层的模拟更加准确，并将结果推广到了大气边界层。

1.1 壁面模型介绍
模拟流体流动在工程设计和分析中具有重要的应用价值。

因为层流的流动特征在时间和空间上的尺度非常小，而湍流的流动特征相对层流而言较大，所以模拟层流到湍流的转换过程是很困难的，这对动态流动的精确模拟产生了很大的阻碍。

在直接数值模拟的过程中，不需要添加其他模型就可以解决所有尺度的流动，但如果要模拟整个流动过程，所需的网格点非常多。

为了减少直接数值模拟的计算量，发展出了一种新的数值模拟方法——大涡模拟错误!未找到引用源。

，通过对湍流进行低通滤波，计算量显著减小，但是经过过滤，大涡模拟消除了许多小尺度结构。

从物理和工程学的角度来看，高频信息对实际问题的重要性不大，但是，它所携带的物理信息对流动的发展有重要影响，所以要对高频信息建立模型来模拟其对大涡模拟的影响，这种模型叫做亚格子模型。

目前已经发展出了许多有效的模型和方法错误!未找到引用源。

应用亚格子模型后，大涡模拟可以准确地应用于多种流动状态。

但是在近壁区域，因为流动是粘性的，所以亚格子模型在近壁区域不太准确，除此之外，这个区域的流动结构趋于各向异性，而亚格子模型主要用来对各向同性涡进行建模，各向同性涡只代表一小部分的总能量流，他们不能精确代表壁面附近的湍流
应力，所以要解决近壁剪切应力产生的涡旋结构，所需要的网格点非常多，这就使得大涡模拟的近壁分辨率与直接数值模拟一样高。

为了模拟高雷诺数的湍流边界层，壁面模型逐渐被发展起来，壁面模型给大涡模拟提供了边界条件，所以大涡模拟不需要直接解决近壁区域。

当不直接计算靠近壁面的区域时，计算量显著减少。

一个典型的壁面模型错误!未找到引用源。

是在固体壁面用近似边界条件代替标准的无滑移速度边界条件，所以在不计算壁面区域的情况下，大涡模拟可以精确计算外部流动。

壁面附近除了垂直壁面方向大的速度梯度外，也包含许多有条纹结构，这些结构对湍动能和剪切应力的产生及输运是非常重要的，因为近似边界条件对这些结构有影响，所以近似边界条件对大涡模拟非常重要。

近似边界条件错误!未找到引用源。

有许多优势，其中一个就是壁面应力，为了更精确的模拟流动结构，必须保证近似边界条件的准确性。

另外，通过已知的平均速度廓线可以使壁面应力与流动状态直接相关，利用这种关系，已经发展出了许多模型。

1.2 代数壁面模型
壁面模型主要分为三种：第一种为代数模型，代数模型在大涡模拟和壁面应力之间使用了一个简单的关系式；第二种为双层模型错误!未找到引用源。

，在双层模型中，壁面应力通过近壁动力学来确定；第三种为基于控制的壁面模型错误!未找到引用源。

，在这种壁面模型中引入了控制器，通过壁面应力的输入来控制大涡模拟。

这些模型共享的一个额外特点是：当垂直壁面的速度为零时，仅仅规定了与壁面速度分量平行的壁面应力。

用这种方法纯粹地从物理推导中很难确定法向速度，因为它的分量和垂直壁面方向的导数在壁面上为零。

一个额外的困难就是，如果法向速度非零，必须使穿过壁面的净质量输运为零，这就意味着，不可能根据局部大涡模拟数据来确定法向速度，需要从壁面上获得其他信息，所以在下面的讨论中假设法向速度为零。

本文中使用的壁面模型是代数应力模型，即壁面应力模型，下面将壁面应力模型做一个简要介绍。

壁面应力模型错误!未找到引用源。

是解决近壁区域的第一类方法，这种方法用壁面应力代替了传统的无滑移边界条件，所以不需要直接解决靠近壁面的湍流运动。

Deardorff （1970）首先应用了这个模型，他把这个模型应用于有限雷诺数平面槽道流的大涡模拟中：
2222211u u y y z
κ∂∂=-+∂∂ （1）
2222w w y x ∂∂=∂∂ （2）
其中u 和w 分别为流向和展向的速度分量，1y 为远离壁面第一个网格点的坐标，
κ为冯卡门常数，这些边界条件是加在速度的二阶导数上的，与此同时，利用这些边界条件可以在壁面画出一条对数廓线，当这个边界条件和壁面条件结合时，u 和w 上的条件提供了模拟所需要的全部边界条件，使用这个模型，Deardorff 可以计算平面槽道中的流动，但平均统计数据和实验数据并不是很一致，这不仅仅归因于壁面模型，还与外部的网格分辨率有关。

Schumann 使用有限体积法和交错网格，计算了有限雷诺数的湍流流动，实现了大涡模拟平面槽道流的标准壁面应力模型错误!未找到引用源。

，在计算过程中假设壁面应力与第一个内部网格点上的速度方向相同，从而确定了壁面应力，特别地，使用了以下模型：
1211()()()
w w t u u y y u y ττνν++<>∂=+=∂<> （3） 13()w t w w y y
τννν∂=+=∂ （4） 其中，<>代表水平平均，ν为分子粘度，t ν为涡流粘度，w τ<>代表平均流向
壁面应力。

通过假设边界层处于平衡状态，反复迭代，使1y 处（槽道内层的第一
个网格点）的水平平均流向速度满足壁面的对数律。

相比Deardorff 的粗分辨率计算，这个模型在槽道流中获得了比较好的结果。

现在对这个模型已经提出了一些改进，例如Piomelli （1989）提出的，在计算中考虑靠近壁面涡旋结构的倾斜；Grotzbach 使用了这种模型在壁面上施加了热量通量，计算过程中涉及到热传递。

像在之前提到的，壁面模型对于模拟环境流动有重要作用，在环境流动中，壁面应力根据局部和瞬时的对数廓线来定义（Mason 和Callen ），除对数律以外，其他平均速度廓线也被用来计算壁面应力，例如，在Werner 和Wengle （1991）的工作中使用了近壁线性廓线，并且用这些廓线预测得到的结果与Schumann 和Piomelli 得到的结果相同。

结合Mason 和Callen 的壁面模型，Mason 和Thomson 提出了一个随机后散射模型错误!未找到引用源。

，这个模型给近壁区的纳维斯托克斯方程施加了一个随机力，分析了从小尺度到大尺度过程中能量后向散射的影响，通过调整这个力的振幅，显著改善了Mason 和Callen 得到的平均速度廓线，Mason 、Thomson 和
Piomelli在报告中指出：这个力“打碎”了大尺度结构，并且对速度场的影响不大，即没有明显改善外流中的平均速度，而且，从瞬时流中可看出，产生的流动结构和边界层上的结构并不相符。

除此之外，目前还没有方法选择随机力的振幅，因此这个壁面模型不能广泛应用于其他情况，但是证明了，为了更好的预测平均速度廓线，必须对标准壁面模型进行修正。

1.3 本文研究内容
本文主要应用了Schumann和Werner提出的壁面模型来模拟壁面对外部流动的影响，在Schumann模型中，Grotzbach把Schumann的方法推广到了压力梯度不已知的情况。

在本文中，首先在低雷诺数情况下对不加壁函数的程序用直接数值模拟进行验证，引入亚格子模型后,同时对大涡模拟做了验证,并将两种数值模拟得到的结果与Lee和Moser的结果进行对比，然后讨论了不同亚格子模型对湍流边界层模拟结果的影响，最后将两种壁函数加到现在的程序中，检验不同壁函数、不同网格对湍流边界层结果的影响。

分别画出平均速度廓线、速度脉动和雷诺应力随湍流边界层高度的变化，与Lee和Moser得到的结果进行对比，并将结果推广到大气边界层，得出结论。

第二章对纳维斯托克斯方程的隐式速度解耦过程
2.1 介绍
随着直接数值模拟和大涡模拟这两种数值模拟方法越来越进步，出现了很多求解不可压缩纳维斯托克斯方程的有效数值算法，其成功的核心是对耦合的不可压缩动量方程和连续性方程解耦。

对文献的精读发现，之前的许多方法使用了半隐式的方案，即把隐式方案应用于粘性条件，显式方案应用于非线性对流条件，时间步通过CFL数控制。

Choi 和Moin在分步法的基础上采用了一个完全隐式的方法，首先对纳维斯托克斯方程在时间上离散，然后进行空间离散，用这种方法得到的中间速度分量是耦合的，后来使用牛顿迭代方法得到了中间速度分量。

为了防止迭代过程，Rosenfeld提出了一个非耦合的隐式解算器错误!未找到引用源。

，他设计了三个时间步的线性化方案，这个方案需要n-1步和n步的速度场来得到n+1步的速度，在不忽略时间二阶精度和稳定性的基础上，控制方程被解耦。

在最近的研究中，Kyoungyoun Kim，Seung-Jin Baek 和Hyung Jin Sung错误!未找到引用源。

对解决不可压缩湍流的纳维斯托克斯方程发展出了一种有效的数值计算方法，这个算法提出了一个新的隐式速度解耦过程，采用完全隐式的时间推进，在块LU分解和近似分解的基础上，速度项和压力项被解耦，同时保留了时间二阶精度，另外，由于隐式的对流条件，中间速度是耦合的，所以重点放在了对中间速度的解耦上，这就需要对第n个时间步的速度近似分解，这些解
耦过程同样保留了时间二阶精度。

本文数值模拟的过程中也用到了这种解耦过程，第二部分将会对目前的数值解耦方法及方程的近似分解过程做一个简要的介绍，在第三部分，将会把这个解耦过程应用于槽道流，并用直接数值模拟进行验证，画出结果进行比较分析。

2.2 数值方法
不可压缩的无量纲纳维斯托克斯方程为：
1,(1,2,3)Re i i i j j i j j
u u p u u i t x x x x ∂∂∂∂∂+=-+=∂∂∂∂∂ (1) 0i i
u x ∂=∂ (2) 其中，i x 是笛卡尔坐标，i u 是每一个方向相应的速度分量，Re 是雷诺数，所有
的分量都用特征长度进行了无量纲化。

在第n+1/2个时间步上，对上述两个方程进行空间和时间离散，方程可以写成如下的形式：
111/2111(()())()22Re n n n n n n n u u H u H u Gp Lu Lu mbc t ++++-++=-+++∆(3)
10n Du cbc +=+ (4)
其中，L 代表离散的拉普拉斯算子，H 代表离散对流算子，G 代表离散梯度算子，D 代表离散散度算子。

必须说明的是，边界条件已经被合并到方程(3)和(4)中，空间离散算子L 、H 、G 和D 定义在交错网格上，并且用二阶中心有限差分来定义，变量1n u +和1/2n p +定义在内部节点上，而不是交错网格上的边界点，边界上已知的速度分量和速度的边界条件被分别加在了mbc 和cbc 上。

在交错网格上，压力被定义在立方体中心，速度分量被定义在正交平面上，并且动量方程的离散化需要边界上的压力或压力梯度。

对于时间离散，使用了一个完全隐式的时间离散，因为对流条件和粘性条件是隐式的，所以产生了非线性的方程。

在最近的研究中，非线性条件被线性化，同时保留了时间二阶精度：
11112()n n n n n n n n i j i j i j i j u u u u u u u u O t ++++=+-+∆ (5)
通过这个线性化，对流项的线性算子N 可以被定义为：
111(()())2
n n n Nu H u H u ++=+ (6) 注意到线性化算子N 中包含第n 步的速度。

通过使用对流算子N ，离散方程(3)和(4)可以以矩阵的形式写出：
100n A G r mbc u D cbc p δ+⎛⎫⎛⎫⎛⎫⎛⎫=+ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭⎝⎭ (7)
11[()]2Re
A I t N L t =+∆-∆ 1/2112Re n n n r u Gp Lu t -=
-+∆ 1/21/2n n p p p δ+-=-
在上述公式中，通过反演方程(7)的系数矩阵，得到了下一个时间步上的值
1n u +和1/2n p +，因为(7)的系数矩阵大而且稀疏，所以方程(7)的求解过程很难，换句话说，因为在动量方程和连续性方程中，速度和压力是相关的，所以方程
(7)不能被直接求解。

对方程(7)使用块LU 分解：
1000n A I tG r mbc u D tDG I cbc p δ+⎛⎫⎛⎫⎛⎫⎛⎫⎛⎫∆=+ ⎪ ⎪ ⎪ ⎪ ⎪-∆⎝⎭⎝⎭⎝⎭⎝⎭
⎝⎭ (8) 这个方程和(7)是不同的，作为近似因子分解结果，把G p δ近似成了tAG p δ∆。

在最近的研究中，压力p 通过p δ来表达。

上面近似分解的误差为：
2()0tMG p O t δ⎛⎫∆∆=
⎪⎝⎭，其中12Re M N L =- (9) 方程(8)也可以被写为： *00A r mbc u D tDG cbc p δ⎛⎫⎛⎫⎛⎫⎛⎫=+ ⎪ ⎪ ⎪ ⎪-∆⎝⎭⎝⎭⎝⎭
⎝⎭ (10) 1*0n I tG u u I p p δδ+⎛⎫⎛⎫⎛⎫∆= ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭
(11) 其中，*u 是1n u +的中间速度，通过以下操作，可以得到更加简化的表达式：
*Au r mbc =+ (12)
*tDG p Du cbc δ∆=- (13)
1*n u u tG p δ+=-∆ (14)
1/21/2n n p p p δ+-=+ (15) 因为在方程(3)和(4)中已经应用了边界条件，所以不需要对边界条件进行其他特殊处理，并且在方程(12)—(15)解耦速度和压力的过程中，同样保留了时间二阶精度。

然后，通过使用*u δ，将近似分解用于速度分量*
u 上，所以方程(12)可以被写为： *n A u Au r mbc δ=-++ (16)
**n u u u δ=- (17)
方程(16)可以以矩阵形式给出：
*11112131*21222322*313233331u I tM tM tM R tM I tM tM u R t tM tM I tM R u δδδ⎛⎫⎛⎫⎛⎫+∆∆∆ ⎪ ⎪ ⎪∆+∆∆= ⎪ ⎪ ⎪∆ ⎪ ⎪ ⎪ ⎪∆∆+∆⎝⎭⎝⎭⎝⎭ (18)
当动量方程通过半隐式方法离散时，方程(18)中的非对角线子阵,i j M （i j ≠
）为零，然而，在最近的全隐式方法中，,()i j M i j ≠
不再为零，因为非线性的对流条件为隐式条件，这就意味着*1u δ、*2u δ和*3u δ是完全耦合的。

通过对系数矩阵(18)进行近似分解，中间速度
分量可以用第n 个时间步的速度进行解耦，在近似分解的过程中，时间二阶精度被保留了。

*11112131*21222322*333132330010000u I tM tM tM R I tM I tM I tM u R t I tM tM tM I R u δδδ⎛⎫⎛⎫⎛⎫+∆∆∆⎫⎛ ⎪ ⎪ ⎪⎪ ∆+∆∆= ⎪ ⎪ ⎪⎪ ∆ ⎪ ⎪ ⎪⎪ +∆ ⎪∆∆⎭⎝⎝⎭⎝⎭
⎝⎭ (19) 像之前所描述的，当前的速度解耦过程只需要第n 个时间步上的速度，在三
个时间步的方案中，第n-1步和第n 步的速度都需要才能保证时间的二阶精度。

方程(19)的二阶误差项如下：
**11122111332***211222113322233***311223113332233()tM M u tM M u O t tM M u tM M u tM M u tM M u tM M u tM M u δδδδδδδδ⎛⎫∆+∆ ⎪∆=∆+∆+∆ ⎪ ⎪ ⎪∆+∆+∆⎝⎭
引入新的变量**1u δ、**2u δ后，*u 可以通过下面的式子得到：
**11111()I tM u R t
δ+∆=∆ (21) ****22222111()I tM u R M u t
δδ+∆=-∆ (22) ******33333113221()I tM u R M u M u t
δδδ+∆=--∆ (23) ****22233
u u tM u δδδ=-∆ (24) *****11122133u u tM u tM u δδδδ=-∆-∆
(25) **,(1,2,3)n i i i u u u i δ=+= (26)
由以上的式子可以看到，不用求解方程(12)中的大矩阵，就可以得到中间速度。

小结：①首先用(21)—(26)的速度解耦过程得出*u ；
②从方程(13)中解出p δ；
③从方程(14)计算得出下一步的速度1n u +。

2.3 直接数值模拟的验证
使现在的解耦方法保留时间二阶精度是很重要的，首先将2.2中介绍的程序应用到三维周期性的槽道流中，控制方程为不可压缩的纳维斯托克斯方程，边界条件为壁面无滑移条件，建立针对控制方程的数值离散化方法，本程序中使用了有限差分法错误!未找到引用源。

：Crank-Nicolson method ，在解耦过程中保持了恒定的质量流率，雷诺数的值是4400，计算域为:02,02,0X Y Z ≤≤∏≤≤≤≤∏,为了使计算结果更加准确，所取的网格点为：3128，0.01t ∆=。

下图为得到的壁面应力随时间的变化情况：
图2.1 直接数值模拟壁面应力图
图 2.1为用直接数值模拟计算的壁面应力随时间的演化过程。

由图可以看出，在无量纲时间400之后，数据趋于稳定，所以在计算流向平均速度、脉动和雷诺应力时，从40000步开始取数据，接下来画出流向速度随壁面高度的变化，在画图过程中不仅在同一水平高度进行了平均，而且将稳定后的各个时间段上的值做了平均错误!未找到引用源。

，因为展向速度和法向速度为零，所以没有画出，横纵坐标都经过了无量纲化，下图中实线表示实际模拟计算得到的结果，而虚线代表Lee 和Moser 错误!未找到引用源。

得到的结果：
图2.2 直接数值模拟流向速度图
图2.2为用直接数值模拟计算的流向速度随时间的演化过程。

可以看出，计算得到的结果与正确结果基本重合，可见，在直接数值模拟中流向速度得到了很好的验证。

下图为脉动''u u <>、''v v <>、''w w <>和雷诺应力''u v <>随壁面高度的变化，在平均上与流向速度的处理方式一样，实线代表实际模拟计算得
到的结果：
图2.3 直接数值模拟流向脉动图图2.4 直接数值模拟法向脉动图
图2.5 直接数值模拟展向脉动图图2.6 直接数值模拟雷诺应力图
图2.3、图2.4、图2.5、图2.6分别为用直接数值模拟计算的流向、法向、展向脉动和雷诺应力随时间的演化过程。

由上面四幅图同样可以看出，计算得到的结果与Lee和Moser的结果基本重合，所以，在直接数值模拟中脉动和雷诺应力也得到了很好的验证。

从以上关于流向速度、脉动和雷诺应力的图中可归纳出，现在的程序应用于直接数值模拟是有效的，并且在计算过程中得到了比较好的结果。

第三章 大气边界层在高雷诺数下的模拟结果
3.1 不同亚格子模型对大气边界层数值模拟结果的影响
在不可压缩湍流中，亚格子涡粘模型采用分子粘性的形式错误!未找到引用源。

，即：
123
ij ij t ij kk S τνδτ-=+
(1) 在公式(1)中，t ν为亚格子涡粘系数， _ij S 是可解尺度的变形率张量。

可以看出，在涡粘模型中，加上涡粘系数后，就可以使用纳维斯托克斯方程直接数值计算，其中，涡粘系数是需要对其封闭的，可以有不同的形式。

将亚格子涡粘模型(1)代入到大涡数值模拟的控制方程中，可以得到大涡数值模拟控制方程错误!未找到引用源。

：
_()()3i j j i kk j i j j i
u u u u u p t x x x x x τνν---
--∂∂∂∂∂∂+=-+++∂∂∂∂∂∂ ， 1Re t νν=+ (2) _0i i
u x ∂=∂ (3) 可以看出，代入亚格子涡粘模型和涡扩散模型后，湍流控制方程的变化只是在分子扩散系数上加亚格子涡粘系数，下面给出亚格子涡粘系数的两种形式。

3.1.1 Smagorinsky 模式
在大涡模拟中将雷诺平均混合长度模式推广，得到原始的Smagorinsky 模式，混合长度模式的涡粘公式错误!未找到引用源。

为：
'2
t u u l l y ν∂<>∝∝∂ 令l =∆，那么在二维空间中有：
2()2ij ij u S S y
∂<>=<><>∂ 将上式推广到三维空间中可得到如下关系： 1
2
2(2)t ij ij l S S ν∝<><>
同样，令l =∆，将平均运算改为过滤运算，并引入模型常数后，可将亚格子涡粘系数写成以下形式：
__
222t S ij ij C S S ν=∆<> (4)
在实际模拟计算中发现模型常数明显偏大,Smagorinsky 模式的致命缺点是耗散比较大，下面利用近壁阻尼函数，对Smagorinsky 系数进行修正，从而减少耗散错误!未找到引用源。

对Smagorinsky 系数S C 做如下修正：
0.18[1exp(/)]s C y A ++=--
式中/y yu τν+=，26A +=，[1exp(/)]f y A ++=--称为近壁阻尼函数。

3.1.2 动力模式
涡粘错误!未找到引用源。

和涡扩散模式在工程计算中有普遍的应用。

但因其导出过程以均匀各向同性湍流为基础，实际流动中的复杂湍流并不满足这些条件，所以通过对亚格子模型的改进，使这种模式适应局部湍流的结构，下面介绍的动力模式是经过改进的亚格子涡粘模型。

动力模式在模拟过程中使用了两次过滤，假定过滤过程是线性的，且有12∆>∆，由1∆、2∆过滤产生的亚格子应力和连续两次过滤产生的亚格子应力
可得到如下关系：
13ij kk ij D ij L L C M δ-=
因为上式为超定方程，所以不能直接计算系数D C ,但有以下几种方法来解决超定
问题错误!未找到引用源。

第一个为变形率张量收缩法，即将上式两边同乘以可解尺度的变形率张量。

但是实际计算表明，用这种方法计算的模式系数很不规则，而且计算过程中的稳定性较差；第二种为最小误差法，令这个超定方程两边的平方差最小，即：
2
103ij kk D ij D L L C M C ⎧⎫∂--=⎨⎬∂⎩⎭ 可得到：
ij ij D ij ij M L C M M = (5)
在式(5)中，ij L 和流动的形态有关，ij M 和模式的形式有关。

最小误差法比变形
率张量收缩法已经有很大改进，但是仍有缺陷：①模式系数D C 可能会出现负值，
导致计算结果发散；②分母可能很小，也导致计算结果发散。

为了克服计算上的
困难，采用平均系数法，对(5)式右端的分子和分母分别求系综平均，得到模式系数为：
ij ij
D
ij ij
M L
C
M M
<>
=
<>
(6)
3.1.3 大涡模拟的验证
将大涡模拟应用到低雷诺数的三维周期性槽道流中，使用的亚格子模型为动力模式，雷诺数的值为4400，因为大涡模拟不直接计算近壁区域，所以增大了网格的分辨率错误!未找到引用源。

，模拟过程中使用的网格点数为：3
65，计算域和直接数值模拟相同:02,02,0
X Y Z
≤≤∏≤≤≤≤∏，0.01
t
∆=。

下图为得到的壁面应力随时间的变化情况：
图3.1 大涡模拟壁面应力图
图3.1为用大涡模拟计算的壁面应力随时间的演化过程。

可见，在无量纲时间200之后，数据趋于稳定，所以在计算平均速度以及脉动和雷诺应力时，从20000步开始取数据，接下来画出流向速度随壁面高度的变化，处理方法同直接数值模拟：
图3.2 大涡模拟流向速度图
图3.2为用大涡模拟计算的流向速度随时间的演化过程。

因为未加壁函数，所以第一个网格点在近壁区域，大涡模拟取得网格点比直接数值模拟少，得到的结果没有直接数值模拟精确，计算结果与正确的结果之间出现了一定的偏差，但是，因为大涡模拟所取的网格点比较稀疏，加之壁面对结果的影响比较大，可看出，得到的结果是在误差范围内的。

下图为脉动''u u <>、''v v <>、''w w <>
和雷诺应力''u v <>随壁面高度的变化，在平均上与速度的处理方式一样，先对
同一个时间点同一高度上的值进行平均，然后在各个时间上平均，实线代表模拟计算得到的结果：
图3.3 大涡模拟流向脉动图 图3.4 大涡模拟法向脉动图
图3.5 大涡模拟展向脉动图 图3.6 大涡模拟雷诺应力图 图3.3、图3.4、图3.5、图3.6分别为用大涡模拟计算的流向、法向、展向脉动和雷诺应力随时间的演化过程。

所以，在大涡模拟中，速度脉动和雷诺应力在误差范围内也得到了相对较好的验证。

同时也证明了现有程序对大涡模拟的有效性。

3.1.4 两种亚格子模型的模拟结果
在3.1.1和3.1.2中分别介绍了Smagorinsky 模式和动力模式，并且利用近壁阻尼函数对Smagorinsky 模式进行了改进，动力模式采用涡粘模式表达式，改进之后的模型在减少计算量的基础上有效的模拟了壁面对湍流边界层的影响，在模拟湍流边界层的过程中雷诺数增加到了15000，网格数为365，计算域不变 ，0.005t ∆=，缩短了时间间隔，首先给出两种亚格子模型下流向速度随壁面高度的变化趋势，并将结果进行对比：
图3.7 不同亚格子模型模拟流向速度图
图3.7为用不同亚格子模型计算的流向速度随时间的演化过程。

用两种亚格子模型得到的结果与Lee和Monser错误!未找到引用源。

得到的结果相比，有一定的偏差，但这个偏差是在误差范围内的，所以得出：在不加壁函数时，这两种亚格子模型对湍流边界层的模拟是有效的。

下面给出速度脉动和雷诺应力随壁面高度的变化：
图3.8 不同亚格子模型模拟流向脉动图图3.9 不同亚格子模型模拟法向脉动图
图3.10 不同亚格子模型模拟展向脉动图图3.11 不同亚格子模型模拟雷诺应力图
图3.8、图3.9、图3.10、图3.11分别为用不同亚格子模型计算的流向、法向、展向脉动和雷诺应力随时间的演化过程。

从以上几幅图中可以看出用两种亚格子模型模拟得到的结果非常接近，但与Lee和Moser的结果出现了一定的偏差，原因是：Lee和Moser是用直接数值模拟对整个槽道模拟得到的结果，而这两种亚格子模型是在半槽道中模拟，并且在湍流边界层的模拟中，强制假定上边界法向速度的边界条件为零，这些可能导致误差比较大，但是可以看出，偏差是。