Model Hierarchy for the Shape Optimization of a Microchannel Cooling System

We model a microchannel cooling system and consider the optimization of its shape by means of shape calculus. A three-dimensional model covering all relevant physical effects and three reduced models are introduced. The latter are derived via a homogenization of the geometry in 3D and a transformation of the three-dimensional models to two dimensions. A shape optimization problem based on the tracking of heat absorption by the cooler and the uniform distribution of the flow through the microchannels is formulated and adapted to all models. We present the corresponding shape derivatives and adjoint systems, which we derived with a material derivative free adjoint approach. To demonstrate the feasibility of the reduced models, the optimization problems are solved numerically with a gradient descent method. A comparison of the results shows that the reduced models perform similarly to the original one while using significantly less computational resources.


Introduction
For small devices, such as chemical microreactors and electronic equipment, it is critical to have a cooling system that is able to absorb a lot of heat over a small surface area since the performance of these devices is directly related to their operating temperature. For this purpose, cooling systems based on microchannels have been used (see, e.g., [24,32,42] and the references therein). Their heat transfer coefficient is large due to the high specific surface area of the microchannels and, as they are rather small, they do not increase the overall size of the heat-emitting device too much. Furthermore, the effectiveness of the cooling system critically depends on the uniform distribution of the coolant among all channels. Otherwise, localized zones of high temperature, so-called hot spots, can occur and potentially damage the device.
The purpose of this paper is to investigate the shape optimization of such a microchannel cooling system. In the literature, this task already received a lot of attention, e.g., in [3,13,19,36,47,51,54,60]. In most of these publications, geometrical properties, such as the shape of the boundary or the topology of the geometry, are parametrized and the optimization of the design is carried out over these parameters, leading to finite dimensional optimization problems. This approach suffers from the obvious drawback that only shapes representable by the parametrization can be reached in the optimization. A more general approach for design optimization consists of using shape or topological sensitivity analysis. These techniques are based on the so-called shape derivative, which measures the sensitivity of a shape due to infinitesimal deformations, and the topological derivative, which measures the sensitivity of a geometry with respect to the insertion of an infinitesimally small hole, see, e.g., [14,61] for shape calculus and [45] for topological sensitivity analysis. In recent years these techniques have been applied to many industrial problems, e.g., the shape design of polymer spin packs [27,[37][38][39], electric motors [20,21], acoustic horns [5,57], automobiles [18,46,48], aircrafts [41,55,56] or pipe systems [25,28,58]. To the best of our knowledge, the optimization of a microchannel cooling system by means of shape calculus has only been investigated in our earlier work [6], where we rigorously analyzed the theoretical aspects of this problem.
To model the cooling system mathematically, we introduce the following models: A three-dimensional model representing the most important physics as well as three reduced models. One of them is a threedimensional porous medium model based on a homogenization of the domain, similar to those used in, e.g., [12,33,34]. For the other two reduced models we use a dimension reduction technique similar to the one of [7] that transforms the previous, three-dimensional models to two-dimensional ones. A numerical comparison of the reduced models with the original one shows that they capture the most important physical effects properly, while reducing the computational resources needed considerably.
For the shape optimization, we introduce a cost functional based on the absorption of heat and the uniform distribution of flow through the microchannels. This is adapted to all reduced models using analogous techniques. Subsequently, we present the shape derivatives and adjoint systems for all optimization problems, which we obtained using the material derivative free Lagrange approach of [62]. We use these to solve the optimization problems numerically with the help of a gradient descent method, which is again adjusted to all models. A comparison of the optimized geometries suggests that the models give similar results for the shape optimization problem while being substantially more efficient.
This paper is structured as follows: We introduce our three-dimensional model of the cooling system as well as the shape optimization problem in Section 2. In the three subsequent sections we then present the reduced models respectively. Additionally, the corresponding adaptations for the shape optimization problem are discussed. The dimension reduction technique is introduced and used to derive a twodimensional model in Section 3. Section 4 provides the details for the three-dimensional porous medium model of the cooler, to which the dimension reduction approach is then applied in Section 5, so that we also get a two-dimensional porous medium model. The accuracy of the reduced models is compared numerically in Section 6. Finally, we discuss the implementation as well as the results of the numerical shape optimization in Section 7, where we again focus on comparing the different models to each other.

Problem Formulation
First, we describe the geometry of the cooling system and the three-dimensional mathematical model as well as the shape optimization problem. Subsequently, we give a short introduction to shape calculus and then present the shape derivative and adjoint system for the optimization problem.

Description of the Geometry.
To model the cooling system we consider only the domain of the coolant Ω ⊂ R 3 , which is assumed to be encased by metal. This metal conducts the heat from the heat source, to which the entire cooler is attached, to the coolant. For simplicity, we call Ω the geometry or domain of the cooling system throughout the rest of this paper. The domain Ω has the structure Ω =Ω × (0, h), whereΩ ⊂ R 2 is a two-dimensional domain and h > 0 is the constant height of the geometry, i.e., Ω is an extrusion ofΩ along the z-axis. The boundary of Ω is denoted by Γ and is divided into three parts: The inlet Γ in , where the coolant enters the system, the wall boundary Γ wall , where heat transfer from the heat source to the cooling system takes place, and the outlet Γ out , where the cooling liquid leaves the domain. The boundaries inherit the structure of Ω, i.e., Γ in =Γ in × (0, h), The structure of Γ wall comes from the fact that the planes { z = 0 } and { z = h }, i.e., the top and bottom of the cooler and part of Γ wall , cannot be represented as boundaries ofΩ. Of course, it holds thatΓ =Γ in ∪Γ wall ∪Γ out . Additionally, we denote by Ω mc the subdomain corresponding to the microchannels, whose wall boundaries are denoted by Γ mc ⊂ Γ wall . This situation is depicted in Figure 1, where the domainΩ is shown. Note, that as  the height of the geometry h is very small compared to the diameter ofΩ, we only show slices through Ω at z = h /2 for all three-dimensional problems. Finally, we remark that such cooling systems can be manufactured, e.g., using wet chemical etching [35] or 3D-printing [50].

Mathematical Model.
As we only investigate the steady state of the system, we model all physical processes using stationary equations. The viscous dominated flow of the cooling fluid, which is caused by the slow flow velocities and the small size of the geometry, is modeled by the Stokes equations where u and p denote the fluid velocity and pressure, respectively. Furthermore, µ denotes the dynamic viscosity of the coolant and ∂ n v is the normal derivative of a function v given by ∂ n v = Dv n, where Dv is the Jacobian of v and n is the outward unit normal on Γ. We model the flow of the coolant into Ω with the fixed inflow condition u = u in on Γ in . Here, u in is the inflow velocity which is chosen to be in the direction of the inward facing normal and for the sake of simplicity we assume that u in ∈ H 1 (Ω) 3 with u in = 0 on Γ wall , as in our theoretical analysis [6]. For the wall boundary Γ wall we use the no-slip condition and on the outlet Γ out we use the do-nothing condition that models the unimpaired outflow of the fluid (see, e.g., [29]). The temperature of the coolant changes due to both conduction and convection. This is modeled by the convection-diffusion equation in Ω, where T denotes the fluid's temperature and u is its velocity, i.e., the solution of (2.1). Furthermore, κ, ρ, and C p denote the coolant's thermal conductivity, density, and specific heat capacity, respectively. The temperature of the inflowing fluid is fixed to T in on Γ in by a Dirichlet condition. Analogously to before, we assume that T in ∈ H 1 (Ω) (cf. [6]). The heat transfer over the wall boundary is modeled by the Robin boundary condition on Γ wall , where α denotes the heat transfer coefficient between coolant and wall. For simplicity, we assume that the heat source generates a temperature distribution that yields a constant temperature T wall on Γ wall , but non-constant temperature distributions could also be considered. The behavior of the temperature at Γ out is modeled by a homogeneous Neumann condition. The values of the physical parameters for our particular problem setting are given in Table 1. They are chosen to resemble a viscous coolant that could be used for the application in a chemical microreactor. The heat transfer coefficient is chosen so that the amount of energy absorbed by the system resembles the amount of heat generated by the chemical reaction. In particular, we consider a setting similar to the one described in [10], where the Sabatier process in microreactors is investigated. As the notion of strong solutions of PDEs is usually too strict for the means of PDE constrained optimization problems (see, e.g., [26,65]), we consider equations (2.1) and (2.2) in their weak form. To this end, we define the (affine) Sobolev spaces The weak form of (2.1) and (2.2) is then given by If Ω has a Lipschitz boundary and the data is sufficiently smooth, it is not too difficult to show that problem (2.4) has a unique weak solution that depends continuously on the data (see, e.g., [6,16,53]).

The Optimization Problem.
Our goal is to improve the cooling system by means of shape optimization, i.e., we want to optimize it by changing its shape without altering its topology. To do so, we assume that we can only change the shape of Γ wall and that Γ in and Γ out remain fixed. This is reasonable since the cooling system is connected to other parts via in-and outlet. Furthermore, we also do not change the underlying structure of the geometry since we assume that the height h of the cooler is fixed, and we also only allow the microchannels to change in length.
The quality of the cooling system is measured by the following criteria. First, it should absorb a specific amount of heat so that a finely detailed cooling of the heat source is possible. Second, we want to find a geometry that distributes the coolant uniformly to the microchannels as this helps preventing hot spots. To model this, we define the cost functional as For the term J 1 we have the following considerations. The heat flux on Γ wall is given by −κ∇T as the fluid velocity vanishes there. Thus, the overall energy entering the domain is given by Γ wall κ∇T · n ds. Using the Robin boundary condition on Γ wall (cf. (2.2)), we can express the amount of energy entering the domain via Q(Ω, U ). Therefore, by minimizing J 1 we try to get a cooling system that absorbs a particular amount of heat, given by Q des . This is particularly useful for cooling devices that rely on a specific operating temperature, such as chemical reactors.
The term J 2 has the purpose of minimizing the distance of the fluid velocity u to some desired velocity u des in L 2 (Ω mc ), where Ω mc is the domain corresponding to the microchannels (cf. Section 2.1). In particular, we choose u des as the velocity profile corresponding to the case of uniformly distributed flow through each channel, which can be computed as Poiseuille flow for a channel with square base (see, e.g., [11,Chapter 3.4]). This term of the cost functional leads to a uniform flow distribution which helps to minimize the occurrence of hot spots.
Finally, the term J 3 is a so-called perimeter regularization (see, e.g., [59,62]), that penalizes the increase of surface area. This is used since it can grant the existence of a minimizer, as discussed in [1].
As we assume that the state system (2.4) has a unique weak solution for every domain Ω with a Lipschitz boundary, we now introduce the so-called reduced cost functional j by j(Ω) = J(Ω, U (Ω)), (2.7) where U (Ω) denotes the solution of (2.4) in Ω. With this, it is evident that (2.6) is equivalent to the reduced optimization problem min where the PDE constraint is formally eliminated.

Shape Calculus.
To solve the shape optimization problem numerically, we calculate the shape derivative of the reduced cost functional (2.7), which is then used in a gradient descent method (see also Section 7). For a detailed introduction to shape calculus we refer to the monographs [14,61]. To compute the shape derivative, we utilize the so-called speed method, which is also known as velocity method. It uses the flow of a vector field V to deform a domain Ω ⊂ D, where D ⊂ R d is a so-called hold-all domain, i.e., a domain that contains all domains under consideration, and d denotes the dimension of the geometry. For a given vector field V ∈ C k 0 (D; R d ) with k ≥ 1, i.e., the space of all k-times continuously differentiable functions from D to R d having compact support in D, the evolution of a point x ∈ Ω under the flow of V is given by the solution of the ODĖ We know that this system has a unique solution x(t) for t ∈ [0, τ ] if τ > 0 is sufficiently small. Hence, we define the flow of V as the diffeomorphism Now we define the shape derivative as in [62], where we write 2 D := { Ω | Ω ⊂ D } for the power set of D.
if it exists. The function J is called shape differentiable at Ω ∈ S if the above limit exists for all V ∈ C ∞ 0 (D; R d ) and if the mapping , is linear and continuous. Then, dJ(Ω)[V] is called the shape derivative of J at Ω in direction V.

Formal Shape Optimization.
Our goal is to calculate the shape derivative of the reduced cost functional (2.7), which is done with an adjoint approach (see, e.g., [14,62] for more details). The rigorous verification of the assumptions for this approach is beyond the scope of this paper and can be found in our earlier work [6]. We define the Lagrangian associated to (2.7) by where U = (u, p, T ) ∈ U and P = (v, q, S) ∈ P. In particular, we have L(Ω, U (Ω), ψ) = j(Ω) for all ψ ∈ P as U (Ω) is the solution of the state system (2.4). Further, for a sufficiently smooth vector field V with associated flow Φ t , we define Ω t := Φ t (Ω). With this, we introduce the associated shape Lagrangian by which corresponds to the original Lagrangian L on the perturbed domain Ω t . We denote by U t = U (Ω t ) the solution of the state system on the domain Ω t and define U t = U t • Φ t . With this, a similar argumentation as in [63] shows that we have G(t, U t , ψ) = j(Ω t ) for all ψ ∈ P, because Φ t is a diffeomorphism. Hence, we can calculate the shape derivative of the reduced cost functional using Additionally, under the assumptions of [62, Theorem 3.1] it holds that where U is the solution of the state system, i.e., (2.4), and P is the solution of the adjoint system It is straightforward to see that this adjoint system is equivalent to is the weak solution of the state system (2.4). To calculate the shape derivative using (2.9), we use the transformation formula to rewrite G(t, U, P ), where we pull-back the integrals to the domain Ω. To this end, we define denotes the inverse of DΦ t , and not the Jacobian of Φ −1 t . Additionally, we note that for all sufficiently smooth functions f and vector fields V for t ∈ [0, τ ] with τ > 0 being sufficiently small (see, e.g., [14,62]). With this, we rewrite G as follows (2.11) Further, due to [14,62] we have the following derivatives is the symmetric part of the Jacobian of V and div Γ denotes the tangential divergence defined as div Γ (V) = div (V) − DV n · n. Finally, differentiating (2.11) w.r.t. t and using (2.9) proves the following theorem (cf. [6]).

Theorem 2.2. Under the assumptions of [62, Theorem 3.1] the shape derivative of
where (u, p, T ) = U (Ω) is the solution of the state system (2.4) and (v, q, S) = P (Ω) is the solution of the adjoint system (2.10).

The Dimension Reduction Technique
For our first reduced model, we introduce a dimension reduction technique, similar to the idea discussed in [7] for the topology optimization of Stokes flow which was also used, e.g., in [22,66]. The technique transforms the three-dimensional state system to a two-dimensional one, making its numerical solution considerably easier. Afterwards, we discuss this technique in the context of the shape optimization problem.

Description of the Model.
The main idea of this model is to exploit the structure of our domain, i.e., we use that Ω =Ω × (0, h). As this model is based on the weak formulation of the PDEs, we now introduce the following two-dimensional (affine) function spaces in analogy to the ones given in (2.3) The key observation for the dimension reduction of the Stokes equation is the following. For viscous flow between parallel, infinite planes we get a parabolic velocity profile in analogy to two-dimensional Poiseuille flow (cf. [11]). This situation is not very different from our setting. We still consider the flow of a viscous fluid between two parallel plates, albeit the geometry is more complex. Due to this observation we make the following assumption: there is no fluid velocity in z-direction and the x and y components of the velocity are parabolic in z, i.e., we describe the fluid velocity The factor 4 /h 2 is chosen so thatũ represents the maximum fluid velocity. Furthermore, we see that the no-slip boundary conditions at { z = 0 } and { z = h } are satisfied by (3.1). For the fluid temperature and pressure we assume that they are constant along the z-direction, i.e., wherep ∈P andT ∈W . For the pressure this is reasonable as we do not have a flow in z-direction (cf. (3.1)) and, hence, we easily get ∂ z p = 0 from the Stokes equation (2.1). For the temperature the situation is more complicated and, in general, it is not constant in z. However, as our problem is convection-dominated and there is no convection in z-direction due to (3.1), we observe that the variation in temperature is significantly higher in the x-y-plane than it is in z-direction, and the assumptions (3.2) are justified. Naturally, we assume thatũ in andT in are obtained from u in and T in analogously to (3.1) and (3.2). As we want to derive a two-dimensional weak formulation, we introduce the functionŝ in analogy to (3.1) and (3.2). Using these functions in the weak form (2.4) and performing the integration w.r.t. z yields the twodimensional system 3) Note, that we useT in = T in for the inflow temperature as this is assumed to be constant. For the inflow velocity we choose a parabolic profile like in (3.1) so that we are compatible with both the dimension reduction and the no-slip boundary condition on the top and bottom of the geometry.

