EFFECTS OF PHYSICOCHEMICAL PARAMETERS ON THE OPTIMIZED TOPOLOGY OF CATALYST IN MICRO REACTORS

Fluid flow, mass transport and catalytic reaction in micro-reactors with catalyst are studied. Optimized structures of the catalyst generating the maximum reaction rate are explored by numerically solving the reactive transport processes combined with topology optimization schemes in COMSOL 5.3a. Effects of several physicochemical parameters including pressure drop, reaction rate constant, permeability of the porous catalyst and radius of the reactor are investigated. Distributions of the optimized structures, velocity, pressure and concentration are discussed in detail. Results show that the four physicochemical parameters have significant effects on the optimized structures of the catalyst obtained, indicating that generally there does not exist a single structure that can result in the maximum averaged reaction rate under different conditions. Three indicators are evaluated, including the averaged reaction rate, conversion ratio and the volume of catalyst. It is found that change trends of different indicators with different physicochemical parameters sometimes are opposite, indicating great challenging for multiple objective optimization.


INTRODUCTION
Reactive transport processes in porous media are widely encountered in many scientific and engineering problems, such as fuel cells, batteries, filters, packed-bed reactors, micro-fluidic reactors, and so on (Chen et al., 2019;Chen et al., 2018a;Krishna and Wesselingh, 1997;Pardo et al., 2014;Sunden and Yuan, 2010;Thundil Karuppa Raj and Thulasi, 2018). There exists strong interaction between the complex structures of porous media and the multiple sub-processes involved including fluid flow, mass transfer and chemical reactions. Optimizing the porous structures and porous media distributions to improve the performance of related systems post a great challenge.
With the rapid development of theory and technology, improving energy efficiency and reducing energy consumption (Sun et al., 2018) are continuously required for applications of the mass transfer process in not only theoretical research but also engineering practice (Sun et al., 2019). To this end, on the one hand, the reactant and catalyst with excellent reactivity and low cost should be identified and fabricated. On the other hand, deep understanding of multiple reactive transport processes in porous media and effects of porous structures are of great importance for improving the performance of porous media and reducing the cost. For instance, Chen et al. (2018a) studied the diffusion reaction process of oxygen in reconstructive catalyst layer and found out that the micro-structure significantly affected the total reaction rate. Okkels and Bruus (2007) found that using the optimized distribution of active porous material could increase the reaction rate apparently. Consequently, it is of vital significance to optimize structures of the equipment which can satisfy the demands of the equipment.
Recognizing the significance of the structure of the key parts in facilities, plenty of efforts have been devoted to gain the optimal structure. Conventional design methods concerned more about changing the geometric parameters or shapes on the basis of the original structures. * Corresponding author. Email: lichennht08@mail.xjtu.edu.cn Nevertheless, the artificial design is so limited that it is not clear whether some changes of simple parameters applied on that basis, which also can be regarded as a restriction, can lead to the best performance for the devices under certain conditions. Therefore, topology optimization methods emerge as required.
Topology optimization is a mathematical method that optimizes the material layout within a given design space, for a given set of loads, boundary conditions and constraints with the goal of maximizing the performance of the system , also can be treated as a kind of inverse design. Compared with size and shape optimization, the biggest advantage of topology optimization is that no presupposition needs to be taken on the geometry. Instead, only design domain, boundary conditions and optimizing objective are required to be configured (Andreasen et al., 2009). That means the structure in the specific domain is completely unknown before optimization, which indicates that there is no restriction of prescribed shape and size. Therefore, it can get some structures which are strange or unforeseeable but optimal.
Since 1988, topology optimization have been widely adopted to solve structural optimization problems (Bendsøe, 1988). Initially, studies on topology optimization are usually concentrated on structural mechanics problems. Suzuki and Kikuchi (1994) verified the realization method and specific application conditions of topology optimization. They found that the extended method can solve any type of linearly elastic structures. Sigmund (2001) developed a 99 line code of topology optimization in MATLAB to obtain the maximum stiffness. Subsequently, topology optimization has been widely studied and extended to other fields such as fluid flow field (Borrvall and Petersson, 2003;Gersborg-Hansen et al., 2005), heat transfer field (Bruns, 2007;Gersborg-Hansen et al., 2006;Iga et al., 2009), coupled field (Alexandersen et al., 2014;Yoon, 2010), and so on. For optimization problem including mass transport with chemical reactions, there are some studies in the literature using topology optimization method. For instance, Okkels and Bruus (2007) optimized distributions of porous catalyst in micro-reactors to obtain higher reaction rate using topology optimization. Yaji et al. (2017) conducted 2D optimization design of negative electrode flow fields in vanadium redox flow batteries to improve its performance. Behrou et al. (2019) designed flow fields of proton exchange membrane fuel cells, and compared the performance of the topology-optimized, parallel and serpentine flow fields. Deng et al. (2018) provided a topology design of the electrode-pair in electroosmotic micro-mixer in order to achieve the complete mixing of two fluids. Hence, topology optimization can facilitate to design and optimize the structure of microfluidic devices.
The present study aims to enhance reactive transport processes in micro reactors. In practice, the micro reactor with different size may be operated under different conditions. Therefore, in the present study, based on the work of Okkels and Bruus (2007), effects of several parameters are further studied including pressure drop, reaction rate constant, permeability of the porous catalyst and radius of the reactor. In Section 2, the physicochemical problem studied are introduced. In Section 3, the topology optimization problem is presented. Then in Section 4, effects of different factors in the optimized structures are discussed in detail. Finally, in Section 5, some important conclusions are drawn.

