Heat release rate in MATLAB

[Maciej Mikulski]

 Hi,

 I'm new with CANTERA and I'm running it as matlab toolbox. I'm trying to start with CANTERA for CH4 combustion using GRI30 mechanism.

How to calculate total heat release rate from the reactions in matlab? I can find creationRates of the species.

 

Appreciate any advice

[Thomas Fiala]

Hi Maciej,

In which context do you want to compute the heat release rate?

For one-dimensional flames, maybe these threads using the Python module will help you:

https://groups.google.com/forum/?fromgroups#!searchin/cantera-users/heat$20release/cantera-users/Qr24bK7Qg3k/yebXO8JxPj0J
https://groups.google.com/forum/?fromgroups#!searchin/cantera-users/heat$20release/cantera-users/nrHehTJVhGA/AOrOpJHTxDkJ

The functions described there use the Python module, but the underlying functions probably have the same names in the Matlab toolbox.

Thomas

[Maciej Mikulski]

Thank You Tomas,

I've managed with the heat release, acording to the script You have send. I have a nother question though. After calculating, heat release I'm calculating mass fractions change

dY/dt from the solutions in GRI30 (doty function). How to calculate mas fraction after time step dt. Should I do with from the definition of the derivative like Y(t+dt)=(dY/dt)*dt+Y(t), or is there a dedicated fuction for this?

Maciej     

[Nick Curtis]

Hi Maciej,
In order to integrate the mass fractions / temperature / pressure over time, you should look into creating a Reactor that fits the model you want to simulate (e.g. constant pressure, volume, etc.)
Check out the examples on this page, especially the reactor1.m and reactor2.m files

Nick

[Maciej Mikulski]

Hi Nick, thank you for the advice,

 I expected that this is done by the reactors solutions. But how about using it in your own "reactor". I've tried to implement the kinetics in my own, simple zero-d CI engine model. The energy equations are solved by iteration, step, by step in Crank angle (CA) domain (but this is just scaling from time domain).

 My first, general idea to implement Cantera kinetics to the code was rather simple:

 1. input parameters to cantera: 

set(Gas,'Pressure',P(i),'Temperature',T(i),'Y',Y(i)) 

 2. Calculate using CANTERA solutions:

-internal energy of the system U=cv_mole(Gas).*moleFractions(Gas)*N     %N is a total number of moles

-heat release rate from combustion  : dQ/da=(30/(pi*n)*enthalpies_RT(Gas).*netProdRates(Gas)*R*T(i)*V(i);   %(30/(pi*n) is just scaling from t->CA domain.

-change of mass fractions dY/da=(30/(pi*n)*ydot(Gas)

-mass fractions after calculation step  Y(i+1)=(dY/da)*da+Y(i)

 3. Calculate using own model state parameters after calculation step (taking to account internal energy U, heat release rate dQ/da from the above calculations) 

P(i+1), T(i+1), V(i+1)

 4. Repeat from point 1 with new values of P, T, Y.

 The algorithm goes well if I set dY/da=0 arbitrarily (i.e. for compression), otherwise after few steps i get negative values of fractions Y and as a consequence negative temperatures (the calculation collapses).

 I suspect that what I’m doing wrong is using set(Gas,'Pressure',P(i),'Temperature',T(i),'Y',Y(i)) command every time to set new parameters in the loop step. Probably I cannot do it simply by this command because it changes other parameters of the gas composition. 

 Could you advise if my suspicions are correct and advise how to solve this problem in a different way.

 Regards,

Maciek  

[Nick Curtis]

Hi Maciek,

The best part about Cantera is that all this can be accomplished very easily!
Basically you'll want to proceed as follows:

1. Install the gas at its initial state in an IdealGasReactor (see here)
2. Create a reservoir, the state of the gas in this reservoir will not change (and doesn't matter) (see here)
3. Install a Wall between the reactor and the Reservoir
4. Derive a function for the velocity of your piston, in terms of time (use the same transform you had earlier). 
   
5. Create a Func object that describes your velocity function.  Note that you can create composite Func objects using the plus, rdivide and times functions (described on the Func page)
6. Use this Func object to set the wall velocity
7. Install the whole thing in a reactor net, e.g. if your IdealGasReactor is named r:
   

   network = ReactorNet({r});

Finally you simply need to step the ReactorNet between your initial and ending crank angle using the step() or advance() functions

Nick