Application to the Shape Optimization Problem.
We now describe the application of the dimension reduction technique to the shape optimization problem (2.6). For this, we again use the parabolic profile (3.1) for the velocity as well as the constant profile (3.2) for pressure and temperature. Note, that we also use the parabolic profile for u des so that this state is reachable by the dimension reduction model. Substituting these into the cost functional (2.5) and carrying out the integration over z yields the dimension-reduced cost functional Note, that the terms Ω 2 dx and Ω 2α(T wall −T ) dx arise due to the integration over lid and bottom of the three-dimensional domain. The corresponding reduced optimization problem reads To derive the shape derivative for this model, we have to make an additional assumption on the vector fields we use to deform the domain. As stated in Section 2.3, we want to keep the height of the cooling system fixed. Therefore, we restrict ourselves to deformations generated by vector fields whose z-component vanishes, i.e., . Any other vector field with non-vanishing z-component only introduces a reparametrization along the z-axis that cannot change the geometry as h is fixed. Using this, we apply the dimension reduction technique to the shape Lagrangian defined in (2.8) and, after lengthy calculations, end up with the shape derivative is the weak solution of (3.3) and (ṽ,q,S) =P (Ω) solves the adjoint system for allṼ = (ũ,p,T ) ∈P.
Remark. Note, that the only simplification made for the dimension-reduction model is that we suppose that the state variables are of the form (3.1) and (3.2). Hence, we could derive the shape derivative given in (3.5) also by applying the dimension-reduction technique directly to the adjoint system (2.10) as well as to the shape derivative (2.12). Thus, the optimization commutes with the dimension reduction technique.

