张琨， 陈昊天， 邓康耀， 胡志龙， 钱跃华， 刘博

(1.上海交通大学 动力机械及工程教育部重点实验室， 上海 200240； 2.上海船舶设备研究所， 上海 200031； 3.中国船舶集团有限公司第七一一研究所， 上海 201108； 4.中船动力研究院有限公司， 上海 200120)

近年来，随着内燃机强化程度的提高，进排气的流速更高和能量更大，特别是排气压力波变化梯度加大，造成压力波在进排气管中的传播过程更加复杂[1]。为了更加准确地预测进排气系统的能量损失和转化过程，进排气系统仿真需有更高精度的计算模型[2]。但是，高精度的计算格式通常需要牺牲计算效率，因此内燃机性能仿真在不降低计算效率的基础上提高进排气流动模型精度至关重要。

Benson等[3-4]开展了内燃机工作过程仿真，引入特征线方法为进排气管道一维非定常流动计算方法，该方法较简便，经改进后能达到较高的精度，但流量误差较大，不能计算带有EGR多种混合成分的流动计算，因此在随后的应用中逐渐被淘汰[5]。随着对称或迎风激波捕捉数值方法用于求解双曲守恒方程组的兴起，其逐渐取代特征线法成为一维气体动力学计算的主流格式[6]，并逐步发展出了Lax-Friedrichs格式、Lax-Wendroff格式、MacCormack格式、ROE格式、TVD格式、ENO格式[7]等多种高精度计算格式。上述迭代格式在内燃机进排气流动仿真过程中，通常网格都是均匀网格，提高计算精度需要增加网格密度或者采用高精度的求解格式，这将导致计算效率降低。同时，实际内燃机进排气过程中，不同区域不同时刻对于网格分辨率的需求是不一样的，例如在排气压力波峰面经过的区域，需要采用较高分辨率捕捉压力从而提高仿真精度，而在其他区域可以采用较粗的网格。因此，传统迭代均匀网格的使用难以实现计算精度与计算效率的平衡。

为了克服传统内燃机进排气一维非定常流动不能实时调整网格密度的缺点，本文在均匀网格TVD格式基础上，考虑空间位置对通量限制器中的Sweby因子的影响，推导建立适用于内燃机进排气流动计算的非均匀网格TVD迭代格式。而后通过捕捉管道中参数梯度变化较大的接触面位置，对该区域实施网格加密，从而实现网格自适应加密和合并，建立可变网格一维非定常流动数值仿真计算方法。最后通过建立激波管和内燃机整机算例对比变网格与传统计算方法的计算精度和计算效率。

1 非均匀网格一维非定常流动模型

1.1 基本控制方程及网格划分

考虑截面变化、有摩擦、传热的一维非定常流动控制方程为[1]：

Wt+F(W)x=S(x,W)

(1)

其中：

一维非定常流动方程组的Jacobi矩阵为[7]：

(2)

Jacobi矩阵J(W)由特征值组成的对角矩阵为Λ，右特征向量矩阵为P，左特征向量矩阵为Q。

图1 有限体积法非均匀网格划分示意Fig.1 Schematic diagram of finite volume method for heterogeneous meshing

1.2 源项处理

对于带有源项的一维非定常流动控制方程组式(1)，在从n时刻推进至n+1时刻时，可将其分裂为2步求解[8]：

Wt+F(W)x=0

(3)

(4)

为了得到二阶精度的解，采用二阶Runge-Kutta方法：

K1=ΔtS(tn,Wn)

(5)

K2=ΔtS(tn+Δt,Wn+K1)

(6)

(7)

1.3 非均匀网格TVD格式

对于不含源项的一维非定常流动控制方程组(3)，因传统的TVD格式是针对标量方程的，因此首先需要对控制方程组运用特征理论解耦。

定义新的特征变量:

U=QW

(8)

将方程组(3)左侧乘以Q，可以得到:

(QW)t+QJP(QW)x=0

(9)

结合式(8)得:

Ut+ΛUx=0

(10)

定义新的变量:

Fu=QF

(11)

根据Harten推导的证明，对于迎风迭代格式:

(12)

(13)

式中：h(φ(rj+1/2))为对角矩阵；r为Sweby因子，其含义为网格边界的上游梯度和下游梯度的比值。式(12)和(13)为TVD格式。在满足稳定性准则前提下，如果通量限制器φ(r)满足:

