Florian Dugast

Yann Favennec

Christophe Josset
email: christophe.josset@univ-nantes.fr

Reactive uid ow topology optimization with the multi-relaxation time lattice Boltzmann method and a level-set function

Keywords: Topology optimization, adjoint-state formulation, reactive uid ow, lattice Boltzmann method, level-set method

Reactive fluid flow topology optimization with the multi-relaxation time lattice Boltzmann method and a level-set function

Florian Dugast, Yann Favennec, Christophe

Introduction

The development of ecient reactors is needed in many chemical process and energy applications: chemical synthesis, catalytic combustion, and electro-chemical conversion, to cite but a few. The growth of renewable energy (wind, solar, or marine) implies an improvement of the electricity storage technologies. Electrolysers, or redox-ow batteries, can be used for the rst stage of the electro-chemical conversion, and fuel cells for the reverse reactions. One of the principal characteristics to improve is the power density (kW/m 3 , kW/m 2 or kW/kg) that represents the "good use" of the constitutive materials. All these systems can be seen as electrochemical reactors, where reactants are supplied to the reaction place by convection and diusion mechanisms while by-products (thermal energy or new species) need to be removed (to avoid thermal problems or to leave free access to the reaction location). In the case of conventional vanadium redox-ow batteries, the reaction occurs on the material surface, whereas in the case of fuel cells, the species transport is more complicated. For example, in the case of polymer exchange membrane fuel cells (PEMFC), the reactants are usually transported into the gas channels and are diused through a thin porous medium, in order to reach the catalyst layers (anode or cathode).

In this context, the eciency of the reaction depends on the uid ow channels geometry inside the domain. Some geometries have been tested in literature to this matter (serpentine [START_REF] Messaggi | Analysis of ow eld design on vanadium redox ow battery performance: Development of 3D computational uid dynamic model and experimental validation[END_REF], inter-digitated [START_REF] Houser | Architecture for improved mass transport and system performance in redox ow batteries[END_REF], etc.), but there are still discussions about which conguration is optimal, or the best, and benchmark parameters are not well-established [START_REF] Cervantes-Alcalá | Flow distribution and mass transport analysis in cell geometries for redox ow batteries through computational uid dynamics[END_REF]. Structural optimization techniques provide a promising method to overcome the strong limitation of determining a priori a system's geometry. Previous methods required this to be specied to allow optimization through testing relevant parameters (e.g. channel dimensions), signicantly limiting the structural congurations possible.

Doing so, the uid/solid distribution is to be determined via optimization methods. Among the three main optimization methods, namely the size optimization, the shape optimization, and the topology optimization [START_REF] Sigmund | Topology optimization approaches[END_REF], the latest is associated to the highest number of degrees of freedom. With topology optimization, very dierent geometries can be dealt with and obtained, leading to high levels of geometrical complexity, and with the ability of creating holes inside the domain of interest (this is not possible with either optimization strategies of size and shape optimization).

The basic idea of topology optimization is to optimize a material allocation problem where the material can have dierent properties. One of the rst applications concerns the structural mechanical optimization, with compliance minimization [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF][START_REF] Sigmund | Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, meshdependencies and local minima[END_REF]. In that case, the material is either void or solid, but it can also be uid or solid, in order to build uid channels, for example.

Topology optimization is well developed for uid ow designs [710], but most often for solving pressure drop minimization problems.

One challenge of topology optimization is to develop multi-physics applications. As the complexity of the physics increases for these problems, topology optimization can be helpful to suggest innovative geometries that could not be intuitively thought. Although some topology optimization algorithms have been proposed for some heat transfer applications [1115] , the development of such methods for chemical reaction is relatively scarce. Non exhaustively, the following works are of great interest. In their seminal publication, Okkels et al. [START_REF] Okkels | Scaling behavior of optimally structured catalytic microuidic reactors[END_REF] extended the work of Borrvall et al. [START_REF] Borrvall | Topology optimization of uids in Stokes ow[END_REF] in the case of advection-diusion-reaction problem. They worked on a catalytic microreactor for topology optimization and compared their results with empirical geometries, like a uniform or a membrane reactor. They also studied the inuence of the relevant physical parameters on the optimized geometries.

Schäpper et al. [START_REF] Schäpper | Topology optimized microbioreactors[END_REF] dealt with topology optimization of bio reactors. They could show that the production of protein could be improved with topology optimization compared to a conventional reactor with a uniform distribution of biomass. A compromise is found for the glucose owrate, sucient enough for the reaction to happen and limited to avoid a strong ethanol production, which decreases the enzym activity. Yaji et al. [START_REF] Yaji | Topology optimization for the design of ow elds in a redox ow battery[END_REF] worked on vanadium redox ow batteries to provide a sucient supply of species to the carbon ber electrode, under a pressure drop constraint. The electrode is modeled as a porous medium with Darcy's law and the generation rate depends on the local velocity. The resulting optimal structure is an interdigitated ow eld with strong tortuosity. This is due to the porous structure of the solid electrode, encouraging the ow through it in order to enhance the reaction. Kim et al. [START_REF] Kim | Topology optimization of gas ow channel routes in an automotive fuel cell[END_REF] worked on PEM fuel cells optimization in order to maximize the reaction rate between hydrogen and oxygen. Their topology optimization method has been applied to dierent geometry initializations such as centered inlet and outlet, multi-terminal inlet and outlet or serpentine channels.

In these three studies: [1618], the catalyst for the reaction is located directly on the porous support, whereas in [START_REF] Kim | Topology optimization of gas ow channel routes in an automotive fuel cell[END_REF], the reaction takes place in the uid channels. In this last application, the specicity of a PEM fuel cell with dierent layers is taking into account: the optimization deals with the best way to feed reactant into a planar electrode and the solid has two functions: electrode supporting and uid routing. The work presented in this article will follow the same approach. All the previously cited works rely on the SIMP modelization regarding the design material. This one is represented by a scalar eld varying continuously from 0 (solid) to 1 (uid). Intermediate states are penalized during the optimization process in order to promote a clear transition between the solid and the uid. However, gray-scale eects occur, and, as mentioned by Makhija et al. [START_REF] Makhija | Topology optimization of multi-component ows using a multi-relaxation time lattice Boltzmann method[END_REF], the porosity model does not prevent diusion into and through solid regions, excluding low Peclet number simulations, when advection and diusion are of the same order of magnitude. In order to overcome this issue, a level-set function (LSF) was used in this study for the uid/solid representation. As the level-set function is continuous, it provides a continuous uid/solid transition as required by the optimization process. However it is possible to use the zero contour of the LSF in order to get a clear uid/solid interface. This latter feature is particularly interesting to obtain accurate boundaries for the forward problem. Many dierent methods can fall under a level-set appellation. Diering from the original work of Allaire [START_REF] Allaire | A level-set method for shape optimization[END_REF], the evolution of the level-set function will not be realized with the Hamilton-Jacobi equation combined with shape derivatives. A similar approach to that of Yamada [START_REF] Yamada | A topology optimization method based on the level set method incorporating a ctitious interface energy[END_REF] was applied here, with no advection in the LSF evolution equation but only a source term depending on the design sensitivities. With this method, the LSF does not need to be re-initialized during the optimization process and the creation of discontinuity in the material distribution is also possible. A smooth Heaviside function is proposed for the calculation of the cost function gradient, and a threshold is then applied to compute the forward problem, with discrete 0-1 boundaries.

Next, all the works previously cited for the topology optimization of chemical process were based on a nite element strategy to solve both the Navier-Stokes equations and the advection-diusion-reaction equation. As an alternative to ordinary discretization methods (nite element method (FEM), nite volume method (FVM)), the lattice Boltzmann method (LBM) may be very ecient for solving that kind of problem. Originally, this method has been used for uid ow problems [START_REF] Succi | The lattice Boltzmann equation: A new tool for computational uid-dynamics[END_REF][START_REF] Zou | On pressure and velocity ow boundary conditions and bounceback for the lattice Boltzmann BGK model[END_REF] in which both the displacement and the collision of particles are solved via the Boltzmann equation, involving a (often simplied) collision operator. The macroscopic quantities for the ow eld (density, pressure, velocity) can then be recovered by the moments of the distribution function involved in the Boltzmann equation.

