Code covered by the BSD License  

Highlights from
Spark Ignition Engine GUI

image thumbnail

Spark Ignition Engine GUI

by

 

01 Apr 2013 (Updated )

computer program for thermodynamic engine model for a spark ignition engine.

eng_sim.m
%% Engine Data
Rc = str2num(get(handles.cr_stroke_ratio, 'String')); %2;                 % Connecting rod/stroke length ratio
%rc = str2num(get(handles.cmp_ratio, 'String')); %11;                % Compression Ratio
Vm = 0.00001 * str2num(get(handles.max_cyl_vol, 'String')); %0.00055;           % Maximum Cylinder Volume (m^3)
%ths = (pi/180) * str2num(get(handles.ign_onset, 'String')); %0.88 * pi;        % Ignition Onset (158.4 Deg)
thb = (pi/180) * str2num(get(handles.brn_dur, 'String')); %0.33 * pi;        % Burn Duration (59.4 Deg)
%rpm = str2num(get(handles.eng_rpm, 'String')); %6000;             % Engine
%(rpm)
ws = 2 * pi * rpm/60;   % Engine (rad/s)
% valve Data, exhaust
thevo = (pi/180) * str2num(get(handles.exh_vlv_open, 'String')); %2 * pi;         % Exhaust valve opens (rad)
thevc = (pi/180) * str2num(get(handles.exh_vlv_clsd, 'String')); %3 * pi;         % Exhaust valve closed (rad)
eportd = str2num(get(handles.exh_vlv_port_dia, 'String')); %0.04;          % Exhaust port diameter (m)
estemd = str2num(get(handles.exh_vlv_stem_dia, 'String')); %0.015;         % Exhaust valve stem diameter (m)
elm = str2num(get(handles.exh_vlv_max_lift, 'String')); %0.035;            % Exhaust valve max lift (m)
%%%%%%%%%%%%%%%%%
% Valve Data Intake
thivo = (pi/180) * str2num(get(handles.inlt_vlv_open, 'String')); %3 * pi;         % Inlet valve opens (rad)
thivc = (pi/180) * str2num(get(handles.inlt_vlv_clsd, 'String')); %4 * pi;         % Inlet valve closed (rad)
iportd = str2num(get(handles.inlt_vlv_port_dia, 'String')); %0.04;          % Inlet port diameter (m)
istemd = str2num(get(handles.inlt_vlv_stem_dia, 'String')); %0.015;         % Inlet valve stem diameter (m)
ilm = str2num(get(handles.inlt_vlv_max_lift, 'String')); %0.035;            % Inlet valve max lift (m)
%%%%%%%%%%%%%%%%%
% Intake Exhaust state
Pi = str2num(get(handles.intk_press, 'String')); %100;               % Intake Pressure (kPa)
Pe = str2num(get(handles.exhst_press, 'String')); %150;               % Exhaust Pressure (kPa)
Ti = str2num(get(handles.inlet_temp, 'String')); %320;               % Inlet temperature (K)
%%%%%%%%%%%%%%%%%
% Unburned Fuel-air equation of state, CH4, phe = 1
Ru = str2num(get(handles.ubair_gas_const, 'String')); %0.2968;            % unburned fuel - air gas constant (kJ/kgK)
Cpu = str2num(get(handles.ubair_cp, 'String')); %1.022;            % unburned fuel - air Cv (kJ/kgK)
ku = str2num(get(handles.ubair_sp_ht_ratio, 'String')); %1.29;              % unburned fuel - air specific heat ratio
hfu = str2num(get(handles.ubair_zero_deg_enth, 'String')); %-692;             % Zero degree enthalpy for unburned fuel - air(kJ/kg)
Cvu = Cpu - Ru;         % unburned fuel - air Cp (kJ/kgK)
%%%%%%%%%%%%%%%%%
% Burned Fuel-air equation of state
Rb = str2num(get(handles.bair_gas_const, 'String')); %0.2959;            % burned fuel - air gas constant (kJ/kgK)
Cpb = str2num(get(handles.bair_cp, 'String')); %1.096;            % burned fuel - air Cv (kJ/kgK)
kb = str2num(get(handles.bair_sp_ht_ratio, 'String')); %1.27;              % burned fuel - air specific heat ratio
hfb = str2num(get(handles.bair_zero_deg_enth, 'String')); %-3471;            % Zero degree enthalpy for burned fuel - air(kJ/kg)
Cvb = Cpb - Rb;         % burned fuel - air Cp (kJ/kgK)
%%%%%%%%%%%%%%%%%
% uf - air bf - air mixture data
dcv = Cvb - Cvu;        %
dhf = hfb - hfu;        %
dR = Rb - Ru;           %
%%%%%%%%%%%%
% Initial condition
xr = 0.1;
P1 = 50;
T1 = 600;
for i = 1:3
    th(1) = 0;
    T(1) = T1;
    P(1) = P1;
    m = P1 * Vm/(((1-xr)*Ru+xr*Rb)*T1);
    V(1) = Vm;
    W(1) = 0;
    mc(1) = m;
    y(1) = 0;
    dy(1) = 0;
    Qt = 0;
    %% Initial mass in the cylinder
    %% Numerical Integration parameters
    N = 200;
    dth = thevo/N;
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Numerical Integration of the 1st law until exhaust valve opens
    for j = 2:N+1
        th(j) = j * dth;
        thm = th(j) - dth;
        mc(j) = m;
        [v,dv] = sz_volm(thm,Rc,rc,Vm);
        [x,dx] = sz_comb(thm,xr,ths,thb);
        q = -(dcv*T(j-1)+dhf)/(Cvu+x*dcv);
        s = -((1-x)*Ru+x*Rb)/(Cvu+x*dcv);
        T(j) = T(j-1) + (q*dx + s*T(j-1)*dv/v)*dth;
        W(j) = W(j-1) + P(j-1)*dv*dth;
        [v,dv] = sz_volm(th(j),Rc,rc,Vm);
        V(j) = v;
        P(j) = m*((1-x)*Ru+x*Rb)*T(j)/v;
    end
    PowecycleEff = W(j)/(m*(1-xr)*(hfu-hfb));
    %PowecycleEff;
    j = N+1;
    Avalve(N+1) = 0;
    %%%%%%
    %% Exhaust process
    K = 101 + N;
    dth = (3*pi - thevo)/100;
    for j = N+2:K
        th(j) = th(j-1) + dth;
        Ae = valve(th(j-1),thevo,thevc,eportd,estemd,elm);
        Avalve(j) = Ae;
        mf = sz_massf(P(j-1),T(j-1),Pe,kb,Rb,Ae,ws);
        [v,dv] = sz_volm(th(j-1),Rc,rc,Vm);
        mc(j) = mc(j-1) - mf*dth;
        T(j) = T(j-1) -(kb -1) *(dv/v+mf/mc(j-1))*T(j-1)*dth;
        W(j) = W(j-1) + P(j-1)*dv*dth;
        [v,dv] = sz_volm(th(j),Rc,rc,Vm);
        V(j) = v;
        P(j) = mc(j) * Rb * T(j)/V(j);
        K = j;
    end    
    K = j;
    thivo = th(K);
    mr = mc(K);
    %%%%%%%%%%%%%%%%%%%%
    %% Intake Process
    L = round(((4*pi-th(K))/dth));
    for j = K+1:K+L+1
        th(j) = th(j-1) + dth;
        Ai = valve(th(j-1),thivo,thivc,iportd,istemd,ilm);
        Avalve(j) = Ai;
        mf = sz_massf(Pi,Ti,P(j-1),ku,Ru,Ai,ws);
        [v,dv] = sz_volm(th(j-1),Rc,rc,Vm);
        x = mr/mc(j-1);
        a = Cvu + x * dcv;
        b = (Ru + x * dR)/a;
        T(j) = T(j-1) + (-b*dv*T(j-1)/v+mf*(Cpu*Ti-Cvu*T(j-1))/(mc(j-1)*a))*dth;    % Temperature
        mc(j) = mc(j-1)+mf*dth;
        W(j) = W(j-1) + P(j-1)*dv*dth;
        [v,dv] = sz_volm(th(j),Rc,rc,Vm);
        V(j) = v;
        x = mr/mc(j);
        P(j) = mc(j) * (x *Rb+(1-x)*Ru)* T(j)/v;
    end
    %%% Outputs
    xr = x; %mr/mc(j);
    T1 = T(j);
    P1 = P(j);
    Work = W(j);
    Eff = Work/((hfu-hfb)*(mc(j)-mr));
    Power = rpm*Work/120;
    assignin('base','W',W);
    assignin('base','V',V);
    assignin('base','P',P);
    assignin('base','mc',mc);
    assignin('base','th',th);
    assignin('base','T',T);
end

Contact us