小型强化内燃机燃烧室爆震流场及波致疲劳模拟

王志坚， 张佳明， 吕文芝， 黄国龙， 孟庆斌

(1.内燃机可靠性国家重点实验室，山东 潍坊 261061； 2.潍柴动力股份有限公司，山东 潍坊 261061； 3.北京科技大学 机械工程学院，北京 100083； 4.教育部流体与材料相互作用重点实验室，北京 100083)

随着环境和能源问题在世界范围得到广泛关注，以及越来越严格的排放法规限制，内燃机具备了高效和低污染的特性，内燃机小型化技术得到了越来越多的认同[1-4]。内燃机缸内燃烧过程的不稳定性增加，爆震倾向增加，爆震发生时最大缸压可达20 MPa以上，远高于常规工况，降低燃油经济性的同时内燃机机构件产生严重破坏。活塞处于一个高速、高压和高温的恶劣工作环境，在这种高频、高振幅的压力震荡激励下，往往会导致活塞失效，因此研究内燃机爆震发生特性及活塞强度校核对小型强化内燃机的安全性及可靠性有着重要的理论价值和工程价值[5-7]。对爆震特性的研究主要通过仿真和试验的手段，爆震发生时燃烧室内的流场规律，对活塞强度的研究主要针对热冲击、热疲劳以及机械疲劳等方面[8-12]，因此对于小型强化后缸压激励作用以及爆震发生后缸压激励作用下活塞的受力规律和强度校核研究不足，很少有学者研究爆震激励下活塞的寿命[13-14]。

本文对缸内发生爆震时的冲击波进行研究，将爆震压力作为激励施加到活塞顶面，分析活塞在不同工况下的温度场、热-机耦合强度以及疲劳寿命，阐明爆震形成活塞破坏性作用的机理，同时为小型强化内燃机活塞的设计和分析提供理论基础和工程依据。

1 爆震数值模拟原理与实验验证

1.1 小型强化内燃机热力学模型

在内燃机中，活塞顶面、汽缸壁面、缸盖底面和进排气门共同围成燃烧室，将该模型简化为一个热力学系统进行研究。缸内工质的状态由能量守恒方程、质量守恒方程及理想气体状态方程把整个工作过程联系起来。

气体状态方程、质量守恒方程及能量守恒方程分别为：

pV=mRT

(1)

(2)

(3)

(4)

式中：p为缸内气体压力，Pa；V为气缸工作容积，m3；m为缸内气体质量，kg；R为气体常数，J/(mol·K)；T为缸内气体温度，K；mb为缸内已燃烧的质量，kg；me为流入气缸内的空气质量，kg；ma为流出气缸的质量，kg；Qb为燃烧放出的热量，kJ；Qw为壁面传热量，kJ；h为气缸内工质的比焓，kJ/kg；h2为进气门前比焓，kJ/kg；Cv为比定容热容，kJ/(kg/K)；αm为瞬时过量空气系数；φ为曲轴转角。

联立式(2)～(4)可得：

(5)

由此可求出气缸内工质的温度变化，再通过式(1)气体状态方程便可求出缸内压力变化。

1.2 增压模拟及爆震特征分析

增压技术的基本原理是通过提高进气密度，来提高进入燃烧室的空气量，同时提高每个循环的供油量，使缸内燃烧更为剧烈，以此提高缸内的平均有效压力，从而提高发动机输出功率。

通过改变气缸内的空燃比可以模拟内燃机的增压过程，将自然进气工况缸内空燃比设置为14∶1，一级增压工况为20∶1，二级增压工况为28∶1，利用式(1)、(5)得到缸内压力变化曲线如图1(a)所示，自然吸气、一级增压、二级增压3种工况下时缸压最大压力约为6、9和12 MPa。随着小型强化内燃机强化程度的提高，爆震的发生几率也随之提高。如图1(b)所示，当爆震发生时，缸内压力首先产生一个巨大的振荡，振幅约10 MPa，紧接着是振幅逐渐衰减的压力波动，振荡频率约10 kHz，最后缸压趋于平稳，可见，该爆震波作用于燃烧室零部件上，会对构件疲劳性能产生影响。因此，通过对二级增压时内燃机最大压力值赋予振荡衰减函数来得到爆震缸压曲线。

图1 不同工况下缸内压力曲线Fig.1 Cylinder pressure curves under different conditions

1.3 爆震实验验证

