发动机缸内工作过程计算

发动机工作过程数值计算作业

缸内封闭过程数值计算

学    院：汽车学院

专    业：动力机械及工程

*    *****

学    号：**********

*    *****

二〇一三年五月

气缸封闭过程的数值计算

发动机缸内工作过程的数值计算，是以内燃机缸内各工作阶段的物理模型为基础， 通过微分方程来对其各阶段工作过程进行数学描述， 然后通过程序编写求解微分方程， 得到缸内温度、压力等参数随曲轴转角的变化曲线。

一、基本热力学模型

图1  发动机缸内热力系统

在简化假设的基础上，取气缸为一个热力学系统，如图1所示。这个热力系统包括了质量交换项， 如排气 dmA，进气 dmE，喷入气缸内瞬时燃料质量 dmB ； 与能量交换项，如焓变， 功， 燃烧放热等。 图中 T 、P 、V、 m 及 u 分别为缸内瞬时气体温度、压力、体积、燃料质量及比内能。

二、简化假设

1.假定工质为理想混合气体；

2.假定缸内各处温度、压力及混合气浓度均匀；

3.用纯空气的气体常数代替混合气气体常数；

4.假定扫气完全，即不考虑残余废气；

5.不计漏气损失，并假定只有在燃烧始点才有燃油喷入气缸；

6.按代用燃烧规律进行喷油，并认为着火延迟等于零；

7.假定放热率为100%。

三、基本方程

1、气体状态方程

PV=mRT

2、压缩期()

能量方程：＝（－p）

质量方程：＝0

3、燃烧期()

能量方程：＝（＋－p－u－m）

质量方程：＝

4、膨胀期()

能量方程：＝（－p）

质量方程：＝0

四、其他计算公式

1、气缸瞬时工作容积：

=

式中：-气缸行程容积

      VL-发动机排量

-压缩比

-曲柄连杆比

-活塞行程

-气缸直径

2、气缸容积随曲轴转角变化率：

3、单位曲轴转角的传热量：

式中：-活塞

-气缸盖

-气缸套

-传热表面积

      ；；

-传热表面平均温度

；；

-换热面平均的瞬时传热系数

利用Woschni公式：

式中：-活塞行程容积 

-气体速度系数

   压缩、膨胀阶段  

   换气阶段        

-燃烧室形状系数，直喷式燃烧室=3.24X10-3

-活塞平均速度，计算得6.9m/s

、、分别为压缩始点的气缸压力、容积、温度。分别取值为102kPa, 0.L, 375K

-发动机倒拖时的气缸压力,取为100kPa

4、代用放热规律：

累积放热量：

瞬时放热率：

式中：-燃烧效率,取为100%

-燃料低热值，取为42500kJ/kg 

-燃烧品质指数,取为1.0

-燃烧持续角，取为80°CA

-循环喷油量,为3.5421X10-6kg

-燃烧始点的曲轴转角，为338°CA

5、工质内能：

u=u(,T ）

   式中：-瞬时过量空气系数

1）压缩阶段，取。

2)燃烧阶段：

3)膨胀阶段：

         =const

              式中：-理论空气量，取14.3kg/kg

    -喷入燃料量随曲轴转角的变化函数

                    -每循环喷油量

6、等容比热：

7、比内能对过量空气系数的偏导

  8、燃烧期混合气的质量：

mH=mL+mB

五、龙格—库塔法的计算步骤

        已知  ，

步骤一：

步骤二：

步骤三：

步骤四：

步骤五：

步骤六：

步骤七：

步骤八：

步骤九：

步骤十：

步骤十一： 

步骤十二：

步骤十三：

六、发动机基本参数

发动机形式：四冲程、直喷、水冷、单缸柴油机

   发动机型号：ZH105W

压缩比：＝16.5

排量：=0.96L

活塞行程：s＝115mm

缸径：D＝105mm

循环喷油量：3.5421X10-5kg

燃烧持续角：80°CA

曲柄半径和连杆长度比 ＝0.3

喷油提前角：

配气相位：

    进气提前（上止点前）12°CA

    进气迟后（下止点后）38°CA

    排气提前（下止点前）55°CA

    排气迟后（上止点后）12°CA

