XFEM level set-based topology optimization for turbulent conjugate heat transfer problems

Solving conjugate heat transfer design problems is relevant for various engineering applications requiring efficient thermal management. Heat exchange between fluid and solid can be enhanced by optimizing the system layout and the shape of the flow channels. As heat is transferred at fluid/solid interfaces, it is crucial to accurately resolve the geometry and the physics responses across these interfaces. To address this challenge, this work investigates for the first time the use of an eXtended Finite Element Method (XFEM) approach to predict the physical responses of conjugate heat transfer problems considering turbulent flow. This analysis approach is integrated into a level set-based optimization framework. The design domain is immersed into a background mesh and the geometry of fluid/solid interfaces is defined implicitly by one or multiple level set functions. The level set functions are discretized by higher-order B-splines. The flow is predicted by the Reynolds Averaged Navier–Stokes equations. Turbulence is described by the Spalart–Allmaras model and the thermal energy transport by an advection–diffusion model. Finite element approximations are augmented by a generalized Heaviside enrichment strategy with the state fields being approximated by linear basis functions. Boundary and interface conditions are enforced weakly with Nitsche’s method, and the face-oriented ghost stabilization is used to mitigate numerical instabilities associated with the emergence of small integration subdomains. The proposed XFEM approach for turbulent conjugate heat transfer is validated against benchmark problems. Optimization problems are solved by gradient-based algorithms and the required sensitivity analysis is performed by the adjoint method. The proposed framework is illustrated with the design of turbulent heat exchangers in two dimensions. The optimization results show that, by tuning the shape of the fluid/solid interface to generate turbulence within the heat exchanger, the transfer of thermal energy can be increased.


Introduction
Solving conjugate heat transfer problems is relevant for various engineering applications. For systems operating at large or small scales, such as gas turbines or micro-electronic modules, an efficient thermal management is crucial to ensure reliable functionality. Heat exchangers are often added to systems to provide an efficient temperature control. They are made of highly conductive materials and are cooled or heated through natural or forced convection, i.e., a moving fluid is used to transport the heat to or away from the systems. By modifying the layout of heat exchangers, the flow path and the flow conditions can be altered which in turn influences the heat transfer between fluid and solid. Therefore, the performance in terms of heat transfer and dissipative losses can be improved by optimizing the geometry of the heating or cooling devices.
Since the seminal work of Bendsøe and Kikuchi (1988), topology optimization has become a popular tool to systematically address design problems and generate layouts with enhanced performance under specific requirements, see Sigmund and Maute (2013) and Deaton and Grandhi (2014). Originally introduced for structural designs, topology optimization was first applied to fluid problems by Borrvall and Petersson (2003). Considering Stokes flows, they introduced a porosity variable to represent the distribution of fluid and solid regions through the design domain. The method was extended to laminar Navier-Stokes models by Gersborg-Hansen et al. (2006). Following these pioneering works, various flow design problems have been solved by topology optimization, see Alexandersen and Andreasen (2020) for a comprehensive review.
While most designs are optimized assuming laminar flow, a wide range of fluid systems operate in the turbulent regime. Accurately modeling turbulence at an acceptable computational cost remains challenging, see Spalart (2000). High fidelity models, such as Direct Numerical Simulations (DNS) or Large Eddy Simulations (LES), are time-consuming and integrating them into an optimization loop remains impractical. Lower fidelity models based on time averaged velocity and pressure fields, such as the Reynolds Averaged Navier-Stokes (RANS) model, are commonly used. To represent the Reynolds stresses, they are generally coupled to eddy viscosity models, such as the Spalart-Allmaras (SA), the k − , or the k − model, see Spalart and Allmaras (1992) and Wilcox (1993). To date, the design of fluid systems for turbulent flows is mainly based on RANS models.
Working with finite volumes, Othmer (2008) proposed a continuous adjoint approach with respect to density design variables and studied turbulent flows in ducts under the assumption of frozen turbulence, i.e., the variations of the turbulence variables are not accounted for in the sensitivity analysis. Zymaris et al. (2009) extended the continuous adjoint approach with respect to shape parameters to account for the turbulence variations modeled by the SA equation. They showed that, considering frozen turbulence, sensitivities can exhibit erroneous signs. Philippi and Jin (2015) developed an adjoint sensitivity method where the turbulent flow in porous media is represented by a k − model. They used finite volumes to design channels for minimal mechanical energy loss under the assumption of frozen turbulence. Papoutsis-Kiachagias and Giannakoglou (2016) leveraged the continuous adjoint approach to perform shape and density-based topology optimization on aero-and hydrodynamic industrial applications. Working with finite elements, Yoon (2016) proposed a density-based optimization framework for the design of channels with minimal turbulent energy dissipation. The framework relies on a discrete adjoint sensitivity analysis of the RANS equations closed with the SA model. The approach was further extended to the k − turbulence model in Yoon (2020). Dilgen et al. (2018a) designed channels and manifolds for turbulent flows with minimal power dissipation using a density-based topology optimization approach. They used a finite volume approach and automatic differentiation to perform the discrete adjoint sensitivity analysis of the RANS equations with different eddy viscosity models. Sá et al. (2021) tackled the specific case of rotating flows with an adapted version of the SA model and using a density-based approach.
While all aforementioned studies rely on density-based approaches, other optimization techniques have also been applied to turbulent flow design problems. Kubo et al. (2021) proposed a level set-based approach with Ersatz material and finite volumes. They represented the turbulence with the k − and k − models under the assumption of frozen turbulence. Picelli et al. (2022) used a binary topology optimization approach in combination with a geometry trimming procedure to design channels for minimum turbulent energy dissipation with the k − and k − models. Alonso et al. (2022) extended the framework to the Wray-Agarwal turbulence model. To improve the predictive capabilities of the RANS equations, Hammond et al. (2022) used data-driven turbulence modeling and designed fluid systems with turbulent flows under the frozen turbulence assumption.
Along with the development of topology optimization approaches for fluid problems, multi-physics fluid problems gained in popularity, in particular conjugate heat transfer problems. Initial works on density-based methods considering conjugate heat transfer by Dede (2009) and Yoon (2010) focused on the design of channels for forced convection with laminar flows. Over the years, the complexity has increased to cover different heat transfer phenomena, different flow conditions, and to include multiple fluids. The interested reader is referred to Dbouk (2017) for an extensive overview.
As turbulence allows for increased heat transfer between fluid and solid phases as well as for an improved heat transport through the fluid, several research efforts have been devoted to the development of topology optimization frameworks for heat transfer with turbulent flows. Kontoleontos et al. (2013) extended the continuous adjoint approach of Zymaris et al. (2009) to account for turbulent heat transfer and generated channel layouts maximizing the energy transfer for a constrained pressure loss using density-based topology optimization. Koga et al. (2013) designed heat sinks for maximum heat transfer and minimum pressure loss considering Stokes flow discretized by finite volumes. In a post-design stage, they carried out validation simulations and considered turbulent flows with a k − model. Pietropaoli et al. (2017) used densitybased topology optimization and a discrete adjoint approach to design internal channels heated at the outer walls. Turbulence was described by a k − model under the assumption of frozen turbulence. Dilgen et al. (2018b) extended their work on turbulent flow topology optimization to include thermal transport and design heat sinks for forced convection with turbulent flow modeled by the k − model. Lee et al. (2020) developed a simplified non-exact adjoint sensitivity analysis for the design of aero-thermal systems including internal and external turbulent flows based on the SA model and finite volumes. Yaji et al. (2020) proposed a multi-fidelity design approach for thermal turbulent flow optimization problems. Zhao et al. (2021) used a multi-layer thermo-fluid model for the design of turbulent forced convective heat sinks using the SA model. Ghosh et al. (2022) focused on the design of internal cooling ducts for gas turbine blades using finite volumes and under the frozen turbulence assumption.
All previously cited works on turbulent conjugate heat transfer rely on density-based methods. As the key physics phenomenon, i.e., the heat transfer, happens at the fluid/solid interface and as the flow fields are strongly influenced by the shape of the fluid/solid interface, it is crucial that both the geometry and the physics responses across the interface are accurately resolved. Combining the level set method (LSM) with immersed finite element methods (IFEM) offers an elegant approach to tackle such problems. Several fluid design problems were studied using level set-based topology optimization with IFEM, such as hydrodynamic Boltzmann transport in , fluid-structure interactions in Jenkins and Maute (2016), heat transfers via natural convection in Coffin and Maute (2016), and species transport in Villanueva and Maute (2017). So far, design problems focusing on turbulent conjugate heat transfer have not been addressed using such techniques.
In this paper, a level set-based topology optimization framework for turbulent conjugate heat transfer problems is proposed. To obtain a crisp description of the fluid/solid interface, the design domain is immersed into a background mesh, and the design geometry is defined implicitly by one or multiple level set functions (LSF). For the first time, an eXtended Finite Element Method (XFEM) approach is used to predict the physical responses of the systems. To model turbulent conjugate heat transfer, the flow is described by the RANS equations for incompressible flows and stabilized by the Streamline Upwind Petrov-Galerkin (SUPG) and the Pressure Stabilizing Petrov-Galerkin (PSPG) formulations. Turbulence is described by the one-equation SA model. As the turbulence explicitly depends on the distance to the closest wall, a distance field is constructed based on the heat method, see Crane et al. (2017). The thermal energy transport in the fluid phase is predicted by an advection-diffusion equation. The temperature field in the solid phase is governed by a linear diffusion equation.
To discretize and solve conjugate heat transfer problems, an XFEM approach is adopted. This approach is similar to the XIGA approach presented in Noël et al. (2022) but is restricted to first-order approximations for the state variable fields. Higher-order B-spline functions are used for the geometry representation as was done in Noël et al. (2020). A generalized Heaviside enrichment strategy with multiple enrichment levels is used to ensure independent approximation on each connected fluid or solid subregion. Boundary and interface conditions are enforced weakly via Nitsche's method, see Nitsche (1971). Additionally, the face-oriented ghost stabilization is used to mitigate numerical instabilities resulting from the creation of small integration subdomains, see Burman (2010) and Burman and Hansbo (2014).
Design problems are solved by a gradient-based optimization algorithm. The gradients of the objective and constraint functions are evaluated by the adjoint method. The optimization framework is applied to the design of heat exchangers considering turbulent flows to enhance heat transfer and thermal management with acceptable pressure losses.
The remainder of the paper is structured as follows. Section 2 is devoted to the geometry representation of the designs using one or several LSFs. In Sect. 3, the basic concepts of the XFEM approach such as the enrichment strategy, the face-oriented ghost stabilization, and the numerical integration scheme are briefly described. Section 4 details the physics model used to predict the solution of turbulent conjugate heat transfer problems. The handling of the governing equations using the XFEM approach is detailed. The ability of the proposed approach to solve turbulent conjugate heat transfer problems is demonstrated in Sect. 5 with benchmark problems, namely the backward facing step and the flow over a thick conducting plate. The formulation of the optimization problems is discussed in Sect. 6. In Sect. 7, the proposed optimization framework is illustrated with two dimensional heat exchanger design problems. Variations of the design problem are considered in terms of the outer geometry of the heat exchangers, of the imposed flow conditions, and of the formulation of the optimization problems. Finally, Section 8 draws conclusions on the developed framework and proposes extensions of the framework for future research.

