摩托发动机冷却风扇流场的数值模拟

2019-10-18谭礼斌袁越锦冷小丽余千英

陕西科技大学学报 2019年5期

谭礼斌,袁越锦,王萍,冷小丽,唐琳,黄灿,余千英

(1.陕西科技大学机电工程学院,陕西西安 710021；2.重庆隆鑫通用动力股份有限公司技术中心基础研究所CFD研究室,重庆 400039 )

0 引言

工程机械如输送泵、挖掘机、压路机等、轻工机械如造纸机、多烘缸纸张干燥机、通用发电机组等设备使用工况都比较恶劣，且运行时间特别长.长时间的运行工作容易造成机械产品灰尘累积、温度高等问题，从而使产品寿命降低.因此，对工程机械、轻工机械使用的发动机的冷却系统的研究越来越受到人们的重视[1].其中，冷却风扇是发动机冷却系统中最核心部件之一，其工作性能的好坏直接影响发动机的工作性能及使用寿命[2].冷却风扇的主要功能是以风扇叶片的旋转控制冷却风流量，合理地实现发动机周围各高温部件的冷却，从而使发动机能够在最佳状态下运行，功率最大化[3].冷却风扇在发动机不同转速下的流量是评估风扇性能的一个重要指标，合理的风扇风量值才能更好地保证发动机设备总体的工作性能及其运行可靠性.

采用试验测试的方法可以获取冷却风扇性能，为风扇设计与优化提供参考，但由于试验周期较长、成本高，在产品前期开发及性能评估应用越来越少.随着计算流体动力学仿真(CFD:Computational Fluid Dynamics)技术和高性能计算机的不断发展与运用，基于CFD数值仿真方法来预测冷却风扇性能已成为旋转机械产品行业的趋势.国内外学者对冷却风扇内部流场及其性能分析进行了大量的仿真模拟及实验研究[4-18]，如陈慧敏等[19]利用ANSYS Workbench软件对风扇内部流场进行了数值模拟，分析了静压、功率、效率与风扇风量间的关系；张代胜等[20]利用流体分析软件Fluent对冷却风扇进行了流场和气动噪声数值模拟分析，得出在降低噪声的同时必须保证冷却风扇的性能;Pogorelov等[21]利用商业流体分析软件CFX研究了叶片叶顶间隙尺寸对风扇流场的影响；Gallonia等[22]利用ANSYS Fluent研究了叶片形状和数量对某电动摩托冷却用径向风扇的工作性能的影响，为叶片优化提供了初步建议.由此可见，依据数值仿真分析的仿真预测结果来进行产品性能的评估和优化方案的寻找已经在旋转机械产品领域得到广泛的运用.

目前国产摩托车采用的强制风冷发动机在发动机工作性能检测工作中都普遍存在热负荷过高的问题，因此，摩托发动机冷却风扇性能的评估对发动机散热性能的评测极其重要[23].摩托车发动机冷却风扇结构简单，体积较小，工作性能的好坏对发动机的冷却效果具有极大的影响.目前对摩托车发动机冷却风扇性能测试实验平台搭建及仿真分析对比的相关研究较少.采用实验分析研究冷却风扇的性能可以评估冷却风扇对发动机的冷却效果，同时利用实测的参数来标定仿真计算模型，可以验证仿真计算模型的准确性，对冷却风扇结构的改进和优化具有非常重要的指导意义.因此，本文以某摩托车发动机冷却风扇为研究对象，搭建摩托车发动机冷却风扇的测试平台；同时，基于计算流体动力学方法(CFD)创建与测试平台等比例的流体计算域物理模型，利用商业流体分析软件STAR-CCM+11.06对冷却风扇的内部流场特性进行数值模拟分析，预测其在运行工况下的风扇性能，并与实验测试的风扇风量值进行对比验证分析，验证CFD数值仿真的可靠性，同时评估实验测试数据的可靠性.本文研究结果可为后续冷却风扇的结构改进及性能提升提供相应的优化思路和理论参考.

1 冷却风扇实验测试

某型发动机冷却风扇实验测试在某司工程研究中心完成，测试方法和测试步骤完全按照GB/T 1236-2000《工业通风机用标准化风道进行性能试验》执行[24].冷却风扇实验测试设备图如图1所示.

