柴油机缸内三维动网格和流场数值模拟程序开发

刘晓刚，明平剑，张格健，张文平

(哈尔滨工程大学 动力与能源工程学院，黑龙江 哈尔滨 150001)



柴油机缸内三维动网格和流场数值模拟程序开发

刘晓刚，明平剑，张格健，张文平

(哈尔滨工程大学 动力与能源工程学院，黑龙江 哈尔滨 150001)

现较为成熟的CFD仿真软件均为国外所有，国内缺乏相应的知识产权。为此，本文自主研发了一套仿真程序，对内燃机缸内三维工作过程进行数值模拟。本文建立了动态层网格模型和匹配的任意边界模型，采用有限体积方法离散流场控制方程，通过SIMPLEC算法求解NS方程，基于FORTRAN语言开发了一套内燃机缸内流场仿真软件。对内燃机的缸内空气纯压缩和膨胀过程做了仿真分析，计算结果表明，缸压曲线、缸内温度曲线与应用商业软件PUMPLINX、FLUENT计算结果基本吻合,证明动网格模型和流场数值模拟方法正确，本文的程序有进一步市场应用的前景。

内燃机；运动网格；CFD；数值模拟；程序开发；软件仿真

内燃机缸内流场数值模拟主要是借助CFD手段预测缸内的流动形式，从而在设计阶段可预先改进进气形式或燃烧室形状，形成更好的缸内气流状态以提高燃烧效率和降低污染物排放，在现代内燃机设计和预先研究中发挥了越来越重要的作用[1]。而目前国内的内燃机设计上主要采用市面上较为流行的商用软件，如KIVA系列软件[2]、AVL Fire软件、STAR-CD软件、ANSYS Fluent软件和新兴的CONVERGE软件等，购买和升级费用昂贵，没有独立的知识产权，一直受制于人，阻碍了中国内燃机行业的发展。

哈尔滨工程大学自主研发的基于非结构化网格的通用输运方程求解软件[3-5](general transport equation analysis,GTEA)采用有限体积方法离散流场控制方程，通过SIMPLEC算法求解NS方程和分离式变量求解方法，目前可以较为精确的求解二维和三维、层流和湍流的流场流动问题。本文在该软件的基础上，根据内燃机的气缸往复直线运动形式建立了动态层网格模型和匹配的任意边界模型，并实现了其与流场求解器的耦合计算，形成了一套内燃机缸内流场专业化仿真软件。

1 数学模型

数值计算中数学模型可分为计算域的几何信息模型和数学控制方程模型。

1.1 运动网格模型

动态层模型[5]是利用增加或删减网格层来适应边界的运动，该方法对网格的划分及边界的运动形式要求较高，边界运动主要为往复式直线运动，网格需要在动边界运动方向上分层。程序实现时首先要建立结构体数组来保存每层网格的单元和节点信息，以便进行网格层间拓扑信息的查找，之后再根据计算域是膨胀还是压缩来做网格层的加减，如下图1所示，只对最下层的网格做加减层操作，以最大程度保持网格原有质量。通过网格节点坐标替换及网格标记阻断相结合的方法来满足增减网格层时网格拓扑结构的变化，减层时1、2、3、4单元被阻断，也就是不参与计算，F、G、H、I、J节点替换A、B、C、D、E节点坐标；加层时增加1、2、3、4单元，增加A、B、C、D、E节点，并替换F、G、H、I、J节点坐标。

图1 增减网格层时层动边界网格拓扑结构变化Fig.1 Topological structure change of the moving boundary mesh when the grid layer is changed

匹配的任意边界网格模型[6]，即商用软件中的r任意边界模型，是指由于相对运动等原因在同一位置有两个边界存在，再通过其上的边界面找到对应关系以传递流场中的信息，以按内部单元面处理。由于动态层网格在三维情况下只适用于规则的六面体或者三棱柱网格单元，且本文只对初始网格进行编号，所以要将整个计算域划分为气缸区域和燃烧室区域，在两个区域结合处采用任意边界面的处理方法。本文为处理方便在网格划分时将这两个边界上的网格一一对应(在Gambit软件中选中Link mesh face选项)，这样两个边界面上的所有网格几何信息全部相同，单元面间的相交判断就可以通过边界单元面的中心坐标是否相同来实现，进而确定边界单元面间的对应关系和所属单元。

1.2 流场控制方程

流体的运动要满足基本守恒定律，其控制方程组是各种CFD软件求解的基本数学模型，通用输运方程的积分形式可写为

(1)

