> restart:#"m16_p31"

A thermodynamic simulation of a reciprocating engine cycle is to be performed using a Weibe-function to model the combustion process, with typical parameters for a SI-ICE: a=5, m=3,  bs=-10º,  bs=-70º. The model is to be applied to a Yamaha YZ250FN four-stroke 250 cc motorcycle, single cylinder, spark-ignition engine, with the following dimensions: bore D=0.0770 m, stroke L=0.0526 m, with a compression ratio r=12.5. Only the compression and expansion strokes are to be considered, with the crank angle,  , as the independent variable, and only at the steady-state corresponding to full load at 8500 rpm, i.e., assume the cylinder initially filled with an stoichiometric octane/air mixture at ambient conditions (unit volumetric efficiency). Perform:

a) Compute the displacement volume, clearance volume, trapped mass of air and fuel, and chemical energy released per cycle.

b) Compute gas volume and energy release as a function of   and make a plot.
c) Establish the energy balance for the gas, and get a single differential equation to compute the pressure profile p( ).
d) Solve numerically for p( ), T( ) and dWsh/d  (the contribution to the shaft work at each step). Plot them.
e)      Compute the shaft power and the energy efficiency of the engine.

Data:

 > su1:="Aire":fuel:=C8H18;dat:=[D=0.0770*m_,L=0.0526*m_,Vd=249e-6*m_^3,r=12.5,N=(8500/60)/s_,a=5,m=3,theta_bs=-10*Pi/180,theta_bd=60*Pi/180,b=3]:evalf(%,3);

Sketch and nomenclature:

Fig. 1. Crank angle diagram showing enclosed volume in m3 (red line), its angle derivative (gree line), and valve controls and burning timing (in black).

D=bore

L=stroke

Vcl=clearance volume

Vd=displacement volume

b=Lrod/(L/2) (of order 2)

ivc=intake valve closes

evo=exhaust valve opens

bs=burning startup

bd=burning duration

Eqs. const.:

Slider and Crank-Shaft mechanism (a typical connecting rod / crank shaft length ratio =3 has been used):

 > eqCS:=V=Vcl+Vd*(sqrt(b^2-sin(theta)^2)+1-cos(theta)-b)/2;eqCS_diff:='V_theta=diff(rhs(eqCS),theta)';eqWeibe:=Qacc=Qcycle*(1-exp(-a*((theta-theta[bs])/theta[bd])^m));eqWeibe_diff:=Qacc_theta='diff(rhs(eqWeibe),theta)';eqWeibe_diff:=Qacc_theta=diff(subs(theta[bs]=theta_bs,theta[bd]=theta_bd,rhs(eqWeibe)),theta):

a)Compute the displacement volume, clearance volume, trapped mass of air and fuel, and chemical energy released per cycle.

b) Compute gas volume and energy release as a function of   and make a plot

 > Qacc_otto:=piecewise(theta<0,0,Qcycle);Qacc_:=piecewise(thetatheta_bs+theta_bd,Qcycle,rhs(eqWeibe));Qacc_theta_:=diff(Qacc_,theta):Qacc_theta_:=piecewise(thetatheta_bs+theta_bd,0,rhs(eqWeibe_diff)):uc:=180/Pi;plot(subs(Qcycle=Qcycle_,theta[bs]=theta_bs,theta[bd]=theta_bd,Vcl=Vcl_,dat,SI0,theta=theta_deg/uc,{Qacc_otto,Qacc_,Qacc_theta_}),theta_deg=-Pi*uc..Pi*uc,Qdot_and_Qacc=0..2000,color=black);

c) Establish the energy balance for the gas, and get a single differential equation to compute the pressure profile p( ).

The real combustion process is assimilated to a heat transfer process (i.e. DUtherm+DUchem=W+0 to DUtherm=W+Q).

d) Solve numerically for p( ), T( ) and dWsh/d  (the contribution to the shaft work at each step). Plot them.

Direct plot

 > dsol1 := dsolve({eqDp_,p(-Pi)=1e5}, numeric, p(theta),range=-Pi..Pi);with(plots):odeplot(dsol1,theta=-Pi..Pi);dsol1(0),p(0)=op(2,op(2,dsol1(0))):evalf(%);

But, before detailed plots are presented, the ideal Otto cycle is worked out for comparison with Weibe model.

Ni=Number of intervals for discretization of crank-angle. For debugging: #print(evalf([theta_(i),p_(i),V_(i),T_(i),Wu_(i)]))

 > Ni:=30:Dtheta:='2*Pi/Ni';for i from 1 to Ni do theta_(i):=-Pi+2*Pi*i/Ni;p_(i):=op(2,op(2,eval(subs(theta=theta_(i),dsol1(theta)))))*Pa_; V_(i):=subs(Vcl=Vcl_,dat,theta=theta_(i),rhs(eqCS)); T_(i):=subs(Adat,dat,p_(i)*V_(i)/(ma_*R));Wu_(i):=subs(dat,(p_(i)-p0)*subs(eqCS_diff,theta=theta_(i),V_theta));od: plot(subs(SI0,{potto,[seq([theta_(i),p_(i)],i=1..Ni)]}),theta=-Pi..Pi);plot(subs(SI0,{Totto,[seq([theta_(i),T_(i)],i=1..Ni)]}),theta=-Pi..Pi);plot(subs(SI0,{Wu_otto,[seq([theta_(i),Wu_(i)],i=1..Ni)]}),theta=-Pi..Pi);plot(subs(eqCS,Vcl=Vcl_,dat,SI0,{[V,potto,theta=-Pi..Pi],[seq([V_(i),p_(i)],i=1..Ni)]}),0..subs(dat,SI0,Vcl_+Vd));loglogplot(subs(eqCS,Vcl=Vcl_,dat,SI0,{[Totto,potto,theta=-Pi..Pi],[seq([T_(i),p_(i)],i=1..Ni)]}),'T_K'=300..5000,'p_Pa'=1e5..2.5e7);

HINT: if no fuel is added (put 0*mF), the plots would correspond to isentropic compression and expansion, with maximum values of p(0)=p0*r^gamma and T(0)=T0*r^(gamma-1).

e) Compute the shaft power and the energy efficiency of the engine.