图1 冷却风扇实验测试装置图

将用于测试的冷却风扇样件安装在图1所示的测试工装上(按照冷却风扇大小制定等比例的安装工装)，并连接到风扇测试风洞装置上进行测试.测试过程中给定风扇输入电源为12 V，冷却风扇的转速通过红外转速仪测定.测试不同风压时风扇的流量、转速、功率，通过数据采集系统与计算机连接，实时采集测试结果.风压通过测定喷嘴前后的静压值来设定，静压值通过内外压差变送器(日本横河公司生产)来控制，压差变送器测试精度为±0.175%FS,控制精度为±2 Pa.通过多喷嘴流量测试系统可以获得冷却风扇不同风压下的流量大小.其中，单个喷嘴的流体流量(风量)与喷嘴流量系数、喷嘴喉部面积、喷嘴前后静压差、流体密度有关系，其中喷嘴流量系数与流体物理属性有关[24].通过单个喷嘴的流体流量按下式计算：

(1)

式(1)中：Q为通过喷嘴的流体流量，m3/h；C为喷嘴流量系数，根据雷诺数可以获得流量系数[24]；A为喷嘴喉部面积m2；Δp为喷嘴前后的静压差，Pa；ρ为喷嘴喉部的流体密度，kg/m3.

冷却风扇测试风洞原理简图如图2所示.冷却风扇测试过程及评估标准依据中国国家标准测试流程一一执行，获取的冷却风扇风量实测结果具有一定可靠性，用于与仿真计算结果对比验证分析，验证基于CFD数值仿真方法预测冷却风扇工作性能的可靠性.表1为测试得到的冷却风扇在不同风压下对应的转速及冷却风扇流量值.

图2 风扇测试风洞原理简图

表1 冷却风扇风洞实验测试结果

2 物理模型

某发动机冷却风扇实物图及等比例绘制的3D模型图如图3所示.将等比例绘制的3D模型在ANSYS SPACECLAIM 2016软件中进行前处理及流体计算域物理模型的创建，完成后将计算域物理模型导入通用流体分析软件STAR-CCM+11.06，采用多面体网格和边界层网格类型对风扇计算域模型进行初步网格划分，获得相应的计算域网格模型图(如图4所示).该计算域模型主要分为进口半球面区域、旋转流体域区域、出口区域(非旋转域)三个部分.考虑到发动机冷却风扇计算域模型中不同区域的流场变化不同，旋转区域的流场变化较剧烈，因此采用分区域划分网格的方法对冷却风扇计算域模型进行网格划分，旋转流体域的网格尺寸较小，进口半球面区域的网格尺寸较大，非旋转区域的网格尺寸最大.依据网格参数控制策略经验，设置初始网格尺寸并执行网格划分，网格划分完成后整个计算域模型的初始网格数量约为590万.流体计算域物理模型按照实验测试装置进行等比例绘制搭建，其中非旋转域半球形面的直径约为风扇叶片外径的4倍.

(a)实物图 (b)3D模型图3 冷却风扇模型示意图

图4 冷却风扇流体计算域物理模型示意图

3 数学模型

计算流体动力学基本控制方程包括连续性方程(质量守恒方程)、Navier-Stokes方程(动量守恒方程)、能量守恒方程[25,26].本文假设冷却风扇的旋转气流流动状态为稳态流动，流体介质为不可压缩流体，整个流动过程模拟中不考虑热量的交换，因此本文冷却风扇内部流场的数值模拟只需求解流体连续性方程和动量方程.在发动机冷却风扇数值模拟研究中，k-ε模型具有较高求解精度和收敛速度的特点而被广泛应用[27].因此，本文选定STAR-CCM+中Realizablek-ε湍流模型进行摩托发动机冷却风扇内部流场特性的数值模拟计算.依据上述计算模型的简化，冷却风扇数值模拟过程中涉及的相关数学方程表达如下：

(1)连续方程(质量守恒定律)

(2)

(2)Navier-Stokes方程(动量守恒定律)

(3)

式(3)中：ui、uj是平均速度分量，m/s，xi、xj为坐标分量，m；p是流体微元体上的压力，Pa；μeff是湍流有效黏性系数，Pa·s.

