Numerical Analysis and Multi-Objective Optimization of Tesla Valve Cold Plates for Lithium-Ion Battery Thermal Management

1 Introduction

Lithium-ion batteries have been widely adopted in electric or hybrid electric vehicles in recent years, owing to their several advantages, including high energy density, long cycle life, and good safety performance [1]. However, as the demand for extended driving range in electric vehicles continues to rise, the increasing number of batteries within the pack and their increasingly compact arrangement make it difficult to effectively dissipate the heat generated during operation. This readily leads to excessively high battery temperatures, impairing their operational performance and, in severe cases, potentially triggering thermal runaway [2,3]. The ideal operating temperature range for traction batteries is 20°C to 40°C, with an optimal point around 30°C [4]. Consequently, employing a high-performance battery thermal management system is crucial to ensuring the safe, reliable, and long-term operation of traction batteries [5,6,7].

Current heat dissipation methods for traction batteries in electric vehicles primarily include air cooling, liquid cooling, phase change material (PCM) cooling, heat pipe cooling, direct refrigerant cooling, and hybrid cooling [8,9,10]. Among these, air cooling offers structural simplicity but struggles to meet basic cooling demands under the high energy density conditions of modern battery packs. PCM cooling can passively absorb heat through phase transition; for example, Harish et al. [11] demonstrated that sodium sulphate decahydrate (SSD) reduces peak pack temperature from 50°C to 34°C at 4C discharge in an 18,650 battery pack. However, PCM cooling still suffers from inherent drawbacks such as low thermal conductivity, encapsulation difficulties, large volume requirements, and poor thermal recovery after melting. Heat pipe cooling demonstrates effective heat transfer performance; however, it involves complex structures and higher costs. Hybrid cooling and direct refrigerant cooling face challenges related to high energy consumption and cost. In contrast, liquid cooling presents more significant advantages in terms of heat exchange efficiency, structural simplicity, and cost-effectiveness [12,13].

To advance research on battery thermal management, we consider the Tesla valve [14]—a unique check valve with no moving parts internally, originally patented by Nikola Tesla in 1920. Compared to conventional valves, the Tesla valve features a simpler structure, contains no moving components, facilitates large-scale manufacturing, and offers a longer service life. Its distinctive bidirectional flow characteristics enable effective flow rectification, while the multi-bend channel design generates significant turbulence, inducing boundary layer disturbance and enhancing heat transfer [15]. Leveraging these characteristics, numerous researchers have explored the application of Tesla valve flow channels in battery thermal management.

Research on the heat transfer and flow characteristics of Tesla valves forms the foundation of this work. Mohammadzadeh et al. [16] investigated the influence of Tesla valve size (from single-stage to four-stage configurations) and Reynolds number (Re) on valve performance. The results indicated that a two-stage Tesla valves layout represents the optimal choice, and the variation pattern of the diodicity (Di) with Re remains unaffected by valve size. Meanwhile, Thompson et al. [17] confirmed via 3D-Computational fluid dynamics (CFD) simulations that for multi-stage Tesla valves (MSTVs), a significant enhancement in diodicity can be achieved by increasing the number of stages and reducing the inter-valve spacing when Re > 50. Furthermore, the performance advantage of MSTVs intensifies with increasing Reynolds number. Qian et al. [18] optimized the thermal dissipation performance of a Tesla valve microchannel heat sink via numerical simulations. Their study investigated the influence of the number of Tesla valve stages, the outer radius of the curved channels, and the orientation of the Tesla valve units on the thermal and flow characteristics. Nobakht et al. [19] investigated the diode mechanism and its influencing factors in microvalves of different structures. It was found that the primary sites for the diode mechanism occur at the T-junctions and side channels; furthermore, eliminating the inlet and outlet channels significantly enhances diodicity. Monika et al. [20] numerically investigated a four-channel MSTV liquid cooling plate, finding that larger channel width/outer radius and smaller valve angle reduce battery temperature rise significantly, with minimal pressure drop penalty from radius and angle changes.

Meanwhile, Zhang et al. [21] developed a three-dimensional parametric model of a Tesla valve, with their study demonstrating that Tesla valves with high aspect ratios exhibit superior flow rectification capability. Nguyen et al. [22] systematically investigated the flow and rectification characteristics of Tesla valves under inertial flow conditions. Through numerical analysis, they revealed that the diodic effect of the Tesla valve initiates abruptly at Re ≈ 200, accompanied by laminar-turbulent transition, and pulsatile flow can effectively enhance its flow rectification performance. Truong et al. [23] systematically proposed, for the first time, a comprehensive optimization methodology for the Tesla valve. Through numerical simulations, they established inverse and proportional relationships between key geometric parameters—specifically, the angle and straight-segment length—and the Reynolds number (Re).