式中:ρ是流体密度，U和Ub分别是流体和网格面运动速度，图1中红色运动边界处的边界单元面需要计算其Ub，本文采用上一时间段平均移动速度计算，n是单元面的外法线方向单位向量，V和A分别为网格单元体积和单元面的面积，φ为通用形式控制方程的求解变量，Γφ为对应的扩散系数，Sφ为相应的源项。各方程的具体形式如表1所示。

表1 扩散系数和源项

表中，μeff为有效粘性系数，求解层流问题时μeff=μl，求解湍流时μeff=μl+μt，H为流体总焓，λ为导热系数，cp为定压比热，Prt为湍流普朗特数，τ为流体所受切应力，k-ε系列湍流方程同样可以写为上述形式，本文算例没有采用湍流方程，这里不再赘述。

2 数值方法

2.1 控制方程离散

时间项的离散计算类似于一阶常微分方程，本文采用欧拉隐式格式对其进行计算，近似可得离散形式为

(2)

式中，Δt为计算时间步长、(ρVPφ)o为上一时刻值，F(φ,t)为通用输运方程中除时间项外的所有其他项当前时刻值之和，所以这种格式称作为隐式格式。

通过对对流项的离散，其可以近似等于各单元面上的对流通量之和：

(3)

式中：nf为当前单元所含有的单元面个数，fj为通过该单元第j个单元面的流量，当流速与单元面外方向n相同时为正，否则为负。单元面上变量值φj的取值取决于对流项采用的格式，常采用的一阶迎风格式：

(4)

则离散后的对流项最终可写为

(5)

扩散项的离散要稍微复杂一些，将变量梯度分为法向n和交叉向τ两部分，略去中间计算过程，最终可得：

(6)

式(6)中右端第一项为法向扩散部分，第二项为交叉向扩散部分，近似离散为

(7)

(8)

(9)

通过有限体积法离散流体通用控制方程，最终可得代数方程组目标形式为

(10)

式中：φP和φPj分别为当前单元和第j个相邻单元的变量值，ap和aj分别为代数方程中当前单元和相邻单元变量系数，bp为代数方程源项，本文通过代数多重网格和共轭梯度法求解代数方程组。

2.2 时间项与对流项特殊处理

为方便程序计算，要对时间项和对流项进行特殊处理。采用一阶迎风格式离散对流项，可进一步写成如下形式：

(11)

可知，当fj>0时，fj=fout；当fj<0时，fj=fin。同样经过上述处理，忽略源项，离散后的连续方程可写为

(12)

将其代入式(11)，对流项最终可得：

(13)

式中:fin只能等于零或-fj，所以可让fin=-max(0, fj)。

离散后的时间项可进一步整理得：

(14)

最后将时间项和对流项合并在一起，可得：

(15)

根据上述各项对代数方程组的贡献，经整理各系数可表示为

这样就会对每个方程形成代数方程组，为使方程求解稳定，本程序采用求解小量形式进行求解，式(10)可变为

(16)

式中：φ'=φ-φo。GTEA采用分离式变量求解方法，即逐一求解各个方程组，本文采用的共轭梯度法和代数多重网格法来求解代数方程组，进而求得各变量。

2.3 非结构化网格SIMPLEC算法

动量预测是SIMPLE系列算法[7]的第一步，其主要思想是根据初始给定的压力p*和ρ*，通过求解动量方程得出预测步的U*。将动量方程中的压力梯度项提出进行单独处理，可得离散后动量方程的代数方程形式为

(17)

则可通过上式求得预测步的各单元中心的速度。

压力修正的主要思想是利用动量预测得出的单元中心的速度，再用求解动量方程相似的形式来计算网格面上的速度，即动量加权或Rhie-Chow插值的方法；并采取压力、速度、密度修正的方法，用压力修正量表示速度和密度的修正量，最后修正单元面的流量，将其代入连续方程。SIMPLEC算法得到的压力修正量和速度修正量之间的关系为

(18)

省略中间计算,最后得出压力修正方程的形式如下

(19)

动量修正的主要目的是通过压力修正方程求出的压力修正量来修正速度和压力。虽然速度修正量与压力修正量之间的关系是通过动量方程求出，但通过压力修正方程之后，修正后的速度是严格满足连续方程的，修正后的压力也同样满足动量方程。程序中修正的公式为

(20)

(21)

而最后修正的密度是通过状态方程和已修正后的速度计算得出。

图2所示为GTEA软件的计算流程图。

图2 程序计算流程图Fig.2 Program calculation flow chart

3 算例验证

为验证本文内燃机缸内动网格模型及其流场计算的正确与否，先后设计了以下的动态层网格模型流场算例和匹配的任意边界网格模型流场算例，并将计算结果与商业软件的计算结果进行了对比。

