COMSOL 6.4 - Multicomponent Tubular Reactor with Isothermal Cooling

Multicomponent Tubular Reactor with Isothermal Cooling

This example uses the Chemical Reaction Engineering Module to study an elementary, exothermic, irreversible reaction in a tubular reactor (liquid phase, laminar flow regime). The reactor uses a constant temperature cooling jacket to keep its temperature down. The steady-state behavior of the reactor is investigated. The reaction kinetics and physical properties of the species are modeled with the Chemistry interface, which is available in the Chemical Reaction Engineering Module.

The Model Definition section provides a general description of the complete reactor model, whereas the Modeling Instructions details how to set up and solve the model.

Model Definition

Reaction

The reaction is a reversible liquid phase conversion of propylene oxide and water into propylene glycol in the manner of

(1)

In the species named the prefix “p” stands for propylene. Water is in excess and is modeled as a solvent. The reaction kinetics is first order with respect to the concentration of propylene oxide:

(2)

Geometry

Figure 1 illustrates the tubular reactor geometry.

Figure 1: Model geometry for the 2-dimensional rotationally symmetric models.

The system is described by a set of differential equations on a 2D surface that represents a cross section of the tubular reactor in the rz-plane. That 2D surface’s borders represent the inlet, the outlet, the reactor wall, and the centerline. The reactor model uses the following differential equations: mass balances for the species, a heat balance, combined with momentum and mass balance for the fluid flow. Due to rotational symmetry, the software need only to solve these equations for half of the domain as shown in Figure 1.

Model Equations

You describe the mass balances and heat balances in the reactors with partial differential equations (PDEs). The equations are defined as follows.

Mass Balance

(3)

where Di denotes the diffusion coefficient of species i, Ci is the concentration, u the flow velocity vector, and Rx is the reaction rate.

Using a two-dimensional axially symmetric geometry, and a constant diffusion coefficient, this equation corresponds to

(4)

where u and w denote the radial and axial velocity component respectively.

The mass balance equations are set up and solved for using the Transport of Diluted Species interface.

Mass Balance Boundary Conditions

•

Inlet (z = 0)

•

At the wall (r = R)

The boundary condition selected for the outlet states that convection dominates transport out of the reactor. Thus this condition keeps the outlet boundary open and does not set any restrictions on the concentration.

•

Outlet (z = L)

where L denotes the length of the reactor.

Energy Balance Inside the Reactor

(5)

where ρ denotes the density, CP equals the heat capacity, k is the thermal conductivity, T the temperature, and ΔHRx is the reaction enthalpy.

Using a two-dimensional axially symmetric geometry, and a constant thermal conductivity, this equation corresponds to

(6)

The energy balance is solved for by the Heat Transfer in Fluids interface.

Energy Balance Boundary Conditions

•

Inlet (z = 0)

•

At the wall (r = R)

where Ta denotes the constant temperature in the cooling jacket, and Uk is an overall heat transfer coefficient.

Mass Balance

As for the mass balance, choose the boundary condition at the outlet for the energy balance such that it keeps the outlet boundary open. This condition sets only one restriction, that the heat transport out of the reactor is convective.

•

Outlet (z = L)

Momentum Balance

The fluid flow is modeled with the Laminar flow interface that solves the Navier–Stokes equations computing the velocity and pressure. At the inlet, a fully developed laminar flow profile with an average flow velocity is prescribed. At the outlet, the pressure is prescribed.

Model Parameters

Below is a list of the model’s input data. You define them either as constants or as expressions involving other constants. In defining each parameter in COMSOL Multiphysics, for the constant’s , use the left side of the equality in the following list and use the value on the right side of the equality for the that defines it. Type the unit inside brackets, like this [mol/m^3].

The constants in the model are:

•

Activation energy, E = 75362 J/mol

•

Frequency factor, A = 16.96E12 1/h

•

Overall heat transfer coefficient, Uk = 1300 W/(m2·K)

•

Thermal conductivity of mixture, ke = 0.559 W/(m·K)

•

Inlet temperature, T0 = 312 K

•

Inlet temperature of the coolant, Ta0 = 273 K

•

Heat of reaction, ΔHRx, dHrx = -84666 J/mol

•

Total flow rate, v0 = 0.1[mol/s]/cpoxide0

•

Average flow velocity, u0 = v0/(pi*Ra^2)

•

Concentration of propylene oxide at inlet, cpoxide0 = rho_poxide/M_poxide/9[1]

•

Concentration of water, cwater0 = rho_H2O/M_H2O*(7/9)[1]

•

Molar heat capacity of water, cpm_H2O = 74.5 J/(mol·K)

•

Reactor radius, Ra = 0.1 m

•

Reactor length, L = 1 m

•

Molar weight of propylene oxide, M_poxide = 58.095 g/mol

•

Molar weight of water, M_H2O = 18 g/mol

•

Molar weight of propylene glycol, M_pglycol = 76.095 g/mol

•

Density of propylene oxide, rho_poxide = 830 kg/m3