Multiple studies have demonstrated that multi-objective optimization is an effective approach to enhance the performance of battery thermal management systems (BTMS). Dong et al. [24] combined the NSGA-II algorithm with neural networks to conduct multi-objective optimization on BTMS. Using channel spacing, width, angle, and mass flow rate as design variables under a maximum temperature constraint, their optimization results showed that, compared with the initial model, the maximum temperature difference was reduced by 14.90%, the heat transfer coefficient was increased by 35.79%, and the pressure drop was decreased by 68.33%, effectively enhancing the heat dissipation performance. Wang et al. [25] proposed a structural optimization scheme for battery packs with bidirectional unequal spacing. This method integrated neural networks and the particle swarm optimization (PSO) algorithm to perform collaborative optimization of the system, aiming at reducing the maximum temperature rise and energy consumption. After optimization, the system energy consumption was reduced to 41.19% of the original design, effectively improving energy efficiency and heat dissipation performance. Zhang et al. [26] developed a BTMS coupled with phase change materials (PCMs) and microchannels. Multi-objective optimization was carried out using the response surface methodology (RSM) and the NSGA-II algorithm. Despite a slight increase in the core temperature, the comprehensive performance of the system was significantly enhanced: better temperature uniformity, a substantial reduction in pressure drop, and a smaller system volume. These findings consistently indicate that multi-objective optimization methods can systematically balance the conflicting performance indicators in BTMS and achieve improvements in overall performance.

While recent studies have made valuable progress in understanding the effects of individual geometric parameters of Tesla valves on cold plate performance, several critical gaps remain unaddressed. For instance, Porwal et al. [27] investigated the influence of shunt angle and stage number on diodicity, but their analysis was limited to flow characteristics without considering heat transfer performance under battery thermal management conditions. Zhao et al. [28] examined multistaged Tesla valves and established power-law correlations for pressure drop, yet the coupled effects of multiple geometric parameters (e.g., asymmetry ratio and branch channel width) on both thermal and hydraulic performance were not systematically explored. Furthermore, existing optimization studies on Tesla valve cold plates [29] have predominantly employed single-factor or orthogonal experiments without adopting surrogate-assisted multi-objective optimization frameworks, thus failing to reveal the trade-off between temperature uniformity and pressure drop. In the context of battery thermal management, a systematic analysis of the synergistic effects of shunt angle, number of unit pairs, channel asymmetry ratio, and branch channel width is still lacking.

To address this research gap, this paper will delve into the impacts of shunt angle, number of unit pairs, channel asymmetry ratio, and branch channel width on the battery thermal management system. Furthermore, centering on the multi-objective optimization requirements of the system, a comprehensive design method is proposed, and the specific process is as follows: First, parametric modeling constructs a design space incorporating key variables including shunt angle (30° to 50°), number of unit pairs (3 to 7 pairs), channel asymmetry ratio (0 to 0.5 h ), and branch channel width (2 to 4 mm); Computational fluid dynamics (CFD) simulations are then employed to obtain flow and heat transfer characteristics for various Tesla valve cold plate configurations; Finally, by integrating Kriging surrogate models with the NSGA-II genetic algorithm and targeting minimization of maximum temperature difference and pressure drop as optimization objectives, Pareto frontiers sets for cold plate performance are systematically identified.

2 Model Description

2.1 Geometric Model

This study employs the Tesla valve as the structural basis for the cold plate, with the valve system primarily comprising the main channel and the branch channel (see Fig. 1). The unidirectional flow characteristic originates from significant flow distribution disparities when fluid traverses the branch channels in forward versus reverse directions [30]. Under constant flow rate conditions, reverse flow experiences preferential fluid entry into branches due to inertia; this branch flow subsequently collides with main channel flow at downstream bifurcations, generating complex flow fields characterized by localized turbulence. These fields increase the heat transfer rate but also induce intense vortex-driven energy dissipation and elevated local pressure losses.

Figure 1: Schematic diagram of the Tesla valve structure.

The Tesla valve parameters in the present study are illustrated in Fig. 2a [31]. The inlet cross-section is square with side length D = 3 mm. The main channel length L 1 and the branch channel length L 2 are both 20 mm.

To ensure physical rationality and consistency with existing research, all geometric parameters are selected within ranges validated by published literature. Specifically, the included angle α between the main channel and branch channels is set to 30°, 35°, 40°, 45°, and 50°. This range is determined based on Reference [23]: angles smaller than 30° would lead to excessive flow resistance in the branch channels, drastically reducing forward flow efficiency; conversely, angles larger than 50° would orient the branch channels nearly perpendicular to the main flow direction, weakening collision-induced dissipation during reverse flow and impairing unidirectional flow performance. As such, the range of 30° to 50° encompasses the typical engineering operating scope for Tesla valves.

The width of the branch channels d is set to 2 mm, 3 mm, and 4 mm, corresponding to ratios of 0.67 to 1.33 relative to the main channel width D = 3 mm. As documented in Reference [19], this ratio range optimizes the balance between the inertia of branch flow and the collision effect within the main channel: excessively narrow branch channels cannot generate effective flow disturbances, whereas overly wide ones cause excessive fluid mixing between the main and branch flows, leading to a sharp increase in pressure drop.

Each unit pair consists of two adjacent upper and lower Tesla valves. Different cold plate configurations are designed with varying numbers of such unit pairs n and inter-unit spacing h, as detailed in Table 1. The number of unit pairs n is chosen from 3 to 7, as validated in Reference [17]. As n increases from 3 to 7, both pressure drop and heat transfer performance exhibit nonlinear changes, which fully capture the transition of flow disturbance from the developing stage to the fully developed stage. The inter-unit spacing h is defined as the horizontal offset between adjacent Tesla valve units; this parameter is normalized by the unit period to form an asymmetry ratio η , set within the range of 0 to 0.5. As reported in Reference [29], this ratio range induces substantial variations in the fluid distribution ratio between the main and branch channels, directly influencing the generation and dissipation of flow turbulence. It therefore fully reflects how structural asymmetry impacts flow and heat transfer characteristics.

