The Cut-Cell Method for the Conjugate Heat Transfer Topology Optimization of Turbulent Flows Using the “ Think Discrete–Do Continuous ” Adjoint

: This paper presents a topology optimization (TopO) method for conjugate heat transfer (CHT), with turbulent flows. Topological changes are controlled by an artificial material distribution field (design variables), defined at the cells of a background grid and used to distinguish a fluid from a solid material. To effectively solve the CHT problem, it is crucial to impose exact boundary conditions at the computed fluid–solid interface (FSI); this is the purpose of introducing the cut-cell method. On the grid, including also cut cells, the incompressible Navier–Stokes equations, coupled with the Spalart–Allmaras turbulence model with wall functions, and the temperature equation are solved. The continuous adjoint method computes the derivatives of the objective function(s) and constraints with respect to the material distribution field, starting from the computation of derivatives with respect to the positions of nodes on the FSI and then applying the chain rule of differentiation. In this work, the continuous adjoint PDEs are discretized using schemes that are consistent with the primal discretization, and this will be referred to as the “ Think Discrete–Do Continuous ” (TDDC) adjoint. The accuracy of the gradient computed by the TDDC adjoint is verified and the proposed method is assessed in the optimization of two 2D cases, both in turbulent flow conditions. The performance of the TopO designs is investigated in terms of the number of required refinement steps per optimization cycle, the Reynolds number of the flow, and the maximum allowed power dissipation. To illustrate the benefits of the proposed method, the first case is also optimized using a density-based TopO that imposes Brinkman penalization terms in solid areas, and comparisons are made.

Irrespective of the approach used to compute the optimal material distribution, the latter is usually introduced into the flow equations through the so-called Brinkman penalization terms, which attempt to drive the flow solution to zero in the solidified parts of the domain.This leads to the imposition of the wall boundary conditions only in a weak sense, which, among other outcomes, may allow some flow leakage into the solidified domain, making TopO prone to inaccuracies and potentially leading to sub-optimal solutions.Converting the TopO solution to a body-fitted grid and successively performing shape optimization can mitigate this effect [27], without being able, however, to make topological changes in this second stage.Attempts to impose boundary conditions in a stronger sense in TopO, using the finite element/volume methods (FEM/FVM), have been made.Regarding the former, the XFEM approach was introduced in [28] to perform TopO based on the level-set approach for 2D laminar flow problems.Recently, it has been utilized to increase the accuracy of TopO solvers in laminar [29] and turbulent flows [20], including heat transfer.In FVM methods for TopO, [30] utilized a ghost-cell immersed boundary method to apply wall boundary conditions for turbulent flow problems, without having to modify the background mesh.In [31], dealing with laminar flow TopO problems, a cut-cell approach was used to introduce new boundary faces at the FSI, allowing the solution of the flow equations on a body-fitted grid and the imposition of boundary conditions in the strong sense.The present paper utilizes a similar methodology and extends it to tackle TopO problems under turbulent flow conditions, including CHT.The flow solver utilizes also adaptive mesh refinement (AMR) to better resolve the flow and thermal boundary layers around the FSI.AMR has been used recently to enhance the accuracy of TopO using the density [21] and the level-set approach [22], but not in the context of TopO with immersed boundary methods, as in this work.The imposition of exact boundary conditions and the AMR affect not only the accuracy of the flow solver but also the optimized geometries designed by it; see Section 4.
The majority of works on TopO rely on gradient-based optimization methods, supported by continuous [10,32] or discrete [11,33] adjoints, to compute the sensitivity derivatives of the objective and constraint functions with respect to the design variables.In a continuous adjoint, the adjoint PDEs must be discretized and numerically solved; devising appropriate discretization schemes is a challenge affecting the accuracy of the computed sensitivity derivatives [34].In a continuous adjoint, it is important to mention the ease of implementation, the physical insight into the adjoint terms/equations, and the lower computational cost and memory footprint of the adjoint code.On the other hand, the discrete adjoint equations are always consistent with the discretized flow equations [33] and may also retain the convergence characteristics of the flow solver [35]; however, the discrete adjoint is known for its larger memory requirements and/or computational costs [36,37].Recently, the present group of authors developed an adjoint methodology named the "Think Discrete-Do Continuous" (TDDC) adjoint, leading to discretization schemes for continuous adjoint equations that are consistent with those of the primal solver, drawing inspiration from the discrete adjoint.Practically, the new TDDC adjoint bridges the gap between the two adjoint approaches, combining the best of both worlds.Thus, using the TDDC adjoint, one enjoys the low cost and memory requirements of the continuous adjoint and the consistency of the discrete adjoint.The TDDC approach was briefly introduced in [38] in shape optimization and is further elaborated upon herein, constituting the second main focus of this article.
This article is structured as follows.Section 2 describes the developed cut-cell approach for CHT TopO problems, including the utilized AMR approach and the primal and TDDC adjoint equations.Section 3 presents the methodological differences between the developed cut-cell TopO approach and a density-based one, with Section 4 showcasing these differences in CHT TopO examples.Here, the importance of applying proper boundary conditions and using AMR in TopO becomes clear.Finally, Section 5 summarizes the findings of this article.
The methods described in this article are implemented in the in-house variant of the adjointOptimisation library of OpenFOAM, developed and made public by the present group of authors.Using the proposed method, the objective function values recorded during (and, of course, at the end) of the TopO loop are accurate, thus eliminating the need to re-evaluate the optimized geometry on a body-fitted grid.In addition, the TDDC adjoint ensures that consistent gradients are computed throughout the optimization.

