*** empty log message ***
This commit is contained in:
Dave Goodwin
parent
a7c2084a14
commit
ad791e12d5
18 changed files with 138 additions and 62 deletions
|
|
@ -1,5 +1,5 @@
|
|||
function mdot = massFlowRate(f)
|
||||
function mdot = massFlowRate(f, time)
|
||||
% MASSFLOWRATE - mass flow rate in kg/s
|
||||
%
|
||||
mdot = flowdevicemethods(21, f.index)
|
||||
mdot = flowdevicemethods(21, f.index, time);
|
||||
|
||||
|
|
|
|||
|
|
@ -4,20 +4,23 @@ function x = ReactorNet(reactors)
|
|||
% A ReactorNet object is a container that holds one or more
|
||||
% Reactor objects.
|
||||
%
|
||||
% See also: Reservoir
|
||||
%
|
||||
if nargin == 1
|
||||
else
|
||||
error('wrong number of arguments');
|
||||
error('wrong number of arguments to ReactorNet constructor');
|
||||
end
|
||||
|
||||
x.index = reactornetmethods(0);
|
||||
x.index = reactornetmethods(0,0);
|
||||
if x.index < 0
|
||||
error(geterr);
|
||||
end
|
||||
x = class(x,'ReactorNet');
|
||||
|
||||
% add reactors
|
||||
unfinished
|
||||
nr = length(reactors);
|
||||
for i = 1:nr
|
||||
addReactor(x,reactors{i});
|
||||
end
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,17 +1,17 @@
|
|||
function advance(r, tout)
|
||||
function advance(n, tout)
|
||||
% ADVANCE - Advance the state of the reactor network in time.
|
||||
%
|
||||
% Method advance integrates the system of ordinary differential
|
||||
% equations that determine the rate of change of the volume, the
|
||||
% mass of each species, and the total energy for each reactor. The
|
||||
% integration is carried out from the current reactor time to time
|
||||
% integration is carried out from the current time to time
|
||||
% 'tout.' (Note 'tout' is an absolute time, not a time interval.)
|
||||
% The integrator may take many internal time steps before reaching
|
||||
% tout.
|
||||
%
|
||||
% for i in 1:10
|
||||
% tout = 0.1*i
|
||||
% advance(r, tout)
|
||||
% advance(n, tout)
|
||||
% ...
|
||||
% <add output commands here>
|
||||
% ...
|
||||
|
|
@ -19,4 +19,4 @@ function advance(r, tout)
|
|||
%
|
||||
% See also: ReactorNet/step
|
||||
%
|
||||
reactormethods(8, reactornet_hndl(r), tout);
|
||||
reactornetmethods(8, reactornet_hndl(n), tout);
|
||||
|
|
|
|||
|
|
@ -9,7 +9,7 @@ function t = step(r, tout)
|
|||
%
|
||||
% Method 'step' takes one internal time step and returns. This
|
||||
% can be useful when it is desired to resolve a rapidly-changing
|
||||
% solution in the output file.
|
||||
% solution.
|
||||
%
|
||||
% This method can be used as follows:
|
||||
%
|
||||
|
|
@ -24,5 +24,5 @@ function t = step(r, tout)
|
|||
%
|
||||
% See also: Reactor/advance
|
||||
%
|
||||
t = reactormethods(21, reactor_hndl(r), tout);
|
||||
t = reactornetmethods(21, reactornet_hndl(r), tout);
|
||||
|
||||
|
|
|
|||
|
|
@ -1,6 +1,13 @@
|
|||
function v = critTemperature(a)
|
||||
% CRITTEMPERATURE - Critical temperature [K].
|
||||
function v = minTemp(p)
|
||||
% MINTEMP - Minimum temperature.
|
||||
%
|
||||
% The parameterizations used to represent the temperature-dependent
|
||||
% species thermodynamic properties are generally only valid in some
|
||||
% finite temperature range, which may be different for each species
|
||||
% in the phase. This method returns the lowest temperature at which
|
||||
% the parameterizations are valid for all species in the phase.
|
||||
%
|
||||
% The critical temperature is the temperature at the critical point
|
||||
%
|
||||
v = thermo_get(p.tp_id,19);
|
||||
% See also: maxTemp
|
||||
%
|
||||
v = thermo_get(p.tp_id,16);
|
||||
|
||||
|
|
|
|||
|
|
@ -1,7 +1,7 @@
|
|||
function ignite2(g)
|
||||
% IGNITE2 Zero-dimensional kinetics: adiabatic, constant volume.
|
||||
%
|
||||
% This example illustrates how to use function 'reactor' for
|
||||
% This example illustrates how to use function 'reactor_ode' for
|
||||
% zero-dimensional kinetics simulations with arbitrary heat flux
|
||||
% and volume vs. time. Here a constant-volume, adiabatic simulation
|
||||
% is conducted by setting vdot and q to zero.
|
||||
|
|
|
|||
|
|
@ -2,26 +2,29 @@ function periodic_cstr
|
|||
%
|
||||
% Periodic CSTR
|
||||
%
|
||||
% This example illustrates a CSTR with steady inputs but periodic interior state.
|
||||
% A stoichiometric hydrogen/oxygen mixture is introduced and reacts to produce water.
|
||||
% But since water has a large efficiency as a third body in the chain termination reaction
|
||||
% This example illustrates a CSTR with steady inputs but periodic
|
||||
% interior state. A stoichiometric hydrogen/oxygen mixture is
|
||||
% introduced and reacts to produce water. But since water has a
|
||||
% large efficiency as a third body in the chain termination reaction
|
||||
%
|
||||
% H + O2 + M = HO2 + M
|
||||
%
|
||||
% as soon as a significant amount of water is produced the reaction stops. After enough time has
|
||||
% passed that the water is exhausted from the reactor, the mixture explodes again and the
|
||||
% process repeats. This explanation can be verified by decreasing the rate for reaction 7 in
|
||||
% file 'h2o2.xml' and re-running the example.
|
||||
% as soon as a significant amount of water is produced the reaction
|
||||
% stops. After enough time has passed that the water is exhausted from
|
||||
% the reactor, the mixture explodes again and the process
|
||||
% repeats. This explanation can be verified by decreasing the rate for
|
||||
% reaction 7 in file 'h2o2.cti' and re-running the example.
|
||||
%
|
||||
% Acknowledgments: The idea for this example and an estimate of the conditions needed to
|
||||
% see the oscillations came from Bob Kee, Colorado School of Mines
|
||||
% Acknowledgments: The idea for this example and an estimate of the
|
||||
% conditions needed to see the oscillations came from Bob Kee,
|
||||
% Colorado School of Mines
|
||||
%
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
||||
|
||||
help periodic_cstr
|
||||
|
||||
% create the gas mixture
|
||||
gas = IdealGasMix('h2o2.xml');
|
||||
gas = IdealGasMix('h2o2.cti');
|
||||
|
||||
% pressure = 60 Torr, T = 770 K
|
||||
p = 60.0*133.3;
|
||||
|
|
@ -86,8 +89,9 @@ downstream = Reservoir(gas);
|
|||
v = Valve;
|
||||
install(v, cstr, downstream);
|
||||
setValveCoeff(v, 1.0e-9);
|
||||
% mdot (kg/s) = coeff * (delta P, Pascals)
|
||||
|
||||
% create the network
|
||||
network = ReactorNet({cstr});
|
||||
|
||||
% now integrate in time
|
||||
tme = 0.0;
|
||||
|
|
@ -97,12 +101,13 @@ n = 0;
|
|||
while tme < 300.0
|
||||
n = n + 1;
|
||||
tme = tme + dt;
|
||||
advance(cstr, tme);
|
||||
advance(network, tme);
|
||||
tm(n) = tme;
|
||||
y(1,n) = massFraction(cstr,'H2');
|
||||
y(2,n) = massFraction(cstr,'O2');
|
||||
y(3,n) = massFraction(cstr,'H2O');
|
||||
end
|
||||
clf
|
||||
figure(1)
|
||||
plot(tm,y)
|
||||
legend('H2','O2','H2O')
|
||||
|
|
|
|||
|
|
@ -51,6 +51,7 @@ disp(['CPU time = ' num2str(cputime - t0)]);
|
|||
% plot results
|
||||
|
||||
clf;
|
||||
%figure(1);
|
||||
subplot(2,2,1);
|
||||
surf(xo2,t,pr);
|
||||
xlabel('Elemental O/(O+H)');
|
||||
|
|
|
|||
|
|
@ -33,18 +33,21 @@ env = Reservoir(a);
|
|||
w = Wall;
|
||||
install(w,r,env);
|
||||
|
||||
% set expansion parameter. dV/dt = K(P_1 - P_2)
|
||||
% set expansion parameter. dV/dt = KA(P_1 - P_2)
|
||||
setExpansionRateCoeff(w, 1.0e6);
|
||||
|
||||
% set wall area
|
||||
setArea(w, 1.0);
|
||||
|
||||
t = 0;
|
||||
% create a reactor network and insert the reactor:
|
||||
network = ReactorNet({r});
|
||||
|
||||
t = 0.0;
|
||||
dt = 1.0e-5;
|
||||
t0 = cputime;
|
||||
for n = 1:100
|
||||
t = t + dt;
|
||||
advance(r, t);
|
||||
advance(network, t);
|
||||
tim(n) = time(r);
|
||||
temp(n) = temperature(r);
|
||||
x(n,1:3) = moleFraction(gas,{'OH','H','H2'});
|
||||
|
|
|
|||
|
|
@ -22,12 +22,15 @@ set(gas,'T',1001.0,'P',oneatm,'X','H2:2,O2:1,N2:4');
|
|||
% create a reactor, and insert the gas
|
||||
r = Reactor(gas);
|
||||
|
||||
% create a reactor network and insert the reactor
|
||||
network = ReactorNet({r});
|
||||
|
||||
t = 0;
|
||||
dt = 1.0e-5;
|
||||
t0 = cputime;
|
||||
for n = 1:100
|
||||
t = t + dt;
|
||||
advance(r, t);
|
||||
advance(network, t);
|
||||
tim(n) = time(r);
|
||||
temp(n) = temperature(r);
|
||||
x(n,1:3) = moleFraction(gas,{'OH','H','H2'});
|
||||
|
|
|
|||
|
|
@ -15,6 +15,12 @@ pause
|
|||
surfreactor;
|
||||
disp('press any key to continue');
|
||||
pause
|
||||
periodic_cstr;
|
||||
disp('press any key to continue');
|
||||
pause;
|
||||
rankine(300.0, 2.0*oneatm, 0.8, 0.7);
|
||||
disp('press any key to continue');
|
||||
pause;
|
||||
prandtl1(0);
|
||||
disp('press any key to continue');
|
||||
pause
|
||||
|
|
|
|||
|
|
@ -12,6 +12,9 @@ gas = importPhase('ptcombust.cti','gas');
|
|||
% set the initial conditions
|
||||
set(gas,'T',t,'P',oneatm,'X','CH4:0.01, O2:0.21, N2:0.78');
|
||||
|
||||
% The surface reaction mechanism describes catalytic combustion of
|
||||
% methane on platinum, and is from Deutschman et al., 26th
|
||||
% Symp. (Intl.) on Combustion,1996, pp. 1747-1754
|
||||
surf = importInterface('ptcombust.cti','Pt_surf', gas);
|
||||
setTemperature(surf, t);
|
||||
|
||||
|
|
@ -39,22 +42,24 @@ setKinetics(w, surf, 0);
|
|||
|
||||
% set the wall area and heat transfer coefficient.
|
||||
setArea(w, 1.0e-4);
|
||||
setHeatTransferCoeff(w,0.1);
|
||||
setHeatTransferCoeff(w,1.0e1); % W/m2/K
|
||||
|
||||
% set expansion parameter. dV/dt = K(P_1 - P_2)
|
||||
setExpansionRateCoeff(w, 1.0);
|
||||
% set expansion rate parameter. dV/dt = KA(P_1 - P_2)
|
||||
setExpansionRateCoeff(w, 1.0e2);
|
||||
|
||||
network = ReactorNet({r});
|
||||
|
||||
t = 0;
|
||||
dt = 0.1;
|
||||
t0 = cputime;
|
||||
names = {'CH4','O2','CO','CO2','H2O'};
|
||||
x = zeros([100 5]);
|
||||
names = {'CH4','CO','CO2','H2O'};
|
||||
x = zeros([100 4]);
|
||||
for n = 1:100
|
||||
t = t + dt;
|
||||
advance(r, t);
|
||||
advance(network, t);
|
||||
tim(n) = t;
|
||||
temp(n) = temperature(r);
|
||||
pres(n) = pressure(r);
|
||||
cov(n,:) = coverages(surf)';
|
||||
x(n,:) = moleFraction(gas,names);
|
||||
end
|
||||
|
|
@ -66,11 +71,15 @@ plot(tim,temp);
|
|||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
subplot(2,2,2);
|
||||
plot(tim,pres);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Pressure (Pa)');
|
||||
subplot(2,2,3);
|
||||
semilogy(tim,cov);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Coverages');
|
||||
legend(speciesNames(surf));
|
||||
subplot(2,2,3);
|
||||
subplot(2,2,4);
|
||||
plot(tim,x);
|
||||
xlabel('Time (s)');
|
||||
ylabel('Mole Fractions');
|
||||
|
|
|
|||
|
|
@ -1,3 +1,5 @@
|
|||
const double Undef = -999.123;
|
||||
|
||||
void reportError();
|
||||
|
||||
void checkNArgs(const int n, const int nrhs);
|
||||
|
|
|
|||
|
|
@ -110,7 +110,8 @@ extern "C" {
|
|||
case FUNC_CLASS:
|
||||
funcmethods(nlhs, plhs, nrhs, prhs); break;
|
||||
default:
|
||||
mexErrMsgTxt("unknown class");
|
||||
mexPrintf("iclass = %d",iclass);
|
||||
//mexErrMsgTxt("unknown class");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -43,7 +43,7 @@
|
|||
iok = flowdev_install(i, int(v), m);
|
||||
break;
|
||||
case 3:
|
||||
iok = flowdev_setSetpoint(i, v);
|
||||
iok = flowdev_setMassFlowRate(i, v);
|
||||
break;
|
||||
case 4:
|
||||
iok = flowdev_setParameters(i, 1, &v);
|
||||
|
|
@ -54,6 +54,9 @@
|
|||
case 6:
|
||||
iok = flowdev_ready(i);
|
||||
break;
|
||||
case 7:
|
||||
iok = flowdev_setMaster(i, int(v));
|
||||
break;
|
||||
default:
|
||||
mexErrMsgTxt("unknown job parameter");
|
||||
}
|
||||
|
|
@ -70,10 +73,7 @@
|
|||
else if (job < 40) {
|
||||
switch (job) {
|
||||
case 21:
|
||||
r = flowdev_massFlowRate(i);
|
||||
break;
|
||||
case 22:
|
||||
r = flowdev_setpoint(i);
|
||||
r = flowdev_massFlowRate(i, v);
|
||||
break;
|
||||
default:
|
||||
mexErrMsgTxt("unknown job parameter");
|
||||
|
|
|
|||
|
|
@ -4,8 +4,6 @@
|
|||
#include "../../../clib/src/ct.h"
|
||||
#include "ctmatutils.h"
|
||||
|
||||
//const double Undef = -999.123;
|
||||
|
||||
void reactormethods( int nlhs, mxArray *plhs[],
|
||||
int nrhs, const mxArray *prhs[] )
|
||||
{
|
||||
|
|
@ -46,18 +44,18 @@
|
|||
case 4:
|
||||
iok = reactor_setInitialVolume(i, v);
|
||||
break;
|
||||
// case 5:
|
||||
//iok = reactor_setInitialTime(i, v);
|
||||
//break;
|
||||
case 5:
|
||||
iok = reactor_setInitialTime(i, v);
|
||||
break;
|
||||
case 6:
|
||||
iok = reactor_setThermoMgr(i, int(v));
|
||||
break;
|
||||
case 7:
|
||||
iok = reactor_setKineticsMgr(i, int(v));
|
||||
break;
|
||||
//case 8:
|
||||
//iok = reactor_advance(i, v);
|
||||
//break;
|
||||
case 8:
|
||||
iok = reactor_advance(i, v);
|
||||
break;
|
||||
case 9:
|
||||
iok = reactor_setEnergy(i, int(v));
|
||||
break;
|
||||
|
|
@ -76,9 +74,9 @@
|
|||
|
||||
else if (job < 40) {
|
||||
switch (job) {
|
||||
//case 21:
|
||||
//r = reactor_step(i, v);
|
||||
//break;
|
||||
case 21:
|
||||
r = reactor_step(i, v);
|
||||
break;
|
||||
case 22:
|
||||
r = reactor_time(i);
|
||||
break;
|
||||
|
|
@ -116,3 +114,4 @@
|
|||
return;
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -17,7 +17,9 @@
|
|||
|
||||
double r = Undef;
|
||||
double v = Undef;
|
||||
double v2 = -1.0;
|
||||
if (nrhs > 3) v = getDouble(prhs[3]);
|
||||
if (nrhs > 4) v2 = getDouble(prhs[4]);
|
||||
|
||||
// constructor
|
||||
if (job == 0) {
|
||||
|
|
@ -45,10 +47,16 @@
|
|||
break;
|
||||
case 4:
|
||||
iok = reactornet_addreactor(i, int(v));
|
||||
break;
|
||||
break;
|
||||
case 5:
|
||||
iok = reactornet_setInitialTime(i, v);
|
||||
break;
|
||||
case 6:
|
||||
iok = reactornet_setMaxTimeStep(i, v);
|
||||
break;
|
||||
case 7:
|
||||
iok = reactornet_setTolerances(i, v, v2);
|
||||
break;
|
||||
case 8:
|
||||
iok = reactornet_advance(i, v);
|
||||
break;
|
||||
|
|
@ -61,4 +69,30 @@
|
|||
if (iok < 0) reportError();
|
||||
return;
|
||||
}
|
||||
|
||||
// options that return a value of type 'double'
|
||||
|
||||
else if (job < 40) {
|
||||
switch (job) {
|
||||
case 21:
|
||||
r = reactornet_step(i, v);
|
||||
break;
|
||||
case 22:
|
||||
r = reactornet_time(i);
|
||||
break;
|
||||
case 23:
|
||||
r = reactornet_rtol(i);
|
||||
break;
|
||||
case 24:
|
||||
r = reactornet_atol(i);
|
||||
break;
|
||||
default:
|
||||
mexErrMsgTxt("unknown job parameter");
|
||||
}
|
||||
plhs[0] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL);
|
||||
double *h = mxGetPr(plhs[0]);
|
||||
*h = r;
|
||||
if (r == Undef) reportError();
|
||||
return;
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -63,7 +63,7 @@
|
|||
iok = wall_setExpansionRateCoeff(i, v);
|
||||
break;
|
||||
case 10:
|
||||
iok = wall_setExpansionRate(i, int(v));
|
||||
iok = wall_setVelocity(i, int(v));
|
||||
break;
|
||||
case 11:
|
||||
iok = wall_ready(i);
|
||||
|
|
@ -73,6 +73,9 @@
|
|||
m = getInt(prhs[4]);
|
||||
iok = wall_setkinetics(i, n, m);
|
||||
break;
|
||||
case 13:
|
||||
iok = wall_setEmissivity(i, v);
|
||||
break;
|
||||
default:
|
||||
mexErrMsgTxt("unknown job parameter");
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue