diff --git a/samples/matlab_experimental/PFR.m b/samples/matlab_experimental/PFR.m new file mode 100644 index 000000000..918ceb0cf --- /dev/null +++ b/samples/matlab_experimental/PFR.m @@ -0,0 +1,168 @@ +% General_Plug_Flow_Reactor (PFR) - to solve PFR equations for reactors +% +% This code snippet is to model a constant area and varying area +% (converging and diverging) nozzle as Plug Flow Reactor with given +% dimensions and an incoming gas. The pressure is not assumed to be +% constant here, as opposed to the Python Version of the +% Plug Flow Reactor model. +% +% The reactor assumes that the flow follows the Ideal Gas Law. +% +% The governing equations used in this code can be referenced at: +% *S.R Turns, An Introduction to Combustion - Concepts and Applications, +% McGraw Hill Education, India, 2012, 206-210.* +% +% The current example is written for methane combustion, but can be readily +% adapted for other chemistries. +% +% Developed by Ashwin Kumar/Dr.Joseph Meadows (mgak@vt.edu/jwm84@vt.edu) on 3-June-2020 +% Research Assistant/Assistant Professor +% Advanced Propulsion and Power Laboratory +% Virginia Tech + +%% Clear all variables, close all figures, clear the command line: +clear all +close all +cleanup +clc + +help PFR + +% Temperature of gas, in K +T0 = 1473; +% Pressure of gas, in Pa +P0 = 4.47*101325; + +% Equivalence Ratio +Phi = 0.2899; + +% Import the gas phase, read out key species indices: +gas_calc = Solution('gri30.yaml'); +ich4 = gas_calc.speciesIndex('CH4'); +io2 = gas_calc.speciesIndex('O2'); +in2 = gas_calc.speciesIndex('N2'); +nsp = gas_calc.nSpecies; +x = zeros(nsp,1); + +% Change the below values for different Phi values of methane Combustion +x(ich4,1) = Phi; +x(io2,1) = 2.0; +x(in2,1) = 7.52; + +% Set the initial state and then equilibrate for a given enthalpy and pressure: +gas_calc.TPX = {T0, P0, x}; +gas_calc.equilibrate('HP'); + +%% Calculation of properties along the reactor length +% The Dimensions and conditions of the reactor are given below + +% Inlet Area, in m^2 +A_in = 0.018; +% Exit Area, in m^2 +A_out = 0.003; +% Length of the reactor, in m +L = 1.284*0.0254; +% The whole reactor is divided into n small reactors +n = 100; +% Mass flow rate into the reactor, in kg/s +mdot_calc = 1.125; + +% Flag to indicate whether the area is converging, diverging, or constant: +% k = -1 makes the solver solve for converging area. +% k = +1 makes the solver solve for diverging area. +% k = 0 makes the solver solve for constant cross sectional area +if A_in > A_out + k = -1; +elseif A_out > A_in + k = 1; +else + k = 0; +end + +dAdx = abs(A_in-A_out)/L; +% The whole length of the reactor is divided into n small lengths +dx = L/n; + +x_calc = 0:dx:L; +nsp = gas_calc.nSpecies; + +% Initialize arrays for T, Y, and rho at each location: +T_calc = zeros(length(x_calc), 1); +Y_calc = zeros(length(x_calc), nsp); +rho_calc = zeros(length(x_calc), 1); + +T_calc(1) = gas_calc.T; +Y_calc(1,:) = gas_calc.Y; +rho_calc(1) = gas_calc.D; + +for i = 2:length(x_calc) + + %Solver location indicator + fprintf('Solving reactor %d of %d\n', i, length(x_calc)) + +%-------------------------------------------------------------------------- +%------The values of variables at previous location are given as initial--- +%------values to the current iteration and the limits of the current------- +%--------------reactor and the gas entering it are being set--------------- +%-------------------------------------------------------------------------- + inlet_soln(1) = rho_calc(i-1); + inlet_soln(2) = T_calc(i-1); + inlet_soln(3:nsp+2) = Y_calc(i-1, :); + limits = [x_calc(i-1), x_calc(i)]; + gas_calc.TDY = {T_calc(i-1), rho_calc(i-1), Y_calc(i-1, :)}; + options = odeset('RelTol', 1.e-10, 'AbsTol', 1e-10,... + 'InitialStep', 1e-8, 'NonNegative', 1); +%-------------------------------------------------------------------------- +%-------------------------------------------------------------------------- + % These values are passed onto the ode15s solver + [~,y] = ode15s(@PFR_Solver, limits, inlet_soln, options,..., + gas_calc, mdot_calc, A_in, dAdx, k); + + T_calc(i) = y(end, 2); + rho_calc(i) = y(end, 1); + Y_calc(i,:) = y(end, 3:nsp+2); +end + +A_calc = A_in + k.* x_calc * dAdx; +vx_calc = zeros(length(x_calc), 1); +R_calc = zeros(length(x_calc), 1); +M_calc = zeros(length(x_calc), 1); +P_calc = zeros(length(x_calc), 1); +for i=1:length(x_calc) + % The gas is set to the solved property values at each location + gas.TDY = {T_calc(i), rho_calc(i), Y_calc(i, :)}; + % Velocity is calculated from Mass flow rate, Area and Density + vx_calc(i) = mdot_calc./ (A_calc(i) * rho_calc(i)); + % Specific Gas Constant + R_calc(i) = gasconstant() / gas_calc.meanMolecularWeight; + % Mach No. is calculated from local velocity and local speed of sound + M_calc(i) = vx_calc(i) / gas_calc.soundspeed; + % Pressure is calculated from density, temperature and gas constant + P_calc(i) = rho_calc(i) * R_calc(i) * T_calc(i); +end + +%% Plotting +plot(x_calc, M_calc) +xlabel('X-Position (m)') +ylabel('Mach No') +title('Mach No. Variation') +figure(2) +plot(x_calc, A_calc) +xlabel('X-Position (m)') +ylabel('Area (m^2)') +title('Reactor Profile') +figure(3) +plot(x_calc, T_calc) +xlabel('X-Position (m)') +ylabel('Temperature') +title('Temperature Variation') +figure(4) +plot(x_calc, rho_calc) +xlabel('X-Position (m)') +ylabel('Density (kg/m^3)') +title('Density Variation') +figure(5) +plot(x_calc, P_calc) +xlabel('X-Position (m)') +ylabel('Pressure (Pa)') +title('Pressure Variation') diff --git a/samples/matlab_experimental/PFR_Solver.m b/samples/matlab_experimental/PFR_Solver.m new file mode 100644 index 000000000..270c23268 --- /dev/null +++ b/samples/matlab_experimental/PFR_Solver.m @@ -0,0 +1,53 @@ +function F = PFR_Solver(x, soln_vector, gas, mdot, A_in, dAdx, k) + % This function defines the spatial derivatives for an ideal gas plug-flow + % reactor, where the cross-sectional area and pressure are allowed to vary, + % axially. The model is set up by the example file 'PFR.m', + % which points the integrator to this function. The integrator integrates the + % derivatives spatially, to solve the density, temperature, and species mass + % fraction profiles as a function of distance x. + + rho = soln_vector(1); + T = soln_vector(2); + Y = soln_vector(3: end); + + if k==1 + A = A_in+k * dAdx * x; + elseif k==-1 + A = A_in + k * dAdx * x; + dAdx = -dAdx; + else + A = A_in + k * dAdx * x; + end + + % the gas is set to the corresponding properties during each iteration of the ode loop + gas.TDY = {T, rho, Y}; + + MW_mix = gas.meanMolecularWeight; + Ru = gasconstant; + R = Ru / MW_mix; + nsp = gas.nSpecies; + vx = mdot / (rho * A); + P = rho * R * T; + + gas.basis = 'mass'; + MW = gas.MolecularWeights; + h = gas.enthalpies_RT.*R.*T; + w = gas.netProdRates; + Cp = gas.cp; + %-------------------------------------------------------------------------- + %---F(1), F(2) and F(3:end) are the differential equations modelling the--- + %---density, temperature and mass fractions variations along a plug flow--- + %-------------------------reactor------------------------------------------ + %-------------------------------------------------------------------------- + F(1) = ((1-R/Cp)*((rho*vx)^2)*(1/A)*(dAdx)... + + rho*R*sum(MW.*w.*(h-MW_mix*Cp*T./MW))/(vx*Cp) )... + / (P*(1+vx^2/(Cp*T)) - rho*vx^2); + + F(2) = (vx*vx/(rho*Cp))*F(1) + vx*vx*(1/A)*(dAdx)/Cp... + - (1/(vx*rho*Cp))*sum(h.*w.*MW); + + F(3:nsp+2) = w(1:nsp).*MW(1:nsp)./(rho*vx); + + F = F'; + +end diff --git a/samples/matlab_experimental/equil.m b/samples/matlab_experimental/equil.m new file mode 100644 index 000000000..1d4016650 --- /dev/null +++ b/samples/matlab_experimental/equil.m @@ -0,0 +1,63 @@ +function equil(g) +% EQUIL A chemical equilibrium example. +% +% This example computes the adiabatic flame temperature and equilibrium +% composition for a methane/air mixture as a function of equivalence ratio. + help equil + + if nargin == 1 + gas = g; + else + gas = Solution('gri30.yaml'); + end + + nsp = gas.nSpecies; + + % find methane, nitrogen, and oxygen indices + ich4 = gas.speciesIndex('CH4'); + io2 = gas.speciesIndex('O2'); + in2 = gas.speciesIndex('N2'); + + nPhis = 50; + phi = linspace(0.2, 2.70, nPhis); + tad(nPhis) = 0; + xeq(nsp, nPhis) = 0; + + for i = 1:nPhis + x = zeros(nsp, 1); + x(ich4, 1) = phi(i); + x(io2, 1) = 2.0; + x(in2, 1) = 7.52; + gas.T = 300; + gas.P = 101325; + gas.X = x; + gas.equilibrate('HP'); + tad(i) = gas.T; + xeq(:, i) = gas.X; + end + + % make plots + clf; + subplot(1, 2, 1); + plot(phi, tad); + xlabel('Equivalence Ratio'); + ylabel('Temperature (K)'); + title('Adiabatic Flame Temperature'); + + subplot(1,2,2); + semilogy(phi, xeq); + axis([phi(1), phi(50), 1.0e-14, 1]); + %legend(speciesName(gas,1:nsp),1); + j = 10; + for k = 1:nsp + text(phi(j), 1.5*xeq(k, j), gas.speciesName(k)) + j = j + 2; + if j > 46 + j = 10; + end + end + xlabel('Equivalence Ratio'); + ylabel('Mole Fraction'); + title('Equilibrium Composition'); + +end diff --git a/samples/matlab_experimental/flame.m b/samples/matlab_experimental/flame.m new file mode 100644 index 000000000..8372dde3e --- /dev/null +++ b/samples/matlab_experimental/flame.m @@ -0,0 +1,73 @@ +function f = flame(gas, left, flow, right) +% FLAME - create a flame object. +% + % Check input parameters + if nargin ~= 4 + error('wrong number of input arguments.'); + end + + if ~gas.thermo.isIdealGas + error('gas object must represent an ideal gas mixture.'); + end + if ~isInlet(left) + error('burner object of wrong type.'); + end + if ~isFlow(flow) + error('flow object of wrong type.'); + end + + flametype = 0; + if isSurface(right) + flametype = 1; + elseif isInlet(right) + flametype = 3; + end + + % create the container object + f = Stack([left flow right]); + + % set default initial profiles. + rho0 = gas.D; + + % find the adiabatic flame temperature and corresponding + % equilibrium composition + gas.equilibrate('HP'); + teq = gas.T; + yeq = gas.Y; + + z1 = 0.2; + mdot0 = left.massFlux; + mdot1 = right.massFlux; + t0 = left.T; + if flametype == 0 + t1 = teq; + mdot1 = -mdot0; + else + t1 = right.T; + end + f.setProfile(2, {'u', 'V'}, [0.0 1.0 + mdot0/rho0 -mdot1/rho0 + 0.0 0.0]); + f.setProfile(2, 'T', [0.0 z1 1.0; t0 2000.0 t1]); + + for n = 1:gas.nSpecies + nm = gas.speciesName(n); + if strcmp(nm, 'H') || strcmp(nm, 'OH') || strcmp(nm, 'O') ||... + strcmp(nm,'HO2') + yint = 1.0*yeq(n); + else + yint = yeq(n); + end + if flametype == 3 + y1 = right.massFraction(n); + else + y1 = yeq(n); + end + f.setProfile(2, nm, [0, z1, 1 + left.massFraction(n), yint, y1]); + end + + % set minimal grid refinement criteria + f.setRefineCriteria(2, 10.0, 0.8, 0.8); + +end diff --git a/samples/matlab_experimental/ignite.m b/samples/matlab_experimental/ignite.m new file mode 100644 index 000000000..296598c88 --- /dev/null +++ b/samples/matlab_experimental/ignite.m @@ -0,0 +1,103 @@ +function plotdata = ignite(g) +% IGNITE Zero-dimensional kinetics: adiabatic, constant pressure. +% +% This example solves the same problem as 'reactor1,' but does +% it using one of MATLAB's ODE integrators, rather than using the +% Cantera Reactor class. +% + + help ignite + + if nargin == 1 + gas = g; + else + gas = Solution('gri30.yaml'); + end + + % set the initial conditions + + gas.TPX = {1001.0, oneatm, 'X','H2:2,O2:1,N2:4'}; + gas.basis = 'mass'; + y0 = [gas.U + 1.0/gas.D + gas.Y']; + + time_interval = [0 0.001]; + options = odeset('RelTol',1.e-5,'AbsTol',1.e-12,'Stats','on'); + + t0 = cputime; + out = ode15s(@reactor_ode,time_interval,y0,options,gas,@vdot,@area,@heatflux); + disp(['CPU time = ' num2str(cputime - t0)]); + + plotdata = output(out,gas); + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % the functions below may be defined arbitrarily to set the reactor + % boundary conditions - the rate of change of volume, the heat + % flux, and the area. + + % Rate of change of volume. Any arbirtrary function may be implemented. + % Input arguments: + % t time + % vol volume + % gas ideal gas object + + function v = vdot(t, vol, gas) + %v = 0.0; %uncomment for constant volume + v = 1.e11 * (gas.P - 101325.0); % holds pressure very + % close to 1 atm + end + % heat flux (W/m^2). + function q = heatflux(t, gas) + q = 0.0; % adiabatic + end + % surface area (m^2). Used only to compute heat transfer. + function a = area(t,vol) + a = 1.0; + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + % Since the solution variables used by the 'reactor' function are + % not necessarily those desired for output, this function is called + % after the integration is complete to generate the desired + % outputs. + + function pv = output(s, gas) + times = s.x; + soln = s.y; + [~,n] = size(times); + pv = zeros(gas.nSpecies + 4, n); + + gas.TP = {1001.0, oneatm}; + + for j = 1:n + ss = soln(:,j); + y = ss(3:end); + mass = sum(y); + u_mass = ss(1)/mass; + v_mass = ss(2)/mass; + gas.Y = y; + gas.UV = {u_mass, v_mass}; + + pv(1,j) = times(j); + pv(2,j) = gas.T; + pv(3,j) = gas.D; + pv(4,j) = gas.P; + pv(5:end,j) = y; + end + + % plot the temperature and OH mass fractions. + figure(1); + plot(pv(1, :), pv(2, :)); + xlabel('time'); + ylabel('Temperature'); + title(['Final T = ' num2str(pv(2, end)) ' K']); + + figure(2); + ioh = gas.speciesIndex('OH'); + plot(pv(1,:),pv(4+ioh,:)); + xlabel('time'); + ylabel('Mass Fraction'); + title('OH Mass Fraction'); + end +end diff --git a/samples/matlab_experimental/isentropic.m b/samples/matlab_experimental/isentropic.m new file mode 100644 index 000000000..ffc8791d4 --- /dev/null +++ b/samples/matlab_experimental/isentropic.m @@ -0,0 +1,59 @@ +function isentropic(g) +% ISENTROPIC isentropic, adiabatic flow example +% +% In this example, the area ratio vs. Mach number curve is +% computed for a hydrogen/nitrogen gas mixture. +% + help isentropic + + if nargin == 1 + gas = g; + else + gas = Solution('gri30.yaml'); + end + + % set the stagnation state + gas.TPX = {1200.0, 10.0 * oneatm, 'H2:1,N2:0.1'}; + gas.basis = 'mass'; + s0 = gas.S; + h0 = gas.H; + p0 = gas.P; + + mdot = 1; % arbitrary + + mach = []; + a = []; + i = 1; + amin = 1.e14; + + % compute values for a range of pressure ratios + for r = 0.005:0.0025:0.995 + p = p0*r; + + % set the state using (p,s0) + gas.SP = {s0, p}; + + h = gas.H; + rho = gas.D; + + v2 = 2.0*(h0 - h); % h + V^2/2 = h0 + v = sqrt(v2); + a(i) = mdot/(rho*v); % rho*v*A = constant + + if a(i) < amin + amin = a(i); + end + mach(i) = v/gas.soundspeed; + i = i + 1; + end + + a = a/amin; + + % plot results + + clf; + plot(mach,a); + ylabel('Area Ratio'); + xlabel('Mach Number'); + title('Isentropic Flow: Area Ratio vs. Mach Number'); +end diff --git a/samples/matlab_experimental/lithium_ion_battery.m b/samples/matlab_experimental/lithium_ion_battery.m new file mode 100644 index 000000000..8140afcde --- /dev/null +++ b/samples/matlab_experimental/lithium_ion_battery.m @@ -0,0 +1,142 @@ +% This example file calculates the cell voltage of a lithium-ion battery +% at given temperature, pressure, current, and range of state of charge (SOC). +% +% The thermodynamics are based on a graphite anode and a LiCoO2 cathode, +% modeled using the 'BinarySolutionTabulatedThermo' class. +% Further required cell parameters are the electrolyte ionic resistance, the +% stoichiometry ranges of the active materials (electrode balancing), and the +% surface area of the active materials. +% +% The functionality of this example is presented in greater detail in the +% reference (which also describes the derivation of the +% BinarySolutionTabulatedThermo class): +% +% Reference: +% M. Mayur, S. C. DeCaluwe, B. L. Kee, W. G. Bessler, “Modeling and simulation +% of the thermodynamics of lithium-ion battery intercalation materials in the +% open-source software Cantera,” Electrochim. Acta 323, 134797 (2019), +% https://doi.org/10.1016/j.electacta.2019.134797 + +% ----------------------------------------------------------------------------- +% Input +% ----------------------------------------------------------------------------- + +clear all +close all +cleanup +clc + +% Operation parameters +SOC = 0:0.02:1; % [-] Input state of charge (0...1) (can be a vector) +I_app = -1; % [A] Externally-applied current, negative for discharge +T = 293; % [K] Temperature +P = oneatm; % [Pa] Pressure + +% Cell properties +inputFile = 'lithium_ion_battery.yaml'; % Cantera input file name +R_elyt = 0.0384; % [Ohm] Electrolyte resistance +S_ca = 1.1167; % [m^2] Cathode total active material surface area +S_an = 0.7824; % [m^2] Anode total active material surface area + +% Electrode balancing: The "balancing" of the electrodes relates the chemical +% composition (lithium mole fraction in the active materials) to the macroscopic +% cell-level state of charge. +X_Li_an_0 = 0.01; % [-] anode Li mole fraction at SOC = 0 % +X_Li_an_1 = 0.75; % [-] anode Li mole fraction at SOC = 100 % +X_Li_ca_0 = 0.99; % [-] cathode Li mole fraction at SOC = 0 % +X_Li_ca_1 = 0.49; % [-] cathode Li mole fraction at SOC = 100 % + +% ----------------------------------------------------------------------------- +% Calculations +% ----------------------------------------------------------------------------- + +% Calculate mole fractions from SOC +X_Li_an = (X_Li_an_1-X_Li_an_0)*SOC+X_Li_an_0; % anode balancing +X_Li_ca = (X_Li_ca_0-X_Li_ca_1)*(1-SOC)+X_Li_ca_1; % cathode balancing + +% Import all Cantera phases +anode = Solution(inputFile, 'anode'); +cathode = Solution(inputFile, 'cathode'); +elde = Solution(inputFile, 'electron'); +elyt = Solution(inputFile, 'electrolyte'); +anode_interface = Interface(inputFile, 'edge_anode_electrolyte', anode, elde, elyt); +cathode_interface = Interface(inputFile, 'edge_cathode_electrolyte', cathode, elde, elyt); + +% Set the temperatures and pressures of all phases +anode.TP = {T, P}; +cathode.TP = {T, P}; +elde.TP = {T, P}; +elyt.TP = {T, P}; +anode_interface.TP = {T, P}; +cathode_interface.TP = {T, P}; + +% Calculate cell voltage, separately for each entry of the input vectors +V_cell = zeros(length(SOC),1); +phi_l_an = 0; +phi_s_ca = 0; +for i = 1:length(SOC) + % Set anode electrode potential to 0 + phi_s_an = 0; + + % Calculate anode electrolyte potential + phi_l_an = fzero(@(E) anode_curr(phi_s_an, E, X_Li_an(i), anode, elde,... + elyt, anode_interface, S_an) - I_app, phi_l_an); + + % Calculate cathode electrolyte potential + phi_l_ca = phi_l_an + I_app*R_elyt; + + % Calculate cathode electrode potential + phi_s_ca = fzero(@(E) cathode_curr(E, phi_l_ca, X_Li_ca(i), ... + cathode, elde, elyt, cathode_interface,... + S_ca) - I_app, phi_s_ca); + + % Calculate cell voltage + V_cell(i) = phi_s_ca - phi_s_an; +end + +% Let's plot the cell voltage, as a function of the state of charge: +figure(1); +plot(SOC*100, V_cell, 'linewidth', 2.5) +ylim([2.5, 4.3]) +xlabel('State of charge / %') +ylabel('Cell voltage / V') +set(gca, 'fontsize', 14) + + +%-------------------------------------------------------------------------- +% Helper functions +% ----------------------------------------------------------------------------- + +% This function returns the Cantera calculated anode current (in A) +function anCurr = anode_curr(phi_s, phi_l, X_Li_an, anode, elde, elyt, anode_interface, S_an) + % Set the active material mole fraction + anode.X = ['Li[anode]:' num2str(X_Li_an) ', V[anode]:' num2str(1-X_Li_an)]; + + % Set the electrode and electrolyte potential + elde.setElectricPotential(phi_s); + elyt.setElectricPotential(phi_l); + + % Get the net reaction rate at the anode-side interface + % Reaction according to cti file: Li+[elyt] + V[anode] + electron <=> Li[anode] + r = anode_interface.rop_net; % [kmol/m2/s] + + % Calculate the current. Should be negative for cell discharge. + anCurr = r*faradayconstant*S_an; % +end + +% This function returns the Cantera calculated cathode current (in A) +function caCurr = cathode_curr(phi_s, phi_l, X_Li_ca, cathode, elde, elyt, cathode_interface, S_ca) + % Set the active material mole fractions + cathode.X = ['Li[cathode]:' num2str(X_Li_ca) ', V[cathode]:' num2str(1-X_Li_ca)]; + + % Set the electrode and electrolyte potential + elde.setElectricPotential(phi_s); + elyt.setElectricPotential(phi_l); + + % Get the net reaction rate at the cathode-side interface + % Reaction according to cti file: Li+[elyt] + V[cathode] + electron <=> Li[cathode] + r = cathode_interface.rop_net; % [kmol/m2/s] + + % Calculate the current. Should be negative for cell discharge. + caCurr = r*faradayconstant*S_ca*(-1); % +end diff --git a/samples/matlab_experimental/periodic_cstr.m b/samples/matlab_experimental/periodic_cstr.m new file mode 100644 index 000000000..1d7e90ccb --- /dev/null +++ b/samples/matlab_experimental/periodic_cstr.m @@ -0,0 +1,105 @@ +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.yaml' 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 = Solution('h2o2.yaml'); + + % pressure = 60 Torr, T = 770 K + p = 60.0 * 133.3; + t = 770.0; + + gas.TPX = {300, p, '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. + gas.TP = {t, 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. + cstr.setInitialVolume(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; + w.install(cstr, env); + w.area = 1.0; + w.setHeatTransferCoeff(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 / gas.P) * ( gas.T / 273.15)); % m^3/s + mdot = gas.D * vdot; % kg/s + mfc = MassFlowController; + mfc.install(upstream, cstr); + mfc.setMassFlowRate(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; + v.install(cstr, downstream); + v.setValveCoeff(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; + network.advance(tme); + tm(n) = tme; + y(1,n) = cstr.massFraction('H2'); + y(2,n) = cstr.massFraction('O2'); + y(3,n) = cstr.massFraction('H2O'); + end + clf + figure(1) + plot(tm,y) + legend('H2','O2','H2O') + title('Mass Fractions') +end diff --git a/samples/matlab_experimental/prandtl1.m b/samples/matlab_experimental/prandtl1.m new file mode 100644 index 000000000..8423b081e --- /dev/null +++ b/samples/matlab_experimental/prandtl1.m @@ -0,0 +1,76 @@ +function prandtl1(g) +% PRANDTL1 Prandtl number for an equilibrium H/O gas mixture. +% +% This example computes and plots the Prandtl number for a +% hydrogen / oxygen mixture in chemical equilibrium for P = 1 +% atm and a range of temperatures and elemental O/(O+H) ratios. + + help prandtl1 + + if nargin == 1 + gas = g; + else + gas = GRI30('Mix'); + end + + pr = zeros(31,31); + xh2 = zeros(31,31); + visc = zeros(31,31); + lambda = zeros(31,31); + t = []; + xo2 = []; + io2 = gas.speciesIndex('O2'); + ih2 = gas.speciesIndex('H2'); + + minT = gas.minTemp; + maxT = gas.maxTemp; + dT = (maxT - minT)/30.0; + + t0 = cputime; + for i = 1:31 + t(i) = minT + dT*(i-1); + for j = 1:31 + xo2(j) = 0.99*(j-1)/30.0; + x = zeros(gas.nSpecies, 1); + x(io2) = xo2(j); + x(ih2) = 1.0 - xo2(j); + gas.TPX = {t(i), oneatm, x}; + equilibrate(gas, 'TP'); + visc(i,j) = gas.viscosity; + lambda(i,j) = gas.thermalConductivity; + gas.basis = 'mass'; + pr(i,j) = visc(i,j) * gas.cp / lambda(i,j); + x = gas.X; + xh2(i,j) = x(ih2); + end + end + disp(['CPU time = ' num2str(cputime - t0)]); + + % plot results + + clf; + %figure(1); + subplot(2,2,1); + surf(xo2,t,pr); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('Prandtl Number'); + + subplot(2,2,2); + surf(xo2,t,xh2); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('H_2 Mole Fraction'); + + subplot(2,2,3); + surf(xo2,t,visc); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('Viscosity'); + + subplot(2,2,4); + surf(xo2,t,lambda); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('Thermal Conductivity'); +end diff --git a/samples/matlab_experimental/prandtl2.m b/samples/matlab_experimental/prandtl2.m new file mode 100644 index 000000000..9ca49bbc7 --- /dev/null +++ b/samples/matlab_experimental/prandtl2.m @@ -0,0 +1,74 @@ +function prandtl2(g) + % PRANDTL2 Prandtl number for an equilibrium H/O gas mixture. + % + % This example does the same thing as prandtl1, but using + % the multicomponent expression for the thermal conductivity. + % + help prandtl2 + + if nargin == 1 + gas = g; + else + gas = GRI30('Multi'); + end + + pr = zeros(31,31); + xh2 = zeros(31,31); + visc = zeros(31,31); + lambda = zeros(31,31); + t = []; + xo2 = []; + io2 = gas.speciesIndex('O2'); + ih2 = gas.speciesIndex('H2'); + + minT = gas.minTemp; + maxT = gas.maxTemp; + dT = (maxT - minT)/30.0; + + t0 = cputime; + for i = 1:31 + t(i) = minT + dT*(i-1); + for j = 1:31 + xo2(j) = 0.99*(j-1)/30.0; + x = zeros(gas.nSpecies,1); + x(io2) = xo2(j); + x(ih2) = 1.0 - xo2(j); + gas.TPX = {t(i), oneatm, x}; + equilibrate(gas, 'TP'); + visc(i,j) = gas.viscosity; + lambda(i,j) = gas.thermalConductivity; + gas.basis = 'mass'; + pr(i,j) = visc(i,j) * gas.cp / lambda(i,j); + x = gas.X; + xh2(i,j) = x(ih2); + end + end + disp(['CPU time = ' num2str(cputime - t0)]); + + % plot results + + clf; + subplot(2,2,1); + surf(xo2,t,pr); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('Prandtl Number'); + + subplot(2,2,2); + surf(xo2,t,xh2); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('H_2 Mole Fraction'); + + subplot(2,2,3); + surf(xo2,t,visc); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('Viscosity'); + + subplot(2,2,4); + surf(xo2,t,lambda); + xlabel('Elemental O/(O+H)'); + ylabel('Temperature (K)'); + zlabel('Thermal Conductivity'); +end diff --git a/samples/matlab_experimental/rankine.m b/samples/matlab_experimental/rankine.m new file mode 100644 index 000000000..005ff67b1 --- /dev/null +++ b/samples/matlab_experimental/rankine.m @@ -0,0 +1,63 @@ +function [work, efficiency] = rankine(t1, p2, eta_pump, eta_turbine) + % This example computes the efficiency of a simple vapor power cycle. + help rankine + + % create an object representing water + w = Water; + + % start with saturated liquid water at t1 + w.setState_Tsat(t1, 1.0); + p1 = w.P; + + % pump it to p2 + basis = 'mass'; + pump_work = pump(w, p2, eta_pump); + h2 = w.H; + p2 = w.P; + + % heat to saturated vapor + w.setState_Psat(p2, 1.0); + h3 = w.H; + + heat_added = h3 - h2; + + % expand adiabatically back to the initial pressure + work = expand(w, p1, eta_turbine); + + % compute the efficiency + efficiency = (work - pump_work)/heat_added; +end + + +function w = pump(fluid, pfinal, eta) + % PUMP - Adiabatically pump a fluid to pressure pfinal, using a pump + % with isentropic efficiency eta. + + fluid.basis = 'mass'; + h0 = fluid.H; + s0 = fluid.S; + fluid.SP = {s0, pfinal}; + h1s = fluid.H; + isentropic_work = h1s - h0; + actual_work = isentropic_work / eta; + h1 = h0 + actual_work; + fluid.HP = {h1, pfinal}; + w = actual_work; +end + + +function w = expand(fluid, pfinal, eta) + % EXPAND - Adiabatically expand a fluid to pressure pfinal, using a + % turbine with isentropic efficiency eta. + + fluid.basis = 'mass'; + h0 = fluid.H; + s0 = fluid.S; + fluid.SP = {s0, pfinal}; + h1s = fluid.H; + isentropic_work = h0 - h1s; + actual_work = isentropic_work * eta; + h1 = h0 - actual_work; + fluid.HP = {h1, pfinal}; + w = actual_work; +end diff --git a/samples/matlab_experimental/reactor1.m b/samples/matlab_experimental/reactor1.m new file mode 100644 index 000000000..0f751e5fc --- /dev/null +++ b/samples/matlab_experimental/reactor1.m @@ -0,0 +1,85 @@ +function reactor1(g) +% REACTOR1 Zero-dimensional kinetics: adiabatc, constant pressure. +% +% This example illustrates how to use class 'Reactor' for zero-dimensional +% kinetics simulations. Here the parameters are set so that the reactor is +% adiabatic and very close to constant pressure. + help reactor1 + + if nargin == 1 + gas = g; + else + gas = GRI30('None'); + end + + P = 101325.0; + % set the initial conditions + gas.T = 1001.0; + gas.P = P; + nsp = gas.nSpecies; + xx = zeros(nsp, 1); + xx(1) = 0.285; + xx(4) = 0.142; + xx(48) = 0.573; + gas.X = xx; + + % create a reactor, and insert the gas + r = IdealGasReactor(gas); + + % create a reservoir to represent the environment + a = Solution('air.yaml','air','None'); + a.P = P; + env = Reservoir(a); + + % Define a wall between the reactor and the environment and + % make it flexible, so that the pressure in the reactor is held + % at the environment pressure. + w = Wall; + w.install(r,env); + + % set expansion parameter. dV/dt = KA(P_1 - P_2) + w.setExpansionRateCoeff(1.0e6); + + % set wall area + w.area = 1.0; + + % create a reactor network and insert the reactor: + network = ReactorNet({r}); + + nSteps = 100; + tim(nSteps) = 0; + temp(nSteps) = 0; + x(nSteps,3) = 0; + t = 0.0; + dt = 1.0e-5; + t0 = cputime; + for n = 1:nSteps + t = t + dt; + network.advance(t); + tim(n) = network.time; + temp(n) = r.T; + x(n,1:3) = gas.moleFraction({'OH','H','H2'}); + end + disp(['CPU time = ' num2str(cputime - t0)]); + + clf; + subplot(2,2,1); + plot(tim,temp); + xlabel('Time (s)'); + ylabel('Temperature (K)'); + subplot(2,2,2) + plot(tim,x(:,1)); + xlabel('Time (s)'); + ylabel('OH Mole Fraction (K)'); + subplot(2,2,3) + plot(tim,x(:,2)); + xlabel('Time (s)'); + ylabel('H Mole Fraction (K)'); + subplot(2,2,4) + plot(tim,x(:,3)); + xlabel('Time (s)'); + ylabel('H2 Mole Fraction (K)'); + clear all + cleanup + +end diff --git a/samples/matlab_experimental/reactor2.m b/samples/matlab_experimental/reactor2.m new file mode 100644 index 000000000..7553b796c --- /dev/null +++ b/samples/matlab_experimental/reactor2.m @@ -0,0 +1,60 @@ +function reactor2(g) +% REACTOR2 Zero-dimensional kinetics: adiabatic, constant volume. +% +% This example illustrates how to use class 'Reactor' for +% zero-dimensional kinetics simulations. Here the parameters are +% set so that the reactor is adiabatic and constant volume. + + help reactor2 + + if nargin == 1 + gas = g; + else + gas = GRI30('None'); + end + + % set the initial conditions + gas.TPX = {1001.0, oneatm, 'H2:2,O2:1,N2:4'}; + + % create a reactor, and insert the gas + r = IdealGasReactor(gas); + + % create a reactor network and insert the reactor + network = ReactorNet({r}); + + nSteps = 100; + tim(nSteps) = 0; + temp(nSteps) = 0; + x(nSteps,3) = 0; + t = 0; + dt = 1.0e-5; + t0 = cputime; + for n = 1:100 + t = t + dt; + network.advance(t); + tim(n) = network.time; + temp(n) = r.T; + x(n,1:3) = gas.moleFraction({'OH','H','H2'}); + end + disp(['CPU time = ' num2str(cputime - t0)]); + + clf; + subplot(2,2,1); + plot(tim,temp); + xlabel('Time (s)'); + ylabel('Temperature (K)'); + subplot(2,2,2) + plot(tim,x(:,1)); + xlabel('Time (s)'); + ylabel('OH Mole Fraction (K)'); + subplot(2,2,3) + plot(tim,x(:,2)); + xlabel('Time (s)'); + ylabel('H Mole Fraction (K)'); + subplot(2,2,4) + plot(tim,x(:,3)); + xlabel('Time (s)'); + ylabel('H2 Mole Fraction (K)'); + clear all + cleanup +end diff --git a/samples/matlab_experimental/reactor_ode.m b/samples/matlab_experimental/reactor_ode.m new file mode 100644 index 000000000..48ee5e729 --- /dev/null +++ b/samples/matlab_experimental/reactor_ode.m @@ -0,0 +1,51 @@ +function dydt = reactor_ode(t, y, gas, vdot, area, heatflux) +% REACTOR ODE system for a generic zero-dimensional reactor. +% +% Function REACTOR evaluates the system of ordinary differential +% equations for a zero-dimensional reactor with arbitrary heat +% transfer and volume change. +% +% Solution vector components: +% y(1) Total internal energy U +% y(2) Volume V +% y(3) Mass of species 1 +% .... +% y(2+nsp) Mass of last species +% + + [m,n] = size(y); + dydt = zeros(m,n); + + for j = 1:n + this_y = y(:,j); + int_energy = this_y(1); + vol = this_y(2); + masses = this_y(3:end); + + % evaluate the total mass, and the specific internal energy and volume. + total_mass = sum(masses); + u_mass = int_energy/total_mass; + v_mass = vol/total_mass; + + % set the state of the gas by specifying (u,v,{Y_k}) + gas.Y = masses; + gas.UV = {u_mass v_mass}; + p = gas.P; + + % volume equation + vdt = feval(vdot, t, vol, gas); + + % energy equation + a = feval(area, t, vol); + q = feval(heatflux, t, gas); + udt = -p * vdt + a * q; + + % species equations + ydt = total_mass * gas.ydot; + + % set up column vector for dydt + dydt(:,j) = [udt + vdt + ydt' ]; + end +end diff --git a/samples/matlab_experimental/surfreactor.m b/samples/matlab_experimental/surfreactor.m new file mode 100644 index 000000000..87c71340c --- /dev/null +++ b/samples/matlab_experimental/surfreactor.m @@ -0,0 +1,102 @@ +% SURFREACTOR Zero-dimensional reactor with surface chemistry +% +% This example illustrates how to use class 'Reactor' for +% zero-dimensional simulations including both homogeneous and +% heterogeneous chemistry. + +help surfreactor + +clear all +close all +cleanup +clc + +t = 870.0; +gas = Solution('ptcombust.yaml','gas'); + +% set the initial conditions +gas.TPX = {t, oneatm, '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.yaml','Pt_surf', gas); +surf.th.T = t; + +nsp = gas.nSpecies; +nSurfSp = surf.th.nSpecies; + +% create a reactor, and insert the gas +r = IdealGasReactor(gas); +r.setInitialVolume(1.0e-6) + +% create a reservoir to represent the environment +a = Solution('air.yaml','air','None'); +a.TP = {t, oneatm}; +env = Reservoir(a); + +% Define a wall between the reactor and the environment and +% make it flexible, so that the pressure in the reactor is held +% at the environment pressure. +w = Wall; +w.install(r, env); + +A = 1e-4; % Wall area + +% Add a reacting surface, with an area matching that of the wall +rsurf = ReactorSurface(surf, r, A); + +% set the wall area and heat transfer coefficient. +w.area = A; +w.setHeatTransferCoeff(1.0e1); % W/m2/K + +% set expansion rate parameter. dV/dt = KA(P_1 - P_2) +w.setExpansionRateCoeff(1.0); + +network = ReactorNet({r}); +% setTolerances(network, 1.0e-8, 1.0e-12); + +nSteps = 100; +p0 = r.P; +names = {'CH4','CO','CO2','H2O'}; +x = zeros([nSteps 4]); +tim = zeros(nSteps); +temp = zeros(nSteps); +pres = zeros(nSteps); +cov = zeros([nSteps nSurfSp]); +t = 0; +dt = 0.1; +t0 = cputime; +for n = 1:nSteps + t = t + dt; + network.advance(t); + tim(n) = t; + temp(n) = r.T; + pres(n) = r.P - p0; + cov(n,:) = surf.coverages'; + x(n,:) = gas.moleFraction(names); +end +disp(['CPU time = ' num2str(cputime - t0)]); + +clf; +subplot(2,2,1); +plot(tim,temp); +xlabel('Time (s)'); +ylabel('Temperature (K)'); +subplot(2,2,2); +plot(tim,pres); +axis([0 5 -0.1 0.1]); +xlabel('Time (s)'); +ylabel('Delta Pressure (Pa)'); +subplot(2,2,3); +semilogy(tim,cov); +xlabel('Time (s)'); +ylabel('Coverages'); +legend(speciesNames(surf)); +subplot(2,2,4); +plot(tim,x); +xlabel('Time (s)'); +ylabel('Mole Fractions'); +legend(names); +clear all +cleanup