让边界层的求解飞起来吧

LBM与流体力学 2020年6月5日浏览：2582 评论：1

上一篇的结尾答应大家要在这一篇里面介绍壁面函数。

一说到函数，估计又有读者要吐槽，你们是不是又打算列一堆公式，然后告诉我们，这就是壁面函数，我到现在连拉普拉斯算子和哈密顿算子还没搞清楚呢。额，好吧，我承诺，这篇文章不放任何的公式，给大家一起慢慢白话边界层和壁面函数。只是，大家不要嫌弃我们太啰嗦。

那我们就从边界层的定义说起吧。

边界层的概念

通常情况下，对于牛顿流体的流动，可以认为N-S方程精确描述了运动的规律。然而迄今为止，对于N-S方程，人们只得到了几种简单流动的解，一般流动情况都是解不出来的。大家知道，我们平时接触到的大部分流动的雷诺数都是很大的，也就是说粘性力的影响应该很小。如果忽略流体的粘性，则N-S方程就变成了欧拉方程。按照这种想法，通过欧拉方程直接求解高雷诺数的流动应该没有问题。可是，通过欧拉方程得出的解和实际流动经常完全不同。比如大家现在看到的这个圆柱绕流的例子，忽略粘性的结果和实际流动大相径庭。

那么问题到底出在哪儿了呢？难道雷诺数不是代表惯性力与粘性力之比？

作为现代流体力学之父，咱们广义的祖师爷普朗特在1904年提出了边界层的概念，并系统性的回答了刚才的问题。普朗特通过观察和推导得出：物体表面存在一个很薄的剪切层，即边界层。边界层之外的流动基本可以看做是无粘的，而紧挨着固体表面的流体分子被吸附在表面上。于是，在固体表面附近，存在着一个法向速度梯度很大的区域，这就是边界层。对于高雷诺数的流动，边界层是很薄的，壁面附近的区域速度梯度很大，与速度梯度成正比的粘性力也非常大，这就是粘性对于高雷诺数流动不可忽略的原因。

边界层的类型

大家知道，流体的流动状态会伴随着雷诺数的变化而出现明显的差异，主要包含层流和湍流两种状态。那么，相应的边界层也会呈现两种完全不同的状态，借用NASA的官图，不同于层流边界层的一目了然，湍流边界层内部是复杂的 瞬态流动，其速度分布的影响因素除了粘性力外，还有速度脉动带来的雷诺应力。我们在实际工作中碰到的绝大部分工程问题都是湍流，也是 CFD软件的主攻方向。如果不激活软件中的转捩模型，默认计算的都是全湍流——而大部分工程问题也是可以忽略转捩的影响。所以我们本节讨论的内容限定为湍流边界层。

湍流边界层的直观认识

尽管湍流及其边界层内的速度呈现明显的脉动特性，但是如果用平均速度的观点来看湍流边界层的话，大致符合下图所示的速度分布。湍流边界层内的速度分布极其复杂，可以按照受壁面影响的程度为内层和外层。 内层又包含线性底层、过渡层和对数层（因为这里的速度分布满足对数率），其中 线性底层+过渡层又合称为粘性底层，这部分区域的粘性力影响很大，尤其是线性底层的流动基本是层流。而湍流边界层的外层，主要受主流区的影响。由于流动的非定常性，外层和内层之间没有固定的界限，所以对数层也称交叠层，同时属于边界层的内层和外层。

另外上图给出了上述各层对应的Y+的值（为了兑现本文不出现公式的承诺，这儿我就假装大家都知道Y+和U+的概念啦）。其中有几个非常重要的数值，大家可以记住，通常情况下，线性底层位于Y+小于5的区域，而对数层则大致位于Y+为30到500之间的区域。

湍流边界层内湍动能的传播

接下来从物理层面理解一下近壁区内湍动能的生成、传播和耗散。下图显示了湍流边界层内的流动特征。最下面的层为粘性底层，这里发生强烈的剪切， 是湍动能的主要生成区；中间为对数层， 湍动能的生成约等于耗散；最外层则大致为湍流边界层的尾流层，这里的流动已经比较接近之前介绍VLES的文章中介绍的Cascade， 该区域的湍动能主要体现为耗散，生成率很小。由此可知，在湍流边界层内部，湍动能由内向外逐步传递，由小涡传递给大涡，称之为能量的反向传递。

让边界层的求解飞起来吧的图4

看到这里，相信大家对于边界层已经有了一个比较清晰的认识了。那么我们在实际的CFD计算里面如何进行求解呢？

01 近壁区的限制

对于大家熟知的RANS方法中的k-epsilon模型，是针对充分发展的湍流才是有效的，只能用于求解湍流核心区域的流动，即较高雷诺数下的流动。而对于边界层，由于壁面的阻碍作用，沿壁面法向的流态变化很大。因此湍流边界层无法用高雷诺数的RANS模型直接求解。

02 如何解决

方法通常有二：

要么在近壁区修改湍流模型，使用低雷诺数的湍流模型来求解粘性底层和对数层；比如叶轮机械领域常用的SST模型，它结合了自由流中的 k-epsilon 和靠近壁面处的 k-omega 模型，要求网格划分到粘性底层，其Y+一般要求小于5；

要么使用半经验公式，也就是壁面函数法，将近壁区的物理量和湍流核心区的物理量直接联系起来。壁面的网格只需画在壁面函数起效的区域，也就是上图中的对数层，其外侧边界Y+约为500这个量级，对于近壁面网格也有相应的要求。当然有些改进的壁面法则已对Y+不那么敏感。

两种方法的优劣：

如果采用第一种方法，给定壁面无滑移条件，即可直接求解，精度高但计算量大；如果采用壁面函数法，由于使用的是一套半经验公式，要求流动处于准平衡态，不适用于流动容易失稳的区域，比如逆压梯度较大的区域或者大分离区，如下图翼型上侧的吸力面区域。

03 对壁面函数的改进

基于上述标准壁面函数的问题，很多软件推出了相应的解决办法。某些RANS软件使用的“非平衡壁面函数”就是其中之一，在标准壁面函数的基础上增加了压力梯度修正项，把应用范围拓宽到分离、再附着等复杂的问题上。

04 RANS壁面函数的实施

在RANS求解器中，计算过程通常是这样的：赏花归去马如飞——去马如飞酒力微——酒力微醒时已暮——醒时已暮赏花归，通过不断迭代获得结果。噢Sorry，计算过程其实是这样的：根据初始流场计算出壁面Y+值——套用壁面函数求得壁面剪切力——计算壁面附近的湍动能生成和耗散——求解湍动能和湍流耗散率的输运方程——更新全场的湍流度——回到起始的步骤，如此迭代求解。遇到较大的逆压梯度，则需在壁面函数中加入压力修正项。

05 PowerFLOW的不同之处

PowerFLOW的高级壁面模型（ABLM）虽然也是考虑压力修正的壁面函数，但由于LBM方法本身是显式的时间推进算法，每个时间步内并没有迭代，所以从第一步的初始流场出发，每时每刻根据壁面外侧的速度和壁面函数联立实时求得壁面剪切力、壁面滑移速度等参数，整个过程也是时间推进的。对于使用了壁面函数的计算而言，最重要的就是准确计算考虑边界层效应的Uslip，尽管边界层并没有直接求解。这句话是不是听上去有点拗口，那就再读一遍。LBM方法的巧妙之处就在于尽管使用了壁面函数，粒子与壁面的碰撞还是能准确的计算，进而提供一个更准确的Uslip。于是ABLM借助LBM+VLES再度自由翱翔。