Porous Medium Model
Both previously described models completely resolve all microchannels which complicates a physically meaningful discretization of the geometry. To circumvent this, we now introduce our second reduced model, for which we model the microchannels as a porous medium. This allows us to use a much simpler geometry which substantially eases the numerical solution of the equations. For a detailed introduction to porous medium models and their application to microchannels we refer, e.g., to [30,43] and [12,33,34], respectively. 4.1. Description of the Geometry for the Darcy Model. First, we introduce the structure of the domain Ω por that we consider for the Darcy model. It consists of the two subdomains Ω por f = Ω \ Ω mc , the domain without the microchannels that remains unaltered, and Ω por d , which is the homogenized geometry corresponding to Ω mc . The boundary of Ω por is denoted by Γ por and consists of the following parts: For the in-and outlet we have Γ por in = Γ in and Γ por out = Γ out , and the wall boundary is given as Γ por wall = Γ wall \ Γ mc . Additionally, we obtain a new boundary Γ por d that is the outer boundary of Ω por d , and the interfaces between Ω por f and Ω por d are denoted by Γ por fd . Naturally, the geometry of the Darcy model has a similar structure to the geometry of the full model, i.e., we have Ω por =Ω por × (0, h). An analogous decomposition also holds for the two subdomains as well as the boundaries (cf. Section 2.1). The corresponding two-dimensional domainΩ por depicting the above situation can be seen in Figure 2.