(14)

对于对角矩阵的第k行元素，均匀网格的Sweby因子定义[9]为:

(15)

图2为均匀网格与非均匀网格的通量求解示意图，网格间距离对参数的梯度变化率计算有较大影响，非均匀网格的网格间距离不同，因此非均匀网格和均匀网格的梯度变化率求解方法是不一样的。

图2 不同网格划分通量求解示意Fig.2 Diagram of flux solution for different mesh division

对于非均匀网格，需考虑空间位置对Sweby因子的影响，因此，非均匀网格的Sweby因子定义为:

(16)

按照上述定义，原TVD格式定义中的式(13)也需要相应修正为:

(17)

当采用均匀网格时，非均匀网格计算式(16)和式(17)退化为均匀网格计算式(15)和(13)。

将式(12)和(17)左乘矩阵P可以获得原始变量的迭代方程，即：

(18)

(19)

式(18)和(19)为考虑管道面积变化的一维非定常流动TVD格式，公式中单元网格内的参数用平均值表示，下标为j+1/2或j-1/2的交界面参数根据单元网格中的参数通过Roe平均获得。

2 管道变网格数值仿真方法

一维非定常流动数值仿真过程中为了保证计算稳定性，迭代步长需要满足CFL准则：

(20)

从式(20)可以看出，管道流动的迭代步长是由网格尺度Δxj、网格速度|uj|和音速aj决定的。分母|uj|+aj是根据流动状态迭代计算出的，其数值大小不可控，而分子网格尺度可以调整，通过调整分子网格尺度Δxj的大小，可以实现迭代步长的调整。一般情况下，管道数值仿真的网格尺度是定值，设置较大的网格尺度Δxj可以增加迭代步长，减少计算所需时间。但是，当沿管长方向管道参数变化较剧烈时，仿真结果不能准确反映各网格点的参数变化情况，仿真结果存在失真现象。如果设置较小的网格尺度，则会增加数据存储和计算时间，大大降低计算效率。通过接触断面区域采用局部加密网格，可以有效提高仿真精度和效率。但在内燃机实际运行过程中，由于接触断面随压力波传播沿管道运动，需要实时捕捉接触断面的位置，并对该区域实施局部加密。因此，为了满足计算精度和计算效率的平衡，在非均匀网格TVD格式基础上，设计管道变网格计算方法。网格划分初始为较粗的网格，计算过程中根据各网格间物理量波动情况，自适应实现局部网格尺度的加密或合并。

2.1 基本原理

为了表征管道内参数变化，定义光滑度指数为：

(21)

式中：φ为密度、速度、浓度等物理量，一维非定常流动基本方程组中，每个方程中均包含密度参数，可以选取密度梯度为光滑度指标参数；Δxj-1,j表示网格j-1与网格j中点之间的距离，下标ref表示该参数为参考值。光滑度指数σ大表示该区域参数变化剧烈，需要提高网格分辨率，即进行网格加密；当光滑度指数σ较小时，参数变化平缓，已经加密的网格可以进行合并操作，减少数据存储空间及计算时间，提高计算效率。

1)当σj≥ε时，当前网格需进行加密操作，将当前网格沿中点一分为二，变为2个网格，由于计算采用的是有限体积法，在粗网格中添加细网格后，能够保证网格交界面上的流量守恒；

2)当σj<ε时，当前网格需进行合并操作，将当前加密网格与其邻近网格合并为一个网格。

有时，对于粗网格进行一次加密操作后仍不能满足计算精度的要求，需要进行多级嵌套加密，加密或合并的阈值ε需要根据加密级数做相应变化，即对于第k级加密：

εk=ε/2k

(22)

综上所述，变网格加密计算基本步骤为：

1)参数输入，给出网格尺度和初值、阈值；

2)初始化，划分网格并给各网格参数赋初值；

3)根据当前时刻各网格参数及阈值，确定该网格是否需要加密或合并，按照前面制定的网格加密和合并实施方案，计算各网格参数；

系统库文件提供用户程序，可直接使用API接口实现使用注册、控制密钥、加解密功能，内部实现与内核驱动之间的数据交互，支持设备复用功能[17]。

4)采用非均匀网格TVD求解格式，迭代下一时刻各网格点参数；

5)判断计算是否完成，如果完成则退出计算，如果没完成则进入第3)步继续计算。

2.2 数据结构设计