Since the pioneering work of Pingen [START_REF] Pingen | Topology optimization of ow domains using the lattice Boltzmann method[END_REF], adjoint-state LBM methods have also been developed for uid ow identication [START_REF] Klemens | Solving uid ow domain identication problems with adjoint lattice boltzmann methods[END_REF][START_REF] Klemens | Cfd-mri: A coupled measurement and simulation approach for accurate uid ow characterisation and domain identication[END_REF] or topology optimization [START_REF] Dugast | Topology optimization of thermal uid ows with an adjoint lattice boltzmann method[END_REF][START_REF] Yaji | Topology optimization using the lattice Boltzmann method incorporating level set boundary expressions[END_REF][START_REF] Liu | Discrete adjoint sensitivity analysis for uid ow topology optimization based on the generalized lattice Boltzmann method[END_REF]. It is also to be noticed that the LBM computational time of both forward and adjoint-state problems can be signicantly reduced with GPU parallel computing [3032]. Further, with the use of a uniform grid, the remeshing process may not be necessary along with the optimization iterations. This is another clear advantage of the LBM over FEM and FVM: it is indeed very convenient in topology optimization where the complex uid/solid geometry evolves with the iterations of the optimization process. Also, the LBM boundary conditions are easy to implement (bounce-back at no-slip walls for example) making the LBM a good candidate if uid/solid porous structures need to be computed, which can occur in topology optimization with the creation of complex geometries. Generally speaking, the LBM is convenient for multiphase ows, as the interface can be implicitly tracked with multiple distribution functions [START_REF] Liu | Multiphase lattice Boltzmann simulations for porous media applications[END_REF]. Though being beyond the scope of this article, such physics will be used in oncoming research.

In order to compute reactive ows, the concentration eld is represented by a passive scalar quantity which is transported by the uid ow. This quantity is computed in the LBM framework by a double distribution function approach [START_REF] Yan | Numerical simulation of heat transfer and uid ow past a rotating isothermal cylinder A LBM approach[END_REF][START_REF] Li | Lattice Boltzmann models for the convection-diusion equation: D2Q5 vs D2Q9[END_REF], more stable than multi-speed models [START_REF] Mcnamara | Stabilization of thermal lattice Boltzmann models[END_REF]. Another distribution function is thus introduced to compute this concentration, in addition to the rst distribution function used within the ordinary LBM related to the uid ow only. Compared to the previous work done in LBM topology optimization in thermal uid ow problems [START_REF] Dugast | Topology optimization of thermal uid ows with an adjoint lattice boltzmann method[END_REF], a source term is introduced in this distribution function to represent the consumption of the reactant during the chemical reaction. Reaction laws are very numerous and can be highly nonlinear and multi-parameters dependent (for example Arrhenius law for the temperature dependence). In optimization literature, generally, the reactant equation has a linear dependence with the concentration. At the inlet, the reactant concentration is maximum, before decreasing inside the domain. The chemical reaction is then considered optimized when the reactant is fully consumed at the outlet. Mathematically speaking, this optimization problem consists in minimizing the mean concentration at the outlet. Obviously, one can seek not only for the full use of the reactant but also for the best use of reactant: many reaction problems require an homogeneous reaction rate, in order to avoid either over heating for exothermic reactions or catalyst aging. In this case the optimization is multi-objectives (for example searching for an entire consumption of the reactant and an homogeneous reaction rate at the same time) and requires specic numerical methods; this problematic is beyond the scope of this publication.

In this article, a topology optimization method for an advection-diusion reaction problem using the lattice Boltzmann method and a level-set function is presented. A multiple relaxation time operator is used to compute uid ows with Reynolds numbers up to 1,000 and for dierent Schmidt numbers (ratio between the momentum diusion and the mass diusion) [START_REF] Makhija | Topology optimization of multi-component ows using a multi-relaxation time lattice Boltzmann method[END_REF]. As this study depends on several physical parameters (ow regime, reactant mass diusivity, reaction rate), their inuence on the forward problem and on the optimized geometry is also investigated. A limitation of the pressure drop is also introduced in the optimization problem as a constraint.

The paper is organized as follows. Section 2 introduces the main features of the topology optimization method including the update of the geometry with a level-set function. Section 3 presents the background theory for LBM and its application for an advection-diusion reaction equation. This forward problem is eventually written down concisely in residual form, and the cost function to be minimized is also rewritten in terms of lattice Boltzmann variables. Section 4 presents the derivation of the cost function gradient through the adjoint-state method. The concise residual form is again given. Section 5 then presents the numerical examples related to the topology optimization of an advection-diusion reaction problem. The geometry of the problem is a 2D square cavity with the inlet close to the top-left corner and the outlet close to the bottom-right corner; the reaction occurs in the entire domain.

This academic test case is a variation of the 2D cavity studied by Okkels and Bruus [START_REF] Okkels | Scaling behavior of optimally structured catalytic microuidic reactors[END_REF], but the methodology presented in the current paper can be extended to real-life applications and 3D congurations. First, a study about the inuence of the Reynolds number is presented with Reynolds numbers from 10 to 1,000. This range is chosen to explore solutions at suciently high Reynolds numbers to have a benet in stability with the MRT-LBM model used here, while focusing on the laminar regime to maintain the validity of the model. Then the inuence of the other physical parameters (Peclet or Damkholer numbers) is also discussed. All the obtained results demonstrate that the proposed numerical method is robust and ecient for solving topology optimization problems of reactive uid ows, in various operating conditions. Some discussions and conclusions are nally drawn in section 6.

Denition of the optimization problem

The physics of advection-diusion-reaction

The chemical species are transported by the uid ow and the reaction is taken into account by a source term involved in the equation of the species concentration. This problem is written with the following macroscopic equations ∀x ∈ D, D being an open bounded set of R 2 :

∇ • u = 0, (1)
(u • ∇) u + 1 ρ ∇p -ν∇ 2 u = 0, (2)
u • ∇c A -D∇ 2 c A + s = 0. (3)
The source term s involved in eq. ( 3) is expressed as a function of the concentration of the species of concern,

c A = [A] = n A V .
One possibility is to express a linear dependency between s and c A , as in [START_REF] Kim | Topology optimization of gas ow channel routes in an automotive fuel cell[END_REF]. Elsewhere, in [START_REF] Schäpper | Topology optimized microbioreactors[END_REF] and [START_REF] Yaji | Topology optimization for the design of ow elds in a redox ow battery[END_REF], the reaction rate is limited by the ethanol production, or by the velocity eld, meaning that the reaction rate reaches a maximum value k for a specic concentration. Many other models can be used depending on the chemical reactions, using linear, quadratic, cubic, logarithmic, or exponential relations, to cite but a few [START_REF] Ponce Dawson | Lattice boltzmann computations for reaction-diusion equations[END_REF][START_REF] Weimar | Nonlinear reactions advected by a ow[END_REF]. With a linear source term, s = k c A , when introducing a given characteristic length of the computational domain L, and a characteristic velocity U into the initial set of equations eqs.

∇ • u = 0, (4)
u • ∇ u + ∇ p - 1 Re ∇ 2 u = 0, (5)
u • ∇ c A - 1 Pe ∇ 2 c A + Da c A = 0. (6)
Here tilde values X stand for dimensionless version of X and the three dimensionless numbers Re, Sc and Da are dened as follow:

the Reynolds number, Re = U L ν , represents the ratio between inertial and viscous forces involved in a uid ow. It can help to distinguish ow regimes, from laminar to turbulent ows;

the Schmidt number, Sc = ν D , is the ratio of the kinematic viscosity (momentum diusivity) over the mass diusivity;

the Damkohler number, Da = kL U , represents the ratio of the reaction rate over the transport phenomena rate in the system.

Alternatively, the Peclet number can be used, Pe = Re × Sc. It represents the ratio between the transfer by convection and the transfer by diusion.

In uid ow topology optimization problems, the principal eect of modifying the topology is the modication of the velocity distribution. It means that for the optimization to be eective, the convection process has to play a major role for the concentration distribution when compared to the diusion process. Otherwise, the inuence of the advection process on the eciency of the chemical reaction would be too small. So a high Peclet number is chosen here to respect this condition. However, diusion dominant processes can also be optimized using this method but the problem will be simpler by removing eqs. ( 1) to (2) to the set of equations to be solved. About the reaction rate, a compromise is to be found: if the reaction rate represented by the Damkohler number is too high, the reactant is consumed too fast in the domain, and, there is no need of using any topology optimization algorithms.

