Design a cascade of two fixed-bed catalytic reactors

I have been trying to work on this problem for couple of days now, but i don't seem to understand it. if you look at the problem it requires some conversion, i am not familiar with how i can implement it using Matlab. Can anyone help?

Problem

The aim is to design a cascade of two fixed-bed catalytic reactors to produce hydrogen from carbon monoxide and water:

CO+ H2O <---> CO2 + H2

The reactors are adiabatic and the temperature of the intermediate stream (from the 1st to the 2nd reactor) is reduced by 130 K. (Note: the catalyst undergoes significant deactivation when the temperature rises above 475 °C). Furthermore, the temperature at the exit from the first floor must be 10 K below its maximum equilibrium temperature. Knowing that the composition of the output stream must be less than 0.6% (mol) CO on a dry basis,

calculate:

a) The temperature and conversion profiles within each reactor. Graphically present the result.

b) The mass of catalyst on each floor.

c) The diameter and length of each reactor. Also calculate the cooling power of the exchanger involved.

Problem details:

Feed molar flow rate (dry basis) = 3.5 mol/s

Steam/CO molar ratio in feed = 6

To=300 °C

Po 2760 kPa

Pmax=70 kPa

Pb = 1100 kg/m3

Pp = 2000 kg/m3

Cp 2.12 J/(gK)

H1 = -39900 J/kmol

-rco = k (Pco*PH2O – ((Pco2*PH2)/ Keq))

Keq = 0.008915exp(4800/T)

K (mol /skg bar2) = 0.779exp(-4900/T)

The code i have written uptil now is:

function reactor_design()

    % Given conditions
    T0 = 300;               % Initial temperature in Celsius
    P0 = 2760;              % Initial pressure in kPa
    Pmax = 70;              % Maximum pressure drop in kPa
    Pb = 1100;              % Density of the bulk in kg/m^3
    Pp = 2000;              % Density of the catalyst in kg/m^3
    Cp = 2.12;              % Heat capacity in J/(gK)
    H1 = -39900;            % Enthalpy change in J/kmol
    feed_flow_rate = 3.5;   % Feed molar flow rate in mol/s
    steam_CO_ratio = 6;     % Steam/CO molar ratio in feed
    
    % Other constants
    R = 8.314;              % Gas constant in J/(molK)
    
    % Equilibrium constant equation
    Keq = @(T) 0.008915 * exp(4800 / T);
    
    % Rate constant equation
    k = @(T) 0.779 * exp(-4900 / T);
    
    % Initial conditions
    T1_initial = T0;        % Initial temperature of the first reactor
    
    % Calculate equilibrium temperature for the first reactor
    T1_equilibrium = fzero(@(T) Keq(T) - (P0 / Pmax), T0);
    
    % Solve for temperature profile in the first reactor
    [T1, ~] = ode45(@temperature_profile, [0 1], T1_initial);
    
    % Intermediate stream temperature
    T_intermediate = T1(end) - 130;
    
    % Solve for temperature profile in the second reactor
    [T2, ~] = ode45(@temperature_profile, [0 1], T_intermediate);
    
    % Display results
   figure;
   subplot(2, 1, 1);
    plot(T1, 'LineWidth', 2);
    title('Temperature Profile in the First Reactor');
    xlabel('Reactor Length');
    ylabel('Temperature (°C)');

   subplot(2, 1, 2);
    plot(T2, 'LineWidth', 2);
    title('Temperature Profile in the Second Reactor');
    xlabel('Reactor Length');
    ylabel('Temperature (°C)');
    

    % Temperature profile differential equation
    function dTdt = temperature_profile(~, T)
        dTdt = -H1 / (feed_flow_rate * Cp) + (T1_equilibrium - T) * k(T) * P0 / (Pb * R * T);
    end

end

Answer 1

[Hassaan]

The following is an extended version of your provided script with additional functionalities and comments explaining each step. Remember that while this code is a strong starting point, actual reactor design can be complex, and additional considerations based on real-world constraints might be necessary.

reactor_design()