Geometry representation
The LSM was introduced by Osher and Sethian (1988) to efficiently track propagating fronts. Soon after its introduction, the LSM has been used in combination with shape and topology optimization, see Sethian and Wiegmann (2000), Wang et al. (2003), Allaire et al. (2004), andvan Dijk et al. (2013) for an overview. The method describes geometries implicitly by LSFs. The iso-level 0 of a LSF ( ) defines an interface or a boundary Γ ± that separates the domain Ω into two subdomains, Ω + and Ω − , such that: This paper follows the work by Vese and Chan (2002) on a multi-phase level set approach. One or multiple LSFs i ( ) with i = 1, … , n are used to describe the design geometry. With n LSFs, 2 n subregions can be identified within the domain Ω . Each subregion is characterized by a unique combination of the LSF signs, describing whether a point lies inside, outside, or on the interface. Each subregion is associated with the fluid, the solid, or the void phase, and their respective properties. An illustration of the geometry description of a fluid/solid/void problem is presented in Fig. 1. The design domain Ω is defined by the dashed line. The fluid, the solid, and the void domains are represented in blue, grey, and white and are denoted by Ω f , Ω s , and Ω v , respectively. The symbols Γ f and Γ s denote the interfaces of the fluid and solid domains with the void domain, respectively. The fluid/solid interface is denoted by Γ fs , i.e., The level set description of the design domain geometry for a fluid/solid/void problem is shown in Fig. 2. Eight LSFs in green, 1 ( ), … , 8 ( ) , are used to define the boundaries of the design domain that is fully immersed in a background analysis domain. The fluid/solid interfaces within the design domain are described by one LSF 9 ( ) in red. The minus sign associated to each isoline i ( ) = 0 indicates the region where i ( ) < 0 .
For numerical analysis, each LSF i ( ) is discretized on a mesh and approximated as: where B k ( ) are B-spline basis functions and ( i ) k are the coefficients associated with the LSF i ( ) , similarly to Noël et al. (2022). In contrast to advancing the LSFs in the optimization process through the solution of the Hamilton-Jacobi equation, see approaches by Sethian and Wiegmann (2000), Wang et al. (2003), and Allaire et al. (2004), the coefficients of the LSFs are expressed as explicit functions of the design variables . They are updated by mathematical programming methods driven by shape sensitivities.

XFEM formulation
IFEMs alleviate the need for generating conforming analysis meshes. IFEMs using enrichment functions enable the representation of discontinuities within a mesh element by augmenting the standard finite element approximation with enrichment functions. The XFEM was proposed by Moës et al. (1999) and Belytschko and Black (1999) to model crack propagation without remeshing. Over the years, the method has been successfully applied to various types of problems including the multi-phase problems considered here, see for example Burman and Hansbo (2014), Schott and Wall (2014), and Schott et al. (2015).
In this paper, an XFEM approach is adopted, which allows for accurate and efficient resolution of multi-phase problems with evolving interfaces and boundaries. The approach is similar to the XIGA proposed in Noël et al. (2022) but restricts the approximation of the state variable fields to first-order functions. The main ingredients of the XFEM approach are briefly recalled in this section. First, the enrichment strategy is described in Subsect. 3.1. Subsection 3.2 introduces the face-oriented ghost stabilization used to mitigate instabilities resulting from the generation of small integration subdomains. Finally, the numerical integration scheme used to accommodate elements occupied by multiple phases is explained in Subsect. 3.3.

Enrichment strategy
In this paper, the enrichment strategy described in  is extended to enrich each individual basis function separately based on the topology of the phase layout within the basis function support as described in Noël et al. (2022). Considering a multi-phase problem, a state field ( ) is approximated as: where K is the number of background basis functions and L k is the number of separate connected subregions {Ω k } L k =1 in the support of the background basis function B k . The coefficients a k are the degrees of freedom (DOFs) associated with the basis function B k and phase subregions Ω k . The indicator function k ( ) determines whether a point belongs to a phase subregion Ω k : The selection of the enrichment levels is illustrated for a fluid/solid problem in Fig. 3. A basis function B k spans both the fluid and the solid phases and its support is delimited by a red-dashed line. Five separate connected material subregions Ω k exist within the basis support. Each subregion is occupied by one and only one phase, such that are occupied by the fluid phase and two subregions Ω =4 k and Ω =5 k are occupied by the solid phase. As five subregions Ω k exist within the basis support, five enrichment levels L k = 5 are necessary.

