基于STAR-CCM+的通机消声器护罩温度场数值模拟

谭礼斌, 袁越锦, 刘小强, 黄 灿, 余千英, 杨 东, 唐 琳

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

0 引言

通机(General Purpose Engine)，也称为通用发动机，常指车用及特殊用途以外的发动机[1].通机主要结构包括缸头、缸头盖、箱体、消声器、空滤器、化油器、起动器、风扇、风扇罩及消声器罩等部件.其中，消声器护罩具有防护消声器的作用，以避免消声器损坏，同时防止人与消声器直接接触，避免人体接触后灼伤皮肤；此外，护罩通过与消声器间的间隙及长条形孔或圆形孔进行散热，可以降低护罩表面的温度，因此消声器护罩的合理设计，保证护罩表面温度的均匀分布是非常重要的.消声器护罩一般通过成型钢板冲压制成，相关法规对消声器护罩表面温度具有严格的要求，因此消声器护罩在设计与开发中为了满足国家法规和外观设计的要求，往往需要经过多轮修模或重新开模才能最终得到合理的方案，严重影响了开发进度，浪费人力物力.随着CFD仿真技术的发展，利用数值模拟的方法进行各个工业领域产品的虚拟仿真分析及虚拟开发已经逐渐成为行业的趋势[2-4].基于计算流体力学的思想，采用商业流体分析软件对通机整机产品进行流场分析及消声器护罩辐射传热分析，可迅速获得整机流场分布性能及消声器护罩表面温度分布，提高开发成功率并减小实验次数，对优化消声器护罩设计、提升产品性能都具有重要的指导意义.

对于消声器的数值模拟和辐射传热领域的研究，许多科研工作者采用ANSYS、Fluent、Virtual Lab等软件进行了深入的研究[5-7].如张智辉等[8]采用Fluent软件对消声器内部流场进行了数值模拟分析，评估消声器流场对其消声性能的影响.张磊[9]采用数值模拟的思想，应用Matlab编程研究了消声器穿孔管结构参数对消声器性能的影响，为消声器的设计与优化提供了思路.姜晓光等[10]采用Fluent对柴油机缸内辐射换热特性进行了分析，验证了辐射换热计算的可行性.目前，用于流体分析的软件主要有Fluent、CFX、STAR-CCM+等，其中，意大利西门子公司旗下的通用流体分析软件STAR-CCM+具有较高的集成度，是真正的几何处理、网格划分、求解器、后处理完全一体化环境的软件，高度集成，且具有自带的多面体网格技术，可直接生成高质量多面体.多面体在促进计算收敛及求解精度方面都有较好的优势.STAR-CCM+具有Thin Mesher网格技术(薄壁层网格)，对流固耦合等共轭换热问题，能够自动探测固体薄件，划分为高质量、均匀分布的类似于棱柱一样的网格，同时能大大控制网格数量，保证能量方程计算能够得到很好的收敛.因此，在汽车、发动机、通用机械等领域，STAR-CCM+得到了业界广泛的使用及认可.

本文以某型通机整机为研究对象，基于STAR-CCM+ 11.06流体仿真分析平台，对该型号通机不同起动器状态下的整机流场及消声器护罩辐射传热进行数值模拟分析，获取消声器护罩表面的温度分布，与实验进行对标验证，建立消声器护罩辐射传热分析方法，为后续消声器护罩前期设计与开发提供分析支撑和理论依据.

1 物理模型

某型通机(三种不同起动器状态)的整机三维模型采用Unigraphics NX 8.0依据实物按照1∶1比例(Siemens PLM Software Company，Italy)建模获得(如图1所示).采用STAR-CCM+ 11.06(Siemens PLM Software Company，Italy)对整机三维模型进行几何清理并网格划分，获得通机整机的网格模型(如图2所示)和通机整机外流场计算域模型(如图3所示).三种起动器状态下的通机网格模型采用相同的网格控制策略进行网格划分，采用相同的边界条件设置进行整机消声器护罩的辐射传热分析.图4表示消声器护罩表面提取的13个温度测量点(实验和仿真测量点).常用的温度测量方法有红外热成像仪、点温计、温度传感器等[11,12]，本文选用点温计进行消声器护罩温度的测量.该通机三种起动器状态的测试环境为室温28 ℃，通机空载运行，风扇转速为3 000 RPM.