The fixed overall dimensions of the cold plate are shown in Fig. 2b.

Figure 2: Geometric configuration parameters: (a) Tesla valve; (b) Cold plate.

2.2 Mathematical Model

2.2.1 Fluid Dynamic Model

In this work, water is selected as the coolant which is treated as an incompressible, homogeneous fluid with constant thermophysical properties. The Reynolds number (Re), a key parameter for distinguishing between its laminar and turbulent flow regimes, which is defined as: Re=ρυμ·2aba+b(1) where, ρ , υ and μ represent the coolant density (kg/m3), fluid velocity (m/s), and dynamic viscosity (Pa·s), respectively; a and b denote the length and width of the flow channel cross-section (m), respectively.

For the 3 mm × 3 mm square channel adopted in this study, the critical velocity between laminar and turbulent flow is calculated to be approximately 0.77 m/s.

In subsequent calculations, the fluid velocity within the cold plate exceeds 0.77 m/s, resulting in a Re greater than 2300, which indicates that the flow enters the transitional or turbulent regime. Considering that the Tesla valve channel contains disturbed structures such as bends and bifurcations, which can induce turbulence at Re > 2300, this study employs the Realizable k- ε model for numerical simulation.

The governing equations of the model include the incompressible continuity equation, modified Navier–Stokes equation, and four auxiliary equations for the turbulent flow energy k and turbulent velocity ε , expressed as follows [32]: ρ∇·u=0(2) ρu·∇u=∇·−PI+μ+μT∇uT+F(3) ∂ρcpT∂t+∇·ρcpu→T=∇·λ∇T+φ+S(4) ρu·∇k=∇u+μTσk∇k+Pk−ρε(5) ρu·∇ε=∇·μ+μTσε∇ε+ρC1Sε−ρC2ε2k+vε(6) μT=ρCμk2ε,Pk=μT∇u:∇uT(7) where, ρ represents the coolant density; u denotes the fluid velocity vector. P represents the fluid static pressure; I denotes the unit tensor; μ stands for the dynamic molecular viscosity; μ T corresponds to the turbulent eddy viscosity; F represents the external volume force term. t represents the time; c p denotes the constant-pressure specific heat capacity of fluid; T stands for the fluid temperature; λ corresponds to the fluid thermal conductivity; φ represents the viscous dissipation heat source term; S denotes the volumetric internal heat source term. k represents the turbulent kinetic energy; σ k denotes the turbulent Prandtl number for turbulent kinetic energy; P k stands for the production term of turbulent kinetic energy; ε corresponds to the turbulent dissipation rate. σ ε denotes the turbulent Prandtl number for turbulent dissipation rate; C 1 , C 2 and C μ stand for empirical constants of the standard k - ε turbulence model; υ corresponds to the kinematic molecular viscosity.

In this study, the external force term is set to F = 0 . The thermophysical parameters ( ρ , c p , λ , μ ) of the cooling medium are provided in Table 2. The parameters of the standard k - ε turbulence model are uniformly set as follows:

To ensure the numerical stability and convergence of the solution, all simulations are performed using the ANSYS Fluent double-precision solver with a pressure-based solution type. The Coupled algorithm is adopted for pressure-velocity coupling, and the second-order upwind scheme is employed for the spatial discretization of all governing equations to ensure high numerical accuracy.

For convergence criteria, the residual tolerance for the energy equation is set to 1 × 10−8, while the residuals for the continuity, momentum, and turbulence equations are set to 1 × 10−5. The solution is deemed converged when all residuals fall below the specified thresholds, and no further significant variations are observed in the monitored maximum temperature and pressure drop across the cold plate.

2.2.2 Thermal Model of Battery

In this study, the positive and negative electrodes of the lithium-ion cell are neglected, and the internal material of the battery is assumed to be isotropic, homogeneous, and uniformly distributed. The thermophysical properties including density, specific heat capacity, and thermal conductivity are treated as constants (independent of temperature). Based on these assumptions, the three-dimensional unsteady thermal conduction differential equation governing the heat transfer in the battery in a Cartesian coordinate system is expressed as follows: ρbcb∂Tb∂t=λx∂2Tb∂x2+λy∂2Tb∂y2+λz∂2Tb∂z2+qv(8) where, ρ b , c b , q v and T b denote the battery density (kg/m3), specific heat capacity (J/(kg·K)), volumetric heat generation rate (W/m3), and temperature (K), respectively; λ x , λ y and λ z represent the thermal conductivities of the lithium-ion cell along the x, y and z directions (W/(m·K)), respectively; and t is the time (s).

The volumetric heat generation term q v in Eq. (9) is calculated according to the battery heat generation equation proposed by Bernardi et al. [33], defined as: qv=1vbI2R−ITb∂UOCV∂Tb(9) where, V b is the battery volume (m3); I is the charge/discharge current (A); ∂ U O C V ∂ T b is the entropy coefficient, with a typical range of 1~2.8 × 10−4 V/K; and R is the internal resistance of the battery ( Ω ).

2.2.3 Performance Evaluation Indicator

