An adjoint-based solver with adaptive mesh refinement for efficient design of coupled thermal-fluid systems

A multi-objective continuous adjoint strategy based on the superposition of boundary functions for topology optimization of problems where the heat transfer must be enhanced and the dissipated mechanical power controlled at the same time, has been here implemented in a finite volume (FV), incompressible, steady flow solver supporting a dynamic adaptive mesh refinement (AMR) strategy. The solver models the transition from fluid to solid by a porosity field, that appears in the form of penalization in the momentum equation; the material distribution is optimized by the method of moving asymptotes (MMA). AMR is based on a hierarchical nonconforming h-refinement strategy and is applied together with a flux correction to enforce conservation across topology changes. It is shown that a proper choice of the refinement criterium favors a mesh-independent solution. Finally, a Pareto front built from the components of the objective function is used to find the best combination of the weights in the optimization cycle. Numerical experiments on two- and three-dimensional test cases, including the aero-thermal optimization of a simplified layout of a cooling system, have been used to validate the implemented methodology.

approach is often time-consuming and does not ensure the optimum design. For this reason, optimization techniques such as shape optimization (SO) 1,2 and topology optimization (TO) 3,4 are usually adopted. The first family includes all the methods that involves the deformation and movement of an existing fluid-solid interface while in topology optimization material is added to a domain. The shared aim is to reach an optimal configuration that minimizes an objective function while respecting the problem constraints. Topology optimization works particularly well to produce preliminary results that involve large modifications of the computational domain: it has been first introduced for structural optimization 5 and, from that moment, it evolved incessantly. 6 Different techniques have been adopted to optimize TO problems, such as genetic algorithms 7 and neural networks, 8 but, given the usually high number of control variables, 9 gradient-based approaches are preferred for the majority of the applications. In particular, the gradient of the objective function can be computed by solving the adjoint equations (a.e.) of the problem, making this operation substantially independent of the number of control variables. Two different strategies are commonly employed to solve adjoint equations, 10 namely the continuous and the discrete adjoint methods. In the continuous adjoint method (the one selected in this work), the adjoint equations are derived and then discretized; in the discrete adjoint method, adjoint equations are directly derived from the algebraic equations coming from the discretization of the original problem. The application of adjoint equations to fluid mechanics optimization dates back to References 11 and 12. A TO finite element (FE) framework to optimize the dissipated power in a Stokes flow was presented for the first time in Reference 13; a penalization term representing the inverse permeability field regulates the material distribution within the domain, allowing to recover the Brinkman-type form of the Darcy equation for porous media from the momentum equation. Several works followed, extending the previous approach to laminar incompressible Navier-Stokes equations discretized through the finite volume method (FV). [14][15][16] Adjoint equations formulation for TO was shortly after applied to conjugate heat transfer (CHT) problems. 3,4 In Reference 3, multi-objective optimization was performed to minimize the mean temperature and total fluid power dissipated in a system, using FE for the space discretization and the method of moving asymptotes (MMA) 17 to update the design variables. In Reference 4, a FE solver combined with the solid isotropic material with penalization (SIMP) method to interpolate the material properties from solid to fluid and the MMA was proposed. Turbulence was added to a topology optimization problem applied to an incompressible ducted flow in Reference 16, while constant eddy viscosity (CEV) was used in Reference 18. In both cases, the eddy viscosity was independent on the design variables; this is commonly known as the assumption of "frozen turbulence". The same approach for turbulence was later extended in Reference 19, where the Spalart-Allmaras model was applied to perform RANS simulations of incompressible flows including heat transfer with constant temperature in the solidified regions. The CEV hypothesis was relaxed for shape optimization problems in References 20 and 21. Thermal optimization of complex systems requires the use of a thermal objective function, but in fluid dynamics, this approach often leads to nonrealistic flow patterns involving dead ends and broken tracks. 22 To fix this, researchers proposed the use of multi-objective optimization (MOO) and a twofold objective function including the weighted contribution of thermal and pressure functions. 3,23 The competing optimization of multiple objective functions and constraints is a common problem in aerodynamic shape and topology optimization and can be faced with different approaches. In Reference 24, the constrained multiobjective design problem is posed as a single membership function fuzzy optimization. For each step of the optimization, only the minimum membership function is updated; the multi-objective problem becomes therefore a single-objective optimization which allows a violation of the constraints by a reduction of the cost function. In References 25 and 26, the authors identified a set of optimal points (the Pareto frontier) by an optimization procedure carried out with different combinations of the two objective functions. The application of topology optimization to industrial problems often involves three-dimensional geometries, as apparent by the published literature. In Reference 27, the authors proposed an optimization method for two and three-dimensional steady-state CHT problems formulated as a heat transfer maximization to obtain high-performance cooling devices employing level set boundary expressions and a Tikhonov-based regularization scheme. In Reference 28, TO solver was developed for twoand three-dimensional turbulent complex systems where the adjoint solver is obtained via automatic differentiation. In Reference 29, a Lagrangian optimization approach to minimize pressure losses and maximize the heat transfer in a three dimensions Conjugate Heat Transfer problem with the steepest descent method was used.
The common approach to model the material transition in topology optimization for fluids, also adopted in this work, consists of including the porosity field in the momentum equation through the Darcy term ( )u, so that the velocity is penalized in the solid regions. This approach poses significant challenges from a numerical standpoint: the penalization term is a function of the volume-averaged cell velocity, that varies with the discretization in space. Variations in the mesh resolution and topology correspond to different values of the penalization term in the momentum equations and, in turn, to different flow solutions and rates of convergence. Finally, the penalization term causes discontinuities in the primary flow quantities and, therefore, in the calculation of the sensitivity of the objective function. The coupling of adaptive mesh refinement (AMR) with TO potentially represents a very interesting solution to limit the solution dependency on the mesh of such a family of problems, maintaining a low computational time and good accuracy. The method consists of increasing iteratively the grid density in those regions where it is needed, independently of the mesh resolution of the initial grid. Some examples of coupling between TO and AMR exist in structural optimization, with applications, for example, to structures made of discrete geometric components to ensure well-defined sensitivities in the used geometry projection method. 30 In Reference 31, a similar method to save computational resources in the topology optimization of large-scale stress-constrained problems is presented. AMR is well-established in fluid dynamics and from the first publications on the topics in the early 1980s 32,33 it experienced significant growth and is nowadays widely adopted. An exception to the common use of AMR in CFD is represented by topology optimization, where few studies can be found. In Reference 34 the authors proposed for the first time the AMR to obtain a high resolution of the interface for TO problems, applying the algorithm to 2D laminar cases. In Reference 35, a script for the solution of topology optimization problems for Stokes flow in dynamic unstructured meshes is proposed. The complexity of the test has been later incremented in Reference 36, where a discrete adjoint for the topology optimization of compressible and incompressible flows is discussed. For the sensitivity computation, the authors used a finite difference scheme with a grouping algorithm which, at the price of an approximate and user-dependent calculation, allowed them to face arbitrarily complex problems including AMR. To the authors' knowledge, AMR has not been still applied to continuous topology optimization of thermal problems, where the material distribution and the position of the fluid-solid interfaces continuously change. In this work, AMR is applied within a thermal adjoint solver to ensure proper grid resolution at the fluid-solid interfaces while preserving a coarse grid where possible.