从上面的叙述不难看出，本文采用的网格加密生成方法为二分法，针对这种类型的数据生成方法，选用二叉树的树形存储结构用于计算过程的数据生成、存储计算和数据删除。如图3所示，对于父网格j，实施加密后，在第二级分裂为2个子网格，网格编号为2j-1和2j。为了保证数据结构的完整性和易操作性，在进行合并操作时，只能将同一父网格下的2个子网格合并，即将第2层的2j-1和2j节点合并为第1层的节点j。

2.3 网格数据转换

图3 变网格二叉树网格划分及数据存储结构示意Fig.3 Diagram of variable mesh binary tree mesh division and data storage structure

网格合并过程中，父网格的参数根据2个子网格的平均值求得，即:

(23)

与网格合并过程类似，网格加密过程中，子网格的参数根据父网格及其邻近网格的参数插值求得，加密过程子网格插值公式为:

(24)

(25)

根据网格加密插值式(24)和(25)，反推网格合并过程，可以得到与网格合并相同的式(23)，保证了网格加密和合并过程的一致性。

3 变网格数值仿真应用

3.1 激波管算例

一维非定常流动激波管问题具有精确解，可以仿真激波管中的气体流动验证非均匀网格TVD格式的准确性。建立如图4所示激波管算例，管道长1 m，中间由膜片分割，膜片两侧网格赋初值如图所示，管道内为理想气体，绝热系数为1.4，气体初始流速均为0 m/s。

图4 管道激波管算例示意Fig.4 Diagram of calculation case of pipeline shock wave tube

3.1.1 均匀网格不同网格密度对比

上述算例中，对比设置不同网格数(50、100、200)时，管道内激波发生后t=0.5 ms时压力、密度、速度和温度分布情况，如图5所示。从图5中可以看出，在不同网格数情况下都能拟合出沿管道方向物理参数的变化，参数变化平缓的区段，其数值仿真结果较好，激波面仿真拟合效果较差。网格数越多，管道计算结果分辨率越高，模拟结果与精确解的拟合度越高。随着网格密度的增加，仿真激波发生后t=0.5 ms所需的CPU时间分别为0.025、0.036、0.064 s，因此，随着网格数增加，CPU计算所需时间增加，并且呈线性增长趋势。

图5 不同网格数管道激波管算例t=0.5 ms时参数分布Fig.5 The parameter distribution of shock wave tube with different mesh numbers under TVD format (t=0.5 ms)

3.1.2 非均匀网格与均匀网格对比

通过分析上述激波发生后t=0.5 ms时的算例参数分布，可以看出，该时刻管道中流动状态发生较大变化的区域为以下3个区域：

A区域：0.1～0.3 m，该区域速度状态由0至250 m/s连续变化，为左行波到达区域；

B区域：0.5～0.7 m，该区域密度温度发生急剧变化，为发生接触间断的区域；

C区域：0.7～0.8 m，该区域速度状态由250～0 m/s变化，为右行波到达区域。

存在物理参数急剧变化的区域是仿真计算产生误差较大的区域，可以通过对局部区域网格加密以达到提高计算精度的目的。因此，针对激波管算例，将上述区域网格加密，管道非均匀网格的网格尺度分布示意如图6所示。在A区域(0.1～0.3 m)，网格尺度ΔxA为0.01 m，与网格数为100的均匀网格网格密度相同；在B区域和C区域(0.5～0.8 m)，网格尺度ΔxB,C为0.005 m，与网格数为200的均匀网格网格密度相同；其他区域网格尺度为0.02 m，与网格数为50的均匀网格网格密度相同。采用非均匀网格激波管算例进行划分的网格总数为110个。

激波管算例中，计算上述非均匀网格管道内工质激波发生后t=0.5 ms时压力、密度、速度和温度分布情况，并将结果与均匀网格(网格数100和200)方案进行对比，压力和密度结果如图7所示。在左行波区域A，由于非均匀网格110方案的网格密度与均匀网格100方案相同，所以2个方案在该区域各状态参数分辨率相同，模拟结果仿真精度相当，但仿真精度都略小于均匀网格200的算例；在接触间断区域B和右行波区域C，由于非均匀网格110方案的网格密度与均匀网格200方案相同，所以两方案在该区域各状态参数分辨率相同，模拟结果仿真精度相当。结合表1所示的非均匀网格和均匀网格激波管算例误差可以看出，均匀网格100的仿真精度最低，最大误差为4.13%；采用与均匀网格100网格数量相当的非均匀网格方案，不同参数的预测精度都有一定提高，非均匀网格的误差均在3.1%以内。综上所述，采用非均匀网格方案后，在网格数增加不多的前提下，整体仿真精度较网格数同量极的均匀网格方案有大幅提高，局部高密度区域可以达到与2倍均匀网格密度相当的精度。