采用高温电容式压力传感器检测爆震激励信号，其工作原理为由传感器膜片直接感知缸内压力，在膜片与介质之间形成微小的位移，该位移随着被测压力的增加而增大，而传感器的电阻也随之变化，高温电容式压力传感器的结构如图2所示。

图2 高温压力传感器结构示意Fig.2 Structure of high-temperature capacitive pressure sensor

利用测量电路来检测电阻的变化，并产生与爆震压力相对应的输出信号。由于压电晶体的原始输出信号非常微弱，因此将电荷放大器与压力传感器共同放置在气缸顶部以增强信号。然后再并联上旋转编码器，捕捉50个周期的缸压曲线，即为爆震激励信号，图3为测量流程图。

图3 爆震激励压力测量流程Fig.3 knock excitation pressure measurement process

实验采用RoverK16型汽油机，气缸容积为1 397 cm3，冲程79 mm。实验得到小型强化内燃机爆震时缸内压力随时间的变化，每隔3个循环出现一次爆震现象，如图4(a)所示。

图4 缸内压力与爆震单周期缸内压力曲线Fig.4 Cylinder pressure and knock single-cycle cylinder pressure with the time diagram

由图4(b)爆震单周期曲线可知，当爆震发生时，首先产生一个巨大的振荡，振幅可达近10 MPa，缸内压力最大可以达到18 MPa，紧接着是振幅逐渐衰减的小振幅高频率振荡，最后缸压趋于平稳，该规律与数值模拟得到的爆震缸压曲线一致，因此数值模拟得到的缸内压力测试数据得到验证。

2 爆震仿真分析

2.1 有限元模型

有限元模型的建立分为活塞、缸盖及汽缸套部分，综合考虑网格划分及实际计算，对模型的火花塞、气门以及部分倒角进行了简化处理，其他部分严格按照图纸进行三维建模。活塞外径为74 mm，高度为66.1 mm，销孔直径为15 mm，材料为灰铸铁，热膨胀系数为1.1×10-5/℃。实体模型如图5所示。

图5 小型强化内燃机活塞三维实体模型Fig.5 3D solid model of downsized engine cylinder

流体区域网格划分共生成36 431个节点，197 375个单元。固体区域的活塞为中心对称结构，因此取1/4模型进行研究，固体区域网格划分共生成388 474个节点，246 345个单元。网格划分完成后导入ICEM CED 平台进行Quality检验，检验结果显示Quality值均在0.2以上，满足应用要求。流体区域及固体区域网格划分如图6所示。

图6 小型强化内燃机活塞有限元模型网格划分Fig.6 Meshing of downsized engine cylinder finite element model

2.2 边界条件及求解模型

为了使研究爆震冲击波的传播过程计算简便，本文使用欧拉方程作为数值模拟的控制方程。采用FLUENT对缸内爆发压力进行CFD模拟计算，由于研究的是瞬态的、高速可压流体的流动问题，所以采用的数值模拟方法是基于密度的、瞬态的方法。又由于要捕捉冲击波这种非连续性的流场，因此选择使用三阶的MUSCL计算格式。

本文借鉴文献[7]的爆震模拟方法，应用破膜法模拟燃烧室中爆震冲击波的产生，即在燃烧室自燃发生区域创造高压区，使其与周围环境形成强大压差，从而形成冲击波。研究主要针对爆震初期压力剧烈波动的时段，这段持续时间很短，只有3°曲轴转角左右，因此在该时段内假设避免绝热，活塞静止不动。研究的锥顶形燃烧室如图7所示，缸内压力赋值为压缩行程终了时的6 MPa，在燃烧发生的高压区赋值为100 MPa；时间补偿估算为：

Δt=l/c

(6)

式中：l为网格的特征尺度；c为冲击波的波速，约为2 000 m/s。因此估算其时间步长约为10-8s。

图7 锥顶形燃烧室模型Fig.7 Cone-top combustion chamber model

在模拟开始之前，在燃烧室底部从中心点向外缘依次取点作为监测点，用于观察不同位置的超压分布，测点位置如图8所示。

图8 监测点位置示意Fig.8 Monitoring point location diagram

2.3 爆震结果分析

爆震发生初期的几次压力振荡对燃烧室零部件的破坏最为严重，该阶段持续期仅为3°曲轴转角左右，对应时间为0.6 ms，因此将模拟的时间设置为0.85 ms，可以保证监测到燃烧室内压力振荡较为剧烈的阶段。燃烧室底面对应活塞顶面，在燃烧室底面选取11个监测点对应活塞顶面10个圆环位置，活塞顶面圆环位置示意如图9所示，监测点8～11的压力波动曲线如图10所示。

