From 8043b28315f4af8d1c5edc30eeae6e00d35b2840 Mon Sep 17 00:00:00 2001 From: Su Sun <56568705+ssun30@users.noreply.github.com> Date: Fri, 21 Jan 2022 16:55:08 -0500 Subject: [PATCH] Added 1D/CounterFLowDiffusionFlame --- .../1D/CounterFlowDiffusionFlame.m | 205 ++++++++++++++++++ 1 file changed, 205 insertions(+) create mode 100644 interfaces/matlab_experimental/1D/CounterFlowDiffusionFlame.m diff --git a/interfaces/matlab_experimental/1D/CounterFlowDiffusionFlame.m b/interfaces/matlab_experimental/1D/CounterFlowDiffusionFlame.m new file mode 100644 index 000000000..141b3ecac --- /dev/null +++ b/interfaces/matlab_experimental/1D/CounterFlowDiffusionFlame.m @@ -0,0 +1,205 @@ +function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidizer) + % Create a counter flow flame stack. + % :param left: + % Object representing the left inlet, which must be created using + % function 'Inlet'. + % :param flow: + % Object representing the flow, created with function + % 'AxisymmetricFlow'. + % :param right: + % Object representing the right inlet, which must be created using + % function 'Inlet'. + % :param tp_f: + % Object representing the fuel inlet gas, instance of class + % 'Solution', and an ideal gas. + % :param tp_o: + % Object representing the oxidizer inlet gas, instance of class + % 'Solution', and an ideal gas. + % :param oxidier: + % String representing the oxidizer species. Most commonly O2. + % :return: + % Instance of 'Stack' representing the left inlet, flow, and right + % inlet. + + %% Check input parameters + + if nargin ~= 6 + error('CounterFlowDiffusionFlame expects six input arguments.'); + end + if ~isIdealGas(tp_f) + error('Fuel gas object must represent an ideal gas mixture.'); + end + if ~isIdealGas(tp_o) + error('Oxidizer gas object must represent an ideal gas mixture.'); + end + if ~isInlet(left) + error('Left inlet object of wrong type.'); + end + if ~isFlow(flow) + error('Flow object of wrong type.'); + end + if ~isInlet(right) + error('Right inlet object of wrong type.'); + end + if ~ischar(oxidizer) + error('Oxidizer name must be of format character.'); + end + + %% Get the density of both fuel and oxidizer streams. + % To be used in determining velocity of each stream. Also get the + % temperature of both inlet streams. + + rhof = tp_f.D; + rho0 = tp_o.D; + tf = left.T; + tox = right.T; + + %% Find the species index of the oxidizer. + % To be used in determining initial strain rate. + + ioxidizer = tp_o.speciesIndex(oxidizer); + + %% Calculate the stoichiometric mixture fraction. + % Needed for determining location of flame edges and composition. + % elMoles function used to calculate the number of moles of C, H, and O + % atoms in the fuel and oxidizer streams: + % + % elMoles = elementalMassFraction/element atomic weight. + % + % From this, the stoichiometric Air/Fuel ratio can be determined. + % 1 Mole of O needs 2 Moles of C and 0.5 Moles of H for stoichiometric + % conditions. The stoichiometric mixture fraction, Zst, is then + % calculated. + + sFuel = elMoles(tp_f, 'O')- 2*elMoles(tp_f, 'C')- 0.5*elMoles(tp_f, 'H'); + sOx = elMoles(tp_o, 'O')- 2*elMoles(tp_o, 'C')- 0.5*elMoles(tp_o, 'H'); + phi = sFuel/sOx; + zst = 1.0/(1.0 - phi); + + %% Compute the stoichiometric mass fractions of each species. + % Use this to set the fuel gas object and calculate adiabatic flame + % temperature and equilibrium composition. + + spec = tp_f.speciesNames; % Get all of the species names in gas object. + nsp = tp_f.nSpecies; % Get total number of species in gas object. + % Get the current mass fractions of both fuel and inlet streams. + yox = tp_o.Y; + yf = tp_f.Y; + ystoich_double = zeros(1, nsp); % Create empty vector for stoich mass frac. + + for n = 1:nsp + % Calculate stoichiometric mass fractions. + ystoich_double(n) = zst*yf(n) + (1.0 - zst)*yox(n); + % Convert mass fraction vector to string vector. + ystoich_str{n} = num2str(ystoich_double(n)); + % Convert string vector to cell with SPECIES:MASS FRACTION format. + y_stoich{n} = [spec{n}, ':', ystoich_str{n}]; + end + + % Initialize stoichiometric mass fraction cell with first SP:Y value. + ystoich = [y_stoich{1}]; + for i = 2:nsp + % Update cell to have format similar to N2:Yst,O2:Yst,... + ystoich = [ystoich ',', y_stoich{i}]; + end + + % Set the fuel gas object as stoichiometric values and use equilibrate + % function to determine stoichiometric equilibrium temperature and mass + % fractions. + tp_f.TPY = {tf, tp_f.P, ystoich}; + tp_f.equilibrate('HP'); + teq = tp_f.T; + yeq = tp_f.Y; + + %% Estimate the strain rate. + % Based on the inlet stream velocities and determine initial 'guess' + % for mixture fraction based on mass flux ratio. + + zz = flow.gridPoints; + dz = zz(end) - zz(1); + mdotl = left.massFlux; + mdotr = right.massFlux; + uleft = mdotl/rhof; + uright = mdotr/rho0; + a = (abs(uleft) + abs(uright))/dz; + diff = tp_f.mixDiffCoeffs; + f = sqrt(a/(2.0*diff(ioxidizer))); + x0num = sqrt(uleft*mdotl)*dz; + x0den = sqrt(uleft*mdotr) + sqrt(uright*mdotr); + x0 = x0num/x0den; + + %% Calculate initial values of temperature and mass fractions. + % These values to be used for energy equation solution. Method is based + % on the Burke-Schumann model. + + nz = flow.nPoints; + zm = zeros(1, nz); + u = zeros(1, nz); + v = zeros(1, nz); + y = zeros(nz, nsp); + t = zeros(1, nz); + + for j = 1:nz + x = zz(j); + zeta = f*(x - x0); + zmix = 0.5*(1.0 - erf(zeta)); % Mixture fraction in flame. + zm(j) = zmix; + u(j) = a*(x0 - zz(j)); % Axial velocity. + v(j) = a; % Radial velocity. + if zmix > zst + for n = 1:nsp + y(j,n) = yeq(n) + (zmix - zst)*(yf(n) - yeq(n))/(1.0 - zst); + end + t(j) = teq + (tf - teq)*(zmix - zst)/(1.0 - zst); + else + for n = 1:nsp + y(j,n) = yox(n) + zmix*(yeq(n) - yox(n))/zst; + end + t(j) = tox + zmix*(teq - tox)/zst; + end + end + zrel = zz/dz; + + %% Create the flame stack. + % Set the profile of the flame with the estimated axial velocities, + % radial velocities, temperature, and mass fractions calculated above. + flame = Stack([left flow right]); + flame.setProfile(2, {'u', 'V'}, [zrel; u; v]); + flame.setProfile(2, 'T', [zrel; t] ); + for n = 1:nsp + nm = tp_f.speciesName(n); + flame.setProfile(2, nm, [zrel; transpose(y(:,n))]) + end +end + +%% Define elMoles function + +function moles = elMoles(tp, element) + % Determine the elemental moles in a gas object per unit mass. + + % Check input parameters + if nargin ~= 2 + error('elMoles expects two input arguments.'); + end + if ~isIdealGas(tp) + error('Gas object must represent an ideal gas mixture.'); + end + if ~ischar(element) + error('Element name must be of format character.'); + end + + % Calculate the moles per mass of mixture of an element within a gas + % object. The equation used is: + % + % elMoles = elMassFrac/Mel + % + % where elMassFrac is the elemental mass fraction within the gas object + % using the elementalMassFraction function; Mel is the atomic mass of + % the element. + + elMassFrac = tp.elementalMassFraction(element); + eli = tp.elementIndex(element); + M = tp.atomicMasses; + Mel = M(eli); + moles = elMassFrac/Mel; +end