图6 管道非均匀网格划分方案示意Fig.6 Diagram of non-uniform meshing scheme for pipeline

均匀网格100、均匀网格200、非均匀网格110等3种方案在仿真激波发生0.5 ms过程消耗CPU时间分别为0.036、0.064、0.040 s。非均匀网格110与均匀网格100方案的网格数相当，仿真耗时较为接近。均匀网格200方案由于网格数较多，计算量较大，仿真耗时是非均匀网格110方案的1.6倍。

图7 非均匀和均匀网格激波管算例t=0.5 ms时参数分布对比Fig.7 Comparison of parameter distribution between the non-uniform mesh and the uniform mesh solution format shock wave tube cases (t=0.5 ms)

表1 非均匀和均匀网格激波管算例误差

综上所述，采用非均匀网格，与同等均匀网格密度方案相比，能够实现在消耗相同CPU时间基础上，提高仿真精度1.03个百分点；与加密均匀网格方案相比，能够在满足仿真精度相当的基础上，降低计算耗时37.5%。

3.1.3 均匀网格与变网格对比

搭建相同的管道激波管算例，变网格格式初始网格数设置为50，阈值设置为5，加密最高级数设置为3，即最多能加密至200个网格。对比均匀网格(网格数为50和200)和变网格(初始网格数50)等不同数值仿真方法的仿真精度和仿真效率。

变网格格式与均匀网格格式数值仿真结果如图8所示。从图中可以看出，与均匀网格数为50求解格式相比，变网格格式仿真精度有明显提高，尤其是在参数数值变化较为剧烈的激波接触面，变网格格式更加逼近精确解。结合表2变网格与均匀网格激波管算例仿真误差可以看出，变网格求解格式对于压力、密度的数值仿真精度能够达到与均匀网格200相当的水平，各参数误差均在2.2%以内。

图8 变网格与均匀网格求解格式管道激波管算例t=0.5 ms时参数分布对比Fig.8 Comparison of parameter distribution between variable grid and uniform grid solution format shock wave tube cases (t=0.5 ms)

表2 不同网格求解格式管道激波管算例误差

图9所示为使用变网格格式求解管道激波管算例，在t=0.5 ms时各加密级网格分布图。沿管道长度方向生成了3类网格密度级别不同的区域，其中第1级网格的网格密度与初始网格相同，网格离散长度为0.02 m；第2级网格经历了一次网格加密过程，网格离散长度为0.01 m；第3级网格的网格密度最密集，相应的网格离散长度也最小，其离散长度为0.005 m。结合图8(b)可以发现，在密度变化梯度较大的3个区域，均生成了网格最密集的第3级网格，而在密度参数变化平缓的区域，仍然保留为初始的第1级网格，第1级和第3级网格间有少量第2级网格过渡。在初始网格数为50的基础上，计算过程中实施了3级网格加密，仿真结束时，沿管道长度方向有第1级网格17个，第2级网格22个，第3级网格88个，3级网格总计有127个，与仿真精度相当的均匀网格数为200的算例相比网格数减少了73个，相应的数据存储量也降低了36.5%。

图9 变网格求解格式管道激波管算例t=0.5 ms时各加密级网格分布Fig.9 Mesh distribution of each encryption layer in the variable-scale solution format pipeline shock wave tube cases

均匀网格50、均匀网格200、变网格50计算激波管算例所需的CPU时间分别为0.025、0.064、0.039 s。变网格求解格式所消耗的CPU时间比均匀网格数50略高。与均匀网格数200的算例相比，变网格求解格式所消耗的CPU时间仅是均匀网格数200的60.94%。采用变网格求解格式能够在保持数值仿真精度不降低的前提下，提高计算效率。

3.2 内燃机排气管内的流动

为实现内燃机整机性能仿真，除了进排气管道流动仿真模型外，还需其他组件和边界模型。参照相关文献，其他组件和边界模型选用已经广泛使用的模型，如表3所示。为便于实现与进排气流动模型的耦合，边界模型均基于特征线方法建立。

