Topology optimization falls within the realm of non-parametric optimization methodologies, with its primary objective being the realization of designs by tackling a numerically discretized representation of the design challenge.
This study focuses on employing topology optimization to simultaneously enhance heat transfer efficiency and reduce the generated EMI by minimizing heat compliance and maximizing the VSWR within a structural layout. Typically, topology optimization involves discretizing design variables (\({\mathbf{x}}\)) using finite element methods. In our research, we utilized linear triangular elements for electromagnetic simulations and bilinear rectangular structured mesh for heat transfer and convection simulations. A value of 1 signifies a solid element, whereas 0 represents a void element, as depicted in Eq. (1).
$${\mathbf{x}} = \left\{ {\begin{array}{ll} 1 & \quad {Design \; materials} \\ 0 & \quad {Void} \\ \end{array} } \right.$$
(1)
In this section, we delve into the mathematical modeling of maximizing heat transfer while addressing heat compliance. Following that, we’ll explore the mathematical modeling of electromagnetic propagation phenomena. Finally, we will discuss the optimization algorithm employed in our study.
Mathematical modelling for maximizing heat conductivity for solid structures
For thermal dissipater optimization, the problem is reciprocal. Not optimizing original design domain, with ROI as heat source and boundaries as heat flow, the approach minimizes heat compliance. It uses distributed heat source and adiabatic boundaries, except ROI at zero. Enhancing conductivity minimizes thermal energy by reducing compliance. The analysis considers heat conduction and transfer from ROI to boundaries with convection. Only steady-state is examined. The equation is:
$$ {\mathbf{Q}} = {\mathbf{KT}} $$
(2)
where, \({\mathbf{Q}}\) is the nodal heat load, \({\mathbf{T}}\) nodal temperature, and is the \({\mathbf{K}}\) is heat conductivity matrix. Steady-state heat conduction topology optimization aims to identify an optimal structural layout that minimizes heat generation within a specified solid material volume. This optimization involves material distribution that seeks to lower the structure’s temperature, thus enhancing thermal dissipation. As such, minimizing the heat potential capacity (i.e., heat compliance \(C\)) has chosen as the objective function (Eq. 3), with design variables determined by element state variables, i.e., density variables that indicate the presence (1) or absence (0) of materials.
$$ C = {\mathbf{Q}}^{T} {\mathbf{T}} $$
(3)
Moreover, in this research, we are considering the heat convection. The heat convection governing equation is the relation between the heat convection load \({\mathbf{Q}}_{Conv}\) and the relative difference of the ambient temperature above the surface of the solid design domain \({\mathbf{T}}_{\infty }\) and the solid design domain \({\mathbf{T}}_{s}\).
$$ {\mathbf{Q}}_{Conv} = A{\mathbf{k}}_{Conv} ({\mathbf{T}}_{s} – {\mathbf{T}}_{\infty } ) $$
(4)
where \({\mathbf{k}}_{Conv}\) is the heat convection coefficient, and \(A\) is the solid surface that is in contact with the void that is filled with air. Therefore, at the free surfaces of the solid materials \(\Omega_{s}\), the convection is taking place, so the element on the solid’s boundary is needed to be updated during the optimization process using the term \({\mathbf{K}}_{Conv}^{i}\) such that:
$$ {\mathbf{K}}_{Conv}^{i} = \int_{{\Omega_{s} }} {{\mathbf{N}}^{T} {\mathbf{Nk}}_{Conv} d\Omega_{s} } $$
(5)
where \({\mathbf{N}}\) is the shape function matrix for the heat distribution field on the surface. Moreover, the corresponding heat flux vector load update is formulated as follows:
$$ {\mathbf{Q}}_{Conv}^{i} = \int_{{\Omega_{s} }} {{\mathbf{N}}^{T} {\mathbf{k}}_{Conv} {\mathbf{T}}_{\infty } d\Omega_{s} } $$
(6)
In this research, the heat conductivity maximization objective function is taking the form:
$$ \begin{aligned} & find\quad {\mathbf{x}} \\ & \min :\;C({\mathbf{x}}) = {\mathbf{Q}}^{T} {\mathbf{T}} \\ & s.t.\;\left\{ \begin{gathered} {\mathbf{K}}({\mathbf{x}}){\mathbf{T}} = {\mathbf{Q}} \hfill \\ \int_{\Omega } {{\mathbf{x}}d\Omega } \le v,\;{\mathbf{x}} = 0|1,\;\forall {\mathbf{x}} \in \Omega \hfill \\ \end{gathered} \right. \\ \end{aligned} $$
(7)
where \({\mathbf{x}}\) is the binary design variable controlling whether ith element is solid or void ((0) or (1)). \(v\) is the volume fraction to be satisfied for the design domain \(\Omega\).
In order to extend the applicability of sensitivity analysis for multiscale analysis, we’ve expressed all equations using Newton notation to depict the derivatives.
The derivative of heat compliance (i.e., \(\dot{C}\)) is taking the general form as:
$$ \dot{C} = {\dot{\mathbf{Q}}}^{T} {\mathbf{T}} + {\mathbf{Q}}^{T} {\dot{\mathbf{T}}} $$
(8)
where \({\dot{\mathbf{Q}}}\), \({\dot{\mathbf{T}}}\), are the first order derivative of the heat load and temperature with considering Newton notation. \({\dot{\mathbf{Q}}}\) can be obtained by deriving equilibrium equation for heat conduction, such that:
$$ {\dot{\mathbf{Q}}} = {\dot{\mathbf{K}}\mathbf{T}} + {\mathbf{K}} {\dot{\mathbf{T}}} $$
(9)
So, the term \({\dot{\mathbf{T}}}\) can be expressed as:
$$ {\dot{\mathbf{T}}} = {\mathbf{K}}^{ – 1} \left( {{\dot{\mathbf{Q}}} – {\dot{\mathbf{K}}\mathbf{T}}} \right) $$
(10)
Consequently, the sensitivity for a general compliance minimization problem can be formulated as follows:
$$ {\mathbf{C}} = 2{\dot{\mathbf{Q}}}^{T} {\mathbf{T}} – {\mathbf{T}}^{T} {\dot{\mathbf{K}}\mathbf{T}} $$
(11)
In case of fixed loading condition, the first term of Eq. (19) (i.e.,\({\dot{\mathbf{Q}}}\)) will vanish leaving the sensitivity to become:
$$ \dot{C} = – {\mathbf{T}}^{T} {\dot{\mathbf{K}}\mathbf{T}} $$
(12)
It is apparent from Eq. (12) that the minimization of heat compliance is exhibiting the characteristics of a self-adjoint problem.
EMI objective function
EMI is electromagnetic disturbance from devices, disrupting function and possibly exceeding human exposure limits for specific absorption rate (SAR). SAR measures human body energy absorption from electromagnetic fields, like mobile phones. EMI is simulated using waveguide propagation and frequency behavior. Engineers refer to the source as RF circuit with fixed wavelength. In simulations, the metallic structure (e.g., heat sink) serves as principal RF circuit (Fig. 1). When leaked frequency aligns with heat sink inductance and air capacitance, it becomes an unintended antenna. Evaluating its emissivity/absorptivity is vital. Antenna wave propagation is intricate due to factors like medium conductivity, obstacles, reflections, and wave traits. Precise electromagnetic wave measurement poses challenges. Near-field and far-field scanning, simulations, common methods, have limitations in capturing wave behavior.
The equivalent electromagnetic circuit of heat sink.
When the leaked frequency aligns with the resonance frequencies of the metallic heat sink’s inductance and the surrounding air’s parasitic capacitance, the heat sink may inadvertently function as an antenna. In such cases, it becomes crucial to assess the metallic structure’s emissivity/absorptivity for electromagnetic waves. Comprehending wave propagation from an antenna entails considering various factors, including medium conductivity, obstacles, reflections, and wave characteristics. However, accurately measuring electromagnetic waves poses challenges, with techniques like near-field and far-field scanning and simulations commonly employed, albeit each with its limitations in capturing wave behavior. To facilitate electromagnetic wave analysis within and beyond the metallic object, researchers have introduced various parameters, including admittance, impedance, scattering transfer, scattering, and hybrid parameters. These parameters play crucial roles in RF component design and evaluation, with a notable emphasis on scattering parameters (\(S_{ij}\))36. \(S_{ij}\) have emerged as the most popular and widely utilized parameters for waveguide RF assessment, owing to their effectiveness in characterizing the behavior of electromagnetic waves within the metallic structure and their interaction with surrounding elements. In RF circuits, it is often described using wave quantities which are the incident wave (Wx) and the reflected wave (Wr) within the RF circuit within the Euler operator (\(e^{i\theta }\)).
$$ {\text{Wx}} = E_{input} e^{i\theta } $$
(13)
$$ {\text{Wr}} = E_{reflected} e^{i\theta } $$
(14)
where \(E_{input}\) and \(E_{reflected}\) are the electric potential of the incident and reflected wave. The ratio of the reflected wave is the complex reflection coefficient (Γ), which denominates the relation between the reflected and the incident waves:
$$ \Gamma = \frac{Wr}{{Wx}} $$
(15)
Therefore, by connecting the device under testing (DUT), to the VNA device, with single port or two ports; it is easy to evaluate the ability of DUT to omit the targeted RF practically. Moreover, using numerical method for evaluating \(S_{ij}\) is highly recommended.
Numerically, the scattering parameters are evaluated using finite element as the following:
$$ {\mathbf{S}}_{ij} = \frac{1}{2}{\mathbf{E}}^{T} {\mathbf{K}}_{A} {\mathbf{E}} $$
(16)
where \({\mathbf{E}}\) is the electric filed, and \({\mathbf{K}}_{A}\) is the global antenna finite-element electromagnetic stiffness matrix. Reducing the potential for the optimized topology heat sink to act as an antenna for leaked electromagnetic power has been achieved through EMI minimization efforts. This is accomplished by evaluating the Voltage Standing Wave Ratio (VSWR)37. VSWR is taking the form:
$$ VSWR = \frac{{1 + \left( {\left| {E_{input} } \right|/\left| {E_{reflected} } \right|} \right)e^{i\theta } }}{{1 – \left( {\left| {E_{input} } \right|/\left| {E_{reflected} } \right|} \right)e^{i\theta } }} $$
(17)
where \(\left| {E_{input} } \right|\) is the absolute supplied electric potential on the antenna, \(\left| {E_{reflected} } \right|\) is the reflected electric potential to the source, and \(\theta\) is the phase difference between the applied and reflected electric potentials. As such, the VSWR can be rewritten in terms of the \(S_{ij}\) as:
$$ VSWR = \frac{{1 + |S_{11} |}}{{1 – |S_{11} |}} $$
(18)
Realizing high heat conductivity with low EMI structure using the MASTO method
MASTO improves topology optimization’s efficiency, speed, and quality. It combines YUKI’s hybrid optimization and gradient descent. MASTO defines binary design area by YUKI. It evaluates and splits area into local and exploration regions. It divides variables into sets without adaptive functions. It increases global solutions38,39.
Various optimization methods are employed to attain optimal solutions and minimize suboptimal outcomes. In topology optimization, vital for sizable problems, heuristic methods are crucial. They trim computational costs and enhance optimization stability. Yet, optimal design hinges on factors like search increment choices and initial algorithm point. The latter is vital, especially in intricate functions like maximizing heat conductivity in point-to-volume approaches14,34,40. Metaheuristic algorithms excel in general-purpose optimization, not requiring gradients for practical deployment. Their multipoint search distinguishes them from gradient descent, enabling extensive solution space exploration and global extrema identification. Metaheuristics’ multipoint strategy benefits initial point selection in topology optimization, yielding optimal designs with reduced computational costs34. In this study, our aim is to attain multiphysics optimization of structures that exhibit both elevated heat conductivity and minimal EMI concurrently. To achieve this dual objective, we adopt an optimization methodology centred on the reduction of the ratio between heat resistance and the VSWR inherent in the design configuration, as illustrated in Eq. (19).
$$ \begin{aligned} & find\;{\mathbf{x}} = [{\text{x}}_{i = 1..N} ] \\ & \min :\;f\left( {\mathbf{x}} \right) = \frac{C}{VSWR} \\ & {\text{s}} .t.\;\left\{ {\int\limits_{\Omega }^{{}} {{\mathbf{x}}d\Omega } \le V,\;{\mathbf{x}} = 0|1\forall {\mathbf{x}} \in \Omega } \right. \\ \end{aligned} $$
(19)
To enhance the efficiency of topology optimization process, we employ a method that explores multiple initial search points34,41. MASTO optimization method starts with initializing candidate solutions from the YUKI algorithm, forming initial strip lines in the design domain (Fig. 2). YUKI explores solution spaces to balance exploration—searching diverse regions—and exploitation, refining solutions in promising areas. This blend prevents local optima traps and improves solutions near high-potential regions. Synergizing these approaches expedites the search for structures with desired attributes like exceptional heat conductivity and minimal EMI. Optimization employs in-house MATLAB codes integrating partial differential equation and antenna toolboxes. MASTO maximizes a predefined objective function, focusing on the heat compliance to Voltage Standing Wave Ratio (VSWR) ratio. Its framework (Fig. 3) begins with population generation, computing optimal and mean design variables. Local search boundaries are set, defining variable ranges, using stochastic processes for exploration or focus. Exploration introduces random perturbations within boundaries to generate new solutions, while focus adjusts variables towards optima. Iterative processes continue until convergence, driven by fitness criteria evaluation. Sensitivity analysis measures responses to variable modifications. Gradient descent updates variables iteratively, aligning them with the steepest objective function descent. Heat compliance assessment tests design resilience iteratively until algorithmic completion. Supplementary materials provide additional details on the MASTO method.
Optimization process with MASTO.
The flowchart of MASTO optimization method.
