1. Introduction
Since its introduction in [
1] for use in solid mechanics and structures, TopO has been widely adopted in a variety of scientific fields, targeting the optimal material layout that minimizes the given objective function under certain design constraints. TopO formulations for fluid mechanics include phase-field approaches [
2], level set methods [
3,
4,
5,
6], porosity/density-based methods [
7,
8], etc., with applications to a wide range of problems, including Stokes flows [
7], laminar [
8] and turbulent flows [
9,
10,
11], unsteady flows [
12,
13], reactive flows [
14], natural convection problems [
15,
16], etc. The use of TopO for fluid problems with CHT first appeared in [
17,
18] and still remains an active area of research [
19,
20,
21,
22], including also the design of bi-fluid heat exchangers [
23,
24,
25]. An overview of the state-of-the-art TopO methods for fluid flow problems, including heat transfer, was recently published [
26].
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.
2. The Cut-Cell TopO Method
The solution of the TopO problem starts with the generation of a background grid for the analysis domain,
. 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
correspond to solids and those with
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
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,
, and solid,
, domains. The governing flow equations, without the Brinkman terms usually met in TopO, are solved only in
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
and
. The continuous adjoint problem is, then, formulated for the objective function
to be minimized and each flow-related constraint
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.
2.1. 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
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.
2.2. 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 (
). Then, cells with
values in the user-defined
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.
and
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
, where
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
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
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,
should hold. Here,
denotes the interpolated
value from adjacent cells to vertex
i. In a grid with hexahedra, this leads to the formation of an
linear system of equations of the form
. By inverting
A, the unknown coefficients can be expressed as a linear combination of
, i.e.,
, 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
is the set of cells sharing vertex
p and
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
.
2.3. 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
iso-surface (or isoline in 2D). An edge connecting nodes
A and
B with values
is intersected at
Based on Equation (
6), the derivative of
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
.
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
and
. 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
iso-surface; see
Figure 2c. Surface areas, centers, and normals for the newly added faces are computed using the
coordinates. The same quantities are re-computed for intersected faces belonging to
or
using the defined
nodes on these faces and the nodes of the background face with
or
, 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
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
, by discarding all cells belonging to
, 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
, with the exception of the temperature equation, which is solved simultaneously in both
and
, i.e., over the entire domain
.
2.4. 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.
2.4.1. 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
is the velocity.
is the stress tensor, where
is the bulk and
the eddy viscosity. The later is modeled by the Spalart–Allmaras turbulence model [
50] equation
solved to compute
. In Equation (
10),
and
are the production and dissipation terms of the model, respectively. The turbulence model functions are given by
and its constants are
,
,
,
,
,
,
, and
. Finally, the eddy viscosity is computed as
. 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
, with
, 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.
2.4.2. Temperature Equation
The temperature field in
is modeled by the advection–diffusion equation,
Here,
is the fluid density, and
is the fluid-specific heat capacity.
is the effective thermal conductivity, which includes a laminar
and a turbulent
part, where
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
, where
h is a constant heat transfer coefficient and
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 , heat conduction () is modeled by the same equation, by replacing with (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
(
) to be available at the grid faces. For the internal faces of
and
, the linear interpolation of
from cells to faces is used. However, special treatment is required at the FSI to accurately compute the diffusive fluxes
. There, the heat flux conservation property implies that
, where
,
are the outgoing heat fluxes from the fluid and solid parts. In discrete form,
where
is the temperature at the FSI;
,
denote the temperature values at the adjacent fluid and solid cell centers, respectively; and
,
are the distances of the cell centers in
and
, adjacent to the FSI, from the latter. By solving Equation (
13) for
and substituting it back to either
or
, the heat flux at the interface can be written in terms of
and
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].
2.5. 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
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
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 (
) and the outlet (
), 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 (
) drop between the inlet and outlet,
which should not exceed a threshold value. This constraint is expressed as
where
is the power dissipation computed for the starting geometry of the TopO and
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 (
) of the analysis domain (
) volume, as follows:
2.6. 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,
, starts by augmenting it with the field integrals of the governing equations multiplied by the adjoint variables to form the Lagrangian,
, i.e.,
where
q is the adjoint pressure,
is the adjoint velocity,
is the adjoint turbulence variable,
is the adjoint distance, and
is the adjoint temperature.
To compute the total variation in
due to changes in the shape of the interface,
is differentiated with respect to the coordinates
of the points on the FSI, with
. Following the terminology used in shape optimization, these derivatives practically stand for the sensitivity map of
[
55] on the FSI computed by the cut-cell method. Since the coordinates
depend on the cell center
values, the chain rule of differentiation is applied to compute the derivatives of
with respect to
.
To derive the adjoint equations and their boundary conditions,
is expanded using Gauss’ theorem and the identity
[
56], relating spatial derivatives to total derivatives with respect to
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.
is written as
The undefined terms in Equation (
21) will become clear as the analysis proceeds.
To avoid the computation of
,
,
,
,
in the fluid 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
where
is the adjoint stress tensor and
is the adjoint heat source term. Terms
,
, and
emerge from the differentiation of the production and dissipation terms in Equation (
10) with respect to
,
, and
y and can be found in [
58]. Eliminating the same derivatives along the boundaries of
and
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
to compute
. Equation (
14) is used to evaluate the adjoint thermal diffusive flux at the FSI, by replacing
T with
. 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
obtained by differentiating the coordinates of face centers along the FSI and the cut-cell centers with respect to
; the latter are given by closed-form expressions that can be found in [
43].
Having computed
, the chain rule of differentiation helps in computing
at cell
P as
where the first summation is over the
nodes along the FSI, whereas the second is over all
nodes at the edges cut by the FSI. The derivatives
,
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
(see
Section 2.1) is applied to
to yield the derivatives with respect to
at each cell
P of the background mesh, viz.
, where
is the cell volume.
2.7. “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
(Equation (
12)). Integrating the convection term of Equation (
12) over the finite volume
, using the Gauss divergence theorem and a second-order upwind interpolation, yields
where
is the normal to the face
f convecting velocity computed by the Rhie–Chow interpolation [
59],
are the face areas, and
,
are the upwind interpolation weights computed based on
. Moreover,
is the vector connecting the cell centroid to the face center. In Equation (
31), summation is performed over the faces
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.,
where
are linear interpolation weights, computed using geometric distances. Here,
,
denote the first- and second-order convecting fluxes at the faces.
The discrete form of the adjoint convection term
integrated over a control volume is
where
and
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,
, the only difference being that the upwind interpolation weights,
,
, have switched places. This suggests that a downwind scheme should be employed for the interpolation of
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
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
). 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.
3. 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
, where
stands for the impermeability value of the solid material and takes high enough values to ensure that
inside
. Its value is selected based on the non-dimensional Darcy number (set here to
) expressing the ratio of viscous to porous forces as follows
where
L is a characteristic length.
is a convex function [
7] and
a steepening parameter, which penalizes intermediate density values, driving them closer to 0. In the examples of
Section 4, whenever
denTopO is utilized,
. 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
. For instance, for
in
, 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
and
, using the
field, as
, with
. The interpolation scheme used for
k allows intermediate values to appear. Then, Equation (
12) computes
T over
using
as the effective thermal conductivity. Note that
denTopO is unable to accurately interpolate
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),
in cut-cell-based TopO takes on non-zero values at cells in the vicinity of the FSI. In
denTopO, the derivative of
with respect to
at cell
P is given by
These take on non-zero values everywhere in
, allowing thus more drastic topological changes between optimization cycles.
4. 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
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
), where a solid material (in gray color) can be added or subtracted. For both cases, a Cartesian background grid is generated, with
K square cells of edge length equal to
. 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 kg/m3, specific heat capacity J/kg·K, thermal conductivity W/m·K, and bulk viscosity m2/s. Moreover, . The solid material corresponds to aluminum with W/m·K. The fluid enters the domain at the temperature of K; the exit pressure is set to zero. The case Reynolds number is defined as , in terms of the inlet velocity magnitude . Unless specified otherwise, the inlet velocity is m/s, or .
In CHT-A, the sought solid domain is heated by a constant volumetric heat source
kW/m
2. TopO seeks to minimize the average temperature over
(
, Equation (
15)). In CHT-B, the volumetric heat source term applied over
is given by
, where
W/m
2K and
K. Since
Q is not fixed, the goal is to maximize the heat transferred to the coolant by maximizing
(Equation (
16)). In both cases, the fluid volume (area in 2D) inside
should be less than
of
(
in Equation (
19)). Power dissipation for CHT-A is controlled by setting
in Equation (
18), whereas three different
values (
,
, and
) 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
and
, respectively, instead of the recommended (in [
39]) values of
and
. In addition,
(as opposed to the default value of
) is used to increase the convexity of the approximated
,
, and
.
4.1. 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),
. As
denTopO cannot trace the FSI, there is no clear boundary separating
and
. Thus,
(Equation (
15)) and
(Equation (
19)) are recast as follows
where integrations over
(instead of
and
) 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,
, 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
value and the imposed constraints. At the beginning of
denTopO, both
and
are violated. The first ∼
cycles are spent to meet the inequality constraints; then, the MMA goes on to reduce
to
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
K, with
. This means that, when re-evaluated with the cut-cell primal solver, the
denTopO solution still has a lot of room before approaching the
threshold, justifying the poor thermal performance. In terms of execution time, the cut-cell flow and adjoint solvers are ∼
and ∼
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
occupies ∼
of
. 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 re-evaluation 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.
4.2. 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
nodes controls the
values in cells lying within its boundaries; the latter coincide with the boundaries of
(
Figure 7). The
values at the control points are the design variables of the problem, instead of the
values of the background grid cells. Since
,
, this leads to 325 design variables in total. A set of parametric coordinates (
) is assigned to each cell
P, the center coordinates
of which lie within the control lattice; the parametric coordinates at
P are computed as
,
and the
value as
where
,
are the second-degree NURBS basis functions. Taking the derivative of Equation (
39) with respect to
yields,
.
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
, whereas
is in gray. The derivatives of
with respect to the
cell values are computed by the adjoint method, as described in
Section 2.6. The chain rule of differentiation transforms
into
. 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
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.
4.3. 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,
has been seeded with multiple square solid structures of size
; 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
K, ∼
K, and ∼
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
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
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
and
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
constraint is not met. Irrespective of the violation or not of this constraint, solution R2 outperforms R0 and R1 in terms of
. 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 m/s, m/s, and m/s, corresponding to , 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 values at cells adjacent to the FSI equal to 3, 5, and 8, respectively.
The optimized solutions for the three
are presented in
Figure 10. The average temperature values in
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
of the initial design (
Figure 8) increases as
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
, larger solid structures, aligned with the flow, are created. By increasing
, 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
K, 378 K, and
K for
, 2500, and 5000, respectively. These observations are in agreement with those in [
19,
20].
Figure 11 shows the convergence of the normalized
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
of
and, thus,
is violated. As the solid material is added to meet
, more heat is given to the system and this leads to higher
values. Once
has been satisfied, an almost constant amount of heat,
W, is imposed in all remaining cycles. After satisfying both
and
, the MMA tries to minimize the value of
; see
Figure 11. Since a further decrease in the fluid’s volume and
would be counterproductive with respect to the cooling of CHT-A, both
and
are practically satisfied as equality constraints (
).
4.4. 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
; see
Figure 4b. Compared to CHT-A, the solid material is heated up according to the following volumetric heat source term:
. TopO seeks to maximize the heat transferred to the fluid. The latter is expressed by
(Equation (
16)), which, after taking into consideration the energy conservation, is equivalent to
. 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
values (Equation (
18)), namely
,
, and
. The background grid is refined twice in each TopO cycle. The optimized solutions are presented in
Figure 12. Increasing
, and thus the allowed power dissipation, leads to the formation of more flow channels for enhanced heat transfer. The average temperature inside
for the three designs shown in
Figure 12 is
,
, and
. As expected, lower values of
are obtained by increasing
as this allows more heat to be exchanged.
5. 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 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 are discarded, the solution cost of the flow and adjoint equations is reduced (by up to ∼ 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 .
To facilitate TopO, the continuous adjoint method is derived for the CHT problem and computes the sensitivity map of 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 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.