强吸气旋转圆筒壁面湍流边界层建模及计算1)

Ievgen Mochalin 林静雯 蔡建程,2) Volodymyr Brazhenko 鄂世举

*(浙江师范大学工学院,浙江金华 321004)

†(浙江省城市轨道交通智能运维技术与装备重点实验室,浙江金华 321005)

引言

旋转可渗透圆柱外的流动是在Taylor-Couette 流动的基础上叠加径向穿透流动[1-2],它具有很多应用,旋转过滤是一个典型例子[3-5].含杂质溶液通过旋转圆筒膜可以大为提高过滤效率,并且旋转膜表面的固相沉积率低不容易阻塞,因此有学者专门分析旋转过滤圆柱表面的流体运动[6-10].

旋转柱体的冷却是Taylor-Couette 流动的另一个研究方向,F´enot 等[11]综述了转子-静子系统中的同心圆筒间隙内的不可压缩流动与换热问题研究.Mochalin 等[12]研究了外筒静止、内筒旋转工况下的间隙流动,发现内筒做成完全多孔或是条状多孔可以使换热效率提高1.5 倍.因此,研究旋转多孔圆筒表面附近的流动很有意义.

旋转圆筒外的流动由于受到离心不稳定的作用,会导致二次涡流即Taylor 涡的产生.对于同心不可渗透圆筒间的流动,在较低的旋转速度时就能发生流动不稳定性,体现为轴对称Taylor 涡的产生,随着旋转速度的提高,流动模式会改变[13].对于经典Taylor-Couette 流动,相比于由离心不稳定弓起的Taylor 涡产生,两圆筒间的流动由层流到湍流的转捩发生较迟.通过多孔圆筒弓入径向穿透流,可以使间隙流动变得稳定,从而Taylor 涡的产生大为延后[14-15].Mochalin 和Khalatov[16]的研究表明强吸气穿透流可以防止Taylor 涡的产生直到非常大的旋转速度,并且发现旋转吸气圆筒壁面弓起的旋转流体运动集中在其壁面附近,而外面静止圆筒对里边旋转圆筒附近的流动影响很小,所以旋转流动可以视为旋转多孔圆筒表面的边界层流动.边界层厚度依赖于旋转速度和径向穿透流强度.

研究表明,多孔壁面径向流可以使边界层从层流转捩到湍流发生在离心不稳定(Taylor 涡) 之前,这就为边界层湍流强度控制提供了新思路.旋转壁面流动具有很多应用价值,比如相分离、表面冷却等[12,16].所以,有必要研究旋转多孔圆筒表面的湍流边界层流动.

Mochalin 和Khalatov[16]通过求解Reynolds 平均Navier-Stokes 方程和Reynolds 应力湍流模型,研究了同心圆筒间的流动.梁田等[17]则利用Spalart-Allmaras 湍流模型计算边界层抽气,进行叶栅角区分离流动控制的研究.在这些研究中,近壁面流动使用非常稠密的低Reynolds 数计算网格,计算量大.本文的研究目标是对于强吸气旋转壁面流动使用简单而可靠的预测方法,并提供计算步骤.其研究思路是使用恰当的代数湍流模型,基于详细数值模拟结果进行代数湍流模型参数的调整和校准,使其适用于一定工况范围.

1 强吸气壁面旋转流动速度分布的解析表达式

1.1 基本假设与无量纲参数化

本文所研究的同心圆筒间流动为半径R外的多孔圆柱壁外的稳定流动(在离心不稳定发生之前).旋转流体运动主要集中在厚度为δ 的边界层内.多孔壁面吸气,假设均匀吸气速度分布Vs.旋转流动认为在轴向不变化,即Vz≡0,∂/∂z≡0,且轴对称即∂/∂φ ≡0(φ 为方位角).

旋转流动可使用Reynolds 平均Navier-Stokes(NS) 方程联合连续方程进行描述.显然,使用柱坐标系(r,φ,z)更为方便,在柱坐标系下稳态不可压缩N-S方程

式中,Vr,Vφ分别为径向和周向的平均速度,ρ 为流休密度,νe为流体的有效运动黏度,它等于分子黏性ν和湍流黏性νt之和,νe=ν+νt.

连续方程为

积分上述方程,注意到壁面Vr(R)=Vs条件,得到

利用式(2),方程(1)可以转化成

在本文中,边界层的计算在径向划分为多个微元段Δr,在每一段中νe认为是常数.遵循这一假设,式(4)中最后一项在每个微元段可以忽略.进一步对式(4)中最后一项在整个边界层内都忽略,这是附加的假设.其有效性将从本文方法的结果与详细的数值仿真结果相比较看出.利用式(3) 以及上述假设,积分式(4),得到

上式中的积分常数C可由旋转壁面的边界条件

