# Topology Optimisation of Conformal Cooling Channels in Injection Moulding: A Scaling Strategy for Forced Convection in High Reynolds Number Flows

Yupeng Sun, Aniket Ghosh Dastidar, Alban Agazzi, Ronan Le Goff, Joe Alexandersen

► To cite this version:

Yupeng Sun, Aniket Ghosh Dastidar, Alban Agazzi, Ronan Le Goff, Joe Alexandersen. Topology Optimisation of Conformal Cooling Channels in Injection Moulding: A Scaling Strategy for Forced Convection in High Reynolds Number Flows. 2025. hal-05000026

HAL Id: hal-05000026

https://hal.science/hal-05000026v1

Preprint submitted on 21 Mar 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Topology Optimisation of Conformal Cooling
Channels in Injection Moulding: A Scaling
Strategy for Forced Convection in High Reynolds
Number Flows

Yupeng Sun¹, Aniket Ghosh Dastidar², Alban Agazzi²,
Ronan Le Goff², Joe Alexandersen¹*

¹*Institute of Mechanical and Electrical Engineering, University of Southern Denmark, Campusvej 55, Odense M, 5230, Denmark.

²Centre Technique Industriel de la Plasturgie et des Composites
(CTIPC), 2 Rue Pierre et Marie Curie, Bellignat, 01100, France.

*Corresponding author(s). E-mail(s): joal@sdu.dk;

**Abstract**

Conformal cooling channels for injection moulds play a vital role in improving the efficiency and quality of the moulding processes by enabling more uniform heat removal. This study develops an integrated framework for optimising such cooling channels by topology optimisation, addressing key challenges in heat transfer under turbulent flow conditions. A scaling strategy is proposed to approximate high velocity turbulent heat transfer using a lower velocity laminar flow model. The scaling strategy overcomes the limitations of traditional Darcy potential flow models by incorporating inertial effects, which are critical for designing complex cooling layouts. Verification against high-fidelity Reynolds-Averaged Navier-Stokes simulations confirms the effectiveness of the scaling model in delivering reliable and physically accurate optimisation results. The optimisation process explores multiple performance and manufacturability constraints. Pressure drop constraints reveal the trade-offs between cooling efficiency and fluid resistance. Temperature uniformity constraints significantly reduce surface temperature deviations, particularly in critical regions such as corners and protrusions. Length scale constraints ensure manufacturability by regulating feature sizes, balancing structural simplicity with thermal performance. By combining these innovations, the proposed methodology establishes a robust foundation for designing advanced conformal cooling channels, offering practical solutions for thermal management in industrial applications.
# 1 Introduction

## 1.1 Motivation

In manufacturing forming processes such as injection moulding (Wang et al, 2024a; Jahan et al, 2019; Kariminejad et al, 2024), thermoforming (Tomasoni et al, 2020), die casting (Karakoc et al, 2022; Andronov et al, 2024), and hot stamping (Yu et al, 2024), effective cooling plays a pivotal role in determining cycle times, productivity, and product quality (Dimla et al, 2005). Despite its importance, the design of cooling systems is often sidelined during tool development, with priority given to mechanical and kinematic functionalities of moulds. Moreover, traditional cooling system designs largely depend on the expertise of designers, lacking systematic approaches to achieve optimal solutions. While certain processes, such as hot stamping (Yu et al, 2024), focus on minimising maximum mould temperatures, they often overlook the significance of temperature homogeneity, which is essential for plastic injection moulding to avoid residual stresses and warping defects on the part (Torres-Alba et al, 2023). An optimal cooling system should cool parts as efficiently and uniformly as possible, all while ensuring manufacturability (Marc-Étienne et al, 2024). Achieving this balance requires a multi-disciplinary approach, blending thermal management and manufacturing constraints. However, designing such systems is time-consuming, highlighting the demand for computational methods to automate this task. Topology optimisation emerges as a powerful computational tool to address these challenges. By optimising the cooling path, topology optimisation offers a systematic method to meet performance objectives, including efficient heat removal and temperature uniformity, while adhering to manufacturability constraints.

## 1.2 Literature

Effective cooling system design is essential for injection moulding, as it directly influences cycle time, productivity, and product quality (Dimla et al, 2005). Traditional straight-drilled cooling channels have long been the dominant solution due to their simplicity and cost-effectiveness in manufacturing (Khan et al, 2014). However, their inherent geometric constraints often result in inefficient cooling performance, particularly in moulds with complex geometries, non-uniform wall thickness, or curved surfaces. These limitations can lead to localised hotspots, prolonged cooling cycles, and inconsistent part quality, significantly impacting productivity and product reliability (Park and Dang, 2017; Sachs et al, 2000).

The advent of advanced manufacturing techniques, such as metal additive manufacturing, has enabled the development of conformal cooling (CC) channels that closely follow the contours of the moulded part (Xu et al, 2001; Jahan et al, 2016; Wei et al, 2020). Compared to traditional channels, CC channels offer enhanced cooling
efficiency, improved temperature uniformity, and reduced cycle times, with reported reductions of up to 70% in some cases (Feng et al, 2021; Nguyen et al, 2023). These advantages make CC channels particularly beneficial for applications requiring high precision, such as optical components and complex-shaped moulds. Moreover, CC channels align with Industry 4.0 principles, promoting efficiency and interconnectedness in manufacturing (Santana et al, 2017). Despite their potential, designing CC channels remains a significant challenge due to the complexity of balancing thermal performance, flow efficiency, and manufacturability. Parametric approaches and heuristic algorithms have been employed to optimise the layouts of CC channels, often focussing individually on thermal performance or flow characteristics (Agazzi et al, 2013; Wu and Tovar, 2018). However, these methods frequently neglect critical factors such as inertial effects, pressure drop constraints, or manufacturability considerations, limiting their practical applicability. For example, particle swarm optimisation and other heuristic methods can optimise channel placement but struggle to generate branched networks necessary for efficient cooling (Storti and Sobotka, 2024). Similarly, studies exploring channel shapes and surface properties highlight trade-offs between thermal performance and structural integrity, but often lack comprehensive integration of these factors (Kanbur et al, 2022; Park et al, 2020; Hanzlik et al, 2024).

Topology optimisation has emerged as a transformative tool for addressing the challenges associated with CC channel design. Unlike traditional size or shape optimisation, which rely on predefined geometries, topology optimisation allows for the free distribution of material within a design space, generating innovative, non-intuitive designs. This capability makes topology optimisation particularly well-suited for designing CC channels, where complex geometries and multiple requirements must be considered simultaneously. Topology optimisation has been applied to diverse heat transfer phenomena, including heat conduction (Li et al, 1999), natural convection (Alexandersen et al, 2014), and forced convection (Dede, 2009; Yoon, 2010). Various numerical approaches, such as density-based (Dilgen et al, 2018; Pan et al, 2022; Zhao et al, 2021), level-set (Kambampati and Kim, 2020; Coffin and Maute, 2016), and MMC methods (Yu et al, 2019; Pan et al, 2022), have been adopted to address these problems, each offering unique advantages and limitations. In addition, it has found success in industrial scenarios such as motor housings (Zhou et al, 2016), LED lamps (Alexandersen et al, 2018), cold plates (Dienemann et al, 2022), underscoring its flexibility and effectiveness in complex design challenges. For a more complete overview of the field, the reader is to the recent review papers of the field (Dbouk, 2017; Alexandersen and Andreasen, 2020; Fawaz et al, 2022).

In the specific context of CC channel design, topology optimisation has shown remarkable potential. Jahan et al (2016) applied topology optimisation and additive manufacturing techniques to redesign traditional injection moulds with porous structures and CC channels, achieving enhanced cooling uniformity and reduced material usage. Building on this work, Jahan et al (2019) introduced a coupled thermal-fluid topology optimisation framework incorporating the Navier-Stokes (NS) and convection-diffusion equations, optimising material distribution of CC channels to balance flow resistance, heat conduction, and forced convection. Experimental validation demonstrated that the optimised designs achieved superior heat dissipation and
cooling performance compared to conventional moulds. Li et al (2018) utilised the boundary element method to optimise thin-wall injection moulds, focusing on minimising maximum melt temperature and its standard deviation, while accounting for geometric and physical constraints. Kambampati and Kim (2020) applied level-set topology optimisation to 2D and 3D cooling channel designs, integrating a Darcy flow model to approximate high Reynolds number flows efficiently. By constraining channel width, their designs achieved strong agreement between Darcy and high-fidelity Reynolds-Averaged Navier-Stokes (RANS) models, effectively capturing turbulent flow dynamics. Wang et al (2024b) proposed an origami-inspired strategy to simplify 3D CC channel design problems by applying a 2D topology optimisation framework. Their approach significantly reduced computational complexity while maintaining effective heat transfer and flow performance. By analysing 2D simplification styles, thermal load distribution, and Reynolds number effects, their work provided insight into optimising channel geometry for axisymmetric castings and other complex shapes. Marc-Étienne et al (2024) demonstrated the effectiveness of density-based topology optimisation for designing 3D CC channels. Their methodology incorporated manufacturability constraints based on density gradients, enabling more feasible designs. By minimising temperature standard deviation, they achieved uniform cooling, with experimental results showing a 43% reduction in standard deviation and improved heat dissipation compared to traditional designs. Navah et al (2024) advanced the field with a density-based framework tailored for laminar Navier-Stokes flow. Their study validated porous model results against body-fitted solutions and examined the sensitivity of optimal designs to hyper-parameters such as fluid volume constraints and filter radius. By comparing domain-averaged and cavity-surface-averaged temperature objectives, they highlighted how different design paradigms influenced cooling performance.

Despite these advances, several challenges remain in the application of topology optimisation to industrial CC channel design. Current methods often rely on simplified physical models, such as Darcy flow or laminar NS equations, which struggle to accurately describe physical phenomena in high Reynolds and Peclet number flows. While high-fidelity models like RANS offer greater accuracy, their computational cost is prohibitive for routine optimisation tasks. Furthermore, many studies lack comprehensive exploration of trade-offs between different functional objectives, such as manufacturability, temperature uniformity, and pressure drop constraints. Addressing these limitations requires innovative approaches to balance computational efficiency, physical fidelity, and practical constraints, ensuring optimised designs can meet the demands of real-world applications.

## 1.3 Contributions

This work aims to contribute to the design and optimisation of CC channels for injection moulding applications through topology optimisation. The key contributions of this study are:

* A scaling strategy is proposed as a way to circumvent the need to solve the full turbulent flow and heat transfer equations. Like the commonly used Darcy potential flow model (Zhao et al, 2018), this approach aims to mitigate the computational
costs associated with turbulent flow simulations. However, it offers an advantage
by to a certain extent accounting for inertial effects, improving its applicability to
industrial flow scenarios. Its performance and limitations are thoroughly analysed
and validated through comparisons with a full turbulent RANS model.

* This study proposes a framework to generate concepts for CC channels using topology optimisation. By systematically combining and analysing multiple functional requirements, including pressure drop, cooling efficiency, temperature uniformity, and feature size constraints, the approach offers a practical framework for design. The results demonstrate promising potential for effective and manufacturable solutions in real-world injection moulding applications.

The study restricts itself to 2D problems and serves as a starting point to present the approach.

## 1.4 Paper layout

The remainder of this paper is structured as follows. Section 2 introduces the injection mould cooling model, corresponding governing equations and proposed scaling strategy for forced convection topology optimisation. Section 3 describes the framework of the topology optimisation method, including the problem formulation, interpolation schemes, and continuation strategies. This section also covers the implementation of pressure, volume, temperature uniformity, and length scale constraints. Section 4 presents a series of numerical case studies to investigate the effects of pressure constraints, temperature uniformity constraints, and length scale constraints on structural performance. Finally, Section 5 summarises the key findings and conclusions of this work.

# 2 Theory

## 2.1 Problem setup

This study focuses on the optimisation of CC channels within injection moulding systems using topology optimisation techniques. The injection moulding system under investigation, illustrated in Fig. 1, comprises three main components: the polymer part ($\Omega_p$), the injection mould ($\Omega_s$), and the cooling channels ($\Omega_f$). The mould is divided into upper and lower sections to enable separation from the polymer part post-process, while the cooling channels are optimised within the light blue design region shown in the figure. The dimensions of each component are detailed in Fig. 1(b) and Tab. 1, and the material properties for the polymer, mould, and cooling medium are summarised in Tab. 2. These parameters reflect typical conditions encountered in industrial injection moulding applications.

Given the inherent complexity of coupled thermal-fluid dynamics in injection cooling systems, a number of simplifications have been employed for this initial study. The study is restricted to a two-dimensional and steady-state model. Although all practical CC design problems are three-dimensional, the complicated interactions between constraints and parameters is more easily handled in two dimensions. Importantly, the proposed optimisation framework can and will be extended to three-dimensional
<figure><img src="image_2.png" /><figcaption>(a) Geometric dimensions.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Boundary conditions.</figcaption></figure>

**Fig. 1:** Geometric dimensions and boundary conditions of the injection moulding cooling model.