The Cut-Cell TopO Method
The solution of the TopO problem starts with the generation of a background grid for the analysis domain, Ω d .In most cases, a typical Cartesian grid suffices.An artificial "density" field, α, stored at the cell centers of the background grid, corresponds to the design variables to be updated in each optimization cycle.The α field is transformed into an auxiliary material distribution field, β (see Section 2.1), indicating the spatial distribution of a fluid or solid.Areas with β ≈ 1 correspond to solids and those with β ≈ 0 to fluids.Cells with intermediate values appear close to the FSI.These cells are refined (h-refinement; see Section 2.2) to effectively model the flow phenomena close to the wall.Then, the β w = 0.5 iso-surface, acting as the interface of a fluid and solid, is computed.Its intersections with the background grid cells give rise to cut cells separating the initial grid into fluid, Ω F , and solid, Ω S , domains.The governing flow equations, without the Brinkman terms usually met in TopO, are solved only in Ω F by applying proper boundary conditions along its wall boundaries, as if a body-fitted grid was used.The solution of the temperature equation takes place over both Ω F and Ω S .The continuous adjoint problem is, then, formulated for the objective function J to be minimized and each flow-related constraint g i to compute their derivatives with respect to the coordinates of grid nodes that lie at the FSI.By applying the chain rule of differentiation, these are transformed into derivatives with respect to the design variables α to be used to update the latter using the Method of Moving Asymptotes (MMA) [39].The MMA creates and solves a convex optimization sub-problem in each TopO cycle, by creating convex approximations to the objective and constraint functions based on their computed gradients.The various steps of the TopO algorithm are presented schematically in Figure 1 and are further analyzed below.

Filtering of the Design Variables
To avoid grid-dependent solutions and eliminate the well-known checkerboard pathologies [40], a filtering scheme for the α field is utilized.This computes an auxiliary smoothed density field, β, by solving a Helmoltz-type filter PDE [41], where R F (β, α) denotes the filtering operation and r is the radius, usually defined as 5 to 10 times the average edge length in the background grid.In addition, the usage of filtering "broadens" the effect of the design variables on the computed FSI and accelerates the convergence of the TopO.

Grid Refinement
In each optimization cycle, the grid close to the FSI is adequately refined in order to accurately simulate the flow and thermal boundary layers.Cells in the vicinity of the FSI are refined using an h-adaptivity approach.In what follows, the cells of the coarse grid to be refined are referred to as "parent" cells, whereas those resulting from the splitting of the "parent" cells are the "offspring".Note that this work exclusively makes use of block-structured grids consisting purely of hexahedral cells (square elements of the same size/area in 2D).
In the beginning of this step, all background grid cells are set to have a uniform refinement level (L 0 ).Then, cells with β values in the user-defined [β min , β max ] range are selected for refinement.This range defines the bandwidths of cells close to the FSI to be refined, even though the β field is not a distance field.β min = 0.3 and β max = 0.7 are used in the examples presented here in order to refine the cells in a broad band across the FSI.Then, cells are added to or subtracted from the previously selected set in order to meet the grid regularity requirements.In particular, cells are marked or unmarked for refinement in such a way that keeps the level difference of adjacent cells equal to 1, at most; this is imposed for stability reasons.The cells to be refined are split into 4 (in 2D) or 8 (in 3D) offspring cells, each having a refinement level of L P + 1, where L P is the level of their parent cell.An octree (in 3D) or quadtree (in 2D) data structure is used to store the information regarding parent and offspring cells.The cell selection and refinement process described above is repeated several times in each optimization cycle, until a user-defined maximum refinement level L max is reached.At the end of each refinement step, the β field is interpolated to the offspring cells, and the selection criterion for new refinement cells is applied again.In this work, a trilinear polynomial of the following form is constructed for each refined parent cell in terms of the coordinates (x, y, z) in the local coordinate system and is used to compute the β values at the offspring cell centers.At each vertex i of the parent cell, βi =⃗ p T i •⃗ c should hold.Here, βi denotes the interpolated β value from adjacent cells to vertex i.In a grid with hexahedra, this leads to the formation of an 8 × 8 linear system of equations of the form A •⃗ c = ⃗ β.By inverting A, the unknown coefficients can be expressed as a linear combination of βi , i.e., c i = α ij βj , and Equation (2) can be written as Equation (3) requires the β values at grid vertices prior to the refinement.These are computed by interpolating β from cells to vertices, according to where Cells(p) is the set of cells sharing vertex p and ω pC = are inverse-distance interpolation weights.Differentiating Equation (4) with respect to β at any cell Q yields if Q belongs on the list of cells sharing p as a vertex, or 0 otherwise.This will later be used as part of the chain rule to compute δJ δα .

Generation of Cut Cells
The β field, interpolated at the refined grid vertices though Equation (4), is used to compute the FSI and create body-fitted grids at each optimization cycle.The FSI reflects the intersection of grid edges by the β w = 0.5 iso-surface (or isoline in 2D).An edge connecting nodes A and B with values βA < β w < βB is intersected at Based on Equation ( 6), the derivative of x S i with respect to β at any node M is computed as where δ is the Kronecker delta.Equation ( 7) is also part of the chain rule computing δJ δα .
Cut cells are created by connecting the computed nodes along the edges intersected by the iso-surface to define new internal grid faces.The node connectivity is computed as in [42,43].The new internal faces split the initial cells into sub-cells belonging to Ω F and Ω S .This is shown schematically in Figure 2, where the introduced internal faces correspond to the dark-gray edges.In practice, a cell might be intersected more than once by the β w iso-surface; see Figure 2c.Surface areas, centers, and normals for the newly added faces are computed using the x S i coordinates.The same quantities are re-computed for intersected faces belonging to Ω F or Ω S using the defined x S i nodes on these faces and the nodes of the background face with β < β w or β > β w , respectively.Finally, the cell centers and volumes of the so-formed cut cells are computed.
In practice, cut cells of small size can be generated.Their presence next to cells of much larger size might cause numerical instabilities during the flow solution, leading to convergence issues.To deal with this, each small cell is merged with a neighboring larger cell if its volume is less than 1% of that of the latter [44], as shown in Figure 3.More than one small cell can be merged with the same neighboring cell, by also updating the grid connectivity information.Further discussion of this matter is beyond the scope of this work and can be found in many papers on conventional cut-cell methods [45][46][47].Compared to level-set-based methods [48], the approach proposed here to track the FSI is simpler, since it avoids the re-initialization of the signed distance field in each optimization cycle by solving, for instance, the Hamilton Jacobi equation.Since the solution of the flow (see Section 2.4) and adjoint (see Section 2.6) equations takes place only in Ω F , by discarding all cells belonging to Ω S , the computational burden is reduced compared to that in standard density-based TopO.In other words, the systems of equations emerging from the discretization of the governing PDEs have as many cells as those in Ω F , with the exception of the temperature equation, which is solved simultaneously in both Ω F and Ω S , i.e., over the entire domain Ω d .