确定,其中Ω 为壁面旋转角速度.式(5)中的速度梯度可用壁面处的摩擦应力

来表示.

联合式(5)、式(6)和式(8),得到

为了得到无量纲表达,使用圆筒半径R、壁面旋转速度W和动压作为长度、速度和压力的参考量,通过下式定义无量纲量

其中,ξ 为无量纲径向坐标,ur,uφ为无量纲径向和周向速度,q为无量纲压力.从而得到无量纲化的式(9)为

式中,uw=Vw/W为无量纲摩擦速度,us=Vs/W为无量纲吸气速度.式(11)两边同除以RW,得到

为求解式(12),将旋转圆柱壁面法向的边界层分成足够小间隔Δξi,在每个间隔内湍流黏性看成为常数,式(12)简化为

积分后得到

积分常数C由下列边界条件得到

式中ξ0i,uφi分别为无量纲径向距离和每个径向间隔Δξi的周向速度;对于第一步,由式(6)、式(10)可知ξ01=1,uφ1=1.把式(14)代入式(13),得到常数C的定义

因此,在每个间隔内ξ ∈[ξ0i,ξ0i+Δξi],周向速度可由式(13)、式(15)得到

1.2 圆筒表面摩擦应力

壁面吸气使得两个圆柱面间流动的Taylor 涡产生大为滞后.这意味着径向Reynolds 数(Rer)足够大时,两圆柱面间流动为层流流动,其周向流动速度分布可表示为[16]

上式沿径向求导,得到

由该速度分布,进而可以得到摩擦系数

表1 给出了所考虑的旋转圆筒壁面摩擦系数.表中的数值仿真按文献[16] 中的方法得到,可以看出由式(18) 得到的结果与数值仿真结果非常接近,误差在1%之内.表1 中的第1、3 行的Reφ,Rer组合,其流态为层流,其余为湍流.基于此,式(18)可以应用于旋转同心圆筒间的湍流流动.所以,有如下关系式

利用式(19)可以化简式(16),得到

这里顺便指出,已有的对吹吸气槽道湍流的研究表明吹气可以使壁面阻力系数下降,而吸气会使壁面阻力系数上升[18-19],而高频周期吹气扰动在狭缝下游产生明显的减阻效果[20],另外通过在平板壁面施加不同频率振幅的周期性扰动,能进行湍流边界层的主动控制减阻[21].表1 的第2 行和第3 行表明,固定旋转速度(Reφ一定)的前提下,摩擦系数随着吸气强度Rer的增加而增加.而随旋转速度Reφ的提高,摩擦系数下降.所以当Reφ和Rer同时提高时,摩擦系数有可能下降.

表1 强吸气旋转圆筒表面摩擦应力Table 1 Friction coefficient values on the surface of a rotating cylinder with strong suction

2 代数湍流模型

2.1 湍流黏性的基本关系式

壁面吸气可以延迟边界层从层流到湍流的转捩[22-25].本文研究边界层除吸气外还涉及离心力和旋转对称等因素影响.文献[16] 中的数值仿真诚然可以考虑这些因素,但计算量过大.本文提供更为快捷的代数湍流型模型,以求解旋转边界层.思路如下:寻找一个合适的代数湍流模型,根据新增加的影响因素进行改进,并依据特定的情况用实验或计算数值结果进行参数校准.

本文依据Cebeci-Smith(CS)湍流模型[26],校准数据使用文献[16] 的结果.CS 模型能考虑压力梯度、壁面吸(吹)气、可压缩性和低Reynolds 数效应等因素的影响.对每个影响因素,分别弓入一个因子(经验参数或是含经验参数的表达式),然后耦合所有的因子.

平板流动的经典CS 两层湍流模型可以表述如下

式中νti,νto为湍流边界层内层与外层的湍流运动黏性系数;ym为νti=νto时的壁面法向距离;Vx,Vy为流向和横向的速度;lm为Prandtl 混合长度;y为离壁面的法向距离;Vm为边界层厚度δ 内的最大速度;δ*为边界层位移厚度.无量纲距离y+定义为

经验参数

式(22)可从下式

广义的CS 模型[26]可以考虑如下因素:趋于边界层外界线时湍流强度的衰减、外层低Reynolds 数对湍流涡的影响、壁面吸(吹)气、压力梯度、可压缩性.本文不考虑最后两个因素,仅考虑前3 者.

边界层外层界线附近湍流强度的衰减,由Klebanoff间歇系数来表示

用它乘式(24)右端.

低Reynolds 数的影响通过改变式(24) 右端的Clauser 系数k实现

式中尾迹参数Π0由Schlichting[24]建议为

式中

δ**为动量厚度

广义CS 模型中通过改变式(23)参数A*以考虑壁面吸气的影响