图9 活塞顶面圆环位置示意Fig.9 Piston top ring position diagram

如图10(a)所示，在0.06 ms时刻，位置8出现一个明显的高幅值压力振荡，幅值达到了33 MPa。如图10(b)～(d)所示，位置9、10、11的压力最大值依次达到了40、70、210 MPa，而且可以明显观察到4次振幅逐渐减小的压力波动，位置11压力最大值在210 MPa左右，远远超过了燃烧室中心高压区域的初始赋值100 MPa。爆震冲击波抵达外缘狭小区域后，由于空间缩小，缸盖与活塞顶部之间的间隙减小，爆震冲击波在缸盖与活塞之间的反射次数也成倍上升，因而导致在端部尖角处的振荡频率急剧增大，进而导在活塞外缘的压力急剧上升，进而在活塞外缘出产生剧烈的压力振荡。同一时刻不同位置的压力值不同，从活塞顶面圆心处向外缘依次产生压力振荡，这符合爆震波在燃烧室内的传播规律。

图10 不同位置的压力波动曲线Fig.10 Pressure fluctuation curves of different positions

由图11可见，随着监测位置从圆心向外缘移动，压力最大值先小幅减小然后大幅升高，并在位置11达到最大。

图11 不同监测位置最大值Fig.11 Maximum pressure at different monitoring positions

2.4 爆震热机耦合强度以分析

1)活塞热-机耦合边界条件。

内燃机运行时活塞主要受到热载荷和机械载荷作用，热载荷使活塞形成温度梯度，从而产生不同的位移形变及热应力应变，机械载荷主要考虑活塞在燃气爆发压力下的载荷，本文采用热机耦合模块计算耦合工况下的活塞强度。

内燃机稳定运行后，活塞温度场基本保持不变，因此采用以对流换热系数作为边界的第3类热边界条件来计算活塞温度场，如图12所示。由于本文主要研究爆震初期破坏性最强的几次反射，在该时间段内，活塞在上止点附近速度较慢，可以忽略活塞的运动，因此通过约束活塞销孔部位来模拟活塞销对活塞的约束力，具体方法是在活塞销孔上施加一个法向约束，限制其轴向运动，其他方向运动自由。活塞温度分布云图如图13所示，活塞顶面中心温度最高为358 ℃，温度由中心向边缘逐渐降低。活塞头部沿轴向出现温度分层，第1环岸与第2环岸间的温差约为50 ℃，第2环岸与第3环岸间的温差约为35 ℃，第3环岸与第4环岸间的温差约为30 ℃，活塞裙部包括活塞内腔温度为102 ℃，活塞顶部与底部温度相差约256 ℃。

基于上述原因,本文提出了一个线性卡尔曼滤波器.该滤波器与MDS相结合来跟踪传感器网络节点的位置,给出了一种在测量对应的成对距离时寻找节点变化位置的有效方法,并据此设计了KL滤波器,以实现对室内机器人在传感器网络下的精确定位.

图12 活塞热机耦合边界条件Fig.12 Piston thermo-mechanical coupling boundary conditions

图13 活塞温度分布云图Fig.13 Temperature contours of piston

2)活塞顶面强度。

以热-机耦合边界条件为基础，将数值模拟及仿真得到的气缸压力施加到活塞顶面，以活塞顶面中心点为圆心生成10个同心圆环，将压力波动曲线施加在对应圆环区域内，即可得到不同工况下活塞的应力和变形，活塞顶面11个监测点的最大值应力和最大变形值如图14所示。应力值由活塞顶面中心向外缘先平稳升高然后逐渐下降，各个监测点的应力值随增压强度的提升不断增大，爆震发生时产生的应力值明显高于其他工况的应力值，但变化趋势基本不变。由图14(b)可知，活塞在自然吸气、一级增压、二级增压工况时最大变形分布趋势基本一致，均为逐渐减小然后趋于平稳，且变形值随增压强度的提升不断增大，爆震发生时产生的变形值明显高于其他工况，且变化趋势与其他工况不同，为先减小后逐渐升高达到最大值，分析原因为爆震发生时产生的压力并非均匀施加到活塞顶面，而是活塞顶面边缘处受到的压力最大，靠近顶面中心处压力变小。

图14 活塞顶面不同位置不同工况的应力值和变形值Fig.14 Stress and deformation values for different working conditions at different positions on the top of piston

3)活塞整体强度。