(3)k-ε两方程湍流模型

(4)

式(4)中：Gk为速度梯度产生的湍动能项，Gb为浮力产生的湍动能项，YM表为脉动扩张项，C1ε,C2ε,C3ε为经验常数，σk,σε分别为与湍动能k和耗散率相对应的Prandtl数，Sk和Sε为用户自定义的源项.

4 边界条件

冷却风扇流体计算域物理模型中主要分为风扇旋转域和非旋转域两个计算域，计算前需要将边界类型及边界数值赋予指定的边界，从而可以真实还原冷却风扇的测试工况，进行冷却风扇实验测试结果与仿真计算结果的对比验证分析.本文的边界条件设置如下：

(1)风扇旋转域：针对风扇的数值模拟，STAR-CCM+中通常有滑移网格(Sliding Mesh)、风扇动量源(Fan Momentum Source)和运动参考系法(Moving Reference Frame，简称MRF)，其中滑移网格用于瞬态模拟，动量源法用于无风扇实体的模拟[11,28,29].本文有风扇实体模型，且为稳态数值模拟，因此本文采用运动参考系法对某发动机冷却风扇性能进行数值模拟分析.通过MRF方法实现动-静面间的数据传递，给旋转叶片区域(旋转域)设定旋转速度和旋转方向，即可实现风扇的运转.冷却风扇的旋转坐标及旋转方向如图5所示，风扇的旋转速度为实验测试中风扇实际工况下的运行转速.

图5 冷却风扇旋转坐标及旋转方向示意图

(2)非旋转域入口边界：设置为滞止入口(stagnation inlet)，相对大气压力的总压为0 Pa.

(3)非旋转域出口边界：设置为压力出口(pressure outlet)，相对大气压力的静压值按照实验测试环境下各工况的测试风压值进行设定.

(4)固壁边界：发动机冷却风扇的叶片、轮毂面和非旋转区域(除去入口、出口)与流体域所有的接触面都设置为壁面(WALL)边界，采用绝热无滑移壁面条件，计算域壁面参数采用 STAR-CCM+推荐的Two-Layer All Y+Wall Treatment ，对近壁区进行处理.

本文采用迭代仿真的方法进行冷却风扇内部稳态流动特性的数值模拟分析，即采用场函数(Field Function)进行变工况边界条件的设置，计算冷却风扇在不同风压不同转速下对应的风量值.具体设置如下：

(1)风扇旋转速度边界场函数设置如下：${Iteration}<=5000?4213:(${Iteration}<=10000?4167:(${Iteration}<=15000?4172:(${Iteration}<=20000?4196:(${Iteration}<=25000?4194:(${Iteration}<=30000?4026:4000)))))

(2)压力出口边界条件场函数设置如下：${Iteration}<=5000?0:(${Iteration}<=10000?9.9:(${Iteration}<=15000?19.9:(${Iteration}<=20000?39.8:(${Iteration}<=25000?50.9:(${Iteration}<=30000?60:80)))))

整个冷却风扇仿真计算中，共计算六个工况点.每个工况点迭代求解步为5 000步，六个工况点的计算总步数为30 000步，即完成冷却风扇在0 Pa、9.9 Pa、19.9 Pa、39.8 Pa、50.9 Pa、60 Pa六个不同风压下进出口流量的计算.计算中采用STAR-CCM+的Auto-Save功能或Solution history对每个计算工况点的计算结果进行存储，存储设置为每5 000计算步存储一次Auto-Save File或Solution History.

5 模拟结果及分析

5.1 网格无关性研究

冷却风扇计算域网格的划分方法和网格参数控制策略对计算域模型的网格数量和网格质量有着直接的联系，而网格质量地好坏对数值仿真预测的结果具有一定的影响[30].采取分区域网格的划分方法和合适的网格参数对主要流动区域进行网格加密细化处理的策略相结合可提高网格质量，减小因计算域网格质量而引起的仿真计算误差.因此，本文对19.9 Pa、39.8 Pa两个风压下的冷却风扇计算模型进行网格无关性验证研究，以确定合适的网格划分策略.本文相关体加密示意图如图6所示.非旋转域的网格对计算结果的影响不大，对数值计算结果会产生相应影响的主要是旋转主流区域和风扇叶片网格.因此，本文非旋转域体网格加密尺寸设置为16 mm，过渡区域为保证较好的网格过渡，block加密的基本尺寸设置为8 mm.