Face-oriented ghost stabilization
Enriched IFEMs may suffer from numerical instabilities when boundaries and interfaces intersect the mesh such that small integration subdomains are generated. In this situation, the contributions of the basis functions approximating the state variable fields might vanish and lead to a poorly conditioned system of equations. To mitigate this issue, the face-oriented ghost stabilization penalizes the jump in state field gradients across so-called ghost facets, as proposed in Burman (2010) and Burman and Hansbo (2014). In this work, the approach of Noël et al. (2022) is adopted.
The set of ghost facets F G is the set of facets that lie inside the design domain and that are intersected by interfaces and boundaries. A ghost facet F is shared between two adjacent elements Ω + F and Ω − F , and is characterized by a normal F , defined so as to align with the outward normal to Ω + F and F = + F = − − F . The elements adjacent to the facet F are subdivided by the phase layout into N + F and N − F connected subdo- respectively to all of ℝ n d with n d the spatial dimension. The jump is defined similarly for the trial field . The operator k n (•) is the k th order normal derivative operator and n (•) = ∇(•) ⋅ F with ∇(•) the spatial derivative operator. The ghost penalty parameters G , used in this paper, are further described in Subsect. 4.6.  Figure 4 illustrates the subdivision of Ω + F and Ω − F into connected subdomains Ω + F,i and Ω − F,j and the associated indices M + F,i and M − F,j for a fluid/solid problem. The considered ghost facet F and the adjacent background elements are marked using a red box. The background element Ω − F is occupied by three connected subdomains Ω − F,1 , Ω − F,2 , and Ω − F,3 while the background element Ω + F is divided into two connected subdomains Ω + F,1 and Ω + F,2 . As Ω − F,2 and Ω + F,1 have the same material index M − F,2 = M + F,1 and meet along the facet, the jump in the associated field gradients is penalized across the facet. The same holds for the pairs Ω − F,1 and Ω + F,2 , and Ω − F,3 and Ω + F,2 .

Numerical integration
The XFEM approach allows for the presence of multiple solid, fluid, and void phases within a single mesh element. The integration of the weak form of the governing equations is performed separately on each subdomain. For this purpose, a conforming integration mesh is built for each intersected element. A quadrangular intersected element is subdivided into triangular integration elements that align with the boundaries and interfaces described by the LSFs. It should be noted that the cost for generating the conforming XFEM models is negligible in comparison to the cost of the forward and sensitivity analyses. Standard quadrature is then performed on each triangular integration element and non-intersected background element. The reader is referred to Villanueva and Maute (2014) and Noël et al. (2022) for further details on the subdivision strategy.

Incompressible Navier-Stokes equations
In this work, the fluid flow is governed by the incompressible Navier-Stokes equations. The residuals of the strong form of the momentum equilibrium equations Ω and of the incompressibility condition r p Ω are given as: where is the fluid velocity, p is the pressure, and is the constant fluid density. The Cauchy stress tensor is denoted by ( , p) and is computed as:

Fig. 4
Ghost stabilization on a fluid/solid problem where is the identity matrix, is the fluid dynamic viscosity, and ( ) is the strain rate tensor defined as: The residuals of the weak form of the momentum equations and of the continuity condition within the fluid domain are given as: where and are the velocity field and test function, and p and p are the pressure field and test function, respectively.
Boundary conditions on the fluid velocity field are imposed weakly via Nitsche's formulation, see Nitsche (1971) and Bazilevs and Hughes (2007). To prescribe a velocity D on Γ f D , the following contributions are added to the velocity and pressure residuals: where N and p N are parameters used to obtain a symmetric ( N = 1 , p N = 1 ) or a skew-symmetric ( N = −1 , p N = −1 ) Nitsche's formulation. In this paper, a skew-symmetric formulation is adopted for both contributions. The Nitsche's penalty parameter is denoted N and is defined following Schott et al. (2015) as: where N is a constant parameter chosen to achieve a desired accuracy in satisfying the boundary conditions. The first term accounts for viscosity-dominated flows and h is the length of the interface within an intersected element, while the second term accounts for convection flows. The operator ‖ • ‖ ∞ is the infinite norm operator, such that At the inflow, imposing inlet velocity conditions might cause inflow instabilities. To alleviate this issue, an additional up-winding term is added to the residual Γ D : where up is the up-winding parameter and is set to up = 1 , as proposed in Schott and Wall (2014). The pressure p N can be imposed on the fluid boundaries Γ f N . The weak form of the residual for the Neumann boundary condition is evaluated as:

Spalart-Allmaras turbulence model
In this paper, heat transfer with turbulent flows is considered. Turbulence introduces additional stresses in the fluid due to the presence of turbulent eddies in the flow. These eddies generate additional shear stress, i.e., the so-called Reynolds stress t . This stress is modeled based on the Boussinesq eddy viscosity assumption as: where t is the turbulent dynamic viscosity. The stress tensor takes the form: The RANS system of equations is closed using the turbulence model of Spalart and Allmaras (1992) and the turbulent dynamic viscosity is evaluated as: where ̃ is the so-called modified viscosity and is the additional state variable introduced in the SA model. The residual of the strong from of the SA equation is given as: where P is the production term, D the wall destruction term, and K the diffusion coefficient. The modified velocity ̃ is defined as: where the constants c b2 and are set to c b2 = 0.622 , and = 2∕3.
The original model presented in Spalart and Allmaras (1992) only admits non-negative values for the modified viscosity ̃ given non-negative boundary and initial conditions. However, coarse grids and transient states may yield negative values. To handle these situations, a continuation model for negative turbulence values was proposed in Allmaras et al.
(2012). The modified model recovers the original SA model for positive ̃ values, while a negative value of the modified viscosity ̃ produces zero eddy viscosity, i.e., t = 0 . The modified version of the model is adopted in this paper and is described hereunder in the context of the XFEM. The turbulent dynamic viscosity t is evaluated as: where is the kinematic viscosity such that = , and c v1 is a constant set to c v1 = 7.1.
The production term P is defined as: where f t2 is a parameter that controls the use of the f t2 term, i.e., f t2 = 1 to account for and f t2 = 0 to exclude the f t2 term from the formulation. The f t2 term was initially introduced to control the laminar regions and to render solutions with ̃= 0 stable. Rumsey (2007) showed that the f t2 term can be omitted for reasonably high Reynolds number. In this paper the influence of omitting the f t2 term on the optimization results is investigated. The constant parameter P is used to enhance the numerical behavior of the model when ̃ becomes negative and should be set such that P ≥ 1.0 , see Anderson et al. (2019). In this paper, the parameter is set to P = 10.0 . The constants c b1 , c t3 , and c t4 are set to c b1 = 0.1355 , c t3 = 1.2 , and c t4 = 0.5 , respectively. The wall destruction term D is defined as: where d is the distance to the closest wall. The constants and c w3 are set to = 0.41 and c w3 = 2 . The intermediate variable g is computed as: , and =̃, where c w2 and r lim are constants set to c w2 = 0.3 and r lim = 10.0. The diffusion coefficient K is defined as: where c n1 is a constant set to c n1 = 16.
The value of the modified vorticity S should always be positive and greater than 0.3 S for physically relevant configurations. Numerically, the model allows for zero or negative values of S . To ensure the non-negativity of the modified vorticity S , the following modifications were proposed in Allmaras et al. (2012): where the magnitude of the vorticity S is evaluated as: and the intermediate variable S as: Finally, the residual of the weak form of the SA equation is given as: where ̃ and ̃ are the modified viscosity field and test function, respectively.
As the flow at the fluid/solid interface is close to zero, i.e., it is laminar, the viscosity is prescribed to zero, i.e., ̃= 0 . A prescribed modified viscosity ̃D is imposed weakly via Nitsche's formulation as: where Ñ is a parameter used to obtain a symmetric ( Ñ = 1 ) or a skew-symmetric ( Ñ = − 1 ) Nitsche's formulation. In this paper, a skew-symmetric formulation is adopted. The Nitsche's penalty parameter Ñ is defined as: where Ñ is a constant penalty parameter chosen to achieve a desired accuracy in satisfying the boundary conditions. At the inflow, enforcing inlet viscosity conditions might cause instabilities. Similarly to Eq. (14), an additional upwinding term is added to the residual R̃Γ D : where the up-winding parameter ũ p is set to ũ p = 1.0.

