function [Tout, Pout, rhoOut, gammaOut, VsoundOut, Mout, noEvapOut, Urs, X, phi, aeroLoadFrac, dropVolFrac, absCrossR5] = aerofrosh(varargin) % % When using this code for research, please cite: % M. F. Campbell, D. R. Haylett, D. F. Davidson, and R. K. Hanson. % AEROFROSH: A shock condition calculator for multi-component-fuel aerosol-laden flows. % Shock Waves 26:4 (2016) 429-447. % % "AEROsol FROsen SHock solver" % PURPOSE: % * Calculate conditions behind incident and reflected shocks in the % aerosol shock tube % INPUTS: % * Uis = incident shock velocity at endwall (m/s) % * T1 = temperature in region 1 (K) % * P1 = pressure in region 1 (torr) % * absorbR2 = absorbance (-ln(I/Io) in region 2 @ wavelength of absCross % * absorbR5 = absorbance (-ln(I/Io) in region 5 @ wavelength of absCross - NOT REQUIRED % * bathGasSpec = cell array of strings of gases in bath gas % e.g.: {'AR','N2','O2'} % * XbathGas = vector of mole fractions of gases in bath gas % e.g.: [0.01,0.79,0.20] % * fuelSpec = cell array of strings of liquids in fuel mixture % e.g.: {'C11H22O2','C12H26'} % * Xfuel = vector of mole fractions of liquids in fuel mixture % e.g.: [0.995, 0.005] % * vpfVals = vector containing vapor pressures of fuel components @ T1 (Torr) % * rfVals = vector containing liquid-phase densities of fuel components @ T1,P1 (g/cm3) % * hfVals = vector containing liquid-phase enthalpies of fuel components @ T1,P1 (J/kg) % * absCrossPolyMatrix = Matrix containing 4th order polynomial of % absorption cross sections for fuels. Each row in the matrix % represents one component of the fuel mixture. % Units produced by evaluation of polynomial should be m^2/mol. % e.g. [x1*T2^3 + x2*T2^2 + x3*T2 + x4 --> C11H22O2 % x1*T2^3 + x2*T2^2 + x3*T2 + x4]; --> C12H26 % * nasaFile = string name of *.dat file containing NASA polynomial % coefficients for molecules in shock tube gas mixture. % Should be in the same folder as "aerofrosh.m" function % e.g.: 'sandia.dat' % * Ltube = optical path length through shock tube (shock tube width) (meters) % e.g.: 0.1000 % OUTPUTS: % * Tout = [T2 T5] % * T2 = temperature in region 2 (K) % * T5 = temperature in region 5 (K) % * Pout = [P2 P5] % * P2 = pressure in region 2 (Pa) % * P5 = pressure in region 5 (Pa) % * rhoOut = [rho1 rho2 rho5] % * rho1 = density in region 1 (kg/m3) % * rho2 = density in region 2 (kg/m3) % * rho5 = density in region 5 (kg/m3) % * gammaOut = [gamma1 gamma2 gamma5] % * gamma1 = cp1/cv1 in region 1 % * gamma2 = cp2/cv2 in region 2 % * gamma5 = cp5/cv5 in region 5 % * VsoundOut = [Vsound1 Vsound2 Vsound5] % * Vsound1 = speed of sound in region 1 (m/s) % * Vsound2 = speed of sound in region 2 (m/s) % * Vsound5 = speed of sound in region 5 (m/s) % * Mout = [Mis Mrs] % * Mis = Mach number of incident shock % * Mrs = Mach number of reflected shock % * noEvapOut = [P2noEvap T2noEvap rho2noEvap] % * P2noEvap = pressure without droplet evap in region 2 (Pa) % * T2noEvap = temperaure without droplet evap in region 2 (K) % * rho2noEvap = density without droplet evap in region 2 (kg/m3) % * Urs = velocity of reflected shock at endwall (m/s) % * X = fully evaporated gas mole fractions % * phi = equivalence ratio of mixture % * aeroLoadFrac = molar percent of fuel loading from aerosol droplets % * dropVolFrac = fraction of volume in region 1 occupied by droplets % * absCrossR5 = absorption ((Io-I)/Io) cross section of fuel mixture % (m^2/mol) at T5, assuming complete evaporation in Region 2. % To verify complete evaporation in Region 2, a visible laser should % be directed through the tube. Its Mie Scattering will go to zero % as the droplets evaporate. % DEVELOPMENT HISTORY: % * Originally written by Dr. David Davidson using Mie Scattering in % Region 1 for fuel loading measurement. % * Updated by Daniel Haylett to use laser absorption in Region 2 % for fuel loading measurement. % * Improved and transcribed into MATLAB code by Matthew Campbell. % VERSION NUMBER: % * 0.5: 11/01/2010 - Demo release % * 1.0: 05/01/2011 - Initial release; Mis accounts for droplets % * 1.1: 05/28/2012 - Check for convergence % * 1.2: 05/29/2012 - Changed criteria for viable shock from Mrs to Urs % * 1.3: 06/02/2012 - Better SANDIA.DAT file reading routine % * 1.4: 05/11/2016 - Fix phi calculation. Change only affects oxygenated fuels. %% PARSE INPUT % where are we? fprintf('Parsing Input\n') % Grab information from "varargin" for i=1:2:length(varargin) switch varargin{i} case 'Uis' Uis = varargin{i+1}; % m/s case 'P1' P1 = varargin{i+1}; % torr case 'T1' T1 = varargin{i+1}; % K case 'absorbR2' absorbR2 = varargin{i+1}; % -ln(I/Io) in Region 2 (just before reflected shock) case 'absorbR5' absorbR5 = varargin{i+1}; % -ln(I/Io) in Region 5 (just after reflected shock) case 'bathGasSpec' bathGasSpec = varargin{i+1}; % cell array of species in bath gas case 'XbathGas' XbathGas = varargin{i+1}; % vector of bath gas mole fractions case 'fuelSpec' fuelSpec = varargin{i+1}; % cell array of species in fuel mixture case 'Xfuel' Xfuel = varargin{i+1}; % vector of fuel component mole fractions case 'vpfVals' vpfVals = varargin{i+1}; % vector of fuel component vapor pressures @ T1 (Torr) case 'rfVals' rfVals = varargin{i+1}; % vector of fuel component liquid-phase densities @ T1,P1 (g/cm3) case 'hfVals' hfVals = varargin{i+1}; % vector of fuel component liquid-phase enthalpies @ T1,P1 (J/kg) case 'absCrossPolyMatrix' absCrossPolyMatrix = varargin{i+1}; % matrix of absorption cross section polynomial; gives units of m^2/mol case 'nasaFile' nasaFile = varargin{i+1}; % string name of *.dat file containing NASA polynomial coefficients case 'Ltube' Ltube = varargin{i+1}; % optical path length through shock tube in meters otherwise % do nothing end end %% CHECK INPUTS % where are we? fprintf('Checking Input\n') % if vector is column, invert it so that it is a row vector. if size(bathGasSpec,1)>size(bathGasSpec,2) bathGasSpec=bathGasSpec'; end if size(XbathGas,1)>size(XbathGas,2) XbathGas=XbathGas'; end if size(fuelSpec,1)>size(fuelSpec,2) fuelSpec=fuelSpec'; end if size(Xfuel,1)>size(Xfuel,2) Xfuel=Xfuel'; end if size(vpfVals,1)>size(vpfVals,2) vpfVals=vpfVals'; end if size(rfVals,1)>size(rfVals,2) rfVals=rfVals'; end if size(hfVals,1)>size(hfVals,2) hfVals=hfVals'; end % check absorption cross section matrix if size(absCrossPolyMatrix,2) ~= 4 error('Absorption cross section polynomial matrix should have 4 columns and each row should represent one fuel component'); end % check gas vectors if (size(bathGasSpec,1) ~= size(XbathGas,1)) || (size(bathGasSpec,2) ~= size(XbathGas,2)) error('Bath gas species and mole fraction vectors must have same size') end % check liquid vectors if (size(fuelSpec,1) ~= size(Xfuel,1)) || (size(fuelSpec,2) ~= size(Xfuel,2)) error('Fuel species and mole fraction vectors must have same size') end % check temperature if T1<273.15 error('Temperature should be in Kelvin, not Celsius!'); end % ensure that mole fractions sum to 1.0 if sum(XbathGas)-1.0 > 0.000001 error('Bath gas mole fractions do not sum to 1.0'); end if sum(Xfuel)-1.0 > 0.000001 error('Fuel component mole fractions do not sum to 1.0'); end %% INITIAL COMPUTATIONS % where are we? fprintf('Performing Initial Computations\n') % create a cell array of all the species names (bath gases, then liquids) speciesVec = [bathGasSpec fuelSpec]; % make it UPPERCASE speciesVec = upper(speciesVec); % remember how many components there are in the fuel mixture numFuelComps = size(fuelSpec,2); % Transform pressure units P1torr = P1; % Torr P1 = P1torr*(1/760)*(101325/1); % Torr --> Pa % The velocity of the test gas relative to the shock in region 1 (u1) is % equal to the incident shock velocity (Uis) u1 = Uis; % m/s %% CONSTANTS % where are we? fprintf('Defining Constants\n') % elemental molecular weights eleMolWts = { % g/mol (cell array) 'AC' 227 'AG' 107.8682 'AL' 26.9815386 'AM' 243 'AR' 39.948 'AS' 74.92160 'AT' 210 'AU' 196.966569 'B' 10.811 'BA' 137.327 'BE' 9.012182 'BH' 262 'BI' 208.98040 'BK' 247 'BR' 79.904 'C' 12.0107 'CA' 40.078 'CD' 112.411 'CE' 140.116 'CF' 251 'CL' 35.453 'CM' 247 'CO' 58.933195 'CR' 51.9961 'CS' 132.9054519 'CU' 63.546 'DB' 262 'DY' 162.500 'ER' 167.259 'ES' 252 'EU' 151.964 'F' 18.9984032 'FE' 55.845 'FM' 257 'FR' 223 'GA' 69.723 'GD' 157.25 'GE' 72.64 'H' 1.00794 'HE' 4.002602 'HF' 178.49 'HG' 200.59 'HO' 164.93032 'HS' 265 'I' 126.90447 'IN' 114.818 'IR' 192.217 'K' 39.0983 'KR' 83.798 'LA' 138.90547 'LI' 6.941 'LR' 262 'LU' 174.9668 'MD' 258 'MG' 24.305 'MN' 54.938045 'MO' 95.96 'MT' 266 'N' 14.0067 'NA' 22.98976928 'NB' 92.90638 'ND' 144.242 'NE' 20.1797 'NI' 58.6934 'NO' 259 'NP' 237 'O' 15.9994 'OS' 190.23 'P' 30.973762 'PA' 231.03588 'PB' 207.2 'PD' 106.42 'PM' 145 'PO' 209 'PR' 140.90765 'PT' 195.084 'PU' 244 'RA' 226 'RB' 85.4678 'RE' 186.207 'RF' 261 'RH' 102.90550 'RN' 222 'RU' 101.07 'S' 32.065 'SB' 121.760 'SC' 44.955912 'SE' 78.96 'SG' 263 'SI' 28.0855 'SM' 150.36 'SN' 118.710 'SR' 87.62 'TA' 180.94788 'TB' 158.92535 'TC' 98 'TE' 127.60 'TH' 232.03806 'TI' 47.867 'TL' 204.3833 'TM' 168.93421 'U' 238.02891 'V' 50.9415 'W' 183.84 'XE' 131.293 'Y' 88.90585 'YB' 173.054 'ZN' 65.38 'ZR' 91.224}; % universal gas constant Ru = 8.314471; % J/mol*K %% READ SANDIA THERMODYNAMIC DATA FILE, COMPUTE SPECIES MOLECULAR WEIGHTS % where are we? fprintf('Reading Sandia Thermo Data\n') % Read file (see subfunction below) [coefVec, nameEle, numEle] = readSandia(speciesVec, nasaFile); % where are we? fprintf('Computing Species Molecular Weights\n') % Find molecular weight of each species, gas and liquid specMolWts = getSpecMolecWts(nameEle, numEle, eleMolWts); % g/mol %% DETERMINE MOLAR AND MASS AVERAGE PROPERTIES FOR FUEL MIXTURE % where are we? fprintf('Computing Average Properties Of Fuel Mixture\n') % Molecular Weight fuelMolWtAvg = specMolWts((end-numFuelComps+1):end)*Xfuel'; % (g/mol,i)*(mol,i/mol,total) = g/mol % Vapor Pressure % use Raoult's Law of ideal mixtures: P,a = x,a*Psat,a % compute the vapor pressure for each component in solution % then sum these partial pressures up to find the total vapor pressure % according to Dalton's law of partial pressures. vpfPartialPressure = vpfVals.*Xfuel; % Torr vpfAvg = sum(vpfPartialPressure); % Torr % Density % assume 1 mol of liquid fuel. the total mass of this 1 mol of liquid fuel % is the fuel's average molecular weight, and the mass of each component is % the mole fraction of that component multiplied by its molar mass. the % total volume is the mass of each component divided by its density. the % density is found by dividing total mass by total volume. totalMass = fuelMolWtAvg*1; % g/mol * 1 mol = g massEachComponent = specMolWts((end-numFuelComps+1):end).*Xfuel; % g/mol, note the ".*" operator for element-by-element multiplication. massFracs = massEachComponent/totalMass; % g/g rfAvg = 1/(sum(massFracs./rfVals)); % g/mol % Specific Enthalpy % assume 1 mol of liquid fuel. find the mass fraction of each component. % then do a mass fraction weighted summation to find the specific enthalpy. hfAvg = hfVals*massFracs'; % (J/kg)*(g/g) = J/kg %% DETERMINE GAS PHASE MOLE FRACTIONS IN REGION ONE % region 1: bath gas + liquid vapor pressure; % region 2: bath gas + liquid vapor pressure + evaporated droplets; % region 5: assume to be same as region 2 (assume no fuel decomposition in region 2); % where are we? fprintf('Determining Region 1 Gas Mole Fracs\n') % assemble mole fractions, accounting for vapor pressure of fuel mole fraction of gas species % the measured P1 includes the vapor pressure; therefore the partial pressures of bath gas % species are less than simply their bath gas mole fraction times total pressure. X1onlyGas = XbathGas*(P1torr-vpfAvg)/P1torr; % add in fuel component mole fractions (derived from vapor pressure) % note that X1onlyGas is the mole fractions IN THE GAS PHASE, NOT the % overall mole fractions in region 1 (which account for both liquid and gas) X1fuelMoleFracs = vpfPartialPressure/P1torr; X1onlyGas = [X1onlyGas X1fuelMoleFracs]; % the ratio of vapor-phase fuel components due to vapor pressure is NOT the same as the % molar fuel ratio in the liquid. This composition is computed below, though the program % uses another method to solve the mass balance equations. fuelFracsVapor = X1fuelMoleFracs/sum(X1fuelMoleFracs); % for reference only %% DETERMINE GAS-ONLY PROPERTIES IN REGION ONE % Weighted sum by the mole fractions in the gas; this is the region 1 value % multiply then add mole fractions by species molec weights gasMolWtNoDrops = X1onlyGas*specMolWts'; % g/mol % compute the enthalpy of gas without drops [CpNoDrops, hNoDrops] = getCph_T(T1, X1onlyGas, gasMolWtNoDrops, Ru, coefVec, speciesVec); % J/kg*K, J/kg % compute ideal gas constant for gas-part of mixture in region 1 Rgas1 = (Ru/gasMolWtNoDrops)*1000; % (J/mol*K)*(mol/g)*(1000 g/1 kg) = J/kg*K % compute gas specific volume in region 1 v1noDrops = Rgas1*T1/P1; % (J/kg*K)*K/Pa = m3/kg %% SOLVE FOR INCIDENT SHOCK CONDITIONS % At this point we have a vector of bath gas with the vapor pressure of the % fuel mixed in. However, we also need to account for the change in mole % fraction which occurs when the aerosolus fuel evaporates. This loop uses % the absorption in Region 2 to compute the new mole fractions. It also % contains an inner loop which uses Newton-Raphson to solve for T2 and P2. % Cm is the mole fraction of moles of molecules in the liquid phase, % ie, in the droplets. Cm = mol_droplets / (mol_droplets + mol_gas). % where are we? fprintf('Solving for Incident Shock Conditions\n') % First guess is that the aerosol added no fuel, ie just vapor pressure % ***In general, Xguess is the mole fractions after complete evaporation*** Xguess = X1onlyGas; % Set Cm initially to be zero, ie there were no moles in liquid phase Cm = 0; % Iteration parameters for Mole Fraction loop maxTimeMF = 10; % second maxIterMF = 1000; % iterations tolMF = 1e-6; % convergence tolerance stopperMF = false; % Iteration loop iterMF = 0; startTimeMF = tic; while stopperMF == false; % inc the counter iterMF = iterMF+1; % CONDITIONS FOR MIXTURE % compute the molecular weight of the mixture once it is fully evaporated gasMolWt = Xguess*specMolWts'; % g/mol % gas constant for fully evaporated mixture % this is defined in regions 2 and 5 where there is only gas. % it is undefined for the total mixture in region 1 where there are both gas and droplets. % however, it IS fully defined for the gas-portion of the region 1 mixture (computed above). Rgas2 = (Ru/gasMolWt)*1000; % (J/mol*K)*(mol/g)*(1000 g/1 kg) = J/kg*K % Assume 1 mole of liquid and gas combined in the pre-evaporation mixture in region 1. % Then the mass and volume of each can be determined. massLiq = Cm*1*fuelMolWtAvg*(1/1000); % (mol frac)*(1 mol)*(g/mol)*(1 kg/1000 g) = kg massGas = (1-Cm)*1*gasMolWtNoDrops*(1/1000); % kg volLiq = Cm*1*fuelMolWtAvg*(1/rfAvg)*(1/1e6); % (mol frac)*(1 mol)*(g/mol)*(cm3/g)*(1 m3/1e6 cm3) = m3 volGas = (1-Cm)*1*Ru*T1/P1; % mol*(J/mol*K)*K/Pa = m3 % compute density and specific volume of total mixture in region 1 rho1 = (massLiq+massGas)/(volLiq+volGas); % kg/m3 v1 = 1/rho1; % m3/kg % Adjust the enthalpy for the gas in region 1 by subtracting out the enthalpy % of vaporization and adding in the enthalpy that the evaporated fuel would have at T1. % This is equivalent to h_liq: [h_liq = h_vap - h_fg] % Add in the enthalpy via mass fractions % Recall that "gasMolWt" is the molecular weight of the gas including evaporated liquid massFracLiq = (Cm*fuelMolWtAvg)/gasMolWt; % (molfrac*g/mol)/(molfrac*g/mol) = no units massFracGas = ((1-Cm)*gasMolWtNoDrops)/gasMolWt; if abs(massFracLiq+massFracGas-1) > 0.000001 error('Mass fractions not summing') end h1 = hNoDrops*massFracGas + hfAvg*massFracLiq; % J/kg % The sound speed required is the FROZEN sound speed, discussed in Temkin 2000. % It is derived from the gaseous sound speed. % specific heat at constant volume CvNoDrops = CpNoDrops-Rgas1; % J/kg*K % ratio of specific heats gamma1 = CpNoDrops/CvNoDrops; % density of gas only, without drops rho1noDrops =1/v1noDrops; % kg/m3 % speed of sound in gas only, without drops Vsound1gas = sqrt(gamma1*P1/rho1noDrops); % m/s % volume fraction of droplets (droplet volume concentration) dropVolFrac = volLiq/(volLiq+volGas); % droplet mass concentration dropMassConc = dropVolFrac*(rfAvg/rho1noDrops)/(1+(1-rho1noDrops/rfAvg)*dropVolFrac*(rfAvg/rho1noDrops)); % frozen sound speed (sound speed corrected for droplets) Vsound1 = Vsound1gas*(1-dropMassConc)^(1/2)/(1-dropVolFrac); % compute mach number of incident shock Mis = u1/Vsound1; % no units % SOLVE IDEAL SHOCK EQNS % Before Newton-Rhapson guess using ideal relations: % P2/P1 = ((2*gamma1*Mis^2-(gamma1-1))/(gamma1+1)) % T2/T1 = (gamma1*Mis^2-(gamma1-1)/2)*((gamma1-1)/2*Mis^2+1)/(((gamma1+1)/2)^2*Mis^2) % gamma = cp1/(cp1-1) --> assumes driven gas is mostly Ar with Cp ~ 2.5 % Mis = U1/Vsound1 % compute ideal pressure in region 2 P2ideal = P1*((2*gamma1*Mis^2-(gamma1-1))/(gamma1+1)); % compute ideal temperature in region 2 T2ideal = T1*(gamma1*Mis^2-(gamma1-1)/2)*((gamma1-1)/2*Mis^2+1)/(((gamma1+1)/2)^2*Mis^2); % K % specific volume v2ideal = Rgas2*T2ideal/P2ideal; % m3/kg % SOLVE NORMAL SHOCK EQNS % give it the ideal values to help pinpoint the solution [P2, T2, v2, Cp2, h2] = incidentShockNR(... P2ideal, T2ideal, v2ideal, Rgas2, ... P1, v1, h1, u1, ... Xguess, gasMolWt, coefVec, speciesVec, Ru); % ACCOUNT FOR OBSERVED ABSORPTION % compute absorption coefficient at T2 for all fuel components for i=1:numFuelComps absCrossPoly = absCrossPolyMatrix(i,:); % grab the coefficients in the row crossR2vals(i) = polyval(absCrossPoly,T2); % (m^2/mol) evaluate the polynomial at T2 end % compute the molar average absorption coefficient % for the first iteration, base this on the mole fractions in % the liquid fuel so if the vapor pressure is ~0 the algorithm % works. After that, base it on the vaporized mole fractions if iterMF == 1; fuelEvapMoleFracs = Xfuel; else tempTot = sum(Xguess((end-numFuelComps+1):end)); % sum of mole fractions of fuel fuelEvapMoleFracs = Xguess((end-numFuelComps+1):end)/tempTot; % divide to find fuel component mole fractions end crossR2avg = crossR2vals*fuelEvapMoleFracs'; % m^2/mol % compute molar concentration of fuel (this is for all evaporated fuel combined) nFuel = absorbR2/(crossR2avg*Ltube); % mol/m^3 =(no units)/(m^2/mol*m) % compute molar concentration of gas in shock tube at T2, P2 nTot2 = P2/(Ru*T2); % mol/m^3 % compute mole fraction of fuel in region 2 XevapFuelNew = nFuel/nTot2; % compute the old mole fraction of fuel total % this is the sum of the mole fractions of individual fuel components XevapFuelOld = sum(Xguess((end-numFuelComps+1):end)); % CHECK RESULTS AND DECIDE WHAT TO DO % check variable checkMF = abs((XevapFuelOld-XevapFuelNew)/XevapFuelNew); % decide what to do if checkMF < tolMF; stopperMF = true; % exit the loop elseif toc(startTimeMF) > maxTimeMF error('Timed out on MF loop') elseif iterMF > maxIterMF error('Itered out on NR loop') else % update Cm (mole fraction of molecules in liquid) % Assume 1 mol total of gaseous mixture prior to evaporation, % then realize that the total moles of gas mixture increase only % due to the addition of evaporated fuel droplets. XnonFuelPreEvap = sum(X1onlyGas(1:end-numFuelComps)); % = 1-sum(X1onlyGas(end-numFuelComps+1)) XnonFuelPostEvap = 1-XevapFuelNew; molAdded = (XnonFuelPreEvap-XnonFuelPostEvap)/(XnonFuelPostEvap); % Cm is the fraction of molecules (liquid and gas % combined) that are in the liquid phase in region 1 Cm = molAdded/(1+molAdded); % mol evaporated divided by total mols (1 mol gas mixture + mol evaporated) % solve for new post-evaporation mole fractions % Again assume 1 mol total of gaseous mixture prior to evaporation. % Find the moles of non-fuel bath gas components present. molesBathGasComponents = X1onlyGas(1:end-numFuelComps)*(1.0); % (mol/mol)*mol = mol % Find the moles of each fuel component prior to evaporation % in the gas phase due to their vapor pressures. molesFuelComponentsVapor = X1onlyGas((end-numFuelComps+1):end)*(1.0); % (mol/mol)*mol = mol % Find the moles of each fuel component added via aerosol evaporation % assume that the fuel nebulizes in proportion to its mole fractions % ie mole fractions in the aerosol drops are identical to the mole fractions in the liquid molesFuelComponentsAerosol = Xfuel*molAdded; % we now know how many moles of each species there are, and from which source. % add the vapor and aerosol fuel mole fractions molesFuelComponents = molesFuelComponentsVapor+molesFuelComponentsAerosol; % assemble new mole factions vector Xguess = [molesBathGasComponents molesFuelComponents]/(1+molAdded); % check to make sure mole fractions are consistent if abs(XevapFuelNew-sum(Xguess(end-numFuelComps+1:end))) > 0.000001 error('Mole fractions not summing') end end end % ensure that mole fraction of fuel did not decrease % note that minimum mole fraction of fuel is that of its vapor pressure if sum(Xguess((end-numFuelComps+1):end)) < sum(X1fuelMoleFracs) error('Xfuelabs(%4.3s) < Xfuelvp(%4.3s)',sum(Xguess(end-numFuelComps+1:end)),sum(X1fuelMoleFracs)) end % compute gas velocity in region 2 u2 = u1*(v2/v1); % m/s % compute gamma in region 2 Cv2 = Cp2 - Rgas2; % J/kg*K gamma2 = Cp2/Cv2; % compute density in region 2 rho2 = 1/v2; % kg/m3 % assign mole fractions X = Xguess; % compute speed of sound in region 2 Vsound2 = sqrt(gamma2*P2/rho2); % m/s % check to make sure shock conservations equations are solved tolCheck = 1e-6; eq1 = ((rho1*u1)-(rho2*u2))/(rho1*u1); % (kg/m3)*(m/s) eq2 = ((P1+rho1*u1^2)-(P2+rho2*u2^2))/(P1+rho1*u1^2); % Pa + (kg/m3)*(m2/s2)*(1 Pa / 1 kg/ms2) eq3 = ((h1+(1/2)*u1^2)-(h2+(1/2)*u2^2))/(h1+(1/2)*u1^2); % J/kg*(1 kgm2/s2 / 1 J) + m2/s2 if abs(eq1)>tolCheck error('Continuity not holding over incident shock') elseif abs(eq2)>tolCheck error('Momentum not holding over incident shock') elseif abs(eq3)>tolCheck error('Energy not holding over incident shock') end %% SOLVE INCIDENT SHOCK CONDITIONS WITHOUT EVAPORATION % where are we? fprintf('Solving for Incident Shock with No Evaporation\n') % It is often helpful to know the shock conditions in Region 2 prior to % evaporation in order to estimate the effect that the aerosol had on the % shock wave. None of the values computed in this section are necessary % for the computation of the reflected shock conditions. % Assume that the aerosol is not present and essentially compute the % properties using only the test gas mixture (including fuel vapor) % Initial guess: % * Use P2ideal, T2ideal, v2ideal, because these will be closer to % frozen shock values (see Guha 1992) % Gas properties: % * Use Rgas1, P1, v1noDrops, h1noDrops, X1onlyGas, gasMolWtNoDrops [P2noEvap, T2noEvap, v2noEvap, dummy1, dummy2] = incidentShockNR(... P2ideal, T2ideal, v2ideal, Rgas1, ... P1, v1noDrops, hNoDrops, u1, ... X1onlyGas, gasMolWtNoDrops, coefVec, speciesVec, Ru); rho2noEvap = 1/v2noEvap; %% SOLVE REFLECTED SHOCK CONDITIONS % where are we? fprintf('Solving for Reflected Shock Conditions\n') % Before Newton-Rhapson guess using ideal relations: % alpha=(gamma+1)/(gamma-1) % P5/P2=(alpha+2-P1/P2)/(1+alpha*P1/P2) % T5/T1=(P5/P2)*(alpha+P5/P2)/(1+alpha*P5/P2) % gamma=cp1/(cp1-1) --> assumes driven gas is mostly Ar with Cp ~ 2.5 % ideal gas constant in region 5 is the same as for region 2 because % molecular species and mole fractions are the same. Rgas5 = Rgas2; % compute alpha alpha = (gamma2+1)/(gamma2-1); % compute ideal pressure in region 2 P5ideal = P2*(alpha+2-P1/P2)/(1+alpha*P1/P2); % Pa % compute ideal temperature in region 2 T5ideal = T2*(P5ideal/P2)*(alpha+P5ideal/P2)/(1+alpha*P5ideal/P2); % K % specific volume v5ideal = Rgas5*T5ideal/P5ideal; % m3/kg % Solve the Reflected Shock Equations [P5, T5, v5, Cp5, h5] = reflectedShockNR(... P5ideal, T5ideal, v5ideal, Rgas5, ... P2, v2, h2, u2, u1, ... X, gasMolWt, coefVec, speciesVec, Ru); % compute density in region 5 rho5 = 1/v5; % kg/m3 % compute gas velocity in region 5 % given the assumption that the gas is stagnated after the reflected shock % in region 5, u5 (the velocity of the test gas in region 5 relative to the % reflected shock) is equal and opposite to the velocity of the reflected % shock (Urs). u5 = v5*(u1-u2)/(v2-v5); % m/s Urs = u5; % m/s % compute gamma in region 5 Cv5 = Cp5 - Rgas5; % J/kg*K gamma5 = Cp5/Cv5; % compute density in region 5 rho5 = 1/v5; % kg/m3 % compute speed of sound in region 5 Vsound5 = sqrt(gamma5*P5/rho5); % m/s % Mach number of reflected shock u2dash = u5+u1-u2; % m/s Mrs = u2dash/Vsound2; % no units % check answers % T5 will always be greater than T2, no mater how much evaporation % P5 will always be greater than P2 % Urs will always be less than Uis % however, Mrs may not be less than Mis; this is because as droplet % volume fraction increases, T2, R2 (often), and gamma2 decrease, % leading to a decrease in Vsound2 and an increase in Mrs. if (T5<=T2) || (P5<=P2) || (Urs>=Uis) error('Impossible reflected shock conditions') end % check to make sure shock conservations equations are solved tolCheck = 1e-6; eq1 = ((rho2*u2dash)-(rho5*u5))/(rho2*u2dash); % (kg/m3)*(m/s) eq2 = ((P2+rho2*u2dash^2)-(P5+rho5*u5^2))/(P2+rho2*u2dash^2); % Pa + (kg/m3)*(m2/s2)*(1 Pa / 1 kg/ms2) eq3 = ((h2+(1/2)*u2dash^2)-(h5+(1/2)*u5^2))/(h2+(1/2)*u2dash^2); % J/kg*(1 kgm2/s2 / 1 J) + m2/s2 if abs(eq1)>tolCheck error('Continuity not holding over reflected shock') elseif abs(eq2)>tolCheck error('Momentum not holding over reflected shock') elseif abs(eq3)>tolCheck error('Energy not holding over reflected shock') end %% COMPUTE ALF, DVF, PHI % where are we? fprintf('Solving for ALF, DVF, and Phi\n') % compute molar fraction of fuel loading due to aerosol % Cm = mole fraction just liquid fuel % sum(X(...)) = total mole fraction of fuel (gas + liquid) aeroLoadFrac = Cm/sum(X((end-numFuelComps+1):end)); % function to compute stoichiometry (see below) phi = getStoich(X, speciesVec, nameEle, numEle); %% COMPUTE ABSORPTION CROSS SECTION IN REGION FIVE % where are we? fprintf('Solving for Absorption Coefficient in Region Five\n') % decide whether to do this part or not if exist('absorbR5','var') % if the user gave an absorption value in R5 % Assume that we have complete evaporation in R2 (this % can be verified by seeing Mie scattering of a visible laser go to % zero in Region 2). Then, we can compute the absorption cross % section for T5 using this mole fraction. Note that the accuracy % of the R5 absorption cross section depends on the accuracy of the % R2 absorption cross section! % compute molar concentration in shock tube at T5, P5 nTotR5 = P5/(Ru*T5); % mol/m3 % Grab the fuel mole fraction in R5 computed via iteration XfuelR5iter = sum(X((end-numFuelComps+1):end)); % compute fuel particle density in R5 (mole frac * total) nFuelR5iter = XfuelR5iter*nTotR5; % mol/m^3 % compute absorption cross section % this is the molar average cross section for the fuel mixture absCrossR5 = absorbR5/(nFuelR5iter*Ltube); % m^2/mol; Ltube defined above else % output an obviously wrong answer absCrossR5 = -999; end %% ASSEMBLE OUTPUT VECTORS % where are we? fprintf('Assembling Output Vectors\n') % output vectors Tout = [T2 T5]; % K Pout = [P2 P5]; % Pa rhoOut = [rho1 rho2 rho5]; % kg/m3 gammaOut = [gamma1 gamma2 gamma5]; VsoundOut = [Vsound1 Vsound2 Vsound5]; % m/s Mout = [Mis Mrs]; % m/s noEvapOut = [P2noEvap T2noEvap rho2noEvap]; % Pa, K, kg/m3 function specMolWts = getSpecMolecWts(nameEle, numEle, eleMolWts) % This function computes the species molecular weights % Find molecular weight of each species (gas and liquid) (g/mol) % rows of nameEle represent species % columns of nameEle are the elements composing that species specMolWts = zeros(1,size(nameEle,1)); % row vector; each entry is mol wt of the species gotIt = false; % used to determine if an element was found for i=1:size(nameEle,1) % outer loop over all rows of nameEle for j=1:size(nameEle,2) % inner loop over all columns of nameEle tempEl = nameEle{i,j}; if isempty(tempEl) continue; % if nothing there, skip it. else for k=1:size(eleMolWts,1) %loop over rows of eleMolWts if strcmp(eleMolWts{k,1},tempEl) % add to the total weight; % the amount to add is the molecular weight of the % element times the number of that element present % in the species specMolWts(i) = specMolWts(i)+eleMolWts{k,2}*numEle(i,j); gotIt = true; % we found the element break; % stop looking, we got it; end end if gotIt == false error('Atomic weight of %1s not found',tempEl); end end gotIt = false; end end function [Cp, h] = getCph_T(T, X, gasMolWt, Ru, coefVec, speciesVec) % This finds the specific heat and enthalpy at the specified temperature % loop over all species to compute for each species first HoverRT = zeros(1,size(coefVec,3)); CpoverR = zeros(1,size(coefVec,3)); for i=1:size(coefVec,3) % loop over 3rd dimension of coefVec % grab temperature limits Tlow = coefVec(3,1,i); Tcom = coefVec(3,2,i); Thigh = coefVec(3,3,i); % check temperature limits if T < Tlow %fprintf('%3s: T (%1.2f K) < Tlow (%1.0f K)\n',speciesVec{i},T,Tlow); elseif T > Thigh %fprintf('%3s: T (%1.2f K) > Thigh (%1.0f K)\n',speciesVec{i},T,Thigh); end % decide which set of data to use if T < Tcom % low temp set j = 2; else % high temp set for T >= Tcom j = 1; end % grab temperature coefficients a1 = coefVec(j,1,i); a2 = coefVec(j,2,i); a3 = coefVec(j,3,i); a4 = coefVec(j,4,i); a5 = coefVec(j,5,i); a6 = coefVec(j,6,i); a7 = coefVec(j,7,i); % compute Cp/R CpoverR(i) = a1 + a2*T + a3*(T^2) + a4*(T^3) + a5*(T^4); % compute h/RT HoverRT(i) = a1 + a2*T/2 + a3*(T^2)/3 + a4*(T^3)/4 + a5*(T^4)/5 + a6/T; end % now compute mixture properties Cp = (X*CpoverR')*Ru*(1000/gasMolWt); % (no units)*(J/mol*K)*(1000 g/kg)*(mol/g) = J/kg*K h = (X*HoverRT')*T*Ru*(1000/gasMolWt); % J/kg function [P2, T2, v2, Cp2, h2] = incidentShockNR(P2ideal, T2ideal, v2ideal, Rgas2, P1, v1, h1, u1, Xguess, gasMolWt, coefVec, speciesVec, Ru) % This function solves for the conditions behind the incident shock % Solve the Incident Shock Equations % P*v = Rgas*T % rho1*u1 = rho2*u2 % P1+rho1*u1^2 = P2+rho2*u2^2 % h1+(1/2)*u1^2 = h2+(1/2)*u2^2 % by minimizing the following pair of equations (f1, f2) % f1 = (P2/P1-1)+(u1^2/(P1*v1))*(v2/v1-1) % f2 = (h2-h1)/((1/2)*u1^2)+(v2^2/v1^2-1) % use Newton-Raphson Method % Assign ideal to first guess % Note that we do not use the P2 and T2 from the last iteration % of the mole fractiion loop, because this prevents the new % Newton-Raphson loop from being ill-conditioned as we changed % the mole fractions. A conservative method of solving the % equations is to change the mole fractions, but start from the % ideal P2 and T2 each time. P2guess = P2ideal; T2guess = T2ideal; v2 = v2ideal; % NR iteration parameters maxTimeNR = 10; % seconds maxIterNR = 100; % iterations tolNR = 1e-6; % convergence tolerance stopperNR = false; % NR loop iterNR = 0; startTimeNR = tic; while stopperNR == false; % inc the counter iterNR = iterNR+1; % compute cp and h for gas in region 2 (see subfunction below) [Cp2, h2] = getCph_T(T2guess, Xguess, gasMolWt, Ru, coefVec, speciesVec); % J/kg*K, J/kg % compute f1 f1 = (P2guess/P1-1)+(u1^2/(P1*v1))*(v2/v1-1); % compute f2 f2 = (h2-h1)/((1/2)*u1^2)+((v2/v1)^2-1); % compute derivatives for jacobian % note that dh/dP = 0 (ideal gas) % note that dT/dP = 0 (independent variables) df1dP2 = (1/P1)+(u1^2/(P1*v1^2))*(-v2/P2guess); % df1/dP2 df1dT2 = (1/P1)*(0)+(u1^2/(P1*v1^2))*(v2/T2guess); df2dP2 = (2/u1^2)*(0)+(2*v2/v1^2)*(-v2/P2guess); df2dT2 = (2/u1^2)*(Cp2)+(2*v2/v1^2)*(v2/T2guess); % form Jacobian matrix [df1/dP2 df1/dT2; df2/dP2 df2/dT2] J = [df1dP2 df1dT2; ... df2dP2 df2dT2]; % matrix multiplication to solve for new T2 and P2 guesses % using matlab backslash operator rather than inverse saves time! % general method: [P2new T2new]' = [P2old T2old]' - inv(J)*[f1old f2old]' PTold = [P2guess T2guess]'; fvec = [f1 f2]'; PTnew = PTold - J\fvec; P2new = PTnew(1); T2new = PTnew(2); % check results - the T2 and P2 we guess should give the same T2,P2 checkP = abs((P2guess-P2new)/P2new); % percent change from new checkT = abs((T2guess-T2new)/T2new); % decide whether to keep going or quit if (checkP < tolNR) && (checkT < tolNR) stopperNR = true; % exit the loop elseif toc(startTimeNR) > maxTimeNR; error('Timed out on NR loop') elseif iterNR > maxIterNR error('Itered out on NR loop') else % reassign the variables for the next iteration P2guess = P2new; T2guess = T2new; % recompute specific volume v2 = Rgas2*T2guess/P2guess; % m3/kg % prevent the algorithm from going the wrong way while P2guess < 0 P2guess = mean([P2ideal P2guess]); % average back to ok range end while T2guess < 0 T2guess = mean([T2ideal T2guess]); end end end % save results P2 = P2guess; T2 = T2guess; function [P5, T5, v5, Cp5, h5] = reflectedShockNR(P5ideal, T5ideal, v5ideal, Rgas5, P2, v2, h2, u2, u1, X, gasMolWt, coefVec, speciesVec, Ru) % This function solves for the conditions behind the reflected shock % Solve the Reflected Shock Equations % P*v = Rgas*T % rho2*u2dash=rho5*u5 % P2+rho2*u2dash^2=P5+rho5*u5^2 % h2+(1/2)*u2dash^2=h5+(1/2)*u5^2 % u2dash=u5+u1-u2 % u2dash=(u1-u2)*(rho5/(rho5-rho2)) % U5=0 (gas relative to lab frame behind reflected shock is stagnant) % by solving the following pair of equations % f3=(P5/P2-1)+(u1-u2)^2/(P2*(v5-v2)) % f4=(h5-h2)/((1/2)*(u1-u2)^2)+(v5+v2)/(v5-v2) % use Newton-Raphson Method % Assign ideal to first guess P5guess = P5ideal; T5guess = T5ideal; v5 = v5ideal; % NR iteration parameters maxTimeNR = 10; % second maxIterNR = 100; % iterations tolNR = 1e-6; % convergence tolerance stopperNR = false; % NR loop iterNR = 0; startTimeNR = tic; while stopperNR == false; % inc the counter iterNR = iterNR+1; % compute cp and h for gas in region 5 (see subfunction below) [Cp5, h5] = getCph_T(T5guess, X, gasMolWt, Ru, coefVec, speciesVec); % J/kg*K, J/kg % compute f3 f3 = (P5guess/P2-1)+(u1-u2)^2/(P2*(v5-v2)); % compute f4 f4 = (h5-h2)/((1/2)*(u1-u2)^2)+(v5+v2)/(v5-v2); % compute derivatives for jacobian % note that dh/dP = 0 (ideal gas) % note that dT/dP = 0 (independent variables) df3dP5 = (1/P2)+(-(u1-u2)^2/(P2*(v5-v2)^2))*(-v5/P5guess); % df3/dP5 df3dT5 = (1/P2)*(0)+(-(u1-u2)^2/(P2*(v5-v2)^2))*(v5/T5guess); df4dP5 = (1/(1/2*(u1-u2)^2))*(0)+(-2*v2/(v5-v2)^2)*(-v5/P5guess); df4dT5 = (1/(1/2*(u1-u2)^2))*(Cp5)+(-2*v2/(v5-v2)^2)*(v5/T5guess); % form Jacobian matrix [df1/dP2 df1/dT2; df2/dP2 df2/dT2] J = [df3dP5 df3dT5; ... df4dP5 df4dT5]; % matrix multiplication to solve for new T5 and P5 guesses % using matlab backslash operator rather than inverse saves time! % general method: [P5new T5new]' = [P5old T5old]' - inv(J)*[f3old f4old]' PTold = [P5guess T5guess]'; fvec = [f3 f4]'; PTnew = PTold - J\fvec; P5new = PTnew(1); T5new = PTnew(2); % check results - the T2 and P2 we guess should give the same T2,P2 checkP = abs((P5guess-P5new)/P5new); % percent change from new checkT = abs((T5guess-T5new)/T5new); % decide whether to keep going or quit if (checkP < tolNR) && (checkT < tolNR) stopperNR = true; % exit the loop elseif toc(startTimeNR) > maxTimeNR; error('Timed out on NR loop') elseif iterNR > maxIterNR error('Itered out on NR loop') else % reassign the variables for the next iteration P5guess = P5new; T5guess = T5new; % recompute specific volume v5 = Rgas5*T5guess/P5guess; % m3/kg % prevent the algorithm from going the wrong way while P5guess < 0 P5guess = mean([P5ideal P5guess]); % average back to ok range end while T5guess < 0 T5guess = mean([T5ideal T5guess]); end end end % save results P5 = P5guess; % Pa T5 = T5guess; % K function phi = getStoich(X, speciesVec, nameEle, numEle) % This function computes the stoichiometry of the shock % phi = (nFuel/nO2)_actual / {nFuel/nO2)_stoich % = sum{xFueli*(Nci+(1/4)Nhi+Nsi-(1/2)Noi)}/xO2 % Stoich = C-->CO2, H-->H2O, S-->SO2 % Inert species list: these are not included in the calculation of phi % N2, AR, etc. need not be listed here because they have no C, H, or S. inertReactantSpecies = {'CO2','H2O','SO2'}; % Find the index of oxygen in the species vector. % If oxygen is not present, or if it is listed more than once, phi = -999 tempArray = strncmp(speciesVec,'O2',2); % vector of 0's with a 1 where O2 is. if sum(tempArray)==0 % oxygen not listed phi = -999; return elseif sum(tempArray)>1 % oxygen listed more than once fprintf('Error: O2 listed more than once in test gas species vector\n') phi = -999; return end oxygenIndex = find(tempArray); % gives index where 1 is present % Compute numerator of phi fraction phiNumerator = 0; for i=1:length(X) % check if this is an inert species speciesName = speciesVec{i}; tempArray = strncmp(inertReactantSpecies,speciesName,length(speciesName)); if sum(tempArray)>0 % this species is inert and does not contribute to phi continue end % check if this is O2 if i==oxygenIndex % oxygen is not considered to be a fuel, it is an oxidizer continue end % determine number of C, H, S, and O atoms in this species nC = 0; nH = 0; nS = 0; nO = 0; for j=1:size(nameEle,2) % looping over columns of nameEle (rows are species, columns are elements) % grab the name of the element (ie, 'H', 'C', 'AR', 'HE') tempEl = nameEle{i,j}; % decide what to do if isempty(tempEl) continue; % if nothing there, skip it. elseif strcmp(tempEl,'C') nC = numEle(i,j); elseif strcmp(tempEl,'H') nH = numEle(i,j); elseif strcmp(tempEl,'O') nO = numEle(i,j); elseif strcmp(tempEl,'S') nS = numEle(i,j); else continue % not one of the atoms we care about ('AR', 'HE', etc.) end end % Compute this part of the numerator and add it to the sum phiNumerator = phiNumerator + X(i)*(nC + (1/4)*nH + nS - (1/2)*nO); end % Compute denominator phiDenominator = X(oxygenIndex); % Compute phi phi = phiNumerator/phiDenominator; % if rich, phi>1, if lean, phi<1 function [coefVec, nameEle, numEle] = readSandia(speciesVec, nasaFile) % This function reads the thermodynamic data file % See the accompanying document explaining the SANDIA.DAT file % there are two fits for data: a high and a low temperature fit % Each fit has 7 coefficients; we need to reassemble these. % Note: we are not reading in the date, and also there is a value in % the 5th position in line 4 which is the enthalpy of formation divided % by the gas constant R; this can be computed using the equations from % the coefficients already read into the file. % for more information, see the following resources: % http://garfield.chem.elte.hu/burcat/THERM.DAT % http://www.galcit.caltech.edu/EDL/public/formats/chemkin.html % http://www.galcit.caltech.edu/EDL/public/formats/nasaold.html % http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19780009781_1978009781.pdf % make a checklist vector to ensure we got everything checkList = zeros(1,length(speciesVec)); % open the file fileID = fopen(nasaFile); % read the first line firstLine = fgets(fileID); % is this the right starter file? tempIndex = strfind(firstLine,'THERMO'); % index where 'THERMO' begins if isempty(tempIndex) % if 'THERMO' does not occur in first line error('%1s is empty',nasaFile) end % read the next line secondLine = fgets(fileID); % find the general temperature limits genTlow = str2num(secondLine(1:10)); % general low temp limit (K) genTcom = str2num(secondLine(11:20)); % general common temp (K) genThigh = str2num(secondLine(21:30)); % general high temp limit (K) % read the file closeNow = false; while (closeNow ~= true) % reset variables elements = {''}; numAtoms = []; Tlow = 0; Tcom = 0; Thigh = 0; TlowWarn = false; TcomWarn = false; ThighWarn = false; % get the first line of the species L1 = fgets(fileID); % check for ending endText = upper(strtrim(L1(1:3))); if strcmp(endText,'END') closeNow = true; % have we found everything? elseif sum(checkList) == length(checkList); % all 1's so this sums to the length if full! closeNow = true; % if not at the end, keep reading the file else % grab the species name speciesName = upper(strtrim(L1(1:18))); % remove white space, make UPPERCASE % check the name; sometimes there is other text beside it refIndex = strfind(speciesName,'REF'); % find where 'REF' occurs if ~isempty(refIndex) % if 'REF' occured, then... speciesName = strtrim(speciesName(1:refIndex-1)); % take a substring end % grab the elements and number of atoms of each elpos = [25, 30, 35, 40, 74]; for i=1:5 % generate indicies tmpIdx1 = elpos(i)+0; %element char 1 tmpIdx2 = elpos(i)+1; %element char 2 tmpIdx3 = elpos(i)+2; %number of elem digit 1 tmpIdx4 = elpos(i)+4; %number of elem digit 3 % grab the element name and the number of atoms tmpEl = upper(strtrim(L1(tmpIdx1:tmpIdx2))); % UPPERCASE tmpNum = str2num(L1(tmpIdx3:tmpIdx4)); % make sure the element is there % not all compounds have 5 elements % if it's there, store it. if not, throw it away. check = strcmp(tmpEl,''); % true if nothing, false if something % some files have a zero in column 74 % this throws off the algorithm so check for it check2 = strcmp(tmpEl,'0'); if (check == false) && (check2 == false) elements{1,i} = tmpEl; numAtoms(1,i) = tmpNum; end end % grab the phase phase = upper(L1(45)); % UPPERCASE letter % grab the low, middle, and high temperatures (K) Tlow = str2num(L1(46:55)); Thigh = str2num(L1(56:65)); Tcom = str2num(L1(66:73)); % check the inputs if sum(size(Tlow))>2 || sum(size(Tcom))>2 || sum(size(Thigh))>2 error('SANDIA.DAT file format issue for %7s; Check for tabs instead of spaces.',speciesName) end % if we did not see Tlow etc, use the general value if isempty(Tlow); Tlow = genTlow; TlowWarn = true; % flag that temperature was replaced end if isempty(Tcom); Tcom = genTcom; TcomWarn = true; % flag that temperature was replaced end if isempty(Thigh); Thigh = genThigh; ThighWarn = true; % flag that temperature was replaced end % get line two L2 = fgets(fileID); % grab the temperature coefficients here Ha1 = str2num(L2(1:15)); % high temperature coefficient a1 Ha2 = str2num(L2(16:30)); Ha3 = str2num(L2(31:45)); Ha4 = str2num(L2(46:60)); Ha5 = str2num(L2(61:75)); % get line three L3 = fgets(fileID); % grab the temperature coefficients here Ha6 = str2num(L3(1:15)); Ha7 = str2num(L3(16:30)); La1 = str2num(L3(31:45)); % low temperature coefficient a1 La2 = str2num(L3(46:60)); La3 = str2num(L3(61:75)); % get line four L4 = fgets(fileID); % grab the temperature coefficients here La4 = str2num(L4(1:15)); La5 = str2num(L4(16:30)); La6 = str2num(L4(31:45)); La7 = str2num(L4(46:60)); % if we want this species, store it, otherwise disregard for i=1:length(speciesVec); % if the same name as one of our wanted species, % and if it is a gas % and if we have not grabbed it yet if strcmp(speciesName,speciesVec{i}) && strcmp(phase,'G') && checkList(i)==0; % coefficeint and temperature matrix coefVec(1,:,i) = [Ha1 Ha2 Ha3 Ha4 Ha5 Ha6 Ha7]; coefVec(2,:,i) = [La1 La2 La3 La4 La5 La6 La7]; coefVec(3,:,i) = [Tlow Tcom Thigh 0 0 0 0]; % element name and number vectors for j=1:length(elements) nameEle(i,j) = elements(1,j); % use parentheses not curly brackets numEle(i,j) = numAtoms(1,j); end % display warnings about changing temperatures if TlowWarn fprintf('WARNING: Defaulting Tlow to %4.2f K for %7s\n',genTlow,speciesName); end if TcomWarn fprintf('WARNING: Defaulting Tcom to %4.2f K for %7s\n',genTcom,speciesName); end if ThighWarn fprintf('WARNING: Defaulting Thigh to %4.2f K for %7s\n',genThigh,speciesName); end % check off that we got it checkList(i) = 1; end end end end % close the file fclose(fileID); % did we get everything? if sum(checkList) ~= length(checkList) for i=1:length(checkList) if checkList(i) == 0 error('Species not in %5s: %5s',nasaFile,speciesVec{i}) end end end