七、程序框图

八、源程序（matlab）

1、压缩过程

自定义ysgc函数并生成m文件，源程序如下：

function f=ysgc(fai,T)

P1=102000;%压缩始点气缸内气体压力，单位pa

V1=0.000;%压缩始点气缸内气体体积(计算得)，单位m3

T1=375;%压缩始点气缸内气体的温度（计算得），单位K

R=8.314;%气体常数，单位J/mol.k

ML=29;%空气摩尔质量，单位g/mol

nL=P1*V1/R/T1;%空气摩尔量，单位mol

mL=nL*ML;%压缩始点气缸内气体质量，单位g

lmd=10000;%过量空气系数

D=0.105;%气缸直径，单位m

c1=2.29;%气体速度系数

cm=6.9;%活塞平均速度（计算得），单位m/s

c2=3.24*10^-3;%燃烧室形状系数

Vh=0.96/16.5*15.5/1000;%行程容积，单位m3

P0=100*10^3;%发动机倒拖时的气缸压力,单位pa

ksai=16.5;%压缩比

lmds=0.3;%曲柄连杆比

V=Vh/2*(2/(ksai-1)+1-cos(fai)+1/lmds*(1-sqrt(1-lmds^2*sin(fai)*sin(fai))));%气缸瞬时工作容积，单位m3

A1=pi*D^2/4;%活塞传热表面积，单位m2

A2=1.5*A1;%气缸盖传热表面积，单位m2

A3=4*V/D;%气缸套传热表面积，单位m2

Tw1=560;%活塞表面的平均温度，单位K

Tw2=548;%气缸盖表面的平均温度，单位K

Tw3=548;%气缸套表面的平均温度，单位K

n=1800;%发动机的转速，单位r/min

P=nL*R*T/V;%气缸内气体瞬时压力，单位pa

ref=0.130*T^(-0.53)*(P/100000)^0.8*D^(-0.2)*(c1*cm+c2*Vh*T1*(P-P0)/P1/V1)^0.8;%换热面平均的瞬时传热系数，单位KJ/m2.s.K

dQw_dfai=1/6/n*ref*(A1*(Tw1-T)+A2*(Tw2-T)+A3*(Tw3-T));%单位曲轴转角的传热量，单位KJ

cv=0.14455*(-3*(0.0975+0.0485/lmd^0.75)*(T-273)^2*10^(-6)+2*(7.768+3.36/lmd^0.8)*(T-273)*10^(-4)+(4.6+46.4/lmd^0.98)*10^(-2));%定容比热，单位KJ/Kg.K

dV_dfai=Vh/2*(sin(fai)+lmds*sin(fai)*cos(fai)/(sqrt(1-lmds^2*sin(fai)*sin(fai))));%气缸容积随曲轴转角的变化率

f=1000/mL/cv*(dQw_dfai-P/1000* dV_dfai);

调用ysgc函数并用龙格库塔法法解微分方程，源程序如下：

[fai,T]=ode45('ysgc',[(180+38)/180*pi:5/180*pi:(360-22)/180*pi],375)

2、燃烧过程

自定义rsgc函数并生成m文件，源程序如下：

function f=rsgc(fai,T)

P1=102000;%压缩始点气缸内气体压力，单位pa

V1=0.000;%压缩始点气缸内气体体积，单位m3

T1=375;%压缩始点气缸内气体的温度，单位K

D=0.105;%气缸直径，单位m

R=8.314;%气体常数，单位J/mol.k

ML=29;%空气摩尔质量，单位g/mol

nL=P1*V1/R/T1;%空气摩尔量，单位mol

mL=nL*ML;%压缩始点气缸内气体质量，单位g

c1=2.29;%气体速度系数

cm=6.9;%活塞平均速度，单位m/s

c2=3.24*10^-3;%燃烧室形状系数

Vh=15.5*0.96/16.5/1000;%行程容积，单位m3

yinu=1;%燃烧效率

m=1;%燃烧品质指数

L0=14.3;%柴油理论空燃比

mB0=0.03542;%每循环喷油量，单位g

Hu=42500;%燃料低热值，单位KJ/Kg