Description of the Model.
To couple the fluid equations on Ω por f , where we have the usual Stokes system as in Section 2, and Ω por d , where we consider the porous medium, we choose the Brinkman equation (cf. [43]). This allows an implicit coupling through transmission conditions that vanish in the weak formulation of the problem (cf. (4.5)). In particular, the Brinkman equation, without the transmission conditions, reads where K denotes the permeability tensor and µ denotes the effective viscosity. Here, u describes the averaged (or Darcy) velocity and not a physical fluid velocity. The permeability tensor K already models the friction arising from the no-slip boundary condition in the microchannels completely so that we must not prescribe a no-slip condition on Γ por u · n = 0 and µ ∂ n u × n = 0 on Γ por d so that the fluid cannot leave the domain and does not experience additional friction. For the permeability tensor we observe that the flow through the channels is basically one-dimensional and in y-direction. To model this anisotropy, K has the structureK · diag(ε, 1, ε), wherê K is the permeability in y-direction and 1 ε > 0 is a relaxation parameter needed for the invertibility of K (cf. (4.1)). We computeK numerically using Darcy's law (see, e.g., [30,43]) and choose ε = 1e−5. Further, as 1 /K is very large, we neglect the influence of the effective viscosity and set µ = µ as originally suggested in [8].
To describe the transmission conditions for the coupling of (4.1) with the Stokes system, we denote by v f and v d the restriction of a function v on Ω por f and Ω por d , respectively. Then, the coupling conditions are given by the continuity of both the velocity and the normal component of the viscous stress tensor, i.e., wheren denotes the outer unit normal to Ω por f . To summarize, the Darcy model for the fluid reads For the modeling of heat transfer in the porous medium we modify the so-called local thermal nonequilibrium approach that can be found, e.g., in [43,Chapter 2.2]. This model uses two equations, one for the liquid and one for the solid phase. However, we assume that the solid's temperature is constant, in analogy to the constant wall temperature T wall (cf. Section 2), and are left with the equation for the fluid phase. In particular, the heat transfer equation in Ω por d reads −∇ · (ϕκ∇T ) + ρC p u · ∇T + h fs (T − T wall ) = 0 in Ω por d , ϕκ ∂ n T = 0 on Γ por d , where ϕ is the porosity of the porous medium, i.e., the fraction of total volume occupied by the fluid, and h fs denotes the interfacial heat transfer coefficient. The latter can, again, be computed numerically using the formulas given in, e.g., [49,64]. The homogeneous Neumann boundary condition in (4.4) is used for a similar reason as the slip boundary condition in (4.1). The heat transfer between solid and fluid phase is completely contained in the term h fs (T − T wall ) and, therefore, we must not have an additional heat source on the boundary Γ por d . The parameters for the porous medium model which we computed numerically using a single microchannel as representative elementary volume are given in Table 2. κ ∂nT f − ϕκ ∂nT d = 0 on Γ por fd . Note, that demanding the continuity of the heat flux' conductive part over the interface is already sufficient for the continuity of the entire flux, as we get the continuity of its convective part directly from   × P por × W por 0 , where we assume that u in ∈ H 1 (Ω por ) 3 with u in = 0 on Γ por wall and u in · n = 0 on Γ por d as well as T in ∈ H 1 (Ω por ) in analogy to Sections 2 and 3. Finally, the weak form of the Darcy model is given by  The difference now lies in the interpretation of the velocity. Whereas we have a physical velocity in Ω mc and, thus, also prescribe a physical one for u des , the fluid velocity in Ω por d is an averaged one. Hence, the desired velocity u por des has to be the mean of u des over Ω por d . The perimeter regularization is modeled analogously to before, i.e., we use Finally, the cost functional J por for the Darcy model reads J por (Ω por , U ) = λ 1 J por 1 (Ω por , U ) + λ 2 J por 2 (Ω por , U ) + λ 3 J por 3 (Ω por , U ). As before, we define the reduced cost functional by j por (Ω por ) = J por (Ω por , U (Ω por )), where U (Ω por ) denotes the solution of (4.5) on Ω por and the reduced optimization problem is given by min Ω por j por (Ω por ).
Applying the same techniques as in Section 2.4 we derive the following shape derivative for j por where (u, p, T ) = U (Ω por ) is the solution of (4.5) and (v, q, S) = P (Ω por ) is the solution of the following adjoint system