PHYSICOCHEMICAL MODEL
Following the work of Okkels and Bruus (2007), the computational domain for the micro-reactor is a rectangle with an inlet and an outlet as shown in Fig. 1. Due to the symmetry, only the top half of the domain is investigated. The reactive transport processes taking place inside the domain can be described as follows. The bulk fluid with a dilute species A flows into the domain through the left inlet Гin. Then chemical reaction, which consumes the dilute species, occurs in the reaction domain Ω which is the domain to be designed, by catalytic action when the species comes into contact with the porous catalyst. Eventually, the surplus reactant flows out of the domain through the right outlet Гout. The distribution of the catalyst, or the catalyst topology in domain Ω, is required to be designed to meet objective functions, such as maximum reaction rate.

Fig. 1 Geometry model of micro-reactor
From the above description, it can be found multiple processes, including fluid flow in both clear region and porous region, mass transport and chemical reaction simultaneously take place in the domain. In the present study, numerical studies will be conducted to investigate above reactive transport processes and to optimize the catalyst distributions. Before establishing the mathematical models, the following assumptions are employed: • The bulk fluid is assuming as an incompressible liquid with fixed density ρ and viscosity μ; • A steady-state, laminar and isothermal flow is hypothesized in the micro-reactor.

•
The consumption of the dilute species has no effect on the fluid flow. In other words, the fluid flow process and mass transport process are one-way coupled.

Governing equations and boundary conditions
The governing equations for the reactive transport processes includes mass, momentum and concentration conservation equations. They are briefly introduced as follows.
The mass conservation equation for steady state incompressible flow is as follows: The momentum equation considering the existence of the porous catalyst is as follows (Yaji et al., 2017): where ε, the design variable of the topology optimization problem, denotes the existence of the porous catalyst. Its values as unity indicate the clear fluid flow region and zero for the porous catalyst region, and that between zero and unity means the local computational mesh is partially occupied by the catalyst. The last term on the right hand side of Eq.
(2) is the Brinkman friction term, which represents the effects of the porous catalyst on the bulk fluid flow. In this term, α(ε) is the inverse permeability expressed as follows (Olesen et al., 2006): where q is a dimensionless parameter. The term ( ) help to make the boundary smoother between catalyst region and the clear flow region. To better understand Eq. (3), the relationship between Fig. 2 under different q. It can be found that for infinite q, ( ) (3) is the maximum inverse permeability, which is related to the Darcy number Da as follows (Olesen et al., 2006): where L is the characteristic length of the computational domain. The relationship between the permeability of the porous catalyst k and αmax can be easily obtained based on Eqs. (4-5): Obviously, based on Eqs. (3-6), it can be found that when ε is equal to 1, which means there is no catalyst, the Brinkman term is equal to 0. Conversely, if ε is equal to 0, which indicates the local region is occupied by the porous catalyst, α(ε) is equal to αmax. The governing equation for the concentration of the dilute species A is as follows (Yaji et al., 2017): where c is the concentration. D(ε) is the diffusivity, which is calculated by the following formula: where D is the binary diffusivity of the species in the clear fluid region, while Deff is the effective diffusivity in the porous catalyst. The latter one is related to the porosity e and tortuosity τ of the porous catalyst: In the present study, the Bruggeman equation is adopted to determine Deff, where τ is taken as e -0.5 . The volumetric source term ϕ(ε) in Eq. (7) is employed to represent the local reaction rate in the porous catalyst Eq. (10): where ka refers to reaction rate constant (Chen et al., 2018b;Kim and Sun, 2012). The boundary conditions are as follows. For the fluid flow, pressures are prescribed at the inlet and outlet of the reactor, with nonslip conditions at the fluid-solid boundaries. For the mass transport, concentration is given at the inlet while a fully developed boundary condition is adopted at the outlet.