Motivation of this research
The research described in this paper is motivated by the need to design and optimize the thermal management of battery-powered electric propulsion systems for Automotive and Aerospace, a promising technology to limit pollutants emission. Major thermal issues of high-power density machines involve the design of the thermal path in presence of power electronic semiconductors, nonisotropic electromagnetic materials and thermal insulation. In these situations, the main output required by the designer is the optimal positioning of the heat sources 37 and the use of CFD simulation for the optimal design of the cooling system of the electronic boards, as discussed in this work.
In the design of the cooling systems, a limitation of the maximal dissipated mechanical power is applied, while the heat flux between the solid and the fluid must be maximized. The Mach number of the carrier fluid (either liquid or air) is low; also, no phase change usually occurs in the pipes of such systems. We can therefore treat the flow as incompressible. The continuous adjoint approach, which is capable of computing the objective function sensitivities with respect to the design variables, is proved to outperform most of the relevant methods available (direct sensitivity analysis, finite differences, or the complex variable approach) for optimization. For this reason, it has been selected to solve the family of problems discussed in this paper. The effectiveness of the optimization method is evaluated through (a) required performance of the system, in terms of dissipated mechanical power and capability to remove the heat; (b) reliability and speed of calculation of the solution method to reach the optimal solution.

Goals and highlights
The goal of the present work is to develop an efficient continuous thermal adjoint solver for the multi-objective optimization of coupled thermal-fluid problems, where penalization is used to model material transition. The coupling with an adaptive mesh refinement (AMR) strategy is presented and its application to mitigate typical problems of topology optimization, as the high computational expense and the mesh dependency, is discussed. We describe the development of an incompressible, single-phase, multidimensional Finite Volume solver, that has been implemented as an open-source C++ code in the OpenFOAM Technology. Two-and three-dimensional numerical experiments are used for validation. The field of application for which the solver has been conceived characterizes some of its features that are here summarized: -a semi-implicit method for pressure linked equations (SIMPLE) 38 is applied to couple pressure and velocity in the primal and in the adjoint problem; -the objective function is defined by a weighted combination of recoverable thermal and dissipated power through the boundaries; -a density-based approach for the optimization problem is combined with the rational approximation of material properties (RAMP) method 39 to model the porosity and material properties between the fluid and the solid; -sensitivity is computed by solving the adjoint equations and used with the method of moving asymptotes 17 to update the design variables. The MMA method is implemented as a dynamic C++ library that is linked to the solver; -the continuous topology optimization solver has been extended to support the adaptive mesh refinement (AMR); implementation details, including the handling of flux correction to enforce mass conservation across topology changes are discussed, as well as the impact of AMR on the optimal solution of the adjoint method as different refinement criteria are used.
These features are included in the solution algorithm of Figure 2 and will be discussed in detail in the next sections.

Paper structure
The paper is organized as follows: in Section 2, the fundamental theory of the continuous adjoint method for topology optimization of a coupled thermal-fluid system is described; in Section 3, the derivation of the solution strategy for the optimization problem is presented, while in Section 4 the details of the numerical implementation of the solver are discussed. Discretization methods applied to the algebraic operators are described in Section 5. Results are presented in Section 6: the verification of the proposed numerical method without the use of AMR is carried out in a two-dimensional setup, then the use of AMR on three-dimensional cases is compared against the solution on the same geometry without mesh refinement for different weights of the objective function. Finally, a test case based on a real geometry is simulated. Conclusions are drawn in Section 7.