Governing Equations
The governing equations for the flow and CHT problem are presented; their boundary conditions are summarized in Appendix A along with those of the adjoint problem.The flow solver (and its adjoint) is implemented in the in-house adjointOptimisation library [49] of OpenFOAM; the latter utilizes a finite-volume discretization and a cell-centered storage strategy.Details regarding the numerical solution of the governing PDEs follow.

Flow Equations
The flow is governed by the Navier-Stokes equations for incompressible fluids, which, in residual form, are written as where p is the pressure divided by the constant fluid density and v i is the velocity.
is the stress tensor, where ν is the bulk and ν t the eddy viscosity.The later is modeled by the Spalart-Allmaras turbulence model [50] equation solved to compute ν.In Equation (10), P( ν, y) = c b1 S and D( ν, y) = c w1 f w ν y 2 are the production and dissipation terms of the model, respectively.The turbulence model functions are given by , g = r+c w2 r 6 − r , r = ν Sκ 2 y 2 and its constants are and c v1 = 7.1.Finally, the eddy viscosity is computed as ν t = f v1 ν.To avoid extreme grid stretching close to the wall, the wall function technique, as in Spalding's formula [51], is used.
The distances y from the nearest wall, including the FSI, are computed via the solution of the Hamilton-Jacobi equation [52], reading in which the diffusion term has been multiplied by (1+ϵ), with ϵ = 0.1, for stability reasons.Equations ( 8)-( 10) are solved in a segregated manner based on the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) [53], to handle the pressure-velocity coupling.Second-order upwind schemes are used to discretize the convection terms.Second-order central schemes are used for Laplace terms, using an explicit correction term to take into account the grid's non-orthogonality.Finally, gradients are computed by applying Gauss' theorem over control volumes; linear interpolation from cell to face centers is used to evaluate the field's values at their boundaries.

Temperature Equation
The temperature field in Ω F is modeled by the advection-diffusion equation, Here, ρ F is the fluid density, and c F p is the fluid-specific heat capacity.k e f f = k F + k t is the effective thermal conductivity, which includes a laminar k F and a turbulent where Pr t is the turbulent Prandtl number.The assumption that the fluid's thermophysical properties do not depend on the temperature is made.In Equation ( 12), Q is a volumetric heat source term that can be either constant or given in the form of h T re f − T , where h is a constant heat transfer coefficient and T re f a reference temperature.This difference in the form of Q (and the corresponding physics) is the main distinction between the two cases presented in Section 4.
In Ω S , heat conduction (R T S = 0) is modeled by the same equation, by replacing k e f f with k S (solid thermal conductivity), while considering zero velocity, guaranteed by the exact boundary conditions over the FSI, due to the use of cut cells.
The discretization of Equation ( 12) using finite volumes assumes integration over the boundaries of each cell, which requires k f e f f ( f = f ace) to be available at the grid faces.For the internal faces of Ω F and Ω S , the linear interpolation of k e f f from cells to faces is used.However, special treatment is required at the FSI to accurately compute the diffusive fluxes q f .There, the heat flux conservation property implies that q f = −q f F =q f S , where q f F , q f S are the outgoing heat fluxes from the fluid and solid parts.In discrete form, where T f is the temperature at the FSI; T F , T S denote the temperature values at the adjacent fluid and solid cell centers, respectively; and ∆ F f , ∆ S f are the distances of the cell centers in Ω F and Ω S , adjacent to the FSI, from the latter.By solving Equation ( 13) for T f and substituting it back to either q f F or q f S , the heat flux at the interface can be written in terms of T F and T S as Equation ( 14) is consistent with the way in which the boundary conditions are imposed at the interface between adjacent domains in CHT simulations based on a partitioned approach [54].

Objective Functions and Constraints
In this work, TopO is applied to improve the performance of cooling devices.In the presented examples (see Section 4), a volumetric heat source term Q is applied over Ω S and this causes the T over the solid to increase.Due to the high thermal conductivity of the solid material, the applied heat is then transferred to the fluid/coolant and is removed through the domain outlet.Depending on the way in which heat is provided to the system, the TopO aims either to minimize the average temperature over Ω S or to maximize the amount of heat transferred from the solid to the fluid.The latter is the difference in the convecting energy fluxes evaluated at the inlet (S I ) and the outlet (S O ), i.e., Note that, due to the way in which the normal vectors are defined, the two integrals on the right-hand side of Equation ( 16) have opposite signs.An increase in cooling performance usually comes at the expense of higher power consumption in order to maintain the flow.The energy loss due to friction forces at the solid walls is given by the volume flowrate-weighted total pressure (p t = p+ 1 2 v 2 k ) drop between the inlet and outlet, which should not exceed a threshold value.This constraint is expressed as where P 0 L is the power dissipation computed for the starting geometry of the TopO and π P L is the percentage value of the threshold.A second inequality constraint is imposed on the fluid volume so as not to exceed a target percentage (π V ) of the analysis domain (Ω d ) volume, as follows:

Continuous Adjoint Method for Sensitivity Analysis
To update the design variables α in each optimization cycle, the derivatives of the objective and the constraint functions with respect to α are computed using the continuous adjoint [10].The formulation of the adjoint problem for a general objective function or constraint, F , starts by augmenting it with the field integrals of the governing equations multiplied by the adjoint variables to form the Lagrangian, Ł(≡ F ), i.e., where q is the adjoint pressure, u i is the adjoint velocity, ν a is the adjoint turbulence variable, y α is the adjoint distance, and T α is the adjoint temperature.
To compute the total variation in F due to changes in the shape of the interface, Ł is differentiated with respect to the coordinates x s m of the points on the FSI, with δŁ δx s m = δF δx s m .Following the terminology used in shape optimization, these derivatives practically stand for the sensitivity map of F [55] on the FSI computed by the cut-cell method.Since the coordinates x s m depend on the cell center β(α) values, the chain rule of differentiation is applied to compute the derivatives of F with respect to α.
To derive the adjoint equations and their boundary conditions, δŁ δx s m is expanded using Gauss' theorem and the identity δ , relating spatial derivatives to total derivatives with respect to x s m and holds for any quantity Φ.This mathematical development for CHT problems (governed by Equations ( 8)-( 12)) can be found in [57] (therein, however, for shape optimization) and is not repeated here in the interest of brevity.δŁ δx s m is written as The undefined terms in Equation ( 21) will become clear as the analysis proceeds.
To avoid the computation of δv i where is the adjoint stress tensor and Q α is the adjoint heat source term.Terms C S , C ν , and C y emerge from the differentiation of the production and dissipation terms in Equation ( 10) with respect to S, ν, and y and can be found in [58].Eliminating the same derivatives along the boundaries of Ω F and Ω S gives rise to the adjoint boundary conditions, summarized in Appendix A.
The solution of the adjoint problem (Equations ( 22)-( 26)) starts by first solving the adjoint temperature equation (Equation ( 26)) over the entire Ω d to compute T α .Equation ( 14) is used to evaluate the adjoint thermal diffusive flux at the FSI, by replacing T with T α .Then, the adjoint mean flow (Equations ( 22) and ( 23)) and turbulence model (Equation ( 24)) PDEs are solved in a segregated manner.In this work, the discretization of the adjoint equations, described in Section 2.7, is consistent with the one used for the primal problem (as in the discrete adjoint), leading to the computation of accurate sensitivity derivatives.
After satisfying the field adjoint equations and their boundary conditions, the remaining terms in Equation ( 21) are the shape sensitivity derivatives along the FSI, where Note that Equation ( 27) requires the grid sensitivities δx k δx s m obtained by differentiating the coordinates of face centers along the FSI and the cut-cell centers with respect to x s m ; the latter are given by closed-form expressions that can be found in [43].
where the first summation is over the N nodes along the FSI, whereas the second is over all NP nodes at the edges cut by the FSI.The derivatives δx sn m δ βQ , δ βQ δβ P are obtained through Equations ( 5) and (7), respectively.Using the proposed method, non-zero derivatives with respect to β are computed only at cells close to the FSI.From this point of view, the proposed TopO method resembles level-set -based rather than the standard density-based TopO.As a final step, the filtering operation R F Ψ β , δF δβ (see Section 2.1) is applied to δF δβ to yield the derivatives with respect to α at each cell P of the background mesh, viz.δF δα P = Ψ β P Ω P , where Ω P is the cell volume.

"Think Discrete-Do Continuous" Adjoint Discretization
To compute the gradient of the discretized objective with respect to the design variables as accurately as in the discrete adjoint, the continuous adjoint PDEs (Equations ( 22)-( 26)) should be discretized in a manner consistent with the primal equations.In this work, this is achieved by the "Think Discrete-Do Continuous" (TDDC) adjoint, firstly introduced by the present group of authors in [38] for use in aerodynamic shape optimization.It is now extended to TopO with CHT.The derivation of the TDDC adjoint starts from the derivation of a (hand-differentiated) discrete adjoint.
The method is exemplified through the derivation of the consistent discretization of the convection term in the adjoint temperature equation (Equation ( 26)), which emerges from the differentiation of term ρ F c F p v j ∂T ∂x j (Equation ( 12)).Integrating the convection term of Equation ( 12) over the finite volume Ω P , using the Gauss divergence theorem and a second-order upwind interpolation, yields where v PN ⟨n⟩ is the normal to the face f convecting velocity computed by the Rhie-Chow interpolation [59], ∆S f are the face areas, and w PN U , w PN U = 1 − w PN U are the upwind interpolation weights computed based on v PN ⟨n⟩ .Moreover, d P f j is the vector connecting the cell centroid to the face center.In Equation (31), summation is performed over the faces F(P) of the finite volume.The spatial gradients of T are computed based on Gauss' theorem and linear interpolation from cells to face centers, i.e., ∂T ∂x j P = 1 where are linear interpolation weights, computed using geometric distances.Here, Φ T,PN FO , Φ T,PN SO denote the first-and second-order convecting fluxes at the faces.The discrete form of the adjoint convection term −ρ F c F p where FO and SO have the same meaning as before.According to the proposed TDDC adjoint, the adjoint fluxes in Equation ( 33) should be computed as follows for the discrete form of Equation ( 33) to be consistent with its primal counterpart (Equation ( 31)).The first-order adjoint flux in Equation ( 34) is similar to the corresponding first-order primal one, Φ TiPN FO , the only difference being that the upwind interpolation weights, w PN U , w PN U , have switched places.This suggests that a downwind scheme should be employed for the interpolation of T α at the faces in conformity with the upwind propagating information in the adjoint problem.Additional differences can be seen in the second-order fluxes, in which the adjoint gradients are multiplied by the linear interpolation weights and are projected onto the face normal vectors, rather than the face-to-cell distance ones.
The spatial gradients in Equation ( 34) are discretized as follows ∂T α ∂x j P where the fluxes in both terms are computed in a non-conservative manner (i.e., fluxes computed on the same face are different when seen from the two cells adjacent to that face, due to the presence of d P f j ).Despite this, the consistent second-order adjoint flux (Equation ( 34)) is conservative; this is in line with the fact that the adjoint energy equation (Equation ( 26)) is conservative too.
The full derivation of Equations ( 34) and ( 35) is given in Appendix B.