Statistical indicators
In this work, three statistical indicators are employed to measure the performance of the micro-reactor with different distributions of the porous catalyst. The first one is the conversion ratio of reactant X: where Fin and Fout refer to the flow of reactant with concentration c through inlet and outlet, respectively. And they are defined as follows: where nin and nout are the normal direction of the boundary Гi and Гout, respectively. The second indicator adopted is the total catalyst volume fraction Ctotal, which is defined as follows: where s is the area of the design domain.
The third indicator is the averaged reaction rate φ , which is defined as follows (Okkels and Bruus, 2007):

TOPOLOGY OPTIMIZATION CONFIGURATION
In the present study, the distribution of the catalyst inside the design domain Ω is optimized to maximize the third indicator, namely the averaged reaction rate. The topology optimization problem is described as follows: It is worth mentioning that the problem defined by Eq. (16) is an optimization problem. This is because if the entire domain is free of the catalyst, φ will be zero. On the other hand, if the domain is full of the catalyst, the flow resistance is extremely high, leading to a small amount of the dilute species supplied into the domain, and thus φ will be low.
Therefore, there exists an optimal distribution of the catalyst, which leads to the maximum φ .The above topology optimization problem is implemented in the software COMSOL Multiphysics 5.3a. In COMSOL, the "Laminar Flow" module is configured to solve the mass and momentum governing equations, while the "Transport of Diluted Species" module is adopted to solve the mass transport problem of the dilute species.
Particularly, the "Optimization" module is activated in which the design variable and optimizing objective for the topology optimization can be defined. The optimization algorithm called SNOPT is adopted for the optimization process, which is especially effective for nonlinear problems in which functions and gradients are expensive to evaluate. For details of the SNOPT, one can find in the work of (Gill et al., 2005).

RESULTS AND DISCUSSION
In this section, first the fluid flow and mass transport in the optimized structures of the base case are analyzed in detail. Then, the effects of several factors including pressure drop, reaction rate, and permeability of the porous catalyst and the size of the domain are studied in detail. Values of the physiochemical variables for the base case are listed in Table 1.
First of all, the grid independence is conducted for the base case (listed in Table 1). Grid number from 2500 to 75000 is studied, and two important variables, namely conversion ratio as defined in Eq. (11) and maximum reaction rate as defined in Eq. (15) are evaluated. As shown in Fig. 3, when the grid number increases from 2500 to 15000, the maximum average reaction rate rises rapidly while the conversion ratio drops dramatically. After the grid number of 16000, the change of two variables becomes slightly. The relative error between the grid number of 16346 and 73789 is 0.2% and 0.11% for the conversation ratio and maximum reaction rate, respectively. For the balance between computational resources and simulation accuracy, the grid number of 16346 is chosen.
To validate the proposed approach, a comparison between the work by Okkels and Bruus (2007) and this work is shown in Fig. 4. In this verification, the computational domain, the initial and boundary conditions, as well as the values of all the variables are the same. In detail, the Darcy number is 1×10 -4 , the binary diffusivity is 3×10 -8 m 2 ·s -1 , and the reaction rate constant is 1 s -1 . Moreover, the pressure drop between inlet and outlet is 0.25 Pa. It can be noticed that the distributions of catalyst in these two works are basically similar with some differences in details. The reason for these differences may be the different computational method. In the study of Okkels and Bruus (2007), MATLAB and COMSOL were combined to solve this optimization problem, while in the present study only COMSOL with "Optimization" module is applied for calculating.

