基于大涡模拟研究分形栅格激发湍流的传热特性

2022-10-18段成栋魏超中国矿业大学徐海学院江苏徐州221008

化工管理 2022年28期

段成栋，魏超(中国矿业大学徐海学院，江苏徐州 221008)

0 引言

不同湍流状态下的流动对换热效果影响显著，例如：湍流强度、压降、均匀各向同性等湍流参数对于换热效果具有重要作用。湍流对圆柱体传热的影响在很多工程应用中有重要意义[1]。在实际工程中，通常使用圆柱形几何体在流体和壁之间进行热交换，例如管壳式热交换器、压水反应堆、水-空气散热器等[2]。

1935年，G.I Taylor在风洞中放置了多排栅格，当均匀气流垂直流过栅格时，会发生不规则扰动。这种不规则扰动向下游移动过程中，由于没有外部干扰，逐渐演化为均匀各向同性湍流。相似的均匀各向同性湍流已经在风洞或水洞中得到了广泛的研究。2007年，Hurst等[3]使用各种多尺度(分形)栅格在风洞中产生湍流，发现复杂的多尺度栅格结构会极大地影响湍流行为，并通过实验方法探索了分形维数Df≤2、栅格有效尺度Meff和网格厚度比tr=tmax/tmin。这三个因素对分形栅格激发湍流的流动机制、各向同性、均匀性和衰减特性有很大影响。

2016年，G.Melina等[4]研究了不同几何形状栅格激发的湍流场中心线上的流动：分形栅格(FSG17)、单一方形栅格(SSG)、规则栅格(RG60)。结果表明，在湍流产生区，单一栅格结构产生的湍流强度大于其他两种情况，这是由于较强的涡旋脱落效应，而分形栅格抑制了涡旋脱落。2015年Torrano等[5]将实验测量数据与基于RANS的栅格生成湍流衰减数值研究进行了比较。结果表明，基于RANS的数值模拟捕捉湍流流场大尺度特征的能力在可接受的精度范围内。R.J.Hearst等[6]设计了一种特殊的分型栅格结构，还增加了分形栅格实验的下游测量范围。结果表明，特殊栅格产生的湍流行为与之前分形栅格附近处观察到的类似，远下游湍流性能则与规则栅格产生的类似。

上述结论为研究湍流参数对传热效果的影响提供了理论依据和实验依据。Smith等[7]提出了一个理论模型，假设在前驻点附近，雷诺应力与自由流中的湍流强度Tu和与壁面的距离成正比。他们的模型得到了实验结果的支持，实验结论表明前驻点处的弗罗斯林数与湍流参数TuRe0.5成正比，其中Re是基于圆柱体直径D和平均流向速度的雷诺数。Kestin等[8]以及Lowery等[9]分别测量了风洞中栅格产生的湍流场中圆柱体的传质和传热。这两项研究都将NuFsp/Re0.5的值与TuRe0.5的二次多项式函数相关联。此外，几项实验表明，流动中的积分长度标度Lu也非常重要[10]。Van Der Hegge-Zijnen[11]报道，在Re和Nu相同条件下，当0＜LU/D＜1.6时，Nu随LU/D增加；当LU/D＞1.6时，Nu随LU/D减小，其中“L”定义为湍流的长度尺度。最近，为了强化传热，人们广泛研究了流场中的涡旋结构，并通过在流场中添加涡流发生器实现了良好的传热效果。

随着计算机技术的进步，数值模拟已成为研究湍流的重要手段。常用的数值模拟方法有三种：直接数值模拟(DNS)、雷诺平均Navier-Stokes方程模拟(RANS)、大涡模拟(LES)。DNS可以获得湍流统计的深度数据，包括壁面附近的高阶量，而它只能用于求解非常有限类型的湍流。RANS能以最经济的方法正确预测高雷诺数的复杂流动，但结果是时间平均的，不能反映流场的瞬时信息。这类模型的典型例子是k-ε或k-ω模型。最近，许多雷诺应力模型都是通过修正湍流的压力-应变关系发展起来的[12-15]。LES直接预测大尺度涡旋运动，对较小的尺度(亚格子尺度)和更接近各向同性的涡旋[16]，采用亚格子尺度(SGS)模型。不仅可以获得详细的流场动态信息，而且可以节省计算成本。因此，大涡模拟得到了越来越广泛的应用。

鉴于上述分析，文章采用大涡模拟的方法研究了三种不同截面栅格产生的湍流以及恒温圆柱在流场中的传热。三种湍流发生装置的横截面形状不同：三角形、倒三角形和正方形。通过分析流场中速度、湍流强度、压降和努塞尔数等瞬时和统计量，揭示了不同扰动下流体运动和传热机理，以及湍流特性与传热之间的关系。

2 数值模拟方法

2.1 计算域与网格划分

基于LES研究流体通过分形栅格产生的不同流场对圆柱体产生不同的热交换效应。如图1所示，具体的模拟方案如下：在长方体风洞中距离入口0.1 m处放置三种不同形状的分形栅格，第四种不放置任何扰动作为参照。在所有工况中的同一位置各放置一个直径为D=0.05 m、长度为0.2 m、温度为323 K的恒温圆柱体。为了消除尺寸因子的不可比性，计算域基于换热圆柱直径进行无量纲化。无量纲后的计算域为Lx×Ly×Lz=24D×4D×4D。