图6 不同区域网格加密示意图

STAR-CCM+11.06中自带的网格生成工具可生成多面体网格、四面体网格、切割体网格三种非结构化网格，也可生成薄板网格、拉伸网格、周期网格等2.5D结构化网格.本文对比分析多面体网格、四面体网格、切割体网格三种非结构化网格划分方法对冷却风扇风量计算结果的影响(风扇叶片基本尺寸设置为1 mm).图7表示不同网格类型对冷却风扇风量计算结果的影响.从图7可以看出，不同网格类型下冷却风扇的风量计算结果差异不大.三种网格类型的网格数量和计算时间排序为：多面体网格<切割体网格<四面体网格.因此，选择多面体网格作为后续网格数量无关性研究和分析计算.

网格数量无关性研究则是通过改变旋转域加密尺寸来研究不同网格数量对冷却风扇风量计算结果的影响.图8所示为不同旋转域网格加密尺寸下风扇叶轮局部网格示意图.图9为不同网格数量下冷却风扇风量计算结果.当网格数量超过600万(旋转域网格加密尺寸1 mm，叶尖加密尺寸0.25 mm时，网格数量约为640万)，冷却风扇风量的仿真计算值几乎没变化.

综合网格类型和网格数量的研究结果，考虑到计算资源和计算精度问题，本文选择多面体网格类型进行冷却风扇流体计算域物理模型的网格划分.冷却风扇计算域模型的网格参数控制策略为：非旋转域网格尺寸16 mm，过渡区域网格尺寸8 mm，旋转域网格尺寸1 mm，叶尖加密尺寸0.25 mm，边界层层数6层，边界层增长比1.5，风扇边界层厚度0.4 mm，其余边界层厚度1 mm.

图7 不同网格类型下冷却风扇风量计算结果对比

图8 冷却风扇叶片网格局部图

图9 不同网格数量下冷却风扇风量计算结果对比

5.2 仿真收敛性判定

仿真计算结果收敛与否要通过查看通过STAR-CCM+的Reports功能监测的相关参数量进行判断.本文重点关注的是冷却风扇的风量值，因此监测了冷却风扇进口、出口流量随迭代步数的变化情况及残差曲线图.图10～11所示为风压19.9 Pa下仿真计算文件的残差曲线图和质量流量监测曲线图，从图中可以看出，冷却风扇风量值已经趋于平稳，仿真预测的冷却风量值为180.24 m3/h，与实验测量值间的误差约为3.5%；且各个残差曲线都低于10-6以下，收敛性较好，表明本文采用的网格及每个工况点的计算迭代步数5 000步可以获得一个较准确的仿真预测值.

图10 残差曲线图

图11 质量流量监测曲线图

5.3 冷却风扇内部流场计算结果

5.3.1 速度场分布

图12～13为冷却风扇叶片矢量图和表面速度分布云图.从图中可以看出，空气从半球面入口(虚拟大气入口环境)被冷却风扇吸入后，流经风扇叶片吸风面，经过风扇叶片旋转的作用，气流流动方向发生改变，呈现螺旋形流动；气流随叶片旋转，且在冷却风扇叶尖处存在最大的旋转气流速度.经过风扇叶片旋转加速后的气流流入出口域自由方腔中，最后从出口边界面流出.冷却风扇叶尖旋转气流速度越大，叶片受到的气流冲击效应越大，长时间运行会发生气流冲击腐蚀现象，导致叶片叶尖及其附近区域发生严重变形或腐蚀断裂；同时，叶尖处高压侧气流与低压侧气流的交互作用导致叶尖涡流的形成[31].

(a)风扇叶片吸风面速度矢量图

(b)风扇叶片背风面速度矢量图图12 冷却风扇叶片速度矢量图

(a)风扇叶片吸风面表面速度分布云图