Base case
Fig. 5 displays the distributions of catalyst, reactant concentration, velocity and pressure in the optimized structure for the base case. The structure is described as a gray-scale image, in which the color of each mesh stands for the value of volume fraction of pore. Black and white parts represent the porous catalyst and fluid region, respectively. It can be found that the optimized distributions of catalyst (black regions) are quite complex and non-intuitive. The optimized structure has two large pores, with one connecting the inlet and the outlet while the other is a dead-end pore. Clear observation found that there are slightly coarse boundaries between pores and the porous catalyst, which will result in difficulties for the subsequent processing and manufacturing.  Correspondingly, the velocity and pressure fields show complex distributions, which are significantly affected by the local catalyst distribution. The velocity in the porous region is relatively low due to the local low permeability. For the concentration of species A, it can be found that the concentration decreases rapidly at the fluid-catalyst interface due to the chemical reaction. Inside the region of the catalyst, the local concentration can be as low as zero. The existence of the pore penetrating the domain serves as a highway for mass transport, and the reactant is not completely consumed at the outlet. Finally, for the based case, the three statistical indicators are calculated. The maximum average reaction rate, conversion ratio and the total catalyst volume fraction are 0.04 mol·m -3 ·s -1 ,74.26% and 35.99%, respectively.

Effects of the pressure drop
In the practical application, the pressure drop between the inlet and outlet can be adjusted to improve the reactor performance. In this section, distributions of the catalyst, concentration, velocity and pressure under different pressure drop ∆p of 0.1 Pa, 0.25 Pa, 1 Pa and 2 Pa are investigated and displayed in Fig. 6(a-d), respectively. As can be seen from Fig. 6, different pressure drops lead to different distributions of catalyst in the reaction domain. For an extremely low-pressure drop, the dominant transport mechanism for mass transport is diffusion. As pressure drop increases, the pores penetrating the reactor that tend to diminish. For ∆p as 0.1 Pa, there are two connected pores from the inlet to the outlet, and the geometry of the catalyst region is almost symmetry along the x direction. As ∆p increases, the catalyst region expands, causing the close of the bottom connected pore and the shrink of the upper connected pore. For ∆p greater than 1 Pa, it can be found that the originally connected pore has been completely blocked. In addition, as ∆p increases, the velocity increases, and the velocity and pressure fields develop according to the catalyst distribution.
For the mass transport, for ∆p as 0.25 Pa, species transport through the two connected pores is quickly depleted inside the catalyst. As the connected pores are closed or become narrow, the species has to penetrate into the catalyst through two narrow slots, as shown in Fig. 6(c). Note that although the catalyst topology is almost the same in Fig. 6(c) and 6(d). The concentration contour with a higher concentration in Fig.  6(d) is larger than that in Fig. 6(c). This is because a higher pressure gradient facilitates the species to transport deeply into the catalyst. Further, it can be found that for higher pressure, the species is not completely assumed in the reactor and a considerable amount of species flows out of the domain. To better understand the effects of pressure drop, Fig. 7 depicts the conversion ratio, maximum average reaction rate and total catalyst volume fraction under different pressure drop. It can be found that as pressure drop increases, the amount of the catalyst with the optimized topology also increases. The change trend of the conversion ratio and the maximum average reaction rate is reverse. It can be found that as pressure drop increases, the conversion ratio decreases, because velocity in the domain increases which reduces the residence time of the species contacting with the catalyst. This is in agreement with the increasing concentration at the outlet under a higher pressure drop shown in Fig. 6. The maximum average reaction rate increases as the pressure drop increases, as more species are supplied into the domain for chemical reaction and the region with higher concentration expands.

Fig. 7 Different statistical indicators under different pressure drop between inlet and outlet.
From the above discussion, it can be found that as the pressure drop increases, the change treads of the three indicators are not consistent. For some processes, higher conversion ratio is required, while under some scenarios the maximum average reaction rate is desirable. In addition, the increase of the catalyst layer volume fraction also leads to higher cost. In the present study, only the maximum average reaction rate is adopted as the optimized objective. For further optimization of multiple objectives, the optimization will be much more complicated, and it is highly possible that the requirement of the multiple objectives cannot be satisfied simultaneously.
Besides, it also can be found that the optimized structures are different under the operating parameters, such as the pressure drop. This indicates that in practice, once the optimized structure has been determined, it is desirable to operate the reactor under the designed operating parameters, and otherwise the reactor performance will be greatly deteriorated.

