Tubular Reactor with Nonisothermal Cooling Jacket

In this simple example, study an elementary, exothermic, and irreversible reaction in a tubular reactor (in a liquid phase and laminar flow regime). The reactor keeps its temperature down via a cooling jacket. In this model, you investigate the reactor’s steady-state behavior.

The Model Definition section provides a general description of the complete reactor model, whereas the Modeling Instructions detail how to set up and solve a nonisothermal reactor model that accounts for the cooling jacket.

This model is based on the example in Ref. 1.

Model Definition

Reaction

The reaction is a conversion of species A, B, and C in liquid.

(1)

A is the notation for propylene oxide, B is water, and C is propylene glycol. The reaction kinetics are first order in regard to the concentration of A.

Geometry

Figure 1 illustrates the model geometry. We assume that the variations in angular direction around the central axis are negligible, and therefore the model can be axisymmetric.

Figure 1: Geometry for the two-dimensional rotationally symmetric model.

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. The borders of the surface represent the inlet, outlet, reactor wall, and centerline. Assuming that the diffusivity for the three species is of the same magnitude, you can model the reactor using three differential equations: one mass balance for one of the species (as noted in the next section, mass balances are not necessary for the other two species); one heat balance for the reactor core; and one heat balance for the heating jacket. Due to rotational symmetry, the software only needs to solve these equations for half of the domain shown in Figure 1.

Model Equations

You describe the mass balances and heat balances in the reactors with partial differential equations (PDEs) as given in the Transport of Diluted Species and Heat Transfer in Fluids interfaces, while one ordinary differential equation (ODE) is required for the heat balance in the cooling jacket. The latter equation is set up with a Coefficient Form Boundary PDE. The equations are defined as follows.

Mass Balance for Species A

(2)

where Dp denotes the diffusion coefficient, CA is the concentration of species A, U is the flow velocity, R gives the radius of the reactor, and rA is the reaction rate.

In this example, we assume that the species A, B, and C have the same diffusivity, which implies that we must solve only one material balance. We know the other species’ concentrations through stoichiometry.

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

(3)

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

Energy Balance Boundary Conditions

•

Inlet (z = 0)

•

At the wall (r = R)

where Ta denotes the constant temperature in the cooling jacket.

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: the heat transport out of the reactor is purely convective.

•

Outlet (z = L)

Energy Balance of the Coolant in the Cooling Jacket

Here, we assume that only axial temperature variations are present in the cooling jacket. This assumption gives a single ODE for the heat balance:

where Tj is the coolant temperature, mc is the mass flow rate of the coolant, CPc represents its heat capacity, and Uk gives the heat transfer coefficient between the reactor and the cooling jacket. You can neglect the contribution of heat conduction in the cooling jacket and thus assume that heat transport takes place only through convection.

Boundary Condition for the Cooling Jacket

You can describe the cooling jacket with a 1D line. Therefore, you need only an inlet boundary condition.

•

Inlet (z = 0)

Model Parameters

You can define the model’s input data either as constants or as logical expressions. To define the constant’s name, use the left-hand side of the equality in the following list (for example, Diff, for the diffusivity of all species). To define the expression, use the value on the right-hand side of the equality (for example, 1E-9, for Diff).

•

Activation energy, E = 75362 J/mol

•

Frequency factor, A = 16.96E12 1/h

•

Diffusivity of all species, Diff = 1E-9 m2/s

•

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

•

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

•

Inlet temperature, T0 = 312 K

•

Inlet temperature of the coolant, Ta0 = 277 K

•

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

•

Total flow rate, v0 = 0.1 m3/s

•

Mass flow rate of coolant, mc = 0.1 kg/s

•

Concentration A at inlet, cA0 = 1587 mol/m3

•

Concentration B at inlet, cB0 = 43210 mol/m3

•

Concentration methanol at inlet, cMe0 = 1587 mol/m3

•

Heat capacity at inlet, Cp0 = (146.54*cA0_po+75.36*cB0+81.095*cMe0)/rho0 J/(kg·K) (here, the numerical factors are the molar specific heat values in the unit J/(mol·K))

