Topology Optimization of Cooling Channels with Conjugate Heat Transfer under Non-Uniform Heat Sources

1  Introduction

Driven by the new industrial revolution and the global energy transition, thermal management has evolved from a traditional supporting function into a core technology for ensuring the safe and efficient operation of advanced systems [1–5]. Concurrently, rapid advancements are driving electronic components and mechanical structures toward miniaturization, integration, and higher power density, which significantly increases thermal loads. For instance, the localized heat flux in advanced integrated circuits has risen sharply, making traditional air-cooling techniques increasingly inadequate [6,7]. In hypersonic vehicles, extreme aerodynamic heating can cause localized temperature rises which may lead to material degradation or structural failure if uncontrolled [8,9]. Similarly, the risk of thermal runaway in new energy vehicle power batteries remains a critical challenge for both safety and driving range [10,11].

Efficient thermal management is therefore critical. Convective heat transfer, as a fundamental mode of heat transfer, dominates industrial energy transport processes owing to its high energy transport capacity [12–15]. This mechanism plays a fundamental role in the design of cooling channels. In practice, heat is typically removed via forced convection by fabricating internal channels within a heated structure and circulating a coolant. However, improving heat exchange performance often comes at the cost of a significantly increased pressure drop. Consequently, balancing heat transfer efficiency, structural weight, pumping power, and temperature uniformity represents a complex multiphysics optimization problem in engineering design [16–18].

Traditional cooling channel designs, including straight channels, serpentine passages, and pin-fin arrays, have been extensively investigated and widely adopted in engineering practice. However, these conventional configurations exhibit several inherent limitations, particularly under non-uniform heat source conditions. First, the design space is constrained by predefined geometric configurations, permitting only parametric optimization of channel size or shape without altering the overall topology [19–21]. This limitation reduces the ability to achieve globally optimal thermal-fluid performance. Second, conventional microchannel designs often suffer from localized hot spots under non-uniform heat flux distributions due to their limited adaptability to spatially varying thermal loads [22,23]. Finally, these designs are typically developed through empirical correlations and iterative experimental testing, which are time-consuming and do not necessarily yield optimal solutions [18,24]. These limitations highlight the need for more flexible and systematic design approaches capable of automatically generating complex channel networks tailored to specific thermal conditions.

To overcome these limitations, topology optimization provides a powerful design alternative for cooling channels [25–31]. This approach generates an optimal distribution of material within a prescribed design domain, with objectives such as maximizing heat transfer or minimizing the mean temperature under constraints like volume fraction or permissible pressure drop. Since its development via the homogenization method [32–34], various numerical approaches, including density-based methods [35–41], the level-set method [42–47], and evolutionary structural optimization [48–51], have been developed. Topology optimization has been notably successful in structural lightweight design and is increasingly applied to fluid-thermal coupled systems.

In fluid topology optimization, Borrvall and Petersson [52] established the foundational framework by introducing porous medium theory into fluidic design. They achieved a unified mathematical representation of fluid and solid materials across the design domain by incorporating a Brinkman penalty term into the Stokes equations. This approach formed the standard framework for subsequent density-based fluid topology optimization. Later, Gersborg-Hansen et al. [53] and Olesen et al. [54] extended the governing equations to the steady-state Navier-Stokes equations, incorporating convective inertial terms. Their works enabled the application of topology optimization to flow problems at moderate Reynolds numbers. However, at higher Reynolds numbers, the nonlinearity of the flow field becomes significantly enhanced, and numerical instabilities become a primary challenge. To address this, subsequent research has focused on developing numerical stabilization techniques and robust solution algorithms [55–59].

While fluid topology optimization primarily seeks to minimize flow resistance, effective thermal management requires enhancing heat transfer. Therefore, integrating the energy equation to construct a fluid-solid coupled conjugate heat transfer (CHT) topology optimization framework is a logical necessity [60–68]. Yoon [69] and Dede [70] performed pioneering work on the topology optimization of forced convection cooling structures. By employing thermal compliance or maximum temperature as the objective function, they obtained innovative heat sink designs under multiphysics coupling conditions. Furthermore, Koga et al. [71] utilized the ability of the level set method to capture sharp fluid-solid interfaces, successfully designing efficient active cooling channels. Their work demonstrated the considerable potential of topology optimization for thermal management applications.

Despite considerable progress, existing CHT topology optimization still faces significant challenges, primarily due to the conflict between accurately modeling turbulent flow and managing computational cost. In practical engineering applications, coolant flows are typically turbulent at high Reynolds numbers, and the laminar flow assumption leads to a marked degradation in predictive accuracy [72]. Although Dilgen et al. [39,73] incorporated turbulence models into topology optimization, the strong nonlinearity and high sensitivity of these models lead to complex sensitivity analysis and poor convergence. More critically, full three-dimensional (3D) turbulent CHT simulations entail a high number of degrees of freedom, rendering even a single simulation iteration computationally expensive [74–77]. Consequently, driving hundreds or thousands of optimization iterations directly with such high-fidelity models remains computationally prohibitive for practical engineering purposes.

To address these challenges from a modeling perspective, existing CHT topology optimization approaches can be broadly classified into full-order models (FOMs) and reduced-order models (ROMs). FOMs are capable of accurately resolving detailed flow and heat transfer physics, and transient full-order simulations can further capture time-dependent thermal responses under realistic operating conditions. However, these approaches are computationally intractable for topology optimization, as each design iteration requires solving large-scale nonlinear systems with millions of degrees of freedom. This issue becomes even more pronounced when considering transient processes and complex geometries, as highlighted in recent studies [66,78]. To alleviate this limitation, ROM strategies have been increasingly adopted. These approaches significantly reduce computational cost by simplifying the governing equations or reducing the dimensionality of the problem. However, this efficiency gain often comes at the expense of physical fidelity, particularly in accurately capturing turbulent flow behavior and thermal effects [79].

Therefore, pseudo-three-dimensional (pseudo-3D) models have been developed to address this computational challenge [80–83]. These models reduce the problem from three to two dimensions by incorporating thickness-direction integral averaging and equivalent source terms into the two-dimensional (2D) governing equations, thereby significantly improving computational efficiency while retaining key 3D effects. Conventional pseudo-3D models typically assume laminar or Darcy flow. Their successful engineering application thus hinges on calibrating their parameters so that the model’s macroscopic behavior accurately represents that of actual turbulent flows.

In practical engineering, heat flux distributions are often non-uniform. For example, localized hot spots on aerospace vehicle skins can experience significantly higher heat fluxes than surrounding areas. This non-uniformity presents a challenge for cooling channel topology optimization [22,84,85]. The optimization must automatically adjust to the heat source distribution, generating denser channel networks in high-flux regions to enhance localized cooling. Simultaneously, under a limited total flow rate, coolant must be allocated efficiently between high- and low-heat-flux regions to achieve temperature uniformity. However, while dense branching channels improve heat dissipation, they also increase flow resistance. Therefore, the optimization must effectively eliminate hot spots while maintaining the total system pressure drop within allowable engineering limits.

To address the aforementioned challenges, this study develops a practical fluid-solid CHT topology optimization method for cooling channel design under non-uniform heat source conditions. We first construct a pseudo-3D CHT model based on Darcy’s law and then calibrate its permeability parameters using high-fidelity Reynolds-averaged Navier-Stokes (RANS) simulations to match its macroscopic behavior to real turbulent flow. Multi-condition topology optimization is performed for various inlet-outlet configurations under non-uniform heat sources. Finally, the optimized 2D topologies are reconstructed into 3D solid models and validated through numerical simulations to assess their performance. This systematic evaluation demonstrates the advantages of the optimized configurations over traditional straight channels in reducing maximum temperature, controlling pressure drop, and improving temperature uniformity.

The remainder of this paper is organized as follows. Section 2 presents the governing equations and the physical model. Section 3 describes the calibration of the permeability parameters. In Section 4, the optimization results under multi-condition are presented and analyzed. The optimized designs are validated through 3D simulations with performance comparisons. Finally, Section 5 summarizes the key findings and discusses future work.

2  Theoretical Analysis and Optimization Formulation

2.1 Physical Problem Description

This work addresses a topology optimization problem associated with an aircraft fuel tank’s active cooling system. We developed a hierarchical physical model to accurately capture the conjugate heat transfer process while maintaining computational efficiency. As illustrated in Fig. 1, the model cross-section features a typical sandwich-type composite layout.

Figure 1: Schematic of the sandwich-type composite structure in the aircraft fuel tank cooling system.

