边界层理论与壁面方法详解

1.      问题

 边界层处理是所有CFD工程师都面临的一个难题，

 第一层取多厚？

    y+值是多少？

 做多少层？

 层与层间的增长比例？

 边界层总厚度多少？

 与试验值不符合，要不要再加密一下……种种问题不一而足。

 边界层理论在所有流体力学著作里都有涉及，最专业的著作当属德国航空专家H. Schlichting的《Boundary-Layer Theory》，这本专著已经出到第八版了，足见其影响力，目前亚马逊有售。本文主要内容都是笔者的工程经验，理论部分主要取自上述专著，难免有疏漏或不当之处，望方家斧正。

2.      理论

 边界层的概念是1904年德国著名的力学家普朗特在海德尔堡第三届国际数学家学会上宣读的“关于摩擦极小的流体运动”的论文中首先提出的。他根据理论研究和实际观察，证实了对于水和空气等粘性系数很小的流体，在大雷诺数下绕物体流动时，粘性对流动的影响仅限于紧贴物体壁面的薄层中，而在这一薄层外粘性的影响很小，完全可以忽略不计。普朗特把这薄层称为边界层，或称附面层。从边界层内的流动过渡到外部流动是渐变的，所以边界层的厚度δ通常定义为从物面到约等于99%的外部流动速度处的垂直距离，它随着离物体前缘的距离增加而增大。根据雷诺数的大小，边界层内的流动有层流与湍流两种形态。一般上游为层流边界层，下游从某处以后转变为湍流，且边界层急剧增厚。层流和湍流之间有一过渡区，称为转捩（Transition）现象，NASA AmesCenter曾在导弹风洞试验中复现这一现象，如图1所示。

图1 风洞试验中的层流-湍流转捩现象

图1很直观的给大家展现了层流到湍流的转捩过程，我们对湍流区域的边界层进一步放大，可以归纳出另一张经典的湍流边界层分区，如图2所示。

图2 湍流边界层分区

图2使用归一化坐标而非物理坐标y，这是因为，使用物理坐标时，随着雷诺数Re的变化，其速度值是变化的，而研究中更需要寻找出不变性，通过归一化处理后，消除了Re显示变化影响（吸收到摩擦速度和粘性长度的定义里），因此得到了普适的结果。

下面对图2的物理量进行解释，横坐标y+的定义为

纵坐标

 其中y是离壁面的距离，u_t可以理解为壁面附近流体的剪切速度，v是运动粘度， u代表主流速度。细心的读者应该可以注意到，这y+其实是个雷诺数的公式，所以y+是y处漩涡的典型雷诺数，也反映了粘性影响随y的变化。根据y+值的大小，可以将边界层内分割为三个区域：

1)     粘性底层（Viscous sublayer）：𝑦+ < 5，粘性底层内，u+ = y+，由于横坐标y+采用对数坐标，所以在图2上显示为曲线；

2)     过渡区（Buffer layer）：5 ≤ 𝑦+ ≤ 30；

3)     对数区（Log-law region）：30 < 𝑦+ < 𝑓(𝑅e)，对数区内u+ = 2.5*ln（y+）+ 5.45，在对数坐标下为一直线。该区域的上限即边界层总厚度，随Re变化，一般认为是99%的外部流动速度处，定义为δ_tot，对于平板流动，Schlichting给出了经验估算方法，

3.   STAR-CCM+壁面模型

根据壁面y+值及边界层网格分布，发展出了不同的壁面模型，如图3所示。

图3  不同y+壁面算法示意

STAR-CCM+壁面模型的特点是自动及全面，如图4所示，具体应用技巧如下：

1)     Low y+ walltreatment：y+值小于1, 一般推荐25+层，其中3+层位于粘性底层内，5+层位于过渡区，5+层位于对数区内，直接求解近壁面速度分布，主要用于外气动计算中的流动分离，转捩现象；

2)     High y+wall treatment：y+值位于30至150之间，5层左右边界层，主要用于无流动分离的情况；

3)     默认选项为All y+ wall treatment，顾名思义，不需要用户选择，可以自动根据来流速度，边界网格计算y+值，并自动切换壁面算法，使用时注意y+<5或30<y+<150，尽量避开过渡区。

图4 STAR-CCM+壁面模型

错误的试用壁面模型可能使计算结果严重失真，图5为一管道流动案例，分别采用y+=1（红色实线），y+=20（蓝色实线）获得的计算结果与DNS（Direct Numerical Simulation，可认为真值，图中圆点）结果对比，可以看到，在横坐标y+值小于200的区间内，y+=1的U+与软流粘性系数计算结果与DNS结果吻合度极高，而蓝色的y+=20位于过渡区内，现有的壁面模型处理该区间都有缺陷，因此在计算中应尽量规避。

图5 不同y+模型计算结果对比

4.      STAR-CCM+边界层网格处理

STAR-CCM+边界层网格生成方法非常简单，默认设置如图6所示，边界层定义方法采用Stretch factor，需要定义的参数有：Numberof Prism Layer （边界层总层数），Prism layer stretching（相邻层间的比例），Prism LayerTotal Thickness（边界层总厚度），边界第一层网格厚度通过上述三个参数计算得出。

图6 Stretch Factor方法

如果想直接定义底层厚度，只需要在Prism Layer Mesher的属性栏里更改生成方法，相应的参数变更为Prism Layer Near WallThickness（底层厚度），层间增长率通过计算得出。

图7 wallThickness方法

以NASA CRM飞机模型巡航状态工况为例（Ma=0.85，Re=5000000，特征弦长Lc=7m），若采用图7的wallThickness方法，需要确定三个主要参数：

1. Prism Layer Near WallThickness（底层网格厚度）

可通过下列公式计算y+=1条件下的底层网格厚度，为0.037mm

2. Prism Layer TotalThickness（边界层总厚度）

利用

计算得出

3. Number of PrismLayer （边界层总层数）

外气动计算一般推荐层间增长率在1.2~1.3之间，此处取增长率1.2，可计算得到大约需要35层边界层。

文章最后展示一个很常见的误区，许多工程师熟悉边界层理论，也按照要求做好了网格，计算结果却不理想，很多时候都是由图8所示的问题引起，忽视了边界层最外层与主流网格的过渡，导致从一个很小的边界层网格直接过渡到主流区的大网格，体积变化率（Volume Change，衡量网格质量标准之一）可达上千。正确的做法是右图所示，通过计算，设置边界层总厚度为面最小网格尺寸1~2倍，即可获得适当的过渡比例。

图8 边界层与主流网格过渡

免责声明：本文系网络转载或改编，未找到原创作者，版权归原作者所有。如涉及版权，请联系删