Dimension Reduction Applied to the Darcy Model
For our final reduced model, we apply the dimension reduction technique of Section 3 to the porous medium model derived in the previous section.

Description of the Model.
For the application of the dimension reduction approach, we now use two different profiles for the fluid velocity.
In Ω por f we have a physical velocity and, hence, we use a parabolic profile as in (3.1). In contrast, in Ω por d we have an averaged velocity and, therefore, we assume to have a constant profile for the velocity there. I.e., we use Note, that we use the scaling factor of 6 /h 2 so thatũ = [ũ 1ũ2 ] T represents the mean fluid velocity in both parts of the domain, making a coupling with transmission conditions similarly to (4.2) possible. For the pressure and temperature we use the same approach as in (3.2) and assume that they are constant in z.
We introduce the (affine) Sobolev spaces whereũ in ∈ H 1 (Ω por ) 2 withũ in = 0 onΓ por wall andũ in · n = 0 onΓ por d as well asT in ∈ H 1 (Ω por ), similarly to before. Proceeding analogously to Section 3, we get the following weak form for the 2D Darcy model for allṼ = (ṽ,q,S) ∈P por , where we rescaled the equations to obtain the continuity of the pressure over the interfaceΓ por fd .

Application to the Optimization Problem.
For the shape optimization problem we again proceed as before and consider the cost functional The corresponding (reduced) shape optimization problem reads miñ Ω porj por (Ω por ) =J por (Ω por ,Ũ (Ω por )), whereŨ (Ω por ) denotes the solution of (5.1). As in Section 3, we assume that the z-component of the vector field V vanishes (cf. (3.4)) in order to calculate the shape derivative, which then reads where (ũ,p,T ) =Ũ (Ω por ) is the solution of (5.1) and (ṽ,q,S) =P (Ω por ) solves the adjoint system for allŨ = (ũ,p,T ) ∈P por .

Numerical Comparison of the Reduced Models
After introducing all reduced models in Sections 3 to 5, we now investigate how well they approximate the original model from Section 2. To distinguish between the models we use the following naming: we call the model from Section 2, i.e., the three-dimensional model without further approximations, the "full 3D" model. In analogy, we call the model from Section 3, i.e., the one arising from applying the dimension reduction technique to the full 3D model, the "full 2D" model. The models of Section 4 and 5 are then called "Darcy 3D" and "Darcy 2D" model, respectively.