CS 模型中并没有考虑流线曲率的影响,流线曲率会产生离心力.本文利用Bradshaw[27]的思想,这一思想被许多学者所采纳[28-30],它把流线曲率的影响与浮力效应相类比,弓入Richardson 数Ri(表示由于离心力和速度梯度产生的湍流能量之比) 到湍流黏性系数中

式中,νt0为不考虑流线曲率的湍流黏性,β,m为经验常数.β,m的选择依赖于用于校准的剪切流,这将在下文中简述.根据Bradshaw[27]的建议,对于具有周向和轴向分量的旋涡流动

所以,在广义CS 模型中,考虑到趋于边界层外界线时湍流强度的衰减、外层内低Reynolds 数对湍流涡的影响、流线曲率、壁面吸(吹)气等方面后,式(22)和式(24)改写为

式(23)中的A*由式(34)给出,而式(38)中的k由式(30)～式(33)得到而不是式(27)中的常数0.016 8.式(38)中的Klebanoff函数Fkl由式(29)给出.

利用式(3)和式(10),上式可以写为

对于壁面法向距离y及无量纲壁面距离y+式(26)有

由式(23)、式(34)、式(39)、式(40),边界层内层湍流黏性式(37)可写为

经验参数β,m,C1的选择在2.3 节讨论.

对于式(36)的Richardson 数,考虑到ξ ≈1,有

对所研究的问题Vm=W,所以由式(38)可以得到边界层外层的黏性系数

由式(29)给出Klebanoff间歇系数Fkl.

由于研究对象为旋转流动,定义位移厚度应使用相对速度W-Vφ,式(25)可写为

动量厚度

2.2 外层湍流黏性公式的修正

本文研究的旋转圆管壁面边界层与平板边界层有区别,其轴对称性使计算需特殊处理.对于平板边界层,边界层从平板前缘点开始,向下游(纵向)发展.下游各横截面上边界层厚度、位移厚度和动量厚度由上一横截面上信息计算得到,这是因为平板边界层沿纵向是抛物型偏微分方程.对于本文中的旋转轴对称边界层,由于初值未知,因此需要使用迭代方法计算,即先设置初始猜想值,通过迭代收敛求解边界层.

下面进行迭代求解的说明.图1 示意了静止壁面和旋转壁面边界层速度分布.实线下的面积记为S1,S2,虚线下的面积记为如果黏性增加,则将导致层与层之间的动量交换.对于静止的平板边界层:.

图1 静止壁面(左)和旋转壁面(右)边界层:1–νt 小;2–νt 大Fig.1 Typical velocity profiles in the boundary layers on a stationary surface(Left)and a rotating cylinder surface(Right):1–less νt values;2–greater νt values

而对于旋转边界层

边界层位移厚度的计算依赖于式(12)所描述的uφ速度分布.由式(19)=us可知无量纲位移厚度主要依赖于us和Reφ.对R=0.11 m 以空气为介质和R=0.055 m 以水为介质的情况,在适用于低Reynolds数的计算网格上使用Reynolds 应力湍流模型进行数值仿真,得到不同us,Reφ组合的边界层厚度(见表2),基于此建立(us,Reφ)的关联公式.

对表2 中的数据分析拟合,得到无量纲位移厚度

表2 旋转吸气圆筒表面无量纲边界层厚度Table 2 Values of the dimensionless boundary layer thicknesses on the surface of a rotating suction cylinder

无量纲厚度

通过上两式得到相应边界层厚度的初始估计值.

式(49)～式(52)作为轴对称旋转边界层迭代求解的初值.

2.3 流线曲率和壁面吸气所对应的经验参数选取

由流线曲率以及相应的离心力对湍流影响由Bradshaw 修正式(35) 得到,结合到式(41) 和式(44)式中进行内层和外层湍流黏性修正.对于经验参数m,取m=1.β 依赖壁面曲率、吸气率和旋转速度.从黏性底层到湍流强烈的外层,β 的变化与形状因子有关,详细的流场数值模拟表明如下的β 取值较合理

吸气对边界层的影响通过式(41)中的经验参数A*修正.A*指数形式的式(42)中系数C1在平板边界层的CS 模型中取C1=11.8.对于本文研究的旋转边界层,C1与边界层中的湍流程度有关,而与湍流程度密切相关.对于层流边界层,=2.湍流黏性为=0,从式(41)和式(42)可知C1→∞;对于本文研究的已经发展的湍流边界层≈1.17,C1≈11.通过详细数值模拟的流场数据的校核标定,表明C1的较好近似

这与上述两个极限情况吻合良好.

注意到系数β 和C1校核与形状因子的确定是耦合的:一方面β 和C1依赖于形状因子的值,而反过来β 和C1值的选择又影响形状因子的值.基于已验证详细的数值模拟结果,计算形状因子值,建立标定所需数据库.对于不同Reφ和us的组合,使用所建立的近似方法进行不同β 和C1参数组合的试算,求得相应的值.然后根据梯度下降法改变β 和C1取值,使详细数值模拟得到的与基于代数湍流模型方法得到的之差在容许误差之内.有了β 和C1离散形式参数值集后,对其进行拟合,得到最小二乘拟合函数β()、C1(),其结果见式(53)和式(54).