The outer layer is composed of an advanced thermal insulation material to withstand high external temperatures from the hot outside airflow encountered during flight. Beneath this lies a cast aluminum skin, which provides structural integrity and high thermal conductivity. The attenuated heat flux then conducts through the aluminum skin to reach the inner cooling channel surface. The innermost layer is the fuel distribution layer, which is the design domain for topology optimization in this study (denoted as Ωd). In this layer, the coolant (the fuel) and the solid channel walls form interpenetrating phases. A pressurized flow drives the fuel through the embedded channel network, absorbing the heat conducted from the outer skin via convection at the fluid-solid interfaces, thereby regulating the temperature and maintaining the overall thermal balance of the structure.

Solving 3D CHT problems with topology optimization presents a fundamental trade-off between accuracy and computational cost. To significantly reduce the computational cost while preserving solution accuracy, a “pseudo-3D” modeling strategy is employed. This approach comprises two key aspects. First, the flow field is simplified to two dimensions (2D) by assuming the coolant flow is predominantly in-plane, neglecting the through-thickness velocity component. This allows the flow equations to be solved in a 2D domain. Second, the thermal field is treated in a pseudo-3D manner. Although the geometry is reduced to 2D, heat transfer in the thickness direction is retained. The external heat flux passes through the insulation and aluminum skin before entering the channel layer. By incorporating equivalent source and conduction terms into the energy equation, the model captures this through-thickness heat flow. This maintains critical 3D thermal effects, such as wall-to-fluid conduction, rather than simplifying the thermal load to a purely in-plane 2D source.

Following these assumptions, the design domain Ωd is defined as a 2D region comprising two design-variable-dependent subdomains: the fluid domain Ωf and the solid domain Ωs. Its boundary consists of the inlet Γin, the outlet Γout, and the solid walls Γwall, as shown schematically in Fig. 2.

Figure 2: Schematic of the 2D design domain (Ωd), its subdomains (fluid Ωf and solid Ωs), and boundaries (inlet Γin, outlet Γout, and walls Γwall).

2.2 Fluid Flow Governing Equations

Conventional fluid topology optimization typically employs the Navier-Stokes (N-S) equations. For steady, incompressible, Newtonian flow, the momentum conservation equation is given by:

ρ(u⋅∇)u=−∇p+μ∇2u+f(1)

where ρ is the fluid density, u the velocity field, p the pressure field, μ the dynamic viscosity, and f an external body force. The convective term (u⋅∇)u introduces nonlinearity. However, in topology optimization, the flow channel layout undergoes significant topological changes, including channel disconnection, merging, and formation. Concurrently, the velocity approaches zero in solid or low-permeability regions, often leading to numerical singularities and spurious oscillations.

To address these numerical challenges, the Darcy flow model from porous media theory is adopted to approximate the macroscopic flow behavior. This model provides a numerically stable and simplified framework for representing fluid flow in topology optimization. It establishes a linear relationship where the fluid velocity within the porous medium is proportional to the pressure gradient and inversely proportional to the dynamic viscosity:

uf=−κμ∇p(2)

where uf is the fluid velocity and κ is the permeability of the medium. While permeability is generally a second-order tensor, it reduces to a scalar under the assumption of isotropy.

It is important to emphasize that the Darcy model is not intended to directly resolve the underlying turbulent flow physics in high-Reynolds-number regimes. Instead, it is employed here as an effective macroscopic resistance model that establishes a simplified relationship between pressure gradient and flow rate. In this context, the permeability parameter does not necessarily correspond to a physical porous medium, but rather serves as an equivalent parameter to represent the overall flow resistance. The nonlinear inertial effects inherent in turbulent flow are implicitly incorporated into this parameter through calibration against high-fidelity simulations, as detailed in Section 3. In this way, the proposed approach avoids directly solving the turbulent governing equations while still capturing their macroscopic influence. It bridges the gap between the linear Darcy model and high-Reynolds-number turbulent flows.

For topology optimization, the permeability κ is the key parameter to implicitly distinguish between fluid and solid regions. In regions defined as flow channels, κ is assigned a very large value (κ→∞), which greatly reduces flow resistance and approximates unobstructed fluid motion. Conversely, in solid regions, κ is set to approach zero (κ→0), rendering the flow resistance effectively infinite and thereby forcing the fluid velocity uf to near zero. This method uses a spatially varying permeability field κ(x) to implicitly define the fluid-solid interface, where x denotes the position vector. It avoids the computationally expensive and algorithmically complex explicit interface tracking required in traditional moving boundary formulations.

For an incompressible fluid, mass conservation is expressed by the continuity equation:

∇⋅uf=0(3)

Substituting the Darcy model from Eq. (2) into Eq. (3) yields a governing equation for the pressure field:

∇⋅(−κμ∇p)=0(4)

Eq. (4) is a second-order linear elliptic partial differential equation. Compared to the N-S equations, which include nonlinear convective terms, this formulation possesses significantly more favorable numerical properties. First, for any κ>0, the equation stays positive definite, guaranteeing stable convergence even with distorted meshes or poor initial guesses. Moreover, upon discretization, it yields a symmetric positive-definite linear system that can be solved efficiently using standard iterative or direct solvers.

Well-posed boundary conditions are imposed based on the physical conditions. At the inlet, a normal velocity (Neumann) condition is specified:

n⋅uf=uin,x∈Γin(5)

where uin is the prescribed average inlet velocity.

A zero-pressure (Dirichlet) condition is set at the outlet:

p(x)=0,x∈Γout(6)

This defines the pressure field relative to the outlet, simplifying the subsequent calculation of pressure drop.

All remaining boundaries are treated as impermeable walls. The no-penetration condition is enforced by setting the normal velocity to zero, corresponding to a homogeneous Neumann condition:

n⋅uf=0⟹n⋅(−κμ∇p)=0,x∈Γwall(7)

2.3 Heat Transfer Governing Equations

The heat transfer within the design domain is a CHT process. The steady-state energy conservation equation for the fluid is given by:

ρcpuf⋅∇Tf−∇⋅(k∇Tf)=Q(8)

where cp is the specific heat capacity at constant pressure, k is the thermal conductivity, Tf is the fluid temperature, and Q is the volumetric heat source representing heat flux from the channel walls.

Substituting the Darcy velocity field from Eq. (2) into Eq. (8) yields the explicit governing equation for the coupled fluid-thermal problem:

ρcp(−κμ∇p)⋅∇Tf−∇⋅(k∇Tf)=Q(9)

In the pseudo-3D model, the volumetric heat source term Q maps the actual 3D thermal load onto the 2D domain. Specifically, it provides an equivalent representation of the surface heat flux applied to the external boundaries of the structure. In practical aerodynamic heating scenarios, thermal loads are imposed as surface heat fluxes on the windward and leeward surfaces. To incorporate this effect into the reduced-order formulation, this boundary heat flux is mathematically transformed into an equivalent internal volumetric heat source. By integrating the heat transfer process along the thickness direction, the surface heat flux is converted into an equivalent volumetric heat source according to Q=Qf/t or Q=Qb/t, where Qf and Qb are the prescribed surface heat fluxes on the windward and leeward surfaces, respectively, and t is the effective thickness of the layer. This approximation implies a uniform distribution of heat across the thickness, which is consistent with the pseudo-3D modeling framework. While this approach neglects detailed boundary-layer effects near the wall, it preserves the global thermal behavior.

Moreover, aerodynamic heating on the fuel tank surface is highly non-uniform in actual flight conditions. For instance, the windward side experiences significantly higher heat flux than the leeward side. To capture this spatial non-uniformity in the optimization, Q is defined as a spatially varying function Q(x). This spatially varying heat load drives the optimization toward generating channel structures with asymmetric or locally intensified features, thereby concentrating thermal transport where it is most needed.

The following thermal boundary conditions are applied. At the fluid inlet, the coolant enters at a prescribed constant temperature Tin, giving the Dirichlet condition:

Tf(x)=Tin,x∈Γin(10)

The lateral walls and the outlet are assumed to be adiabatic, resulting in a homogeneous Neumann condition:

n⋅∇Tf=0,x∈Γwall∪Γout(11)

2.4 Topology Optimization Formulation

In this work, the material layout within the design domain Ωd is optimized via topology optimization. A continuous design variable field γ(x)∈ [0, 1] represents the material distribution, where γ= 0 corresponds to pure fluid, γ= 1 to pure solid, and intermediate values 0 <γ<1 represent intermediate states.

A mapping must be defined between the material properties and the design variable γ. Scalar properties such as thermal conductivity, density, and specific heat capacity are interpolated linearly:

k(γ)=kf+γ(ks−kf)(12)

ρ(γ)=ρf+γ(ρs−ρf)(13)

cp(γ)=cpf+γ(cps−cpf)(14)

where the subscripts f and s denote the fluid and solid properties, respectively.

