A Multiscale 3D Packed Bed Reactor

The packed bed reactor is used in heterogeneous catalytic processes and is one of the most common reactors in the chemical industry. Its basic design is a column filled with porous catalyst particles, and in some cases the reactor also has a specially designed bottom plate through which the reaction mixture enters. The catalyst particles can be contained within a supporting structure, such as tubes or channels, or they can be packed in one single compartment in the reactor.

Figure 1: An example of the macroscale (bed volume with entry holes) and the microscale (pellet) of a packed bed reactor.

The bed with the packed catalyst particles makes the modeling of mass transport and reactions in the reactor a challenge. The challenge is that species transport and reaction occur in dimensions of different orders of magnitude:

•

In the macropores between the dumped pellets, and

•

inside the catalyst pellets in micropores.

As such, the problem is regarded as a multiscale problem. The Reactive Pellet Bed feature, available with the Transport of Diluted Species interface, is dedicated to these multiscale problems.

The structure between particles in the bed is described as a macroporous material of meter dimensions. The particle radii are often in the order of 1 mm. The pores inside the catalyst particles form the microscale structure of the bed. The pore radii in the particles are often between 1 and 10 microns. There are two porosities that are important: bed porosity (macroscale) and pellet porosity (microscale). Sometimes such models are called double-porosity models.

When a pressure drop is applied across the bed, flow and convection of the fluid is initiated in the bed. The transport of chemicals inside the pellets are dominated by diffusion.

This model is an extension to the 1D example, Packed Bed Reactor, which contains more complex reactions.

Model Definition

A model geometry made up of one eighth of the reactor in Figure 1 can be used due to symmetry. The geometry is shown in Figure 2.

Figure 2: The packed bed reactor simulation geometry. Symmetry observations enables modeling of 1/8 of the true geometry. The results will be expanded to the true geometry with aid of a sector data set.

The reversible catalytic chemical reaction occurs inside a pellet. The reactant species A and B forms a product C:

A + B ↔ 2C

The reaction kinetics are assumed to be equimolecular and are set up with the Chemistry interface. The automatic reaction rate can thus be used and has the following form:

where k is the rate factor (SI unit: m3/(mol·s)) with the superscripts f and r denoting the forward and reverse reaction, respectively. ci is the concentration (SI unit: mol/m3) of species i. The forward reaction constant is defined with the inbuilt Arrhenius expression and the reverse is computed with the equilibrium constant of the reaction.

The mass transport of the reacting species in the reactor is modeled with the Transport of Diluted Species interface, which accounts for diffusion, convection, and reaction in diluted solutions. The species are assumed to be diluted in water.

The reaction inside the pellets is added to the mass balances in the Transport of Diluted Species interface with the Reactive Pellet Bed feature. This feature has a predefined extra dimension (1D) on the normalized radius (r = rdim/rpe) of the pellet particle. The mesh on the extra dimension has a default of 10 elements with a cubic root sequence distribution. If spherical pellets are selected, the following spherical diffusion/reaction equation is set up and solved along the pellet radius for each species i:

(1)

Here, r is a dimensionless radial coordinate that goes from 0 (center) to 1 (pellet surface), rpe is the pellet radius, and N the number of pellets per unit volume of bed. The advantage of formulating Equation 1 on a dimensionless 1D geometry is that the pellet radius can be changed without changing the geometry limits.

Dpe is an effective diffusion coefficient (SI unit: m2/s) and Rpe, i is the reaction source term (SI unit: mol/(m3·s)). Note that the latter term is taken per unit volume of porous pellet material.

At the pellet-fluid interface, a film condition assumption is made. The flux of mass across the pellet-fluid interface into the pellet is possibly rate determined by the resistance to mass transfer on the bulk fluid side. The resistance is expressed in terms of a film mass transfer coefficient, hDi, such that:

(2)

where Ni, inward is the molar flux from the free fluid into a pellet and has the unit moles/(m2·s). The mass transfer coefficient to calculated automatically as described in the section Theory for the Reactive Pellet Bed in the Chemical Reaction Engineering Module User’s Guide.

The pressure drop in the reactor is also accounted for and is modeled with the Darcy’s Law interface.

In Table 1 the model parameters are tabulated.

Table 1: Summary of input data.
Property	Value	Description
Hr	1 [m]	Height of the packed bed reactor
Rr	0.2 [m]	Radius of packed bed reactor
ρb	0.51 [g/cm3]	Density of packed bed
ρpe	0.68 [g/cm3]	Density individual pellet
εb	1-ρb/ρpe	Macroscale porosity (of bed)
εpe	0.70 (-)	Microscale porosity (of pellet)
rpe	0.5 [mm]	Pellet radius (spherical shape)
Dpe,A	1.5e-9 [m2/s]	Diffusion coefficient of A in pellet
Dpe,B	2e-9 [m2/s]	Diffusion coefficient of B in pellet
Dpe,C	0.5e-9 [m2/s]	Diffusion coefficient of C in pellet
A	2e12 [m3/(mol s)]	Frequency factor reaction
E	75000[J/mol]	Activation energy reaction
Keq0	1000	Equilibrium reaction constant
kappa	1.88e-10[m2]	Permeability of Bed
CA,in	1[mol/m3]	Inlet concentration A
CB,in	1[mol/m3]	inlet concentration B
CC,in	0[mol/m3]	inlet concentration C
DA	1e-8 [m2/s]	Diffusion coefficient of A in bed
DB	1.5e-8 [m2/s]	Diffusion coefficient of B in bed
DC	0.5e-8 [m2/s]	Diffusion coefficient of C in bed
pDarcy	0.4 [atm]	Inlet pressure offset