Advection-diffusion equation
The heat transfer in the fluid is described by an advection-diffusion equation where the velocity is governed by the flow model. The residual of the strong form of the advection-diffusion equation is given as: where c p is the heat capacity, and Q is the body heat load. The diffusive heat flux is denoted by and is computed as: where is the isotropic diffusion tensor and = with the conductivity. The presence of turbulence results in an additional turbulent heat flux t defined as: where t is the turbulent isotropic diffusion tensor and t = t with t , the turbulent conductivity. Finally, the diffusive heat flux takes the form: The turbulent conductivity t is evaluated as: where Pr t is the turbulent Prandtl number and is set to Pr t = 0.9. The residual of the weak form of the advection-diffusion equation is given as: where T and T are the temperature field and test function, respectively.
For heat conduction in the solid, a linear diffusion model is used and can be obtained by omitting the advection term in Eqs. (32) and (53) and replacing the conductivity of the fluid, = f , with the one of the solid, = s .
Boundary conditions on the temperature field are prescribed weakly via Nitsche's formulation: where T D is a prescribed temperature, T N is a parameter used to obtain a symmetric ( T N = 1 ) or a skew-symmetric ( T N = − 1 ) Nitsche's formulation. A skew-symmetric formulation is implemented in this paper. The Nitsche's penalty parameter T N is expressed as: where T N is a constant parameter chosen to achieve a desired accuracy in satisfying the boundary conditions.
The continuity of temperature and heat flux at the fluid/ solid interface Γ fs is imposed weakly via Nitsche's formulation as follows: where T I is a parameter used to obtain a symmetric ( T I = 1 ) or a skew-symmetric (  where meas (Ω f ) and meas (Ω s ) are the surface area of fluid and solid within the intersected element. The penalty parameter T I is evaluated as: where meas (Γ fs ) is the length of the interface within the intersected element.

Wall distance through the heat method
The one-equation SA turbulence model requires the distance to the closest wall d. In this paper, the wall distance is constructed using the heat method, see Crane et al. (2017). First, a transient heat conduction problem for the additional state variable field is solved over the design domain. The residual of the strong form of the conduction equation is written as: where , c p , and are the density, the specific heat capacity, and the isotropic diffusion tensor of the heat method. Zero initial boundary conditions, i.e., | |t=0 = 0.0 , are imposed on the design domain Ω and a fixed temperature, i.e., | |Γ fs = 1.0 , is enforced at the interface Γ fs between fluid and solid. The equation is solved in one single time step. The time step size is selected following the work by Geiss et al. (2019), i.e., just large enough to yield non-zero gradient over the design domain.
The distance D is constructed by solving a Poisson's equation with a body heat load dependent on the normalized gradient of the additional temperature field . The strong form of the residual for the Poisson's equation is expressed as: A distance D | |Γ fs = 0.0 is enforced on the interface Γ fs between fluid and solid.
The weak form of the residuals in Eqs. (43,44) are formulated similarly to the advection-diffusion equations presented in Subsect. 4.3. It should be noted that building a wall distance field using the heat method requires the introduction of two additional state variable fields and D and the solution of two additional scalar partial differential equations.
In the vicinity of and at the walls, the obtained distance field D might exhibit values close to machine precision and suffer from small spurious node-to-node variations. This issue results from the weak imposition of zero distance boundary conditions within the XFEM approach. The small variations might be amplified through the formulation of the SA source terms and lead to large nodeto-node variations of the production and wall destruction coefficients, which in turn jeopardize the stability of the SA model. To mitigate this issue and regularize the distance field, a L 2 projection scheme is used. First, a logarithmic function of the square of the distance D is projected to a intermediate distance variable P , such that: where P and P are the intermediate distance field and the test function, respectively. Finally, the wall distance d is computed from the intermediate distance variable P as:

Subgrid stabilization
Spurious node-to-node velocity oscillations can arise from the convective terms in the incompressible Navier-Stokes equations, the SA turbulence equation, and the advectiondiffusion equation, see Eqs. (7,18,32). Additionally, using approximations of the same order for both the velocity and the pressure fields may lead to spurious pressure oscillations. To mitigate these issues, the incompressible Navier-Stokes equations are augmented by the SUPG and the PSPG formulations of Tezduyar et al. (1992). The SUPG formulation introduced in Franca et al. (1992) is used to augment both the SA turbulence equations and the advection-diffusion equation.
The contribution of the SUPG/PSPG formulations to the weak form of the residuals of the incompressible Navier-Stokes equations is given as: Following the approach of Taylor et al. (1998) and Whiting and Jansen (2001), the stabilization parameters SUPG and p PSPG are defined as: where C I is a constant set to 36.0 for linearly interpolated finite elements, see Whiting (1999). The operator tr (•) evaluates the trace as tr (•) = • ii and the contravariant metric tensor is expressed as: with are the parametric coordinates in the background element.
The weak form of the residual for the turbulence equation is augmented with the SUPG formulation as: The stabilization parameter S UPG is evaluated following Tezduyar and Osawa (2000) as: where ‖ • ‖ is the Frobenius norm and h̃ is a length measure defined as: where | • | is the absolute value operator and B i are the basis functions used to approximate the velocity field.
Finally, the weak form of the residual for the advectiondiffusion equation is augmented with the SUPG formulation as: The stabilization parameter T SUPG is formulated, similarly to Eq. (50), as: where the length measure h is defined as:

Face-oriented ghost stabilization
To counteract the effect of numerical instabilities related to small integration subdomains, the face-oriented ghost stabilization presented in Subsect. 3.2 is adopted in this work. The formulation of the residual of the weak form of the stabilization for a generic state field is given as: where  (2014), the associated ghost penalty parameter ,visc G is given as: where ,visc G is a constant penalty parameter, h F is the length associated with ghost facet F , and p is the order of the considered approximation, as introduced in Eq. (5).
A convective ghost penalty term ,conv F G acting on the first order gradient of velocity field only, as defined in Eqs. (56, 5) with p = 1 , is added to the residual Ω as proposed in Schott and Wall (2014). The associated ghost penalty parameter ,conv G is given as: where ,conv G is a constant penalty parameter. To control pressure instabilities, a ghost penalty term R p F G is added to the pressure residual R p Ω , as proposed in Schott et al. (2015), and the associated ghost penalty parameter p G is defined as: where both viscous and convective flows are handled by the first and second term respectively and p G is a constant penalty term. Instabilities in the modified viscosity field ̃ are mitigated by adding a ghost stabilization term to the viscosity residual R̃Ω . The associated ghost penalty parameter is given as: where G is a constant penalty parameter.
Finally, a ghost penalization term is added to the temperature residual R T Ω and the ghost penalty parameter T G is defined as: where T G is a constant penalty parameter. Similar ghost penalization terms are introduced to stabilize the additional temperature , distance D , and approximate distance P fields.

Validation examples
To validate the physics model presented in Sect. 4, two wellknown academic examples are reproduced using our XFEM approach and are compared against numerical results from the literature. All verification examples are analyzed with the flow channels immersed in a background mesh to demonstrate the accuracy of the proposed numerical models for settings encountered in the optimization studies of Sect. 7.
Studying turbulent flows, the physics responses are often presented in terms of normalized quantities. In the following examples, the normalized distance to the wall y + , the normalized velocity u + , and the normalized temperature T + are used and are expressed as: The friction velocity u and the friction temperature T are defined as: where w , T w , q w are quantities evaluated at the wall and are the wall shear, the wall temperature, and the wall heat flux, respectively. A normalized distance to the wall of approximately y + ≈ 1 is recommended to accurately resolve boundary layer phenomena in turbulent flows.
The numerical studies presented in the paper are performed with a parallelized implementation of the XFEM approach into an in-house C++ code. The forward analysis uses a staggered solution strategy, i.e., the problem is solved successively for the heat method state variables, for the distance variables, followed by the flow and SA state variables, and finally the temperature variables. The flow-SA equations are solved staggered using a pseudo time stepping with a local CFL strategy and an increasing time step size are used to reach steady state, see Ceze and Fidkowski (2013) and Witherden et al. (2017). The fluid and SA subproblems are solved by Newton's method. All linear and linearized systems of equations are solved by MUMPS, see Amestoy et al. (2001). Gauss quadrature rules are used to integrate the set of discretized governing equations on each integration subdomains. Working in two dimensions, a 2×2 -point integration rule is used for quadrangular integration elements, a 7-point integration rule is used for triangular integration elements, and a 2-point integration rule is used for interface line elements.