However, a strongly penalized interpolation scheme is required for the permeability κ to drive the design toward a clear solid/fluid distribution. This penalization encourages binary designs, i.e., clear distinctions between flow channels (κ→∞) and solid regions (κ→0). A convex penalization function is adopted:

κ(γ)=κf+γpκ(κs−κf)(15)

with a penalization factor pκ=3, which significantly penalizes the permeability of intermediate material. Theoretically, κs should be zero for solid regions. However, this would make the system matrix singular, causing numerical instability. In practice, a small positive value is used for κs to maintain numerical robustness while effectively blocking flow.

The sensitivity filtering technique based on partial differential equation is employed. The filtered design variable γ~ is obtained by solving the following Helmholtz-type partial differential equation:

−Rmin2∇2γ~+γ~=γ, in Ωd(16)

where Rmin represents the filter radius, which effectively defines the minimum feature size of the design. This parameter implicitly enforces a minimum length scale in the design, effectively suppressing checkerboard patterns and imposing a lower bound on channel width. However, filtering typically produces gray, blurred interfaces. To recover a crisp binary distribution, a hyperbolic tangent projection is applied to γ~:

γ¯=tanh⁡(β(γ~−η))+tanh⁡(βη)tanh⁡(β(1−η))+tanh⁡(βη)(17)

where β controls the projection slope and η is the threshold. In the present study, β= 2 and η= 0.5 are employed. This projection sharpens the intermediate gray regions, yielding a manufacturable fluid-solid interface.

The filtering and projection steps not only improve numerical stability but also play a crucial role in controlling the geometric characteristics of the optimized topology. Specifically, the filter introduces a minimum length scale, which suppresses variations in the design field and prevents the formation of checkerboard patterns and disconnected microstructures. It can be interpreted as an approximate lower bound for the channel width, which is critical for both numerical stability and manufacturability. The projection operation further promotes the discreteness of the material distribution by enforcing a clear distinction between fluid and solid regions. These combined effects are clearly reflected in the optimized results, which are characterized by smooth boundaries, clear connectivity, and physically consistent flow pathways.

The objective is to design a flow channel layout that maximizes heat dissipation efficiency. The problem is mathematically formulated as minimizing the average temperature subject to multiple constraints:

F⁡i⁢n⁢d:γi∈[0,1],i=1,2,…,NM⁢i⁢n⁢i⁢m⁢i⁢z⁢e:Φ=∫ΩdT⁡(γ)dΩΩdS⁢u⁢b⁢j⁢e⁢c⁡t⁢t⁢o:{Governing Eqns:{∇⋅uf=0uf=−κ⁡(γ)μ∇pρ⁡(γ)⁢cp⁡(γ)⁢uf⋅∇T−∇⋅(k⁡(γ)∇T)=QVolume Constraint:∫Ωd(1−γ)⁢d⁢Ω≤Vf⋅ΩdPressure Drop:h1=1|Γi⁢n|⁢∫Γi⁢np⁢d⁢Γ−h1∗≤0Temp. Limits:h2=Tmax−h2∗≤0Temp. Uniformity:h3=(Tmax−Tmin)−h3∗≤0(18)

Here, the objective function Φ represents the spatially averaged temperature across the entire design domain Ωd. The optimization is subject to several constraints. First, a volume constraint limits the material usage, where Vf denotes the maximum allowable volume fraction. Second, the pressure drop constraint h1 requires that the average pressure at the inlet boundary Γin does not exceed the specified limit h1∗. Regarding thermal performance, Tmax and Tmin represent the maximum and minimum temperatures within the domain, respectively. The maximum temperature constraint h2 ensures that Tmax does not exceed the critical threshold h2∗ to avoid boiling or thermal decomposition of the coolant. Finally, the temperature uniformity constraint h3 restricts the global temperature difference (Tmax−Tmin) to a prescribed bound h3∗, thereby mitigating thermal stress and maintaining structural safety.

The topology optimization problem defined in Eq. (18) is solved using the Method of Moving Asymptotes. The optimization process is terminated when the relative change in the objective function between successive iterations falls below 10−3 or when the maximum number of iterations (e.g., 300) is reached. The sensitivities of the objective function and constraints with respect to the design variables are computed via the adjoint method. The filtering and projection techniques are incorporated into the sensitivity analysis through the chain rule to ensure consistency between design variables and physical properties.

3  Parameter Calibration for the Darcy Model

Coolant flow in aerospace cooling channels typically occurs at moderately high Reynolds numbers, on the order of 103 to 105, where the flow is turbulent. In this regime, flow resistance is dominated by inertial effects and nonlinear convection. Consequently, the pressure drop ΔP scales approximately with the square of the velocity, i.e., ΔP∝u2. This inherent nonlinearity contributes significantly to the computational cost of solving the N-S equations. In contrast, the Darcy flow model employed here is linear, with flow resistance arising primarily from viscous drag. This yields a linear relationship where the pressure drop is proportional to velocity, i.e., ΔP∝u. Because the model neglects inertial effects, the permeability κ must be carefully calibrated to match the macroscopic flow behavior of actual turbulent conditions. If κ is not properly calibrated, the Darcy model will yield flow fields and pressure drops that deviate significantly from reality, rendering the optimization results of limited practical value for engineering applications.

From this perspective, the calibration process establishes a mapping between the Darcy model and the high-Reynolds-number turbulent flow, ensuring that the reduced-order model reproduces the correct macroscopic behavior. Therefore, in this study, the Darcy model is not employed to capture high-fidelity local turbulent phenomena. Instead, by calibrating the effective permeability parameter against high-fidelity turbulent data, the surrogate model is constrained to match the global flow resistance of the true turbulent system. This ensures that the resulting channel configurations are physically consistent from a system-level energy balance standpoint. The localized non-linear turbulent effects omitted during the optimization phase are subsequently captured and evaluated during the final 3D validation.

Parameter calibration determines the effective permeability values (κf, κs) so that the Darcy model approximately reproduces the macroscopic behavior of turbulent flow. The calibration proceeds in two steps. First, the ratio κf/κs is chosen to prevent fluid flow in solid regions. Second, the absolute value of κf is adjusted to meet the prescribed pressure drop constraints. For reliable calibration, a benchmark case with typical flow features such as bifurcation and bend is selected, as shown in Fig. 3. To obtain reference data for calibration, turbulent flow is simulated using the RANS method coupled with the k−ω turbulence model.

Figure 3: Benchmark model for parameter calibration: (a) geometry and boundary conditions; (b) non-uniform heat source distribution.

In this case study, a non-uniform heat source is applied. The heat flux is higher in the central region and lower near the sides. Along the y-axis, it changes with transitional boundaries at y= 0.3 m and y= 0.7 m. Material thermal properties are given in Table 1, and structural dimensions and boundary condition parameters are listed in Table 2.

Initially, the fluid permeability is set to κf = 1 × 10−8 m2. We analyze the influence of the permeability ratio κf/κs by varying it from 101 to 1010. Table 3 lists the corresponding pressure drops and the average temperatures in the fluid and solid regions. Fig. 4 plots these trends. Fig. 5 shows velocity fields and streamline patterns for permeability ratios of 10 and 108 together with the turbulent flow solution for comparison. It is found that when the permeability ratio is low, the streamlines penetrate regions intended as solid walls (Fig. 5a), resulting in unphysical flow leakage. As the ratio increases, flow in solid regions is effectively suppressed and the pressure drop rises significantly. At κf/κs = 108, streamlines are confined strictly to the prescribed channels (Fig. 5b), ensuring physically valid fluid-solid boundaries. Therefore, a permeability ratio of 108 is selected.

Figure 4: Variation of pressure drop and mean temperature with κf/κs.

Figure 5: Comparison of velocity fields: (a) Darcy model with κf/κs=10; (b) Darcy model with κf/κs=108; (c) turbulent flow solution.

To quantitatively assess the accuracy of the Darcy-based model, a comparison with the turbulent flow solution is conducted. The pressure drop and average temperature for the turbulent flow in Fig. 5c are 41,582 Pa and 378.13 K. The relative deviations of pressure drop and average temperature are both within 5%, which is acceptable for engineering-level predictions. These results demonstrate that the calibrated reduced-order model achieves satisfactory accuracy for engineering applications while reducing the computational cost.

Next, κf is calibrated based on the pressure drop. With the ratio fixed at κf/κs = 108, κf is systematically varied from 10−1 to 10−10 m2. Table 4 lists the resulting pressure drops and average temperatures. The corresponding trends are plotted in the Fig. 6.

Figure 6: Variation of pressure drop and average temperature with κf.

In the double-logarithmic plot (Fig. 6), the pressure drop ΔP follows a straight line of slope −1 against κf. This confirms that κf is the main parameter governing system behavior. Meanwhile, the average temperature changes only slightly as κf varies. This means that once the ratio κf/κs is fixed, the system’s thermal dissipation performance is largely insensitive to the absolute value of κf. In engineering design, the allowable pressure drop is usually set between 0.01 and 0.02 MPa. According to the curve in Fig. 6, a fluid permeability of κf = 2.5 × 10−8 m2 keeps the system pressure drop within this range. Therefore, the calibrated permeability values provide a reliable foundation for the subsequent topology optimization.

The calibration results further demonstrate that, although the Darcy model neglects local inertial effects, it can accurately reproduce the global pressure drop and flow distribution when an appropriate permeability is selected. This indicates that the calibrated Darcy model achieves macroscopic consistency with turbulent flow, which is sufficient for guiding topology optimization.

4  Numerical Designs under Non-Uniform Heat Source

4.1 Optimization Problem Setup

All optimization cases in this section employ the material properties listed in Table 1. The design domain is a 2D rectangular region of size 1.0 m × 1.25 m. Inlet velocity is fixed at uin= 0.35 m/s, and outlet pressure is set to Pout= 0 Pa. All other boundaries (excluding inlets and outlets) are modeled as adiabatic no-slip walls, consistent with the physical conditions.

The non-uniform heat sources are applied to simulate the concentrated heating encountered in aircraft thermal protection systems. The domain is divided vertically into high- and low-heat-flux regions. The central region, y∈ [0.3 L, 0.7 L] is given a high heat flux of Qhigh= 4700 W/m2. The side regions, y∈ [0, 0.3 L) and y∈ (0.7 L, L], receive a lower heat flux of Qlow= 3300 W/m2. This step-wise heat source provides a simple representation of the concentrated heating often seen in aerospace thermal protection systems. Moreover, in the present study, the optimization formulation uses the following parameters: volume fraction Vf= 0.6, maximum temperature h2∗= 373.15 K, temperature uniformity h3∗= 40 K, and pressure drop limits h1∗= 0.01 MPa (double-inlet) or 0.02 MPa (single-inlet).

4.2 Example I: Single-Inlet Single-Outlet (SISO) Symmetric Configuration

First, we consider a SISO optimization case. Fluid enters through a centered inlet on the left boundary and exits through a symmetrically placed outlet on the right boundary. The design domain and boundary conditions are shown in Fig. 7.

Figure 7: Setup for the SISO topology optimization case: (a) dimensions and boundary conditions; (b) applied non-uniform heat source distribution.

Fig. 8 shows the optimized flow channel topology for a SISO symmetric structure under a non-uniform heat source, along with the corresponding pressure, temperature, and velocity fields. The optimized design exhibits a typical hierarchical branching structure (Fig. 8a). Near the inlet, a wide channel maintains a large fluid-filled region to promote bulk heat removal. The flow subsequently splits into multiple finer branches, ensuring efficient heat transport. Within the high-heat-flux region (y∈ [0.3L, 0.7L]), the channels are more concentrated than in the adjacent low-flux zones, creating an increased solid-fluid contact area for heat transfer. This configuration shortens conduction paths and fully utilizes the coolant’s cooling capacity. Toward the outlet, the branches gradually merge into wider collecting channels. The pressure, temperature, and velocity fields in Fig. 8 demonstrate that this channel design effectively enhances heat transfer.

Figure 8: Optimized results for the SISO case: (a) channel layout; (b) pressure field; (c) temperature field; (d) velocity field.

To quantitatively characterize the optimized channel design, several geometric metrics are introduced. The fluid volume fraction represents the proportion of the flow region within the design domain. The wetted perimeter is defined as the total length of the fluid–solid interface. The number of effective channels is estimated based on the connectivity of the flow paths in the optimized topology. These metrics provide a quantitative basis for comparing different designs beyond qualitative descriptions.

The optimized design exhibits a fluid volume fraction of 58%. The number of effective flow channels increases to 14, and the wetted perimeter is increased by 28% compared to the straight channel. This increase in channel number and wetted perimeter indicates that the optimization process promotes the formation of distributed branching structures, particularly in regions with high thermal loads. Moreover, similar quantitative trends are observed in other configurations. As a result, the coolant can more effectively access critical regions, thereby improving temperature uniformity and mitigating thermal hotspots.

4.3 Example II: Single-Inlet Dual-Outlet (SIDO) Asymmetric Configuration

This example examines a configuration with a single centered inlet on the left side and dual outlets on the upper and lower sides, as shown in Fig. 9. The two outlets are placed asymmetrically about the domain centerline. The upper outlet is located nearer the outer boundary, whereas the lower outlet is positioned closer to the centerline. This configuration presents a more complex scenario where the flow must distribute into multiple downstream branches under inherently unbalanced conditions.

Figure 9: Setup for the SIDO asymmetric topology optimization case: (a) dimensions and boundary conditions; (b) applied non-uniform heat source distribution.

Fig. 10 shows the topology optimization results for the SIDO asymmetric flow channel under non-uniform heat source conditions, along with the pressure, temperature, and velocity field distributions. Near the inlet, the channel structure covers a wide region, enhancing heat dissipation. Here, the fluid temperature stays low and can absorb heat efficiently. Further inside, the flow splits into multiple fine channels to ensure uniform flow distribution and effective heat dissipation. During this process, the fluid is heated by the heat source and its temperature gradually rises. Near the outlet, the channels merge again and the fluid temperature reaches its peak. Moreover, channels are more densely distributed in the middle region with higher heat flux. This design effectively enhances heat dissipation. This design promotes rapid heat transfer to the fluid, effectively enhancing dissipation and preventing local overheating. This demonstrates that the optimization algorithm prioritizes local cooling efficiency by maximizing the solid-fluid contact area in critical regions, while simultaneously adapting the global flow distribution to the imposed outlet constraints.

Figure 10: Optimized results for the SIDO case: (a) channel layout; (b) pressure field; (c) temperature field; (d) velocity field.

4.4 Example III: Dual-Inlet Single-Outlet (DISO) Symmetric Configuration

The DISO configuration is widely adopted in engineering to improve flow and temperature uniformity, particularly for wider cooling structures. Here, fluid enters from inlets on both sides and converges at a centered outlet on the right boundary. Fig. 11 shows the topology optimization results for the symmetric DISO structure under non-uniform heat source conditions, together with the distributions of pressure, temperature, and velocity fields. Compared to the single-inlet design, the dual-inlet layout allows coolant to enter from both sides simultaneously, ensuring coolant coverage across the full width of the structure and promoting more uniform and rapid convergence toward the central high-heat region. Within this region, a dense channel network maximizes the solid-fluid contact area, accelerating heat transfer. Near the outlet, the branching channels merge into thicker collectors.

Figure 11: Optimized results for the DISO symmetric case: (a) channel layout; (b) pressure field; (c) temperature field; (d) velocity field.

To examine the influence of inlet placement, we optimize the structures by varying the inlet location relative to the top boundary. Three cases are examined, with the distance l1 from the inlet to the top boundary set to 0.1, 0.25, and 0.4 m, respectively. The corresponding model configurations and heat source distributions are shown in Fig. 12.

Figure 12: Setup for the DISO topology optimization case with various inlet: Setup for the DISO topology optimization case: (a) dimensions and boundary conditions; (b) applied non-uniform heat source distribution.

The optimized topologies and field distributions for different inlet positions are shown in Fig. 13. Table 5 compares the pressure drop, maximum temperature, and maximum flow velocity across the four inlet positions examined. When the inlet is close to the side boundary, the main channel forms near that boundary, forcing the fluid to travel along the edge before turning toward the central hot spot. This increases the flow path length and, consequently, the pressure drop. As the inlet is moved inward, the main channel aligns more directly toward the central high-heat-flux region, yielding a more direct flow path that reduces both flow length and bending losses. Despite these variations in inlet position, all optimized designs maintain a high channel concentration in the high-heat-flux region. This demonstrates that the optimization algorithm is more sensitive to the heat source distribution than to the inlet location, with the primary objective remaining the effective mitigation of hot spots.

Figure 13: Optimized results for the DISO case with varying inlet positions. Rows correspond to different inlet distances l1 from the top boundary: (a–c) channel layouts for l1 = 0.1, 0.25, and 0.4 m; (d–f) corresponding temperature fields; (g–i) corresponding pressure fields; (j–l) corresponding velocity fields.

4.5 Example IV: Dual-Inlet Dual-Outlet (DIDO) Asymmetric Configuration

The DIDO configuration presents the most challenging fluid-thermal topology optimization problem in this study. It involves both multiple inputs and multiple outputs. This offers greater design flexibility but also demands that the optimization algorithm balance flow distribution among multiple branches. The design domain is shown in Fig. 14, where the upper outlet is positioned nearer the outer boundary, while the lower outlet lies closer to the domain centerline.

Figure 14: Setup for the DIDO topology optimization case: (a) dimensions and boundary conditions; (b) applied non-uniform heat source distribution.

Fig. 15 presents the topology optimization results for the DIDO asymmetric structure under non-uniform heat source conditions, together with the corresponding distributions of pressure, temperature, and velocity fields. The optimized design exhibits a complex network that is highly interconnected and hierarchical branching. As the fluid enters from the two inlets, it begins to split immediately. Some of the flow goes straight to the outlets, while the rest converges in the high-heat-flux region. After completing the main heat exchange there, it then splits and exits through the outlets. The channel layout consists of multiple parallel flow paths. This arrangement significantly improves the system’s cooling capacity by balancing flow distribution and maximizing heat transfer area.

Figure 15: Optimized results for the DIDO case: (a) channel layout; (b) pressure field; (c) temperature field; (d) velocity field.

4.6 3D Validation and Comparative Analysis

Although the pseudo-3D Darcy model with parameter calibration captures macroscopic behavior, it has inherent limitations. To validate the reliability and performance benefits of the optimized designs for real engineering applications, we perform a full 3D reconstruction and simulation.

4.6.1 3D Reconstruction of Channel Geometries

To evaluate the performance and physical validity of the channel layouts under realistic operating conditions, the optimized channel designs, including the SISO and DISO configurations are reconstructed into 3D geometries. In addition to the optimized configurations, two conventional channel designs, namely a straight channel and a serpentine channel, are also constructed to serve as baseline references for comparison. The reconstruction models are shown in Fig. 16.

Figure 16: 3D reconstruction of channel geometries. SISO cases: (a) optimized channel; (b) straight channel; (c) serpentine channel. DISO cases: (d) optimized channel; (e) straight channel; (f) serpentine channel.

The reconstruction is carried out through a systematic extrusion-based procedure. First, the optimized design variable field is converted into a binary field using a threshold value of 0.5. This thresholding operation distinguishes fluid regions from solid regions and ensures a clear definition of the channel boundaries. The choice of the threshold value follows common practice in density-based topology optimization and provides a balance between geometric fidelity and numerical stability.

To improve geometric quality, a smoothing operation is applied to the channel boundaries. This step removes jagged edges and staircase-like artifacts introduced by the discretized design field, thereby producing smoother channel walls and more physically realistic geometries. Such treatment is particularly important for avoiding artificial flow separation and numerical instabilities in subsequent simulations.

Subsequently, the resulting binary design is extruded along the out-of-plane direction with a thickness of 0.05 to generate the 3D structure. A solid substrate with a thickness of 0.12 is further added beneath the channel to model the structural wall. This construction preserves the optimized design while enabling a physically consistent conjugate heat transfer simulation.

From an engineering perspective, the reconstructed geometries exhibit favorable manufacturability characteristics. The minimum feature size of the channels is controlled by the filter radius used during the topology optimization process. As a result, excessively thin or disconnected structures are effectively suppressed. Moreover, the generated channel networks consist of continuous and smoothly varying flow paths, which are compatible with modern fabrication techniques such as additive manufacturing. Therefore, the proposed approach not only produces numerically optimized designs but also yields geometries with practical engineering relevance.

4.6.2 Mesh Convergence Study

To ensure the reliability and numerical accuracy of the simulations, a mesh convergence study is performed using the straight channel configuration shown in Fig. 16b as a representative baseline case. The straight channel is selected for this purpose due to its simple geometry and well-defined flow characteristics, which allow for a clear assessment of numerical discretization effects.

Three levels of mesh resolution are considered, corresponding to coarse, medium, and fine meshes. The mesh refinement is performed by systematically increasing the number of elements while maintaining similar mesh quality and distribution. The maximum temperature across the channel is selected as quantitative indicators for evaluating mesh convergence. Table 6 lists the relevant results.

The results show that the differences in maximum temperature between the medium and fine meshes are within 0.3%. These small deviations indicate that further mesh refinement has a negligible effect on the solution, confirming that the numerical results are effectively mesh-independent. Therefore, the medium mesh is adopted for all subsequent simulations, as it provides a good compromise between computational efficiency and solution accuracy.

4.6.3 Comparative Performance Analysis

To comprehensively evaluate the performance of the proposed topology optimization method, the reconstructed optimized channel designs, including the SISO and DISO configurations, are compared with two conventional channel layouts: the straight channel and the serpentine channel.

The straight channel represents the simplest cooling configuration, characterized by a direct flow path from inlet to outlet. The serpentine channel is widely used in engineering applications due to its enhanced heat transfer capability. All channel configurations are subjected to identical boundary conditions, including inlet velocity, outlet pressure, and non-uniform heat flux distribution. This ensures a fair and consistent comparison across different designs. The performance of each configuration is evaluated in terms of temperature distribution, velocity field, and pressure drop, providing a comprehensive assessment of both thermal and hydraulic characteristics.

Fig. 17 presents the distributions of pressure, temperature, and velocity fields for the SISO configuration (Fig. 8a), as well as for the straight and serpentine channel designs. A detailed quantitative comparison is provided in Table 7.

Figure 17: Comparison of flow field performance across optimized, straight and serpentine channel configurations: (a–c) pressure fields; (d–f) temperature fields; (g–i) velocity fields.

Fig. 18 presents the distributions of pressure, temperature, and velocity fields for the DISO configuration with l1 = 0.25, as well as for the straight and serpentine channel designs. A detailed quantitative comparison is provided in Table 8.

Figure 18: Comparison of flow field performance across optimized, straight and serpentine channel configurations: (a–c) pressure fields; (d–f) temperature fields; (g–i) velocity fields.

The comparative results reveal differences in the thermal-fluid performance of the considered channel configurations. For the straight channel, the coolant flows along a single direct path, resulting in limited interaction between the flow and the heated regions. Consequently, significant temperature gradients develop within the domain, and pronounced hot spots appear in regions with high heat flux. The overall thermal performance is inadequate for applications involving non-uniform heat loads.

The serpentine channel improves thermal performance by extending the flow path and increasing fluid mixing. This leads to a reduction in peak temperature. However, the repeated bending of the flow path introduces additional frictional losses and flow resistance, resulting in a substantially higher pressure drop compared to the straight channel. This highlights the classical trade-off between heat transfer enhancement and hydraulic efficiency.

In comparison, the optimized channel designs (SISO and DISO configurations) exhibit better overall performance by effectively balancing these competing factors. The optimized channel networks feature branching and merging structures that enable efficient redistribution of the coolant throughout the domain. As a result, regions with high heat flux receive increased flow supply, leading to a significant reduction in local hot spots and improved temperature uniformity. At the same time, the presence of multiple parallel flow paths reduces the overall flow resistance. Unlike the serpentine channel, where the flow is constrained to follow a single elongated path, the optimized designs allow for smoother and more distributed flow trajectories. This results in a moderate pressure drop that is significantly lower than that of the serpentine configuration while maintaining comparable or even superior thermal performance.

Overall, the results demonstrate that the proposed topology optimization method successfully resolves the trade-off between heat transfer enhancement and pressure loss. This confirms the effectiveness and practical applicability of the proposed method for thermal management under complex operating conditions.

5  Conclusions

This study presents a topology optimization method for cooling channel design under non-uniform heat sources, based on a pseudo-3D fluid-solid CHT model. By combining reduced-order modeling, parameter calibration, multi-condition optimization, and 3D validation, the proposed method addresses the challenge of balancing computational efficiency and physical fidelity in cooling channel optimization. The results demonstrate that the calibrated Darcy-based model serves as an efficient surrogate for capturing the macroscopic behavior of high-Reynolds-number flows. The proposed approach provides a practical balance between computational efficiency and physical accuracy for engineering design.

The main conclusions can be summarized as follows:

(1)   A pseudo-3D Darcy-based CHT model is developed to retain essential 3D thermal effects within a 2D optimization framework. Through systematic calibration against turbulent simulations, the model is shown to reproduce the macroscopic flow resistance and thermal performance of practical turbulent cooling flows.

(2)   Topology optimization under non-uniform heat sources yields channel networks that adapt primarily to the spatial distribution of thermal loading. Dense branching structures naturally emerge in high-heat-flux regions to enhance local heat transfer, while sparser channels are formed elsewhere to limit flow resistance. This behavior is consistently observed across different inlet-outlet configurations, indicating that thermal load plays a dominant role on topology optimization.

(3)   Three-dimensional reconstruction and validation demonstrate that the optimized designs perform better than conventional straight and serpentine channel layouts under identical operating conditions. These improvements confirm that the performance gains predicted by the pseudo-3D optimization are preserved in full 3D simulations.

To further assess the advantages of the proposed method, a comparison with existing topology optimization approaches is discussed. Compared to full 3D topology optimization methods, the present approach significantly reduces computational cost by solving a problem with a calibrated Darcy model. This leads to a substantial reduction in degrees of freedom while maintaining good agreement with high-fidelity 3D simulations. In comparison with level-set methods, the density-based formulation used in this work offers improved numerical robustness and avoids difficulties associated with boundary evolution and re-initialization. Furthermore, unlike conventional design approaches, the proposed method enables the automatic generation of complex branching channel networks, providing greater design flexibility and improved thermal-fluid performance. These comparisons indicate that the approach achieves a favorable balance between computational efficiency, physical accuracy, and design capability.

Overall, the proposed method provides a practical and physically meaningful tool for design of cooling channels under complex, non-uniform thermal environments. Although the present study is based on numerical simulations, the proposed method has been validated through comparisons with high-fidelity results, which demonstrate its accuracy in simulating thermal-fluid behavior. Experimental validation of the optimized channel designs would be a valuable extension of this work. Future work will focus on extending the method to more complex material systems and exploring its integration with advanced manufacturing techniques to further enhance its applicability to real-world thermal management problems.

Acknowledgement: Not applicable.

Funding Statement: This research was supported by the Liaoning Provincial ‘Jiebang Guashuai’ Science and Technology Program (Grant No. 2023JH1/10400048), the National Natural Science Foundation of China (Grant No. 12102079), and the Postdoctoral Fellowship Program of CPSF (Grant No. GZC20230344).

Author Contributions: The authors confirm contribution to the paper as follows: conceptualization, Jingjie He and Xiaopeng Zhang; methodology, Yuhui Jing; software, Yuhui Jing; validation, Jingjie He, and Yuhui Jing; formal analysis, Yuhui Jing; investigation, Yuhui Jing; resources, Xiaopeng Zhang; data curation, Yuhui Jing; writing—original draft preparation, Jingjie He and Yuhui Jing; writing—review and editing, Jingjie He, Yuhui Jing and Xiaopeng Zhang; visualization, Yuhui Jing; supervision, Xiaopeng Zhang; project administration, Jingjie He and Xiaopeng Zhang; funding acquisition, Jingjie He and Xiaopeng Zhang. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Hossain MM, Gu M. Radiative cooling: principles, progress, and potentials. Adv Sci. 2016;3(7):1500360. doi:10.1002/advs.201500360. [Google Scholar] [PubMed] [CrossRef]

2. Murshed SS, De Castro CN, Lourenço M, Lopes M, Santos F. A review of boiling and convective heat transfer with nanofluids. Renew Sustain Energy Rev. 2011;15(5):2342–54. doi:10.1016/j.rser.2011.02.016. [Google Scholar] [CrossRef]

3. Fan S, Duan F. A review of two-phase submerged boiling in thermal management of electronic cooling. Int J Heat Mass Transf. 2020;150(41501):119324. doi:10.1016/j.ijheatmasstransfer.2020.119324. [Google Scholar] [CrossRef]

4. Shahjalal M, Shams T, Islam ME, Alam W, Modak M, Hossain SB, et al. A review of thermal management for Li-ion batteries: prospects, challenges, and issues. J Energy Storage. 2021;39(2):102518. doi:10.1016/j.est.2021.102518. [Google Scholar] [CrossRef]

5. Bell LE. Cooling, heating, generating power, and recovering waste heat with thermoelectric systems. Science. 2008;321(5895):1457–61. doi:10.1126/science.1158899. [Google Scholar] [PubMed] [CrossRef]

6. Li W, Yang F, Alam T, Khan J, Li C. Experimental and theoretical studies of critical heat flux of flow boiling in microchannels with microbubble-excited high-frequency two-phase oscillations. Int J Heat Mass Transf. 2015;88(6):368–78. doi:10.1016/j.ijheatmasstransfer.2015.04.061. [Google Scholar] [CrossRef]

7. Tuckerman DB, Pease RFW. High-performance heat sinking for VLSI. IEEE Electron Device Lett. 1981;2(5):126–9. doi:10.1109/EDL.1981.25367. [Google Scholar] [CrossRef]

8. Peters AB, Zhang D, Chen S, Ott C, Oses C, Curtarolo S, et al. Materials design for hypersonics. Nat Commun. 2024;15(1):3328. doi:10.1038/s41467-024-46753-3. [Google Scholar] [PubMed] [CrossRef]

9. Xu H, Li M-J, Yang X-F, Du Y-X, Bi L-S. Numerical and experimental studies on an improved simulated annealing algorithm for inverse heat flux identification in transpiration cooling thermal protection systems. Int J Heat Mass Transf. 2026;260(1):128404. doi:10.1016/j.ijheatmasstransfer.2026.128404. [Google Scholar] [CrossRef]

10. Xiang X, Yang R, Sun H, Hoffman D. Performance improvement of a hybrid battery thermal management system with pseudo-3D topology optimization. Int J Heat Mass Transf. 2025;240(10):126632. doi:10.1016/j.ijheatmasstransfer.2024.126632. [Google Scholar] [CrossRef]

11. Zhang J, Zhang D, Zheng Y. Multi-objective topology optimization of flow channels in liquid-cooled plates for battery thermal management system. Struct Multidiscip Optim. 2025;68(10):203. doi:10.1007/s00158-025-04156-y. [Google Scholar] [CrossRef]

12. Guo H, Zhang X, Feng L, Wu Y, Wu Y. Experimental Nusselt number correlations for heat transfer of a single spherical particle in turbulent flow. Int J Heat Mass Transf. 2026;256(9):127939. doi:10.1016/j.ijheatmasstransfer.2025.127939. [Google Scholar] [CrossRef]

13. Prajapati YK, Bhadari P. Flow boiling instabilities in microchannels and their promising solutions–a review. Exp Therm Fluid Sci. 2017;88(2):576–93. doi:10.1016/j.expthermflusci.2017.07.014. [Google Scholar] [CrossRef]

14. Niknam SA, Mortazavi M, Li D. Additively manufactured heat exchangers: a review on opportunities and challenges. Int J Adv Manuf Technol. 2021;112(3):601–18. doi:10.1007/s00170-020-06372-w. [Google Scholar] [CrossRef]

15. Hirose KY, Yamazaki W, Vengadesan S. Investigation of thermal management in a differentially heated square cavity through topology optimization. Int J Heat Fluid Flow. 2026;120(5):110365. doi:10.1016/j.ijheatmasstransfer.2025.128205. [Google Scholar] [CrossRef]

16. Webb RL. Performance evaluation criteria for use of enhanced heat transfer surfaces in heat exchanger design. Int J Heat Mass Transf. 1981;24(4):715–26. doi:10.1016/0017-9310(81)90015-6. [Google Scholar] [CrossRef]

17. Zhou L, Gao T, Wang Y, Wang J, Gong J, Li J. Large eddy simulation of enhanced heat transfer in grooved channel with pulsating flow corresponding to hydrodynamic frequency. Int J Heat Mass Transf. 2024;218:124822. doi:10.1016/j.ijheatmasstransfer.2023.124822. [Google Scholar] [CrossRef]

18. Fawaz A, Hua Y, Le Corre S, Fan Y, Luo L. Topology optimization of heat exchangers: a review. Energy. 2022;252(8):124053. doi:10.1016/j.energy.2022.124053. [Google Scholar] [CrossRef]

19. Haertel JHK, Nellis GF. A fully developed flow thermofluid model for topology optimization of 3D-printed air-cooled heat exchangers. Appl Therm Eng. 2017;119(5):10–24. doi:10.1016/j.applthermaleng.2017.03.030. [Google Scholar] [CrossRef]

20. Alexandersen J, Sigmund O, Meyer KE, Lazarov BS. Design of passive coolers for light-emitting diode lamps using topology optimisation. Int J Heat Mass Transf. 2018;122(3–4):138–49. doi:10.1016/j.ijheatmasstransfer.2018.01.103. [Google Scholar] [CrossRef]

21. Li H, Ding X, Meng F, Jing D, Xiong M. Optimal design and thermal modelling for liquid-cooled heat sink based on multi-objective topology optimization: an experimental and numerical study. Int J Heat Mass Transf. 2019;144(8):118638. doi:10.1016/j.ijheatmasstransfer.2019.118638. [Google Scholar] [CrossRef]

22. Zhang T, Jing T, Qin F, Sun X, Li W, He G. Topology optimization of regenerative cooling channel in non-uniform thermal environment of hypersonic engine. Appl Therm Eng. 2023;219(10):119384. doi:10.1016/j.applthermaleng.2022.119384. [Google Scholar] [CrossRef]

23. Hekmatara M, Kharati-Koopaee M. Numerical study of the influence of pin fin arrangement and volume fraction on the heat transfer and fluid flow phenomena within open microchannels. Int Commun Heat Mass Transf. 2024;155(8):107595. doi:10.1016/j.icheatmasstransfer.2024.107595. [Google Scholar] [CrossRef]

24. Liu Y, Chen C, Yuan Y, Yang J, Guo Z, Shi J. Study of microchannel heat transfer characteristics based on topology optimization. Int J Therm Sci. 2025;214(2):109898. doi:10.1016/j.ijthermalsci.2025.109898. [Google Scholar] [CrossRef]

25. Sun S, Liebersbach P, Qian X. 3D topology optimization of heat sinks for liquid cooling. Appl Therm Eng. 2020;178(4):115540. doi:10.1016/j.applthermaleng.2020.115540. [Google Scholar] [CrossRef]

26. Lee JS, Yoon SY, Kim B, Lee H, Ha MY, Min JK. A topology optimization based design of a liquid-cooled heat sink with cylindrical pin fins having varying pitch. Int J Heat Mass Transf. 2021;172(3):121172. doi:10.1016/j.ijheatmasstransfer.2021.121172. [Google Scholar] [CrossRef]

27. Duan Z, Xie G, Yu B, Jin P. Multi-objective topology optimization and thermal performance of liquid-cooled microchannel heat sinks with pin fins. Case Stud Therm Eng. 2023;49(2):103178. doi:10.1016/j.csite.2023.103178. [Google Scholar] [CrossRef]

28. Han H, Han Y, Lin Y, Wang C, Korvink JG, Deng Y. Topology optimization of microchannel heat sinks for laminar flows of thermal–fluid. Appl Therm Eng. 2025;270(14):126153. doi:10.1016/j.applthermaleng.2025.126153. [Google Scholar] [CrossRef]

29. Wang G, Wang D, Liu A, Dbouk T, Peng X, Ali A. Design and performance enhancement of thermal-fluid system based on topology optimization. Appl Math Model. 2023;116:168–86. doi:10.1016/j.apm.2022.11.031. [Google Scholar] [CrossRef]

30. Xia Y, Chen L, Luo J, Tao W. Numerical investigation of microchannel heat sinks with different inlets and outlets based on topology optimization. Appl Energy. 2023;330:120335. doi:10.1016/j.apenergy.2022.120335. [Google Scholar] [CrossRef]

31. Yan K, Wang Y, Yan J. Topology optimization of two fluid heat transfer problems for heat exchanger design. Comput Model Eng Sci. 2024;140(2):1949. doi:10.32604/cmes.2024.048877. [Google Scholar] [CrossRef]

32. Bendsøe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng. 1988;71(2):197–224. doi:10.1016/0045-7825(88)90086-2. [Google Scholar] [CrossRef]

33. Suzuki K, Kikuchi N. A homogenization method for shape and topology optimization. Comput Methods Appl Mech Eng. 1991;93(3):291–318. doi:10.1016/0045-7825(91)90245-2. [Google Scholar] [CrossRef]

34. Ji Y, Dede EM, Lee J. Multiscale topology optimization with explicit de-homogenization for graded pin-fin heat sink design. Appl Therm Eng. 2026;288:129723. doi:10.1016/j.applthermaleng.2026.129723. [Google Scholar] [CrossRef]

35. Bendsøe MP, Sigmund O. Material interpolation schemes in topology optimization. Arch Appl Mech. 1999;69(9–10):635–54. doi:10.1007/S004190050248. [Google Scholar] [CrossRef]

36. Zuo W, Saitou K. Multi-material topology optimization using ordered SIMP interpolation. Struct Multidiscip Optim. 2017;55(2):477–91. doi:10.1007/s00158-016-1513-3. [Google Scholar] [CrossRef]

37. Ramalingom D, Cocquet PH, Bastide A. A new interpolation technique to deal with fluid-porous media interfaces for topology optimization of heat transfer. Comput Fluids. 2018;168(3):144–58. doi:10.1016/j.compfluid.2018.04.005. [Google Scholar] [CrossRef]

38. Xu S, Cai Y, Cheng G. Volume preserving nonlinear density filter based on heaviside functions. Struct Multidiscip Optim. 2010;41(4):495–505. doi:10.1007/s00158-009-0452-7. [Google Scholar] [CrossRef]

39. Dilgen SB, Dilgen CB, Fuhrman DR, Sigmund O, Lazarov BS. Density based topology optimization of turbulent flow heat transfer systems. Struct Multidiscip Optim. 2018;57(5):1905–18. doi:10.1007/S00158-018-1967-6. [Google Scholar] [CrossRef]

40. Hu Z, Zhang H, Wang J, Wang L. Heuristic optimality criterion algorithm for topology optimization of conjugate heat transfer problem. Int J Therm Sci. 2024;200:108949. doi:10.1016/j.ijthermalsci.2024.108949. [Google Scholar] [CrossRef]

41. Li L, Xu J, Li S. Topology optimization of microchannel structures for enhanced heat flow in liquid cooling garments. Phys Fluids. 2024;36(12):127140. doi:10.1063/5.0244742. [Google Scholar] [CrossRef]

42. Wang MY, Wang X, Guo D. A level set method for structural topology optimization. Comput Methods Appl Mech Eng. 2003;192(1–2):227–46. doi:10.1016/S0045-7825(02)00559-5. [Google Scholar] [CrossRef]

43. Allaire G, Jouve F, Toader AM. Structural optimization using sensitivity analysis and a level-set method. J Comput Phys. 2004;194(1):363–93. doi:10.1016/j.jcp.2003.09.032. [Google Scholar] [CrossRef]

44. Challis VJ, Guest JK. Level set topology optimization of fluids in Stokes flow. Int J Numer Methods Eng. 2009;79(10):1284–308. doi:10.1002/nme.2616. [Google Scholar] [CrossRef]

45. van Dijk NP, Maute K, Langelaar M, van Keulen F. Level-set methods for structural topology optimization: a review. Struct Multidiscip Optim. 2013;48(3):437–72. doi:10.1007/s00158-013-0912-y. [Google Scholar] [CrossRef]

46. Noël L, Maute K. XFEM level set-based topology optimization for turbulent conjugate heat transfer problems. Struct Multidiscip Optim. 2023;66(1):2. doi:10.1007/s00158-022-03353-3. [Google Scholar] [CrossRef]

47. Chen D, Hasegawa Y. Level-set based topology optimization in conjugate heat transfer with large solid-to-fluid thermal conductivity ratios. Appl Therm Eng. 2025;279:127626. doi:10.1016/j.applthermaleng.2025.127626. [Google Scholar] [CrossRef]

48. Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct. 1993;49(5):885–96. doi:10.1016/0045-7949(93)90035-C. [Google Scholar] [CrossRef]

49. Huang X, Xie YM. Optimal design of periodic structures using evolutionary topology optimization. Struct Multidiscip Optim. 2008;36(6):597–606. doi:10.1007/S00158-007-0196-1. [Google Scholar] [CrossRef]

50. Picelli R, Vicente WM, Pavanello R. Evolutionary topology optimization for structural compliance minimization considering design-dependent FSI loads. Finite Elem Anal Des. 2017;135(10):44–55. doi:10.1016/j.finel.2017.07.005. [Google Scholar] [CrossRef]

51. Munk DJ, Kipouros T, Vio GA. Multi-physics bi-directional evolutionary topology optimization on GPU-architecture. Eng Comput. 2018;35(3):1059–79. doi:10.1007/s00366-018-0651-1. [Google Scholar] [CrossRef]

52. Borrvall T, Petersson J. Topology optimization of fluids in Stokes flow. Int J Numer Methods Fluids. 2003;41(1):77–107. doi:10.1002/fld.426. [Google Scholar] [CrossRef]

53. Gersborg-Hansen A, Sigmund O, Haber RB. Topology optimization of channel flow problems. Struct Multidiscip Optim. 2005;30(3):181–92. doi:10.1007/s00158-004-0508-7. [Google Scholar] [CrossRef]

54. Olesen LH, Okkels F, Bruus H. A high-level programming-language implementation of topology optimization applied to steady-state Navier-Stokes flow. Int J Numer Methods Eng. 2006;65(7):975–1001. doi:10.1002/NME.1468. [Google Scholar] [CrossRef]

55. Kreissl S, Pingen G, Maute K. Topology optimization for unsteady flow. Int J Numer Methods Eng. 2011;87(13):1229–53. doi:10.1002/NME.3151. [Google Scholar] [CrossRef]

56. Deng Y, Liu Z, Zhang P, Liu Y, Wu Y. Topology optimization of unsteady incompressible Navier–Stokes flows. J Comput Phys. 2011;230(17):6688–708. doi:10.1016/j.jcp.2011.05.004. [Google Scholar] [CrossRef]

57. Nørgaard S, Sigmund O, Lazarov B. Topology optimization of unsteady flow problems using the lattice Boltzmann method. J Comput Phys. 2016;307(2):291–307. doi:10.1016/j.jcp.2015.12.023. [Google Scholar] [CrossRef]

58. Nørgaard SA, Sagebaum M, Gauger NR, Lazarov BS. Applications of automatic differentiation in topology optimization. Struct Multidiscip Optim. 2017;56(5):1135–46. doi:10.1007/s00158-017-1708-2. [Google Scholar] [CrossRef]

59. Sasaki Y, Sato Y, Yamada T, Izui K, Nishiwaki S. Topology optimization for fluid flows using the MPS method incorporating the level set method. Comput Fluids. 2019;188(2):86–101. doi:10.1016/j.compfluid.2019.05.010. [Google Scholar] [CrossRef]

60. Qian X, Dede EM. Topology optimization of a coupled thermal-fluid system under a tangential thermal gradient constraint. Struct Multidiscip Optim. 2016;54(3):531–51. doi:10.1007/s00158-016-1421-6. [Google Scholar] [CrossRef]

61. Zhao X, Zhou M, Sigmund O, Andreasen CS. A poor man’s approach to topology optimization of cooling channels based on a Darcy flow model. Int J Heat Mass Transf. 2018;116(2):1108–23. doi:10.1016/j.ijheatmasstransfer.2017.09.090. [Google Scholar] [CrossRef]

62. Makhija DS, Beran PS. Concurrent shape and topology optimization for steady conjugate heat transfer. Struct Multidiscip Optim. 2019;59(3):919–40. doi:10.1007/s00158-018-2110-4. [Google Scholar] [CrossRef]

63. Subramaniam V, Dbouk T, Harion JL. Topology optimization of conjugate heat transfer systems: a competition between heat transfer enhancement and pressure drop reduction. Int J Heat Fluid Flow. 2019;75(5):165–84. doi:10.1016/j.ijheatfluidflow.2019.01.002. [Google Scholar] [CrossRef]

64. Yu M, Ruan S, Wang X, Li Z, Shen C. Topology optimization of thermal-fluid problem using the MMC-based approach. Struct Multidiscip Optim. 2019;60(1):151–65. doi:10.1007/s00158-019-02206-w. [Google Scholar] [CrossRef]

65. Dong X, Liu X. Multi-objective optimal design of microchannel cooling heat sink using topology optimization method. Numer Heat Transf Part A: Appl. 2020;77(1):90–104. doi:10.1080/10407782.2019.1682872. [Google Scholar] [CrossRef]

66. Ji R, Qu Z, Lei R, Wang H, Meng J, Zhang J, et al. Topology optimization method of cold plates for transient conjugate heat transfer process based on continuous adjoint method. Int J Heat Mass Transf. 2026;259(4):128387. doi:10.1016/j.ijheatmasstransfer.2026.128387. [Google Scholar] [CrossRef]

67. Wang G, Yuan Y, Wang D, Chen L, Xiang S, Chen J. Multi-objective topology optimization of conjugate heat transfer in bi-fluid heat exchanger. Int Commun Heat Mass Transf. 2026;171(2):110053. doi:10.1016/j.icheatmasstransfer.2025.110053. [Google Scholar] [CrossRef]

68. Zhou L, Luo K, Pan H, Liu L, Tian F, Li W. Improvement of topological optimization of turbulent conjugate heat transfer in complex design domains by the initialization layout method based on flow field control. Int J Heat Mass Transf. 2026;260(9–10):128425. doi:10.1016/j.ijheatmasstransfer.2026.128425. [Google Scholar] [CrossRef]

69. Yoon GH. Topological design of heat dissipating structure with forced convective heat transfer. J Mech Sci Technol. 2010;24(6):1225–33. doi:10.1007/s12206-010-0328-1. [Google Scholar] [CrossRef]

70. Dede EM. Optimization and design of a multipass branching microchannel heat sink for electronics cooling. J Electron Packag. 2012;134(4):041001. doi:10.1115/1.4007159. [Google Scholar] [CrossRef]

71. Koga AA, Lopes ECC, Villa Nova HF, Lima CRD, Silva ECN. Development of heat sink device by using topology optimization. Int J Heat Mass Transf. 2013;64:759–72. doi:10.1016/j.ijheatmasstransfer.2013.05.007. [Google Scholar] [CrossRef]

72. Yaji K, Ogino M, Chen C, Fujita K. Large-scale topology optimization incorporating local-in-time adjoint-based method for unsteady thermal-fluid problem. Struct Multidiscip Optim. 2018;58(2):817–22. doi:10.1007/s00158-018-1922-6. [Google Scholar] [CrossRef]

73. Dilgen CB, Dilgen SB, Fuhrman DR, Sigmund O, Lazarov BS. Topology optimization of turbulent flows. Comput Methods Appl Mech Eng. 2018;331(1):363–93. doi:10.1016/j.cma.2017.11.029. [Google Scholar] [CrossRef]

74. Santhanakrishnan MS, Tilford T, Bailey C. Performance assessment of density and level-set topology optimisation methods for three-dimensional heat sink design. J Algorithms Comput Technol. 2018;12:273–87. doi:10.1177/174830181877901. [Google Scholar] [CrossRef]

75. Pietropaoli M, Montomoli F, Gaymann A. Three-dimensional fluid topology optimization for heat transfer. Struct Multidiscip Optim. 2018;59(3):801–12. doi:10.1007/s00158-018-2102-4. [Google Scholar] [CrossRef]

76. Alexandersen J, Sigmund O, Aage N. Large scale three-dimensional topology optimisation of heat sinks cooled by natural convection. Int J Heat Mass Transf. 2016;100(4):876–91. doi:10.1016/j.ijheatmasstransfer.2016.05.013. [Google Scholar] [CrossRef]

77. Lei T, Alexandersen J, Lazarov BS, Wang F, Haertel JHK, De Angelis S, et al. Investment casting and experimental testing of heat sinks designed by topology optimization. Int J Heat Mass Transf. 2018;127(1):396–412. doi:10.1016/j.ijheatmasstransfer.2018.07.060. [Google Scholar] [CrossRef]

78. Yu M, Wang X, Gu J, Ruan S, Li Z, Qian S, et al. A synergic topology optimization approach on distribution of cooling channels and diverse-intensity heat sources for liquid-cooled heat sink. Struct Multidiscip Optim. 2022;65(2):48. doi:10.1007/s00158-021-03113-9. [Google Scholar] [CrossRef]

79. Sun Y, Yao S, Alexandersen J. Topography optimisation using a reduced-dimensional model for transient conjugate heat transfer between fluid channels and solid plates with volumetric heat source. Struct Multidiscip Optim. 2024;67(4):45. doi:10.1007/s00158-024-03760-8. [Google Scholar] [CrossRef]

80. Haertel JHK, Engelbrecht K, Lazarov BS, Sigmund O. Topology optimization of a pseudo 3D thermofluid heat sink model. Int J Heat Mass Transf. 2018;121(4):1073–88. doi:10.1016/j.ijheatmasstransfer.2018.01.078. [Google Scholar] [CrossRef]

81. Zeng S, Kanargi B, Lee PS. Experimental and numerical investigation of a mini channel forced air heat sink designed by topology optimization. Int J Heat Mass Transf. 2018;121(3):663–79. doi:10.1016/j.ijheatmasstransfer.2018.01.039. [Google Scholar] [CrossRef]

82. Zeng S, Lee PS. Topology optimization of liquid-cooled microchannel heat sinks: an experimental and numerical study. Int J Heat Mass Transf. 2019;142:118401. doi:10.1016/j.ijheatmasstransfer.2019.07.051. [Google Scholar] [CrossRef]

83. Yan S, Wang F, Hong J, Sigmund O. Topology optimization of microchannel heat sinks using a two-layer model. Int J Heat Mass Transf. 2019;143(5):118462. doi:10.1016/j.ijheatmasstransfer.2019.118462. [Google Scholar] [CrossRef]

84. Thumma T, Mishra SR, Abbas MA, Bhatti MM, Abdelsalam SI. Three-dimensional nanofluid stirring with non-uniform heat source/sink through an elongated sheet. Appl Math Comput. 2022;421(16):126927. doi:10.1016/j.amc.2022.126927. [Google Scholar] [CrossRef]

85. Wei J, Xu X, Zhao Z, Pan Y, Hao N, Ou B, et al. Topology optimization of microchannel heat sinks for non-uniform integrated chip systems. Int Commun Heat Mass Transf. 2026;171(6):110046. doi:10.1016/j.icheatmasstransfer.2025.110046. [Google Scholar] [CrossRef]