# periodic_cstr.m (Source)

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106``` ```function periodic_cstr % % 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.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 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% help periodic_cstr % create the gas mixture gas = IdealGasMix('h2o2.cti'); % pressure = 60 Torr, T = 770 K p = 60.0*133.3; t = 770.0; OneAtm = 1.01325e5; set(gas,'T', 300.0, 'P', p, 'X', 'H2:2, O2:1'); % create an upstream reservoir that will supply the reactor. The % temperature, pressure, and composition of the upstream reservoir are % set to those of the 'gas' object at the time the reservoir is % created. upstream = Reservoir(gas); % Now set the gas to the initial temperature of the reactor, and create % the reactor object. set(gas, 'T', t, 'P', p); cstr = IdealGasReactor(gas); % Set its volume to 10 cm^3. In this problem, the reactor volume is % fixed, so the initial volume is the volume at all later times. setInitialVolume(cstr, 10.0*1.0e-6); % We need to have heat loss to see the oscillations. Create a % reservoir to represent the environment, and initialize its % temperature to the reactor temperature. env = Reservoir(gas); % Create a heat-conducting wall between the reactor and the % environment. Set its area, and its overall heat transfer % coefficient. Larger U causes the reactor to be closer to isothermal. % If U is too small, the gas ignites, and the temperature spikes and % stays high. w = Wall; install(w, cstr, env); setArea(w, 1.0); setHeatTransferCoeff(w, 0.02); % Connect the upstream reservoir to the reactor with a mass flow % controller (constant mdot). Set the mass flow rate to 1.25 sccm. sccm = 1.25; vdot = sccm * 1.0e-6/60.0 * ((OneAtm / pressure(gas)) * ( temperature(gas) / 273.15)); % m^3/s mdot = density(gas) * vdot; % kg/s mfc = MassFlowController; install(mfc, upstream, cstr); setMassFlowRate(mfc, mdot); % now create a downstream reservoir to exhaust into. downstream = Reservoir(gas); % connect the reactor to the downstream reservoir with a valve, and % set the coefficient sufficiently large to keep the reactor pressure % close to the downstream pressure of 60 Torr. v = Valve; install(v, cstr, downstream); setValveCoeff(v, 1.0e-9); % create the network network = ReactorNet({cstr}); % now integrate in time tme = 0.0; dt = 0.1; n = 0; while tme < 300.0 n = n + 1; tme = tme + dt; 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') title('Mass Fractions') ```