Backward facing step
This subsection focuses on the validation of the XFEM formulation of the incompressible Navier-Stokes and SA models. The energy transport, in the fluid or solid phase, is not considered here. To this end, the backward facing step problem is investigated numerically using the approach proposed in Sect. 4. Following NASA Langley Research Center (2022a), the geometry of the backward facing step is illustrated in Fig. 5, where the dimension H is set to H = 0.0127 m. It should be noted that the geometry is not represented to scale. The fluid is chosen as water and is characterized by a density f = 1.18 kg/m 3 and a kinematic viscosity f = 1.57×10 −5 m 2 /s at atmospheric conditions, and at a reference temperature of 300 K. The boundary conditions are presented in Fig. 5. The inflow conditions are represented in blue. A uniform inlet velocity with u x,in = 41.5 m/s and u y,in = 0 m/s is enforced corresponding to a Reynolds number of Re = 33, 570 . A uniform inlet viscosity is prescribed with ̃i n = 3 f . Additionally, symmetry conditions, in green, are imposed close to the inlet, i.e., a zero tangent velocity u y,sym = 0 m/s and a prescribed viscosity ̃s ym = 3 f . At the walls, colored in black, a no-slip condition is imposed on the velocity, i.e., u x,wall = u y,wall = 0 m/s, and the modified viscosity is set to zero ̃s ym = 0 m 2 /s. At the outlet, marked in red, the pressure is set to zero. The properties used for the u x,wall u y,wallνwall  Table 1. The ability of the proposed method to model turbulent flow is demonstrated by comparison against experimental results from Driver and Seegmiller (1985) and numerical results from NASA Langley Research Center (2022b). This study is performed on a refined mesh around the walls to allow for an accurate resolution of the wall distance field and of the boundary layer phenomena, as shown in Fig. 6. Mesh elements have an initial size h = H∕5 and are refined seven times around the walls leading to a mesh size of h = H∕(5×2 7 ) , 2,727,478 mesh elements, and 7,781,480 DOFs. This refinement leads to a normalized distance to the wall y + ≈ 3 for the first element near the wall.
Different settings of the SA model are considered. In particular, the use of the f t2 term in Eqs. (21,22) and the use of the proposed approximated distance d in Eq. (46) instead of an exact analytical distance function are investigated. The solution obtained using the f t2 term and the exact wall distance is presented in Fig. 7. Figure 8 presents the velocity and turbulence profiles at different locations downstream of the step, i.e., x∕H = 1.0, 4.0, 6.0 and 10.0. Figure 8 shows that, regardless of the SA model configuration, the numerical results produced using the XFEM approach are in good agreement with both the numerical results from NASA Langley   These results indicate that the wall distance field is sufficiently resolved by the approximate wall distance proposed in Eq. (46). Therefore, this approximation will be used in subsequent simulations. Concerning the effect of the f t2 term, numerical results show that the term has a minor influence on the flow conditions at this Reynolds number. Small difference between velocity profiles with f t2 = 0 (solid and dashed red curves) and f t2 = 1 (solid and dashed blue curves) are only observed inside the channel. These results agree with the observations of Rumsey (2007). The influence of the f t2 term on topology optimization results is further studied in Sect. 7.

Laminar and turbulent flow over a thick flat plate
In this subsection, the coupling of the flow model with the advection-diffusion equation is investigated. The problem of a laminar or a turbulent flow over a thick plate with a prescribed temperature at its bottom side is studied numerically using the proposed XFEM approach. Following Vynnycky et al. (1998), the geometry of the computational domain is presented in Fig. 9, where the dimensions are set with H = 1 m. The boundary conditions are illustrated in Fig. 9. The inflow conditions consist of a constant inlet velocity with u x,in = 0.1 m/s and u y,in = 0 m/s and an inlet temperature results from XFEM approach with f t2 term ( f t2 = 1 ) and with exact distance, without f t2 term ( f t2 = 0 ) and with exact dis-tance, with f t2 term ( f t2 = 1 ) and with approximate distance d, without f t2 term ( f t2 = 0 ) and with approximate distance d, numerical results from NASA Langley Research Center (2022b), and experimental results from Driver and Seegmiller (1985) u

Fig. 9
Turbulent flow over a thick flat plate with a prescribed temperature at its bottom side: geometry and boundary conditions T in = 300 K. If a turbulent flow is considered, a uniform inlet viscosity ̃i n = 3 f is prescribed. Symmetry conditions are imposed along the edges in green. A zero tangent velocity u y,sym = 0 m/s and a zero viscosity ̃s ym = 0 m 2 /s are imposed. On the surface of the plate and at the bottom of the fluid domain, represented in yellow, a no-slip condition is imposed, and the modified viscosity is set to zero. The plate is heated to a constant temperature T b = 310 K at its bottom side colored in orange. The outflow conditions consist of a zero pressure at the outlet in red. Finally, the flow and heat transfer are governed by fixed adimensional numbers. For all simulations, a Reynolds number Re = 10, 000 is used. Two different Prandtl numbers, Pr = 0.01 and 100 are investigated.
The solid plate has a density s = 1.0 kg/m 3 and a conductivity s = 100.0 W/mK. The fluid has a density f = 1.0 kg/m 3 , a kinematic viscosity f = u x,in H∕ Re , and a dynamic viscosity f = f f related to the flow conditions. The thermal properties of the fluid are defined by a conductivity f = s ∕ , where is a prescribed conductivity ratio between the solid and the fluid, and a specific heat capacity c f p = f Pr ∕ f J/kgK. When considering turbulence, the SA model is used with f t2 = 0 and with the approximate wall distance. The properties used for the heat method are = 1.0 kg/m 3 , c p = 0.1 J/kgK, and = 1.0 W/mK. Finally, the penalty parameters are set as detailed in Table 2.
The ability of the proposed method to model conjugate heat transfer is demonstrated by comparison against analytical results from Vynnycky et al. (1998) and numerical results from Yau (2016). The XFEM analysis is performed on a mesh refined near the plate to allow for an accurate resolution of boundary layer phenomena. Mesh elements have an initial size h = H∕5 and are refined six times near the plate, as shown in Fig. 10. The resulting mesh has a minimum element size of h = H∕(5×2 6 ) with 143,844 quadrangular elements and 398,592 DOFs. For the turbulent case, this mesh refinement lead to a normalized distance to the wall y + = 0.87. Figure 11 presents the adimensional temperature contours in the fluid and solid domains obtained with the XFEM approach considering a laminar and a turbulent flow with a fixed Reynolds number of Re = 10, 000 and for two Prandtl numbers, Pr = 0.01 and Pr = 100 . The adimensional temperature T adim is defined as follows: For a lower Prandtl number, a low ratio allows for a larger temperature drop within the solid. For a higher Prandtl number, the temperature drop is limited to the solid domain. The temperature profiles observed for the laminar and the turbulent flow are similar suggesting that turbulence has little effect on the thermal response and the problem is convection dominated.
The evolution of the adimensional temperature in terms of the distance x along the plate is shown in Fig. 12. The numerical results obtained with the XFEM approach are compared against studies from Vynnycky et al. (1998) and Yau (2016) for a laminar flow over the plate. The obtained temperature profiles are in good agreement with both analytical results from Vynnycky et al. (1998) and numerical results from Yau (2016). This comparison verifies the ability of the proposed XFEM approach to accurately model conjugate heat transfer.

