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
View Full Document