<table><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Mould length</td><td><i>L</i></td><td>0.3</td><td>m</td></tr><tr><td>Part width</td><td><i>W</i><sub>1</sub></td><td>0.04</td><td>m</td></tr><tr><td>Protrusion width 2</td><td><i>W</i><sub>2</sub></td><td>0.04</td><td>m</td></tr><tr><td>Protrusion width 3</td><td><i>W</i><sub>3</sub></td><td>0.04</td><td>m</td></tr><tr><td>Protrusion width 4</td><td><i>W</i><sub>4</sub></td><td>0.04</td><td>m</td></tr><tr><td>Lower mould height</td><td><i>H</i><sub>1</sub></td><td>0.04</td><td>m</td></tr><tr><td>Protrusion height 2</td><td><i>H</i><sub>2</sub></td><td>0.04</td><td>m</td></tr><tr><td>Protrusion height 3</td><td><i>H</i><sub>3</sub></td><td>0.04</td><td>m</td></tr><tr><td>Protrusion height 4</td><td><i>H</i><sub>4</sub></td><td>0.04</td><td>m</td></tr><tr><td>Upper mould height</td><td><i>H</i><sub>5</sub></td><td>0.04</td><td>m</td></tr><tr><td>Inlet diameter</td><td><i>d</i><sub>c</sub></td><td>0.15</td><td>m</td></tr><tr><td>Offset thickness</td><td><i>d</i><sub>0</sub></td><td>0.15</td><td>m</td></tr><tr><td>Part thickness</td><td><i>d</i><sub>p</sub></td><td>0.15</td><td>m</td></tr></tbody></table>

**Table 1:** Dimensions for injection moulding cooling model.
<table><caption>**Table 2:** Physical parameters and material properties for injection moulding cooling model.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Inlet Velocity</td><td>U<sub>in</sub></td><td>1.35</td><td>m/s</td></tr><tr><td>Inlet temperature</td><td>T<sub>in</sub></td><td>323.15</td><td>K</td></tr><tr><td>Part conductivity</td><td>k<sub>p</sub></td><td>0.1</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Solid conductivity</td><td>k<sub>s</sub></td><td>30</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Fluid conductivity</td><td>k<sub>f</sub></td><td>0.6</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Part heat capacity</td><td>c<sub>pp</sub></td><td>2000</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Solid heat capacity</td><td>c<sub>ps</sub></td><td>460</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Fluid heat capacity</td><td>c<sub>pf</sub></td><td>4580</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Part density</td><td>ρ<sub>p</sub></td><td>1000</td><td>kg/m<sup>3</sup></td></tr><tr><td>Solid density</td><td>ρ<sub>s</sub></td><td>7800</td><td>kg/m<sup>3</sup></td></tr><tr><td>Fluid density</td><td>ρ<sub>f</sub></td><td>1000</td><td>kg/m<sup>3</sup></td></tr><tr><td>Fluid viscosity</td><td>μ</td><td>5.4 × 10<sup>-4</sup></td><td>Pa·s</td></tr><tr><td>Heat source</td><td>Q</td><td>1000</td><td>kW/m<sup>3</sup></td></tr></tbody></table>

problems in future work. Furthermore, the injection moulding cooling process is fundamentally transient, where transient analysis is important for accurately capturing the thermal dynamics of injection cooling (Agazzi et al, 2013). However, in scenarios with relatively stable boundary conditions, it has been observed that the outcomes of steady-state optimisation closely align with those obtained from transient analyses (Sun et al, 2024). Consequently, this study adopts a steady-state approach, striking a balance between computational efficiency and practical accuracy. The following assumptions and simplifications are made:

* A uniform volumetric heat source (Q) applied within the polymer part (Ω<sub>p</sub>). This volumetric heat source is calculated as an equivalent average heat source, derived from the energy required to cool the polymer part from its initial temperature to the target temperature during the injection moulding process.

* Perfect thermal contact is assumed between the mould (Ω<sub>s</sub>) and polymer part (Ω<sub>p</sub>), with heat conduction dominating in the solid regions and convection occurring in the fluid domain (Ω<sub>f</sub>).

The primary objective of this optimisation investigation is to achieve uniform cooling and minimise the temperature of the polymer part, while ensuring the structural feasibility of the optimised design, which are specific as follows.

* Minimising the bulk average temperature in the polymer part (Ω<sub>p</sub>), thereby improving cooling efficiency.

* Enhancing temperature uniformity across the surface of the polymer (Γ<sub>p</sub>) to reduce residual stresses and deformation.

* Reducing pressure drop within the cooling channels (Ω<sub>f</sub>) to ensure operational feasibility.

* Limiting the structural complexity of the optimised result to enhance the manufacturability of the cooling channels.
## 2.2 Governing equations

### 2.2.1 Fluid flow

It is assumed that the coolant flow can be modelled using the steady-state incompressible Navier-Stokes (NS) equations in two dimensions:

$$ \left. \begin{aligned} \rho u_j \frac{\partial u_i}{\partial x_j} - \mu \frac{\partial}{\partial x_j} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) + \frac{\partial p}{\partial x_i} - f_i(\mathbf{x}) &= 0 \\ \frac{\partial u_i}{\partial x_i} &= 0 \end{aligned} \right\} \text{for } \mathbf{x} \in \Omega \quad (1) $$

where $i \in 1, 2$, $\rho$ is the mass density of the fluid, $u_i$ is the velocity component in the $x_i$ coordinate direction, $p$ is the pressure and $f_i(\mathbf{x})$ is a spatially-dependent body force. The equations are active in the entire domain, $\Omega$, in order to facilitate topology optimisation. The body force is defined as a so-called Brinkman penalty term:

$$ f_i(\mathbf{x}) = -\alpha(\mathbf{x})u_i \quad (2) $$

which represents a resistance term. The Brinkman penalty factor, $\alpha$, is a spatially-varying parameter defined as follows for the separated domains introduced in Fig. 1:

$$ \alpha(\mathbf{x}) = \begin{cases} \alpha_{\min} & \text{if } \mathbf{x} \in \Omega_f \\ \alpha_{\max} & \text{if } \mathbf{x} \in \Omega_s \end{cases} \quad (3) $$

Since the considered domain is two-dimensional with a finite out-of-plane thickness, the value in the fluid region is set to $\alpha_{\min} = \frac{10\mu}{h^2}$ (Borrvall and Petersson, 2003; Gersborg-Hansen et al, 2005) in order to include viscous losses due to out-of-plane no-slip conditions, where $h$ is the out-of-plane thickness or height of the flow channel. Ideally, the value in the solid region would be infinite, which would theoretically ensure identically zero velocities inside the solid, effectively negating the solid domain. However, due to numerical reasons, the value is capped at a maximum value set to $\alpha_{\max} = 10^7$ throughout all examples. This maximum value was found to be sufficiently high to ensure negligible velocities in the solid regions.

### 2.2.2 Heat transfer

For the conjugate heat transfer, the steady-state convection-diffusion equation is used:

$$ \rho c_p u_i \frac{\partial T}{\partial x_i} - \frac{\partial}{\partial x_i} \left( k(\mathbf{x}) \frac{\partial T}{\partial x_j} \right) - Q(\mathbf{x}) = 0 \quad \text{for } \mathbf{x} \in \Omega \quad (4) $$

where $c_p$ is the specific heat capacity of the fluid, $T$ is the temperature field, $k(\mathbf{x})$ is the thermal conductivity, and $Q(\mathbf{x})$ is a spatially-varying volumetric heat source. Inside the solid domain ($\mathbf{x} \in \Omega_s$), it is assumed that the velocity is negligibly small and, thus,
the convection term disappears, leaving just conduction. The thermal conductivity is spatially-varying and defined as follows:

$$k(\mathbf{x}) = \begin{cases} k_f & \text{if } \mathbf{x} \in \Omega_f \\ k_s & \text{if } \mathbf{x} \in \Omega_s \end{cases} \qquad (5)$$

where $k_f$ and $k_s$ are the thermal conductivity of the fluid and solid, respectively. The volumetric heat generation is defined as:

$$Q(\mathbf{x}) = \begin{cases} Q_0 & \text{if } \mathbf{x} \in \Omega_p \\ 0 & \text{otherwise} \end{cases} \qquad (6)$$

where $\Omega_p$ is the domain of the polymer part, in which a volumetric heat generation, $Q_0$, is imposed.

## 2.3 Scaling strategies for heat transfer TO

In this study, the flow within the CC channels is characterised by a Reynolds number of approximately 10,000 calculated as $Re = \rho_f U_{in} d_c / \mu$. This value unequivocally places the flow within the turbulent regime, which is a common condition in industrial CC channel applications. Such high Reynolds numbers arise due to the combination of small channel dimensions and elevated flow rates, both of which are essential for achieving efficient convective heat transfer. While accurately capturing turbulent flow dynamics is crucial for effective cooling, directly employing high-fidelity turbulence models, such as RANS, in topology optimisation is computationally prohibitive for many users. To address this, the Darcy potential flow model has been widely used due to its simplicity and efficiency (Zhao et al, 2018). However, the complete lack of inertia and viscosity limits its applicability to scenarios where inertial effects are significant, such as in CC channels operating under turbulent conditions.

To address this limitation, we propose a scaling strategy as a computationally efficient alternative for scenarios where inertia cannot be ignored. By solving the NS equations at a reduced velocity or Reynolds number, this strategy preserves both inertial and viscous effects - albeit inertial effects are scaled down to some extent. The down-scaled velocity field is then rescaled to its original magnitude in the calculation of the heat transfer, providing proper order-of-magnitude of the temperature field and cooling heat fluxes.

### 2.3.1 Velocity and pressure scaling strategies

The core idea of the strategy is to scale down the velocity field to a regime where solving the NS equations is computationally feasible, and then rescale the resulting velocity field to its original magnitude for the heat transfer calculation. The first step is determining the scaling factor $\eta$ based on Reynolds number. For a given optimisation problem, the original Reynolds number $Re_o$ is typically known or can be estimated from the flow conditions ($Re_o = 10,000$ in this work). Next, a scaled Reynolds number $Re_s$ must be chosen, ensuring both numerical stability and computational feasibility.
The chosen Reynolds number should be sufficiently low to guarantee stable solutions for the laminar NS equations, while as close as possible to the original Reynolds number $Re_o$ to maintain similar flow characteristics. In this study, we select $Re_s = 100$, which is commonly used in previous optimisation works in the field. Unless otherwise specified, we adopt $Re_s = 100$ throughout this paper. The scaling factor $\eta$ is then determined as:

$$ \eta = \frac{Re_o}{Re_s} \qquad (7) $$

Once the scaling factor is determined, we define the scaling ratios for the velocity and pressure fields. The scaled velocity field $u_s$ and pressure field $p_s$ are related to the original fields $u_o$ and $p_o$ as follows:

$$ u_s = \frac{u_o}{\eta} \qquad p_s = \frac{p_o}{\eta_p} \qquad (8) $$

where $\eta_p$ is a scaling factor for pressure that depends on the flow regime. For laminar flow, pressure scales linearly with Reynolds number, whereas for turbulent flows, the pressure typically exhibit a quadratic or higher-order dependence on the Reynolds number (Schlichting and Gersten, 2017). Therefore, the specific value of $\eta_p$ must be determined based on the flow regime, and trial calculations are necessary to determine the appropriate scaling in each case.

After defining the scaling factor and obtaining the proportional relationship between velocity and pressure, the boundary conditions of velocity and pressure must be adjusted accordingly. The NS equations (Eq. 1) are then solved under the adjusted boundary conditions to obtain the scaled velocity and pressure field. After solving for the scaled velocity and pressure fields, the velocity field is rescaled back to its original magnitude $u_o$ by applying the inverse scaling relation:

$$ u_o = \eta u_s \qquad (9) $$

This rescaled velocity field $u_o$ is then substituted into the convective term of the heat transfer equation (Eq. 4), as shown:

$$ \rho_f c_p \eta u_s \cdot \nabla T - \nabla \cdot k_f \nabla T - Q = 0 \qquad (10) $$

The resulting convective velocity field now accurately captures the convection effects essential for heat transfer in the topology optimisation process.

### 2.3.2 Comparative verification with RANS and Darcy models

To verify the proposed scaling strategy, we compare the results of the scaling model with those obtained from RANS and Darcy models using two test cases, as illustrated in Fig. 2. The first case replicates a benchmark problem introduced in the pioneering work that first utilised the Darcy model for heat transfer topology optimisation (Zhao et al, 2018). In this example, the flow channel has a uniform width and lacks branching. The second case involves a more complex geometry, where inertial effects are significant
<figure><img src="image_4.png" /><figcaption>**Fig. 2:** Dimensions and boundary conditions of the verification case: (a) The benchmark case proposed in (Zhao et al, 2018); (b) Verification case where inertial effects are influential.</figcaption></figure>

due to a larger fluid domain and multiple branches. The parameters for these two cases are summarised in Tab. 3 and Tab. 4, respectively.