Optimization problems
In this paper, we consider optimization problems that can be formulated as follows: where the vector contains the N s design variables and s and s are the imposed lower and upper bounds, respectively. The vector ( ) collects all state variables such that ( ) = [ p̃T D P ] T . Note that the state variables used in Eq. (65) satisfy the discretized governing equations for a given design . The objective function is denoted by z( , ( )) and g j ( , ( )) is the j th constraint function considered in the problem.  The objective function is defined as: (66) z( , ( )) = w f F( , ( )) + w p P p ( ) + w r P r ( ), where F( , ( )) is the quantity of interest to minimize; P p ( ) and P r ( ) are the perimeter penalty and the regularization penalty, as further defined in Subsect. 6.1. The parameters  Fig. 12 Turbulent flow over a flat plate with a prescribed temperature at its bottom side: adimensional temperature profile over the plate for conductivity ratios = 1, 5, 20 considering laminar and turbulent results from the XFEM approach and laminar analytical results from Vynnycky et al. (1998) and numerical results from Yau (2016) w f , w p , and w r are weights associated with the terms in the objective function. Following Haber et al. (1996), the perimeter control penalization P p is added to the objective function to ensure the well-posedness of the optimization problem and to prevent the emergence of numerical artefacts, such as irregular geometric features:

Regularization
where P 0 is the initial perimeter value.
To avoid spurious oscillations of the level set field that might deteriorate the stability and the convergence of the optimization problem, the regularization scheme proposed in Geiss et al. (2019) is adopted. The regularization promotes the convergence of the level set field to a smooth target field ̃ with advantageous properties, i.e., the convergence to an upper or lower bound B away from the interface and to a uniform gradient along the interface. The scheme is enforced by adding the penalty term P r to the objective function: The weights w and w ∇ are evaluated as: where the parameter P r controls the regions of influence of the regularization, i.e., regions near and away from the interface. The weights w 1 , w 2 control the mismatch between design and target LSF near and away from the interfaces, respectively. The weights w ∇ 1 and w ∇ 2 play the same role for the mismatch between design and target gradients of the LSF. The smooth target level set field ̃ is built from the distance field D , obtained by the heat method in Eq. (44), as: where B is a parameter that is set to 1 or −1 to select the upper bound B for the fluid domain or lower bound − B for the solid domain respectively.

Sensitivity analysis
In this paper, the update of the design variables is performed by mathematical programming techniques driven by shape sensitivity computations. The derivative of a generic objective or constraint function F( , ( )) with respect to a design variable s i is computed as: where the first and the second term in the right-hand side account for the explicit and the implicit dependency on the design variables, respectively. Using the adjoint approach, see Michaleris et al. (1994), the implicit term is evaluated as: where is a collection of the residuals presented in Sect. 4 for the forward analysis and represents the adjoint responses, evaluated through: The semi-analytical approach of Sharma et al. (2017) is used to evaluate the derivatives with respect to the design variables. During the sensitivity analysis, all the terms in the governing equations introduced in the forward analysis in Sect. 4 are accounted for, and all required derivatives with respect to the state variables are derived analytically following the approach of Zymaris et al. (2009) andYoon (2016) to yield consistent sensitivities.

Numerical optimization examples
In this section, the features of the proposed XFEM level set-based optimization framework for solving turbulent conjugate heat transfer design problems are investigated. Two different optimization problem formulations are considered aiming at: (i) minimizing the maximum temperature, and (ii) maximizing the heat transfer within the heat exchangers, while maintaining an acceptable pressure drop. Additionally, two different heat exchanger geometric configurations are studied to demonstrate the flexibility of the approach.
Unless specified otherwise, all subsequent simulations are carried out using the following settings. By default, the f t2 term is omitted in the SA model. For specific configurations the analysis and optimization results are compared when using or omitting the f t2 term. The approximate distance field in Eq. (46) is used for the wall distance and the properties for the heat method are chosen as = 1.0 kg/m 3 , c p = 0.1 J/kgK, and = 1.0 W/mK. The penalty parameters introduced in Sect. 4 are detailed in Table 3. In the following examples, the state variable fields are approximated with linear B-spline functions, while the design variable field is approximated with quadratic B-spline functions, as they promote smoother designs, mitigate the development of spurious features, and limit the need for additional filtering techniques, see Noël et al. (2020). Filtering is achieved through a coarser and higher-order B-spline interpolation for the design variable field.
The set of discretized governing equations and adjoint sensitivity equations are built and solved following the approach described in Sect. 5, except that for the sensitivity analysis, the flow and the SA adjoints are solved for simultaneously using a monolithic formulation. The optimization problems are solved by the Globally Convergent Method of Moving Asymptotes (GCMMA) of Svanberg (2002). The initial, lower, and upper asymptote adaptation parameters in GCMMA are set to 0.5, 0.7, and 1.2, respectively. No inner iterations are used. The optimization problems are considered converged if the design is visually converged, if the objective function stagnates, and if the constraints are satisfied. The weights and parameters introduced in Sect. 6 are given in Table 4.
While adaptive mesh refinement is performed to resolve the boundary layer for the benchmark problems presented in Sect. 5, uniformly refined meshes are used for the optimization problems. This choice is made as the Reynolds numbers considered in this section are lower, and therefore sufficiently fine, uniform meshes can be used. This is verified in Subsect. 7.2 by solving the same optimization problem on a reference and a refined mesh.

Heat exchanger with single outlet
In this subsection, a heat exchanger with a single outlet is optimized for minimal average temperature in the solid subject to a constraint on the power dissipation. The problem is taken from Dilgen et al. (2018b) and the setup is presented in Fig. 13. The fluid is chosen as water and is used to cool a heated solid, here aluminum. The dimensions of the heat exchanger are described in terms of half the inlet size H = 0.1 m. The inlet and outlet channels are extended, see red and blue dashed lines, to reduce the influence of the inlet and outlet conditions on the flow inside the heat exchanger.
The working fluid is water with a density f = 1.18 kg/m 3 and a dynamic viscosity f = 1.77×10 −5 kg/ms. The fluid has a specific heat capacity c f p = 1004.9 J/kgK and a conductivity f = 2.54×10 −2 W/mK. The turbulent Prandtl number is set to Pr t = 0.9 . The solid phase is aluminum with a thermal conductivity s = 237 W/mK.
At the inlet Γ in , boundary conditions are prescribed for the fluid. To obtain a developed flow at the inlet, a power law velocity profile is enforced, following Kubo et al. (2021), such that: The magnitude u in is chosen to yield specific Reynolds numbers Re = f u in H∕ f . A constant viscosity is prescribed at the inlet with a magnitude ̃i n = 3 f . The fluid enters the heat exchanger at a temperature T in = 300 K. At the outlet Γ out , the pressure is prescribed to zero p out = 0 Pa. The solid is heated by a constant volumetric heat load Q = 10.0 kW/m 2 .
To minimize the solid temperature, the objective function is formulated as: u x, in = u in 1 − y H 1 n , n = 2 log 10 Re 10.0 ,

