I am writing a code for and HCCI combustion engine and need to plot the effect of changing equivalence ratio in the range 0.3<phi<0.8. I have it set up in a loop with hold on so that it combines the points onto one graph, but my plot shows up empty?

for i = 1:np
    theta(i) = 180+n*CR; %set crank to go from 180deg-540deg
    V(i) = Vc*(1+0.5*(Rc-1)*(2+1-cosd(theta(i))-(4-sind(theta(i))^2)^0.5));
    n = n+1;
end %solves for volume
P(1) = Pbdc;
T(1) = Tbdc;
gamma = 1.4;
%isentropic relation so solve for T and P using the following
set(gas,'T',T(1),'P',P(1),'X',x);
for phi = 0.3:.05:0.8
for i = 2:360/CR
    t = 0;
    dt = 8.33e-5;
    set(gas,'T',T(i-1),'P',P(i-1));
    r = Reactor(gas);
    network = ReactorNet({r});

      for n = 1:10
          t = t+dt/10;
          advance(network,t);
      end

      temp = temperature(r);
      Pr = pressure(r);
      P(i) = Pr*(V(i-1)/V(i))^gamma;
      T(i) = temp*(V(i-1)/V(i))^(gamma-1);

end %Have T and P, now place them in arrays

      subplot(2,1,1);
      plot(phi,T);
      xlabel('Equivalence Ratio (phi)');
      ylabel('Temperature (K)');
      hold on

      subplot(2,1,2);
      plot(phi,P);
      xlabel('Equivalence Ratio (phi)');
      ylabel('Pressure (Pa)');
      hold on
  end  end

Answer 1

[Geoff Hayes]

Jack - without having access to your data it is difficult to say why you are not seeing anything being shown in your two subplots. If I do something trivial like

 for phi = 0.3:.05:0.8
    T = rand(12,1);
    P = rand(12,1);
    subplot(2,1,1);
    plot(phi,T);
    xlabel('Equivalence Ratio (phi)');
    ylabel('Temperature (K)');
    hold on
    subplot(2,1,2);
    plot(phi,P);
    xlabel('Equivalence Ratio (phi)');
    ylabel('Pressure (Pa)');
    hold on
 end

where T and P are random vectors of twelve elements each ("simulating" your inner for loop) then each iteration of the for loop plots something distinct for each phi. So I don't think that there is anything wrong with the way you have set up the code to plot to each of the subplots.

Put a break point at the line

 plot(phi,T);

and run your code. When the debugger pauses at this line, look at T. How many elements are in this array. There should be 360/CR (what is the value for CR?) Is this the case? Now look at the elements themselves. What is the minimum and what is the maximum value in both arrays?

[Jack Walsh]

here is the data:

gas = importPhase('GRI30.xml');
nsp = nSpecies(gas);
ic3h8 = speciesIndex(gas,'C3H8');
x = zeros(nsp,1);
io2 = speciesIndex(gas,'O2');
in2 = speciesIndex(gas,'N2');
phi = 0.5;
x(ic3h8,1) = 1;
x(io2,1) = 5/phi;
x(in2,1) = (5*3.76)/phi;

%function [V,P,T,theta] = project
Tbdc = 550; %temperature at bottom dead center in degrees Kelvin
Pbdc = 101000; %pressure at bottom dead center in pascals
Rc = 14; %compression ratio
Vd = 0.5; %displacement volume in liters
Vc = Vd/(Rc-1);
rps = 33.3333; %complete cycles per second of the engine
tfr = 1/rps; %time it taks for a complete cycle
dt = tfr/360; %time it takes to rotate one degree
%t = 0; %initial time at bottom dead center
%gas = importPhase('GRI30.xml')
CR =1; %crank resolution
np = 360/CR; %number of points to be iterated through
n = 1;

[Walter Roberson]

We do not have your GRI30.xml and importPhase and various other names you mention are not routines that I recognize.

[Jack Walsh]

Oh I am sorry that is because I am using it in unison with the freeware Cantera. I figured out my problem and have since turned in my repot on it! Thanks for your help!

[sebsatian kluge]

Could you tell us how you made it work?