•

Heat capacity per mass of coolant, Cpc = 4180 J/(kg·K)

•

Density at inlet, rho0 = cA0_po*M_po+cB0_po*M_w+cMe0_po*M_po kg/m3

•

Reactor radius, Ra = 0.1 m

•

Reactor length, L = 1 m

•

Density, propylene oxide, rho_po_p = 830 kg/m3

•

Density, methanol, rho_m_p = 791.3 kg/m3

•

Density, water, rho_w_p = 1000 kg/m3

The following section lists the definitions for the model expressions. Again, to put each expression in COMSOL Multiphysics form, use the left-hand side of the equality (for instance, u0) for the variable’s . Use the right-hand side of the equality (for instance, v0/(pi*Ra^2)) for its .

•

The superficial flow rate is defined according to the analytical expression

which you define in COMSOL Multiphysics as u0 = v0/(pi*Ra^2).

•

The superficial, laminar flow rate

becomes uz = 2*u0*(1-(r/Ra)^2).

•

The conversion of species A is given by

which in COMSOL Multiphysics form is xA = (cA0-cA)/cA0.

•

The concentration of species B is according to

which in COMSOL Multiphysics form becomes cB = cB0-cA0*xA.

•

The concentration of species C is expressed as

which becomes cC = cA0*xA.

•

The rate of reaction takes the following form:

which in COMSOL Multiphysics form is rA = -A*exp(-E/R_const/T)*cA.

•

The heat production term becomes

which in COMSOL Multiphysics form is Q = -rA*dHrx.

Results

The following figures collect the results as shown in Ref. 1.

Surface plots for the surface temperature and conversion are shown in Figure 2 and Figure 4. These show that where the temperature is low, little conversion takes place and vice versa. This is 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 is shown in these figures as well.

Figure 2: Temperature surface.

Figure 3: Radial temperature surface profiles.

Figure 4: Conversion surface.

Figure 5: Radial conversion surface profiles.

Exercises

Try these example exercises with the model to better understand the system:

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

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

Expand the model. Is the convection in the radial direction important?

Reference

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

COMSOL_Multiphysics/Chemical_Engineering/tubular_reactor

Modeling Instructions

When you start COMSOL Multiphysics, 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 if 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 sets up the required mass balance equation for species A.

Click .

In the text field, type cA.

cA is the dependent variable name.

In the tree, select .

Click .

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

Finally, select the to model the cooling jacket. Tc is the temperature of the coolant.

In the tree, select .

Click .

In the table, enter the following settings:


Tj

Click

In the dialog box, type temperature in the text field.

Click

In the tree, select .

Click .

In the window, In the table, enter the following settings:


Source term quantity	Unit
Custom unit	W/m

click

In the tree, select .

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

Click

Geometry 1

Start by defining the reactor geometry. In 2D axisymmetry, 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 0.1.

Click

Click the

button in the toolbar.

Root

Now, move on to define model-specific constants and expressions. You can type in constant names and their values in the section. 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.

Global Definitions

Parameters 1

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 tubular_reactor_parameters.txt.

Just as with the model constants, you will find it convenient to group user-defined expressions in a list. You can type in expressions that contain constants from the list as well as dependent variables that you solve for; for example, cA.

Definitions

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

Variables 1

In the toolbar, click

and choose .

In the window for , locate the section.

Click

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

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 species A. As you click the node, the section of the window will tell you which equations are being solved for. The section 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 input defining the velocity field of the reacting mixture as well as the diffusivity of species A. In this model the velocity field is given by an expression describing a parabolic laminar flow profile. The variable names you type in have previously been defined in the section and the section. The velocity field could also be calculated by COMSOL Multiphysics using a interface.

Transport Properties 1

In the window, under click .

In the window for , locate the section.

Specify the u vector as


0	r
uz	z

Locate the section. In the DcA text field, type Diff.

Initial Values 1

In the window, click .

In the window for , locate the section.

In the cA text field, type cA0.