In the simulation analysis, the Performance Evaluation Criterion (PEC) is ultimately selected to comprehensively evaluate the optimization performance. The specific formula is defined as [34]: PEC=Nue/Nurfe/fr13(10) Nu=hwDhλw(11) hw=QAhTs−Tw(12) where, h w represents the convective heat transfer coefficient of the fluid; D h denotes the characteristic length of the heat transfer surface; λ stands for the thermal conductivity of the fluid. Additionally; Q is the heat transfer rate; A is the convective heat transfer area; while T s and T w correspond to the average temperatures of the tube wall and the coolant during convective heat transfer, respectively.

The surface friction coefficient of the straight channel inside the cold plate is estimated by adopting the Darcy-Weisbach friction loss equation, to reflect the flow resistance within the channel, as shown in: f=ΔP·Dh2ρwυ¯2LC(13) where, the subscripts e and r denote the structure to be evaluated and the reference structure, respectively. Furthermore, f is the friction factor of the channel; Δ p is the pressure drop at the inlet and outlet; υ ¯ is the average velocity of the fluid; is L C the length of the channel.

2.3 Initial and Boundary Conditions

In this study, the Tesla valve inlet is designated as the reverse flow direction. The ambient temperature and initial temperature are both set to 20°C. The coolant inlet is defined as a mass flow inlet with a flow rate of 0.0048 kg/s and a temperature of 20°C. The coolant outlet is specified as a pressure outlet with gauge pressure set to 0 Pa. It is assumed that the internal battery materials are isotropic and homogeneously distributed, with their specific heat capacity, density, and thermal conductivity being constant. Furthermore, uniform internal heat generation within the battery is assumed. A constant heat flux of q = 5000 W/m2 [35] is applied to the bottom surface of the XY-plane to simulate the heat dissipation requirement of the battery under typical high-rate discharge conditions. This heat flux value is determined with sufficient literature support and engineering rationality: it corresponds to the average heat flux density of a typical prismatic LiFePO4 battery under 2C–3C high-rate discharge. According to references [35], the average surface heat flux density of the battery is approximately 4500–5500 W/m2 at 2C discharge. Thus, 5000 W/m2 represents a moderately high discharge condition, which can effectively verify the heat dissipation capability of the cold plate under typical high-load scenarios. The air convective heat transfer coefficient on the outside surface of the plate is set to 5 W/(m2·K) [36]. In the simulations, the effects of gravity are neglected in the numerical simulations.

The thermophysical properties of the coolant, cold plate, and battery are listed in Table 2.

2.4 Grid Independence Verification

The cold plate with the following geometric parameters was selected for the grid independence verification: α = 40 ∘ , n = 3 , η = 0.375 , d = 3 mm. Grid independence was assessed by monitoring the maximum temperature within the cold plate and the pressure drop across the cooling channel.

In this verification, 8 grid schemes with different element sizes were established. Given the stable structure of the solid domain and the complex flow characteristics of the fluid domain, a differentiated meshing strategy was adopted for the two domains: the element size of the fluid domain was set smaller than that of the solid domain, and local mesh refinement was implemented for the key structures in the channel. The detailed parameters are presented in Table 3.

As illustrated in Fig. 3, both T max and Δ p progressively stabilize as the mesh count increases. Using the results obtained with 325,142 mesh elements as the baseline (yielding a maximum temperature of 28.9545°C and a pressure drop of 1459.34 Pa) the relative errors in both T max and Δ p remain below 0.1% when the mesh count is increased to 473,625. Consequently, the mesh configuration with 325,142 elements was adopted as the baseline mesh for all subsequent simulations.

By comprehensively considering both computational accuracy and computational cost, the grid configuration with a solid domain element size of 2 mm, a fluid domain element size of 1 mm, and a total of 325,142 mesh elements was finally determined as the baseline mesh scheme for all subsequent numerical simulations.

Figure 3: Verification of grid independence.

2.5 Verification of the Model

To verify the effectiveness of the numerical simulation methodology adopted in this study, the Tesla valve model investigated by Liu et al. [37] was used for pressure drop validation, and the Tesla valve model reported by Han et al. [38] was employed to verify the temperature difference. The boundary conditions and thermophysical parameters of materials applied in the simulations strictly followed the specifications in the corresponding references.

Fig. 4a,b presents the comparative analysis between the experimentally measured pressure drop under different flow rates and the numerical results, as well as the comparison between the experimentally obtained maximum temperature under various heat fluxes and the simulated data. The figure indicates that the present numerical results are in high agreement with the two sets of experimental data, and the maximum relative errors are all below 3%.

Consequently, the established numerical model and solution methodology exhibit sufficient reliability for application to the Tesla valve cold plate configuration examined in this work.

Figure 4: Verification of the model (a) pressure drop; (b) temperature difference.

3 Analysis of Tesla Valve Cooling Performance

3.1 Orthogonal Experimental Design

Structural factors, including shunt angle α , number of unit pairs n , channel asymmetry ratio η , and branch channel width d , were selected, with maximum temperature difference Δ T max (the temperature difference between maximum temperature on the heating surface and the initial coolant temperature), the average temperature difference Δ T a v e (the temperature difference between average temperature on the heating surface and the initial coolant temperature), and pressure drop Δ p as response variables. An orthogonal experiment was conducted to determine the effects of these four factors on the response variables. The ranges and levels for these four parameters are presented in Table 4.

The orthogonal experiment employed the standard orthogonal array L25 (43 × 31), yielding 25 sets of computational results as presented in Table 5.

