Produced using MATLAB software TR 1D model1 SS TR 1D model1 SS m TR 1D model1 SS TR 1D model1 SS m function imain flag TR 1D model1 SS This program calculates the steady state concentration and temperature profiles in a 1 D tubular reactor for an arbitrary number of species and an arbitrary reaction network The reaction network is specified by the stoichiometric coefficients and the exponential powers to which the concentrations of each species are raised in the rate laws The effective diffusivities for each species and the density and heat capacity of the medium are assumed to be constant The heats of reaction are likewise assumed constant and the temperature dependence of each rate constant is specified by the value of the rate constant at a reference temperature and a constant activation energy The heat transfer coefficient for the cooling jacket is assumed constant Dankwert s boundary conditions are applied at the inlet and outlet A constant superficial velocity obtained from knowledge of the reactor dimensions and volumetric flow rate is used to quantify the convective contribution to the fluxes of each species concentration and the enthalpy PROGRAM INPUT OUTPUT DATA problem dimension data struct ProbDim num species IN INT the number of species num rxn IN INT the number of reactions reactor data struct Reactor len IN REAL the length of the tubular reactor dia IN REAL the diameter of the tubular reactor Qflow IN REAL the volumetric flow rate through the reactor Along with the dimensions of the reactor it defines the superficial velocity used in the convective terms of the species and enthalpy balances Temp cool IN REAL the temperature of the reactor coolant jacket U HT IN REAL the overall heat transfer coefficient of the reactor conc in IN REAL ProbDim num species the concentrations of each species in the reactor inlet Temp in IN REAL the temperature of the reactor inlet volume PROG REAL the volume of the reactor cross area PROG REAL the cross sectional area of the reactor surf area PROG REAL the surface area of the reactor available for heat transfer to the cooling jacket velocity PROG REAL the superficial velocity in the reactor that is included in the convective flux terms physical data struct Physical diffusivity IN REAL num species the constant diffusivities of each species density IN REAL the constant density of the medium Cp IN REAL the constant heat capacity of the medium thermal conduct IN REAL the constant thermal conductivity of the medium thermal diff PROG REAL the constant thermal diffusivity of the medium rxn data struct Rxn stoich coeff IN REAL ProbDim num rxn ProbDim num species the stoichiometric coefficients possibly fractional of each species in each reaction ratelaw exp IN REAL ProbDim num rxn ProbDim num species the exponential power possibly fractional to which the concentration of each species is raised each reaction s rate law is rxn elementary IN INT ProbDim num rxn if a reaction is elementary then the rate law exponents are zero for the product species and the negative of the stoichiometric coefficient for the reactant species In this case we need not enter the corresponding components of ratelaw exp since these are determined by the corresponding values in stoich coeff We specify that reaction number irxn is elementary by setting is rxn elementary irxn 1 Otherwise default 0 we assume that the reaction is not elementary and require the user to input the values of ratelaw exp for reaction irxn k ref IN REAL ProbDim num rxn the rate constants of each reaction at a specified reference temperature T ref IN REAL ProbDim num rxn This is the value of the reference temperature used to specify the temperature dependence of each rate constant E activ IN REAL ProbDim num rxn the constant activation energies of each reaction divided by the value of the ideal gas constant delta H IN REAL num rxn the constant heats of reaction PROGRAM IMPLEMENTATION NOTES Section 1 Method of discretizing PDE s To discretize the partial differential equations that describe the balances on the species concentrations and the enthalpy use the method of finite differences To avoid spurious oscillations when convection dominates and the local Peclet number is greater than two use upwind differencing Implement the finite difference procedure so that the grid point spacing may be non uniform grid data struct Grid num pts PIN INT the number of grid points in the axial direction z POUT REAL Grid num pts the values of the z coordinate at the grid points state data struct State conc POUT REAL Grid num pts ProbDim num species the values of the species concentrations at grid points Temp POUT REAL Grid num pts the values of the temperature at each grid point Section 2 Method of solving for the steady state profiles To solve for the steady state profiles we will use a robust two step procedure We will initially assume that the inlet conditions hold uniformly throughout the reactor As this is likely to be far from the true solution we will first perform a number of implicit Euler time integration steps to get within the vicinity of the stable steady state solution The time integration will proceed until a maximum number of time steps have been performed or until the norm of the time derivative vector falls below a specified value If the time derivative has become sufficiently small we will switch to Newton s method with a weak line search to aid global convergence If one wishes to use only Newton s method to solve for the steady state profile for example to find an unstable steady state then Solver max iter time is set to 0 Otherwise if the maximum number of time integration steps has been performed and the time derivative is still too large the program exits without performing any Newton s method iterations A restart utility will be added so that if convergence is not achieved executing the program again will start from the previously saved results Upon a restart new time step and convergence tolerances are input At each time or Newton s method iteration the values of the concentration and temperatures at each grid point will be constrained to be non negative iflag restart PIN INT This integer flag indicates whether the simulation is a restart of a previous simulation in which only new convergence parameters need be input or is an initial simulation in which all system parameters must be input If iflag restart is non zero then it is a restart if 0 then it is an initial simulation imain flag