•

Density of water, rho_H2O = 1000 kg/m3

•

Density of propylene glycol, rho_pglycol = 1040 kg/m3

•

Reference dynamic viscosity of water (at 293 K), myref_H2O = 1 mPa·s

•

Molar heat capacity of water, cpm_H2O = 75.36 J/mol/K

•

The conversion of poxide is given by

All necessary reaction kinetics and mass transport properties are incorporated into the model with the Chemistry interface.

Results

Surface plots for the surface temperature and conversion are shown in Figure 2 and Figure 4. These figures show that where the temperature is low, little conversion takes place and vice versa. This occurs because the rate of the reaction is temperature dependent. The low temperature closest to the wall is due to the coolant.

Figure 3 and Figure 5 show the temperature and conversion surface profiles at three locations along the length of the reactor. The further along the reactor the reactants travel, the more reaction takes place and the higher the temperature becomes. The impact of the coolant are shown in these figures as well.

Figure 2: Temperature in the reactor.

Figure 3: Radial temperature profiles.

Figure 4: Conversion of propylene oxide in the reactor.

Figure 5: Radial conversion profiles for propylene oxide.

Exercises

Some example exercises below can easily be performed with the model to understand the system better.

How does the thermal conductivity of the mixture affect the temperature distribution?

How does the coolant temperature decrease the mixture temperature at the outlet?

Add more reacting components, for example, a side reaction to the Chemistry interface. How will this affect the results?

Use more complicated reaction kinetics. The reaction can be shifted to a second order reversible reaction.

References

1. S. Fogler, Elements of Chemical Reaction Engineering 4th ed., p. 557, Example 8-12 Radial Effects in Tubular Reactor, Prentice Hall, 2005.

Chemical_Reaction_Engineering_Module/Tutorials/multicomponent_tubular_reactor

Modeling Instructions

Starting you are greeted by the . Here you choose the dimension of your model geometry as well as the physics interfaces required. You can return to the later in the modeling process should you want to expand your model to include additional physics interfaces.

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

This interface can be used to calculate reaction kinetics and thermal and mass transport properties.

Click .

In the tree, select > .

This sets up the required mass balance equation for propylene oxide and propylene glycol. Water is a solvent and is not accounted for here.

Click .

In the text field, type 2.

In the table, enter the following settings:


cpoxide
cpglycol

cpoxide and cpglycol are the dependent variable names, where p stands for propylene.

In the tree, select > .

Selecting this physics interface adds an energy balance to the model.

Click .

Also, set up the interface to describe the fluid flow in the reactor.

In the tree, select > > .

Click .

Click

In the tree, select > .

The analysis type lets you investigate the steady-state behavior of the reactor.

Click

Global Definitions

Parameters 1

Start by adding the . You can type in constant names and their values in the dialog. Note that you can enter units enclosed in brackets after the constant values. This can be very useful as the software is able to keep track of unit consistency throughout the model setup procedure.

In this case, the model parameters are available in a text file that is imported.

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file multicomponent_tubular_reactor_parameters.txt.

Geometry 1

Now, move on to define the reactor geometry. In , the representation of the tubular reactor is reduced to a simple rectangle.

Rectangle 1 (r1)

In the toolbar, click

The geometry is automatically drawn as you leave the node. You can also click the button in the toolbar.

In the window for , locate the section.

In the text field, type Ra.

In the text field, type L.

Click

Add a variable computing the conversion of propylene oxide.

Definitions

Variables 1

In the toolbar, click

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
Xpoxide	(cpoxide0-cpoxide)/cpoxide0		Conversion of propylene oxide

Select the mixture type and enable transport mixture properties to be computed.

Chemistry (chem)

First the reaction kinetics and mixture properties will be set up, this is done in the interface.

In the window, under click .

In the window for , locate the section.

From the T list, choose .

Locate the section. From the list, choose .

Click to expand the section.

Reaction 1

In the toolbar, click

and choose .

The reaction is irreversible and contains three species: propylene oxide (poxide), water, and propylene glycol (pglycol).

In the window for , locate the section.

In the text field, type poxide+H2O=>pglycol.

Click .

Locate the section. From the list, choose .

In the rj text field, type chem.kf_1*chem.c_poxide.

Locate the section. Find the subsection. In the text field, type 1.

Define the reaction expression and use the in-built Arrhenius expression for calculation of the reaction rate constant.

Locate the section. Select the checkbox.

In the Af text field, type A.

In the Ef text field, type E.

Locate the section. From the list, choose .

In the Q text field, type -chem.r_1*dHrx.

Under each species, their respective chemical formulas can be entered. Entering a chemical formula gives the species’ molar mass and enables balancing the reaction.

Species: poxide

In the window, click .

In the window for , locate the section.

Select the checkbox.

In the text field, type C3H6O.

In the ρ text field, type rho_poxide.

Species: H2O

Since the species’ name (H2O) is a chemical formula, the chemical formula field is already filled in. As a result, so is the molar mass.

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the ρ text field, type rho_H2O.