Comparison with Density-Based TopO (denTopO)
A sharp distinction between the proposed TopO workflow and that of the density-based approach (denTopO) is the ability of the former to apply proper boundary conditions along the FSI and generate (and adapt/refine) body-fitted grids.On the other hand, denTopO utilizes Brinkman penalization terms to simulate the effect that solidified areas have on the flow.These terms, expressed in terms of the material distribution field β, are added to Equations ( 9)-(11) (and their adjoints) and effectively treat solid walls as regions of high impermeability.For a flow variable ϕ, the corresponding Brinkman term is given by BS ϕ = β max I(β)ϕ, where β max stands for the impermeability value of the solid material and takes high enough values to ensure that ϕ ≈ 0 inside Ω S .Its value is selected based on the non-dimensional Darcy number (set here to Da = 10 −5 ) expressing the ratio of viscous to porous forces as follows where L is a characteristic length.I(β) = β 1+b I (1−β) is a convex function [7] and b I a steepening parameter, which penalizes intermediate density values, driving them closer to 0. In the examples of Section 4, whenever denTopO is utilized, b I = 10.The main disadvantage of this approach is that, in practice, the imposed penalization terms are not sufficient to impose exact zero values of the corresponding flow field inside Ω S .For instance, for ϕ = v i in BS ϕ , spurious unwanted flow leakage inside the solid areas, hindering the accuracy of the flow solution, usually occurs.Moreover, it does not support the use of the wall function technique in turbulent flows.To solve the CHT problem, the thermal conductivity k at the cell centers is interpolated between k F and k S , using the β field, as k(β) = k F +β p (k S − k F ), with p = 3.The interpolation scheme used for k allows intermediate values to appear.Then, Equation ( 12) computes T over Ω d using k e f f = k(β)+(1−β)k t as the effective thermal conductivity.Note that denTopO is unable to accurately interpolate k e f f at the FSI (Equation ( 14)), which is crucial in simulating heat transfer phenomena.
The second difference between cut-cell TopO and denTopO is related to the computation of the adjoint sensitivities with respect to α.The cut-cell TopO is driven by shape sensitivities, which are, then, transformed into derivatives with respect to α.These are used to update the design variables and, subsequently, the FSI's position.Depending on the r value used in Equation ( 1), δJ δα in cut-cell-based TopO takes on non-zero values at cells in the vicinity of the FSI.In denTopO, the derivative of J with respect to α at cell P is given by These take on non-zero values everywhere in Ω d , allowing thus more drastic topological changes between optimization cycles.

Application in CHT Problems
The proposed method is demonstrated in two different CHT TopO problems.The two topologies, having a single (CHT-A) and a split (CHT-B) outlet, can be found in Figure 4, with dimensions given in terms of the inlet height H = 0.1 m.CHT-A was first defined in [19].The design variables α take on a zero value inside the inlet and outlet channels (Figure 4) and are allowed to vary only in the central part of the geometry (shown in blue color in Figure 4; this stands for the design domain Ω D ), where a solid material (in gray color) can be added or subtracted.For both cases, a Cartesian background grid is generated, with 21.6 K square cells of edge length equal to h = H 20 .Even if both cases are symmetric with respect to the top horizontal line, the entire domain is shown in most figures presenting results.
The working fluid is air, with a constant density ρ F = 1.18 kg/m 3 , specific heat capacity c F p = 1004.9J/kg•K, thermal conductivity k F = 0.0254 W/m•K, and bulk viscosity ν = 1.5 × 10 −5 m 2 /s.Moreover, Pr t = 1.The solid material corresponds to aluminum with k S = 237 W/m•K.The fluid enters the domain at the temperature of T in = 273 K; the exit pressure is set to zero.The case Reynolds number is defined as Re = v in H ν , in terms of the inlet velocity magnitude v in .Unless specified otherwise, the inlet velocity is v in = 0.375 m/s, or Re = 2500.
In CHT-A, the sought solid domain is heated by a constant volumetric heat source Q = 1 kW/m 2 .TopO seeks to minimize the average temperature over Ω S (J T , Equation ( 15)).In CHT-B, the volumetric heat source term applied over Ω S is given by Q = h T re f − T , where h = 10 W/m 2 K and T re f = 373 K. Since Q is not fixed, the goal is to maximize the heat transferred to the coolant by maximizing J H (Equation ( 16)).In both cases, the fluid volume (area in 2D) inside Ω D should be less than 45% of Ω D (π V ≈ 0.49 in Equation ( 19)).Power dissipation for CHT-A is controlled by setting π P L = 0.33 in Equation ( 18), whereas three different π P L values (0.25, 0.33, and 0.5) are used in CHT-B to investigate its effect on the cooling performance (see Section 4.4).Regarding the MMA settings used in both cases, the adaptation parameters controlling the amount of the asymptote decrease or increase in each optimization cycle are set to 0.5 and 1.05, respectively, instead of the recommended (in [39]) values of 0.7 and 1.2.In addition, ρ = 5 × 10 −5 (as opposed to the default value of 1 × 10 −5 ) is used to increase the convexity of the approximated J , g P L , and g V .