Discretization and Numerical Solution of the PDEs.
The computational meshes are generated with the help of FreeCAD 0.16 [52] and Gmsh 4.1.0 [23]. To get a comparable discretization in the x-y plane, we generate the three-dimensional meshes by extruding the ones for the two-dimensional models. Further, for the Darcy models we use a similar discretization on the in-and outlet domains as for the full models, whereas the homogenized part of the domain Ω por d is discretized significantly coarser than Ω mc to reduce computational cost. We discretize all PDEs with the finite element method using FEniCS 2018.1 (cf. [2,40]). For the fluid velocity we use quadratic Lagrange elements, and for both the pressure and temperature we use linear Lagrange elements, i.e, we use the LBB-stable Taylor-Hood elements for the Stokes system (see, e.g., [15,29]). For the temperature we additionally use a streamline-upwind Petrov-Galerkin (SUPG) method to stabilize the convection-dominated system (see, e.g., [9,15,29]).
The linear systems arising from the finite element method are solved with the library PETSc, version 3.10.5 (cf. [4]). For all two-dimensional problems we use the direct solver MUMPS. For the threedimensional Stokes equations we use GMRES as it showed a better convergence than the MINRES method. As preconditioner, we use a FIELDSPLIT method based on the Schur complement, consisting of an ILU preconditioner for the velocity block and the algebraic multigrid preconditioner BOOMERAMG for the pressure block. For the three-dimensional convection-diffusion equation we again use GMRES with an ILU factorization as preconditioner. All linear systems are solved to a relative tolerance of 1e−10.
In Table 3 we compare the computational efficiency of the models w.r.t. the size of the mesh and linear systems, as well as the resources needed for the solution of the corresponding PDEs. We see that all reduced models lead to smaller meshes and, thus, also to smaller linear systems. Especially the dimension reduction yields a significant decrease in system size that comes from the fact that we use     Table 6. Errors in the L 1 (Ω) norm.
quadratic elements for the velocity and that one velocity component vanishes for the 2D models. Further, both the solution time and memory requirements also go down quite a bit from the full 3D model to the 3D Darcy one, and they decrease substantially when considering the two-dimensional models.

Comparison of the Models.
To compare the models we proceed as follows. First, we solve the PDEs on their respective domains. For the full 2D model we then only have to "reverse" the dimension reduction using the profiles given in (3.1) and (3.2) to obtain three-dimensional quantities. For the 3D Darcy model we do not alter the solution in Ω por f , since we did not change the model there, and compute from the mean velocity in Ω por d the corresponding physical velocity as solution of Poiseuille flow in the channels. For the pressure and temperature we interpolate the values on Ω por d to Ω mc as they describe values extended to the larger averaged domain Ω por d . Finally, we combine both approaches for the 2D Darcy model.
We computed both the absolute and relative errors of the reduced models to the full 3D model in three different norms: The L ∞ (Ω)-norm (Table 4), the L 2 (Ω)-norm (Table 5), and the L 1 (Ω)-norm (Table 6), and the relative errors are also shown in Figure 3 in a logarithmic plot. Note, that the relative errors of the temperature are computed w.r.t. the reference temperature T in , i.e., we use where T is the temperature computed by the full 3D model andT is the temperature obtained from a reduced model. The largest difference between all models can be seen for the L ∞ (Ω)-norm. There, the errors in velocity and temperature for the full 2D model are ten times as small as the ones for both Darcy Full 2D Darcy 3D Darcy 2D Figure 3. Comparison of the models' relative errors.
models. This comes from the fact that the full 2D model is posed on the domainΩ that still includes the microchannels, whereas the Darcy models utilize an averaged domain. In particular, the flow of the coolant into and out of the channels is only resolved by the full 2D model, which explains the large L ∞ error for the other ones. However, this is remedied when we consider the L 2 and L 1 norms. There, we observe that the Darcy models are still worse compared to the full 2D model, but the difference of the errors is now considerably smaller. Further, we observe that the pressure is approximated very well by all models since the corresponding relative errors are below 1 % in all considered norms. Altogether, except for the L ∞ norm, all relative errors are below 5 % for all quantities, which indicates that our reduced models work rather well.
In addition to this, we see that both Darcy models have very similar errors in all considered norms. This suggests that most of the error of the 2D Darcy model comes from the modeling as a porous medium, and not from the subsequent dimension reduction. From this, we conclude that we can use the 2D Darcy model instead of the 3D one as their errors are nearly identical, but the 2D model is significantly more efficient (cf. Table 3). Finally, the full 2D model shows the best overall performance of the reduced models due to its combination of efficiency and accuracy.

Numerical Results for the Optimization Problem
Now, we describe the numerical solution of shape optimization problems using (2.6) as an example and then explain the details and modifications we use for the reduced models. Afterwards, we discuss the results obtained from the different models and compare them.