(c)更换起动器和消声器护罩(C)图1 整机三维模型

图2 整机共形网格模型(蓝色显示的为缸 体、箱体及消声器部件等处理为固体部件)

图3 整机外流场计算模型

(a)A型和B型整机消声器护罩温度测量点

(b)C型整机消声器护罩温度测量点图4 整机消声器护罩测量点 (13个)示意图

2 数学模型

流体流动满足三大守恒定律：质量守恒定律、动量守恒定律和能量守恒定律[13].在流体流动处于湍流状态时，整个体系还要遵循湍流运输方程.以上这些守恒定律的数学描述，统称为控制方程.文中选用STAR-CCM+中提供的Realizable k-ε湍流模型进行数值计算.

湍流控制方程为三维不可压缩雷诺时均Navier-Stokes方程：

(1)连续方程[14]

(1)

(2)动量方程(N-S方程)[15]

(2)

(3)能量方程[16]

(3)

式(1)～(3)中：ui、uj是平均速度分量，m/s，xi、xj为坐标分量，m；p是流体微元体上的压力，Pa；μeff是湍流有效黏性系数，Pa·s；T是温度，K；λ为流体换热系数，W/m2·K；Cp为流体比热容，J/(kg·K)；ST是流体内热源和由粘性作用引起流体机械能转变为热能，J.

(4)k-ε湍流模型方程[17]

(4)

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

(5)辐射传热模型

本文采用STAR-CCM+ 11.06进行消声器护罩的辐射传热分析.辐射传热模型选择为Surface to Surface Radiation，类型为Gray Thermal Radiation[18].表面A1和表面A2之间的辐射换热能量的计算公式为[19]：

(5)

式(5)中：φb12为辐射能量，W；X12、X21为角系数；A1、A2为表面积，m2；Eb1、Eb2为辐射能力，W/m2.

通机整机固体计算域模型中任意两个表面，若其表面面积和温度作为已知条件，那么角系数一旦确定，就可以求出辐射换热能量.因此问题的关键为求解并确定角系数[19].STAR-CCM+中可以通过View factor Calculator计算角系数，通过固体热边界条件(温度和换热系数)的输入设置即可计算受热辐射部件的表面温度，获取固体部件的辐射温度场.

3 模型求解

本文某型通机整机消声器护罩辐射传热分析采用稳态分析，湍流模型选择为Realizablek-ε湍流模型，计算域的壁面参数采用STAR-CCM+推荐的Two-Layer All y+Wall Treatment 设置.流体计算域网格采用多面体网格和边界层网格(3层)模型进行网格划分，固体计算域网格采用多面体网格和薄壁层网格(5层)模型进行网格划分，网格细化和网格控制策略采用已建立的网格参数细化方法进行设置[1].生成的流体计算域网格约为550万，固体域计算网格约100万.消声器护罩辐射传热分析的相关边界条件设置如下：

(1)旋转域：风扇旋转速度为3 000 RPM，处理为Moving Reference Frame(MRF)；

(2)流体域：消声器入口流量为1.016 g/s(实测数据)(质量流量入口)，入口温度为280 ℃；空滤器出口流量为1.016 g/s，出口温度为40 ℃；虚拟计算域入口设置为滞止入口(Stagnation Inlet)，出口设置为压力出口(Pressure Outlet)；实验测试环境温度为28 ℃；流体属性选择为标准大气压下理想气体(Ideal Gas)，定压比热容为1 003.62 J/kg·K，导热系数为0.027 W/m·K，动力粘度为1.85×10-5Pa·s.