For the single-channel case, Fig. 3 shows the pressure, velocity, temperature, and convective heat flux distributions for the three models. For the scaling model, the pressure scaling factor $\eta_p$ was calculated to be $\eta_p = 8000$. While the absolute pressure levels are not discussed in detail here, as both the Darcy and scaling models require numerical trial-and-error calibration to align with the RANS results, we focus on the pressure distribution trends. The Darcy model predicts a nearly uniform pressure drop along the flow direction, lacking the localised variations captured by the RANS model. In contrast, the scaling model identifies low-pressure regions at the two channel bends, closely aligning with the RANS trends. The velocity distribution reveals further differences between the models. The RANS results clearly exhibit a gradient velocity, with
<table><caption>**Table 3:** Parameters for verification case (a).</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Side length</td><td>D</td><td>0.5</td><td>m</td></tr><tr><td>Fillet radius 1</td><td>R<sub>1</sub></td><td>0.1</td><td>m</td></tr><tr><td>Fillet radius 2</td><td>R<sub>2</sub></td><td>0.05</td><td>m</td></tr><tr><td>Fillet radius 3</td><td>R<sub>3</sub></td><td>0.00125</td><td>m</td></tr><tr><td>Inlet temperature</td><td>T<sub>in</sub></td><td>0</td><td>K</td></tr><tr><td>Solid conductivity</td><td>k<sub>s</sub></td><td>44</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Fluid conductivity</td><td>k<sub>f</sub></td><td>0.6</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Solid heat capacity</td><td>c<sub>ps</sub></td><td>460</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Fluid heat capacity</td><td>c<sub>pf</sub></td><td>4200</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Solid density</td><td>ρ<sub>s</sub></td><td>7800</td><td>kg/m<sup>3</sup></td></tr><tr><td>Fluid density</td><td>ρ<sub>f</sub></td><td>1000</td><td>kg/m<sup>3</sup></td></tr><tr><td>Fluid viscosity</td><td>μ</td><td>0.001</td><td>Pa·s</td></tr><tr><td>Solid permeability</td><td>K<sub>s</sub></td><td>2.5 × 10<sup>-11</sup></td><td>m<sup>2</sup></td></tr><tr><td>Fluid permeability</td><td>K<sub>f</sub></td><td>2.5 × 10<sup>-5</sup></td><td>m<sup>2</sup></td></tr><tr><td>Heat source</td><td>Q</td><td>200</td><td>kW/m<sup>3</sup></td></tr></tbody></table>

<table><caption>**Table 4:** Parameters for verification case (b).</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Side length</td><td>L</td><td>0.3</td><td>m</td></tr><tr><td>Inlet length</td><td>L<sub>in</sub></td><td>0.04</td><td>m</td></tr><tr><td>Part length</td><td>L<sub>p</sub></td><td>0.15</td><td>m</td></tr><tr><td>Part width</td><td>W<sub>p</sub></td><td>0.032</td><td>m</td></tr><tr><td>Side length 2</td><td>L<sub>2</sub></td><td>0.041</td><td>m</td></tr><tr><td>Side length 3</td><td>L<sub>3</sub></td><td>0.034</td><td>m</td></tr><tr><td>Side width 2</td><td>W<sub>2</sub></td><td>0.01</td><td>m</td></tr><tr><td>Side width 3</td><td>W<sub>3</sub></td><td>0.054</td><td>m</td></tr><tr><td>Fillet radius 1</td><td>R<sub>1</sub></td><td>0.016</td><td>m</td></tr><tr><td>Inlet temperature</td><td>T<sub>in</sub></td><td>323.15</td><td>K</td></tr><tr><td>Part conductivity</td><td>k<sub>p</sub></td><td>0.1</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Solid conductivity</td><td>k<sub>s</sub></td><td>30</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Fluid conductivity</td><td>k<sub>f</sub></td><td>0.6</td><td>W/m·K<sup>-1</sup></td></tr><tr><td>Part heat capacity</td><td>c<sub>pp</sub></td><td>2000</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Solid heat capacity</td><td>c<sub>ps</sub></td><td>460</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Fluid heat capacity</td><td>c<sub>pf</sub></td><td>4580</td><td>J/kg·K<sup>-1</sup></td></tr><tr><td>Part density</td><td>ρ<sub>p</sub></td><td>1000</td><td>kg/m<sup>3</sup></td></tr><tr><td>Solid density</td><td>ρ<sub>s</sub></td><td>7800</td><td>kg/m<sup>3</sup></td></tr><tr><td>Fluid density</td><td>ρ<sub>f</sub></td><td>1000</td><td>kg/m<sup>3</sup></td></tr><tr><td>Fluid viscosity</td><td>μ</td><td>5.4 × 10<sup>-4</sup></td><td>Pa·s</td></tr><tr><td>Heat source</td><td>Q</td><td>1000</td><td>kW/m<sup>3</sup></td></tr></tbody></table>

higher velocities near the channel centre and lower velocities near the walls, reflecting turbulent boundary layer effects. The Darcy model, however, predicts a nearly uniform velocity field, failing to account for inertial and boundary layer effects. The scaling model demonstrates an intermediate behaviour, capturing a velocity gradient but overestimating the centreline velocities while underestimating near-wall velocities. This discrepancy arises from the scaling model's assumption of laminar-like flow, leading to a thicker boundary layer and enhanced momentum transfer in the channel centre. Examining the temperature and convective heat flux distributions provides
<figure><img src="image_5.png" /><figcaption>**Fig. 3:** Comparison of physical fields obtained using different physical models in the single channel case.</figcaption></figure>

additional insights. The scaling model overestimates the temperature around the heat source, while the Darcy model underestimates it. Furthermore, the scaling model exhibits a distinct asymmetry in temperature distribution near the outlet, with higher temperatures on one side. This behaviour is tied to differences in convective heat flux. The Darcy model, with its higher near-wall velocities, overestimates convective effects, resulting in lower temperatures overall. Conversely, the scaling model's thicker
<figure><img src="image_6.png" /><figcaption>Fig. 4: Comparison of velocity and temperature at the cutting line obtained using different models.</figcaption></figure>

boundary layer leads to reduced near-wall convective flux, causing higher near-wall temperatures.

Figure 4 compares velocity and temperature profiles at different cross-sections. The scaling model predicts higher centreline velocities and lower near-wall velocities than the RANS model. Darcy model shows a velocity distribution closer to RANS in magnitude but lacks the ability to represent velocity gradients in the boundary layer. For temperature, the scaling model's peak values are closer to those of the RANS model but show larger deviations near the walls due to the exaggerated boundary layer effect.

For the multi-channel case, Fig. 5 shows the pressure, velocity, temperature, and convective heat flux distributions for the three models. The more complex geometry
<figure><img src="image_7.png" /><figcaption>**Fig. 5:** Comparison of physical fields obtained using different physical models in the multi-channel case.</figcaption></figure>

highlights the critical role of inertia in flow behaviour. For the scaling model, the pressure scaling factor was calculated as $\eta_p = 8600$, while the Darcy model's pressure was scaled by a factor of 38.4 based on the RANS model's pressure drop. The Darcy model predicts a uniform pressure decrease from inlet to outlet, failing to capture the localised variations seen in the RANS results. By contrast, the scaling model successfully captures these variations, reflecting the influence of branching and localised inertial effects. Similar trends are observed in the velocity distribution. The RANS results show that flow is concentrated in the central channels, with low velocities in the outer branches. The scaling model accurately reproduces this behaviour, identifying the near-zero velocities in the outer branches and the flow concentration in the centre. However, the Darcy model predicts close to uniformly distributed velocities,
<figure><img src="image_8.png" /><figcaption>**Fig. 6:** Comparison of velocity and temperature at the cutting line obtained using different models.</figcaption></figure>

incorrectly assuming almost equal flow distribution among all channels. Moreover, the Darcy model fails to capture reversed flow in the outer branches, a feature clearly present in the RANS results and partially reproduced by the scaling model. Temperature and heat flux distributions further underscore the limitations of the Darcy model. Due to its failure to predict reversed flow, the Darcy model underestimates the temperatures in regions influenced by recirculating flow. In contrast, the scaling model captures the reversed flow and its impact on temperature distribution. However, in the central channels, the scaling model underestimates convective effects due to its boundary layer thickness, resulting in higher temperatures compared to the RANS results.

Figure 6 provides a detailed comparison of temperature and velocity profiles at critical cross-sections. In the central channels, the scaling model exhibits higher centreline velocities and larger near-wall velocity gradients than the RANS model. Consequently, convective effects are overestimated, leading to slightly elevated temperatures. The Darcy model, however, produces flat velocity profiles with magnitudes far lower than those of both the RANS and scaling models, failing to capture the boundary layer effects. Near the heat sources in the outer branches, the scaling model successfully captures reversed flow, matching the RANS results although with a significantly too low velocity magnitude. In stark contrast, the Darcy model predicts a uniform, forward velocity field, entirely missing the reversed flow and associated inertial effects. This results in significant inaccuracies in both velocity and temperature distributions.
These two test cases illustrate that when inertial effects are negligible, both the Darcy and scaling models offer valid approximations for turbulent flow. The scaling model better captures flow characteristics, but the Darcy model’s slip boundary condition provides a more accurate representation of near-wall convective effects, potentially leading to better temperature predictions. However, when inertial effects become significant, the Darcy model fails to replicate the observed flow and temperature distributions. The scaling model, by incorporating inertial effects, provides a more reliable approximation. This makes the scaling model particularly advantageous during early optimisation stages, where significant structural changes occur, offering more accurate physical information for topology optimisation.

# 3 Topology optimisation

## 3.1 Optimisation problem formulation

### 3.1.1 Optimisation formulation

The general topology optimisation problem can be formulated as follows:

$$ \text{Minimise: } \phi \qquad (11) $$

$$ \text{Subject to: } g_i \le 0, i \in [1, \dots, N_g] \\ 0 \le \gamma \le 1 $$

The goal is to find the optimal structural distribution that minimises the objective function, $\phi$, while simultaneously satisfying the $N_g$ constraints $g_i$. The objective function $\phi$ and the constraints $g_i$ are both functions of the design and physical fields, whose values are obtained by solving the NS equation (Eq. 1) and conjugate heat transfer governing equation (Eq. 4). Based on the nested analysis and design (NAND) framework, the governing equation is not explicitly listed in the optimisation formulation.

In this study, the optimisation is applied to the cooling channels in injection mould design. The objective is to enhance production efficiency while meeting manufacturing process requirements. Improving the cooling performance of these channels is the primary and most essential task, making it the optimisation objective. Additionally, the uniformity of temperature during the cooling process is another critical metric. Large internal temperature gradients in the part can lead to residual stresses or deformations, which affect product quality. Therefore, temperature uniformity is also considered during the optimisation process. Furthermore, to meet manufacturing constraints and account for the complexity of design optimisation, constraints are imposed on structural feature size, volume fraction, and pressure drop in the cooling channels.

It is important to note that this study adopts a single objective with multiple constraints formulation. Although some studies on heat transfer optimisation employ multi-objective formulation, this often requires scaling different objectives to the same magnitude and carefully selecting weight coefficients, which is more suitable for exploratory research without clear functional requirements. We recommend that for
problems with clear prioritisation of objectives, using a single objective with multi-ple constraints provides a clearer direction for optimisation and yields more desirable designs.

In this research, cooling performance is intuitively expressed as the bulk average temperature, $\bar{T}_{\text{bulk}}$, of the polymer part under given flow and heat transfer boundary conditions. The optimisation objective is to minimise $\bar{T}_{\text{bulk}}$, which is defined as follows:

$$ \phi = \bar{T}_{\text{bulk}} = \frac{\int_{\Omega_p} T \, d\Omega}{\int_{\Omega_p} d\Omega} \qquad (12) $$

Temperature uniformity during the cooling process is represented by the variance in surface temperature, $T_{\text{edge}}$, with the constraint expressed as:

$$ g_1 = \frac{1}{\int_{\Gamma_p} d\Gamma} \int_{\Gamma_p} \left( \frac{T - \bar{T}_{\text{edge}}}{T_{\epsilon}} \right)^2 d\Gamma - 1 \qquad (13) $$

with

$$ \bar{T}_{\text{edge}} = \frac{1}{\int_{\Gamma_p} d\Gamma} \int_{\Gamma_p} T \, d\Gamma \qquad (14) $$

where $\bar{T}_{\text{edge}}$ represents the average surface temperature, $T_{\epsilon}$ represents the upper limit of the temperature fluctuation range, which can be roughly interpreted as $T \in [\bar{T}_{\text{edge}} - T_{\epsilon}, \bar{T}_{\text{edge}} + T_{\epsilon}]$. By adjusting $T_{\epsilon}$, different requirements for temperature uniformity can be achieved. It should be noted that this is a global integral constraint, ensuring the overall level of temperature fluctuation without imposing strict limits on local temperature variations.

The constraint on pressure drop, $\Delta p$, across the cooling channel network is expressed as:

$$ g_2 = \frac{\Delta p}{\Delta p_{\text{max}}} - 1 = \frac{1}{\Delta p_{\text{max}} \int_{\Gamma_{\text{in}}} d\Gamma} \int_{\Gamma_{\text{in}}} p \, d\Gamma - 1 \qquad (15) $$

where $\Delta p_{\text{max}}$ is the upper limit for the pressure drop and $\Gamma_{\text{in}}$ represent the inlet of the cooling channel. The integral over the outlet is left out since $p = 0$ is applied here. This constraint can be adjusted to control different pressure drops between the inlet and outlet of the channel. This constraint not only meets the practical manufactur-ing requirements for the channel's dynamic performance but, from an optimisation perspective, the pressure drop constraint combined with the Brinkman interpola-tion penalty mechanism helps drive the 0-1 distribution of the design field, thereby promoting the formation of the optimal design.

The constraint on volume fraction of the cooling channels, $V_f$, is formulated as:

$$ g_3 = \frac{V_f}{V_{\text{max}}} - 1 = \frac{1}{V_{\text{max}} \int_{\Omega} d\Omega} \int_{\Omega} \bar{\gamma} \, d\Omega - 1 \qquad (16) $$

