COMSOL 6.4 - Optimization of a Catalytic Microreactor

Optimization of a Catalytic Microreactor

In this model, a solution is pumped through a catalytic bed, where a reactant undergoes chemical reaction as it gets in contact with the catalyst. The purpose of the example is to maximize the total reaction rate for a given total pressure difference across the bed. This is achieved by finding an optimal catalyst distribution. The distribution of the porous catalyst determines the total reaction rate in the bed. A large amount of catalyst results in a low flow rate through the bed, while less catalyst gives a high flow rate but low conversion of the reactant.

This modeling example is based on Ref. 1.

This application requires the Optimization Module.

Model Definition

The model geometry is shown in Figure 1. The reactor consists of an inlet channel, a fixed catalytic bed, and an outlet channel.

Figure 1: Model geometry.

The optimal catalyst distribution should maximize the average reaction rate, which is expressed as the integral of the local reaction rate, r (SI unit: mol/(m3·s)), over the domain, Ω:

Assuming a first-order catalytic reaction with respect to the reactant species, the local reaction rate is determined by

(1)

where θ denotes the volume fraction of solid catalyst, c refers to the concentration (SI unit: mol/m3), and ka is the rate constant (SI unit: 1/s).

The mass transport is described by the convection and diffusion equation

where u denotes the velocity vector (SI unit: m/s) and D is the diffusion coefficient (SI unit: m2/s). The Navier–Stokes equations describe the fluid flow:

(2)

The coefficient α(ε) depends on the distribution of the porous catalyst as

(3)

where Da is the Darcy number; L is the length scale (SI unit: m); and q is a dimensionless parameter, the interpretation of which is discussed in the next section.

From Equation 3, the direct conclusion is that when ε equals 1, α equals zero and Equation 2 reduces to the ordinary Navier–Stokes equations. In this case the reaction rate is zero; see Equation 1.

To summarize, the optimization problem is

(4)

where

and physical boundary conditions apply.

Convex Optimization Problems

One of the most important characteristics of an optimization problem is whether or not the problem is convex. This section therefore briefly describes this property. For a more general discussion of the subject, see for example Ref. 2.

A set C is said to be convex if for any two members x, y of C, the following relation holds:

that is, the straight line between x and y is fully contained in C. A convex function is a mapping f from a convex set C such that for every two members x, y of C

(5)

An optimization problem is said to be convex if the following conditions are met:

•

the design domain is convex

•

the objective and constraints are convex functions

The importance of convexity follows simply from the result that if x* is a local minimum to a convex optimization problem, then x* is also a global minimum. This is easily proven by simply assuming that there is a y such that f(y) < f(x*), and then using Equation 5.

This particular optimization problem is nonlinear, because a change in ε implies a change in the concentration, c. Because of this implicit dependence, it is very difficult to determine whether or not the objective is convex. There is therefore no guarantee that the optimal solution you obtain is globally optimal or unique. In the best of cases, running the optimization gives a good local optimum.

The parameter q can be used to smoothen the interfaces between the catalyst and the open channel. To see the effect of this parameter, rewrite Equation 3 as

It follows that when q approaches infinity, α is the (inverse) porosity. On the other hand, lowering the value of q decreases the magnitude of α.

Figure 2: q(1 − ε)/(q − ε) plotted as a function of ε for different values of q.

Figure 2 shows q(1 − ε)/(q − ε) plotted as a function of ε for different values of q. This plot shows that lowering the value of q, increases the convexity of the force coefficient. For a low q value, an increase in ε around 0.5, imposes a small increase of the force coefficient, while for a higher value of q, a change in ε imposes an almost equal change for the whole range. Therefore, for a lower q value, the solution is not sharp at the interfaces. On the other hand, for small values of ε, the force term decreases rapidly when q is small, and thus affects the flow field to a much wider extent. In the limit when q approaches infinity, α as a function of ε is a straight line.

Results and Discussion

Figure 3 shows the velocity field in the empty channel. This is the starting point for the optimization.

Figure 3: Velocity field in the open channel.

Figure 4: Distribution of the porous catalyst seen in black and open channel in white.

Figure 4 shows the distribution of the porous catalyst in black and the open channels in white. This result shows that, optimally, the supply of the reactant should be distributed over a large area of the reactor. Note also that the amount of open channel volume is significant.

Figure 5 shows the concentration distribution in the reactor. This plot shows how the porous catalyst is fed with the reactant through the open channels. The plot naturally resembles that of Figure 4.

Figure 5: Concentration distribution in the reactor after optimization.

Let

where nflow refers to the normal to the boundary ∂Ωi in the flow direction (that is, pointing in to the domain at the inlet and out from the domain at the outlet). Then Fi is a measurement of the flow of the species with concentration c through the boundary ∂Ωi per unit length in the transverse dimension. The conversion, X, of the reactant is defined as

In this case, the conversion of the reactant is around 50%.

Figure 6 shows the velocity field in the reactor. The porous catalyst slows down the flow significantly compared to Figure 3.

Figure 6: Velocity field in the reactor after optimization.

References

1. F. Okkels and H. Bruus, “Scaling Behavior of Optimally Structured Catalytic Microfluidic Reactors,” Phys. Rev. E, vol. 75, pp. 016301 1–4, 2007.

2. S.G. Nash and A. Sofer, Linear and Nonlinear Programming, McGraw Hill, 1995.

Chemical_Reaction_Engineering_Module/Reactors_with_Mass_Transfer/microreactor_optimization

Modeling Instructions

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 .

Click

In the tree, select > .

Click

Root

Load 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 microreactor_optimization_parameters.txt.

Geometry 1

Next, create the geometry. The reactor consists of three domains: the inlet channel, the reacting domain, and the outlet channel (see Figure 1).

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2*L.

In the text field, type L.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 6*L.

In the text field, type 3*L.

Locate the section. In the text field, type 2*L.

Locate the section. Select the checkbox.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 2*L.

In the text field, type L.

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

Click

Click the

button in the toolbar.

The geometry should now look like that in Figure 1.

Define integration couplings to use for calculating the conversion of the reactant.

Add a density topology feature, which can be used to distinguish between free flow and solid regions. This variable will be coupled back to the interfaces later.

Component 1 (comp1)

Density Model 1 (dtopo1)

In the toolbar, click

and choose .

Only the center part of the channel geometry is needed in the optimization, so you only have to define the feature there.

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

In the q text field, type q.

Locate the section. In the θ0 text field, type 1.

Now the design variable used in the optimization is defined. The initial value 1 corresponds to a channel free from porous material.

Definitions

Variables 1

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Select Domain 2 only.

Locate the section. In the table, enter the following settings:


Name	Expression	Unit	Description
phi	k_a(1-dtopo1.theta)c	mol/(m³·s)	Local reaction rate
alpha	(mat1.def.eta(minput.T)/(DaL^2))dtopo1.theta_p		Drag-force coefficient

Objective Function

In the toolbar, click

and choose .

In the window for , type Objective Function in the text field.

In the text field, type obj.

Locate the section. From the list, choose .

Click in the upper-right corner of the section. From the menu, choose > > > .

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Transport of Diluted Species (tds)

Click the

button in the toolbar.

In the dialog, select > in the tree.

In the tree, select the checkbox for the node > .

Click .

In the window, expand the > node, then click .

In the window for , click to expand the section.

From the list, choose , so that the adjoint problem converges without warnings (this step is optional). Alternatively, the warning can be avoided by disabling compensation for nojac terms on the .

Laminar Flow (spf)

In the window, under click .

In the window for , click to expand the section.

Find the subsection. From the list, choose to save computational time.

Volume Force 1

In the toolbar, click