Fig. 13
Heat exchanger with a single outlet: geometry and boundary conditions where T s is the temperature integrated over the solid domain Ω s : The reference temperature T ref is set to the inlet temperature, i.e., T ref = 300 K. The volume V s of the solid domain Ω s is evaluated as: The parameters T s 0 and V s 0 are the integrals of the temperature over the solid domain and the volume of the solid evaluated for the initial design.
The power dissipation constraint is formulated as: where and D in∕out,0 is the dissipated power integrated over the inlet or the outlet for the initial design. The weight w d controls the allowed power dissipation through the heat exchanger. Finally, a volume constraint is enforced on the fluid domain to avoid a trivial solution with no solid and to control the amount of energy injected into the heat exchanger through the solid. The constraint is formulated as: where V f is the volume of the fluid domain Ω f evaluated as in Eq. (76), V is the volume of the design domain, and w v is the allowed fluid volume faction, here set to w v = 0.55.
The initial design is seeded with multiple solid inclusions, as shown in Fig. 14. The level set field is approximated by quadratic B-spline and is discretized on a coarse uniform background mesh with an element size h = H∕5 , leading to 4284 design variables. The state variable fields are discretized on a finer uniform background mesh with an element size h = H∕40 , leading to 251,264 quadrangular elements and 896,292 DOFs including 414,924 flow DOFs, 125,879 temperature DOFs, 251,758 DOFs for the heat method, and 103,731 DOFs for the construction of the approximate distance. This mesh size leads to normalized distances to the wall y + ≈ 2, 5, and 10 for Reynolds numbers Re = 1000, 2500, and 5000, respectively, which is insufficient to fully resolve boundary layer phenomena, but offers a reasonable compromise between accuracy and computational cost. The optimization problem is solved considering three different inlet velocities u in = 0.15 , 0.375, and 0.75 m/s leading to the following Reynolds numbers, Re = 1000 , 2500, and 5000. The power dissipation is limited to a third of the dissipation observed for the initial design, i.e., the weight in Eq. (77) is set to w d = 0.33 . It should be noted that, as the inlet velocity changes, so does the initial power dissipation, and thus the allowed power dissipation in the heat exchanger. Therefore, comparison between the optimized designs and their performances should be made carefully.
The optimized designs and the corresponding level set fields are shown in Fig. 15. For all Reynolds numbers, the optimized design features a network of channels. Owing to the regularization, the level set fields converge to the upper and lower bounds in the bulk of the fluid and solid domains and have close to uniform spatial gradient along the fluid/ solid interface.
The solution fields over the optimized designs are presented in Fig. 16. For different Reynolds numbers, the velocity magnitude, the viscosity ratio, and the temperature over the optimized heat exchanger, as well as the performance of the designs in terms of average temperature in the solid domain T s ∕V s , are provided. Focusing on the geometry of the designs, a tendency to generate smooth straight channels that align with the flow is observed for all considered Reynolds numbers. For lower Reynolds numbers and thus lower power dissipation, smooth wide channels are created. As the Reynolds number is increased and a larger power dissipation is allowed, narrower channels with higher flow velocities are formed. These observations agree with the results presented in Dilgen et al. (2018b).
As expected, the heat transfer is enhanced for higher Reynolds numbers, allowing a larger dissipation in the heat exchanger. This is evident from the decrease of the maximum temperature in the solid phase, i.e., the adimensional maximum temperature reaches T max ∕T in = 2.4755, 1.5528, The contour plots of the viscosity ratio in Fig. 16 show that the optimized designs exhibit low or no turbulence in the channels within the heat exchanger. Interestingly, the maximum viscosity ratio decreases as the Reynolds number increases: = 5.9648 , 4.3823, and 3.7297 for Reynolds numbers Re = 1000 , 2500, and 5000, respectively. As the Reynolds number and the allowed power dissipation are increased, the optimized designs feature narrower channels which impede the development of turbulence.
To separate the influences of increasing the Reynolds number and the allowed power dissipation on the generated design, the optimization problem is solved considering three different constraint weights w d = 0.25 , 0.33, and 0.50 at a fixed inlet velocity u in = 0.375 m/s leading to a Reynolds number Re = 2500 . Figure 17 shows the velocity magnitude, the viscosity ratio, and the temperature over the optimized heat exchanger considering different allowed power dissipation. The performance in terms of average temperature in the solid domain T s ∕V s is provided for each design.
As the allowed pressure drop is increased, a larger number of narrower channels through the heat exchanger is created. A larger pressure drop across the heat exchanger allows for the development of narrow channels with higher flow velocities, and advective energy transport, which in turns enhances the heat transfer and lowers the adimensional maximum temperature: T max ∕T in = 1.5985, 1.5528, and 1.4009 for constraint weights w d = 0.25, 0.33, and 0.50, respectively. It should be noted that the distribution of the turbulence and the maximum viscosity ratios within the heat exchanger = 4.1064, 4.3823, and 3.7496 are similar for all considered pressure drops.
The influence of the f t2 term in the SA model is investigated by solving the optimization problem with and without the f t2 term, i.e., setting f t2 to 1 or 0. The velocity magnitude, the viscosity ratio, and the temperature distributions over the optimized heat exchangers, as well as the performance in terms of the average temperature in the solid domain T s ∕V s , are presented in Fig. 18 for Re = 1000 and w d = 0.50 and in Fig. 19 for Re = 5000 and w d = 0.33 . Reanalyses are performed on the designs, i.e., designs optimized with f t2 = 0 are reanalyzed setting f t2 = 0 and f t2 = 1 , see Fig. 18a, and designs optimized with f t2 = 1 are reanalyzed setting f t2 = 0 and f t2 = 1 , see Fig. 18b.
For a lower Reynolds number in Fig. 18, the designs generated with and without the f t2 term differ only slightly. Performing reanalyses shows that marginally lower mean temperatures over the heat exchanger are obtained when dropping the f t2 term, i.e., T s ∕V s = 1.6774 K for f t2 = 0 versus 1.6779 K for f t2 = 1 when optimizing designs with f t2 = 0 , and T s ∕V s = 1.6880 K for f t2 = 0 versus 1.6886 K for f t2 = 1 for design obtained with f t2 = 1 . Note that when predicting the thermal response with f t2 = 1 , the design optimized without the f t2 term has a slightly lower mean temperature than the design optimized with the f t2 term. This issue is likely explained by the non-convexity of the problems solved.
For a higher Reynolds number in Fig. 19, similar observations can be made and reanalyses lead to T s ∕V s = 1.2063 K for f t2 = 0 versus 1.2064 K for f t2 = 1 when generating designs with f t2 = 0 , and T s ∕V s = 1.2022 K for f t2 = 0 versus 1.2024 K for f t2 = 1 when generating designs with f t2 = 1 . It should be noted that the effect of the f t2 term is less pronounced for Re = 5000 than for the lower Reynolds number case. In conclusion, while the design layout marginally changes, performances across designs optimized with and without the f t2 term differ insignificantly. In agreement with Rumsey (2007), the effect of including the f t2 term becomes negligible as the Reynolds number increases.

Heat exchanger with split outlet
In this subsection, the versatility of the proposed optimization approach is illustrated by varying the design domain, the inflow conditions, and the heat loading. The influence of the formulation of the optimization problem is also investigated. A heat exchanger with a split outlet is optimized for maximum heat transfer under a mass flow constraint, see Fig. 20. Contrary to the single outlet case in Sect. 7.1, the geometry of the design domain prevents the development of a straight channel between the inlet and the outlet. In contrast to the previous example, a constant pressure is imposed at the inlet fixing the allowed pressure drop through the heat exchanger. In this example, the solid is heated up to a prescribed temperature and the design goal is to maximize the heat transfer between fluid and solid. The effects of these variations on the optimization process and the geometry of the optimized designs, i.e., the arrangement of the cooling channels within the heat exchanger, are investigated.
The problem setup is presented in Fig. 20. Water is used to cool heated aluminum. The dimensions of the design domain are described in terms of half of the inlet size H = 0.1 m. The inlet and outlet channels are again extended, see red and blue dashed lines, to mitigate the influence of inlet and outlet conditions on the flow in the heat exchanger. The fluid and solid properties are identical to the ones described for the single outlet heat exchanger in Sect. 7.1.
At the inlet Γ in , the following boundary conditions are prescribed on the fluid flow. A constant inlet pressure p in = fũ2 in is enforced and is defined as a function of a nominal inlet velocity ũ in that corresponds to a desired Reynolds number Re = fũ in H∕ f . The inlet pressure increases with the Reynolds number as the dynamic viscosity remains constant. A constant viscosity is prescribed at the inlet with a magnitude ̃i n = 3 f . The fluid enters the heat exchanger at a temperature T in = 300 K. At the outlet Γ out , the pressure is prescribed to zero p out = 0 Pa. The solid is heated up to a temperature of T ref = 325 K through a temperature dependent volumetric heat load To maximize the heat transfer through the exchanger, the objective function is defined as: where the convective flux Q in∕out at the inlet or the outlet is evaluated as: and Q in∕out,0 is the convective flux evaluated at the inlet or the outlet for the initial design. The mass flux constraint ensures a non-zero flow from the inlet to the outlet and that the solid does not obstruct the heat exchanger. The constraint is formulated as: where the mass flux M in∕out through the inlet or the outlet is evaluated as: and M in∕out,0 is the mass flux computed at the inlet or the outlet for the initial design. The weight w m specifies the allowed mass flux through the heat exchanger and is set to w m = 0.5. The initial design is seeded with multiple solid inclusions, as shown in Fig. 21. The level set field is discretized with quadratic B-spline on a coarse uniform background mesh with an element size h = H∕10 , leading to 16,059 design variables. The state variable fields are discretized on a finer mesh with an element size h = H∕40 , leading to 245,616 quadrangular elements and 875,477 DOFs including 405,640 DOFs for the turbulent flow, 122,809 DOFs for the temperature, 245,618 DOFs for the heat method, and 101,410 DOFs to build the distance field. This mesh size leads to normalized distances to the wall y + ≈ 2 , 5, and 10 for Reynolds numbers Re = 1000 , 2500, and 5000, respectively, which is again insufficient to fully resolve boundary layer phenomena (83) M in∕out = ∫ Γ in∕out ⋅ dΓ, but offers a reasonable compromise between accuracy and computational cost, as further discussed later. The influence of the inlet flow conditions on the optimized designs and their performances is investigated by solving the problem for three different inlet pressures corresponding to Reynolds numbers Re = 1000 , 2500, and 5000. The optimized designs are shown in Fig. 22. For each design, the velocity magnitude, the viscosity ratio, the temperature, and the performance in terms of heat flux ΔQ = Q in − Q out are provided. Contrary to the single outlet case in Sect. 7.1, tortuous channels around multiple solid inclusions and recirculation areas are generated for all Reynolds numbers considered here. Such geometric features increase the area between fluid and solid, and in turn the heat transfer. With the increase in the Reynolds number, the large central channel tends to split into multiple narrow and tortuous channels.
As the Reynolds number increases, the turbulence distribution changes and areas with larger viscosity ratio move from the inlet to the core of the heat exchanger. Note that with increasing Reynolds number, the surface of the channels tends to be less smooth and presents a secondary oscillatory Adimensional temperature T /T in pattern, especially visible for the design at a Reynolds number Re = 5000 . In the following, these patterns are investigated to understand whether they are a numerical artifact or present a physical advantage. To this end, the maximum heat transfer optimization problem with Re = 5000 is solved on a refined mesh with an element size of h = H∕80 and a normalized distance to the wall of approximately y + ≈ 5 . Note that the level set field is represented with the same spatial discretization as used previously, i.e., a coarser mesh with an element size h = H∕10 . The velocity magnitude, the viscosity ratio, and the temperature fields on the optimized designs, as well as the performance in terms of heat flux ΔQ , for the reference and the refined mesh are shown in Fig. 23. Adimensional temperature T /T in Optimized design for β f t2 = 1