(3)辐射边界：辐射边界(发射率、吸收率等)可以在流体侧加载，也可加载于固体侧进行辐射分析.本文辐射边界加载于流体侧.缸体、箱体的辐射发射率设置为0.5；消声器各部件的材料为ST14，辐射发射率设置为0.5.依据《传热学》中“大部分非金属材料辐射发射率都很高，一般在0.85～0.95之间，且与表面状况(包括颜色在内)的关系不大，在缺乏资料时，可近似地取作0.9″，因此，本文在进行整机消声器护罩辐射传热计算中，油箱、风扇罩、空滤器等非金属材料的辐射发射率取值为0.9[20].

(4)固体域：辐射传热分析中缸体、箱体、消声器部件等处理为固体，缸体、箱体的材料为ADC12，密度为2 800 kg/m3，定压比热容为880 J/kg·K，导热系数为190 W/m·K；消声器各部件的材料为ST14，密度为7 830 kg/m3，定压比热容为434 J/kg·K，导热系数为65 W/m·K.固体计算域热边界设置依据燃烧分析结果进行设定.其中，缸体、箱体内部表面(internal surface)热边界设置为：温度80 ℃，对流换热系数100 W/m2·K；缸套及燃烧面的区域划分如图5所示，缸套三个区(Cylinder-up、Cylinder-middle、Cylinder-down)加载的温度分别为200 ℃、225 ℃、250 ℃，对流换热系数都设置为300 W/m2·K；进排气道(Intake Pipe、Exhaust Pipe)、燃烧室面四个区(火花塞Sparkplug、燃烧室Zone1、燃烧室Zone2、燃烧室Zone3)按空载运行状态下加载温度及对流换热系数，加载的温度分别为100 ℃、500 ℃、750 ℃、650 ℃、600 ℃、550 ℃，对应的对流换热系数分别为100 W/m2·K、500 W/m2·K、700 W/m2·K、700 W/m2·K、650 W/m2·K、650 W/m2·K.

图5 热边界加载区域示意图

4 模拟结果及分析

4.1 共形网格与非共形网格计算结果对比

图6表示共形网格模型和非共形网格模型下消声器护罩各测量点温度变化曲线图.将流体计算域模型和固体计算域模型一起划分网格，设置相同的网格基准尺寸和最小尺寸，采用自定义控制对流体域侧禁用薄壁层网格模型(Part Control)，固体域禁用边界层网格模型(Surface Control)，网格划分完成后即形成共形网格.共形网格和非共形网格作为流固耦合计算的两种网格类型，在网格节点上进行数据交互或插值计算时，共形网格模型下，固体计算域的网格与流体计算域的网格完全重合，两者间可以通过重合的网格单元进行更好的数据交互.从图6可以看出，共形网格模型下消声器护罩各测量点温度仿真值与实验值间的误差较小，表明共形网格模型下的温度计算精度更高.

(a)A型起动器下消声器护罩温度变化曲线图

(b)B型起动器下消声器护罩温度变化曲线图

(c)C型起动器下消声器护罩温度变化曲线图图6 消声器护罩测量点温度变化曲线

4.2 边界条件设置对消声器护罩温度计算结果的影响

4.2.1 空滤器流量、消声器流量及消声器出口温度

图7表示不同流量及消声器出口温度状态下消声器护罩测量点温度变化曲线图.空滤器及消声器流量1.016 g/s为实测边界条件换算的流量值，消声器出口温度280 ℃为实测温度值.从图7中可以看出，实测状态下的消声器护罩温度仿真值与护罩温度实验值的变化趋势基本一致.流量及出口温度变化，对消声器护罩温度仿真结果具有一定的影响，流量及出口温度增大，消声器护罩温度仿真值增大.

图7 流量及排温对计算结果的影响

4.2.2 材料辐射发射率的影响

图8表示材料不同辐射发射率下消声器护罩测量点温度变化曲线图.从图8中可以看出，材料辐射发射率减小，消声器护罩温度降低，材料辐射发射率增大，消声器护罩温度增大.整体上分析，三种辐射发射率条件下获得的消声器护罩仿真温度值间存在的差异不大，原因是通机整机的温度不大，通过辐射作用产生的热量就不会存在较明显的差异.本文采用的材料辐射发射率值(固体部件)都为0.5，该辐射发射率下的消声器护罩温度仿真结果间的误差较小.