爆震的发生不仅对活塞顶面的应力及变形产生影响，而且对活塞整体产生影响，图15～18为活塞在不同工况下的应力和变形云图，将不同工况下的最大拉应力、最大压应力、最大变形数据及发生的位置如表1所示。最大变形位置为活塞裙部底端外侧。

图15 自然吸气工况活塞应力及变形云图Fig.15 Stress and deformation contours of piston naturally aspirated condition

图16 一级增压工况活塞应力及变形云图Fig.16 Stress and deformation contours of piston stepped boom condition

图17 二级增压工况活塞应力及变形云图Fig.17 Stress and deformation contours of piston secondary boost condition

图18 爆震工况活塞应力及变形云图Fig.18 Stress and deformation contours of piston knocking condition

表1 不同工况下活塞热机耦合仿真结果Table 1 Thermal-mechanical coupling simulation results of piston under different working conditions

由图15～18所示应力位移云图数据所示，自然吸气、一级增压、二级增压3种工况下最大拉应力值逐渐增大，且均出现在活塞内腔顶部中心处，分析由于活塞头部直接受到缸内压力作用后向内凹陷导致。3种工况下最大压应力值逐渐增大，均出现在活塞内腔与活塞头部交界处，分析由于活塞头部中心即薄弱区域向内凹陷后形成的应力集中区域。爆震工况下最大拉应力为195 MPa，出现在活塞销孔与活塞内腔交界处，最大压应力为118 MPa，出现在活塞环第1环的内侧，活塞应力与变形变化曲线与爆震压力波动曲线基本一致，如图19所示。4种工况下主应力均小于活塞材料灰铸铁的屈服极限280 MPa。

4种工况下最大变形值逐渐增大，均出现在活塞裙部底端外侧，分析由于该位置为活塞薄弱区域，刚度较小，受力之后容易产生变形。爆震工况下最大变形值为0.07 mm，没有超过活塞与气缸套的标准间隙0.2 mm，符合间隙要求。

综上所述，爆震的发生一方面大幅度增大了活塞顶面各个监测点、活塞内部的应力值和变形值，另一方面改变了活塞最大应力位置，即危险点位置发生了改变。活塞模型的强度分析结果显示4种工况下活塞所受热机耦合应力没有超过材料的屈服极限，且变形结果没有超过标准间隙，即满足应用要求。

2.5 爆震波致疲劳仿真分析

爆震发生时燃烧室内产生强烈的爆震波，该爆震波为一种冲击波，作用在活塞上产生应力与变形，冲击波对活塞施加反复的高频的冲击载荷。当活塞受到多次重复冲击载荷作用时，因累积损伤会形成裂纹萌生、扩展直至断裂破坏，该过程为活塞受波致疲劳的基本过程，在该过程中因受到冲击载荷作用而产生的应力与应变为冲击应力与应变。

图19 活塞爆震时应力及变形曲线Fig.19 Deformation and stress curve of piston during knocking

本文使用nCode软件进行疲劳寿命分析，将有限元分析的应力应变场结果导入，材料疲劳特性参数由Ansys-Workbench数据库中导入，然后对载荷谱曲线进行峰谷值抽取，进行雨流计数，设置疲劳寿命预测方法为名义应力法，应力组合方法为SignedVonMise，应力修正方法为Doodman。活塞疲劳寿命计算结果如图20所示。

按照平均时速100 km/h计算，疲劳寿命云图显示循环次数均在4.43×107次以上，最先发生破坏的位置在活塞顶面中心、活塞内腔交界处以及第1气环处，假设爆震每隔3个正常循环发生一次，最终计算得到发动机可连续运行249 398 km，达到规定里程22万km。

3 结论

1)通过数值模拟的方法，得到内燃机在自然吸气、一级增压、二级增压以及爆震时的缸压曲线，得出爆震缸压波动特征，并得到了爆震压力测量实验的验证。

2)爆震冲击波使活塞外缘的压力急剧上升，活塞顶面外缘压力最高，达到210 MPa，远超燃烧室中心高压区域的初始赋值100 MPa。

3)爆震工况下活塞最大拉、压应力分别为195 MPa和118 MPa，较自然吸气、一级和二级增压工况明显增大，同时爆震改变了活塞最大应力的位置。

4)爆震波致疲劳仿真分析显示，爆震工况下最先发生破坏的位置为活塞顶面中心、活塞内腔交界处、第1气环处，计算得到内燃机可以连续运行249 398 km，达到规定里程22万km。