where $V_{\text{max}}$ denotes the upper limit of the volume constraint and $\bar{\gamma}$ is the projected design field as will be outlined in Section 3.2.1. This constraint limits the volume of the
cooling channels during optimisation. Beyond manufacturing requirements, the volume constraint serves to eliminate unnecessary fluid or solid regions in the optimisation result, leading to a cleaner design, as will be demonstrated and discussed in later sections.

Lastly, constraints on the feature sizes are introduced to enhance the manufacturability, imposing limits on the maximum and minimum feature sizes for the channels, as well as the minimum feature size for solid regions. These constraints are expressed as: Maximum length scale constraint:

$$g_4 = \frac{1}{\epsilon_f^{\max} \int_{\Omega} d\Omega} \int_{\Omega} H_p(F_f(\tilde{\gamma}, r), \beta, \eta) d\Omega - 1 \quad (17)$$

Minimum length scale constraints:

$$g_5 = \frac{1}{\epsilon_f^{\min} \int_{\Omega} d\Omega} \int_{\Omega} I_i^f(\min(\tilde{\gamma} - \eta_e), 0)^2 d\Omega - 1 \quad (18)$$

$$g_6 = \frac{1}{\epsilon_s^{\min} \int_{\Omega} d\Omega} \int_{\Omega} I_i^s(\min(\tilde{\gamma} - \eta_d), 0)^2 d\Omega - 1 \quad (19)$$

where $\epsilon_f^{\max}$, $\epsilon_f^{\min}$, and $\epsilon_s^{\min}$ represent the upper bounds of the constraints on maximum channel size, minimum channel size, and minimum solid size, respectively. $H_p(\gamma, \beta, \eta)$ denotes the hyperbolic tangent projection function, $F_f(\gamma, r_{max})$ is the PDE filter function, and $I_i^f$ and $I_i^s$ are indicator functions used to select the inflection points for fluid and solid regions, respectively. The filter functions and projection methods will be detailed in Sect. 3.2.1 and the principles and parameters behind the maximum and minimum size constraints will be further explained in Sect. 3.2.2.

### 3.1.2 Interpolation scheme

The interpolation function plays a crucial role in topology optimisation by coupling optimisation variables with material properties, thereby establishing a mapping between material spatial distribution and physical response. In this study, interpolation is applied to both the Brinkman penalty factor in Eq. 3 and the thermal conductivity in Eq. 5. For the Brinkman penalty factor, we employ the Rational Approximation of Material Properties (RAMP) interpolation scheme (Stolpe and Svanberg, 2001; Alexandersen, 2022). For thermal conductivity, we use an inverse RAMP interpolation function. The functions are defined as follows:

$$\alpha(\tilde{\gamma}) = \alpha_{\min} + (\alpha_{\max} - \alpha_{\min}) \frac{1 - \tilde{\gamma}}{1 + q_{\alpha} \tilde{\gamma}} \quad (20)$$

$$k(\tilde{\gamma}) = k_s - (k_s - k_f) \frac{\tilde{\gamma}}{1 + q_k(1 - \tilde{\gamma})} \quad (21)$$

where $\tilde{\gamma}$ represents the projected design field, while $q_\alpha$ and $q_k$ are the interpolation function curve slopes.
<figure><img src="image_9.png" /><figcaption>Fig. 7: Interpolation functions.</figcaption></figure>

The interpolation functions with different values of $q_α$ and $q_k$ are illustrated in Fig. 7. The figure reveals that the interpolation functions for the Brinkman penalty factor $α$ and thermal conductivity $k$ have opposite convexity. In the above interpolation schemes, the Brinkman penalty factor $α$ corresponding to intermediate design field values is lower than that obtained through linear interpolation, whereas the thermal conductivity $k$ for intermediate design field value is higher than that obtained through linear interpolation. This strengthens the convection, promoting the formation of flow channels in the early stages of optimisation. Furthermore, the lower $α$ makes it easier to satisfy pressure drop constraints, enabling the optimisation algorithm to focus on finding the most efficient channel layout.

As the optimisation progresses, we employ a continuation strategy, gradually reducing $q_α$ and $q_k$ to achieve a final design that meets all constraints. The details of this continuation strategy will be discussed in Sect. 3.2.3.

## 3.2 Topology Optimisation Methodology

### 3.2.1 Filter and projection

In topology optimisation, filtering techniques are introduced to regularise the design field and prevent numerical instabilities such as checkerboard patterns. Additionally, the minimum bandwidth of the filtered design field is controlled by the filter radius, which allows filtering techniques to be used for controlling feature sizes. In this study, we apply the widely-used PDE filter (Lazarov and Sigmund, 2011), expressed as follows:

$$-r^2\nabla^2\tilde{\gamma} + \tilde{\gamma} = \gamma \qquad (22)$$

For clarity and simplicity, the above PDE formulation is reformulated explicitly as $\tilde{\gamma} = F_f(\gamma, r)$. The relationship between the minimum radius $R$ of the PDE-filtered design field and the filter parameter $r$ is as follows (Lazarov and Sigmund, 2011):

$$r = \frac{R}{2\sqrt{3}} \qquad (23)$$
In addition, this work applies a smooth Heaviside projection to the filtered design field $\tilde{\gamma}$ to promote a 0-1 distribution in the optimised design field layout. The expression for the smoothed Heaviside projection is as follows (Wang et al, 2011):

$$ \tilde{\gamma} = H_p(\tilde{\gamma}, \beta, \eta) = \frac{\tanh(\beta\eta) + \tanh(\beta(\tilde{\gamma} - \eta))}{\tanh(\beta\eta) + \tanh(\beta(1 - \eta))} \quad (24) $$

where $\tilde{\gamma}$ denotes the projected design field, $\beta$ represents the slope of the projection, $\eta$ denotes the threshold value for the projection. As $\beta$ approaches infinity, the projection function converges to a 0-1 step function with the step at $\eta$.

### 3.2.2 Length scale constraints

In this study, to meet the characteristic size requirements for practical applications and to enhance the manufacturability of optimised designs, minimum and maximum length constraints are introduced, as shown in Eq. 17, Eq. 18 and Eq. 19.

The minimum length scale is achieved using the geometric constraints proposed by Zhou et al (2015), which, within a three-field density framework, controls the minimum length scale without incurring additional finite element analysis costs. The principle of this method is that the minimum length scale can be implicitly enforced on the blueprint design $\tilde{\gamma}$ (with the projection threshold $\eta = \eta_i$), provided that the topology of all designs remains consistent within the threshold range $\eta \in (\eta_e, \eta_d)$, where $0 < \eta_e < \eta_i < \eta_d < 1$. A sufficient condition for this consistency can be expressed as:

$$ \tilde{\gamma}(\mathbf{x}) \geq \eta_e, \quad \forall \mathbf{x} \in \Omega_1 = \{\mathbf{x} | \tilde{\gamma}(\mathbf{x}) = 1 \text{ and } \nabla \tilde{\gamma} = 0\} \quad (25) $$

$$ \tilde{\gamma}(\mathbf{x}) \leq \eta_d, \quad \forall \mathbf{x} \in \Omega_2 = \{\mathbf{x} | \tilde{\gamma}(\mathbf{x}) = 0 \text{ and } \nabla \tilde{\gamma} = 0\} \quad (26) $$

where $\Omega_1$ and $\Omega_2$ denote the inflection region of the filtered design field within the fluid and solid regions, respectively. This condition ensures that in the inflection region, the filtered design field must be greater than threshold $\eta_e$ or less than threshold $\eta_d$, thus preserving the topology across the threshold range. To identify these inflection regions, the following indicator functions are introduced:

$$ I_i^f = \tilde{\gamma} \cdot \exp \left( \ln(\lambda) \cdot R_{\min}^2 \cdot |\nabla \tilde{\gamma}|^2 \right) \quad (27) $$

$$ I_i^s = (1 - \tilde{\gamma}) \cdot \exp \left( \ln(\lambda) \cdot R_{\min}^2 \cdot |\nabla \tilde{\gamma}|^2 \right) \quad (28) $$

where $\lambda$ controls the slope of the indicator function, allowing for adjustments in precision. In this work, we set $\lambda = 10^{-4}$ and $R_{\min}$ represents the minimum length scale. This indicator function is proposed by Yang et al (2019) based on Zhou's work (Zhou et al, 2015), refining the original function to improve detection accuracy of inflection regions. Using this indicator function, we can further propose the geometric constraints in Eq. 18 and Eq. 19, which equal zero when the conditions in Eq. 25 and Eq. 26 are satisfied and are non-zero otherwise. By integrating these across the entire design domain and constraining it to be less than $\epsilon_f^{\max}$ and $\epsilon_s^{\max}$, the conditions are ensured in a weak sense, thus guaranteeing the minimum length scale. Additionally, according
<table><caption>**Table 5:** Parameters list for continuation scheme</caption><thead><tr><th>Step</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>β</td><td>8</td><td>8</td><td>16</td><td>16</td><td>32</td><td>32</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td><td>64</td></tr><tr><td>q<sub>α</sub></td><td>10<sup>4</sup></td><td>3 × 10<sup>2</sup></td><td>3 × 10<sup>2</sup></td><td>10<sup>2</sup></td><td>10<sup>2</sup></td><td>3 × 10<sup>1</sup></td><td>3 × 10<sup>1</sup></td><td>10</td><td>3</td><td>1</td><td>0</td><td>0</td></tr><tr><td>ε<sub>max</sub></td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>10<sup>-2</sup></td><td>10<sup>-2</sup></td><td>10<sup>-2</sup></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>ε<sub>min</sub></td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>10<sup>-2</sup></td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td><td>10<sup>-5</sup></td><td>10<sup>-6</sup></td><td>10<sup>-6</sup></td><td>10<sup>-7</sup></td><td>10<sup>-7</sup></td></tr><tr><td>ε<sub>s</sub></td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td><td>10<sup>-1</sup></td><td>10<sup>-2</sup></td><td>10<sup>-2</sup></td><td>10<sup>-3</sup></td><td>10<sup>-4</sup></td><td>10<sup>-5</sup></td><td>10<sup>-6</sup></td><td>10<sup>-6</sup></td><td>10<sup>-7</sup></td><td>10<sup>-7</sup></td></tr></tbody></table>

to (Wang et al, 2011), the values of $\eta_e$ and $\eta_d$ influence the relationship between the imposed minimum length $R_{min}$ and the filter parameter $r$. When $\eta_d = 1 - \eta_e = 0.23$, it can be assumed that $R_{min}/2\sqrt{3} = r$, which is adopted in this work.

For the maximum length scale constraint, we apply a morphology-based approach following the work of Lazarov and Wang (2017). Taking channel width as an example, a secondary filtering of the projected design field $\tilde{\gamma}$ is performed using a filter parameter $r_{max} = R_{max}/2\sqrt{3}$. This operation can be expressed as:

$$ \tilde{\gamma} = F_f(\tilde{\gamma}, r_{max}) \qquad (29) $$

Following this filter, an erosion operation is conducted via a projection function:

$$ \tilde{\gamma} = H_p(\tilde{\gamma}, \beta, \eta = 0) \qquad (30) $$

This erosion operator assesses whether the feature size adheres to the maximum constraint $2R_{max}$. For designs meeting the maximum size constraint:

$$ \int_{\Omega} \tilde{\gamma} d\Omega = 0 \qquad (31) $$

By integrating the geometric constraint over the domain and constraining areas that exceed the maximum size, the desired maximum length scale is enforced in a weak sense.

### 3.2.3 Continuation scheme

In this study, a continuation scheme is applied to various parameters, including the slope parameter $\beta$ in the projection function (Eq. 24), the penalty coefficient $q_\alpha$ in the interpolation function (Eq. 20), and the upper bounds $\epsilon_f^{max}$, $\epsilon_f^{min}$, $\epsilon_s^{min}$ for the minimum and maximum length scale constraints (Eq. 17, Eq. 18 and Eq. 19). This scheme ensures a gradual adjustment of these parameters during the optimisation process, enhancing stability and convergence. The specific parameter progression is detailed in Tab. 5.

To balance optimisation stability and computational cost, these parameters are updated every 25 iterations or when the relative tolerance of the optimisation $\epsilon_{con} < 10^{-3}$, triggering the next continuation step. The influence of different $q_\alpha$ and $q_k$ values in the interpolation function is illustrated in Fig. 7. In the continuation scheme, a larger $q_\alpha$ initially implies lower flow resistance in regions with intermediate design field values, encouraging non-solid structures to emerge. As $q_\alpha$ is progressively decreased, the
interpolation function transitions to a linear form, which increases the flow resistance of intermediate regions and helps eliminate intermediate variables. The continuation strategy for $q_\alpha$ and the slope parameter $\beta$ facilitates the search for a superior structural layout and maintains the stability of the solution process. For the thermal conductivity interpolation parameter $q_k$, it was observed that its impact on the optimisation results is relatively small but beneficial to set to a high value $q_k = 10^4$ throughout.

Regarding the upper bounds $\epsilon_f^{\max}$, $\epsilon_f^{\min}$, $\epsilon_s^{\min}$ for the length scale constraints, larger values are initially set to prevent prematurely over-constraining the design process, which can lead to suboptimal or locally trapped solutions. Once the main structural layout emerges, the upper bounds are tightened in subsequent continuation steps, ensuring compliance with the design requirements without compromising performance.

## 3.3 Implementation