Results and Discussion

The following figures display the results at 180 s. Figure 3 shows the velocity distribution in the fluid between the pellets.

Figure 3: Velocity distribution on the macroscale.

Figure 4 shows the macroscale concentration of the reactant A in the bed column fluid. The species is consumed due to the catalytic chemical reaction in the pellets.

Figure 4: Concentration of reactant A.

Streamline plots can be useful to get an understanding of the flow pattern. It can be seen from Figure 5 that no recirculation occurs at the entry holes. The fluid is evenly spread out in the bed chamber as it enters the holes in the bottom place.

Figure 5: The streamlines show how the fluid enters the holes and then spread out in the bed volume as it. The colors of the lines represent the reactant concentration in moles/m3.

A line plot of the concentration in a pellet at a certain position in the bed is interesting in order to understand the local reaction. Figure 6 shows the position at which the pellet line plot is sampled: (x = 0.5, z = 0, y = 0), and Figure 7 is the line plot of both the reactant and the product inside a pellet in the same position.

Figure 6: Coordinate at which the pellet plot is sampled: Centerline of reactor and at 0.5 m height.

Figure 7: Concentration of species A, B and C within the pellet at 0.5 m bed height.

In Figure 8, the concentrations of the species are shown along the reactor height in the center of the geometry. Both the concentrations in the bed and averaged in the pellets are shown and illustrates the local reaction in detail. The species C concentration profiles portray a reaction intense zone within the reactor. A closer look at this zone shows it expanding toward the outlet with time.

Figure 8: Concentration of the species in reactor bed and averaged within the pellets along the reactor height.

Figure 9 shows a 3D concentration plot of the product C within the pellet at the sampled coordinate. It can be seen that the concentration is higher closer to the center of the pellet, where products build up and from where these diffuse into the bulk gas.

Figure 9: Concentration of species C at 0.5 m bed height.

Chemical_Reaction_Engineering_Module/Reactors_with_Porous_Catalysts/packed_bed_reactor_3d

Modeling Instructions

Start by adding the necessary physics interfaces for a 3D model.

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

In the tree, select .

Click .

In the text field, type 3.

In the table, enter the following settings:


cA
cB
cC

In the tree, select .

Click .

Click

In the tree, select .

Click

Geometry 1

Add the model parameters from a text file.

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

Geometry 1

Now create the geometry. You can simplify this by inserting a prepared geometry sequence from a file with prepared geometry selections.

In the window, expand the node.

Right-click and choose .

Browse to the model’s Application Libraries folder and double-click the file packed_bed_reactor_3d_geom_sequence.mph.

In the toolbar, click

Specify the material properties. Some properties can be found in the COMSOL built-in materials, other are manually entered.

Assume the reaction mixture has mainly aqueous properties.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

to close the window.

Chemistry (chem)

Start with the interface and create the needed reaction kinetics expressions by typing in the reaction formulas.

In the window, under click .

In the window for , locate the section.

From the list, choose .

Reaction 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type A+B<=>2C.

Click .

Locate the section. Select the check box.

Select the check box.

In the Af text field, type A.

In the Ef text field, type E.

Locate the section. From the list, choose .

In the Keq0 text field, type Keq0.

The molar masses for the reacting species can be entered for possible future use. For example, if the mass-based Concentrations feature is used in the interface, it can pick up the molar mass values from the node automatically.

Species: A

In the window, click .

In the window for , locate the section.

In the M text field, type Mn_A.

Species: B

In the window, click .

In the window for , locate the section.

In the M text field, type Mn_B.

Species: C

In the window, click .

In the window for , locate the section.

In the M text field, type Mn_C.

The reactive species are diluted in water. For completeness, add the solvent H2O, which does not partake in the reactions. It can be used later if the model is extended.

Species 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type H2O.

Locate the section. From the list, choose .

Locate the section. In the M text field, type Mn_solvent.

Now tell the interface which concentrations to use as input for the rate expressions. Select the pellet concentrations. The entries will at this stage appear yellow since the feature is not yet created.

In the window, click .

In the window for , locate the section.

Find the subsection. In the table, enter the following settings:


Species	Species concentration type	Molar concentration (mol/m^3)
A	Variable	tds.rpb1.cpe_cA
B	Variable	tds.rpb1.cpe_cB
C	Variable	tds.rpb1.cpe_cC
H2O	Solvent	C_solvent