Completely on the other hand, if the reactant is consumed too slowly, then the concentration gradients are too small and the impact of the uid/solid geometry will also be negligible. For these two extreme cases, the optimization problem is insensitive because of the bad sizing of the space domain D

: this one has to be decreased (resp increased) in order to improve the sensitivity of the physical problem to the topology of the domain. In all cases in between, for ordinary reaction rates, the use of topology optimization may be very impact-full.

Setting up the optimization problem: reaction maximization

In topology optimization, one searches for the best material distribution or at least a good enough distribution to satisfy a given objective. The material distribution is represented by the design variable, and the objective is dened in terms of a cost function, J . The update of the geometry will be obtained by the evolution of a level-set function, involving the cost function gradient. This gradient is given by the derivative of the cost function J with respect to the design variable. Section 4 is dedicated to such a derivation.

The objective of the optimization problem is the maximization of the reaction, which is equivalent to the minimization of the product of the velocity with the concentration, at the outlet of the domain. Mathematically, this is written as:

J = 1 |∂D out | ∂Dout u • n c A + D∇c A • n dx, ( 7
)
where n is the outward unit normal vector to the boundary, the rst term of the integral is the convective part of the reactant ux while the second term represents the diusive component. As the Pe number is far greater than one in all simulations, the diusive component is neglected in this work.

A constraint is also to be added to the optimization problem in order to allow a maximum amount of pressure drop during the optimization process.

Such a limitation is useful to prevent from unfeasible designs. Rather than adding inequality constraints on the optimization problem, a penalization term is added to the cost function, as in [START_REF] Allaire | Structural optimization with FreeFem++[END_REF], in structural optimization, for example. The augmented cost function, J + , to be minimized, is thus:

J + = J + J 1 with J 1 = ∆p max exp ∆p ∆p max , (8)
and ∈ R + is a user-dened value controlling the relative importance of the dierent contributions J and J 1 .

Geometry update

The design domain is divided into N elements and there is one design parameter per element. The topology is represented by the signed variable α(x) ∈ R N . Then the level-set method is used for the mapping of the design variable. Such a method has been rstly implemented for the representation of interfaces in multiphase ows [START_REF] Osher | Front propagating with courvature-dependent speed: algorithms based on hamilton-Jacoby formulations[END_REF] and for image segmentation [START_REF] Osher | Geometric Level Set Methods in Imaging, Vision, and Graphics[END_REF].

It has then been used for topology optimization, by Sethian et al. [START_REF] Sethian | Structural Boundary Design via Level Set and Immersed Interface Methods[END_REF],

Wang et al. [START_REF] Wang | A level set method for structural topology optimization[END_REF], and Allaire et al. [START_REF] Allaire | Structural optimization using topological and shape sensitivity via a level set method[END_REF].

The continuous level-set function, is such that its zero contour denes the interface between the two media. The update of the solid/uid distribution α is thus performed in two steps. The rst one is via the evolution of the discrete level-set function Ψ(x) ∈ R N :

Ψ (n+1) (x) = Ψ (n) (x) -P (n) ∇ Ψ J + (n) (x). ( 9
)
The superscript is the iteration count, P is the iteration matrix, ideally a good approximation of the inverse of the Hessian matrix, P ≈ H -1 , and ∇ Ψ J + is the augmented cost function gradient (this vector contains the derivatives of the cost function with respect to all components of Ψ). The second step is the topology update which consists in applying the mapping between the level-set function and the topology.

In a previous study, the mapping from the level set function to the topology was dened as [START_REF] Dugast | Topology optimization of thermal uid ows with an adjoint lattice boltzmann method[END_REF]:

α(Ψ(x)) = 1 2 (1 + signΨ(x)) , (10)
such that

α(x) = 0 if Ψ(x) < 0; 1 if Ψ(x) > 0. (11)
A smooth version of the mapping can be considered, for example [START_REF] Aghasi | Parametric level set methods for inverse problems[END_REF]:

α (Ψ(x)) = 1 2 + 1 π arctan Ψ(x) . (12)
The main interest in using ( 12) rather than [START_REF] Challis | Level set topology optimization of uids in Stokes ow[END_REF] is that the former is a continuous and dierentiable function, with α (Ψ) nite everywhere, as soon as = 0. Besides, if → 0, a clear discontinuous uid/solid interface is obtained, which is useful while solving the forward problem. Indeed, the LBM implementation used here for the uid ow and for the concentration eld does not involve a force term to deal with a Brinkmann penalization for a porous media. As a consequence, each node must be either strictly uid or solid, and no intermediate state is allowed. It means that a threshold is applied after the update of the level-set function to obtain a clear uid/solid interface for the forward problem. The choice of this parameter is discussed in the section dedicated to numerical results. In the following, the subscript is avoided for readability considerations.

Note that other strategies such as the Brinkmann penalization for modeling the porous medium, coupled with the SIMP method could have been used [START_REF] Makhija | Topology optimization of multi-component ows using a multi-relaxation time lattice Boltzmann method[END_REF]. However, the level-set function allows a clearer denition of the topology, and thus minimizes the possible bias involved within the forward model.

Lattice Boltzmann method

Multi-relaxation time model for the uid ow

The lattice Boltzmann method (LBM) is an alternative way of solving uid ow problems compared to classical methods based on the discretization of the NavierStokes equations (nite volume or nite element methods).

This method is based on the calculation of the distribution function f involved in the Boltzmann equation and which gives the probability of nding a particle at the position x, at the time t and at the speed c. Compared to the Boltzmann equation, the velocity vector is discretized along directions in the lattice Boltzmann equation, such as [START_REF] Wolf-Gladrow | Lattice-Gas Cellular Automata and Lattice Boltzmann Models -An Introduction[END_REF][START_REF] Succi | The Lattice Boltzmann Equation for Fluid Dynamics and Beyond[END_REF]:

∂f i ∂t + c i • ∇f i = Ω i (f ) ∀i ∈ J0; I -1K. (13)
The notation f = {f i } I-1 i=0 will be used for the distribution function.

Such a notation will be used for other quantities, as soon as there is no possible ambiguity. In two dimensions, the reference discretization scheme involves I = 9 discrete directions yielding the so-called D2Q9 scheme. The directions and associated weights are given by [START_REF] Qian | Lattice BGK Models for Navier-Stokes Equation[END_REF]:

c = 0 1 0 -1 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 , (14)
w = 4
Also, from the denition of the velocity directions, it is convenient to introduce the opposed velocity directions ī with respect to c 0 , to write the bounce-back boundary conditions:

ī = 0 1 2 3 4 5 6 7 8 = 0 3 4 1 2 7 8 5 6 . ( 16
)
The lattice Boltzmann equation ( 13) is generally solved in two steps: the collision step, and the streaming step. For the collision operator, Ω(f ), the simple model introduced by Bhatnagar-Gross-Krook (BGK) [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF][START_REF] Marié | Etude de la méthode Boltzmann sur Réseau pour les simulations en aéroacoustique[END_REF] which relies on a single-relaxation time (SRT) is often used. However, this is often a source of spurious oscillations and instabilities as the Reynolds number increases. As an example, d'Humières et al. [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF] reported severe oscillations for the pressure eld in the cavity ow for the SRT model at Re=500. So, in order to improve the numerical stability of the LBM, a multi-relaxation time (MRT) operator has further been introduced [START_REF] Humières | Multiple-relaxation-time lattice Boltzmann models in three dimensions[END_REF]. It is written as:

Ω(f ) = -M -1 SM (f -f eq ) , (17)
in which M is the transformation matrix between the moments and the distribution functions:

M =               1 1 1 1 1 1 1 1 1 -4 -1 -1 -1 -1 2 2 2 2 4 -2 -2 -2 -2 1 1 1 1 0 1 0 -1 0 1 -1 -1 1 0 -2 0 2 0 1 -1 -1 1 0 0 1 0 -1 1 1 -1 -1 0 0 -2 0 2 1 1 -1 -1 0 1 -1 1 -1 0 0 0 0 0 0 0 0 0 1 -1 1 -1               , ( 18
)
and S is the diagonal relaxation rate matrix dened by [START_REF] Liu | A High Scalable Hybrid MPI / OpenMP Parallel Model of Multiple-relaxation-time Lattice Boltzmann Method Multiple-relaxation-time Lattice Boltzmann Method[END_REF]:

S = diag(1, s e , s , 1, s q , 1, s q , s ν , s ν ). (19)
Note that s i = 1 corresponds to conserved moments (conservation of mass and momentum), and s ν is related to the uid viscosity:

1 s ν = 3ν + 1 2 . ( 20
)
The other components are to be adjusted in order to improve the numerical stability and are often determined to be close to one through an optimization process [START_REF] Tekitek | Adjoint lattice Boltzmann equation for parameter identication[END_REF][START_REF] Xu | Sensitivity analysis and determination of free relaxation parameters for the weakly-compressible mrt-lbm schemes[END_REF]. In this article the uid used in the cavity is water, so an incompressible model is used for the LBM. The equilibrium distribution f eq i is then given by [START_REF] Zou | A improved incompressible lattice Boltzmann model for time-independent ows[END_REF][START_REF] He | Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation[END_REF]:

f eq i = ω i ρ + ρ 0 3α(c i • u) + 9 2 α(c i • u) 2 - 3 2 αu 2 . ( 21
)
The design-related variable α(ψ(x)) involved in eq. ( 21) has been dened in section 2. The macroscopic quantities, the uid density ρ, the velocity u, and the pressure can further be computed with the moments of f :

ρ(x, t) = 8 i=0 f i (x, t), ρ 0 u(x, t) = 8 i=0 c i f i (x, t), p(x, t) = c 2 s ρ(x, t), (22)
where ρ = ρ 0 + δρ, δρ being the density uctuation and the mean density ρ 0 is chosen to be equal to 1 in our simulations [START_REF] He | Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation[END_REF], and c s = 1/ √ 3 the lattice speed of sound. The uid velocity needs to be very low compared to c s in order to satisfy the low mach number assumption required for the equilibrium distribution function expansion [START_REF] Chen | Lattice Boltzmann method for uid ows[END_REF]. Compared to the classical model, neglecting the density uctuations in terms involving velocity induces a decoupling between density and velocity [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF]. While reducing the errors due to compressibility [START_REF] Mei | Consistent initial conditions for lattice Boltzmann simulations[END_REF], it should be noted that this model is only valid for steady-state ows [START_REF] Dellar | Lattice Boltzmann algorithms without cubic defects in Galilean invariance on standard lattices[END_REF]. Next, the bounce-back boundary condition is applied on all interfaces but on the inlet and outlet, where the boundary conditions proposed by Zou and He [START_REF] Zou | On pressure and velocity ow boundary conditions and bounceback for the lattice Boltzmann BGK model[END_REF] are applied. The user-parameter involved in ( 12) is chosen to tend to zero within this multi-relaxation time lattice Boltzmann model, so that α equals either 0 or 1, and a clear denition of the topology is considered to avoid any bias in the uid ow forward modeling.

Single-relaxation time model for the concentration distribution

Based on a double distribution function approach, a second distribution function, g, similar to f , is introduced to solve the concentration eld [START_REF] Yan | Numerical simulation of heat transfer and uid ow past a rotating isothermal cylinder A LBM approach[END_REF][START_REF] Li | Lattice Boltzmann models for the convection-diusion equation: D2Q5 vs D2Q9[END_REF]:

∂g i ∂t + c i • ∇g i + s i = Ω i ∀i ∈ J0; I -1K, (23)
with the collision operator dened as:

Ω i = - 1 τ g (g i -g eq i ) . ( 24
)
Although the multi-relaxation time collision operator is used for the distribution function f , the simpler single-relaxation time model is used here, with only one relaxation time τ g for the distribution function g. The lattice Boltzmann equation for the concentration is more stable than the one for the uid ow, and the multi-relaxation time model is not mandatory for this distribution function.Indeed, the SRT-BGK model is suitable to solve this problem as long as the diusion is isotropic [START_REF] Chai | A Multiple-Relaxation-Time Lattice Boltzmann Model for General Nonlinear Anisotropic Convection-Diusion Equations[END_REF][START_REF] Rasin | A multi-relaxation lattice kinetic method for passive scalar diusion[END_REF], which is the case here.

The relaxation time τ g is related to the mass diusivity D of the chemical species:

τ g = 3D + 1 2 . ( 25
)
In [START_REF] Succi | The lattice Boltzmann equation: A new tool for computational uid-dynamics[END_REF], the source term involved in the reaction process equation is distributed on all the discrete velocity directions, according to the D2Q9 scheme.

For a linear source term, s = k c A , one can write:

s i = ω i k c A ∀i ∈ J0; I -1K. ( 26
)
More generally, the LBM formulation for a given source model s is relaxed towards the dierent velocities directions as follow:

s i = ω i s. ( 27
)
The equilibrium distribution function g eq i is given by [START_REF] Pingen | Topology optimization for thermal transport[END_REF]:

g eq i = ω i c A (1 + 3α c i • u) . ( 28
)
It can be noticed that the velocity eld is used to compute this equilibrium distribution function, illustrating the coupling between the concentration and the velocity.The concentration c A is computed with the zeroth order of the distribution function g such as [START_REF] Li | Lattice Boltzmann models for the convection-diusion equation: D2Q5 vs D2Q9[END_REF] :

c A (x, t) = 8 i=0 g i (x, t) (29)
The mass ux is written as [START_REF] Yaji | Topology optimization in thermal-uid ow using the lattice Boltzmann 42 method[END_REF]:

q(x, t) = 8 i=0 c i g i (x, t) -c A u (30)
In order to implement the boundary conditions, the unknown distribution functions are assumed to be the equilibrium distribution functions with the concentration c A to be determined [START_REF] Dugast | Topology optimization of thermal uid ows with an adjoint lattice boltzmann method[END_REF][START_REF] Inamuro | A Lattice Boltzmann Method for a Binary Miscible Fluid Mixture and Its Application to a Heat-Transfer Problem[END_REF]. For the prescribed inlet concentration, c A is calculated by [START_REF] Liu | Discrete adjoint sensitivity analysis for uid ow topology optimization based on the generalized lattice Boltzmann method[END_REF], with c A = c Ain . At the outlet, c A is calculated by [START_REF] Obrecht | Multi-GPU implementation of the lattice Boltzmann method[END_REF], with q = 0.

A similar solver has been validated for a uid ow and heat transfer problem by the authors [START_REF] Dugast | Optimisation topologique en convection thermique avec la méthode Lattice Boltzmann[END_REF]. In that case, the passive scalar approach is also used, the only dierence is that the concentration is replaced by the temperature.

Residual form

In order to calculate the cost function gradient, the adjoint-state method will be used, see section 4. In order to set it up, the forward problem, which consists of both the multi-relaxation time model for the uid-ow, and the single-relaxation time model for the concentration distribution, is written down in terms of residuals. We denote R and P the residuals for the Boltzmann state equations, and the boundary conditions, respectively.

Further, we use the superscripts f and g for the uid ow, and for the concentration, respectively. Such residuals have to be written down for each velocity direction. Added to that, the boundary ∂D of the medium D is partitioned with the inlet D in , the outlet D out , and the no-slip boundaries D ns . With such notations, the forward problem is concisely written as:

search (f i , g i )(x, t), i ∈ J0, ..., 8K such that: • ∀x ∈ D : R f i (f , α) = ∂f i ∂t + c i • ∇f i + (M -1 SM (f -f eq )) i = 0, R g i (f , g, α) = ∂g i ∂t + c i • ∇g i + s i + 1 τg (g i -g eq i ) = 0,
• ∀x ∈ ∂D in :

P f 1 (f ) = f 1 -f 3 -2 3 u in = 0, P f 5 (f ) = f 5 -f 7 -1 6 u in -1 2 (f 4 -f 2 ) = 0, P f 8 (f ) = f 8 -f 6 -1 6 u in + 1 2 (f 4 -f 2 ) = 0, P g 1 (g) = g 1 -1 9 χ (1 + 3u in ) = 0, P g 5 (g) = g 5 -1 36 χ (1 + 3u in ) = 0, P g 8 (g) = g 8 -1 36 χ (1 + 3u in ) = 0,
• ∀x ∈ ∂D out :

P f 3 (f ) = f 3 -f 1 + 2 3 η = 0, P f 6 (f ) = f 6 -f 8 + 1 6 η -1 2 (f 4 -f 2 ) = 0, P f 7 (f ) = f 7 -f 5 + 1 6 η + 1 2 (f 4 -f 2 ) = 0, P g 3 (f , g) = g 3 -1 9 ζ (1 -3η) = 0, P g 6 (f , g) = g 6 -1 36 ζ (1 -3η) = 0, P g 7 (f , g) = g 7 -1 36 ζ (1 -3η) = 0,
• ∀x ∈ ∂D ns :