3.1 动态层网格模型验证

1)测试算例

任意选取了一个气缸模型，活塞的半径为0.017m，冲程为0.05m，曲轴转速为1 000r/min，回转半径为0.024 75m，连杆长度为0.1m。计算的压力初值为100 000Pa，温度初值为300K。采用六面体网格单元，共划分了550个网格，初始网格模型如图3所示，计算时间为180°BTDC到180°ATDC，计算得到的平均压力和温度曲线如图4所示。

由图4可知，计算所得的缸内平均压力和平均温度曲线与商业软件PUMPLINX计算所得结果几乎完全吻合。

图3 测试算例网格模型Fig.3 Grid model of test case

图4 气缸平均压力、温度曲线对比Fig.4 Cylinder mean pressure and mean temperature curve comparison

图5 TBD620内燃机气缸模型Fig.5 Internal combustion engine cylinder model of TBD620

2)TBD620气缸算例

TBD620内燃机活塞的半径为0.085m，冲程为0.2m；曲轴转速为1 800r/min，回转半径为0.097 5m，连杆长度为0.35m。计算的压力初值为30 000Pa，温度初值为380K。采用六面体网格单元，共划分了25 752个网格。初始缸内计算域网格模型如图5所示，计算得到的平均压力和温度曲线如图6所示。

由图6可知，计算所得的缸内压力和温度曲线与商业软件FLUENT的计算所得结果几乎完全吻合。

3.2 燃烧室流场仿真验证

为了模拟完整的缸内过程流场仿真，还要加上燃烧室部分，建立的TBD620完整的缸内燃烧室模型如图7所示，气缸部分为六面体网格，燃烧室部分为金字塔和四面体混合单元，结合面处为一一对应的四边形面单元，网格数为10 786。

如图8所示为缸内压缩和膨胀过程所得到的压力和温度曲线与FLUENT计算结果对比图，可看出结果还是基本完全重合，说明了任意边界网格模型与流场耦合正确。

图6 气缸平均压力、温度曲线对比Fig.6 Cylinder mean pressure and mean temperature curve comparison of TBD620

图7 TBD620内燃机气缸Fig.7 Complete model of engine of TBD620

如图9所示为10°BTDC和10°ATDC时的气缸纵截面矢量图和全域速度矢量图的仰视图，可以发现活塞在上行到上止点之前，四周的气体流向燃烧室，会在燃烧室的两个凹坑内形成非常明显滚流漩涡，同理当活塞下行到上止点后，燃烧室内的气体由流出到活塞定于气缸盖的空间，会在燃烧室凹坑内形成两个逆向的滚流漩涡，可见这种气体流动形式可以使燃料与空气混合完全，符合柴油机工作时的实际情况有利于燃料燃烧。

图8 气缸平均压力、温度曲线对比Fig.8 Cylinder mean pressure and mean temperature curve comparison of TBD620

图9 流场速度矢量图Fig.9 Flow velocity vector diagram

4 结论

本文在动态层网格模型和任意边界网格模型的基础上，采用有限体积法离散流场控制方程，SIMPLEC算法求解NS方程，开发了一套三维缸内动网格流场数值模拟程序，即自主程序GTEA。并利用气缸算例和燃烧室算例进行验证，将GTEA仿真模拟所得的缸内平均压力值和缸内平均温度值与商业软件PUMPLINX和FLUENT计算结果对比，发现本文所得的缸内曲线与商业软件所得结果几乎完全吻合，并且给出了TBD620柴油机的缸内流场速度矢量图，就目前来讲，本自主程序可为相关仿真研究提供一定的指导和参考意义。

[1]舒歌群, 马维忍, 许世杰, 等. 喷雾夹角对柴油机性能影响的数值模拟[J]. 工程热物理学报, 2008, 29(7): 1239-1242. SHU Gequn, MA Weiren, XU Shijie, et al. Simulation of the effect of spray angle on diesel performance[J]. Journal of engineering thermophysics, 2008, 29(7): 1239-1242.

[2]AMSDEN A A, O'ROURKE P J, BUTLER T D. KIVA-II: a computer program for chemically reactive flows with sprays[R]. Technical Report LA-11560-MS. Los Alamos, New Mexico: Los Alamos National Laboratory, 1989.

[3]明平剑. 基于非结构化网格气液两相流数值方法及并行计算研究与软件开发[D]. 哈尔滨: 哈尔滨工程大学, 2008. MING Pingjian. Development of numerical modeling for gas-liquid two-phase flows based on unstructured grids and parallel computing[D]. Harbin: Harbin Engineering University, 2008.

