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, Ω. This is equivalent to minimizing the negative of this average reaction rate:

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 check box.

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.

Definitions

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 1 only.

Integration 2 (intop2)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 12 only.

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.

In the toolbar, click

and choose .

Topology Optimization

Density Model 1 (dtopo1)

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 qDarcy 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

Use the integration operators to define the conversion of reactant and the scaled objective function.

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
F_in	intop1(tds.tflux_cx)	mol/(m·s)	Molar flow at inlet
F_out	intop2(tds.tflux_cx)	mol/(m·s)	Molar flow at outlet
X_conv	(F_in-F_out)/F_in		Conversion of reactant

Here, tds.tfluxx_c is the COMSOL Multiphysics variable for the x-component of the total flux.

Variables 2

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	(eta/(DaL^2))dtopo1.theta_p	Pa·s/m²	Drag-force coefficient