Numerical Solution of the Shape Optimization Problem.
In the previous sections we presented the shape derivatives for our optimization problems. However, for numerical purposes this is not yet sufficient as the shape derivative only gives us the sensitivity of the cost functional for a given deformation generated by the flow of a vector field V. In particular, the shape derivative dj(Ω)[V] is a linear functional acting on V. To solve the shape optimization problem numerically, we compute a descent direction V s , i.e., a vector field whose associated flow gives a descent in the cost functional, from the shape derivative. To do so, we choose a symmetric, continuous and coercive bilinear form a : H(Ω)×H(Ω) → R, where H(Ω) denotes some suitable Hilbert space on the current geometry Ω which is specified later on. Then, the shape gradient G w.r.t. a is defined as the solution of the variational equation due to the coercivity of a. Therefore, we use the negative shape gradient in a line search method to solve problem (2.6). After computing the search direction V s = −G, we deform the domain numerically with the so-called perturbation of identity given by Starting from an initial guess, the step size t is accepted if it satisfies the Armijo condition where the last equality holds due to (7.1), and the parameter σ ∈ (0, 1) is chosen to be σ = 1e−4 (cf. [31,44]). If this is not satisfied, the step size is halved and the procedure is repeated. Additionally, we do a mesh quality control based on conditions given in [17] to avoid step sizes that lead to excessively large deformations. We only accept step sizes that satisfy 1 2 ≤ det (I + t DV s ) ≤ 2 as well as t ||DV s || F ≤ 0.3, (7.4) for each element of the mesh, where ||·|| F denotes the Frobenius norm of a matrix. In case the step size is accepted, we double it after the deformation of the domain to get a good initial guess for the next iteration.
For the stopping criterion we make the following considerations. As a is symmetric, continuous, and coercive, we see that it defines an inner product on H(Ω). Hence, for the norm of the shape gradient we use the one induced by this scalar product, i.e., ||G|| H(Ω) := a(G, G) (cf. [17,59]). We denote the shape gradient on the initial geometry Ω 0 by G 0 and terminate the algorithm if the relative stopping criterion ||G|| H(Ω) is satisfied. For all of our numerical experiments we choose the relative tolerance as ε = 1e−3. Additionally, we also stop the optimization after a fixed amount of iterations or in case the conditions (7.3) and (7.4) cannot be satisfied. The general numerical method is summarized in Algorithm 1, where problem (2.6) is used as an example. The only thing left to do before we can apply it is the specification of the bilinear form a for all models, which is done in the subsequent section. Increase the step size for the next iteration: t = 2t 7.2. Choice of the Bilinear Form for the Shape Gradient. The bilinear forms used for computing the shape gradient in (7.1) are based on the equations of linear elasticity which are used, e.g., in [5,27,59]. We modify this approach and use equations of anisotropic, inhomogeneous, linear elasticity for our problems.
Full 3D Model. Let us start with the model from Section 2. We define the Hilbert space as Note, that this choice is consistent for our optimization problem, i.e., it respects the geometrical constraints stated in Section 2.3. For V ∈ H(Ω) we get by the perturbation of identity (7.2) for Ω t = (I+tV)Ω that the boundaries Γ in and Γ out are fixed for all t ∈ [0, τ ]. Furthermore, thanks to the slip condition V · n = 0 on { z = 0 } ∪ { z = h } ∪ Γ mc , the height of the geometry remains fixed and the microchannels can only change their length. For the computation of the shape gradient, we use the bilinear form a : H(Ω) × H(Ω) → R given by Here, λ elas and µ elas are the so-called Lamé parameters and δ ≥ 0 is a damping parameter. Further, C 1 is a numerical constant that leads to an anisotropic strain tensor E(U ), which we choose as C := 1e5. Hence, the z-component of the shape gradient is small in comparison to its other components. This is done to avoid unnecessarily small step sizes due to a deformation in z-direction that does not have an actual effect on the geometry as its height is already fixed. Finally, the term ν(x) models an inhomogeneous stiffness of the geometry which can only be defined after discretization of the geometry with a triangulation T h . Then, it is given by i.e., it is one over the relative (d-dimensional) volume of the considered element. This idea from [5] ensures that large elements have a lower stiffness than small ones since they absorb large deformations better. As in [5] we compute ν once on the initial mesh and do not update it during the optimization process so that elements cannot become arbitrarily stiff. We found that using this approach significantly increased the mesh quality during the optimization and that we did not have to remesh at all.

Full 2D Model.
For the dimension reduction model we proceed similarly to Section 3. We already assumed that the z-component of the deformation vanishes to derive the two-dimensional shape derivative (cf. (3.4)). This is now also used to derive a dimension-reduced formulation of the bilinear form a. To this end, we introduce the Hilbert spacẽ H(Ω) = Ṽ ∈ H 1 (Ω) 2 Ṽ = 0 onΓ in ∪Γ out ,Ṽ · n = 0 onΓ mc .

3D Darcy Model.
For the Darcy model we use the Hilbert space which acts as an equivalent of H(Ω) for the domain Ω por . As before, the perturbation of identity with a vector field from H por leaves both the height h of Ω por and the boundaries Γ por in and Γ por out fixed. We choose the bilinear form and σ(U ) is given in (7.6). As before, the Darcy model coincides with the full model in Ω por f , and in Ω por d we have the following modification. The strain tensor E d in Ω por d now also exhibits a large stiffness in the x-direction. This models the influence of the microchannels on the shape gradient, due to the following reason. As we have lots of narrow channels in Ω mc and use the slip condition V · n = 0 on Γ mc , the geometry of the channels is only allowed to stretch or compress along the y-axis. Therefore, we can neglect the x-component of the shape gradient in Ω mc . The anisotropic strain tensor E d achieves a similar effect by increasing the stiffness in x-direction in Ω por d . 2D Darcy Model. The 2D Darcy model again combines the dimension reduction technique with the porous medium model. As before, we assume that the z-component of the deformation vanishes, and now use the Hilbert spacẽ H por (Ω por ) = Ṽ ∈ H 1 (Ω por ) 2 Ṽ = 0 onΓ por in ∪Γ por out ,Ṽ · n = 0 onΓ por d .
The bilinear form for the 2D Darcy model then reads For all models we choose µ elas = 1 as well as λ elas = 1e−1 for the Lamé parameters. The damping parameter was chosen to be δ = 1e−1. For the numerical solution of (7.1) we again use FEniCS. We discretize the equations using linear Lagrange elements and solve the resulting linear systems using MUMPS for the two-dimensional problems, and a conjugate gradient method preconditioned with BOOMERAMG for the three-dimensional ones.