图14表示冷却风扇流体计算域进口处某截面速度矢量图.从图中可以看出，冷却风扇外圈因风扇罩设计的缺陷导致产生严重回流现象，从而会引起冷却风扇工作效率的降低.后续研究可以考虑进行风扇罩的结构优化设计来提升该冷却风扇的运行效率.

图14 冷却风扇外圈速度矢量图

图15为风扇叶片速度局部放大云图及线探针(Line Probe)示意图.从图中可以看出，风扇叶尖的速度较大.Line Probe由30个探测点组成，可以获得叶片上不同点处的速度曲线图(如图16所示).图16中横坐标表示各探测点到坐标轴X轴的垂直距离.图16(a)中表示随着距离的增大，速度基本上呈现出增加的趋势，表明从叶片尾部到叶尖处的速度逐渐增大.同理，图16(b)中可以反映出叶尖处往叶尾处的速度逐渐减小.通过线探针取点监测叶片速度的方法可以反映叶片表面各点的速度.图17中17个监测点的速度分布基本反映了风扇单个叶片大致的速度变化趋势，该方法也常用于温度场计算的取值计算，可为后续对整体流场的实验验证提供对比参考.

图15 风扇叶片速度云图及Line Probe示意图

(a)Line1上各点速度变化曲线图

(b)Line2上各点速度变化曲线图图16 风扇叶片表面线探针监测点速度分布曲线

图17 风扇叶片表面17个监测点速度分布

5.3.2 压力场分布及噪声源预测

图18为冷却风扇叶片压力分布云图.由风扇叶片吸风面的压力云图可以看出，半径较大处的叶片表面的压力数值较大(压力的绝对值)，沿着叶片半径方向，整个叶片的压力呈现出阶梯式的变化.风扇吸风面存在大量的负压区域，在叶片尾缘与风扇外圈交接处，此处负压数值最大.由风扇叶片背风面的压力云图可以看出，在叶尖和叶片边缘处存在部分负压区域，风扇背风面大部分都是正压区域，即风扇的做功区域，风扇叶片前缘的压力最大.

(a)风扇叶片吸风面压力分布云图

(b)风扇叶片背风面压力分布云图图18 冷却风扇叶片压力分布云图

图19为两个线探针上各点压力变化曲线图.图19(a)中大部分监测点的压力都为负值，图19(b)中大部分监测点的压力都为正值，原因是线探针1各监测点位于风扇吸风面，线探针2各监测点位于风扇背风面.

(a)Line 1上各点压力变化曲线图

(b)Line 2上各点压力变化曲线图图19 风扇叶片表面线探针监测点压力分布曲线

图20中17个监测点的压力分布基本反映了风扇单个叶片大致的压力变化趋势，叶尖处及其附近区域压力为负值.冷却风扇叶片旋转压力场存在的压力脉动与固体壁面间产生相互摩擦作用，从而导致风扇在旋转运行中形成离散噪声，该现象称为风扇的旋转噪声[32].风扇叶尖及其附近区域的压力变化较大，产生的旋转噪声越明显.

图20 风扇叶片表面17个监测点压力分布曲线

通过仿真计算已获得冷却风扇稳态压力场，此时，选择声学模块Acoustic中Broadband Noise Sources模型，就可对冷却风扇面噪声源(Curle surface)和体噪声源(Proudman)进行初步预测.图21表示冷却风扇噪声源预测分布云图.从图中可以看出，叶尖处及其附近区域的噪声源较大，在运行过程中产生的离散噪声或旋转噪声就会越明显.

图21 冷却风扇表面噪声源分布云图

5.4 仿真计算结果与实验对比

图22表示不同风压下冷却风扇风量的数值仿真预测值与实验测试结果的对比及误差分析.从图中可以看出，冷却风扇的数值仿真计算结果与实测结果比较吻合.从仿真值与实测值间的相对误差来看，最大误差约为4.9%，总体平均误差为2.94%，表明本文按照实验测试装置等比例搭建的仿真计算模型可以真实可靠地还原冷却风扇的实验测试，且预测的风扇风量值具有一定的精度.验证了本文搭建的冷却风扇风量数值仿真计算域物理模型和计算方法是有效可行的，可以有效地预测冷却风扇内部流场特性.