//this is a file in the VAc dynamic model software package //this is the main file for the reactor //parent-routine: VAmodel.c //subroutines involve: reactor.c reaction_rate.c Cp.c //Copyright: Rong Chen and Kedar Dave, June 2002 //VERSION 1.0 //if a user wants to build a DLL for the reactor only, follow the instruction below: //1. Add this line to the head of this file: #include "mex.h" //2. Enable mexFunction(), which is the second routine included in this file (currently disabled) //3. Use the following line to generate a DLL file: reactor.dll // mex reactor.c reaction_rate.c Cp.c //4. The name of the generated function is "reactor()". This function should be called in MATLAB. //The plug-flow reactor is modeled as a 10-section distributed model #include #include //#include "mex.h" void Reactor(double *dstatedt, double states[], double time, double Fin, double Yin[], double Tin, double Pin, double Cin, double deltap, double Shell_T, double Coefficient1[], double Coefficient2[], int NR, double tube_L, int tube_number, double tube_diameter, double tube_area, double UA, double cata_bulk_density, double cata_heatcapacity, double cata_porosity, double cata_weight, double E_r1, double E_r2, double MW[], double HVAP[], double HCAPLa[], double HCAPLb[], double HCAPVa[], double HCAPVb[], double T_ref) { /* model - This is a nonlinear model of a catalytic plug flow reactor with */ /* 10 sections. Two exothermic reactios and seven components are */ /* in the system including an inert. Gas input and gas output. */ /* Model assumptions: */ /* 1. Heat is removed by shell steam so it's not an adiabatic reactor */ /* 2. The heat transferred between bed and gas is ignored, therefore T_bed=T_gas like homogeneous */ /* 3. Gas drop is linearized without compensation */ /* 4. Ideal gas law is enforced */ /* This dynamic model is based on the equations (18) to (21) in Luyben's paper */ /* */ //-------------------------------------------------------------------------- //physical properties: MW, HVAP, HCAPLa, HCAPLb, HCAPVa, HCAPVb, T_ref //-------------------------------------------------------------------------- /*output: dstatedt: state derivatives --------------------------------------------------------------------------*/ /*input: states: states time: current running time in minute Fin: reactor feed flowrate Yin: reactor feed composition Tin: reactor feed temperature Pin: reactor feed pressure Cin: reactor feed concentration deltap: total pressure drop Shell_T: shell temperature ------------------------------------------------------------------------*/ /*Equipment Specification Coefficient1: stoichimetric coefficients for reaction 1 Coefficient2: stoichimetric coefficients for reaction 2 NR: section number + 1 tube_L: tube length tube_number: tube number tube_diameter: tube diameter tube_area: tube cross-section area UA: Inverse of the Total Thermal Resistance cata_bulk_density: catalyst bulk density cata_heatcapacity: catalyst heat capacity cata_porosity: catalyst porosity cata_weight: catalyst weight E_r1: heat of reaction for reaction 1 E_r2: heat of reaction for reaction 2 ------------------------------------------------------------------------*/ //local variables //states double C_O2[11], C_CO2[11], C_C2H4[11], C_VAc[11], C_H2O[11], C_HAc[11], Tg[11]; //intermediate variables double C_C2H6[11], CC[11], y[11][7], cp[11][7], P_profile[11], Gas_density[11], FCP[11], CCP[11], v[11], Q_rct[11], r_all[11][2]; double dp; //prossure drop in each section double dz; //section length double G; //specific mass flowrate double avg_MW; //average molecular weight double PR[7]; //partial pressure double Heat_Cap[7]; //partial heat capacity //derivatives double C_O2_t[11], C_CO2_t[11], C_C2H4_t[11], C_VAc_t[11], C_H2O_t[11], C_HAc_t[11], Tg_t[11]; //others int i, j, n; double sum1, temp_T[1], dummy[7], r[2]; /*Calculate pressure for each section, assuming linear drop*/ dp=deltap/(NR-1); P_profile[0]=Pin; for (i=1;i