Based on the 25 experimental cases presented in Table 4, the structural parameters demonstrate a significant influence on battery performance. The difference between the maximum and minimum values is 21.79% for Δ T max , 136.98% for Δ p , and 46.34% for Δ T a v e (these percentages represent the increase from the minimum to the maximum value). Although the average temperature difference ( Δ T a v e ) has a more pronounced impact on overall system performance, the paramount concern for battery safety is local overheating. Therefore, this study will focus on the implications of the maximum temperature difference ( Δ T max ) in subsequent research. Consequently, subsequent range analysis will focus exclusively on Δ T max and Δ p .

3.2 Range Analysis and Pearson Correlation Analysis

Range analysis was applied to determine the influence of the aforementioned four factors on the experimental response, and to analyze the interaction effects among multiple factors on the response pattern. The range is calculated as follows: xij¯=1n∑1nΔTij(14) xij¯=1n∑1nΔpij(15) Ri=maxxij¯−minxij¯(16) where, x i j ¯ is the mean value of the experimental response for factor i at level j , n is the number of levels, R i is the range of factor i , the range of subscripts i , j are given by: i = α , n , η , d , j = 1 , 2 , 3 , 4 , 5 .

The range analysis results for Δ T max and Δ p are presented in Table 6 and Table 7, respectively. For Δ T max , the factors are ranked in descending order of influence magnitude as follows: R η > R n > R d > R α . For Δ p , the influence magnitude ranking is R η > R n > R α > R d . Collectively, the two key factors influencing both Δ T max and Δ p are the number of unit pairs n and the channel asymmetry ratio η . As n increases, Δ T max decreases while Δ p increases; similarly, as η increases, Δ T max decreases while Δ p increases. The range rankings reflect only the main effects. However, simulation results demonstrate that the cooling performance of the Tesla valve cold plate is governed by the synergistic regulation of multiple parameters, including α , n , η and d . Therefore, a multi-objective optimization approach is required to determine the optimal combination of these parameters.

Range analysis only reflects the overall influence of factor level variations on the performance indicators, but cannot quantify the direction and strength of the linear correlation between factors and indicators. Therefore, the Pearson correlation coefficient r between each structural parameter (actual physical values) and performance indicator was further calculated. The calculation formula is as follows [39]: r=∑i=1NXi−X¯Yi−Y¯∑i=1NXi−X¯2·∑i=1NYi−Y¯2(17) where, X i is the actual value of the structural parameter; Y i is the performance indicator value; and X ¯ are Y ¯ their respective sample means; N = 25 is the number of experiments. The closer r is to 1, the stronger the linear correlation; the sign indicates positive or negative correlation.

Based on the 25 sets of experimental data, the values of the Pearson correlation coefficients between each parameter and the performance indicators were calculated, as shown in Table 8.

As shown in Table 6: For Δ T max , n has the largest correlation coefficient (−0.725), followed by η (−0.574); both are negatively correlated, indicating that increasing the number of unit pairs or the asymmetry ratio helps reduce the maximum temperature difference. For Δ p , η has the largest correlation coefficient (0.617), showing a positive correlation; n ranks second (0.479), also positively correlated; meaning that increasing the the number of unit pairs or the asymmetry ratio can increase the pressure drop. For Δ T a v e , n has the largest correlation coefficient (−0.712), followed by η (−0.562); with a negative correlation, indicating that increasing the number of unit pairs or the asymmetry ratio helps reduce the average temperature difference.

The results of the range analysis and the Pearson correlation coefficient analysis are in good agreement: the number of unit pairs and the channel asymmetry ratio are the most critical parameters affecting the performance of the Tesla valve liquid cooling plate.

3.3 Analysis of Cooling Performance under Multi-Factor Coupling Conditions

According to existing research [40], when fluid enters the Tesla valve through the reverse inlet, it splits into two streams: one enters the straight channel while the other flows into and through the branched section. Due to inertia, the fluid collides with the wall of the curved channel, leading to a considerable pressure drop and energy dissipation. Furthermore, the fluid exiting the curved channel undergoes a secondary collision with the fluid in the straight channel, dissipating additional energy. To investigate this flow behavior and the associated energy dissipation, this section examines the turbulent kinetic energy (TKE), which represents the kinetic energy of turbulent fluctuations, along with the corresponding wall heat flux. These parameters serve as effective indicators of turbulence intensity and heat transfer efficiency.

Two Tesla valves are compared in this work. The first one is the Tesla valve with a shunt angle α = 30°, number of unit pairs n = 5, channel asymmetry ratio η = 0.5, and branch channel width d = 3 mm. In this case, the internal flow exhibits significant turbulence enhancement characteristics. The global TKE distribution is shown in Fig. 5a revealing consistently high turbulence intensity levels throughout the entire structure, on the order of 1 × 10−4 to 6 × 10−2 m2/s2.

Notably, this configuration exhibits a pressure drop of 3440.72 Pa, which is markedly higher than that of other reference structures. Analysis of the local turbulent kinetic energy distribution reveals that this unusually high pressure drop stems from its unique flow path design. The upper and lower Tesla valve bifurcation units are directly connected via branch segments. In this compact two-stage Tesla valve structure, strong flow disturbances are generated as the fluid passes through the junction that connects the two bifurcation units along the main flow direction. Specifically, when the fluid transitions abruptly from the outlet of the first-stage bifurcation unit to the inlet of the second-stage unit, its flow direction changes sharply. Moreover, as the number of unit pairs increases, disturbances generated in the preceding units superimpose at the inlets of the downstream units. This repeated structural arrangement significantly amplifies the turbulent kinetic energy.