At this point, right-click the node, select a feature and associate it with the reactor domain. Click the reactor domain to highlight it in red. Right-clicking the same domain changes the color to blue, meaning the domain is selected. The number of selected domains appears in the list of the node. This adds a sink term to the mass balance equations that takes the depletion of species A through chemical reaction into account.

Reactions 1

In the toolbar, click

and choose .

Select Domain 1 only.

In the window for , locate the section.

In the RcA text field, type rA.

Now that the domain equations have been defined for the model, it is time to set the boundary conditions. First, you will select a concentration inflow condition at the inlet boundary.

Concentration 1

In the toolbar, click

and choose .

Select Boundary 2 only.

In the window for , locate the section.

Select the check box.

In the c0,cA text field, type cA0.

Outflow 1

In the toolbar, click

and choose .

Select Boundary 3 only.

Assigning the condition to the outlet boundary imposes −n·D∇c = 0, that is, the transport of mass across the boundary is dominated by convection. Note that the mathematical representation of the boundary conditions is 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 species A. Next, set up the interface.

Heat Transfer in Fluids (ht)

Fluid 1

In the window, under click .

In the window for , locate the section.

Specify the u vector as


0	r
uz	z

In addition to the velocity field, the feature asks for the thermal conductivity, density, and heat capacity of the fluid.

Locate the section. From the k list, choose . In the associated text field, type ke.

Locate the section. From the list, choose .

From the ρ list, choose . In the associated text field, type rho0.

From the Cp list, choose . In the associated text field, type cpm.

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 to the heat balance.

Select Domain 1 only.

In the window for , locate the section.

In the Q0 text field, type Q.

Next, add the boundary conditions specifying a temperature at the inlet, the heat flux between the 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.

In the q0 text field, type -Uk*(T-Tj).

Outflow 1

In the toolbar, click

and choose .

Select Boundary 3 only.

Finally, add an energy balance describing the temperature distribution in the cooling jacket.

Coefficient Form Boundary PDE (cb)

In the window, under click .

In the window for , locate the section.

In the list, choose , , and .

Click

Select Boundary 4 only.

Coefficient Form PDE 1

In the window, under click .

In the window for , locate the section.

In the c text field, type 0.

Locate the section. In the f text field, type 2*pi*Ra*Uk*(T-Tj).

Click to expand the section. Specify the β vector as


0	r
Cpc*mc	z

Dirichlet Boundary Condition 1

In the toolbar, click

and choose .

The only required boundary condition is the temperature at the inlet of the jacket, Ta0.

Select Point 3 only.

In the window for , locate the section.

In the r text field, type Ta0.

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 mesh. The software uses the mesh when applying the finite element method to numerically solve the differential equations. In this particular model, you will create a mapped 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 the sharp concentration and temperature gradients expected when the reactor is run under nonisothermal conditions.

Mapped 1

In the toolbar, click

Distribution 1

Right-click and choose .

First set up 50 vertical mesh lines by selecting the inlet and outlet boundaries and using predefined distribution settings. Then, in the same fashion, set up the horizontal lines to complete the mapped mesh.

Select Boundaries 2 and 3 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 50.

In the text field, type 0.01.

From the list, choose .

Select the check box.

Distribution 2

In the window, right-click and choose .

Select Boundaries 1 and 4 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 200.

In the text field, type 0.01.

From the list, choose .

Select the check box.

In the window, right-click and choose .

The figure below shows the created mesh, containing 10,000 elements.

Study 1

The model will be solved using two study steps. First, the software solves for the mass balance with the temperature kept constant. The computed solution is then used as an initial guess when solving the coupled mass and energy balance equations. This step-wise approach is often useful for tightly coupled equation systems, as a good initial guess helps to improve numerical convergence. It is straightforward to set up the mentioned solver sequence by defining two separate study steps under the node.

The node has a single step set up as a subnode. This study was generated as a consequence of selections in the , initiating the model building process. To set up a two-step solution process, add a second step to the node.

Stationary 2

In the toolbar, click