表3 其他组件模型Table 3 Other component models

搭建HC4132柴油机整机试验台架及测试系统，相关细节见文献[11]，选取HC4132柴油机外特性1 500、1 700、1 900和2 100 r/min等4个工况点，其中2 100 r/min为柴油机标定转速，1 500 r/min为最大扭矩点转速，测量不同转速工况排气压力波以及整机性能指标。搭建HC4132柴油机整机仿真模型，对比管道采用均匀网格和变网格不同转速下排气压力波、整机性能预测精度以及计算效率。均匀网格算例网格数设置为10；管道变网格算例初始网格数设置为5，阈值设置为5，加密最高级数设置为3。

3.2.1 排气压力波对比

HC4132柴油机排气系统采用的是2缸连接一根排气管，即第1缸与第4缸共用一根排气管，第2缸与第3缸共用一根排气管，试验过程中在第3缸排气管处安装传感器用于测量排气压力波变化。对比管道采用均匀网格和变网格不同转速下排气压力波仿真结果与试验结果如图10所示。第3缸排气管测得的压力波，其压力波最大峰值出现在第3缸排气阶段；由于第3缸与第2缸邻近，且两缸连接在同一根排气管上，因此第2缸排气过程对第3缸排气管处压力波影响也较大。不同转速工况，采用变网格计算方案，第3缸排气压力波波峰波谷处预测精度均有明显提高。

图10 不同转速第3缸排气压力波仿真与试验结果对比Fig.10 Comparison of calculated and measured exhaust pressure wave for the third cylinder under different speed

表4为不同转速均匀网格和变网格第3缸排气压力波循环平均误差，均匀网格和变网格循环平均误差低转速工况预测精度低于高转速工况。1 500 r/min转速工况计算误差最大，均匀网格和变网格格式分别为10.62%和5.47%，采用变网格方案排气压力波预测精度最大可提高5.15个百分点。

表4 不同转速不同网格格式排气压力波循环平均误差

3.2.2 整机性能对比

HC4132柴油机主要性能参数不同求解格式仿真结果与稳态试验结果误差如表5所示。从表中可以看出，不同求解格式仿真结果中功率、涡轮进口温度等参数预测较准确，各转速工况误差均在2.82%以内，燃油消耗、中冷前温度、涡轮进口压力等参数预测精度虽然偏低，但最大误差也在4.68%以内。不同管道求解格式整机性能预测精度相当，均匀网格最大误差为4.68%，变网格最大误差为4.14%，采用变网格整机性能预测精度较均匀网格提高0.54个百分点。

表5 HC4132不同求解格式整机性能仿真误差

3.2.3 计算效率对比

程序运行平台处理器主频为4.20 GHz，分别进行4个工况计算，设置仿真20个循环，管道采用不同求解格式进行计算消耗CPU时间如表6所示，表中差值为变网格比均匀网格计算增加的耗时。变网格方案虽然初始网格数较少，但计算过程中需要生成和合并网格造成计算耗时增加，不同转速工况采用变网格方案均略高于均匀网格方案。仿真计算20个循环，变网格方案平均增加耗时2.08 s。

表6 不同求解格式整机性能仿真消耗CPU时间

4 结论

1)激波管算例中，均匀网格TVD格式网格数越多，管道计算分辨率越高，模拟结果与精确解拟合度越高，但随着网格数增加，CPU消耗时间呈线性增长；采用非均匀网格，与同等网格密度方案相比，能够在消耗同等CPU时间基础上，提高仿真精度1.03个百分点，与加密网格方案相比，在满足仿真精度相当的基础上，可以降低计算耗时37.5%。

2)初始网格数为50的变网格计算方法与网格数为200的均匀网格计算方法相比，两者数值仿真精度相当，误差均在2.2%以内，但前者数据存储量比后者降低了36.5%，所消耗的CPU时间仅是后者的60.94%。因此，采用变网格数值仿真方法后，可以在保证数值仿真精度不降低的前提下，降低数据存储量，减少CPU消耗时间，提高计算效率。

3)整机性能仿真中，采用变网格计算20个循环平均增加CPU耗时2.08 s，排气压力波预测精度较均匀网格方案最大可提高5.15个百分点，整机性能预测精度较均匀网格最大可提高0.54个百分点。