This study was fully implemented using COMSOL Multiphysics v6.2, with significant customisation to address the unique requirements of the optimisation problem. For the fluid flow analysis in the cooling channel networks, the weak form of the NS equations was manually implemented using the “Weak form PDE” functionality, supplemented with pressure-stabilising/Petrov-Galerkin (PSPG) and streamline-upwind/Petrov-Galerkin (SUPG) stabilisation terms specialised for topology optimisation (Alexandersen et al, 2014; Alexandersen, 2022). These stabilisation techniques are critical in maintaining numerical stability, especially for cases involving higher values of $\alpha_{max}$, where the default “Laminar flow” NS implementation in COMSOL encounters convergence issues. The heat transfer analysis within both the part and cooling channels was performed using the “Heat Transfer in Solids and Fluids” functionality.

The optimisation problem was solved using the Method of Moving Asymptotes (MMA) (Svanberg, 1987, 2002). Custom interpolation functions (Eqs. 20 and 21), filtering techniques (Eq. 22), and hyperbolic tangent projections (Eq. 24) were implemented within the optimisation framework to ensure smooth and manufacturable designs. Additionally, separate “Coefficient form PDE” nodes were employed to impose morphological constraints (Eq. 17), effectively managing maximum feature sizes in the design domain.

The computational mesh for the optimisation is illustrated in Fig. 8. The design domain was discretised using a quadrilateral mesh with an element size of $h_{mesh} = d_c/6$, supplemented by six layers of boundary layer elements along the walls in case the fluid channels are placed here. For the part geometry, triangular elements with a size of 0.6mm were utilised, while the non-design domain was meshed using the built-in mesh tool with varying element size. Linear Lagrangian discretisation was employed for the velocity, pressure, and design fields, while quadratic Lagrangian discretisation was used for the temperature field to enhance the accuracy of the thermal analysis. To balance computational efficiency and accuracy, the mesh resolution was chosen to be fine enough to capture intricate design features during the optimisation process while maintaining manageable computational costs.

For verification with the RANS model, the built-in “Turbulent Flow” module using the $k-\epsilon$ turbulence model was employed. A “Fine” mesh setting provided by COMSOL was adopted for these simulations, ensuring a high level of accuracy through adaptive
<figure><img src="image_10.png" /><figcaption>Fig. 8: Mesh for the optimisation.</figcaption></figure>

meshing based on the underlying physics and geometry. A mesh independence study
was conducted to validate the accuracy and reliability of the results, confirming that
the chosen mesh parameters were sufficient to resolve key features of the flow and
temperature fields. The detailed results of this study are omitted for brevity, but are
available upon request.

4 Results

## 4.1 Impact of pressure drop constraints on optimised design and performance

In this subsection, we address a relatively simplified case involving the optimisation of cooling channels with the objective of minimising the average bulk temperature of the part, subject to pressure drop constraints (Eq. 15) and volume constraints (Eq. 16). The volume constraint, capped at $V_{max} = 0.1$, ensures compact designs by limiting excessive branching in the cooling channels. Unless otherwise specified, this volume limit is consistently applied across all cases. The optimisation results, including the channel configurations, velocity distributions, and temperature fields, are presented in Figs. 9, 10 and 11, respectively. These figures illustrate the progression of channel designs and cooling characteristics under different pressure drop constraints. It should be noted that since the velocity scaling strategy was employed in this research, it necessitated corresponding scaling of pressure levels. The value of the scaling factor depends on the velocity and pressure levels, ranging from 100 to 10,000 as outlined in earlier sections. To accurately impose pressure drop constraints, multiple trial simulations are needed using the RANS model to determine the appropriate scaling coefficient; this methodology is further detailed in (Geng and Zhou, 2023), which has validated the effectiveness of the strategy. Here, the primary objective is to investigate the relationship between optimised channel configurations and varying constraints while providing insights into scaling factors obtained through this case study; hence, the unscaled pressure drop is presented.
<figure><img src="image_11.png" /><figcaption>**Fig. 9:** Optimised design with varying pressure drop constraints.</figcaption></figure>

Visual analysis of the channel designs in Fig. 9 reveals distinct trends as the pressure drop constraint is progressively relaxed and a increase pressure drop is allowed. For low pressure drop, the optimised channels adopt an umbrella-like configuration in the upper mould, while the lower mould features a straight channel. This configuration minimises flow resistance by providing the shortest path between the inlet and outlet, as shown in Figs. 9(a). With relaxed constraints, the channels conform more closely to the part’s surface, as seen in Fig. 9(b). At sufficiently high pressure drop allowances, the designs transition to fully CC channels, as shown in Figs. 9(c) and (d). Another notable trend is the variation in channel width. For low pressure drops (Fig. 9(a)), the channels are widened<fnref n="1" /> to reduce internal velocities and maintain low pressure drops. Conversely, for higher pressure drops (Fig. 9(d)), the channels narrow, especially near corners, where increased flow velocity (Figs. 10c and 10d) enhances

<footnote n="1">¹For $\Delta p_{max}$ = 0.5Pa, the topology is the same, but with even wider channels. The figure has been left out to reduce space usage.</footnote>
local convective heat transfer. This adaptive behaviour improves cooling performance in areas of insufficient cooling.

Velocity and temperature distributions provide further insights into the performance of the optimised designs, shown in Figs. 10 and 11, respectively. Fig. 10 shows the flow field obtained using both the presented approach (re-scaled values) and the RANS verification model. The velocity field is visualized only in the fluid domain, with negligible values in the solid regions. Overall, the flows look quantitatively very similar. However, there are some important differences, especially in the secondary channels (Fig. 10b) and near channel constrictions (Figs. 10c and 10d). For design (a), characterised by its simple unbranched channel geometry, both the proposed optimisation model and the RANS verification model exhibit nearly identical velocity distributions (Fig. 10a). The primary discrepancy between the two models is localised near boundary regions, where minor velocity variations arise. In contrast, design (b) reveals significant differences in velocity due to the presence of narrow secondary channels (Fig. 10b). For designs (c) and (d), the global velocity distributions show intermediate agreement between verification models. However, magnified views reveal critical boundary-layer discrepancies: the RANS model has larger velocities at the boundaries, while the optimisation model obtains smaller velocities at the boundaries (Figs. 10c and 10d).

Figure 11 illustrates the temperature distributions obtained using the RANS verification model, alongside the relative error between the RANS model and the optimisation model. The relative error is calculated as follows:

$$ERROR = \frac{|T_{OPT} - T_{RANS}|}{T_{RANS} - T_{Ref}} \qquad (32)$$

where $T_{OPT}$ represents the temperature obtained from Brinkman optimisation model, $T_{RANS}$ represents the temperature obtained from RANS model, and $T_{Ref} = 273.15$K. Notably, the observed errors originate from two coupled effects: (i) the Brinkman model's inherent approximations in resolving velocity fields across solid-fluid interfaces, and (ii) discrepancies introduced by scaling strategies that employ laminar model to approximate turbulent flow. From the figures, it is evident that the relative error is predominantly controlled around 10%, with minimal discrepancies observed in the temperature of the polymer part. The majority of the error is concentrated in regions between the cooling channel and the part, where the flow-field interactions are more complex. Since the primary focus of the optimization process is on maintaining the temperature of the polymer part, the observed accuracy of the model is considered acceptable for practical applications. Additionally, a comparison of temperatures obtained during the optimisation process using the scaling model with those predicted by the RANS verification model reveals an unexpected trend. Specifically, the temperatures from the optimisation model are lower than those from the RANS model. This contradicts the findings in earlier sections, where the scaling model was shown to overestimate thermal resistance at boundaries, resulting in higher predicted temperatures compared to the RANS model.

Quantitative comparisons are summarised in Tab. 6, which lists the average bulk temperature, actual pressure drop, and channel volume fraction for each design.
<figure><img src="image_12.png" /><figcaption>**Fig. 10:** Velocity of the optimised design with varying pressure drop constraints. For the optimisation model, only the velocities in the fluid areas are shown, although near-zero velocities are present in the solid.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 11: Temperature of the optimised design with varying pressure drop constraints.</figcaption></figure>
<table><caption>**Table 6:** Comparison of temperature, pressure drop and volume fraction for different optimised designs.</caption><thead><tr><th rowspan="2">Design</th><th colspan="2">Averaged temperature [K]</th><th colspan="2">Pressure drop [Pa]</th><th rowspan="2">Volume fraction [%]</th></tr><tr><th>Optimisation</th><th>RANS</th><th>Optimisation</th><th>RANS</th></tr></thead><tbody><tr><td>Δp<sub>max</sub> = 0.5</td><td>389.16</td><td>392.26</td><td>0.61</td><td>929.41</td><td>15.45</td></tr><tr><td>Δp<sub>max</sub> = 1.0</td><td>388.70</td><td>391.59</td><td>1.00</td><td>922.54</td><td>10.85</td></tr><tr><td>Δp<sub>max</sub> = 2.0</td><td>369.27</td><td>372.69</td><td>2.00</td><td>1666.15</td><td>10.14</td></tr><tr><td>Δp<sub>max</sub> = 3.0</td><td>368.86</td><td>370.97</td><td>3.00</td><td>4398.29</td><td>9.97</td></tr><tr><td>Δp<sub>max</sub> = 4.0</td><td>368.78</td><td>370.73</td><td>4.00</td><td>6020.91</td><td>9.63</td></tr><tr><td>Δp<sub>max</sub> = 5.0</td><td>368.69</td><td>370.56</td><td>5.01</td><td>10271.55</td><td>9.29</td></tr></tbody></table>

Results from both the optimisation and RANS models are included to evaluate the accuracy of the optimisation framework. The data indicates that as the pressure drop constraint is relaxed, the average bulk temperature decreases, a trend consistently observed in both the optimisation and RANS models. This makes sense due to the observed shift to more close conformed channels. Relative errors between the two models are approximately 5% across all designs, indicating a uniform level of agreement. This discrepancy arises from differences in the treatment of thermal and flow resistances, where the optimisation model slightly underestimates convective effects compared to the RANS model. Nonetheless, the optimisation results obtained using the scaling model correctly capture the underlying trends and performance improvements, as confirmed by the RANS simulations. This consistency underscores the robustness of the scaling-based optimisation framework, despite minor prediction discrepancies. Pressure drop comparisons further reveal that the scaling factors between the optimisation and RANS models vary between 1,000 and 2,000 across different designs, which falls within the broader range of 100 to 10,000 discussed previously. This variability highlights the necessity of recalibrating the scaling factor for each design to ensure accurate pressure control. Multiple trial simulations are often required to determine the appropriate scaling coefficient, as detailed in (Geng and Zhou, 2023), where the methodology for achieving this accuracy has been verified.

Figure 12 provides a clearer perspective by plotting the trends of average bulk temperature and channel volume fraction as functions of the actual pressure drop. To better analyse the performance, the results are categorised into three pressure drop ranges: low-pressure drop (900-1600), medium-pressure drop (1600-4300), and high-pressure drop (4300-10,000).

In the low-pressure drop range (900-1600), the average temperature decreases significantly with increasing pressure drop, indicating that this region benefits the most from increasing the pressure drop constraint. However, the temperature levels remain relatively high in this stage, which suggests that pressure drop plays a dominant role in determining the cooling structure. The steep gradient of temperature reduction in this range corresponds to the transition from basic umbrella-like channels to CC branches that begin to closely follow the part surface, leading to substantial cooling performance improvements.

In the medium-pressure drop range (1600-4300), the temperature continues to decrease with increasing pressure drop, albeit at a slower rate compared to the
<figure><img src="image_14.png" /><figcaption>**Fig. 12:** Relationship between average temperature, real pressure drop, and volume fraction.</figcaption></figure>

low-pressure drop range. Both pressure drop and temperature levels fall into an intermediate range, making this stage a balance point for cooling performance and channel volume fraction. In this range, the optimisation process focuses on refining the channel geometry to enhance localized heat transfer while maintaining reasonable channel dimensions. This range demonstrates a practical trade-off between cooling performance and design constraints, offering a viable operational zone for CC channel systems.

In the high-pressure drop range (4300-10,000), the rate of temperature reduction diminishes substantially as the pressure drop increases. While the channel volume fraction continues to decrease, the improvement in cooling efficiency becomes marginal. This diminishing return underscores the nature of Pareto frontiers in optimisation problems, where pushing for further gains in cooling performance becomes increasingly inefficient. For practical CC design, investing in significantly higher pressure drops within this range provides little additional benefit in terms of thermal performance improvement, making it an impractical choice from an engineering perspective.

## 4.2 Temperature uniformity constraints in cooling channel optimisation

### 4.2.1 Impact of varying $T_\epsilon$ on optimised designs

This subsection delves into a more complex scenario, where additional temperature uniformity constraints are incorporated, as formulated in Eq. 13. These constraints are implemented alongside the pressure drop and volume constraints detailed in Sect. 4.1. Specifically, a scaled pressure drop of $\Delta p_{max} = 3$ Pa (corresponding to an actual drop
<figure><img src="image_15.png" /><figcaption>**Fig. 13:** Optimised designs with varying $T_\epsilon$.</figcaption></figure>

of 4398 Pa) is applied, while varying the parameter $T_\epsilon$ to evaluate its impact on the optimised designs.