Effects of the reaction rate constant
In practice, the reaction rate constant will vary under different operating condition, such as different temperature, or different over-potential for electrochemical reaction. In this section, the effects of the reaction rate constant are studied. The distributions of catalyst, concentration, flow rate and pressure drop under different rate constant of 0.025 s -1 , 0.25 s -1 , 1 s -1 and 10 s -1 are presented in Fig. 8.
It can be found that the optimized structures are significantly affected by the reaction rate，showing remarkable alternation of the morphology. As reaction rate constant increases, the volume of the porous catalyst is greatly reduced, and ε becomes higher as indicated by the gray scale image. For the concentration distribution, it can be found that for a low reaction rate constant, even after transporting through the catalyst region, the species is consumed not much, leading to higher concentration even at the reactor outlet, as shown in Fig. 8(a). While for relatively high reaction rate constant, the species is completely consumed after contacting the catalyst.   9 shows the impacts of reaction rate constant on the three indicators. It can be found that both the conversion ratio and the maximum reaction rate increase as the reaction rate constant increases, and the increase is sharp before 5 s -1 . Due to the high chemical reaction rate, the conversion ratio can approach one hundred percent with an extremely low volume of catalyst.

Effects of the permeability
Permeability is an important parameter that describes the capacity of the porous media for fluid flow through. In this section, the effects of permeability of the porous catalyst on the reactive transport processes are studied. Fig. 10 depicts the distributions of catalyst, concentration, flow rate and pressure drop under different permeability of 10 -11 m 2 , 10 -10 m 2 , 10 -9 m 2 and 2×10 -9 m 2 . From Fig. 10, remarkable alternations of the optimized structures of the porous catalyst can be observed. As permeability increases, the catalyst volume goes up and the catalyst tends to fill the entire reactor. For example, when the permeability is 10 -11 m 2 , there are two connecting pores from the inlet to the outlet. But when the permeability is higher than 10 -9 m 2 , there is one connected pore which is partially filled with the catalyst. For the fluid flow, under extremely low permeability, the catalyst region is almost not available to the fluid flow, and thus the velocity in the catalyst region is almost zero. As the permeability increases, the velocity in the porous catalyst region increases. For the case with permeability as high as 2×10 -9 m 2 as shown in Fig. 10(d), the velocity in the catalyst region is quite close to that in the pores.
For the concentration distribution, as the catalyst with low permeability is not available to the fluid flow and mass transport, only the pore-catalyst interface is allowed for chemical reaction, leading to less reactant consumed and thus higher concentration at the domain outlet ( Fig. 10(a)). As the permeability increases, the reactant penetrates into the catalyst region where the reaction takes place, and thus more reactant participates the reaction. Fig. 11 shows the influences of the permeability on the three statistical indicators. As shown in Fig. 11, all three indicators grow up with the permeability. For the conversion ratio, higher permeability of catalyst leads to the increased contact area between catalyst and reactant, leading to a higher conversion ratio. The volume fraction of the catalyst increases, agreeing with that in Fig. 10.  Values of other parameters are the same as the base case. As can be observed from Fig. 12, as the radius increases, the catalyst topology changes significantly. Generally, as the radius increases, more branches of the porous catalyst are generated, leading to a complicated flow network inside the reactor. For the concentration, it can be seen that the concentration at the outlet decreases as the radius increases. This is because larger reactor domain reduces the flow rate and thus prolongs the residence time, which facilitates the reaction of the dilute species. This is supported by the conversion ratio curve in Fig. 13. Intuitively, the larger size of the reactor with higher radius results in more reaction area, and thus might lead to higher maximum reaction rate. However, the fact is that the maximum average reaction rate greatly reduces, from 0.04 mol·m -3 ·s -1 to 0.025 mol·m -3 ·s -1 , while the absolute volume of the catalyst is greatly increased. This result is undesirable, which indicates a huge waste of the catalyst. The undesirable result is caused by the slower flow rate in the larger domain, which decelerates the transport process.

CONCLUSIONS
In the present study, reactive transport processes inside a 2D microreactor are numerically investigated. Topology optimization method is adopted to find the optimized structures of porous catalyst inside the micro-reactor which can generate the maximum reaction rate. Effects of several factors on the optimized structures, pressure, velocity and concentration are explored. The main conclusions are as follows (1) Under different conditions, the optimized structures of porous catalyst are different.
(2) As pressure drop increases, the maximum averaged reaction rate and the total catalyst volume fraction increase while the conversion ratio drops.
(3) As the reaction rate increases, all the indicators can be optimized.
(4) As permeability of the catalyst rises, all the three indicators increase, which means that only the total catalyst consumption is not optimized.
(5) When the radius of reaction domain rises, the conversion ratio and the total catalyst volume fraction can be optimized slightly but the maximum average reaction rate becomes lower.
It is worth mentioning that the boundaries obtained from the simulation are not very smooth. Such rough boundaries of the optimized structures will pose great challenging for practical manufacturing. Therefore, the problem that how to generate boundaries clearer and avoid checkerboard phenomenon at the same time should be addressed.