faivb=(360-22)/180*pi;%喷油始点对应的曲轴转角

drfai=80/180*pi;%燃烧持续角

ksai=16.5;%压缩比

lmds=0.3;%曲柄连杆比

V=Vh/2*(2/(ksai-1)+1-cos(fai)+1/lmds*(1-sqrt(1-lmds^2*sin(fai)*sin(fai))));%气缸瞬时工作容积，单位m3

QB=yinu*mB0/1000*Hu*(1-exp(-6.908*((fai-faivb)/drfai)^(m+1)));%累积放热量，单位KJ

dQB_dfai=6.908*yinu*mB0/1000*Hu/drfai*(m+1)*((fai-faivb)/drfai)^(m+1)*exp(-6.908*((fai-faivb)/drfai)^(m+1));%单位曲轴转角对应的放热量

mB=QB/Hu;%曲轴转角对应fai的时刻,喷入气缸燃料的累积质量，单位Kg

lmd=mL/1000/(mB+0.0000000001)/L0;%过量空气系数

mH=mL/1000+mB;%曲轴转角对应fai的时刻,气缸内混合气的质量，单位Kg

cv=0.14455*(-3*(0.0975+0.0485/lmd^0.75)*(T-273)^2*10^(-6)+2*(7.768+3.36/lmd^0.8)*(T-273)*10^(-4)+(4.6+46.4/lmd^0.98)*10^(-2));%定容比热，单位KJ/kg.K

dV_dfai=Vh/2*(sin(fai)+lmds*sin(fai)*cos(fai)/(sqrt(1-lmds^2*sin(fai)*sin(fai))));%气缸容积随曲轴转角的变化率,单位m3

A1=pi*D^2/4;%活塞传热表面积，单位m2

A2=1.5*A1;%气缸盖传热表面积，单位m2

A3=4*V/D;%气缸套传热表面积，单位m2

Tw1=600;%活塞表面的平均温度，单位K

Tw2=548;%气缸盖表面的平均温度，单位K

Tw3=548;%气缸套表面的平均温度，单位K

n=1800;%发动机的转速，单位r/min

P0=100*10^3;%发动机倒拖时的气缸压力,单位pa

P=mH*1000/29*R*T/V;%气缸内气体瞬时压力，单位pa

ref=0.130*T^(-0.53)*(P/100000)^0.8*D^(-0.2)*(c1*cm+c2*Vh*T1*(P-P0)/P1/V1)^0.8;%换热面平均的瞬时传热系数，单位KJ/m2.s.K

dQw_dfai=1/6/n*ref*(A1*(Tw1-T)+A2*(Tw2-T)+A3*(Tw3-T));%单位曲轴转角的传热量，单位KJ

dmH_dfai=1/Hu*dQB_dfai;%气缸内气体随曲轴转角的变化率,单位Kg

u=0.14455*(-1*(0.0975+0.0485/lmd^0.75)*(T-273)^3*10^(-6)+(7.768+3.36/lmd^0.8)*(T-273)^2*10^(-4)+(4.6+46.4/lmd^0.93)*(T-273)*10^(-2)+1356.8);%比内能，单位KJ/Kg

zdu_zdlmd=0.14455*(0.75*0.0485*lmd^(-1.75)*(T-273)^3*10^(-6)-0.8*3.36*lmd^(-1.8)*(T-273)^2*10^(-4)-0.93*46.4*lmd^(-1.93)*(T-273)*10^(-2));

dlmd_dfai=-mL/1000/L0/(mB+0.0000000001)^2/Hu*dQB_dfai;

f=1/mH/cv*(dQB_dfai+dQw_dfai-P/1000*dV_dfai-u*dmH_dfai-mH*zdu_zdlmd*dlmd_dfai);

调用rsgc函数并用龙格库塔法解微分方程，源程序如下：

[fai,T]=ode45('rsgc',[(360-22)/180*pi:1/180*pi:(360-22+80)/180*pi],843)

3、膨胀过程

自定义pzgc函数并生成m文件，源程序如下：

function f=pzgc(fai,T)

P1=102000;%压缩始点气缸内气体压力，单位pa