Fig. 19
Heat exchanger with a single outlet: optimized designs and performance in terms of average temperature in the solid domain T s ∕V s for a Reynolds number Re = 5000, at a fixed power dissipation with w d = 0.33 While their geometry differs in their detailed layout, both designs present similar turbulence generating features, i.e., multiple narrow, tortuous channels and secondary patterns along the channels. The designs performance differs by less than one percent. The design optimized on the reference mesh achieves a slightly higher heat flux with ΔQ = 25, 959.9 W/m 2 compared to the design optimized on the fine mesh with ΔQ = 25, 029.3 W/m 2 . The close-ups provided in Fig. 23 confirm that the finer mesh allows for an improved resolution of the boundary layer. The velocity profile and the temperature gradient close to the solid/ fluid interface are represented across several elements, which yields a better accuracy. Nonetheless, the reference mesh provides a good compromise between accuracy and computational cost.
To gain insight into the reasons behind the development of oscillatory channel walls, the optimization problem at a fixed Reynolds number Re = 5000 is solved considering three different settings of the SA model namely, the presence of the f t2 term and the turbulent fluid conductivity t . Configuration (a) represents the setting used previously, i.e., the f t2 term is dropped and the turbulent fluid conductivity t is accounted for. In Configuration (b), the f t2 term is accounted for, i.e., f t2 = 1 , as well as the turbulent fluid conductivity t . Finally, for Configuration (c), the f t2 term is included, but the turbulent fluid conductivity is set to zero, i.e., f t2 = 1 and t = 0.
The velocity magnitude, the viscosity ratio, and the temperature distributions over the optimized heat exchangers are presented in Fig. 24. For each configuration of the SA model, reanalyses are performed on each optimized design. The design generated with setting (a) is analyzed using settings (a), (b), and (c), see Fig. 24a. A similar procedure is followed for designs obtained with settings (b) and (c) as shown in Fig. 24b and c, respectively. For each analysis, the value of the heat flux ΔQ = Q in − Q out is provided.
From the results in Fig. 24a and b, similar conclusions to the study on the f t2 term presented in Subsect. 7.1 can be drawn. While the design geometry slightly changes, the performances of the designs generated including or neglecting the f t2 terms only differ marginally. However, omitting the turbulent diffusivity t leads to noticeable drop in heat flux. Figure 24c shows the design generated when omitting the turbulent diffusivity t . In this case, the design only features smooth channel surfaces, but the heat flux is significantly lower than the ones achieved when considering the turbulent conductivity t as a result of increased turbulence. This study suggests that the development of oscillatory shapes along the channels is used to generate an increased turbulence within the heat exchanger, which in turn increases the total thermal conductivity in the fluid and leads to an enhanced heat transfer. Additional studies using different turbulence models are needed to confirm these results.

Conclusions
In this paper, we proposed an XFEM level set-based topology optimization framework for turbulent conjugate heat transfer design problems. To achieve a crisp description of the fluid/solid interface, the geometry of the design is represented by one or multiple LSFs and the physical responses of the design are predicted using an IFEM, here the XFEM approach. The fluid flow is modeled by the RANS equations and the system is closed by the SA turbulence model. As the turbulence equation depends on the distance to the closest wall, a distance field is constructed using the heat method. The heat transfer in the fluid phase is described by the advection-diffusion equation and in the solid by linear isotropic diffusion. Boundary and interface conditions are enforced weakly with Nitsche's method, and the face-oriented ghost stabilization is used to mitigate numerical instabilities associated with small integration subdomains. Verification cases are solved to show that the XFEM approach can reproduce benchmark results using an immersed geometry description. The ability of the proposed optimization framework to solve turbulent conjugate heat transfer problems is illustrated with the design of heat exchangers in two dimensions.
Numerical simulations show that the proposed XFEM approach is able to accurately predict the physics responses in turbulent flows. Using an immersed geometry description and an IFEM, benchmark results for the backward facing step and the turbulent flow over a thick conducting plate are reproduced. Additionally, the influence of the setting of the SA model on the numerical results is investigated and the numerical results with and without specific terms, such as the f t2 term, are in agreement with literature.
Numerical optimization examples demonstrate the ability of the proposed level set-based optimization framework to solve conjugate heat transfer design problems considering turbulent flows. The design of heat exchangers is performed under different flow conditions and for different geometries of the design domains. The optimized layouts achieve enhanced cooling or heat transfer under prescribed conditions. Furthermore, the versatility of the framework is illustrated by solving different optimization problems. The numerical examples show that the formulation of the optimization problem has a large influence on the optimized geometries. The optimization on the split outlet heat exchanger suggests that under particular flow conditions heat transfer can be enhanced through channel wall shapes that generate turbulence, and in turns increase the total conductivity of the fluid.
This work opens up several follow-up research tracks. As shown in the benchmark examples, an increased resolution of the fluid/solid interface is desirable, i.e., y + around 1, to achieve a more precise description of the interface and a higher accuracy of the physics predictions. To this end, different strategies could be investigated: the use of hierarchical B-splines to carry out local mesh refinement, the use of wall functions for an improved description of the near-wall behavior, or the ability to generate a boundary layer mesh with immersed methods. The second optimization problem suggests that, using the SA model, oscillatory channel walls can be used to generate turbulence and increase the heat transfer. Additional studies relying on different turbulence models should be carried out to confirm this result and to understand the influence of turbulence models on the optimized topology and geometry. The extension of the framework to three dimensional optimization problems also presents various challenges. An overall increased efficiency of the computational strategy, including the use of local mesh refinement and iterative solvers, will be necessary to meet the mesh resolution requirements and to handle the sizes of the systems of discretized equations. Finally, including additional physics phenomena in the model, such as fluid/structure interaction or chemical reaction, would allow us to tackle a broader range of applications. otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.