Numerical Results of the Optimal Shape Design Problem.
After giving the details of the numerical solution of the shape optimization problems as well as the choice of the bilinear forms for computing the shape gradient for all models, let us now investigate the results we obtained. We choose the weights for the cost functionals as where Ω 0 denotes the initial geometry. The scaling of the cost functional is chosen so that we weight the functions J 1 and J 2 equally, while having only a slight regularization from J 3 . For the reduced models we choose an analogous scaling. Finally, we discuss the choice of u des and Q des . For the former, we choose the velocity corresponding to the case of uniformly distributed flow among the microchannels, which we compute numerically as Poiseuille flow. For the latter, we have the following considerations.
The cooling system absorbs 6.74 W of heat with its initial shape. For the optimized cooler we want to increase this and, hence, choose Q des = 7.1 W, which corresponds to an increment of about 5 %.  The results of the shape optimization for the full 3D model can be seen in Figure 4. In Figure 4a the initial and optimized geometries are shown. First, we observe that both the initial and optimized geometry are (nearly) point-symmetric to the center of the geometry. For the optimized shape, we see that the in-and outlet domains are pushed to the outside. Additionally, they are dented in the middle, creating a kind of U-shape (on the top). The length of the channels changed accordingly, to balance the pressure differences generated in the in-and outlet domains. The results of this can be seen in Figure 4b, where the mass flow rate of the coolant in the channels is depicted. We see that for the initial geometry the flow through the channels resembles a U-shape. In particular, the outer channels get the most amount of fluid and the middle ones get the least amount. This discrepancy is removed on the optimized domain, where we observe a nearly uniform flow distribution among all channels, achieving one goal of the optimization. Additionally, the heat absorbed by the cooler increases from 6.74 W on the initial geometry to 7.094 W on the optimized one, nearly reaching the desired amount of 7.1 W. Moreover, the size of the geometry did not increase substantially. Hence, the optimized cooler could be used in place of the initial configuration. As the reduced models perform nearly indistinguishably with regards to these two objectives, we do not show the analogous results for them for the sake of brevity.
The history of the optimization process is shown in Figure 5 for all four models, where both the value of the cost functional and the relative norm of the gradient are shown. We observe that the function value does not decrease much further after five iterations for all models which indicates that all geometries converge very quickly to the optimal one. Note, that the cost functional value is comparatively similar for all four models, again suggesting that they perform similarly. We terminated the optimization algorithm after 20 iterations since the norm of the gradient nearly decreased by two orders of magnitude for all     models, which is sufficient for industrial applications. These results indicate that Algorithm 1 converged and found a (local) minimizer. The optimized geometries for all four models can be seen in Figure 6, where additionally the corresponding temperature distribution is shown. To compare them, we also incorporated the boundaries of the optimized geometry for the full 3D model into the plots for the other models. We observe that all four geometries look very similar. In particular, there are only minor differences between the full 3D and the full 2D model. The full 2D model has some slightly longer or shorter channels as well as slight differences in the outer boundaries of the in-and outlet domains. Compared to the optimized geometry of the full 3D model, both Darcy models show larger deviations. The in-and outlet domains for the Darcy models are consistently further to the outside than their equivalents. Additionally, the curves near the in-and outlet are pushed deeper into the geometry. However, these differences are to be expected due to the changes in modeling. In contrast, the region Ω por d behaves similarly to its counterpart, Ω mc . The curves describing the in-and outlets of the channels are very similar between the full 3D model and the Darcy ones (not visualized). Moreover, comparing the optimized geometries of both Darcy models reveals only very subtle differences. These results suggest that the differences in the optimized geometries for the Darcy models mainly arise from the porous modeling and not from the dimension reduction.
In conclusion, all reduced models also work very well in the shape optimization context. Additionally, since the main numerical work of Algorithm 1 consists of solving the state and adjoint systems, the reduced models demand significantly less computational resources. This is depicted in Table 7  the time for the solution of the shape optimization problems as well as the speedup for the reduced models relative to the full 3D model is shown. Again, we observe that the two-dimensional models are over 100 times faster than the full 3D models, making them particularly attractive. Finally, we note that, as for the state system, the full 2D model shows the best performance of all reduced models due to its combination of accuracy and efficiency.

Conclusion and Outlook
In this work, we introduced four different models for a microchannel cooling system using both porous medium modeling and a dimension reduction technique. A numerical comparison showed that all reduced models approximate the original one quite well while requiring substantially less computational resources. Further, we presented a shape optimization problem based on heat absorbed by the cooler and the uniform distribution of coolant among the microchannels which we adapted to all reduced models. For all our models, we presented both the shape derivative and the adjoint systems we derived with a material derivative free Lagrangian approach, which we rigorously analyzed theoretically in our earlier work [6]. We solved the shape optimization problems numerically using a gradient descent method, where we computed the shape gradient with equations of anisotropic, inhomogeneous, linear elasticity. The numerical results for the shape optimization show that all models behave similarly and that they yield similar optimized geometries.
For future work, the shape optimization of a heat source, e.g., a chemical reactor or an electronic device, coupled to the cooling system can be investigated. This would yield more realistic models and optimization problems. Additionally, the models introduced in this work can be used as building blocks for models describing such a coupled system. Finally, considering a topology optimization before the shape optimization could further increase the quality of the geometries.