CHT-A Results Using denTopO
For the purpose of comparison, CHT-A is firstly optimized using the "standard" denTopO approach with Brinkman terms in the flow equations (see Section 2.4).Using the inlet height H as the characteristic length in Equation (36), β max = 150.As denTopO cannot trace the FSI, there is no clear boundary separating Ω F and Ω S .Thus, J T (Equation ( 15)) and g V (Equation ( 19)) are recast as follows where integrations over Ω d (instead of Ω F and Ω S ) are performed by considering the β field.In addition, in Equation ( 12), the imposed heat source term Q is multiplied by β to activate the term in the solidified part only.Note that, in this problem, as the solid volume changes during the optimization, the total amount of heat given to the solid, Ω d βQdΩ, is not fixed.The volume constraint ensures that a minimum amount of heat is given to the system and avoids trivial all-fluid solutions.
Figure 5 shows the convergence of the normalized J T value and the imposed constraints.At the beginning of denTopO, both g V and g P L are violated.The first ∼20 cycles are spent to meet the inequality constraints; then, the MMA goes on to reduce J T to 335.4 K.The optimized distribution of β is shown in the bottom half of Figure 6a.Even though the computed β field in the optimized solution is almost binary, thin gray areas between the fluid and solid appear.The top half of Figure 6a corresponds to the optimized geometry, for which the FSI has been computed by post-processing the β field to generate cut cells.On the background grid that is adapted to the FSI, the optimized solution is re-evaluated using the cut-cell code.This re-evaluation gives J T = 447.6K, with g P L = −0.14.This means that, when re-evaluated with the cut-cell primal solver, the denTopO solution still has a lot of room before approaching the P L threshold, justifying the poor thermal performance.In terms of execution time, the cut-cell flow and adjoint solvers are ∼38% and ∼30% faster than the corresponding Brinkman solvers, for the same amount of SIMPLE iterations, without the usage of AMR, and on the same hardware.This reduction in computational cost is expected, considering that Ω S occupies ∼51% of Ω d .Isolines for the non-dimensional velocity magnitude , the intermediate Spalart-Allmaras model variable χ = ν ν , and the temperature T using the Brinkman and the cut-cell-based flow solvers are shown in Figure 6b-d, respectively.As indicated by the iso-velocity contours of Figure 6b, the reevaluation on the cut-cell grid results in a different distribution of the incoming flow within the various channels.Due also to the inability of the denTopO flow solver to impose accurate boundary conditions on the FSI, significant discrepancies in the computed temperature fields between the two flow solvers can be seen, as in Figure 6d.This study highlights the importance of performing TopO based on solvers able to accurately compute physics.As shown in Section 4.3, a better-performing (in terms of cooling) design can be obtained using the proposed methodology.The top half is the same (optimized using denTopO) solution that has been post-processed to generate cut cells and re-evaluated by imposing accurate FSI conditions.

Verification of Computed Sensitivities Using TDDC Adjoint
In this subsection, the computed sensitivity derivatives using the proposed TDDC adjoint of Section 2.7 are validated against finite differences (FDs) for CHT-A.To reduce the cost of the FDs, which depends directly on the number of design variables, a volumetric NURBS-based parameterization is used.In particular, a 2D lattice of NX × NY nodes controls the α values in cells lying within its boundaries; the latter coincide with the boundaries of Ω D (Figure 7).The α CP values at the control points are the design variables of the problem, instead of the α values of the background grid cells.Since NX = 25, NY = 13, this leads to 325 design variables in total.A set of parametric coordinates (u, v) is assigned to each cell P, the center coordinates (x P , y P ) of which lie within the control lattice; the parametric coordinates at P are computed as u P = x P −x min x max −x min , v P = y P −y min y max −y min and the α value as (39) where U i , V j are the second-degree NURBS basis functions.Taking the derivative of Equation ( 39) with respect to α lm CP yields, δα P The control lattice is shown in Figure 7a, whereas the resulting geometry after transforming the α field into cut cells is presented in Figure 7b; the velocity field is shown inside Ω F , whereas Ω S is in gray.The derivatives of J T with respect to the α cell values are computed by the adjoint method, as described in Section 2.6.The chain rule of differentiation transforms δJ T δα P into δJ T δα lm CP .The so-computed derivatives for the first 100 design variables are plotted in Figure 7c and are in perfect agreement with those computed by the FDs.This verifies the accuracy of the sensitivity derivatives of J when the TDDC adjoint discretization schemes are used.The same level of accuracy is obtained for the derivatives with respect to the remaining design variables, which are not shown here.It should be noted that the lattice shown in Figure 7 is used only in this comparison, to keep the cost of the FDs as low as possible, and not in any other study that follows.

