在工程应用领域，两方程湍流模型仍然有很大的优势[1]．但是，两方程湍流模型在推导中有诸多假设，所以仍存在一定的不足，一直处于不断发展的过程中．目前，两方程湍流模型对于边界上湍流变量的研究还不是很完善，许多计算流体力学代码都指定默认值，并允许用户在需要时进行调整，但没有统一的标准做法．通常入口边界条件的选择对数值模拟结果的影响很大[2-4]，但很少有研究给出不同选择是如何影响模拟结果的分析，也很少有人会注意入口湍流参数及湍流模型封闭系数的选择．有研究指出k-ε湍流模型在流动模拟中还会出现湍流衰减现象[5-7]，即当采用k-ε湍流模型进行数值模拟时，入口边界给定的湍流参数初值会随着下游流动逐渐衰减，当到达研究对象附近时，当地湍流参数的值已经不再是入口边界的初值，导致模拟对象周围的湍流和真实的物理环境中会不一致，进而影响模拟结果的准确性．真实大气边界层中的风廓线都是带有剪切特性的，研究不同入口条件下两方程湍流模型在计算剪切流时湍流动能衰减的规律及原因对于研究大气边界层中风力机周围及风电场内的流动特性有重要意义．这里以中国规范B类风场平均风剖面为入口条件，基于k-ε模型研究了剪切流中湍流动能的衰减规律．首先，通过Matlab ode45函数求解剪切流条件下的湍流变量的输运方程；然后，通过fluent数值模拟，并与Matlab理论数值解、均匀流数值模拟解对比，分析在剪切流中湍流模型、入口边界条件、入流速度梯度和输运方程扩散项对湍流变量衰减的影响规律．1 理论推导方法1.1　标准k-ε湍流模型标准k-ε湍流模型[8]的湍动能k和耗散率ε的输运方程为∂(ρk)/∂t+∂(ρkui)/∂xi=∂∂xjμ+μtσk∂k∂xj+Pk-ρε; (1)∂∂t(ρε)+∂∂xi(ρεui)=∂∂xjμ+μtσε∂ε∂xj+C1εεPk/k-C2ερε2/k, (2)式中：ρ为流体密度；ui为流速分量；t为时间变量；xi和xj为流向分量；μ为分子黏性；μt=ρCμk2/ε为涡黏性系数；Pk=μtS2为由平均速度梯度产生的湍流动能；S=2SijSij为平均应变率张量的模[9]，其中，应变率张量Sij=(∂ui/∂uj+∂uj/∂ui)/2，uj为流速分量；Cμ，C1ε，C2ε，σk，σε为一组模型常数，其默认值为Cμ=0.09，C1ε=1.44， C2ε=1.92，σk =1，σε=1.3．1.2　平均风剖面在大气边界层内常用的风速剖面模型主要有对数律和指数律两种，而对数律模型的适用范围是距地面100 m[10]以下，大多数学者进行数值计算时也都推荐使用指数律[11-13]．采用指数律模型表征平均风剖面，其表达式为u¯(z)/u¯b=(z/zb)a，式中：zb和u¯b分别为参考高度和在该高度处的平均风速；z和u¯(z)分别为任一高度和该高度处的平均风速；ɑ为地面粗糙度指数，随地形地貌而变化．1.3　理论推导过程在剪切流中存在速度梯度，所以方程(1)和(2)的右侧湍流动能生成项不为0；另外，假设不存在湍流变量的梯度(实际的剪切流动中，湍流变量在垂直方向上一般也都有切变，这里为简化模型，以便理论值和数值模拟结果对比，将其假设为常数)，即入口湍流变量不随高度变化，则方程(1)、(2)右侧的扩散项也为0．最后，湍流变量输运方程右侧只剩下生成项和耗散项：u∂k/∂x=Pk/ρ-ε；(3)u∂ε/∂x=C1ε(ε/k)(Pk/ρ)-C2εε2/k．(4)假设在剪切流中，平均风剖面满足中国规范B类风场[14]要求，流向为单一的x方向，则风速为u¯(z)=kuu(z/zG)α，(5)式中：ku为常数，取1.77；u为参考风速，取10 m/s；zG梯度风高度的参考值，取350 m；ɑ为地面粗糙度指数，取0.16．此时，直角坐标系中的应变率张量可写成：Sij=1200∂u/∂z000∂u/∂z00；S≡2SijSij=(∂u/∂z)2；Pk=μtS2=μt(∂u/∂z)2．由于方程(3)和(4)无法求得解析解，因此采用龙格-库塔算法求数值解，通过调用Matlab中ode45函数进行求解．湍动能和耗散率的初始值采用表1中的7种入口条件，垂直高度z取10和100 m作为典型高度分析，表中Iin，μin/μ，kin，εin，μin分别为入口湍流强度、湍流黏性比、湍动能、耗散率和湍流黏性的初始值．湍动能k和耗散率ε分别为k=(3/2)[(I/100)u]2；(6)ε=ρCμ(k2/μ)(μt/μ)-1，(7)式中I为湍流强度．压强为101.325 kPa、温度为20 ℃时，空气的动力黏度和运动黏度分别为μ=17.9×10-6 Pa·s，ν=14.8×10-6 m2/s．10.13245/j.hust.210915.T001表1不同入口边界条件下湍流变量的初始值入口边界条件Iin/%μin/μkin/(m2·s-2)εin/(m2·s-3)μin/(kg·m-1·s-1)10.0010.011.500×10-81.385 82×10-101.79×10-720.0100.011.500×10-61.385 82×10-61.79×10-730.1000.011.500×10-40.013 861.79×10-741.0000.010.015138.582 401.79×10-7510.0001.001.5001.385 82×1041.79×10-5610.00010.001.5001 385.824 021.79×10-4710.000100.001.500138.582 401.79×10-32 数值模拟方法2.1　计算域和网格计算域大小为80 m×200 m，采用结构网格．在剪切流中，近壁面须要进行一定加密．当采用标准k-ε湍流模型时，须保证y＋在30~300范围内，采用标准壁面函数时，应尽量接近30[15]．假定y+值取30，则首层网格高度为1.46×10-3 m．在均匀流条件下，x和y方向网格均匀分布．采用表1中第7种入口条件完成网格无关性验证．剪切流中，垂直方向(y方向)上的网格增长率取1.05，同时验证x方向首层网格高度对不同位置处湍动能衰减的影响．经验证，剪切流条件下，x方向首层网格高度采用0.01 m，y方向首层网格高度采用通过y+计算所得值，总网格数为400×1 000时计算结果比较可靠，计算域和网格如图1．在均匀流条件下，当网格数取400×1 000时计算结果也比较可靠，相应的计算域和网格如图2所示．10.13245/j.hust.210915.F001图1剪切流条件下的计算域和网格10.13245/j.hust.210915.F002图2均匀流条件下的计算域和网格2.2　数值方法采用fluent完成对方程(1)和(2)的求解，并与上文中Matlab的理论数值解进行对比．求解流动方程时，空间离散采用中心节点的有限体积格式，基于SIMPLEC算法求解稳态纳维-斯托克斯方程．对流项离散使用QUICK 格式，其余各项采用二阶迎风格式离散．残差收敛至1×10-6，迭代至收敛．2.3　边界条件入口边界条件取为速度入口，入口速度按式(5)计算，湍流动能和耗散率的初始值分别取表1中的7种入口条件；出口边界条件取为压力出口边界，压力值等于大气压，湍流参数值同入口；顶部边界条件取为对称边界；底部边界条件取为无滑移壁面，粗糙度高度取为零(在真实的剪切流动中，存在粗糙度，此处为了与理论值对比，忽略粗糙度的影响)．均匀流条件下，除入口速度外(入口速度取10 m/s)，其余边界条件与剪切流中一致．3 结果与分析图3和图4分别是采用表1中第3种入口条件计算均匀流和剪切流时湍动能和耗散率的fluent数值模拟结果值在整个流场的分布云图，可以看出：在均匀流和剪切流中，湍流变量随着下游流动均有不同程度的衰减，且在剪切流中，湍流变量的衰减还受高度的影响，即速度梯度的影响，速度梯度越大，衰减越快．下面将着重分析入口湍流强度和湍流黏性比对下游湍动能衰减的影响．10.13245/j.hust.210915.F003图3均匀流中湍动能和耗散率的衰减云图10.13245/j.hust.210915.F004图4剪切流中湍动能和耗散率的衰减云图图5是采用标准k-ε湍流模型时得到的湍动能的fluent数值模拟解、Matlab理论数值解和均匀流数值模拟解的对比及其局部放大图，图中：x为下游距离；纵坐标为流场中当地湍动能(k)与入流湍动能(kin)的比值；图中从上到下的7个分图分别为采用表1中7种入口边界条件计算得到的结果．其中，剪切流的数值模拟解取自图1中参考线A和B位置，参考线A和B分别距流场下边界10和100 m．均匀流的数值模拟解取自图2中参考线A位置，参考线A位于上、下边界的中间．从图5可以看出：对于表1中的后3种入口条件，fluent数值模拟解和Matlab求得的理论数值解有很好的一致性；但对于前4种入口条件，fluent数值模拟解和Matlab求得的数值解偏差较大，这是因为在Matlab的求解过程中只考虑了生成项和耗散项，但实际上在流向上由于湍动能的衰减存在湍流变量的梯度，即在fluent数值求解过程中考虑了扩散项的影响，通过fluent数值解和Matlab理论数值解的对比可以分析湍流变量输运方程扩散项的影响，可见当湍流黏性比较小时，扩散项对湍动能的衰减影响较大．从局部放大图可以看出：当采用第一种入口条件时，剪切流中湍动能的衰减整体比均匀流中小一些，而且随着高度的增大，衰减情况接近于均匀流．尤其在10 m高度处，湍动能反而增加．这种增加趋势和湍动能的生成项有关，湍动能的生成主要取决于涡黏性和速度梯度，在同一种入口条件下，涡黏性不变，越靠近地面，速度梯度越大，导致湍动能生成项比较大，而在100 m高度处，速度梯度很小，所以生成项对湍流衰减几乎没有影响．当湍动能呈现衰减趋势时，速度梯度较大时的衰减速度要比速度梯度小时稍微快一些．当入口湍流黏性比增大时，不同速度梯度下湍流衰减变化不大，而当湍流强度增大时，不同速度梯度下湍流衰减变化较为明显．即在剪切流动中，入口湍流强度对湍流变量输运方程生成项的影响较大，湍流黏性比对湍流变量输运方程生成项的影响较小．10.13245/j.hust.210915.F005图5采用不同入口条件时湍动能数值模拟解和理论数值解的对比同时可以看出：当入口给定均匀的湍流参数时，不同高度的湍动能衰减程度不一样，导致在下游湍动能在各个高度不一样，均匀性不能保持，即单纯的风速剪切和均匀的湍动能是无法共同存在的．当入口湍流强度一样时，入口湍流黏性比越大，速度梯度对湍动能衰减的影响越小；当入口湍流黏性比相同时，入口湍流强度越大，速度梯度对湍动能衰减的影响越大．从整体看，当入口湍流黏性比相同时，入口湍流强度越大，湍动能衰减的越快；当入口湍流强度一样时，入口湍流黏性比越大，湍动能衰减的越慢，这与文献[7]中分析均匀流时的结论一致．对比分析了均匀流和剪切流中湍流动能的衰减规律可得到以下结论．a.在剪切流中，存在湍流动能衰减现象，且当湍流黏性比较小时(如0.01)，扩散项对湍动能衰减的影响较大．b.入口湍流强度对湍流变量输运方程生成项的影响较大，湍流黏性比对湍流变量输运方程生成项的影响较小．当入口湍流强度一样时，入口湍流黏性比越大，速度梯度对湍动能衰减的影响越小；当入口湍流黏性比相同时，入口湍流强度越大，速度梯度对湍动能衰减的影响越大．当湍流强度较小，速度梯度较大时，湍动能随着下游距离反而增加．c.在剪切流中，当入口湍流黏性比相同时，入口湍流强度越大湍动能衰减越快；当入口湍流强度相同时，入口湍流黏性比越大湍动能衰减越慢．