THE CONTINUOUS ADJOINT METHOD FOR TOPOLOGY OPTIMIZATION OF A COUPLED THERMAL-FLUID SYSTEM
The notion of topology optimization was introduced in structural mechanics, where equations for material density were solved to identify areas to add material and increase structural stiffness. The same idea was adapted to fluid dynamics problems later: in fluid problems, the geometry is described by a volume mesh and the porosity is used as a control variable for each cell. The porosity determines which portion of the domain is fluid and which solid, according to the gradient of the cost function. The final porosity distribution determines the optimal topology of the domain. Design modifications are accomplished without changing the mesh with significant savings in the computational effort. The starting point of the continuous adjoint formulation for topology optimization problems consists of computing the derivative of the cost function J with respect to the design variables (known as sensitivity), in a way that is independent of their number. The sensitivity required for a gradient-based optimization is: being b, the vector of the design variables. The optimization procedure employs a gradient-based method to iteratively update the design variables, using the adjoint technique for the computation of the sensitivity. In principle, computing the sensitivities would imply solving the state equations once for each design variable with a significant computational cost. The derivation of continuous adjoint equations through the use of Lagrange multipliers is adopted to make the sensitivity independent of the cardinality of the design variables. 40,41 The main objective of this work is the implementation of a methodology for multi-objective topology optimization of a coupled heat transfer and fluid flow system where the recovered thermal power generated over the boundaries from single or multiple heat sources needs to be controlled. Such problems are subjected to a set of physical constraints R, that can be enforced introducing a Lagrange function L and reformulating the cost function as: where V is the domain and is the vector of Lagrange multipliers, called adjoint variables: being q interpreted as the adjoint pressure, v the adjoint velocity and T a the adjoint temperature, while: is the set of governing equations for the problem, that is completed by their boundary conditions. The governing equations of an incompressible steady flow with heat transfer are written as: -mass conservation: -momentum balance: where ( )u is a Brinkman penalization term, is the porosity, that is updated using the pseudo density as design variable. If tends to zero, the governing equations for the fluid are recovered; for large values of , u tends to zero, and the temperature transport in the control volume is dominated by solid diffusion; -equation of energy, that for incompressible flows is written in the form of temperature transport-diffusion: The solution of the governing/state equations is a constraint for the problem. 16 The resulting optimization problem is then based on the calculation of the sensitivity of the augmented cost function L with respect to the design variables b: 42 The pseudo-density will be used as design variable, b = , as it will be discussed in Section 2.1. In Equation (8), the derivative ∕ b differs from the partial derivative ∕ b because it includes contributions from the domain deformation. If the Leibniz theorem for the differentiation of volume integrals with moving boundaries is applied to Equation (8), it follows: Under the assumption of fixed-boundaries, it follows: Since changes in the design variable force variations of the flow field through the governing equations, the total variation of L also includes contributions from the changes in the primary flow variables u∕ b, p∕ b and T∕ b. Equation (10) can be rewritten in a more convenient form (passages reported in the Appendix) as: so the terms J p , J u , and J T in Equations (13), (14), and (15) vanish. The adjoint equations in the inner domain are therefore independent of the cost function of the problem. This statement is not valid at the boundaries and the use of different cost functions requires different boundary conditions to apply to the same solver.
The solution of the adjoint equation is completed by the set of the boundary conditions: whose general form is: The specific formulation of the adjoint boundary conditions depends on the kind of boundary conditions applied to the flow variables in the primal problem. For the cases considered in this article, these aspects will be discussed in the next sections.

Topology optimization using continuous adjoint
Equation (11) must be implemented and solved into an adjoint CFD solver. The optimization problem consists of minimizing the augmented objective (or cost) function L: where is the nondimensional pseudo density, that is used as design/control variable. In the topology optimization strategy adopted in this work, solid material is added to the fluid domain by a density-based approach. 43 The solution of the design problem should consist in a 0(solid)−1(fluid) field for . If = 1, the governing equations for the fluid are recovered; for = 0, the flow velocity tends to zero and the temperature transport in the control volume is dominated by solid diffusion. On the other hand, this class of problems results unfortunately to be ill-posed 44 and the common approach is to substitute the discrete with a continuous variable, by introducing a penalty that addresses the solution to a binary one. Pseudo-density also appears in the definition of the material properties (density , specific heat Cp, thermal flow diffusivity D, and porosity ), that are calculated by the rational approximation of material properties (RAMP) model: 39 In Equations (24)-(27), the subscript f is used for quantities related to the fluid and s for the solid; the parameter q governing the shape of the functions is set to q = 0.1. 13 In thermal cooling of power electronics of propulsion systems, the maximal total pressure drop in the cooling circuit must be limited while the heat is removed. For the considered problem, the cost function J of Equation (2) is defined as: where J 1 is the mechanical power dissipated by the fluid through the boundaries, while J 2 is the net thermal power recoverable from the domain: 25,26 and In Equation (28), a maximization is desired for the recoverable thermal power only, while J 1 must be minimized. Functions J 1 and J 2 are defined at the boundary S, while their value is null in the inner domain. While J 1 is related to the features of the flow field (Equation 29), J 2 is linked to temperature transport and diffusion (Equation 30): their magnitude can therefore be very different, so it is convenient to normalize them. A possible normalization approach is to compute the optimal values of J 1 and J 2 respectively from the optimization of the fluid problem and the thermal problem, and use these values to scale the corresponding objective function; 26,45 the weighted sum method is finally applied to perform the multiobjective optimization as: being J 1 and J 2 the normalized functions and w ∈ (0, 1) the weighting factor. Once the cost function is defined, a constrained optimization problem, with the constraints being the governing/state equations, is obtained. The optimization problem must respect the following constraints: -maximum quantity of material to be added to the domain, set to avoid trivial full solid solutions: where n c is the total number of control volumes (CVs) in the domain and Φ max the maximum volume fraction that can be solid; -the governing equations for the flow problem (see Equation 4), enforced through Lagrange multipliers method (Section 2).
If the adjoint variables are chosen in a way the adjoint equations and their boundary conditions are equal to zero: the computation of the gradient becomes independent of u∕ b, p∕ b and the sensitivity (Equation 11) reduces to the following algebraic relation: The computational cost for the calculation of the adjoint problem is similar to cost of the solution of the flow equations (primal problem).