Figure 5: Distribution of turbulent kinetic energy and heat flux: (a) TKE; (b) heat flux.

Although the turbulent kinetic energy in this region is considerably higher than that in other configurations, the local fluid maximum temperature difference remains as high as 5.49°C (see Table 4). It can also be observed from Fig. 5b that the high-turbulence-intensity zones formed in local regions of this Tesla valve structure correspond to relatively high local heat flux densities, with a peak value of up to 3.6 × 104 W/m2, which significantly enhances the transient heat transfer capacity. However, the high turbulent kinetic energy is accompanied by substantial energy dissipation and pressure loss, which weakens the effective turbulent perturbation of the thermal boundary layer and leads to highly non-uniform fluid distribution [41]. The majority of the coolant is concentrated in high-perturbation regions, while other areas exhibit inferior heat transfer performance due to insufficient flow kinetic energy. As a result, the overall thermal management efficiency is compromised, and temperature non-uniformity is exacerbated.

The synergistic action of the aforementioned flow characteristics induces a marked surge in turbulent kinetic energy and a sharp intensification of TKE dissipation within the junction region and the immediately adjacent downstream flow field (local TKE peaks reach 6 × 10−2 m2/s2 with a corresponding maximum temperature difference of 5.49°C). This highly disordered flow regime, induced by the uniquely compact coupled configuration, constitutes the fundamental hydrodynamic origin accounting for both the substantially elevated pressure drop relative to alternative Tesla valve cold plate architectures and the limited improvement in heat transfer.

The second one is the Tesla valve cold plate structure with a shunt angle α = 30°, number of unit pairs n = 7, channel asymmetry ratio η = 0.375, and branch channel width d = 4 mm (see Fig. 6a). Exceptional thermal-hydraulic performance was observed in orthogonal testing: a maximum temperature difference as low as 5.08°C (see Table 4) and a pressure drop controlled at 1636.82 Pa. Compared to the high-resistance configuration ( n = 5, d = 3 mm) at the same shunt angle, this structure achieves a balance of low flow resistance and high heat transfer intensity through synergistic optimization of key parameters.

Figure 6: Distribution of turbulent kinetic energy and heat flux: (a) TKE; (b) heat flux.

Compared to the previous Tesla valve cold plate structure, this configuration exhibits a significant reduction in peak turbulent kinetic energy (TKE) throughout the entire domain, accompanied by a progressive downstream attenuation (see Fig. 6a); the TKE within the flow field remains predominantly below 2.24 × 10−2 m2/s2, indicating a dominant low-turbulence flow with ordered structures, which effectively avoids the formation of high-turbulence dissipation zones.

Although this configuration exhibits a lower peak heat flux compared to the previous Tesla valve cold plate structure, its heat flux distribution is more uniform (see Fig. 6b), Indicating that the flow is predominantly governed by moderate-intensity turbulence, thereby avoiding additional energy dissipation caused by localized vortex concentration. This smooth and continuous turbulent structure promotes a more uniform distribution of the coolant within the flow channel, leading to more consistent heat transfer at the walls. Consequently, while maintaining sufficient heat transfer intensity, this configuration effectively mitigates the formation of local hot spots, thereby improving overall temperature uniformity and energy utilization efficiency.

Meanwhile, increasing the branch channel width to 4 mm significantly reduces energy dissipation caused by high-intensity turbulence, effectively lowering the pressure drop. Appropriately decreasing the channel asymmetry factor weakens the strength of separation vortices. Increasing the number of unit pairs promotes spatial dispersion of flow disturbances, preventing local energy concentration and dissipation. Additionally, the expanded heat transfer area contributes to reduced temperature gradients. The synergistic regulation of these multi-parameters enables the configuration to achieve both low pressure drop and low maximum temperature difference.

In conclusion, the magnitude and distribution of turbulent kinetic energy directly affect the uniformity of heat flux distribution and heat transfer efficiency. Excessively high local TKE, while enhancing heat transfer, often comes at the cost of flow field inhomogeneity and increased energy consumption. In contrast, moderate and well-distributed TKE facilitates efficient, uniform, and long-term thermal management performance.

4 Multi-Objective Optimization

To further identify optimal solutions, a multi-objective optimization was performed using Kriging surrogate models based on orthogonal experimental design. With maximum temperature difference Δ T max and pressure drop Δ p as optimization objectives, the structural parameters α , n , η and d were treated as design variables. The optimization workflow is illustrated in Fig. 7.

Figure 7: Optimization procedure flowchart.

4.1 Optimal Latin Hypercube Sampling

During surrogate model optimization, the quality of sampled data directly determines model accuracy and computational efficiency. To balance spatial coverage comprehensiveness with sampling economy, this study employs the Optimal Latin Hypercube Sampling (Optimal LHS) method. Compared to traditional Latin Hypercube Sampling, this approach significantly enhances space-filling uniformity and orthogonality among variables through optimized spatial distribution of sample points [42]. Consequently, it efficiently characterizes design space properties with fewer sampling points.