P f (c i •n<0) (f ) = f (c i •n<0) -f (cī•n>0) = 0, P g (c i •n<0) (g) = g (c i •n<0) -g (cī•n>0) = 0, (31)
with:

η = -ρ out + f 0 + f 2 + f 4 + 2(f 1 + f 5 + f 8 ), χ = 6(c Ain -g 0 -g 2 -g 3 -g 4 -g 6 -g 7 ) (1+3u in ) , ζ = 6(g 1 +g 5 +g 8 ) (1+3η) . (32)

LBM-based cost function

The cost function written in terms of the primal variables is given combining eqs. ( 7) and ( 8), which gives:

J + (u, c A , ∆p) = 1 |∂D out | ∂Dout u • n c A dx + ∆p max exp ∆p ∆p max . (33)
This cost function is to be re-written in terms of the lattice Boltzmann variables, f and g. Taking into account of the following relationships between the primal variables and the lattice Boltzmann variables:

u = 8 i=0 c i f i (x, t) c A = 8 i=0 g i (x, t) ∆p = 1 |∂D in | 1 3 ∂D in 8 i=0 f i dx -1 3 ∂Dout 8 i=0 f i dx , (34)
the cost function that is, nally, to be minimized is:

J + (f , g) = 1 |∂D out | ∂Dout 8 i=0 c i f i • n 8 i=0 g i dx + ∆p max exp 1 3 ∂D in 8 i=0 f i dx |∂D in | ∆p max . (35)
Note that the pressure term at the outlet is not involved here. This one being considered as prescribed, the pressure dierence is only driven by the modication of the inlet pressure during the optimization.

The cost function gradient

In order to derive the full adjoint-state problem needed for the computation of the cost function gradient, the methodology and notations of Gunzburger [START_REF] Gunzburger | Adjoint equation-based methods for control problems in incompressible, viscous ows[END_REF] are introduced:

F (φ, Ψ) = 0 is the forward model gathering all the equations of eq. ( 31); φ = (f , g) is the global forward state;

φ * = (f * , g * ) is the global adjoint-state.
The Lagrange function associated to the optimization problem is written as 1 :

L (φ, Ψ, φ * ) = J + (φ) + F (φ, Ψ), φ * . ( 37
)
The optimization problem is solved searching for the stationary point of the Lagrange function. This point satises:

δL = ∂L ∂φ δφ + ∂L ∂Ψ δΨ + ∂L ∂φ * δφ * = 0, (38)
where δL, δφ, δΨ and δφ * stand for arbitrary variations [START_REF] Gunzburger | Adjoint equation-based methods for control problems in incompressible, viscous ows[END_REF]. Three terms appear in the right hand side of (38):

the rst term, the dierential with respect to the state, ∂L/∂φ, yields the adjoint-state equation;

the second term, the dierential with respect to the design variables, ∂L/∂Ψ, yields the optimality conditions;

1 With the notations given above, this eq. ( 37) can be straitforwardly expanded to:

L(f , g, Ψ, f * , g * ) = J + (f , g) + t f 0 D q={f,g} i=0,...,8 R q i q * i dx + ∂Din q={f,g} i={1,5,8} + ∂Dout q={f,g} i={3,6,7} + ∂Dns q={f,g} ci•n<0 (P q i q * i dx) dt. ( 36
)
the third term, the dierential with respect to the adjoint-state variable, ∂L/∂φ * , yields back the equations of the forward model that are to be satised.

A similar derivation of the adjoint-state problem for a convective uid ow problem has been detailed by the authors [START_REF] Dugast | Topology optimization of thermal uid ows with an adjoint lattice boltzmann method[END_REF]. Based on this method, the adjoint-state is nally written as :

search (f * i , g * i )(x, t), i ∈ J0, ..., 8K, with (f i , g i )(x, t)
given by [START_REF] Delbosc | Real-Time Simulation of Indoor Air Flow Using the Lattice Boltzmann Method on Graphics Processing Unit[END_REF], such that:

• ∀x ∈ D : R * ,f i (f , g, f * , g * , α) = - ∂f * i ∂t -c i • ∇f * i + (M -1 SM (f * -f * ,eq )) i + Q * ,f i = 0, R * ,g i (f , g, g * , α) = - ∂g * i ∂t -c i • ∇g * i + 1 τg (g * i -g * ,eq i ) + Q * ,g i = 0,
• ∀x ∈ ∂D in :

P * ,f 3 (f * , f ) = f * 3 -f * 1 + 1 3|∂D in | exp 1 3 ∂D in 8 i=0 f i dx |∂D in |∆pmax P * ,f 7 (f * , f ) = f * 7 -f * 5 + 1 3|∂D in | exp 1 3 ∂D in 8 i=0 f i dx |∂D in |∆pmax P * ,f 6 (f * , f ) = f * 6 -f * 8 + 1 3|∂D in | exp 1 3 ∂D in 8 i=0 f i dx |∂D in |∆pmax P * ,g 3 (g * ) = g * 3 + 1 6 (4g * 1 + g * 5 + g * 8 ) P * ,g 6 (g * ) = g * 6 + 1 6 (4g * 1 + g * 5 + g * 8 ) P * ,g 7 (g * ) = g * 7 + 1 6 (4g * 1 + g * 5 + g * 8 )
• ∀x ∈ ∂D out :

P * ,f 1 (f * , g * , f , g) = f * 1 -f * 3 + 1 3 (4f * 3 + f * 6 + f * 7 ) + 1 6 ζ 1+3η (4g * 3 + g * 6 + g * 7 ) + c A (g) |∂Dout| = P * ,f 5 (f * , g * , f , g) = f * 5 -f * 7 + 1 3 (4f * 3 + f * 6 + f * 7 ) + 1 6 ζ 1+3η (4g * 3 + g * 6 + g * 7 ) + c A (g) |∂Dout| = P * ,f 8 (f * , g * , f , g) = f * 8 -f * 6 + 1 3 (4f * 3 + f * 6 + f * 7 ) + 1 6 ζ 1+3η (4g * 3 + g * 6 + g * 7 ) + c A (g) |∂Dout| = P * ,g 1 (g * , f ) = g * 1 -1 6 1-3η 1+3η (4g * 3 + g * 6 + g * 7 ) + η |∂Dout| P * ,g 5 (g * , f ) = g * 5 -1 6 1-3η 1+3η (4g * 3 + g * 6 + g * 7 ) + η |∂Dout| P * ,g 8 (g * , f ) = g * 8 -1 6 1-3η 1+3η (4g * 3 + g * 6 + g * 7 ) + η |∂Dout| • ∀x ∈ ∂D ns : P * ,f (c i •n>0) (f * ) = f * (c i •n>0) -f * (cī•n<0) = 0, P * ,g (c i •n>0) (g * ) = g * (c i •n>0) -g * (cī•n<0) = 0. (39)
with:

f * ,eq i = 8 j=0 ∂f eq j ∂f i f * j , Q * ,f i = 8 j=0 -g * j 3 ω j c A τg c j • c i , g eq, * i = 8 j=0 ∂g eq j ∂g i g * j , Q * ,g = 8 j=0 ∂s j ∂g i g * j , η = -ρ out + f 0 + f 2 + f 4 + 2(f 1 + f 5 + f 8 ), ζ = 6(g 1 +g 5 +g 8 ) (1+3η) . ( 40
)
To compute the dierent adjoint-state variables, the steady-state solution of the forward problem is used. The adjoint-states being computed, the cost function gradient is nally computed through

∇ Ψ J + (f , g, f * , g * ) = -α (Ψ) t f 0 8 i=0 ω i f * i M -1 SM 3 c j • u + 9 2 (c j • u) 2 - 3 2 u 2 i dt -α (Ψ) t f 0 8 i=0 3 ω i g * i c A c i • u τ g dt. (41)
The general topology optimization algorithm is detailed in Algorithm 1.

The iterative algorithm involves a criterion based on the stabilization of the cost function value which is detailed in the next section.

Test cases and results

