大气边界层在高雷诺数下的模拟结果

大气边界层在高雷诺数下的模拟结果
3.1 不同亚格子模型对大气边界层数值模拟结果的影响
在不可压缩湍流中，亚格子涡粘模型采用分子粘性的形式[19]，即：
123ij ij t ij kk
S τνδτ-=+
(1) 在公式(1)中，t ν为亚格子涡粘系数， _ij S 是可解尺度的变形率张量。

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

将亚格子涡粘模型(1)代入到大涡数值模拟的控制方程中，可以得到大涡数值模拟控制方程[12]：
_()()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 模式，混合长度模式的涡粘公式[20]为：
'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 系数进行修正，从而减少耗散[19]。

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

3.1.2 动力模式
涡粘[13]和涡扩散模式在工程计算中有普遍的应用。

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

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

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

但是实际计算表明，用这种方法计算的模式系数很不规则，而且计算过程中的稳定性较差；第二种为最小误差法，令这个超定方程两边的平方差最小，即：
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，因为大涡模拟不直接计算近壁区域，所以增大了网格的分辨率[15]，模拟过程中使用的网格点数为：365，计算域和直接数值模拟相同:02,02,0X Y Z ≤≤∏≤≤≤≤∏，0.01t ∆=。

下图为得到的壁面应力随时间的变化情况：
图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[16]得到的结果相比，有一定的偏差，但这个偏差是在误差范围内的，所以得出：在不加壁函数时，这两种亚格子模型对湍流边界层的模拟是有效的。

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

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

3.2 不同壁函数对大气边界层数值模拟结果的影响
3.2.1 本文使用的壁函数
自然界中的流动一般为高雷诺数的流动，模拟时大约70%的网格点位于近壁区域中[21]，所以，如果要完全模拟近壁区域的流动，所需要的网格点非常多，这就失去了大涡数值模拟的优越性。

为了保持大涡模拟的优点又不增加计算量，近壁处常要用壁模型进行近似，在模拟过程中要使第一层网格的坐标位于对数区，从而在第一层网格上构造壁面亚格子模型[17]。

本文首先应用Schumann 提出的壁模型来模拟湍流边界层，并将结果推广到大气边界层的模拟上，要用到的公式在文章的第一章已经做了简要的介绍，其基本思想是假定近壁区存在当地可解尺度速度分布和当地剪应力之间的平衡关系式，第一层网格的垂直距离1y +>30，视流动的雷诺数而定，雷诺数越高，取较大
的1y +值，用流向和展向的可解速度得到相应方向的亚格子应力。

在模拟湍流边界层的过程中也用到了Werner 提出的模型，Werner 建议用幂函数作为平衡层的速度分布，具体公式如下：
2121211112/,0.5/11()(),0.5/2B P P P P w B B B B B B P P P
P P u u A B B A u u A A μδμδτμμμδδδ-+++--⎧≤⎪⎪=⎨⎡⎤-+⎪+>⎢⎥⎪⎣⎦⎩ 式中，P u 为壁面附近第一个网格点P 上与壁面相切的速度，P δ为P 与壁面的垂
直距离，A ，B 为常数，在本文的数值模拟中取A ，B 的值为：A=8.3，1/7B =。

3.2.2 模拟结果
本节为两种不同壁函数对湍流边界层数值模拟结果的影响，在数值模拟的过程中，为模拟高雷诺数下的边界层流动，使雷诺数增加到了187500，因为程序中加了壁函数，所以法向网格在降低计算量和保留计算精度的条件下调整为35，而保持流向网格和展向网格数为65不变，模拟的计算域与直接数值模拟相同:02,02,0X Y Z ≤≤∏≤≤≤≤∏ ，0.01t ∆=，因为是大涡模拟，所
以在计算过程中要注意使第一层网格点130y +>，下图画出在Schumann 壁面模
型和Werner 壁面模型下流向速度随壁面高度的变化，并且与Lee 和Moser 的模拟结果进行了对比：
图3.12 不同壁函数模拟流向速度图
图3.12为用不同壁函数计算的流向速度随时间的演化过程。

从结果中可以看出，模拟的第一层网格点在对数区。

Schumann 壁面模型与Werner 壁面模型模拟得到的结果是不一样的，有一定的偏差，将湍流边界层的模拟结果与Lee 和Moser 相比，Schumann 壁面模型模拟的结果已经很精确了，而Werner 壁面模型
画出的结果偏差很大。

下面画出速度脉动''u u <>、''v v <>、''w w <>和雷诺
应力''u v <>随壁面高度变化的廓线：
图3.13 不同壁函数模拟流向脉动图 图3.14 不同壁函数模拟法向脉动图
图3.15 不同壁函数模拟展向脉动图 图3.16 不同壁函数模拟雷诺应力图 图3.13、图3.14、图3.15、图3.16分别为用不同壁函数计算的流向、法向、展向脉动和雷诺应力随时间的演化过程。

可看出，使用两种壁函数模拟得到的脉动''u u <>、''v v <>、''w w <>和雷诺应力''u v <>在误差范围内是基本吻合的，但是与Lee 和Moser 的结果相比，两种壁函数都有一定的不准确性，这种不准确是由亚格子模型的耗散和壁函数引起的。

因为Werner 壁面模型模拟湍流边界层的结果并不是很好，分析原因可能是因为法向网格数太少引起的，故将法向网格数调整为65和95，保持流向网格和展向网格数不变，向上平移坐标，使第一个网格点在对数区，分析不同网格对模拟结果的影响。

下图为Werner 模式下不同网格引起流向速度的变化：
图3.17 不同网格下的流向速度图
图3.17为Werner 壁面模型在不同网格下计算的流向速度随时间的演化过
程。

可见，增加网格使结果有所改善，但改善的幅度并不是很大，模型的不准确
性还是很明显。

下面为速度脉动''u u <>、''v v <>、''w w <>和雷诺应力
''u v <>在不同网格下的变化：
图3.18 不同网格下的流向脉动图 图3.19 不同网格下的法向脉动图
图3.20 不同网格下的展向脉动图 图3.21 不同网格下的雷诺应力图 图3.18、图3.19、图3.20、图3.21分别为Werner 壁面模型在不同网格下计算的流向、法向、展向脉动和雷诺应力随时间的演化过程。

从结果中可以得出这样的结论：网格点并不是引起结果不准确的关键因素，主要原因还是亚格子模型和壁函数。

所以，接下来的主要工作就是发展更好的亚格子模型和壁函数。