The flow characteristics of the Tesla valve, such as pressure drop, exhibit high sensitivity to parametric variations. Increasing the sample size facilitates capturing nonlinear responses, particularly the interaction effects among the shunt angle, number of unit pairs, channel asymmetry ratio, and branch channel width. Accordingly, 51 sampling points were selected—corresponding to 12.75 times the number of variables—exceeding the empirical lower threshold (5–10 times) [43]. This ensures the surrogate model’s predictive accuracy for coupling effects among α , n , η and d . Simultaneously, the Optimal Latin Hypercube Sampling method avoids local sparsity in the design space by optimizing the minimum inter-sample distance. This sample size maintains computational feasibility for CFD solutions while adequately supporting subsequent model training and validation. The spatial distribution of these points is depicted in Fig. 8. Response values for all sampling points were computed via Fluent fluid dynamics simulations, providing the foundational dataset for surrogate model construction, with results summarized in Table 9.

Figure 8: Spatial distribution of sampling points.

4.2 Kriging Surrogate Model

The Kriging surrogate model is used to perform response surface fitting on the 51 sets of data obtained above. Kriging is a Gaussian process-based interpolation surrogate model proposed by South African geologist Danie G. Krige, later introduced into the field of engineering optimization by Kuhnt et al. [44]. It is a regression algorithm that predicts random processes based on covariance functions. In this study, the Kriging model was constructed with a Gaussian correlation function and a second-order polynomial regression model, which offers advantages such as high robustness and high fitting accuracy compared to other models. In this model, the predicted response at an unsampled location is calculated through the following equation [45]: yx0¯¯=∑i=1nλiyxi(18) where, y x 0 ¯ ¯ is the predicted response at an unsampled location; y x i denotes the observed response value; λ i represents the undetermined weight coefficient.

Thus, the response value at any point in the design space can be obtained once the weight coefficients λ i are determined. The values of λ i must satisfy the following two constraints:

• (1)Unbiased estimationEyx0¯¯−yx0=0,∑i=1nλi=1(19)where, yx0 is the estimated true value of the assumption.
• (2) Minimization of the variance between estimated and true valuesDyx0¯¯−yx0=min(20)Dyx0¯¯−yx0=−∑i=1n∑j=1nλiλjγxi,xj+2∑i=1nλiγxi,x0(21)where, γxi,xj denotes the semi-variance value with the distance between points xi and xj as the parameter; γxi,x0 denotes the semi-variance value with the distance between points xi and x0 as the parameter.

The accuracy of the surrogate model predictions was evaluated using the R 2 error metric, calculated as follows: R2=1−∑i=1nyxi−yxi¯¯2∑i=1nyxi−yxi¯2(22) where, y x i denotes the actual value, y x i ¯ ¯ denotes the predicted value, and y x i ¯ denotes the mean of the observed data points. An R 2 value approaching 1 indicates higher model accuracy. In engineering applications, an R 2 value greater than 0.9 is generally considered [46] to meet application requirements. As shown in Fig. 9, the Kriging surrogate models constructed in this study achieved R2 values of 0.96138 and 0.97214 for Δ T max and Δ p , respectively.

Figure 9: Accuracy analysis figure for the surrogate model: (a) pressure drop; (b) maximum temperature difference. Note: The blue scatter points represent the sample points corresponding to the actual values from CFD simulation and the predicted values from the model; the red solid line is the linear regression fitting line between the predicted values and the actual values; the black dashed line is the y = x ideal prediction reference line.

4.3 NSGA-II Genetic Algorithm

Based on the constructed Kriging surrogate models, a multi-objective optimization of maximum temperature difference and pressure drop was performed using the Non-dominated Sorting Genetic Algorithm II (NSGA-II) [47], which significantly reduces the computational complexity of NSGA. Its introduced elitist strategy enhances the algorithm’s robustness while improving computational efficiency. In this study, the elitist retention strategy adopts Pareto sorting with an elite ratio of 0.2. The convergence criterion is defined as follows: the algorithm is determined to converge when the average change of the Pareto front in 5 consecutive generations is less than 0.01%. This optimization problem can be mathematically formulated as follows:

The parameter settings in this study were configured as follows: population size = 100, number of generations = 20, crossover probability = 0.8, mutation probability = 0.05, crossover distribution index = 10, and mutation distribution index = 20.

For battery thermal management, the primary objective is to keep the maximum temperature difference within a safe range rather than pursuing its absolute minimization. Once this safety requirement is met, the optimization focus should shift to reducing energy consumption, i.e., pressure drop. In the present study, the Tesla valve cold plate inherently exhibits high flow resistance due to its numerous bends and bifurcations, which leads to a generally higher pressure drop compared with traditional straight-channel structures. An excessively high pressure drop directly increases the pumping power consumption of the system, a critical factor restricting the overall performance of the thermal management system. Therefore, to determine the optimal structural parameters of the Tesla valve cold plate, this study employs the linear weighted sum method to perform multi-objective optimization for the two conflicting objectives of pressure drop ( Δ p ) and maximum temperature difference ( Δ T max ). Based on the comprehensive consideration of system energy consumption and temperature uniformity, this study sets the weight coefficients to prioritize pressure drop, with w 1 = 0.55 (corresponding to Δ p ) and w 2 = 0.45 (corresponding to Δ T max ) [48]. This weight allocation aims to give priority to ensuring low operating energy consumption while controlling the maximum temperature difference of the battery pack within an allowable range.