For the validation of the proposed method, the reader is referred to a previous work of the authors [START_REF] Dugast | Topology optimization of thermal uid ows with an adjoint lattice boltzmann method[END_REF] where a thermal uid ow topology optimization case has been computed by following a similar procedure. The obtained results have been then compared with a benchmark test in literature and a study about the grid independency and the sensitivity of the algorithm to dierent initializations have also been conducted.

The geometry and the conguration of the test cases are shown on g. 1.

The 2D square domain is enclosed by no-slip walls. The gray layers, shown The convergence criterion for both the forward LBM problem and the adjoint-state LBM problem has been chosen to be:

|a (n) -a (n-10,000) ∞ a (n-10,000) ∞ < 10 -4 with a (n) ∞ = max q={f,g} j=0,...,8 i=1,...,N q (n) j (x i ) , (42)
in which the innity norm is applied on the vector involving all the nine directions, and all spatial nodes, for both quantities, and the superscript (n) stands for the LBM iteration count. Besides, the convergence criterion of the optimization problem is based on the stabilization of the cost function value along with ten successive iterations:

J (k) -J (k-10) J (k-10) < 10 -4 , (43)
where J (k) is the cost function value at the iteration count k. Calculations were performed on a NVIDIA Quadro K6000 GPU card for taking advantage of the LBM algorithm parallelism. Normally, for the adjoint-states calculation, one needs to store the macroscopic values of the forward problem, at each LBM iteration, which is prohibitive in terms of memory requirement.

Only steady-state problems are considered here, so, as stated by [START_REF] Liu | Discrete adjoint sensitivity analysis for uid ow topology optimization based on the generalized lattice Boltzmann method[END_REF], the solution of the forward problem at nal time can be used as the stationary source term for the whole adjoint-states calculation.

Regarding the reaction source term, instead of the simple linear model introduced previously in the adimensionnalised Damkohler number, a twoparameters model is chosen for a more tunable reaction:

s = k (1 -exp(-rc A )) , (44)
The coecient k controls the maximum reaction rate, and r is a positive user-dened parameter that controls the dependency to the concentration: a high value prescribes a quasi-independent concentration rate whereas a small value tends to a linear model.

The numerical examples presented in this section use the following parameters:

the spatial discretization for the domain is 200 × 200 elements, the reaction rate coecient k is equal to 2 × 10 -5 , the Peclet number Pe is 2000, the initialization starts with the full uid topology, i.e. with Ψ(x) > 0, the iteration matrix P used for the level-set function evolution in eq. ( 9)

is chosen to be a diagonal matrix satisfying P i = 0.01 /α (Ψ i ). The benet of this choice is threefold: i) the derivative α (Ψ) ≥ 0 is no longer involved in the level-set function evolution, so even a discontinuous function α(Ψ) can be used; ii) the user-parameter is not to be pre-assigned to a given value, this parameter being no longer involved anywhere in the optimization process; iii) there is no more bias between the response of the forward model, i.e. the value of the cost function J which would need = 0, and its gradient, ∇ Ψ J , which would need = 0 for its existence condition. Next, the descent parameter 0.01 involved at the numerator of the iteration matrix, has been chosen following Dugast et al. [START_REF] Dugast | Topology optimization of thermal uid ows with an adjoint lattice boltzmann method[END_REF].

the maximal allowed pressure drop is twice the initial pressure drop (calculated at the rst iteration with the full uid geometry), the relaxation times for the MRT are set to s e = s = s q = 0.6 and

s ν = 1/(3ν + 0.5).
the reaction rate coecient r is equal to 10, which corresponds to a Damkholer number Da equal to 0.1.

As the MRT model allows the calculation of higher Reynolds number ows compared to the SRT model, the rst study described hereafter is based on the inuence of the Reynolds number on the optimized topology. Three cases are presented : Re =10, Re =100 and Re =1000. It is to be noticed that although the forward and adjoint-state problems could have been solved with the single-relaxation time model, for both Re=10 and Re=100, the convergence of the LBM solver is not possible without a modication of the relaxation times s 0-6 for Re =1000. This study shows that the introduction of the solid parts breaks the main uid ow path going from the inlet to the outlet via the center of the domain into several uid paths. This is true for the three tested Reynolds numbers. This path division participates to a better distribution of the uid ow on the entire domain, especially on the top-right and bottom-left corners compared to the initial conguration. One can notice that the complexity of the ow paths increases with the Reynolds number. Apart from advection, the forward problem involves also diusion and reaction phenomena. Consequently, the second study deals with the inuence of the mass diusivity and of the reaction rate on the obtained topology. The modication of the mass diusivity aects the balance between advection and diusion and therefore the impact of the uid ow velocity on the concentration. The reaction rate is used to control the speed on the reaction. The inuence of a faster reaction will be studied. With these two studies, the topology optimization method has been tested for dierent ow regimes and reactants. It shows its ability to produce specic optimized topologies for a variety of congurations.

Inuence of the Reynolds number

The rst study deals with the inuence of the Reynolds number on the obtained optimized topology. The objective is the maximization of the reaction within the domain, this measurement being calculated at the outlet, see eq. [START_REF] Li | Lattice Boltzmann models for the convection-diusion equation: D2Q5 vs D2Q9[END_REF]. For this study, the Peclet number is kept constant. As Pe = Re × Sc, it means that the Schmidt number is also modied along with the Reynolds number.

[Re=10] at the iteration 0 (it is the initial guess), T

[Re=10] at the iteration 30, T

[Re=10] at the iteration 50, and after reaching the stabilization of the cost function value,

T * [Re=10]
, at the iteration 900.

From this gure it is seen that the solid parts are introduced only in the center of the domain, but not close to both the inlet and the outlet. That is the reason why the increase in the pressure drop ratio is not much, and the constraint on the pressure drop is by far fullled. In fact, the introduction of the solid medium close to the center of the medium divides the ow in two distinct paths, as can be seen from g. 2, in the middle. This makes such that the uid velocity is greater elsewhere, especially close to the bottom right and top left corners. The eect is that more reactant is consumed in these areas; this participates to the high decrease of the outlet concentration (see the bottom of g. 2). Note that, from the initial guess to the optimized topology, the cost function has decreased from 5.10 -3 (LB unit) down to 1.96.10 -3 (LB unit), i.e. by a factor more than 2.5.

[Re=100] at the iteration 0 (it is the initial guess), T

[Re=100] at the iteration 20, T

[Re=100] at the iteration 50, and after reaching the stabilization of the cost function value,

T * [Re=100]
, at the iteration 390. Comparing g. 3 with g. 2, it can be seen that the decrease of the cost function is more important for Re=100 than for Re=10. Here, the cost function has decreased from 6.1.10 -3 (LB unit) down to 5.96.10 -4 (LB unit), i.e. by a factor greater than 10 .

Conversely to the previous case with Re=10, the optimal uid/solid conguration is here composed of several parts, located near the diagonal of the domain, facing the main uid ow stream. The uid ow is thus mainly directed towards the corners, but also a part of it goes through the center of the domain, as the solid line is not continuous, containing some holes, see the plot in the middle of g. 3.

The outlet concentration for the optimized geometry is very close to zero, as can be seen from the plot in the bottom of g. 3. This means that almost all the reactant has been consumed inside the domain. In that sense, the eciency of the chemical reaction has been highly increased.

Case 3 : Re =1000

In this last test case about the Reynolds number study, the velocity and concentration elds for Re = 1000 are radically dierent to previous cases, with Re = 10 and Re = 100. As a matter of fact, the g. 4 is to be compared to g. 2 and g. 3, middle and left-hand side for the velocity, and bottom and left-hand side for the concentration, at the initial guess, i.e. for the full uid ow case. In fact, when Re = 1000, a re-circulation zone appears around the center of the domain. As a consequence, the concentration is very small there at the center of the vortex, but the concentration remains very high at the outlet.

The uid/solid conguration found by the implemented topology optimization algorithm is, indeed, very far away from an optimized one, as the optimization process is struggling to nd any proper topology. In order to facilitate the convergence of the optimization problem, the initial congura- Right: optimized topology.

tion was set to a uniform distribution of small solid circles inside the domain.

With such an initialization, the re-circulation disappears, and the research of optimized geometries is easier. Note that the cost function value for the full uid geometry was 12.6.10 -3 (LB unit) while it is now 4.5.10 -3 (LB unit) with this initialization, i.e. with circles.