Click to expand the section. In the k text field, type ke.

Click to expand the section. From the list, choose .

In the Cp text field, type cpm_H2O.

Species: pglycol

In the window, click .

In the window for , locate the section.

Select the checkbox.

In the text field, type C3H8O2.

In the ρ text field, type rho_pglycol.

1: poxide + H2O => pglycol

In the window, click .

In the window for , locate the section.

Click in the upper-right corner of the section.

As can be seen, no coefficients appeared in front of any of the species when balancing the reaction, which means that 1 mole of propylene oxide and water is needed to form 1 mole of propylene glycol.

In the window, click .

In the window for , locate the section.

Find the subsection. From the list, choose .

In the table, enter the following settings:


Species	Type	Molar concentration	Value (mol/m^3)
H2O	Solvent	User defined	cH2O0
pglycol	Variable	cpglycol	Solved for
poxide	Variable	cpoxide	Solved for

Locate the section. In the μref text field, type myref_H2O.

In the Tref text field, type Tref_my.

Transport of Diluted Species (tds)

In the next step of the model setup you will specify the parameters and source terms needed for the mass balance equation defined for propylene oxide and propylene glycol. As you click the node the section of the window will tell you which equations are solved for. The shows a list of the geometry domains to which the equations apply. Note that you can change the mass transport mechanisms included in the mass balance equation through selections in the section. This can be done at any time in the modeling process.

Moving on to the node, you are expected to provide diffusivity of propylene oxide and propylene glycol. The variable names you type in have previously been defined in the and lists. Here, it is also possible to couple the interface to other interfaces. In this example, you instead use the node to do this.

Fluid 1

In the window, under > click .

In the window for , locate the section.

In the Dcpoxide text field, type chem.D_poxide.

In the Dcpglycol text field, type chem.D_pglycol.

Reactions 1

In the toolbar, click

and choose .

The reaction rates are selected from the interface.

Select Domain 1 only.

In the window for , locate the section.

From the list, choose .

Inflow 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the c0,cpoxide text field, type cpoxide0.

Selecting the Danckwerts boundary condition is a manner to speed up the computation and improve the solution.

Locate the section. From the list, choose .

Outflow 1

In the toolbar, click

and choose .

Select Boundary 3 only.

Assigning the condition to the outlet boundary imposes −n*Dgrad(c)=0, that is, the transport of mass across the boundary is dominated by convection. Note that the mathematical representation of the boundary conditions are displayed in the section of the window. The boundary conditions for the axis of symmetry as well as the no flux condition for the reactor wall are set by default.

This concludes the definition of the mass balance for propylene oxide and propylene glycol. Now, move on to set up the interface.

Heat Transfer in Fluids (ht)

Fluid 1

The feature asks for the thermal conductivity, density, and heat capacity of the fluid mixture. These are taken from the interface.

In the window, under > click .

In the window for , locate the section.

From the k list, choose .

Locate the section. From the list, choose .

From the ρ list, choose .

From the Cp list, choose .

From the γ list, choose .

Initial Values 1

In the window, click .

In the window for , locate the section.

In the T text field, type T0.

Heat Source 1

In the toolbar, click

and choose .

Add a feature to include the effect of the exothermic reactions as defined in the interface to the heat balance.

Select Domain 1 only.

In the window for , locate the section.

From the Q0 list, choose .

Next, add the boundary conditions specifying a temperature at the inlet, the heat flux between reactor and cooling jacket, and an outflow condition at the outlet.

Temperature 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

In the T0 text field, type T0.

Heat Flux 1

In the toolbar, click

and choose .

Select Boundary 4 only.

In the window for , locate the section.

From the list, choose .

In the h text field, type Uk.

In the Text text field, type Ta0.

Outflow 1

In the toolbar, click

and choose .

Select Boundary 3 only.

Now, move on to the interface.

Laminar Flow (spf)

Fluid Properties 1

In the window, under > click .

In the window for , locate the section.

From the ρ list, choose .

From the μ list, choose .

Assume that the flow has a pattern as it enters the reactor.

Inlet 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

From the list, choose .

Locate the section. In the Uav text field, type u0.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 3 only.

Last, couple the interfaces with the node.

Multiphysics

Nonisothermal Flow 1 (nitf1)

In the toolbar, click

and choose > .

Reacting Flow, Diluted Species 1 (rfd1)

In the toolbar, click

and choose > .

This completes the setup of the physics interfaces. The next step of the modeling process involves meshing.

Mesh 1

Following the steps below you will discretize the geometry with a . The software uses the mesh when applying the finite element method to numerically solve the partial differential equations. In this particular model you will create a mesh. This meshing technique is often a good choice for simple geometries as it allows detailed control over the mesh distribution. The mesh is dense near the reactor inlet and reactor outer wall. This is needed to resolve sharp concentration and temperature gradients expected when the reactor is run under nonisothermal conditions.

Mapped 1

In the toolbar, click