[4]雷国东. 非结构网格FVM在复杂几何结构的湍流反应流计算中的应用研究[D]. 哈尔滨: 哈尔滨工程大学, 2008. LEI Guodong. The application and research of the unstructured grid FVM in the turbulent reaction flows simulation with complex geometries[D]. Harbin: Harbin Engineering University, 2008.

[5]张志荣, 冉景煜, 张力, 等. 内燃机缸内气体CFD瞬态分析中动态网格划分技术[J]. 重庆大学学报: 自然科学版, 2005, 28(11): 97-100, 121. ZHANG Zhirong, RAN Jingyu, ZHANG Li, et al. Techniques of dynamic mesh in the CFD transient analysis of an engine[J]. Journal of Chongqing university: natural science edition, 2005, 28(11): 97-100, 121.

[6]孙华文. 基于非结构动网格的内燃机瞬态流场仿真[D]. 哈尔滨: 哈尔滨工程大学, 2014. SUN Huawen. Numerical simulation of the transient flow in diesel engine based on unstructured dynamic meshes[D]. Harbin: Harbin Engineering University, 2014.

[7]PATANKAR S V, SPALDING D B. A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[J]. International journal of heat and mass transfer, 1972, 15(10): 1787-1806.

[8]刘永丰. 基于非结构网格的缸内两相反应流数值模拟方法研究及软件开发[D]. 哈尔滨: 哈尔滨工程大学, 2013. LIU Yongfeng. Development of numerical simulation for the process of two phases reacting flow in cylinder based on unstructured grids[D]. Harbin: Harbin Engineering University, 2013.

[9]RHIE C M, CHOW W L. Numerical study of the turbulent flow past an airfoil with trailing edge separation[J]. AIAA journal, 1983, 21(11): 1525-1532.

[10]王福军. 计算流体动力学分析-CFD软件原理与应用[M]. 北京: 清华大学出版社, 2004: 1-2, 25-28. WANG Fujun. Computational fluid dynamics analysis-theory and application of CFD software[M]. Beijing: Tsinghua University Press, 2004: 1-2, 25-28.

[11]PERAIRE J, VAHDATI M, MORGAN K, et al. Adaptive remeshing for compressible flow computations[J]. Journal of computational physics, 1987, 72(2): 449-466.

[12]BATINA J T. Unsteady Euler airfoil solutions using unstructured dynamic meshes[J]. AIAA journal, 1990, 28(8): 1381-1388.

Development of numerical simulation program for 3D moving grid and flow field in the cylinder of a diesel engine

LIU Xiaogang, MING Pingjian, ZHANG Gejian, ZHANG Wenping

(College of Power and Energy Engineering, Harbin Engineering University, Harbin 150001, China)

As mature CFD simulation software is currently developed overseas and domestic companies lack intellectual property rights, to ameliorate this problem, Harbin Engineering University developed a set of simulation programs that, after many years of accumulation, have proven to be remarkably successful. To continue conducting a numerical simulation of the three-dimensional working process in the cylinder of an internal combustion engine, a dynamic layer model and a matching model with an arbitrary boundary were firstly established. Moreover, the finite volume method was used to discretize the flow field control equations, and the Navier-Stokes equation was solved using the SIMPLEC algorithm. Based on the FORTRAN language, a set of simulation software of the internal flow field inside the cylinder of a diesel engine was then developed. The air compression and expansion process inside the cylinder were subsequently simulated and analyzed. Moreover, results showed that the cylinder pressure curve and cylinder temperature curve agree well with results of commercial software PUMPLINX and FLUENT, testifying that the moving grid model and the flow field simulation method are correct and that the numerical simulation program has potential as a market application.

I.C. Engine; moving grid; CFD; numerical simulation; program development; software simulation

2015-09-08.

日期：2016-10-12.

国家自然科学基金项目(51206031,51479038);国防预研基金项目(034315003).

刘晓刚(1990-)，男，硕士； 明平剑(1980-)，男，副教授，博士生导师.

明平剑, E-mail: pingjianming@hrbeu.edu.cn.

10.11990/jheu.201509023

TK402

A

1006-7043(2016) 11-1485-07

刘晓刚，明平剑，张格健，等. 柴油机缸内三维动网格和流场数值模拟程序开发[J]. 哈尔滨工程大学学报, 2016, 37(11): 1485-1491. LIU Xiaogang, MING Pingjian, ZHANG Gejian, et al. Development of numerical simulation program for 3D moving grid and flow field in the cylinder of a diesel engine[J]. Journal of Harbin Engineering University, 2016, 37(11): 1485-1491.

网络出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20161012.0927.004.html