The evolution along with the optimization iterations of the cost function and of the pressure drop limitation are shown on g. 5, for this case. At the beginning of the optimization process, the pressure drop increases. This is due to the fact that some solid parts are introduced at the outlet of the domain. This may be seen as a drawback, but in fact, this also participates to the decrease of the concentration eld at the outlet, which is the main objective of the optimization problem. Then, after some iterations, the limitation of the pressure drop included in the augmented cost function helps to prevent from the introduction of more solid parts near this area. The dierent solid parts introduced inside the domain are found to be very small: it is even dicult to see the dierence between the dierent optimization steps, without a close look. But these small modications have a huge impact on the cost function decrease.

As can be seen on g. 5, a uid path has been drawn on the left-hand side of the domain, by removing some solid elements from the initial topology.

The distribution of the reactant is then more uniform, and the reaction rate is more important in this area, which participates to a diminution of the outlet concentration, as can be seen on the plot on the bottom of g. 5. In that case again, the eciency of the chemical reaction has been highly increased.

Crosscheck of the dierent optimized congurations

Based on the Reynolds number test cases, a crosscheck comparison of the dierent optimal topologies has been performed in order to discuss on the results. In this additional study, the uid ow with Reynolds numbers equal to 10, 100, and 1000 is applied on the three previously obtained optimal topologies. The cost function measuring the reaction eciency is then calculated for each of these nine cases. The results are reported in table 1.

The rst row is related to the three dierent found optimal topologies, and the rst column is related to the Reynolds number which is applied. The other numbers are the values of the calculated cost function.

The expected result is that the minimum value of the cost function is to be found for the applied Reynolds number equal to the one related to its optimal geometry. For readability considerations, the minimum value of T *

[Re=10]

T *

[Re=100]

T *

[Re=1000]

Re =10 1.9 0.8 has been actually reached. In fact, it is well known that gradient-type methods, because they are based on a local study of the function to be minimized, usually reach a local minimum next to the initial guess [START_REF] Bruns | Topology optimization of convection-dominated, steadystate heat transfer problems[END_REF]. This is usually considered as the main drawback of these gradient-type algorithms. But, at the same time, compared to gradient-free optimization methods which, indeed, are more likely to nd the global minimum, gradient-type algorithms reach the minima much more eciently. This point makes the dierence.

In fact, the optimization problem treated here involves 200×200 unknowns; such a discretization would be impracticable with any gradient-free optimizer.

This diculty could of course be overcome by modifying the initial topology, many times, but, as the physical problem is complex, it will remain dicult (and likely impossible) to nd the unique global minimum.