function reactor_design()
    % Given conditions
    T0 = 300 + 273.15;       % Initial temperature in Kelvin
    P0 = 2760;               % Initial pressure in kPa
    Pmax = 70;               % Maximum pressure drop in kPa
    Pb = 1100;               % Density of the bulk in kg/m^3
    Pp = 2000;               % Density of the catalyst in kg/m^3
    Cp = 2.12 * 1e3;         % Heat capacity in J/(kgK)
    H1 = -39900;             % Enthalpy change in J/kmol
    feed_flow_rate = 3.5;    % Feed molar flow rate in mol/s
    steam_CO_ratio = 6;      % Steam/CO molar ratio in feed
    Mw_CO = 28.01;           % Molecular weight of CO in g/mol

    % Other constants
    R = 8.314;               % Gas constant in J/(molK)
    
    % Equilibrium constant equation
    Keq = @(T) 0.008915 * exp(4800 / T);
    
    % Rate constant equation
    k = @(T) 0.779 * exp(-4900 / T);

    % Function to calculate the rate of reaction
    function r = rate(T, P_CO, P_H2O, P_CO2, P_H2)
        Keq_T = Keq(T);
        k_T = k(T);
        r = k_T * (P_CO * P_H2O - (P_CO2 * P_H2) / Keq_T);
    end
    
    % Differential equations representing the system
    function dYdt = reactorODE(z, Y)
        % Unpack the variables from Y
        T = Y(1);   % Temperature
        
        % Assume partial pressures based on stoichiometry and conversion for simplicity
        P_CO = P0 * (1 - steam_CO_ratio); % This will change based on your system
        P_H2O = P0 * steam_CO_ratio;      % This will change based on your system
        P_CO2 = P0 * steam_CO_ratio;      % This will change based on your system
        P_H2 = P0 * steam_CO_ratio;       % This will change based on your system

        % Reaction rate
        r = rate(T, P_CO, P_H2O, P_CO2, P_H2);
        
        % Energy balance to find temperature change
        dTdz = (-H1 * r) / (feed_flow_rate * Cp * Mw_CO); % Adjust as per molar flow rates
        
        % Update the system of ODEs
        dYdt = [dTdz;];
    end

    % Initial conditions for the ODEs
    Y0 = [T0;];
    
    % Solve the ODEs for the first reactor
    [z1, Y1] = ode45(@reactorODE, [0, 10], Y0); % 10 is a placeholder for reactor length
    
    % Assume intermediate cooling
    T_intermediate = Y1(end, 1) - 130;
    
    % Update initial conditions for the second reactor
    Y0_intermediate = [T_intermediate;];
    
    % Solve the ODEs for the second reactor
    [z2, Y2] = ode45(@reactorODE, [0, 10], Y0_intermediate); % 10 is a placeholder for reactor length

    % Plotting the results
    figure;
    subplot(2, 1, 1);
    plot(z1, Y1(:, 1) - 273.15, 'LineWidth', 2); % Convert back to Celsius
    title('Temperature Profile in the First Reactor');
    xlabel('Reactor Length (m)');
    ylabel('Temperature (°C)');

    subplot(2, 1, 2);
    plot(z2, Y2(:, 1) - 273.15, 'LineWidth', 2); % Convert back to Celsius
    title('Temperature Profile in the Second Reactor');
    xlabel('Reactor Length (m)');
    ylabel('Temperature (°C)');
end

Key Aspects of the Code:

• Differential Equations: The reactorODE function defines the differential equations for the reactor system. It currently only considers the temperature change due to the reaction enthalpy and the rate of reaction.
• Reaction Rate: The rate function calculates the reaction rate based on the provided kinetics and equilibrium constant.
• Solving ODEs: The ode45 solver is used twice, once for each reactor, to solve the system of ODEs.
• Intermediate Cooling: After solving for the first reactor, the temperature is adjusted to account for the intermediate cooling.
• Plotting: The script plots the temperature profiles for both reactors.

Considerations:

• Partial Pressures: This code assumes partial pressures for simplicity. In a real scenario, these would change with the conversion in the reactor.
• Reactor Dimensions: The actual reactor dimensions (diameter and length) and catalyst mass would need to be calculated based on the volumetric flow rates, desired conversion, and catalyst properties.
• Validation: The model needs to be validated with experimental or literature data to ensure its accuracy.
• Complexity: Real reactor systems might require more complex models considering pressure drops, multiple reactions, heat transfer limitations, and more.

Please review and adapt the code according to the specific details of your system and the reaction mechanism.

------------------------------------------------------------------------------------------------------------------------------------------------

If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.

[HHuzzam]

Thank you that was quiet helpful!