CHT-A Results Using Cut-Cell-Based TopO
Having verified the accuracy of the computed derivatives using the TDDC adjoint, the latter is used by the cut-cell-based TopO to optimize the cooling in CHT-A.To initialize the case, Ω D has been seeded with multiple square solid structures of size H 2 ; see Figure 8a.First, the importance of using the grid refinement strategy of Section 2.2 is demonstrated by solving the TopO problem with 0, 1, and 2 refinement steps.The three solutions will be referred to as R0, R1, and R2, respectively; see Figure 9. Small but significant differences exist among the three optimized geometries.Since the design variables of TopO are defined on the background grid, the unknown optimization variables are the same in all three cases; therefore, the differences in the optimized geometries can be attributed only to differences in the flow solution's accuracy due to the grid refinement.At the end of the refinement process (if applied), the grid used to solve the flow consists of 21.6 K, ∼62 K, and ∼140 K cells for R0, R1, and R2, respectively; a close-up view of the refined grid close to the FSI walls in the central area of the three optimized geometries is given in Figure 9.The maximum y + values computed at the cell centers adjacent to the FSI are also reported; as expected, lower values are computed as the grid resolution close to the FSI increases.According to the J T values listed in Figure 9, each computed using the grid utilized during TopO, the three optimized designs perform equally well in terms of cooling.For a fair comparison, solutions R0 and R1 are re-evaluated after refining the background grid twice using the already optimized β field.The computed J T and ∆p t values are listed in Table 1.The power dissipation values of (the re-evaluated) R0 and R1 solutions exceed the maximum allowed value, which means that the g P L constraint is not met.Irrespective of the violation or not of this constraint, solution R2 outperforms R0 and R1 in terms of J T .The current study confirms the need to refine the grid during TopO.The Reynolds number's effect on the optimized solution is also investigated.The optimization problem is solved for three different inlet velocity values, namely v in = 0.15 m/s, 0.375 m/s, and 0.75 m/s, corresponding to Re = 1000, 2500, and 5000, respectively.Based on the previous study, the grid is refined twice in each cycle of TopO and this leads to the maximum y + values at cells adjacent to the FSI equal to 3, 5, and 8, respectively.
The optimized solutions for the three Re are presented in Figure 10.The average temperature values in Ω S are also given.Note that the solution presented in the middle column of Figure 10 is identical to R2, as in Figure 9.The power dissipation value P 0 L of the initial design (Figure 8) increases as Re increases and so does the maximum allowed power dissipation, according to Equation (18).In all optimized designs, a number of solid bodies, forming a network of flow channels inside the design domain for improved cooling, can be observed.For Re = 1000, larger solid structures, aligned with the flow, are created.By increasing Re, more solid structures of smaller sizes, and thus more and narrower flow channels, are formed.Thus, the heat dissipation is enhanced and this is reflected in the average temperature values being equal to J T = 547.25 K, 378 K, and 340.4K for Re = 1000, 2500, and 5000, respectively.These observations are in agreement with those in [19,20].Figure 11 shows the convergence of the normalized J T value and of the imposed constraints for the three optimizations of Figure 10.In the beginning of the TopO, the fluid volume is more than 45% of Ω D and, thus, g V is violated.As the solid material is added to meet g V , more heat is given to the system and this leads to higher J T values.Once g V has been satisfied, an almost constant amount of heat, Q = 275 W, is imposed in all remaining cycles.After satisfying both g V and g P L , the MMA tries to minimize the value of J T ; see Figure 11.Since a further decrease in the fluid's volume and ∆p t would be counterproductive with respect to the cooling of CHT-A, both g V and g P L are practically satisfied as equality constraints (g V , g P L ≈ 0).Max y + = 10 Max y + = 6.5 Max y + = 5 J T = 378K J T = 378.4KJ T = 378.6K

CHT-B Results Using Cut-Cell-Based TopO
The ability of the proposed method to generate optimal topologies is demonstrated in a second configuration, in which the outlet duct has been re-located to the bottom right of Ω D ; see Figure 4b.Compared to CHT-A, the solid material is heated up according to the following volumetric heat source term: Q = h T re f − T .TopO seeks to maximize the heat transferred to the fluid.The latter is expressed by J H (Equation ( 16)), which, after taking into consideration the energy conservation, is equivalent to Ω S h T re f − T dΩ.As in CHT-A, the initial design has been seeded with multiple square solid structures; see Figure 8b.
The influence of the maximum allowed power dissipation inside the heat exchanger is investigated by performing the optimization for three different π P L values (Equation ( 18)), namely π P L = 0.25, π P L = 0.33, and π P L = 0.5.The background grid is refined twice in each TopO cycle.The optimized solutions are presented in Figure 12.Increasing π P L , and thus the allowed power dissipation, leads to the formation of more flow channels for enhanced heat transfer.The average temperature inside Ω S for the three designs shown in Figure 12 is J T = 327.9K,323.4K, and 318.6K.As expected, lower values of J T are obtained by increasing π P L as this allows more heat to be exchanged.