Select the check box because the is coupled to the feature which is defined in extra dimension.

Click to expand the section. Select the check box.

Continue with the interface to set up the mass transport model.

Transport of Diluted Species in Porous Media (tds)

Fluid 1

In the window, expand the node, then click .

In the window for , locate the section.

From the u list, choose .

Locate the section. In the DF,cA text field, type DA.

In the DF,cB text field, type DB.

In the DF,cC text field, type DC.

From the list, choose .

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the εp list, choose . In the associated text field, type epsilon_b.

Reactive Pellet Bed 1

Add the feature. A predefined extra dimension is attached to this feature. The extra dimension is 1D on the radial coordinate of the pellet particle of which the radius is normalized to 1. The mesh for the extra dimension has a default of 10 elements with a cubic root sequence distribution.

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the rpe text field, type r_pe.

In the εpe text field, type epsilon_pe.

Use a film theory condition to account for any film resistance to mass transfer between the bulk fluid and the pellet. Use spherical pellets.

Locate the section. From the list, choose .

Locate the section. In the Nelem text field, type 6.

Enter the user-defined diffusion coefficients.

Porous Media Diffusion 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

In the Dpeff,cA text field, type DAp.

In the Dpeff,cB text field, type DBp.

In the Dpeff,cC text field, type DCp.

Use the reaction rates calculated in the interface.

Reactions 1

In the window, click .

In the window for , locate the section.

From the Rpe,cA list, choose .

From the Rpe,cB list, choose .

From the Rpe,cC list, choose .

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

Inflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the c0,cA text field, type CA_in.

In the c0,cB text field, type CB_in.

In the c0,cC text field, type CC_in.

Outflow 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Darcy’s Law (dl)

Lastly, enter the model specifications for the interface to compute the convective flow in the reactor.

Fluid and Matrix Properties 1

Take the porosity directly from the feature.

In the window, under click .

In the window for , locate the section.

From the εp list, choose . In the associated text field, type epsilon_b.

From the κ list, choose . In the associated text field, type kappa.

Pressure 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Pressure 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the p0 text field, type p_Darcy.

This completes the setup of the model equations describing the reacting flow and heat transfer in the packed bed reactor. Before solving the problem numerically, the geometry needs to be meshed.

First create a free triangular mesh at the reactor inlet and sweep that mesh along the x direction (the height) of the reactor.

Mesh 1

Free Triangular 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Swept 1

In the toolbar, click

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type 15.

In the text field, type 5.

Click

Since this is a one-way problem, it can be solved in two steps in order to consume less memory: First solve the interface for the velocity, which is a stationary problem. Then solve the interface with a time dependent study step.

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0,10,180).

Locate the section. In the table, clear the check box for .

Stationary

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, clear the check boxes for and .

Right-click and choose .

In the toolbar, click

Root

Click the

button in the toolbar.

In the dialog box, in the tree, select the check box for the node .

Click .

Create views for plotting different angles of the geometry.

Results

Column view

In the window, under right-click and choose .

In the window for , type Column view in the text field.

Pellet view

Right-click and choose .

In the window for , type Pellet view in the text field.

Sector 3D 1

Create a dataset that can be used to plot the column with a sector cut-out for better view.

In the toolbar, click

and choose .

In the window for , locate the section.

In row , set to 1 and to 0.

Locate the section. In the text field, type 8.

From the list, choose .

In the text field, type 5.

Adjust the view angle of the plot with the mouse, then go to the -> under and select the check-box to save the view.

First create Figure 3 showing the velocity distribution in the reactor.

Velocity

In the toolbar, click

In the window for , type Velocity in the text field.

Locate the section. From the list, choose .

Slice 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose .

Locate the section. In the text field, type 8.

In the toolbar, click

Click the

button in the toolbar.

Continue with Figure 4 illustrating the concentration of species A in the reactor.

Concentration, A, Surface (tds)

In the window, click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the toolbar, click

Click the

button in the toolbar.

To create the first xy plot, a new dataset is required. Also, the special syntax comp1.atxd3 is needed to plot variables within the pellet.

Study 1/Solution 1 (3) (sol1)

In the window, under right-click and choose .

In the window for , locate the section.

From the list, choose .

1D Plot Group 10

In the toolbar, click

and choose .

Line Graph 1

Right-click and choose .

Pellet x-y plot

In the window for , type Pellet x-y plot in the text field.

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

In the text area, type Pellet concentrations @ x=0.5[m], y=0, z=0.

Locate the section. From the list, choose .

In the window, under click .

In the window for , type A in the text field.

Locate the section. From the list, choose .

Locate the section. In the text field, type comp1.atxd3(0.5,0,0,comp1.tds.rpb1.cpe_cA).

The syntax comp1.atxd3(0.5,0,0,comp1.tds.rpb1.cpe_cA) means that you visualize the internal pellet concentration 0.5 m from the inlet, and in the center of the column.

Select the check box.

In the associated text field, type pellet (comp1.tds.rpb1.cpe_cA) @ x=0.5[m], y=0, z=0.