图8 材料辐射发射率对计算结果的影响

4.2.3 缸头热边界的影响

图9表示不同缸头热边界设置状态下(缸体中缸套及其燃烧室区域的温度边界及换热系数边界)消声器护罩测量点温度变化曲线图.从图9中可以看出，依据实际热边界进行数值模拟获得的消声器护罩温度结果与实验值间的误差较小.热边界温度减小，辐射传热过程中形成的辐射能量传递减小，消声器护罩温度值降低；热边界温度增大，则消声器护罩温度值增大.

图9 缸头热边界设置对计算结果的影响

4.2.4 重力模型的影响

该型通机整机为垂直轴发动机，因此在计算过程中可以考虑重力模型对消声器护罩温度计算结果的影响.图10表示重力模型设置状态下消声器护罩测量点温度变化曲线图.从图10中可以看出，重力模型的设置对消声器护罩温度值没有显著的影响.因此，后续分析忽略重力模型对消声器护罩温度值的影响.

图10 重力模型对计算结果的影响

4.2.5 边界层网格及薄壁层网格的影响

实际求解计算过程中，网格的细化及其网格质量对求解精度具有一定的影响.本文网格划分已按照建立的通机网格参数控制策略进行相应的网格细化，具有较好的求解精度.而鉴于温度场的流固耦合计算问题中，边界层网格和薄壁层网格对整机流固耦合温度场计算求解精度和效率都有一定的影响，因此在已有的网格控制策略下探究边界层网格及薄壁层网格层数对计算结果的影响是非常重要的[21,22].
图11表示模型边界层网格层数与薄壁层网格层数对消声器护罩温度计算结果的影响.从图11中可以看出，不同薄壁层网格层数下消声器护罩温度仿真值基本一致；边界层网格层数对消声器护罩温度的计算结果具有显著的影响，计算域模型的边界层层数为6层时消声器护罩温度仿真值与实验值更吻合，误差减小.因此，后续分析中采用边界层网格层数6层、薄壁层网格层数5层进行分析，减小仿真值与实验值间的误差.

(a)6层边界层网格局部放大图

(b)8层薄壁层网格局部放大图

(c)消声器护罩各测量点温度变化曲线图11 边界层及薄壁层网格 对计算结果的影响

4.2.6 综合边界因素的影响

图12表示流量为空滤器及消声器流量为1.016 g/s，消声器出口温度280 ℃，边界层网格层数6层，薄壁层网格层数5层，流体域和固体域网格共形下消声器护罩温度仿真值与实验值的对比曲线图.从图12中可以看出，调试后的仿真计算值与实验值的变化趋势更接近，两者间的误差减小.因此，将该设置作为基准设置，进行三种型号起动器状态的(A型、B型、C型)通机消声器护罩辐射传热分析.

图12 消声器护罩温度仿真值 与实验值对比曲线

4.3 整机流场分析计算结果

图13表示不同起动器状态下整机表面风速分布云图.从图13中可以看出，三种起动器状态下的整机表面风速分布趋势基本一致.

(a)A型起动器下整机表面风速分布图

(b)B型起动器下整机表面风速分布图

(c)C型起动器下整机表面风速分布图图13 整机表面风速分布图

4.4 辐射传热分析计算结果

4.4.1 整机温度场分布

图14为不同起动器状态下整机固体部件温度分布云图.三种起动器状态下的固体部件温度场分布趋势基本一致.消声器护罩前侧(正对冷却风扇侧)的温度降低，则背侧(背离冷却风扇侧)的温度较高.造成这种现象的原因是风扇吹出的冷却风较多地冲刷在消声器护罩前侧上，具有较好的冷却作用，而消声器护罩靠近火花塞高温部件，且此侧的冷却风较少，未得到较好的冷却.A型整机火花塞温度仿真结果为180 ℃，实际测量温度为175 ℃，仿真值与实验值相差较小.