图1 计算域

分形栅格的特征由以下参数定义：

—分形结构的维数N，此研究中N=3；

—某一维度分形图案的数量4j；

—栅格长度Lj=0.5jL0和横向厚度tj=0.5jt0，L0=0.5Ly，t0=tmax=0.01 m；

—厚度比tr/tmax/tmin=4，即形成最大正方形的横向厚度与最小正方形的横向厚度之间的比率；

—三棱柱栅格的角度α=53°。

文章重点研究了三棱柱分形网格中湍流的演化规律。

考虑到三角柱相交处的结构网格质量难以保证，在工况1和工况2的整个计算域中使用非结构网格。为了保证更高的计算精度，设置y+=Δy×u∞/ν=1，因此，第一层网格厚度只有0.005 6 mm，y+是从第一层网格的质心到壁面的无量纲距离。在工况1和工况2中，使用MESH生成四面体非结构化网格，选择曲率尺寸函数来控制曲率较大的区域和网格生长，以及表面相邻区域的分布。将最小尺寸、最大面尺寸、最大尺寸分别设置为2.0×10-5m、0.002 m和0.005 m。总节点数为1 629 407，总网格数为8 917 056。

由于工况3和工况4中的几何模型较为简单，因此工况3和工况4中使用了在流向、法向方向和展向方向上网格数为624×200×160的结构化网格。

2.2 控制方程

LES的基本思想是通过过滤将大尺度涡旋与小尺度涡旋分离。为了进行湍流计算，连续性方程和Navier-Stokes方程过滤如下：

式(1)～(3)中：ρ为流体密度；ν为流体的动力黏度；为速度的过滤值；为压力的过滤值；τij为子网格应力项(SGS)，反映了小尺度涡旋对大尺度涡旋的影响。

为了使过滤后的方程封闭，引入以下方程：

式(4)～(7)中：τkk为亚格子应力的各向同性部分，在不可压缩流中通常可以忽略，根据Kolmogorov’s-5/3定律，通常Cs取0.10～0.30之间；Δ≈0.86 mm为滤波宽度，用来区分大尺度和小尺度涡旋的特征长度尺度；μt为SGS黏性项；为平均应变率张量。

2.3 边界条件

LES是一种瞬态湍流模型，将时间项中引入其边界条件。

(1)初始表压值为0 Pa。

(2)所有情况下的入口速度都设置为u∞=5 m/s，雷诺数为Re=u∞l/ν=13 611，其中l是长方体风洞的水力直径。工况1～3中，在速度入口处不增加扰动；工况4中，在入口处设置初始扰动。

(3)在风洞出口边界处使用出流边界条件。

(4)风洞壁面和格栅采用无滑移条件和恒温壁面T=283 K。圆柱也被设置为无滑移壁面，温度为323 K。

在本研究中，LES模拟使用有限体积求解器进行，该求解器采用有限体积法(FVM)在瞬态假设下求解不可压缩Navier-Stokes方程。压力项采用二阶离散格式，动量方程采用有界中心差分离散格式。数值计算方法采用PISO算法，亚格子模型采用Smagorinsky-Lilly模型，且Cs=0.10[17]。

3 结果与讨论

3.1 仿真验证

为了验证本文模拟结果的正确性，用相同的模拟方案对参考文献[4]中研究的RG60流场进行了模拟。在图2中，比较了数值模拟和风洞试验中沿风洞中心线的流向平均速度和湍流强度，结果吻合良好，计算误差小于25%。

图2 沿风洞中心线的平均速度和湍流强度

3.2 流动换热特性

当流体绕圆柱流动时，由于流体与壁面之间的摩擦阻力，圆柱周围形成边界层，边界层区域内的流动为层流。各层之间没有质点掺混，因此主要通过热传导进行传热。边界层分离后，由于宏观上的相互混合，热量被涡旋流动带走，其中的热交换主要是对流换热。

3.2.1 速度分布分析

文章通过对瞬时流向速度云图和数据的分析，得到了湍流场的特征参数和传热效果。由于计算域中的流动在流动开始时尚未完全展开，仅选择0.35～0.55 s时间段的流场数据进行分析。所有工况都使用相同的笛卡尔坐标系，且原点相同，x方向为流向方向，y方向为展向方向，z方向为垂向方向。

定性分析方面，提供了x-y平面(y=0)上的等值面云图以及y-z平面在x/D=2，x/D=10，x/D=12.5,的等值面云图。而在定量分析上，展示了沿着风洞中心线上和圆柱体周围的各变量变化趋势。

图3是t=0.35 s时的速度云图，U是瞬时流向速度。图中显示，当流体通过格栅时，格栅后部会产生一个局部低速区，涡旋逐渐在低速区产生，并沿流向方向交替脱落。对比3种工况，工况2产生的局部低速区域明显大于其他两种情况。从x/D=2平面上看，在栅格下游附近流动速度云图具有非常明显的几何印记，其中工况3最为明显。在工况3和工况4中，流体绕流圆柱时，圆柱后部会出现较为明显的卡门涡街。原因是：在这两种情况下，圆柱体周围环境的流动湍流强度较小，一旦受到扰动就会发生强烈地速度脉动。