Fig. 13 presents the resulting designs for different values of $T_\epsilon$. It can be seen that tightening the temperature uniformity constraint (i.e., reducing $T_\epsilon$) results in noticeable modifications to the channel layout. At $T_\epsilon = 2.0$, the design closely resembles the unconstrained case. When $T_\epsilon$ is reduced to 1.0, the channels deviate from the part surface, forming distinct "Ω"-shaped structures in protruding sections while introducing solid regions near part corners to compress flow channels. For $T_\epsilon = 0.8$, the optimisation produces more intricate features, such as further displacement of channels from the protruding sections, expanded solid regions at the corners, and a newly formed branch near the inlet, significantly altering the topology.

Fig. 14 provides a direct comparison of part surface temperature profiles for designs with different $T_\epsilon$. Figs. 14a and 14b show the surface temperature of the part for the designs with $T_\epsilon = 2$ and $T_\epsilon = 0.8$, while Fig. 14c contrasts the temperature profiles along the part surface for all designs. For the stricter constraint ($T_\epsilon = 0.8$), the optimisation achieves enhanced temperature uniformity across most of the part surface.
<figure><img src="image_16.png" /><figcaption>**Fig. 14:** Comparison of surface temperature profiles for designs with varying $T_e$: (a) Surface temperature contour of design with $T_e = 2$; (b) Surface temperature contour of design with $T_e = 0.8$; (c) Surface temperature deviations (relative to each design's mean contour temperature) along the part.</figcaption></figure>

However, localised regions with slightly elevated temperatures remain, particularly in areas with complex geometric features or abrupt changes in thermal boundary conditions. These elevated temperature regions, commonly referred to as hotspots in industrial applications, are primarily located at corner positions (e.g. C1 and C5) and protruding features (e.g. P1 and P2). Hotspots represent areas of concentrated heat accumulation that are notoriously challenging to mitigate due to geometric limitations and the inherent constraints of thermal diffusion. While localised cooling strategies, such as placing channels closer to these regions, can provide some relief, the fundamental limitations of thermal transport properties mean that complete elimination of hotspots is rarely achievable in practice. Interestingly, the regions between hotspots, such as the areas along the C3-P2-C4 path, benefit the most from the enhanced uniformity introduced by the design. The optimised channel abandons the conformal layout in these regions, instead opting to position channels further away from the part as seen in Figs. 13b and 13c. This strategy leads to a slight increase in overall average temperature but enhances temperature uniformity. It is important to note that the part position markings introduced here (e.g., C1-C5 and P1-P2) will be consistently used throughout the work to facilitate discussions. Unless significant changes occur in the design geometry, these markings will remain the reference framework for analysing temperature distributions and design performance.

The quantitative results in Tab. 7 corroborate the observed trend of improved temperature uniformity with decreasing $T_e$. Both the TO and RANS models demonstrate consistent trends, with parameters responding similarly to variations in $T_e$. However,
<table><caption>**Table 7:** Comparison of temperature, pressure drop and uniformity for different optimised designs.</caption><thead><tr><th rowspan="2">Design</th><th colspan="2">T<sub>bulk</sub> [K]</th><th colspan="2">Pressure drop [Pa]</th><th colspan="2">T<sub>profile</sub> [K]</th><th colspan="2">Std T<sub>profile</sub> [-]</th></tr><tr><th>TO</th><th>RANS</th><th>TO</th><th>RANS</th><th>TO</th><th>RANS</th><th>TO</th><th>RANS</th></tr></thead><tbody><tr><td>Without constraint</td><td>368.86</td><td>370.97</td><td>3.00</td><td>4398.29</td><td>325.69</td><td>327.72</td><td>1.04</td><td>1.52</td></tr><tr><td>T<sub>e</sub> = 2.0</td><td>368.96</td><td>371.25</td><td>3.00</td><td>4343.56</td><td>325.81</td><td>328.01</td><td>1.03</td><td>1.58</td></tr><tr><td>T<sub>e</sub> = 1.0</td><td>369.15</td><td>371.76</td><td>3.00</td><td>3756.82</td><td>326.01</td><td>328.56</td><td>0.91</td><td>1.51</td></tr><tr><td>T<sub>e</sub> = 0.9</td><td>369.09</td><td>372.08</td><td>3.00</td><td>3518.27</td><td>325.96</td><td>328.91</td><td>0.84</td><td>1.51</td></tr><tr><td>T<sub>e</sub> = 0.8</td><td>369.22</td><td>371.84</td><td>3.00</td><td>3807.31</td><td>326.11</td><td>328.72</td><td>0.78</td><td>1.28</td></tr></tbody></table>

discrepancies are evident in the standard deviation of the surface temperature. In the TO model, designs with lower $T_e$ constraints exhibit notably reduced deviations compared to the optimised design without the uniformity constraint, reflecting the effectiveness of the uniformity constraint. In contrast, the RANS model fails to replicate the same degree of reduction, also for designs with $T_e = 1$ and $T_e = 0.9$, where the predicted improvements are less pronounced. These differences may arise from two primary sources: inaccuracies in extracting optimised geometries during the design process and the scaling model's underestimation of convective effects at boundaries. Despite these discrepancies, the overall alignment between the two models regarding performance trends is still encouraging. The scaling model's ability to introduce intricate structural details, leading to valid temperature control, offers a practical pathway for optimised cooling designs.

### 4.2.2 Analysis of structural features in the $T_e = 0.8$ design

To verify the performance of the $T_e = 0.8$ optimised design, this subsection investigates the removal of specific structural features and their influence on temperature performance. Fig. 15a highlights the additional branch and solid features in the optimised design, while Figs. 15b-d show designs with features removed: (b) without feature 1 (lower branch), (c) without feature 2 (upper branch), and (d) without feature 3 (solid region in the channel).

Fig. 16 compares the temperature deviations (relative to each design's mean contour temperature) along the part surface for the original design and the three reference designs. Removing feature 1 significantly increases temperature deviation near C1, indicating its critical role in controlling temperature in this region. The design without feature 2 shows minimal changes in deviation compared to the original design, except for slightly higher values near P1, suggesting that feature 2 primarily controls fluid pressure loss. Conversely, removing feature 3 results in notable deviations at C3, C4, and C5, hotspot areas directly associated with its location, but little impact elsewhere. These observations underscore the localised influence of feature 3 on temperature control, while its effect on overall uniformity is less pronounced.

Quantitative comparisons in Tab. 8 reinforce these findings. All feature-removed designs exhibit higher standard deviations in temperature compared to the optimised design, reflecting reduced uniformity. Interestingly, while removing feature 2 slightly reduces bulk and profile temperatures, it also causes a sharp increase in pressure drop, highlighting its importance in maintaining flow constraints. These results
<figure><img src="image_17.png" /><figcaption>**Fig. 15:** Display of cooling channels with features removed: (a) Optimised design with $T_e = 0.8$; (b) Design with feature 1 removed; (c) Design with feature 2 removed; (d) Design with feature 3 removed.</figcaption></figure>

<table><caption>**Table 8:** Comparison of temperature, pressure drop and uniformity for designs with different features removed using RANS verification model.</caption><thead><tr><th>Design</th><th><i>T</i><sub>bulk</sub> [K]</th><th>Pressure drop [Pa]</th><th><i>T</i><sub>profile</sub> [K]</th><th>Std <i>T</i><sub>profile</sub> [-]</th></tr></thead><tbody><tr><td>Optimised design</td><td>371.84</td><td>3807.31</td><td>328.72</td><td>1.28</td></tr><tr><td>Remove feature 1</td><td>372.54</td><td>3833.35</td><td>329.46</td><td>1.93</td></tr><tr><td>Remove feature 2</td><td>371.73</td><td>8094.87</td><td>328.57</td><td>1.50</td></tr><tr><td>Remove feature 3</td><td>372.25</td><td>3703.66</td><td>329.12</td><td>1.79</td></tr></tbody></table>

demonstrate the intricate balance between structural features and overall design performance, where individual features play distinct roles in temperature control and flow regulation.
<figure><img src="image_18.png" /><figcaption>**Fig. 16:** Comparison of profile temperature along the part surface for different designs using RANS verification model.</figcaption></figure>

### 4.2.3 Comparison between Scaling and Darcy models

To further evaluate the impact of the proposed scaling model and the Darcy model on the optimisation results, this subsection compares the optimised designs obtained under the temperature uniformity constraint ($T_e = 0.8$). The optimisation formulation follows the same setup as in the previous subsection. Two key points are noted here: Firstly, due to the discrepancy between the pressure levels predicted by the Darcy model and the actual pressure levels observed in the RANS model, it is challenging to ensure the imposed pressure constraints are consistent with those of the scaling model. To address this, two optimised designs are presented for the Darcy model, corresponding to pressure drops of $\Delta p_{max} = 6000$ Pa and $\Delta p_{max} = 8000$ Pa. These designs provide a rough evaluation of the differences between the models, with one design exhibiting a slightly lower pressure drop than that achieved with the scaling model and the other slightly higher. Secondly, while attempting to apply the same interpolation and continuation strategies to the Darcy model as used in the scaling model, it was observed that the Darcy model exhibited a higher dependency on penalisation within the interpolation functions. This necessitated modifications to the interpolation function and continuation strategies to ensure the formation of clear fluid channels, the details of which can be found in Appendix A.

Fig. 17 illustrates the optimised designs obtained using the Darcy model. The design corresponding to $\Delta p_{max} = 6000$ Pa features a primary channel layout similar to the umbrella-shaped design observed for $\Delta p_{max} = 2$ Pa in Sect. 4.1. This umbrella-shaped main channel ensures a low overall pressure drop, while the branches aim to reduce localised temperature variations. The design for $\Delta p_{max} = 8000$ Pa closely resembles the layout for $\Delta p_{max} = 5$ Pa in Sect. 4.1, where channels are positioned
<figure><img src="image_19.png" /><figcaption>**Fig. 17:** Optimised designs with Darcy model: (a) Optimised with $\Delta p_{max} = 6000$ Pa; (b) Optimised with $\Delta p_{max} = 8000$ Pa.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 18:** Comparison of profile temperature along the part surface for different designs using RANS verification model.</figcaption></figure>

close to the part surface. However, due to the temperature uniformity constraint, these channels maintain a certain distance from the surface. Compared to the optimised design obtained with the scaling model, the Darcy model's designs tend to position channels closer to the part surface.

Fig. 18 presents the surface temperature deviations of the part, as calculated using the RANS verification model for the different optimised designs. It can be observed that the temperature deviation for the $\Delta p_{max} = 6000$ Pa design from the Darcy model
<table><caption>**Table 9:** Comparison of temperature, pressure drop and uniformity for different optimised designs.</caption><thead><tr><th rowspan="2">Design</th><th colspan="2"><i>T</i><sub>bulk</sub> [K]</th><th colspan="2">Pressure drop [Pa]</th><th colspan="2"><i>T</i><sub>profile</sub> [K]</th><th colspan="2">Std <i>T</i><sub>profile</sub> [-]</th></tr><tr><th>TO</th><th>RANS</th><th>TO</th><th>RANS</th><th>TO</th><th>RANS</th><th>TO</th><th>RANS</th></tr></thead><tbody><tr><td>Darcy &Delta;<i>p</i> = 6000</td><td>369.19</td><td>372.82</td><td>5999.70</td><td>2675.60</td><td>326.10</td><td>329.6</td><td>0.80</td><td>1.78</td></tr><tr><td>Darcy &Delta;<i>p</i> = 8000</td><td>368.73</td><td>371.10</td><td>8000.00</td><td>5842.72</td><td>325.62</td><td>327.93</td><td>0.80</td><td>1.28</td></tr><tr><td>Scaling</td><td>369.22</td><td>371.84</td><td>3.00</td><td>3807.31</td><td>326.11</td><td>328.72</td><td>0.78</td><td>1.28</td></tr></tbody></table>

is significantly higher than that of the other two designs. For the $\Delta p_{max} = 8000$ Pa design and the scaling model's design, the overall deviation levels are comparable, but notable local differences exist. For instance, at corner C1, the scaling model's design exhibits significantly lower temperature deviation. Structurally, this difference corresponds to the scaling model's curved and branched channel leading from the inlet to C1, compared to the straight channel in the Darcy model. A similar phenomenon is observed between C5 and C1 on the lower surface of the part, where the scaling model again shows reduced deviation due to its curved channel design. The root cause of these differences lies in the inability of the Darcy model to account for inertial effects, which the scaling model inherently incorporates.

To facilitate quantitative analysis, Tab. 9 summarises the bulk temperature ($T_{bulk}$), temperature profile ($T_{profile}$), and standard deviation of the temperature for different designs as evaluated by the RANS model. For $T_{bulk}$ and $T_{profile}$, the $\Delta p_{max} = 8000$ Pa design exhibits the lowest values, followed by the scaling model's design and then the $\Delta p_{max} = 6000$ Pa design, showing a direct correlation with pressure drop. Regarding the standard deviation, the $\Delta p_{max} = 8000$ Pa design requires nearly double the pressure drop to achieve the same standard deviation as the scaling model's design, highlighting the advantage of the scaling model in managing complex channel distributions. Additionally, the numerical results obtained using the optimised mesh align well with those obtained using body-fitted RANS meshes, indicating comparable accuracy levels for both models.

## 4.3 Additional length scale constraints

The importance of intricate structural features for ensuring temperature uniformity was demonstrated in the previous section. However, these features often pose significant challenges for practical manufacturing. To address this issue, additional length scale constraints (Eqs. 17, 18, 19) are introduced to simplify structural features or adapt them to specific requirements. In this section, the pressure drop constraint is set as $\Delta p_{max} = 3$ Pa and the temperature uniformity constraint is defined as $T_{\epsilon} = 0.8$. The impact of length scale constraints on structural performance was evaluated through two scenarios: In the first scenario, the minimum channel diameter was constrained to $d_c$, while the maximum diameter was set to $2d_c$. In the second scenario, the channel diameter was made uniform, with both the minimum and maximum diameters constrained to $d_c$.

Figure 19 compares the optimised designs under these constraints with a reference design without length scale constraints. The reference design exhibits channel dimensions ranging between $d_c$ and $2d_c$, except near the protrusions, where dimensions are
<figure><img src="image_21.png" /><figcaption>**Fig. 19:** Optimised designs with different length scale constraints.</figcaption></figure>

smaller than $d_c$. Solid regions within the channels also include features smaller than $d_c$. In contrast, the design with $R_{max} = 2d_c$ (Fig. 19b) eliminates sub-$d_c$ features, ensuring that all structural components respect the defined constraints. Furthermore, the uniform design (Fig. 19c) achieves consistent channel widths of $d_c$, demonstrating the effectiveness of the constraints.

Next, the performance trade-offs caused by the length scale constraints are analysed by comparing the optimised designs against the reference. Fig. 20 shows the velocity distributions obtained from the scaling model during optimisation and validated using the RANS model. In Fig. 20a, the velocity distribution between the optimisation and RANS models is generally consistent, but as noted earlier, the scaling model overestimates boundary effects, leading to underestimation of near-wall velocities and slightly higher central velocities. At channel branches, the optimised model predicts near-zero velocities, while the RANS model shows significantly higher velocities. In Fig. 20b, after applying the constraints, the channel branches widen, mitigating the exaggerated boundary effect in the scaling model and producing velocity distributions closer to the RANS model. Among the three cases, Fig. 20b achieves the best
<figure><img src="image_22.png" /><figcaption>Fig. 20: Velocity comparison for different designs.</figcaption></figure>

agreement between the scaling and RANS results. However, in Fig. 20c, the narrower channel widths amplify boundary inaccuracies, causing the greatest deviation between the scaling and RANS models.

Fig. 21 presents the temperature distributions obtained from the RANS verification model for designs with different length scale constraints, along with the corresponding relative errors in temperature distributions compared to the optimisation model. The relative error is calculated using Eq. 32. From the figure, it is evident that the relative errors align well with the earlier discussions on model accuracy: the errors are primarily concentrated in the mould regions between the cooling channels and the polymer part, whereas the temperature errors within the part itself remain relatively low. This
<figure><img src="image_23.png" /><figcaption>**Fig. 21:** Temperature and relative error distributions of the optimised design with varying length scale constraints.</figcaption></figure>

demonstrates the effectiveness of the scaling strategy in capturing the temperature distributions to an acceptable accuracy.

It is worth emphasising that the introduction of length scale constraints has not negatively impacted the model's accuracy. In fact, the designs with length scale constraints exhibit overall smaller relative errors compared to the reference design. This may be attributed to the constraints limiting the channel diameter, which alters the convective effects within and at the boundaries of the channels. Consequently, the scaling model benefits from these constraints. Nonetheless, the key takeaway is that the
<figure><img src="image_24.png" /><figcaption>Fig. 22: Optimised designs with different length scale constraints.</figcaption></figure>

<table><caption>Table 10: Optimised designs with different length scale constraints.</caption><thead><tr><th></th><th>Design</th><th>Without constraint</th><th>R<sub>max</sub> = 2d<sub>c</sub></th><th>R<sub>max</sub> = d<sub>c</sub></th></tr></thead><tbody><tr><td rowspan="5">RANS</td><td>Bulk temperature [K]</td><td>371.84</td><td>372.20</td><td>372.09</td></tr><tr><td>Pressure drop [Pa]</td><td>3807.31</td><td>3852.16</td><td>6208.43</td></tr><tr><td>Profile temperature [K]</td><td>328.72</td><td>329.02</td><td>328.92</td></tr><tr><td>Thermal uniformity [-]</td><td>1.28</td><td>1.47</td><td>1.64</td></tr><tr><td>Volume fraction [%]</td><td>10.08</td><td>10.02</td><td>7.69</td></tr><tr><td rowspan="4">Optimisation</td><td>Bulk temperature [K]</td><td>369.22</td><td>369.75</td><td>369.65</td></tr><tr><td>Pressure drop [Pa]</td><td>3.00</td><td>2.99</td><td>5.48</td></tr><tr><td>Profile temperature [K]</td><td>326.11</td><td>326.65</td><td>326.55</td></tr><tr><td>Thermal uniformity [-]</td><td>0.78</td><td>0.83</td><td>0.87</td></tr></tbody></table>

inclusion of length scale constraints does not introduce additional errors, underscoring the reliability of the optimisation results.

Fig. 22 presents the temperature deviation profiles, calculated as the absolute difference between local and mean surface temperatures. Overall, the deviations are comparable across all cases, suggesting that the additional length scale constraints do not significantly compromise temperature uniformity. However, in regions such as C4-C5, the designs with constraints exhibit increased temperature deviations compared to the unconstrained design. This is particularly evident for the stricter constraint case ($R_{max} = R_{min} = d_c$), where deviations are higher. The observed increase is likely due to the restricted channel dimensions in these regions, which limit the ability to locally enhance convection by reducing channel sizes. In contrast, the unconstrained design benefits from the flexibility to locally adjust channel dimensions, thereby improving temperature control in critical regions such as C5.
For a quantitative comparison, Tab. 10 summarises the performance metrics of each design under both the optimisation and RANS models. The volume-averaged temperature errors between the two models are small, but larger discrepancies are observed in boundary temperatures and their standard deviations. These deviations arise partly from geometric extraction errors during optimisation and partly from the scaling model's limited accuracy in turbulent boundary simulation. Nonetheless, the errors remain within an acceptable range. Comparing the optimised designs, the reference design achieves the lowest average temperature and standard deviation, indicating that the additional length scale constraints introduce some performance trade-offs. However, these trade-offs are minimal and acceptable considering the improved manufacturability. When comparing the two constrained designs, the uniform-width design achieves lower temperatures and similar uniformity compared to the design with a maximum diameter of $2d_c$. However, the pressure drop for the uniform-width design is significantly higher—almost double—due to the narrower channel dimensions.

## 5 Conclusion

This study has presented a scaling strategy designed for the topology optimisation of cooling channels in injection moulds, effectively reducing the challenges of heat transfer under turbulent flow conditions by approximating it using a laminar flow model. By embedding the scaling model into a conjugate heat transfer framework, the approach strikes a balance between computational efficiency and physical accuracy, overcoming the inherent limitations of the traditional Darcy potential flow model. Verification case studies confirm that the scaling strategy serves as a reliable approximation for tackling industrial problems involving turbulent flow conditions. While discrepancies between the scaling model and the high-fidelity RANS model are observed, particularly near boundaries, the optimisation trends derived from the scaling model remain consistent. The RANS model verification further substantiates the effectiveness of the proposed designs when evaluated in a realistic turbulent flow context.

This work establishes a comprehensive and scalable methodology for the automated conceptual design of CC channels using topology optimisation. The methodology incorporates multiple critical constraints to address both performance and manufacturability:

*   **Temperature Uniformity Constraints:** The inclusion of temperature uniformity constraints led to designs with significantly reduced surface temperature deviations, particularly in thermally sensitive regions such as corners and protrusions, demonstrating the framework's ability to enhance thermal control.

*   **Length Scale Constraints:** Imposing minimum and maximum feature size constraints ensured the manufacturability of optimised designs while maintaining acceptable performance trade-offs. These constraints successfully regulated local structural features, enhancing the feasibility of practical implementation.

In examining the discrepancies between the scaling and RANS models, the study identifies two key sources of error: geometric extraction inaccuracies during optimisation and the scaling model's simplified boundary treatment, which underestimates
convective effects near walls. These errors, while noticeable, do not compromise the overarching trends or the reliability of the proposed framework. Future refinements, such as advanced boundary modelling and improved scaling methods, are expected to bridge the gap between models and further enhance optimisation accuracy.

Looking ahead, this methodology offers considerable potential for extension and refinement. Future work will focus on expanding the framework to three-dimensional cooling channel optimisation, transient heat transfer scenarios, and multi-objective optimisation. Additionally, incorporating advanced manufacturing constraints and leveraging high-fidelity turbulence models will enable more accurate and practical solutions for industrial applications.

**Acknowledgments.** The work by YS and JA was directly funded by Centre Technique Industriel de la Plasturgie et des Composites (CTIPC) through an industrial research collaboration agreement.

## Declarations

**Funding.** The work by YS and JA was directly funded by Centre Technique Industriel de la Plasturgie et des Composites (CTIPC) through an industrial research collaboration agreement.

**Replication of results / data availability.** A COMSOL file generating the results in Fig. 9 is available on GitHub: https://github.com/sdu-multiphysics/conformalcooling

**Author contributions.** YS: Software, Validation, Formal analysis, investigation, Writing - Original Draft, Visualisation. AGD: Methodology, Formal analysis, Writing - Original Draft, Supervision. AA: Conceptualisation, Methodology, Software, Writing - Review & Editing, Supervision, Project administration. RLG: Conceptualisation, Writing - Review & Editing, Supervision, Project administration, Funding acquisition. JA: Conceptualisation, Methodology, Software, Formal Analysis, Investigation, Resources, Writing - Review & Editing, Supervision, Project administration, Funding acquisition.

**Conflict of interest.** AGH, AA, and RLG have a commercial interest in the developed research through their employee CTIPC. YS and JA were funded by CTIPC through an industrial research collaboration agreement. JA is the Lead Guest Editor of the Topical Collection to which this manuscript is submitted.

**Ethics approval and consent.** Not relevant at all.

## Appendix A Details of Darcy model optimisation

Zhao et al (2018) proposed to approximate the turbulent flow as a inviscid flow using the Darcy potential flow model, which is governed by the following equation:

$$ \nabla \cdot \left( -\frac{\kappa}{\mu} \nabla p \right) = 0 \qquad (A1) $$
where $\kappa$ denotes the permeability of the material, and by assigning different permeability values to the solid and fluid regions, it is possible to ensure that impermeable and permeable conditions are imposed in the solid and fluid regions, respectively. The above can be expressed as:

$$ \kappa(\mathbf{x}) = \begin{cases} \kappa_f & \text{if } \mathbf{x} \in \Omega_f \\ \kappa_s & \text{if } \mathbf{x} \in \Omega_s \end{cases} \qquad (A2) $$

where $\kappa_s$ represents the permeability of the solid, $\kappa_f$ represents the permeability of the fluid. Consistent with the work of Zhao et al (2018), the present work makes $\kappa_s = 2.5 \cdot 10^{-11} \text{m}^2$, $\kappa_f = 2.5 \cdot 10^{-5} \text{m}^2$.

Notably, it is observed that during the optimisation process using Darcy's model, for this problem, the traditional RAMP interpolation is not able to obtain an optimised design with clear boundaries, so this study makes a small adjustment to the traditional interpolation so that the penalty curves can further penalise the intermediate variables, thus ensuring a clear configuration. The modified inverse interpolation function is expressed as follows:

$$ \kappa(\bar{\gamma}) = \begin{cases} \kappa_f + (\kappa_s - \kappa_f) \frac{1-\bar{\gamma}}{1+q_{\kappa}\bar{\gamma}} & \text{if } q_{\kappa} > 0 \\ \kappa_s - (\kappa_s - \kappa_f) \frac{\bar{\gamma}}{1-q_{\kappa}(1-\bar{\gamma})} & \text{if } q_{\kappa} \le 0 \end{cases} \qquad (A3) $$

where $q_\kappa$ controls the curvature of the permeability interpolation. The parameters for the continuation strategies are detailed in Tab. A1. Additionally, to achieve a well-optimised design with the Darcy model, the maximum number of iterations per continuation step was set to 100, significantly higher than the 25 iterations required for the scaling model.

<table><caption>**Table A1:** Continuation parameters for optimisation using Darcy model.</caption><thead><tr><th>Step</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th></tr></thead><tbody><tr><td>β</td><td>4</td><td>8</td><td>16</td><td>32</td><td>48</td><td>64</td><td>64</td></tr><tr><td>q<sub>κ</sub></td><td>5</td><td>3</td><td>1</td><td>0</td><td>-1</td><td>-3</td><td>-5</td></tr><tr><td>q<sub>k</sub></td><td>10<sup>4</sup></td><td>10<sup>4</sup></td><td>10<sup>4</sup></td><td>10<sup>4</sup></td><td>10<sup>4</sup></td><td>10<sup>4</sup></td><td>10<sup>4</sup></td></tr></tbody></table>

## References

Agazzi A, Sobotka V, Legoff R, et al (2013) Optimal cooling design in injection moulding process-a new approach based on morphological surfaces. Applied Thermal Engineering 52:170–178. https://doi.org/10.1016/j.applthermaleng.2012.11.019

Alexandersen J (2022) A detailed introduction to density-based topology optimisation of fluid flow problems with implementation in matlab. Structural and Multidisciplinary Optimization 66(1):12. https://doi.org/10.1007/s00158-022-03420-9, URL https://doi.org/10.1007/s00158-022-03420-9
Alexandersen J, Andreasen CS (2020) A review of topology optimisation for fluid-based problems. *Fluids* 5. https://doi.org/10.3390/fluids5010029

Alexandersen J, Aage N, Andreasen CS, et al (2014) Topology optimisation for natural convection problems. *International Journal for Numerical Methods in Fluids* 76(10):699–721. https://doi.org/10.1002/fld.3954

Alexandersen J, Sigmund O, Meyer KE, et al (2018) Design of passive coolers for light-emitting diode lamps using topology optimisation. *International Journal of Heat and Mass Transfer* 122:138–149. https://doi.org/10.1016/j.ijheatmasstransfer.2018.01.103

Andronov V, Pitrmuc Z, Zajic J, et al (2024) Conformal cooling as a support tool for eliminating local defects in high-pressure die casting series production. *Progress in Additive Manufacturing* https://doi.org/10.1007/s40964-024-00721-x

Borrvall T, Petersson J (2003) Topology optimization of fluids in stokes flow. *International Journal for Numerical Methods in Fluids* 41(1):77–107. https://doi.org/10.1002/fld.426

Coffin P, Maute K (2016) Level set topology optimization of cooling and heating devices using a simplified convection model. *Structural and Multidisciplinary Optimization* 53:985–1003. https://doi.org/10.1007/s00158-015-1343-8

Dbouk T (2017) A review about the engineering design of optimal heat transfer systems using topology optimization. *Applied Thermal Engineering* 112:841–854. https://doi.org/10.1016/j.applthermaleng.2016.10.134

Dede EM (2009) Multiphysics topology optimization of heat transfer and fluid flow systems. In: proceedings of the COMSOL Users Conference

Dienemann R, Schewe F, Elham A (2022) Industrial application of topology optimization for forced convection based on Darcy flow. *Structural and Multidisciplinary Optimization* 65(9):265. https://doi.org/10.1007/s00158-022-03328-4

Dilgen SB, Dilgen CB, Fuhrman DR, et al (2018) Density based topology optimization of turbulent flow heat transfer systems. *Structural and Multidisciplinary Optimization* 57:1905–1918. https://doi.org/10.1007/s00158-018-1967-6

Dimla D, Camilotto M, Miani F (2005) Design and optimisation of conformal cooling channels in injection moulding tools. *Journal of Materials Processing Technology* 164-165:1294–1300. https://doi.org/10.1016/j.jmatprotec.2005.02.162

Fawaz A, Hua Y, Corre SL, et al (2022) Topology optimization of heat exchangers: A review. *Energy* 252. https://doi.org/10.1016/j.energy.2022.124053
Feng S, Kamat AM, Pei Y (2021) Design and fabrication of conformal cooling channels in molds: Review and progress updates. https://doi.org/10.1016/j.ijheatmasstransfer.2021.121082

Geng D, Zhou M (2023) Hydraulic pressure control in topology optimization of cooling channels with Darcy flow model. Structural and Multidisciplinary Optimization 66(6):128. https://doi.org/10.1007/s00158-023-03575-z

Gersborg-Hansen A, Sigmund O, Haber RB (2005) Topology optimization of channel flow problems. Structural and Multidisciplinary Optimization 30(3):181-192. https://doi.org/10.1007/s00158-004-0508-7

Hanzlik J, Vanek J, Pata V, et al (2024) The impact of surface roughness on conformal cooling channels for injection molding. Materials 17(11). https://doi.org/10.3390/ma17112477, URL https://www.mdpi.com/1996-1944/17/11/2477

Jahan S, Wu T, Shin Y, et al (2019) Thermo-fluid topology optimization and experimental study of conformal cooling channels for 3d printed plastic injection molds. vol 34. Elsevier B.V., pp 631-639, https://doi.org/10.1016/j.promfg.2019.06.120

Jahan SA, Wu T, Zhang Y, et al (2016) Implementation of conformal cooling & topology optimization in 3d printed stainless steel porous structure injection molds. vol 5. Elsevier B.V., pp 901-915, https://doi.org/10.1016/j.promfg.2016.08.077

Kambampati S, Kim HA (2020) Level set topology optimization of cooling channels using the darcy flow model. Structural and Multidisciplinary Optimization 61:1345-1361. https://doi.org/10.1007/s00158-019-02482-6

Kanbur BB, Zhou Y, Shen S, et al (2022) Metal additive manufacturing of plastic injection molds with conformal cooling channels. Polymers 14(3). https://doi.org/10.3390/polym14030424, URL https://www.mdpi.com/2073-4360/14/3/424

Karakoc C, Dizdar KC, Dispinar D (2022) Investigation of effect of conformal cooling inserts in high-pressure die casting of alsi9cu3. The International Journal of Advanced Manufacturing Technology 121(11):7311-7323. https://doi.org/10.1007/s00170-022-09808-7

Kariminejad M, McAfee M, Kadivar M, et al (2024) Sensorised metal am injection mould tools for in-process monitoring of cooling performance with conventional and conformal cooling channel designs. Journal of Manufacturing Processes 116:25-39. https://doi.org/10.1016/j.jmapro.2024.02.021

Khan M, Afaq SK, Khan NU, et al (2014) Cycle time reduction in injection molding process by selection of robust cooling channel design. ISRN Mechanical Engineering 2014. https://doi.org/10.1155/2014/968484
Lazarov BS, Sigmund O (2011) Filters in topology optimization based on helmholtz-type differential equations. International journal for numerical methods in engineering 86(6):765-781

Lazarov BS, Wang F (2017) Maximum length scale in density based topology optimization. Computer Methods in Applied Mechanics and Engineering 318:826-844.
https://doi.org/10.1016/j.cma.2017.02.018

Li Q, Steven GP, Querin OM, et al (1999) Shape and topology design for heat conduction by Evolutionary Structural Optimization. International Journal of Heat and Mass Transfer 42(17):3361-3371. https://doi.org/10.1016/S0017-9310(99)00008-3

Li Z, Wang X, Gu J, et al (2018) Topology optimization for the design of conformal cooling system in thin-wall injection molding based on bem. International Journal of Advanced Manufacturing Technology 94:1041-1059. https://doi.org/10.1007/s00170-017-0901-1

Marc-Étienne LG, Marjan MZ, Vincent R, et al (2024) Additively manufactured conformal cooling channels through topology optimization. Structural and Multidisciplinary Optimization 67. https://doi.org/10.1007/s00158-024-03846-3

Navah F, Étienne Lamarche-Gagnon M, Ilinca F (2024) Thermofluid topology optimization for cooling channel design. Applied Thermal Engineering 236. https://doi.org/10.1016/j.applthermaleng.2023.121317

Nguyen VT, Minh PS, Uyen TMT, et al (2023) Conformal cooling channel design for improving temperature distribution on the cavity surface in the injection molding process. Polymers 15(13). https://doi.org/10.3390/polym15132793, URL https://www.mdpi.com/2073-4360/15/13/2793

Pan S, Yu M, Li H, et al (2022) An integrated two-step strategy for an optimal design of liquid-cooled channel layout based on thearis density approach. Structural and Multidisciplinary Optimization 65(8):221. https://doi.org/10.1007/s00158-022-03315-9

Park HS, Dang XP (2017) Development of a smart plastic injection mold with conformal cooling channels. Procedia Manufacturing 10:48-59. https://doi.org/10.1016/j.promfg.2017.07.020, 45th SME North American Manufacturing Research Conference, NAMRC 45, LA, USA

Park HS, Dang XP, Nguyen DS, et al (2020) Design of advanced injection mold to increase cooling efficiency. International Journal of Precision Engineering and Manufacturing-Green Technology 7(2):319-328. https://doi.org/10.1007/s40684-019-00041-4
Sachs E, Wylonis E, Allen S, et al (2000) Production of injection molding tooling with conformal cooling channels using the three dimensional printing process. *Polymer Engineering & Science* 40(5):1232–1247. https://doi.org/https://doi.org/10.1002/pen.11251

Santana A, Afonso P, Zanin A, et al (2017) Costing models for capacity optimization in industry 4.0: Trade-off between used capacity and operational efficiency. *Procedia Manufacturing* 13:1183–1190. https://doi.org/https://doi.org/10.1016/j.promfg.2017.09.193, manufacturing Engineering Society International Conference 2017, MESIC 2017, 28-30 June 2017, Vigo (Pontevedra), Spain

Schlichting H, Gersten K (2017) Boundary-Layer Theory. Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/978-3-662-52919-5

Stolpe M, Svanberg K (2001) An alternative interpolation scheme for minimum compliance topology optimization. *Structural and Multidisciplinary Optimization* 22(2):116–124. https://doi.org/10.1007/s001580100129, URL https://doi.org/10.1007/s001580100129

Storti BA, Sobotka V (2024) A numerical framework for three-dimensional optimization of cooling channels in thermoplastic printed molds. *Applied Thermal Engineering* 238. https://doi.org/10.1016/j.applthermaleng.2023.121988

Sun Y, Yao S, Alexandersen J (2024) Topography optimisation using a reduced-dimensional model for transient conjugate heat transfer between fluid channels and solid plates with volumetric heat source. *Structural and Multidisciplinary Optimization* 67(4):45. https://doi.org/10.1007/s00158-024-03760-8

Svanberg K (1987) The method of moving asymptotes—a new method for structural optimization. *International Journal for Numerical Methods in Engineering* 24(2):359–373. https://doi.org/10.1002/nme.1620240207

Svanberg K (2002) A Class of Globally Convergent Optimization Methods Based on Conservative Convex Separable Approximations. *SIAM Journal on Optimization* 12(2):555–573. https://doi.org/10.1137/S1052623499362822

Tomasoni D, Colosio S, Giorleo L, et al (2020) Design for additive manufacturing: Thermoforming mold optimization via conformal cooling channel technology. *Procedia Manufacturing* 47:1117–1122. https://doi.org/https://doi.org/10.1016/j.promfg.2020.04.128, 23rd International Conference on Material Forming

Torres-Alba A, Mercado-Colmenero JM, Caballero-Garcia JdD, et al (2023) Application of new conformal cooling layouts to the green injection molding of complex slender polymeric parts with high dimensional specifications. *Polymers* 15(3). https://doi.org/10.3390/polym15030558
Wang F, Lazarov BS, Sigmund O (2011) On projection methods, convergence and robust formulations in topology optimization. Structural and Multidisciplinary Optimization 43(6):767-784. https://doi.org/10.1007/s00158-010-0602-y

Wang J, Cheng F, Wei Y (2024a) Research on automatic generation technology of double helix-shaped conformal cooling channels for injection mold. The International Journal of Advanced Manufacturing Technology 134(11):5935-5949. https://doi.org/10.1007/s00170-024-14545-0

Wang ML, Zheng LJ, Kang HW (2024b) 3-Dimensional conformal cooling channel design: Origami-inspired topology optimization approach. Applied Thermal Engineering 242:122526. https://doi.org/10.1016/j.applthermaleng.2024.122526

Wei Z, Wu J, Shi N, et al (2020) Review of conformal cooling system design and additive manufacturing for injection molds. Mathematical Biosciences and Engineering 17:5414-5431. https://doi.org/10.3934/MBE.2020292

Wu T, Tovar A (2018) Detc2018-85511 design for additive manufacturing of conformal cooling channels using thermal-fluid topology optimization and application in injection molds

Xu X, Sachs E, Allen S (2001) The design of conformal cooling channels in injection molding tooling. Polymer Engineering & Science 41:1265-1279

Yang K, Fernandez E, Niu C, et al (2019) Note on spatial gradient operators and gradient-based minimum length constraints in SIMP topology optimization. Structural and Multidisciplinary Optimization 60(1):393-400. https://doi.org/10.1007/s00158-019-02269-9

Yoon GH (2010) Topological design of heat dissipating structure with forced convective heat transfer. Journal of Mechanical Science and Technology 24(6):1225-1233

Yu D, Rachik M, Blaise A, et al (2024) Topology optimization designed twisted conformal cooling channel for additive-manufactured hot-stamping tool. Applied Thermal Engineering p 124784. https://doi.org/10.1016/j.applthermaleng.2024.124784

Yu M, Ruan S, Wang X, et al (2019) Topology optimization of thermal-fluid problem using the MMC-based approach. Structural and Multidisciplinary Optimization 60(1):151-165. https://doi.org/10.1007/s00158-019-02206-w

Zhao J, Zhang M, Zhu Y, et al (2021) Topology optimization of planar cooling channels using a three-layer thermofluid model in fully developed laminar flow problems. Structural and Multidisciplinary Optimization 63:2789-2809. https://doi.org/10.1007/s00158-021-02842-1
Zhao X, Zhou M, Sigmund O, et al (2018) A “poor man’s approach” to topology optimization of cooling channels based on a Darcy flow model. International Journal of Heat and Mass Transfer 116:1108–1123. https://doi.org/10.1016/j.ijheatmasstransfer.2017.09.090

Zhou M, Lazarov BS, Wang F, et al (2015) Minimum length scale in topology optimization by geometric constraints. Computer Methods in Applied Mechanics and Engineering 293:266–282. https://doi.org/10.1016/j.cma.2015.05.003

Zhou M, Alexandersen J, Sigmund O, et al (2016) Industrial application of topology optimization for combined conductive and convective heat transfer problems. Structural and Multidisciplinary Optimization 54(4):1045–1060. https://doi.org/10.1007/s00158-016-1433-2