尾迹参数Π0使用经典公式(31),它对壁面曲率和吸气的依赖通过合理选择参数β 和C1体现.

3 算法及验证

3.1 迭代计算过程

边界层流动方程(20)的近似解以及上述定义代数湍流模型的方程,表明可以用迭代法求解吸气旋转圆筒壁面的边界层流动.主要步骤如下.

(1)确定下列参数的初始值:

由式(20)旋转Reynolds 数Reφ;无量纲吸气速度us;径向空间步长kr;径向最大的无量纲尺寸=ξmax;无量纲位移厚度计算迭代容差ε*.

(2)第一步计算:

(a)从式(19)得到无量纲摩擦速度uw;

(c) 由式(51) 和式(52) 预测无量纲边界层厚度和动量厚度和;

(8)下面几步为边界层外层计算循环:

(a)当前步计算初始湍流黏性由式(30)、式(31)、式(43)、式(45)、式(46)、式(49)得到;

(c)当前步结束时计算导数∂uφ/∂ξ:

下一步长

(d) 当达到区域边界时,计算结束,否则回到(a)继续计算.

3.2 迭代计算方法的验证

为了验证上述方法的正确性,对比本文的代数湍流模型和详细流场数值模拟得到的边界层周向速度,其结果如图2 所示.图中也包括由式(17)描述的层流边界层.Reφ和Rer的组合覆盖层流边界层和离心不稳定发生前的湍流边界层.图2 中的Reφ和Rer组合并不是用于校准湍流模型参数的Reφ和Rer组合.

从图2 可以看出,本文的代数湍流模型所预测的旋转圆筒周向速度与Reynolds 应力微分形式湍流模型的详细数值模拟得到的结果几乎一致.并且它也能预测强吸气条件下的层流边界层速度分布,见图2(f).

图3 显示了边界层内的湍流黏性分布,可以看出代数湍流模型的预测结果与Reynolds 应力模型的数值计算结果比较接近.其中图3(a)为层流边界层,理论上湍流黏性为0,本文的代数湍流模型的湍流黏性预测值也非常小.

以上流场结果也证实,对于强吸气旋转圆筒Couette 流动,其主要流动为旋转圆筒壁附近的边界层流动.随着旋转Reynolds 数Reφ的增大,湍流及离心不稳定性增强.吸气量减弱(即Rer减小),边界层增厚,边界层内湍流强度增大.在大Rer情况下,离心不稳定发生大为延后.因此,可以通过合适的Reφ和Rer组合来控制边界层的厚度和湍流强度,具有很好的创新实用价值.

图2 旋转圆筒壁面边界层的平均速度分布:1–Reynolds 应力湍流模型的数值仿真结果;2–本文代数湍流模型结果;3–层流边界层,即式(17)Fig.2 Mean velocity profiles in the boundary layer on the surface of a rotating suction cylinder:1–detailed numerical simulation,RSM turbulence model;2–application of the algebraic turbulence model;3–laminar profile,equation(17)

图3 湍流黏性:1–Reynolds 应力湍流模型的数值仿真结果;2–本文代数湍流模型结果Fig.3 Turbulent viscosity ratio profiles:1–detailed numerical simulation,RSM turbulence model;2–application of the algebraic turbulence model

图3 湍流黏性:1–Reynolds 应力湍流模型的数值仿真结果;2–本文代数湍流模型结果(续)Fig.3 Turbulent viscosity ratio profiles:1–detailed numerical simulation,RSM turbulence model;2–application of the algebraic turbulence model(continued)

4 结论

本文针对吸气旋转壁面的湍流边界层流动建立了代数湍流模型.基于Cebeci-Smith 两层代数湍流模型,通过推导分析,对其进行修正和经验参数的调整,用以考虑离心力场(流线曲率) 和壁面吸气等因素的影响.利用Reynolds 应力湍流模型的详细数值模拟得到的流场,校准Cebeci-Smith 代数模型的经验参数.

建立了壁面吸气条件下轴对称旋转壁面边界层的迭代算法,它容易用计算机程序实现.该迭代算法能用于离心失稳(Taylor 涡)发生前的层流及湍流边界层计算.在旋转多孔圆筒强吸气条件下,离心不稳定在较大的旋转速度(Reφ～105) 条件下发生,通过调整吸气强度可以控制边界层厚度和湍流强度.这为旋转壁面的剪切力、动量传递、传热和传质的控制提供了新思路.本研究可应用于动态旋转过滤、转子轴的热保护、新型旋转热交换器等领域.