Boundary conditions for the adjoint problem
The boundary conditions for the adjoint variables are derived from Equations (20)-(22) when the relation BC a = 0 is enforced. Depending on the type of boundary condition imposed to the primal flow variables, different boundary conditions for the adjoint variables are obtained. 16 -Inlet: velocity u and temperature T are set, so their derivative with respect to is 0. Equations (20)- (22) reduces to: The relations for the boundary conditions can be further manipulated to obtain: -Outlet: at the outlet, zero Neumann conditions are imposed for the velocity and the temperature, while the value of the pressure is fixed. As a consequence, p∕ = 0 and Equations (21)-(22) thus become: From which the boundary conditions for q, for the tangential component of the adjoint velocity v t , and for T a can directly be extracted: The boundary conditions for the normal component of the adjoint velocity v n are obtained from continuity: being ∇ || the in-plane component of the derivatives at the boundary. -Adiabatic wall: the boundary conditions for the pressure and velocity at the wall are as at the inlet boundary (Neumann for p, Dirichlet for u). As a consequence, the boundary conditions for adjoint pressure and velocity are the one expressed in Equations (38)- (40). The normal derivative of the temperature is fixed and u = 0, so the boundary condition for the adjoint temperature reads: -Fixed-temperature wall: if the temperature is fixed at the wall, the derivation of the boundary condition for T a leads to: Boundary conditions for the adjoint pressure and velocity are as in Equations (38)-(40).

DESIGN PERFORMANCE
The design performance of topology optimization by continuous adjoint varies due to different parameters, such as material properties. Also, the use of the penalization term introduces a well-known mesh dependence of the solution: in structural mechanics, it has been related to the nonexistence or nonuniqueness of an analytical solution for the constrained optimization problem 46 and tackled through the application of an implicit filtering scheme for the sensitivity of the objective function. [47][48][49] Alternatively, a local adaptive mesh refinement (AMR) can be used during the simulation.

Implicit filtering scheme for the sensitivity of the objective function
An implicit filtering scheme for the sensitivity of the objective function in the form of a Helmohltz-type PDE 47 has been implemented. Filtered sensitivityf is computed as: where homogeneous Neumann conditions are set at the boundaries: In Equation (50), r is the characteristic length of the isotropic filter,f and f are respectively the filtered and unfiltered sensitivities. Filtered sensitivity is eventually used to update the pseudo density .
F I G U R E 1 Representation of a mesh section for four consecutive optimization steps with adaptive mesh refinement (AMR). The color denotes the number of refinement undergone by each cell: , no AMR; , 1 refinement; , 2 refinements. In iter = 3, AMR allows to ensure a high grid resolution at the fluid-solid interface. when the material distribution is modified during the optimization process. [Colour figure can be viewed at wileyonlinelibrary.com]

Adaptive mesh refinement for topology optimization
Error reduction in CFD simulations could in principle be obtained by increasing the mesh resolution in the regions where the large gradients are expected. In topology optimization, this is not possible: the material distribution is not known a priori and it is updated at run-time, at each step of the adjoint optimization process. With adaptive mesh refinement (AMR), the grid density can be adapted to the material distribution to adequately capture discontinuities in the flow field, where usually steep gradients of the primal and/or of the adjoint variables are found. The method presented in this work is designed to (a) dynamically adapt the mesh resolution in the transition regions from fluid to solid, maintaining an overall coarse grid and a moderate computational effort for the simulation; (b) make the solution of the adjoint solver independent of the initial grid topology. Starting from an initial grid  0 , the applied AMR performs a hierarchical h-refinement strategy 50 that takes advantage of the polyhedral mesh support of the solver. Points and faces are dynamically added locally to the control volumes, without producing modification to the rest of the grid. Cells involved in topological modifications are selected based on an indicator i( ): the mesh is not refined with i min and i max respectively the maximum and minimum threshold for the refinement region. The choice of the refinement indicator is done before the execution of the simulation and depends on the flow properties that should be enhanced: if the goal is to adequately capture the fluid-solid interface, the mesh should be refined according to the quantities that undergo a sharp variation in this region, that is, porosity and velocity. AMR is performed before the solution of the direct and adjoint problems and of the update of the flow properties. In Figure 1, it is shown how AMR is applied near a varying solid/fluid interface.