After obtaining the Pareto optimal solution set (see Fig. 10), all non-dominated solutions are first normalized to eliminate the influence of dimensional differences. The normalization is based on the pressure drop range of approximately 1450–2100 Pa and the maximum temperature difference range of approximately 4.9–5.9°C. Subsequently, the comprehensive evaluation value of each solution is calculated using the following composite objective function [49]: f=w1·Norm(Δp)+w2·Norm(ΔT)(24) where w 1 and w 2 are weight coefficients, subject to w 1 + w 2 = 1 , and w 1 , w 2 > 0 ; and Norm(·) denotes the normalization the objective values. Norm(Obji)=Obji−ObjiminObjimax−Objimin(25) where O b j i is the true value of the objective function, and o b j i max and O b j i min are its maximum and minimum values, respectively, on the entire Pareto frontiers.

The computational and comparative results indicate that the design point with geometric parameters α = 30 ∘ , n = 7 , η = 0 , and d = 4 mm (see Fig. 11) achieves the minimum comprehensive evaluation value. Its performance metrics are Δ p = 1540 Pa (Norm( Δ p ) ≈ 0.138) and Δ T max = 5.42°C (Norm( Δ T max ) ≈ 0.689), yielding a corresponding f ≈ 0.386. Consequently, based on the established weighting criteria, this point is identified as the global optimum.

Figure 10: Pareto frontiers obtained by multi-objective optimization.

Figure 11: Optimized Tesla valve structure.

To validate the reliability of the multi-objective optimization results, simulation verification was performed on the optimized flow channel. As presented in Table 10, the discrepancies between the optimized predicted values and simulation values for both objective functions are minor (maximum relative error < 1.5%). This demonstrates that the Pareto frontiers set possesses high engineering accuracy, meeting the required reliability standards.

The performance results of the optimized structure were compared with those of the initial model ( α = 45°, n = 3, η = 0, d = 2 mm), with the comparative data presented in Table 11. Compared with the initial configuration, the optimized configuration achieves simultaneous thermal uniformity enhancement and pressure drop reduction: the pressure drop Δ p is reduced by 0.36%, the maximum temperature difference Δ T max is decreased by 9.13%, and the average temperature difference Δ T ave is decreased by 15.03%. To further quantify the comprehensive performance improvement, the (PEC) is calculated based on Eq. (10), and the result is 1.18. A PEC value greater than 1 indicates that the optimized structure achieves a synergistic improvement in heat transfer performance with negligible increase in flow resistance, which clearly demonstrates the effectiveness of the optimized design in improving the system’s cooling performance.

5 Conclusion

To ensure the comprehensive performance of lithium-ion batteries, this study focuses on optimizing the thermal efficiency of the Tesla valve cold plate. Orthogonal experimental design was employed to analyze the effects of key structural parameters, including shunt angle, number of unit pairs, channel asymmetry ratio, and branch channel width, on both the pressure drop and temperature uniformity of the cold plate. Building upon these findings, a multi-objective optimization methodology was implemented for the structural design of the cold plate, aiming to simultaneously enhance its thermal management performance and flow characteristics. The principal conclusions are summarized as follows:

• (1)A multi-factor structural model of the Tesla valve cold plate was developed and numerically simulated. The simulation analysis revealed significant interaction effects among the structural parameters-shunt angle, number of unit pairs, channel asymmetry ratio, and branch channel width. Cooling performance varied substantially across different parameter combinations, confirming the necessity of cooperative multi-parameter optimization design.
• (2) Simulation results indicate that variations in the structural parameters induce significant fluctuations in the key metrics: the maximum temperature difference ΔTmax varies by 21.79%, the pressure drop Δp by 136.98%, and the average temperature difference ΔTave by 46.34%. Range analysis for ΔTmax and Δp yielded parameter influence rankings of Rη>Rn>Rd>Rα and Rη>Rn>Rα>Rd, respectively.
• (3) Kriging surrogate models correlating Tesla valve parameters (α,n,η,d) with heat transfer performance metrics (ΔTmax,Δp) were constructed. Multi-objective optimization was subsequently performed using the NSGA-II genetic algorithm, obtaining optimal structural parameters for the Tesla valve cold plate. The optimization results demonstrate that, compared to the initial configuration, the improved flow channel design achieves the dual benefits of enhanced thermal uniformity and reduced pressure drop. Specifically, the pressure drop Δp is lowered by 0.36%, the maximum temperature difference ΔTmax is reduced by 9.13%, and the average temperature difference ΔTave is decreased by 15.03%.

However, the research conclusions in this paper are entirely based on numerical simulations (CFD), lacking experimental validation, and fail to consider the impact of microchannel fabrication accuracy on actual performance. In the model simplification, it is assumed that the battery generates heat uniformly, while the variation of coolant properties with temperature is neglected, thus failing to reflect electrochemical inhomogeneity and viscosity changes at high temperatures. Additionally, the discussion scope of geometric parameters is limited. Meanwhile, combining Tesla valves with nanofluids to explore the regulatory effects of nanoparticles on turbulent transport characteristics, boundary layer perturbation, and local heat transfer enhancement is expected to break through the existing limits of cooling efficiency. Breakthroughs in the aforementioned directions will promote the development of battery thermal management systems and provide core guarantees for the safe application of high-energy-density batteries.