图3 在0.5s时，x-z平面(y=0)和不同横流平面(x/D=2,10,12.5)速度云图

图4 显示了沿流场中心线的流向平均速度分布。对工况1～3，由于栅格阻塞的原因，0≤x/D≤2范围内流向平均速度明显增加。对工况2，平均速度最高可达8.3 m/s，并且所有工况中速度衰减最小。此外，对工况3和工况4，在圆柱后部产生涡旋，并沿流向交替脱落。圆柱后部附近的流向平均速度急剧增加，并伴随着震荡，直到速度恢复到与进口流速相当。

图4 沿中心线的平均速度剖面

3.2.2 湍流强度

图5显示了四种情况下沿中心线的湍流强度(tu)分布。湍流强度由公式(8)计算：

式中：u′为瞬时速度的均方根为时间平均流向速度。

从公式(8)中可以看出，湍流强度与随机湍流波动直接相关。较高的湍流强度意味着更大的随机湍流波动。

如图5所示，工况1～3中圆柱前面区域湍流强度均有不同程度的递增趋势。在第三种工况下，流场产生的最高的湍流强度可达20%。但是，在没有栅格的情况下湍流强度是恒定的。在圆柱后部区域，由于流场受到圆柱的扰动，湍流效果会进一步加强。其中，工况4产生的湍流强度最大约为56%，在工况1中约为48%。而且，在工况3和工况4的圆柱后部附近，湍流强度呈波浪式减小趋势。

图5 沿中心线的湍流强度

3.2.3 压降分析

为了表示流场中圆柱周围的压力分布，文章引入了无量纲压降公式，定义如下：

式中：P为圆柱周围的总压力；Pinlet为入口的总压力；ρ、u分别为入口处的平均速度和密度。

图6显示了气流在圆柱体周围流动时的边界层分离现象。图7显示了圆柱体周围的速度分布。在所有工况下，平均速度沿抛物线分布，先增大后减小，最后趋于稳定。当内压梯度小于零时，边界层内流动加速，流动动能增大。当内压梯度大于零时，边界层中的流动减速，动能迅速耗尽。当边界层中的剪切应力大于上游流体中的剪切应力，就会出现与初始流动方向相反的回流，从而形成了边界层分离。对工况1和工况2，边界层分离发生在φ≈100°；对工况3和工况4，边界层分离发生在φ≈160°。

图6 边界层形成示意图

图7 圆柱周围的速度

在图 8中[18]，展示了圆柱周围 0°≤φ≤180°的压降系数分布。其中工况1和工况2在圆柱周围呈线性下降趋势，对于情况1，压降系数的变化非常小。对情况3和情况4，压降系数的趋势几乎相同。压降系数在φ≈0°时最大，在 0°≤φ≤180°时呈指数下降。φ≈80°压降系数降至最低点，并开始出现回升。在80°≤φ≤100°阶段压降系数略微上升，然后变得稳定。同时，边界层动能耗尽开始出现分离。最终，工况3的最小压降系数维持在-4.5左右，低于工况4的-3.5。

图8 圆柱周围的压降系数分布

3.2.4 涡结构

图9显示了四种情况下的涡结构。在工况1和工况2中，流场中栅格附近和圆柱后部产生了少许的涡结构，而工况3无论在栅格附近还是在圆柱周围都产生了较多的涡旋结构，这是因为工况3的栅格结构前后压强差较大，更容易在流场中诱导涡旋产生。

图9 介于-51 000和-50 000之间涡结构

3.3 换热特性

图10 显示了加热圆柱表面的Nu分布的俯视图、前视图(面向入流方向)、侧视图和后视图。图中可以看出前三种工况在边界层处均有对流换热和热传导现象，而且边界层中的换热效果要强于尾流区域中涡流脱落产生的热对流。这是由于栅格结构在上游流场产生的一定数量的涡结构导致。θ是边界层分离角，即x方向与边界层分离点和圆心之间的线的角度。直观上看，工况1和工况2分离点均发生在，远大于工况2和工况3的θ≈85°。此外，工况1～3中Nu的分布在边界层区域是非常不均匀的，特别是在案例3中，存在相当一部分Nu≥6 000的区域，这些良好的换热效果是由于流场中大量的涡结构存在。

图10 圆柱体的表面努塞尔数

图11 是平面的温度云图。直观上看，图11中得出的结论与图8、图9和图10有很好的匹配。图12显示了圆柱表面0°≤φ≤180°努塞尔数的分布情况。在 0°≤φ≤78°，工况 3的最大Nu约为 6 200，比工况4高约30%，比工况1和工况2高约63%。然而，它在105°≤φ≤180°之间急剧下降，达到最低值1 000后逐渐恢复，最后达到约2 600。这表明，工况3的热交换效果是不均匀的。从面平均努塞尔数上来看，工况3的换热效果是最好的，平均能达到3 870。