SOLVER IMPLEMENTATION
The solution method is shown in Figure 2. The sequential solution of the governing equations for the primal and adjoint problems is carried out on a polyhedral mesh that is dynamically updated during the optimization cycle, then the sensitivity is computed and used to update the design variables by a gradient-based method: 3,13,25 this can be for example the steepest descent method 15 or the method of moving asymptotes. 17 In this work, MMA has been considered for the calculation of the design variable; the algorithm has been implemented in the form of a C++ class to be linked as a dynamic library to the adjoint solver. Reynolds-averaged-Navier-Stokes (RANS) modeling for turbulence has been used for closure; the frozen turbulence hypothesis is applied to simplify the processing of turbulence for the adjoint problem. Either in the primary and the adjoint problem, pressure-velocity coupling is treated by combining mass and momentum equations in a sequential solution, implemented in the form of the SIMPLE (semi-implicit method for pressure linked F I G U R E 2 Solution method for the thermal adjoint solver: once the primal and the adjoint problems are solved, the sensitivity is computed; finally, a gradient-based method is applied to update the design variables before mesh adaptivity is eventually employed.

Solution of the primal flow problem
In the primal problem, the segregated implicit solution of the governing equations (Equations 5,6,7) for a homogeneous, incompressible, steady Newtonian fluid with heat transfer is performed. The conservation of mass, Equation (5) with matrix A including the diagonal components only; the intermediate velocity can be computed treating explicitly the extra-diagonal terms as: The velocity field u * does not satisfy the continuity equation. To enforce mass conservation, a correction equation (written in semi discrete form) is created by substituting velocity fluxes from Equation (54) into the mass conservation (written, applying the divergence theorem, as a sum fluxes evaluated at the cell faces ∇ ⋅ u = ∑ f f = 0) and considering a pressure-correction term: The corrected face fluxed are calculated as follows: Equation (56) is also known as flux corrector equation. Finally, under-relaxed cell-centered pressure is used to correct the intermediate velocity: In Equation (57), p is the under-relaxation factor on pressure. The corrected velocity is treated as a new guess at the next iteration loop, until a converged solution is obtained (see Figure 2). Before passing to the new outer iteration, the transport equation for the flow temperature at the cell center is solved:

Solution of the adjoint problem
The structure of the solution method for the primal and the adjoint problem are similar, since the adjoint equations resembles the primal counterparts ( Figure 3B). In the solution of the adjoint problem, a one-way coupling between the temperature and the pressure-velocity equations applies: the adjoint temperature T a can be determined independently of the adjoint velocity v and pressure q. Equation (15) is therefore solved before the coupled equation for adjoint pressure and velocity: Once the adjoint temperature T a is known, the coupling of Equations (13)- (14) is calculated in a segregated manner. First, an intermediate adjoint velocity v * is calculated in a predictor step, where the adjoint pressure gradient is neglected: Differently from its primal counterparts, the adjoint momentum balance is linear. At the same time, numerical difficulties arises for multiple reasons: the term ∇uv, commonly called Adjoint Transpose Convection (ATC), is an instability source for the solution of the adjoint equations 51 even at moderate Reynolds numbers. A definitive solution to numerical instabilities due to ATC is not available in the literature. In the present work, the ATC is treated explicitly and its value is bounded, if needed. Also the T∇T a term appears as a forcing term in Equation (14). In regions with large adjoint temperature gradients, this term can be large and may slow down the solver convergence.
Once the intermediate value of adjoint velocity v * from Equation (60) is known, it is used to calculate the adjoint pressure: and q ′ is used to correct the adjoint velocity fluxes: The cell-centered adjoint pressure is under-relaxed: and the intermediate adjoint velocity is corrected to provide the updated value v n+1 : The solution method is summarized in Figure 3B. Similarly to what is done for the final flow problem, the corrected adjoint velocity is treated as a new guessed value at the next iteration loop, until a converged solution for the adjoint problem is obtained.

Flux correction with adaptive mesh refinement
A key aspect for a formally correct AMR consists to ensure the fulfillment of the continuity, Equation (5), across the cells where AMR is triggered, which is obtained by a correction applied to the face fluxes. As for any topological mesh change, AMR requires mapping the resolved quantities (physical and material properties, design variables) onto the newly added volumes. 52 The discretization of the operators is thought for a collocated grid arrangement of the variables, which are evaluated at the cell centers, while face values interpolated. As a mesh refinement is triggered across two consecutive iterations n and n + 1, flow quantities are remapped from the old to the new grid. In the Finite Volume method, spatial interpolation of flow quantities does not ensure continuity: remapping applied to surface quantities does not account that interpolation weights are dynamically changing with the grid, so it is not consistent with the remapped cell-centered quantities. For this reason, face fluxes must in principle be corrected as the mesh is changing. In the context of a segregated solver, this can be done by a pressure correction applied to the face fluxes mapped onto the new grid. Before the solution of both the primal and adjoint flow fields, mass conservation is enforced to correct face fluxes across a topology change, by solving a Poisson equation as: where p c is the corrected pressure, u n (x n+1 ) denotes the velocity field computed at the nth outer iteration of the segregated solver, but remapped onto the new mesh (x n+1 ). Equation (65) is completed by the boundary conditions: p c = 0 on permeable boundaries.
The corrected face-interpolated mass fluxes c f are computed as: and continuity is recovered. It can be noticed that the solution strategy for the adjoint equation is similar to the one of the direct equation and that an adjoint continuity equation must be satisfied as well, so the idea of the flux correction should be extended to the adjoint variables. Figure 4 shows the effect of the flux correction: both continuity errors and adjoint continuity errors, computed before the SIMPLE outer iterations, are displayed for a test case in which AMR is performed every 100 iterations. Flux correction drastically limits the peak in the continuity error, that would occur anytime a topological modification is triggered. Despite flux correction does not have a significant impact on the overall convergence rate of the solver, it favors an improved calculation of the initial solution and the convergence of the linear solver.

VARIABLES POSITIONING AND SPATIAL DISCRETIZATION
For the spatial discretization of the developed solver, finite volume schemes based on collocated grid arrangement are used. Primary and adjoint variables are defined at the cell centers and their derivatives are computed as in the following: -Diffusive term (Laplacian) of a quantity Ψ, for example, ∇ ⋅ (Γ∇Ψ): where ∇ n Ψ f is the surface normal gradient of Ψ. The subscript f in Equation (69) indicates the cell-to-face interpolated quantities. In the present work, linear cell-to-face interpolation has been applied: for irregular polyhedral meshes, interpolation is generalized by defining a weight w for each face: where Ψ f is the face-interpolated quantity. Subscripts P and N indicate values at the centers of two neighboring cells.
For the nonorthogonal grid of Figure 5 in a collocated variable arrangement, the surface gradient of a quantity Ψ is decomposed onto an orthogonal part and a (nonorthogonal) correction (see Figure 6): where = 1 n f ⋅d and ∇ n Ψ f is the uncorrected normal gradient from the two values of the two cells sharing the face. The explicit part is computed from Equation (70) as: being w = 0.5 in this work; the normal gradient is computed as: -Gradient terms: we use the Green-Gauss theorem: being V P the volume of the polyhedral cell P, and S f the surface vector of the f th face of the cell.  53 is used to mimic the staggered-grid discretization to prevent checkerboard effects. Using the divergence theorem, the convective terms are rewritten as: The velocity u f is interpolated with the same approach presented in Equation (70), while a second-order central differencing scheme is used for the fluxes.

RESULTS
Simulations have been performed on three different cases to test the numerical properties of the newly implemented solver in terms of ability to: (a) provide results in accordance with the existing literature; 25,26 (b) provide distinct nontrivial optimal solution when different input parameters (e.g., weighting factors and solid material) are imposed; (c) provide results that are independent on the initial grid resolution, thanks to the use of the AMR. Similarly to what done in previous work, 25,26 the simulation of a two-dimensional square channel with single inlet and single outlet subjected to constant wall temperature on top and bottom walls is used to assess the properties of the optimization strategy. The same square channel geometry is then extended to three-dimensions and tested to show the convergence to a unique final configuration when the adaptive mesh refinement is applied. Finally, the robustness of the solver and its application is tested on a geometry of industrial interest.

Two-dimensional square channel
The correct operation of the implemented solver is verified on a simple configuration: a two-dimensional square channel whose side is A= 0.1 m, with single inlet and single outlet and subjected to constant wall temperature on top and bottom walls. The full detail of the geometry studied and the material properties are summarized in Figure 7. The velocity profile at the inlet (whose size is A∕5) is parabolic with a centerline (peak) velocity of the u max = 2.25 ⋅ 10 −2 m s −1 corresponding to Re=30. Being the inverse permeability max = 150 s −1 , the nondimensional Darcy Number Da = ∕ max A 2 = 1 ⋅ 10 −5 is sufficiently low to ensure solid impermeability. 54 Two sets of tests were selected in this work to verify the dependency of the results on: (a) the value of the weighting factor w of Equation (31) within an interval w ∈ [0.1; 0.8] for solid 1 (left column, Figure 8A) and solid 2 (right column, Figure 8B); (b) the value of the thermal diffusivity of the solid material. A 100 × 100 cells grid has been used. For w ≤ 0.6, optimization of the heat recovery is dominant: as apparent from Figure 8A, for solid 1 the porous material collects at the center of the channel and forces the flow to split and skim the hot walls; the fluid velocity is maximum in the region near the walls and a large amount of thermal power is recovered. The material is distributed over striped structures, limiting in that way the heat transferred to the highly diffusive solid. For w > 0.6, the solid material at the center of the channel is progressively removed to minimize the pressure drop across the channel. The recovered thermal power is transferred by conduction in the solid region, that extends from the middle of the channel to the heated walls. A similar analysis is presented with solid 2 of Figure 7, that is characterized by a lower thermal diffusivity D. For w ≤ 0.6, the solid stripes are now replaced by a single region of insulating material, that becomes larger for higher values of the weight w. The single-block still appears in the solution for w ≤ 0.65 but, for w = 0.8, the optimal shape of the solid region drastically changes: the solid material concentrates mostly near the inlet section. It is important to note that the threshold value of w at which the solid material is arranged far from the channel center is different from solid 1 and solid 2, because the thermal diffusivity has a direct influence on the flow temperature and, in turn, on J 2 . Hence, different thermal diffusivities of the solid (D solid1 > D solid2 ) used for the two sets of simulations have a severe impact on the final optimal solution, as apparent in Figure 8A and 8B. A Pareto front (see Figure 9) is built from J 1 and J 2 for different w and provides at the same time the value of both the objective functions; in optimization, this tool supports the designer towards the desired solution. In Figure 9A, for w > 0.65, the Pareto efficiency reaches its maximum for any variation of J 2 ; similarly, for w < 0.04 the recoverable thermal power is not changing significantly with the fluid dissipated power. For 0.04 ≤ w ≤ 0.65, J 1 and J 2 of the objective function contributes to determine the global Pareto efficiency. The optimal solution of the constrained problem is the trade-off solution, that is dependent on the desired enhancement of J 1 and J 2 . Similar considerations can be made for solid 2 when Figure 9B is examined. It is worth noting that the sudden variation of material arrangement observed for solid 2 at w > 0.65 is visible in the Pareto front and corresponds to a significant reduction of the recoverable thermal power. Results of Figure 9 are qualitatively in agreement with, 26 as well as the flow characteristics of the optimal solution. In particular, the Pareto front assumes a convex shape, proving that the obtained solutions are indeed optimal. On the other hand, the optimized shapes of Figure 8 differ from the published solutions; this is not surprising due to the differences of the material properties and optimization parameters. This confirms as the optimization process is strongly influenced by the choice of the input parameters and by the physics of the problem, and possibly affected by local minima problem. 46

Three-dimensional square channel
In this section, simulations on the square channel have been extended to three-dimensions: differently from the two-dimensional case, a uniform velocity profile is imposed at the inlet, (u = 2.25 ⋅ 10 −2 ms −1 ). A fixed temperature is set at the upper and the lower walls, while side walls are adiabatic. The material properties correspond to fluid and solid 1 of Figure 7. The aim of this section is to prove the potential of the adaptive refinement technique applied to the adjoint solver; this is achieved by studying the influence of relevant input parameters on the optimization results. As summarized in Table 1, the weighting factor w, the refinement criterion and eventually the mesh type (static grids vs dynamic/AMR) vary in the different simulations; the results are reported in Figure 11. Results on simulations on fixed size (static) grids

Refinement criterion
are shown in Figure 11A-D: two grids of 20 3 and 80 3 elements are reported for two different weights, w = 0.5 and w = 0.8; the solution of the 80 3 cell mesh will be assumed as reference for comparison with dynamic simulations. For w = 0.5 (i.e., higher relevance of J 2 than J 1 ), the solid material is mostly located at the cavity center and the fluid is forced to change its path and deviate towards the hot walls. As the minimization of the pressure losses becomes more important, the solid material concentrates near the walls and the fluid flows in the central region of the channel; as a consequence, a minor amount of dissipated and recoverable thermal power is observed. System performances are summarized in Figure 10 and are measured in terms of percent variation of the optimized objective functions J OPT with respect to the reference (nonoptimized) solution J REF , that is computed on the 80 3 elements grid: The addition of solid material in the proximity of the channel center for w = 0.5 results in an increase of both the dissipated mechanical power (≃ 16%) and recovered thermal power (≃ 77%), as represented in Figure 10A. For w = 0.8 ( Figure 10B) a lower J 1 (≃ −7.5%) and a limited thermal recovery (ΔJ 2 ≃ 57%) is observed. The overall predicted dissipated power and thermal recovery do not seem very sensitive to changes in the mesh resolution; this does not hold for the distribution of the solid material, as apparent comparing Figure 11A,C. Differences are also visible in Figure 11B,D, proving that the optimal distribution of the solid is sensitive to the mesh resolution. The application of Adaptive Mesh Refinement can potentially make the prediction of the optimal material distribution independent on the initial mesh resolution. As summarized in Table 1, simulations with AMR are performed on two grids of different initial resolution (20 3 and 40 3 elements). Two refinement criteria for AMR have been tested: (a) the normalized gradient of the porosity ∇ N , that should follow the discontinuity represented by the fluid-solid interface, and (b) the magnitude of the normalized gradient of the velocity ∇u N , used to refine the fluid domain. In Figure 11E-H the optimization results for the 20 3 ( Figure 11E,F) and 40 3 ( Figure 11G,H) grids with i = ∇ N are reported. Cells are refined at the beginning of the optimization cycle if ∇ N falls in the interval 0.001 ≤ i ≤ 1. The use of AMR drastically contributes to limit the dependence of the solution on the size of the initial mesh: results from Figure 11E-H closely resemble the results obtained with the static fine grid. At the same time, AMR allows to increase the mesh density in specific regions, while limiting the overall number of elements; this is confirmed by Figure 12, where the number of computational points and the value of the objective function computed for the two AMR cases are compared with the reference cases of Figure 11C,D. For w = 0.5, the cell number (and thus the computational time) is reduced, with respect to the 80 3 grid, of the 68% for case Figure 11E and 64% for case Figure 11G; the objective function is comparable with the reference case Figure 11C. A similar trend is observed for w = 0.8: the reduction of the cells is of 84% for case Figure 11G and 75% for case Figure 11H. Again, the objective function has the same value for cases Figure 11H and case Figure 11D, while an error is introduced for case Figure 11F; this is due to a different prediction of J 1 , as apparent from Figure 10B. If ∇u N is selected as refinement criterion, the mesh topology can be updated multiple times at the same step of the optimization cycle. Additional cases ( Figure 11I-L) have been produced using ∇u N and refining four times the mesh for each optimization. Being u ≃ 0 in the solid region, the mesh is refined in the fluid domain. Solutions obtained with AMR are similar to the reference solution for w = 0.8, while some differences in the distribution of the solid at the center of the channel appear when w = 0.5. The differences in the distribution of porosity is not visible in the plot of the objective function ( Figure 12). The number of cells at the last iteration is lower than in cases Figure 11C,D, by 49% ( Figure 11I) and 45% ( Figure 11K) when w = 0.5, 74% ( Figure 11J) and 66% ( Figure 11L) when w = 0.8. Figure 13 shows the speedup for the cases discussed in Table 1, assuming case c and d as reference for w = 0.5 and w = 0.8 respectively. The computational speedup is: The reduction of the overall mesh size with AMR is only one of the parameters affecting the simulation performance. The speedup obtained with AMR also depends on the selection of the weights w of the objective function and on refinement criteria adopted; the simulation setup influences the location and the size of the refinement regions, the number of refinements that are triggered at run-time and the load unbalance among processors. Moreover, each time AMR is triggered, computational resources are spent in the re-meshing. The combination of all these effects contributes to the overall speedup reported in Figure 13, which shows as a reduction of the computational time is always achieved with AMR.

Thermal management of an electronic board
The thermal adjoint solver has been finally applied to optimize a representative geometry for the cooling system of an electronic board ( Figure 14). The studied domain consists of a thin rectangular cavity with three heat sources A, B, and C at T = T max on the upper wall. A uniform flow at Re = 5000 is set at the inlet boundary Γ in . The − model is applied to model turbulence; the frozen turbulence hypothesis is adopted. The turbulent thermal diffusivity is computed as D T = D + T ∕Pr T , being Pr T = 0.75 the turbulent Prandtl number. AMR on the initial mesh of 4.4 ⋅ 10 5 elements is triggered by monitoring the quantity ∇ N so that, in light of the discussion of Section 6.2, the fluid-solid interface is captured with a fine grid. Two different combinations of the objective functions are studied, corresponding to two different sets of simulations, in the following case a and case b (see Table 2). The minimization of the dissipated power is the main goal in case a, while the maximization of the recovered thermal power is the main goal of case b; w 1 and w 2 of Table 2 are determined with a trial and error procedure and set to values that guarantee a clear contribution of J 2 on the optimal geometry. The results of the two optimizations are reported in Figure 15, where the final material distribution is displayed together with

F I G U R E 11
Optimization solution for three-dimensional grids of different resolution with and without Adaptive Mesh Refinement (cases summarized in Table 1). For each subfigure it is shown: optimal material distribution calculated by the solver (left) and mesh refinement level together with the velocity streamlines, colored accordingly to the normalized temperature fieldT (right). [Colour figure can be viewed at wileyonlinelibrary.com] the computational grid at the channel half-height. The different weights for J 1 and J 2 in the objective function clearly influence the material distribution. The final solution for both the cases is nontrivial and difficult to predict a priori; AMR may therefore represent an effective strategy to adapt the mesh to the material configuration and to concentrate the computational efforts at the fluid-solid interfaces. The optimized material configuration in case a ( Figure 15A) forces the majority of the fluid to turn towards the outlet direction, while a less significant part of the remaining flow is redirected to sources B and C. This solution favors an increment of the recovered thermal power since part of the flow is forced to skim the heat sources, with limited mechanical power losses. In case b (Figure 15B), the flow splits in two directions: part of it reaches directly the outlet, while a part is directed towards the heat sources near the end of the channel. Interestingly, next to the heat sources, the material is arranged in pin-and fin-like structures, which are commonly employed solutions to increment the transferred thermal power in systems as heat exchangers. In that way, the heat transferred from the fluid to the solid is maximized. With this second solution, the extracted mechanical power is incremented but, given the complexity of the layout, the dissipated power is larger than in the previous solution. Results of Figure 15 are consistent with data of Table 2: ΔJ 1 and ΔJ 2 represent the variation of the optimized objective functions J OPT , computed as described in Equation (77). The final solution is related to the choice of the weights w 1 and w 2 : for case a, both J 1 and J 2 are optimized, resulting in a 10% decrease of the dissipated power and a 162% increase of the recovered thermal power with respect to the initial solution. In case b, the increase of J 2 is of the 660%, but this is paid by a small reduction of J 1 : 3.5%. At the end, the best solution is determined by the design requirements (usually the mechanical power available from the pumping system). TA B L E 2 Parameters and results for the optimization of the thermal management system for an electronic board.

F I G U R E 15
Optimal solid material distribution for the two test cases: in Case A the maximum optimization weight is given to J 1 , in Case B the maximum optimization weight is given to J 2 (see Table 2). [Colour figure can be viewed at wileyonlinelibrary.com]

SUMMARY AND CONCLUSIONS
We implemented a continuous thermal adjoint solver within the finite volume method for the optimization of coupled thermal-fluid systems, supporting adaptive mesh refinement (AMR) on polyhedral meshes. Algorithms are developed in the form of an object-oriented C++ code and in a set of dynamically linked libraries, including the implementation of run-time selectable gradient-based methods to update the design variables and models to handle properties interpolation for a porous material. It is shown that the Pareto front built from the optimal values of the objective functions tends to a convex shape; this is used to prove that the solver behaves correctly and to find the best combination of J 1 and J 2 in the optimization cycle. Numerical experiments highlight how different solid distributions may correspond to similar values of the objective function, proving that the solution of a thermal adjoint solver is mesh dependent. The use of AMR has been proved to reduce the mesh dependency: in particular, if a refinement criterium based on the porosity gradient is applied, a unique mesh-independent solution is achieved. A flux correction to enforce mass conservation across the mesh topology changes from AMR has been proposed. Despite in the case of a steady-state solution the initial error caused by mesh remapping cancels out as long as the convergence is reached, flux correction can speed up the convergence, thanks to a more accurate initial solution after the mesh update. A limitation of an AMR-based approach is mostly related to the work per processor and to the locality of the objects that change during computation. For optimal performance, a processors' work redistribution would be suggested; this is done by dynamic re-partitioning, 55 which is natively supported by the method. The ingredients of the proposed solution strategy, taken individually, are not particularly complex. On the other hand, this novel approach includes notable features that contribute to increasing the fidelity of numerical results of the thermal adjoint and, most importantly, that reduce the solution dependence from pre-processing. For this reason, the proposed method looks suitable for applied research and industry.

ACKNOWLEDGMENT
Authors gratefully acknowledge the Laboratory Computing Resource Center (LCRC) of the Argonne National Laboratory (Lemont, US) for the computing resources provided and Vitesco Technologies SAS for providing data and geometries for validation and testing. Algorithmic developments have been implemented and tested by the authors at the Dept. of Aerospace Science and Technology (DAER), Politecnico di Milano and are based on the most recent OpenFOAM Technology available at the time the draft was submitted. Open Access Funding provided by Politecnico di Milano within the CRUI-CARE Agreement.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

APPENDIX . DERIVATION OF THE ADJOINT EQUATIONS FOR THE COUPLED FLUID-THERMAL PROBLEM
Considering that the spatial and partial differentiation ∕ b commute, 56 the terms of Equation (10) can be developed to obtain the adjoint equations (for brevity, integration variables are omitted inside the integral): a Mass conservation: b Momentum balance: considering the convective term )) : that can be rewritten, considering the continuity equation, as: With similar steps, the following relation valid for the second of the convective terms can be obtained, being (∇v) ij = v j The viscous term v ⋅ ( ∇ ⋅ (2 ( u b )) ) (see Equation A2) can be rewritten using index notation as ) v i and, fulfilling a first integration by parts and applying again the divergence theorem, it becomes: Performing a second integration by parts to the last term of Equation (A6), the following result is obtained: ) . (A7) The diffusive term can then be rewritten as: The pressure term is developed as: c Temperature transport diffusion: considering separately the two terms of Equation (A10): that respectively produce: because u i ∕ b x i = 0 due to continuity equation. With similar steps: For what concerns the second term of Equation (A10):