Conclusions
A TopO framework that combines a material distribution field (design variables) and the cut-cell method to tackle laminar and turbulent flow problems with CHT is proposed.Here, the design variables are defined at the cell centers of a background grid and their post-processing yields an auxiliary field β that helps to identify fluid and solid areas, as well as the FSI as the β w = 0.5 iso-surface.The FSI's intersections with the background grid generate cut cells, giving rise to a body-fitted grid on which the flow equations are solved by applying proper boundary conditions, including wall functions.As cells belonging to Ω S are discarded, the solution cost of the flow and adjoint equations is reduced (by up to ∼38% in the examined cases).Additionally, the proposed workflow is able to locally refine the grid close to the FSI to capture the flow physics.The proposed method replaces the Brinkman penalization terms in the flow equations, frequently met in TopO, with the imposition of exact boundary conditions on the FSI and, thus, avoids spurious flow leakage inside Ω S .
To facilitate TopO, the continuous adjoint method is derived for the CHT problem and computes the sensitivity map of J on the FSI.The chain rule is applied on the computed shape sensitivities, yielding the derivatives with respect to the material field.This is the mechanism that allows topological changes to be made.This work also focuses on the derivation of consistent discretization schemes for the adjoint PDEs.This derivation, guided by a hand-differentiated discrete adjoint, under the name "Think Discrete-Do Continuous" or the TDDC adjoint, enables the computation of consistent gradients for the discretized objective function (as in the discrete adjoint) and provides physical meaning to the underlying adjoint PDEs and a low memory footprint (as in the continuous adjoint).
The cut-cell TopO workflow was demonstrated in two 2D cases having a single (CHT-A) and a split (CHT-B) outlet, both operating in turbulent flow conditions.CHT-A was first optimized using the density approach for TopO, relying on the usage of Brinkman penalization terms and interpolation schemes for the thermal conductivity.The re-evaluation of the optimized geometry using the developed cut-cell flow solver showed large discrepancies in the objective and constraint values.In the same case, the accuracy of the proposed TDDC adjoint solver in computing the gradient of J with respect to α was verified against FDs.The use of the cut-cell-based TopO method highlighted the importance of grid refinement to accurately compute the performance during TopO and obtain solutions that truly meet the imposed constraints.CHT-A and CHT-B were optimized for different Re numbers and with different threshold values for the allowed power dissipation, respectively, reaffirming that higher values of both lead to the creation of more flow channels and, eventually, better cooling performance.
The proposed method was demonstrated in 2D cases of academic interest drawn from the literature, to allow for several parametric studies at a low computational cost.The corresponding s/w for 3D flow/optimization problems has already been implemented and its application will be the focus of subsequent research.In Equation (A8), ρ F c F P v PN ⟨n⟩ w PN U T αP +w PN U T αN ∆S f stands for the consistent first-order adjoint convective flux, as given in Equation (34).
The differentiation of the second-order flux term R T α SO,P with respect to T P gives gives rise to the so-called adjoint gradient term given by Equation (35).To proceed with the derivation, the discretization of ∂T ∂x j using Gauss' theorem (Equation ( 32)) is differentiated with respect to T P , and this yields where, considering also that the normal vector n PN j points outwards from Ω P , the identity ∑ f ∈F(P) n PN j ∆S f = 0 has been used to obtain the last expression.From Equation (A9), it can be seen that the term within the summation is identical to the second-order consistent adjoint flux given in Equation (34).

Figure 1 .
Figure 1.Flowchart of the cut-cell-based TopO algorithm.

Figure 2 .
The two-step generation of cut cells.(a) (Step 1): β is interpolated at each vertex q using the β values at the neighboring cells Q. (b,c) (Step 2): The x S i positions where grid edges are cut by the β w = 0.5 iso-surface are computed.Cut cells are generated by connecting these nodes to create new internal faces that split the initial cell into fluid (blue) and solid (red) parts.

Figure 3 .
Figure 3. Merging a small-volume cut cell belonging in Ω F with a neighboring one.
and solid domains, field integrals including them in Equation (21) are eliminated by setting their multipliers to zero, giving rise to the field adjoint equations (a) CHT-A (b) CHT-B Figure 4.A single-outlet (a) and a split-outlet (b) 2D topology [19].The geometry and boundary conditions imposed at the inlet and outlet.The external walls are adiabatic.The available design space, Ω D , where TopO can add or subtract solid material, is the central rectangular area shown in blue.H = 0.1 m.

Figure 5 .
Figure 5. denTopO of a single-outlet heat exchanger (CHT-A): convergence of the normalized objective function J T J 0 T

Figure 6 .
Figure 6.CHT-A results using denTopO.The bottom half of each domain is the outcome of denTopO.The top half is the same (optimized using denTopO) solution that has been post-processed to generate cut cells and re-evaluated by imposing accurate FSI conditions.

Figure 7 .
Figure 7. CHT-A: verification of the TDDC adjoint sensitivity derivatives of J T against FDs for the cut-cell TopO solver.(a) The 25 × 13 lattice controlling the α field.The control points are colored by their α CP values, which take on user-defined values close to 0/1.(b) The channel geometry (Ω F ) created by transforming α CP into cut cells and iso-velocities.(c) Sensitivity derivatives of J T with respect to α CP computed by the TDDC adjoint method and FDs.

Figure 9 .Figure 10 .
Figure 9. CHT-A results using the cut-cell-based TopO: optimized designs at Re = 2500 using no grid refinement (left), one (center), and two grid refinement steps (right).Rows from top to bottom correspond to optimized geometries, grid density over a selected area, and ∥v i ∥ v in and T isolines.

Figure 11 . 0 T
Figure 11.CHT-A: convergence of the normalized objective function J T J 0 T

Figure 12 .
Figure 12.CHT-B results using the cut-cell-based TopO: optimized designs for three different ∆p t threshold values.Rows from top to bottom correspond to ∥v i ∥ v in and T isolines.
p = 0, u i t i = 0 BC u i t i = 0,∂(u i n i ) ∂x j n j = 0 BC p = 0, u i t i = 0 ∂(u i t i ) ∂x j n j = 0, u i n i = 0 y α y α = 0 y α = 0 ofConsistent Discretization Schemes for the Adjoint Convection Term of the Temperature EquationIn this appendix, the formulas for the consistent adjoint fluxes, given by Equation(34), to the second-order upwind discretization of the energy equation convection term are derived.In the discrete adjoint, the discrete residual R T α P of the T α equation at cell P is given by 31) as the discrete residual of the primal equation,The term R T α FO,P , including the derivatives of first-order fluxes, is first developed:⟨n⟩ w PN U (T αP −T αN )∆S f (A7)By adding and subtracting ρ F c F p v PN ⟨n⟩ T αP ∆S f , and using w PN U = 1 − w PN U as well as the mass fluxes over all faces of Ω P that sum up to zero, Equation (A7) becomesR T α FO,P =− ∑ f ∈F(P) ρ F c F p v PN ⟨n⟩ w PN U T αP +w PN U T αN ∆S f (A8)

Table A1 .
Boundary conditions for the primal and adjoint fields.