V1=0.000;%压缩始点气缸内气体体积，单位m3

T1=375;%压缩始点气缸内气体的温度，单位K

R=8.314;%气体常数，单位J/mol.k

ML=29;%空气摩尔质量，单位g/mol

nL=P1*V1/R/T1;%空气摩尔量，单位mol

mL=nL*ML;%压缩始点气缸内气体质量，单位g

mB0=3.5421*10^(-5);%每循环喷油量，单位Kg

L0=14.3;%理论空燃比

lmd=mL/1000/L0/mB0;%过量空气系数

D=0.105;%气缸直径，单位m

c1=2.29;%气体速度系数

cm=6.9;%活塞平均速度，单位m/s

c2=3.24*10^-3;%燃烧室形状系数

S=0.115;%活塞行程，单位m

Vh=pi*D^2*S/4;%行程容积，单位m3

P0=100*10^3;%发动机倒拖时的气缸压力,单位pa

ksai=16.5;%压缩比

lmds=0.3;%曲柄连杆比

V=Vh/2*(2/(ksai-1)+1-cos(fai)+1/lmds*(1-sqrt(1-lmds^2*sin(fai)*sin(fai))));%气缸瞬时工作容积，单位m3

A1=pi*D^2/4;%活塞传热表面积，单位m2

A2=1.5*A1;%气缸盖传热表面积，单位m2

A3=4*V/D;%气缸套传热表面积，单位m2

Tw1=560;%活塞表面的平均温度，单位K

Tw2=548;%气缸盖表面的平均温度，单位K

Tw3=548;%气缸套表面的平均温度，单位K

n=1800;%发动机的转速，单位r/min

P=(mL+1000*mB0)/29*R*T/V;%气缸内气体瞬时压力，单位pa

ref=0.130*T^(-0.53)*（P/1000)^0.8*D^(-0.2)*(c1*cm+c2*Vh*T1*(P-P0)/P1/V1)^0.8;%换热面平均的瞬时传热系数，单位KJ/m2.s.K

dQw_dfai=1/6/n*ref*(A1*(Tw1-T)+A2*(Tw2-T)+A3*(Tw3-T));%单位曲轴转角的传热量

cv=0.14455*(-3*(0.0975+0.0485/lmd^0.75)*(T-273)^2*10^(-6)+2*(7.768+3.36/lmd^0.8)*(T-273)*10^(-4)+(4.6+46.4/lmd^0.98)*10^(-2));%定容比热

dV_dfai=Vh/2*(sin(fai)+lmds*sin(fai)*cos(fai)/(sqrt(1-lmds^2*sin(fai)*sin(fai))));%气缸容积随曲轴转角的变化率

f=1/(mL/1000+mB0)/cv*(dQw_dfai-P* dV_dfai);

调用pzgc函数并用龙格库塔法解微分方程，源程序如下：

[fai,T]=ode45('pzgc',[(360-22+80)/180*pi:5/180*pi:(540-55)/180*pi],997)

九、程序运行结果

1、压缩过程

表1  压缩过程运行数据结果

2、燃烧过程

表2  燃烧过程运行数据结果

3、膨胀过程

表3  膨胀过程运行数据结果

十、用运行数据结果绘制图像

绘制图像源程序：

在matlab命令窗口中输入fai、V、T、P的值；

利用函数plot（fai，V）、plot（fai，T）、plot（fai，P）绘制对应的曲线

气缸瞬时容积随曲轴转角变化曲线见图2、缸内温度随曲轴转角变化曲线见图3、缸内气体压力随曲轴转角变化曲线见图4

图2  气缸瞬时容积随曲轴转角变化曲线

图3  缸内温度随曲轴转角变化曲线

图4  缸内气体压力随曲轴转角变化曲线

十一、小结

本次模拟计算得到的缸内参数随曲轴转角的变化曲线与实际柴油机缸内工作情况基本上相符合，计算结果较为合理。但由于时间和经验有限，计算过程中用于模拟计算的某些参数的取值都是参照前人模拟计算过程中所用参数直接取值的，可能会不够精确，由此可能会导致模拟计算结果与真实情况不能完全贴合。