In further perspectives, the local minimum treatment is an interesting issue to solve, despite the fact that it can be dicult to avoid them while using a gradient-based optimization method. Nevertheless, one possibility would be to use a multi-scale resolution approach. Such a multi-resolution approach, formalized by Liu [START_REF] Liu | A multiresolution method for distributed parameter estimation[END_REF] for solving some ill-posed inverse problems, has been later on successfully used by Dubot et al. [START_REF] Dubot | A wavelet multiscale method for the inverse problem of diuse optical tomography[END_REF] and Liu [START_REF] Liu | A wavelet multiscale method for the inverse problem of a nonlinear convectiondiusion equation[END_REF], for example. In the framework of reactive uid ow topology optimization, the use of a multi-scale resolution approach would enable to perform a scale-byscale optimization on successive convex cost functions, following the work of Chavent [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-by-step guide for applications[END_REF]. Moreover, added to this interesting property, this would also save computational time due to faster convergence. This point will be addressed in a future research.

Following the study of the Reynolds number, the inuence of both the reaction rate and diusion is investigated. Additionally, the Damkholer and Peclet number are modied to determine their eect on the resultant system geometry.

Case 4 : Inuence of the Damkohler number

The reaction rate k is now modied in order to evaluate the inuence of the Damkohler number on the optimal topology. It has been increased from 2 × 10 -5 to 4 × 10 -5 . Though this modication does not have any impact on the velocity, the concentration eld is altered, as the reactant is consumed faster. This eect is clearly visible by comparing gs. 5 and 6.

Concerning the optimization results, one can observe a similar trend between this case and the previous case 3. First, the evolution of the cost function and the pressure drop ratio, represented on the top of g. 6, are similar: one observes an increase and then a decrease of the pressure drop ratio, while the outlet concentration is decreasing. The optimal topology is also comparable to the one of case 3: it implies the same modication on the velocity in which the ow is divided into two main streams just after the inlet. The main dierence between both the previous case 3 and this case 4 is in the magnitude of the outlet concentration. This one is much smaller for the higher Damkholer number, as it was expected.

One can also see from g. 7 (left) that the solid part at the core of the domain is denser for Da=0.1 (in red) than for Da=0.2 (in blue). Please note that, to enhance readability, the solid parts have been represented larger (in this gure only). For Da=0.2, because the reactant is consumed faster (with respect to Da=0.1), it is consequently easier to obtain a small concentration at the outlet, even without the need to alter the velocity eld through the addition of signicant solid material. The comparison of the concentration eld for the two Damkholer numbers is shown on g. 7 (right). The threshold c A ≷ 0.05 is applied to the concentration eld to get a better contrast between the two cases. For the higher Damkholer number, one can see that the concentration becomes very low shortly after the middle of the domain. As the reactant has been almost entirely consumed, the eciency has been highly increased regarding to the cost function. However, the other conclusion that can be made based on this result is that the reactor has been over-sized for this reaction rate. In fact, based on the concentration eld obtained for the optimized geometry, one can conclude that a good eciency could have also been achieved with a smaller reactor. The real cost of the chemical reactor and the dimensions associated with it were not included in the optimization problem; it is nevertheless a useful conclusion from an engineering point of view.

Case 5 : Inuence of the Peclet number

In the previous numerical examples, the Peclet number was equal to 2000.

The advection process played a major role compared to the diusion process.

In this test case 5, the Peclet number is modied to be equal to 100. For the rest, the conguration is kept the same as for the previous case 3, and the Reynolds number is 1000. Note that in this case, the diusion process is no longer negligible in the chemical reaction. Concerning the forward problem, the velocity eld is not altered by this modication. However, important changes are observed in the concentration eld. The pressure drop ratio is important for this case, as can be seen from g. 8, on the top left. This a symmetry of the concentration can be observed, from the inlet down to the outlet. One can notice that, to obtain this behavior, the uid ow is progressively divided into smaller channels from the inlet to the center of the domain. Indeed, the density of solid parts is the most important on the diagonal connecting the bottom-left corner to the top-right corner. With such a pattern, the minimization of the outlet concentration is well achieved.

Overview and analysis of the dierent results

The characteristics of the optimized geometries obtained for all cases are summarized in table 2.

Among all cases, the decrease of the cost function is less important for Re=10 (case 1) and for Pe=100 (case 5). In these two cases, either the advection is the smallest (Re=10, case 1), or the diusion is the greatest (Pe=100, case 5), and the consequences are the same: the impact of the velocity eld on the concentration is reduced, and so is the impact of the topology on the chemical reaction eciency.

For the study on the Damkholer number, the reaction rate is multiplied by 2 (from the test case 3 to the test case 4), and the impact of this modication is consequent on the nal value of the cost function. The initial value was also lower than the other cases, about one order of magnitude compared to the case 3. As it has been stated in section 5.2, the major conclusion regarding the optimized geometry obtained for this case is not only the eciency of the reaction, but rather the over-sizing of the reactor. In g. 7, one can see that there are less solid parts in the core of the domain for Da =0.2 than for Da =0.1. The porosity value in table 2 does not reect this, but this can be explained by the fact that for Da =0.2, the biggest solid part (which is curved and close to the inlet) represents alone 0.37 % of the total porosity.

From table 2, it is seen that the pressure drop ratio slightly increases for Re=10 (case 1) and for Re=100 (case 2), but this pressure drop decreases for Re=1000 (case 3) and Da=0.2 (case 4). In all cases, as the maximal pressure drop is prescribed to twice the pressure drop given for the initial geometry, the initial pressure drop ratio is equal to 0.5. So, for the two rst cases, any added solid within the void medium would increase this pressure drop ratio. However, as the initialization for Re=1000 is realized with a uniform distribution of small circles, the pressure drop can be smaller to the initial one if some circles are removed, especially near the inlet and near the outlet of the domain. This phenomenon happens for Re=1000 (case 3) and for Da=0.2 (case 4). For the test case 5 (study on the Peclet number), the initialization also involves circles, but, as some additional solid parts are introduced, especially near the inlet and near the outlet, this participates to an increase of the pressure drop ratio. As it can be seen on the dierent gures concerning this study, the porosity for a lower Peclet number is higher than for the other cases. The optimization domain is composed of more solid parts to have a stronger control on the velocity eld, useful to counter-balance the higher diusivity, but this is costly in term of pressure drop. One can notice that in all the cases, the optimization process is achieved by the addition of a small amount of materials (98.26 % of minimal nal porosity, except for the low Pe number test). This last conguration is characterized by a lower porosity (97.1 %, which corresponds to three times more added material).

The term porosity not only explains the amount of material but also its structure. Indeed, except at the lowest Re number, where optimal form acts as a prole by splitting the initial major ow, the optimal solid distributions look like porous media. In both the Re=100 case and low Pe number case, a plane structure appears in the diagonal of the cavity, perpendicular to the line joining the entrance and the exit. This structure acts as a perforated plane with adjusted holes in order to control and to balance the uid ow.

In all these congurations, the optimization algorithm has been able to design dierent topologies able to increase the eciency of the chemical reactors, showing its adaptability and usefulness working with several physical parameters involved in a multiphysics problem.

Case # Da Re Pe J * J (0) /J * ∆p ratio Porosity 1 0.

Conclusion

A topology optimization algorithm has been developed for advectiondiusion-reaction equations. The physical goal was the maximization of the reaction within the porous domain. This was measured through the minimization of the concentration ow out of the domain. The evolution of the topology was performed thanks to the evolution of a level-set function and a gradient-type method. Both the forward and the adjoint-state problems have been computed with the lattice Boltzmann method, with a passive scalar approach, and a reaction term for the modeling of the concentration eld.

Though the single-relaxation time collision operator was used for the concentration eld, the ow eld has been conputed with the multi-relaxation time collision operator. This enabled computations with regimes up to Re=1,000.

For both clarity and conciseness considerations, both forward and adjointstate problems have been written down in residual forms.

The numerical applications enabled to exhibit strongly dierent behaviors with respect to the ow regime. Particularly, ows with Re=10, Re=100, and Re=1,000 have been tested. For the specic case of Re=1,000, it has been shown that the presence of the recirculation zone went against a good convergence of the optimization process. As such, a modication of the initial guess topology, with the introduction of dierent small solid zones inside the medium, leads to a better convergence of the optimization problem.

Then, a study of the most important physical parameters (Damkholer and Peclet numbers) involved in the chemical reaction was done. As this problem involves multiple coupled phenomena (advection, diusion and reaction), their inuence on the forward and optimization problem was important to be discussed. First, the reaction rate has been modied for the study on the Damkholer number. Even if the reactant was consumed faster, it was still possible for the algorithm to suggest an ecient design to minimize the cost function. Second, about the Peclet number, the choice was made for the rst numerical examples to work with an advection-dominated problem, compared to mass diusion, with a high Peclet number, in order for the velocity ow to have an important impact on the eciency of the reaction. However, it was shown that the algorithm was also working ne with a lower Peclet number, where diusion was not negligible anymore. In this case, a very regular symmetry has been achieved by the optimized uid/solid distribution, participating to a substantial decrease of the outlet concentration.

A crosscheck has also been performed in order to evaluate the results, and enhance the reliability of the topology optimization method. Although a local minimum was found for Re=10, the comparisons with Re=100 and Re=1000 has been successful as the best performance was achieved with the corresponding optimized geometry.

In further perspectives, the local minimum treatment is an interesting issue to solve, despite the fact that it can be dicult to avoid them while using a gradient-based optimization method. Nevertheless, one possibility would be to use a multi-scale resolution approach. Such a multi-resolution approach, formalized by Liu [START_REF] Liu | A multiresolution method for distributed parameter estimation[END_REF] for solving some ill-posed inverse problems, has been later on successfully used by Dubot et al. [START_REF] Dubot | A wavelet multiscale method for the inverse problem of diuse optical tomography[END_REF] and Liu [START_REF] Liu | A wavelet multiscale method for the inverse problem of a nonlinear convectiondiusion equation[END_REF], for example. In the framework of reactive uid ow topology optimization, the use of a multi-scale resolution approach would enable to perform a scale-byscale optimization on successive convex cost functions, following the work of Chavent [START_REF] Chavent | Nonlinear least squares for inverse problems: theoretical foundations and step-by-step guide for applications[END_REF]. Moreover, added to this interesting property, this would also save computational time due to faster convergence. This point will be addressed in a future research.

( 1 )

 1 to (3), and combining the three physical parameters, namely the uid viscosity, ν, the mass diusivity of the chemical species, D, and the maximal reaction rate, k, one can recover the classical dimensionless version of the initial problem:

on g. 1 ,Figure 1 :

 11 Figure 1: Schematic representation of the optimization problem conguration.

5. 1 . 1 .Figure 2

 112 Figure2presents the results of the optimization process along with the iterations, for the rst case, with Re =10. On the top and left-hand side is presented the evolution of both the cost function value (mean of u × c out ), and of the pressure drop ratio ∆p ∆pmax , along with the iteration count. On the right-hand side, is presented the evolution of the topology: T

Figure 2 :Figure 3

 23 Figure 2: Re=10. Top left: evolution of the cost function and of the pressure drop ratio along with the iteration count. Top right: the topology at iterations 0, 30, 50, and 900. Middle: module of the velocity eld. Left: guess conguration (full uid). Right: optimized topology. Bottom: concentration eld. Left: guess conguration (full uid).Right: optimized topology. 26

Figure 3 :

 3 Figure 3: Re=100. Top left: evolution of the cost function and of the pressure drop ratio along with the iteration count. Top right: the topology at iterations 0, 20, 50, and 390. Middle: module of the velocity eld. Left: guess conguration (full uid). Right: optimized topology. Bottom: concentration eld. Left: guess conguration (full uid).Right: optimized topology. 28

Figure 4 :

 4 Figure 4: Re=1000. Velocity (left) and concentration (right) elds for the full uid topology.

Figure 5 :

 5 Figure 5: Re=1000. Top left: evolution of the cost function and of the pressure drop ratio along with the iteration count. Top right: the topology at iterations 0, 50, 150, and 250. Middle: module of the velocity eld. Left: guess conguration. Right: optimized topology. Bottom: concentration eld. Left: guess conguration. Right: optimized topology.

Figure 6 :

 6 Figure 6: Da = 4 × 10 -5 . Top left: evolution of the cost function and of the pressure drop ratio along with the iteration count. Top right: the topology at iterations 0, 100, 200, and 340. Middle: module of the velocity eld. Left: guess conguration (full uid). Right: optimized topology. Bottom: concentration eld. Left: guess conguration (full uid). Right: optimized topology. 34

2 c

 2 A < 0.05 c A > 0.05 -Da=0.1

Figure 7 :

 7 Figure 7: Comparison of optimization cases for Da=0.1 and Da=0.2: optimal topologies (left) and concentration levels (right)

Table 1 :

 1 Cost function values (×10 -3 ) for dierent congurations the cost function for each Reynolds number is written in bold characters,

	2.6

Table 2 :

 2 Optimization results for dierent congurations

		1	10	2000	1.9	2.60	0.53	98.72%
	2	0.1 100 2000 0.59	10.40	0.61	98.26%
	3	0.1 1000 2000 0.49	9.08	0.33	98.85%
	4	0.2 1000 2000 0.012	62.50	0.306	98.75%
	5	0.1 1000 100	0.72	2.36	1.03	97.1%

Acknowledgments

The authors thank the French Ministry of Higher Education and Research for funding this research. The authors also sincerely thank the reviewers for their highly constructive and wise remarks that enabled us to improve drastically the quality of this paper.

Algorithm 1: General topology optimization algorithm

Input: Level-set function Ψ (0) ; Topology T (0) through α (0) (see eq. ( 10)), along with the smooth version eq. ( 12) ; while criterion [START_REF] Wang | A level set method for structural topology optimization[END_REF] not satised do

Compute the Boltzmann variables f and g solving eq. ( 31); Compute the cost function value J + (f , g) with eq. ( 35); Compute the adjoint Boltzmann variables f * and g * solving eq. ( 39); Compute the cost function gradient ∇ Ψ J + from eq. ( 41); Update of the geometry : actualization of the level-set function with eq. ( 9); return Optimal topology T ( * ) ; Right: optimized topology.