(a)A型起动器下整机温度分布图

(b)B型起动器下整机温度分布图

(c)C型起动器下整机温度分布图图14 整机温度分布云图

4.4.2 消声器护罩温度场分布

图15为不同起动器状态下消声器护罩温度分布云图.从图15中可以看出，不同起动器状态下的消声器护罩温度分布趋势基本相同.消声器护罩正对风扇的那侧温度较低，远离风扇侧温度较高.原因是正对风扇侧具有较多的冷却风冲刷在消声器护罩上，具有一定的冷却作用.三种起动器状态下消声器护罩平均温度在60 ℃～70 ℃之间，温度分布合理，表明该消声器护罩的设计可行.

(a)A起动器状态下护罩温度分布云图

(b)B起动器状态下护罩温度分布云图

(c)C起动器状态下护罩温度分布云图图15 消声器护罩温度分布云图

4.4.3 消声器护罩辐射能量场分布

图16为不同起动器状态下消声器护罩辐射能量分布云图.从图16中可以看出，不同起动器状态下的消声器护罩辐射能量分布趋势基本一致.消声器护罩远离风扇侧的辐射能量较大，在温度场中此侧对应的温度较高.

(a)A起动器状态下护罩辐射能量分布云图

(b)B起动器状态下护罩辐射能量分布云图

(c)C起动器状态下护罩辐射能量分布云图图16 消声器护罩辐射能量分布云图

4.5 实验对标验证

图17为不同起动器状态下消声器护罩各测量点温度仿真值与实验温度值的对比曲线图.从图17中可以看出，消声器护罩温度仿真值与护罩温度实验值的变化趋势基本一致，表明本文数值模拟计算的护罩温度值具有一定的可行性，模拟结果可以有效地预测消声器护罩表面温度分布情况，评估温度值的合理性，可为消声器护罩前期的设计及开发提供分析支撑，减小实验次数及开发成本.

图17中，温度仿真值与实验值具有一定的误差，产生误差的原因可能是三维模型数据的简化处理(未考虑螺钉、线束、小型板件等)、材料辐射发射率的选取、稳态计算结果差异、仿真温度测量点的提取与实验测量点的差异等.实际情况下，材料辐射发射率率是随着温度、材料表面状态及颜色是变化的，整个计算过程依据实际运行工况应为整机系统级瞬态热管理分析计算，且应尽可能将通机整体部件都考虑成固体部件，从而进行整机热管理分析计算(辐射温度场分析)，使数值仿真结果与实验测试结果更接近.

(a)A型起动器下护罩实验温度值与模拟值对比

(b)B型起动器下护罩实验温度值与模拟值对比

(c)C型起动器下护罩实验温度值与模拟值对比图17 某型整机消声器护罩实验 温度值与模拟值对比曲线图

5 结论

(1)采用STAR-CCM+对通机消声器护罩进行温度场的数值模拟，与点温计测量的消声器护罩表面温度的实验值进行了对比研究，验证了本文采用的辐射传热分析方法具有一定的可行性，可有效地预测消声器护罩表面的温度分布情况.

(2)共形网格具有较好的计算求解精度；依据消声器护罩温度仿真值与实验值间的对比分析，得出材料辐射反射率、热边界条件、边界层网格对消声器护罩温度的计算结果具有显著的影响；薄壁层网格、重力模型对温度仿真结果的影响较小.

(3)不同起动器状态下消声器护罩温度和辐射能量分布趋势基本相同，远离风扇侧的辐射能量大，温度较高.消声器护罩温度仿真值与实验值产生误差的原因主要是材料辐射发射率的选取、几何模型的简化及求解模式(稳态)的选取.

(4)采用计算流体力学的数值模拟方法可以有效预测消声器护罩表面温度分布，评估温度值的合理性，可为消声器护罩的设计与开发提供分析支撑，减小实验次数，降低开发成本，达到温度法规要求.