Optimization and Key Factor Analysis of Immersion Cooling Performance for 18650 Lithium-Ion Batteries in a Serpentine Channel Based on Response Surface Methodology

1. Introduction

The battery is one of the core components of a pure electric vehicle. At present, the power batteries of electric vehicles mainly include supercapacitors, nickel-based batteries, lead-acid batteries, lithium-ion batteries, zinc-air batteries, etc. Among them, lithium-ion batteries have become the preferred solution for electric vehicle manufacturers due to their high energy density, light weight, and long service life. However, the performance, lifespan and safety of lithium-ion batteries are very sensitive to temperature. An efficient battery thermal management system should ensure that the battery operates within the optimal working temperature range (15-35 °C) [1,2]and make the temperature distribution on the battery surface uniform. If the temperature of lithium-ion batteries exceeds 50°C, the charging efficiency and service life will be greatly reduced [3]. In the temperature range of 70-100°C, lithium-ion batteries may triggering thermal runaway and other safety accidents [4].

At present, the most commonly used thermal management methods in lithium-ion batteries include air cooling, liquid cooling, phase change materials (PCMs) and heat pipe cooling technologies [5]. The air cooling technique pumps ambient air to purge the battery module to accelerate the heat dissipation. It has the characteristics of simple structure, easy implementation and relatively low cost, but the cooling efficiency is low. Akinlabi and Solyali [6] reviewed the air-cooling battery thermal management systems, and concluded that it is difficult to control the temperature of high-energy-density lithium batteries by air-cooling battery thermal management system (BTMS). The PCMs cooling system is a passive cooling method, and the PCMs used usually has a large latent heat of fusion and an ideal melting point, and can store a large amount of heat generated during the charging and discharging process of the battery. However, the PCMs usually have low thermal conductivity and lead to poor thermal propagation. The PCMs should be modified[7,8,9] or added with other additives [10,11], and used in conjunction with other cooling methods [12,13,14]. Similar to the phase change material cooling system, heat pipe cooling usually needs to be used in combination with other cooling methods [15,16,17], and the heat transfer performance of the heat pipe will be reduced when the vehicle is affected by shock and vibration.

The liquid cooling system has gradually become the hotspot of BTMS. According to the contact type between the coolant and the batteries, it is usually divided into indirect liquid cooling (ICLC) and direct liquid cooling (DCLC, also called as immersion cooling) [18]. Indirect contact liquid cooling usually uses water or a mixture of water and ethylene glycol as the coolant. The coolant flows in channels of pipes, cold plates or water jackets, and removes the heat generated by the battery. Due to the electrical conductivity of the coolant, the indirect cooling system needs to be well sealed to prevent liquid leakage from causing a short circuit. At the same time, the thermal contact resistance between the metal separator and the battery needs to be considered. Although the indirect cooling system has good heat dissipation capacity, the structure of the system is relatively complex, the whole system is very bulky, expensive, and the maintenance cost increases. In the early days of direct contact liquid cooling, deionized water, silicone oil, mineral oil, and refrigerants were commonly used as coolants[19]. Although liquid has a high thermal conductivity compared with air, which can effectively reduce the maximum temperature of the batteries, the direct contact liquid cooling system still greatly relies on the liquid cooling arrangement to enhance the convective heat transfer coefficient, thereby improving the overall temperature uniformity of the battery pack.

In recent years, the immersion cooling system for power battery thermal management has become a new research hotspot. The battery is immersed in a dielectric coolant. Unlike the indirect liquid cooling system, the immersion cooling system has no contact thermal resistance, smaller volume, higher cooling efficiency, and more uniform temperature distribution. To prevent short circuits and electrochemical corrosion between the battery and the working fluid, the coolant should be selected to have insulation properties, non-toxicity, chemical inertness, high thermal conductivity, and flame-retardant characteristics [20,21]. Such fluids typically offer superior thermal performance for BTMS. Nelson and colleagues [22] developed a battery module comprising 12 serially connected cells, incorporating 1 mm wide flow channels filled with silicone oil between adjacent cells. Their experiments demonstrated that silicone oil, serving as a thermal medium, offers significantly better heat dissipation performance than air cooling. Karimi et al. [23] carried out a comparative assessment of water, silicone oil, and air as cooling media, showing that liquid-based cooling can effectively keep battery temperatures within an acceptable operational range, although minor internal temperature variations may still occur. Hirano et al. [24] explored the use of 3M’s Novec7000 for immersion cooling applications. The results indicated that when the battery undergoes 20C high-rate charge and discharge cycles, submerging it in Novec7000 successfully maintains maximum battery temperature T_max at a stable 35°C. Tan et al. [25] utilized HFE-6120 as the working fluid, and conducted CFD simulations to analyze its key operational parameters in DCLC system. Satyanarayana et al. [26] compared the immersion cooling effectiveness of mineral oil and terminal oil under varying discharge rates. Their findings revealed that low-viscosity mineral oil is more efficient in managing T_max, achieving a reduction of up to 49.16% in T_max compared to natural air convection cooling. Patil et al. [27] investigated the direct cooling performance of Li-ion batteries using dielectric fluid immersion cooling combined with tab cooling. By inserting baffles, the disturbance and residence time of the coolant are increased. The results demonstrated significantly improved thermal management, with up to a 46.8% reduction in maximum tab temperature and better overall pack performance compared to conventional cooling methods. Wu et al. [28] compared single-phase and two-phase liquid cooling systems for large lithium-ion batteries, showing that the single-phase system with deionized water and copper foam provides superior thermal performance, achieving lower temperatures and better uniformity during high-rate discharge. Choi et al. [29] proposed a hybrid immersion-cooling structure with one pass partition and graphite fins. The performance of the structure was investigated with variations in key design parameters such as flow direction, battery spacing, pass partitions, and thermal conductive materials. The results show that, compared to conventional structures, this hybrid immersion structure achieves a maximum temperature reduction of 6.7 °C and a temperature difference reduction of 3.0 K under nearly the same weight condition. Zhao et al. [30] developed a hybrid battery thermal management system combining direct liquid cooling with forced air cooling. A jacket was designed outside the battery, with liquid coolant filled between the battery case and the jacket to achieve direct cooling, while air cooling was applied externally. The effects of gap spacing between the battery and liquid-cooling jacket, number of cooling pipelines, liquid flow rate, and fan position on cooling performance were analyzed through numerical simulation. Febriyanto et al. [31] investigated the thermal performance of a Serpentine Channel Immersion Cooling (SCIC) system. The results showed that increasing the coolant flow rate enhances the heat transfer coefficient and reduces battery surface temperature, with a notable balance between heat dissipation and pumping power at 1 LPM (liters per minute).

The above-mentioned studies demonstrate that improvements in the single-phase immersion cooling performance of batteries are mainly achieved through the use of higher-performance coolants or the optimization of flow channel structures. While structural optimization of flow channels has proven effective in enhancing heat dissipation, there remains a significant lack of quantitative research on the relationships between channel design parameters, operational parameters, and thermal response indicators—particularly the maximum battery temperature T_max and pump power consumption P_w. Current research is largely confined to experimental validations or simulation trials based on discrete parameter values, with no universally applicable mathematical correlation model yet established. In this study, 18650-type lithium-ion batteries are selected as the research subject, and HFE-7100 was chosen as the immersion coolant. A battery thermal management model incorporating a serpentine channel configuration was developed using the ANSYS Fluent computational fluid dynamics platform. The focus is on investigating the impact mechanisms of key factors such as partition length, coolant flow velocity, and ambient temperature on thermal performance. To overcome the limitations of traditional single-variable analysis methods, this study employed Response Surface Methodology (RSM) to construct a quadratic polynomial regression model. Experimental conditions were designed using the Box-Behnken Design (BBD) approach, enabling systematic quantification of the contributions of each design variable and their interactions on T_max and P_w. This mathematical model accurately captures the nonlinear variation patterns of thermal responses under multi-parameter coupling conditions, providing a solid theoretical foundation for the optimized design of immersion cooling systems.

2. Methodology

2.1. Physical Model

The battery module’s geometric design is derived from literature [31], as shown in Figure 1. It uses a 290 mm × 120 mm × 75 mm rectangular housing. Twenty-four 18650 LiFePO₄ batteries are connected in series. The individual battery cell model has dimensions of 18 mm diameter × 65 mm height, represented by a cylindrical approximation. The spacing between adjacent columns is 24 mm, while the spacing between rows is 31 mm. The layout follows a six-loop serpentine-channel pattern with 4 batteries per row. All batteries maintain identical electrode orientation. Partitions physically separate battery groups. The module provides 78.6 V nominal voltage. Table 1 summarizes parameters of 18650-type batteries, while Table 2 lists key thermophysical properties of HFE-7100 coolant (data sourced from [31]).

Figure 1. Schematic diagram of the battery module: (a) Three-dimensional view of battery holder design; (b) Flow channel configuration.

Figure 1. Schematic diagram of the battery module: (a) Three-dimensional view of battery holder design; (b) Flow channel configuration.

Table 1. Technical parameters of lithium-ion battery [31].

Parameter	Specification	Parameter	Specification
Type	LiFePO4	Max Charge Current	1C
Nominal Capacity	1800 mAh	Weight	46±1 g
Minimum Capacity	1775 mAh	Diameter	18 mm
Nominal Voltage	3.2 V	Height	65 mm
Internal Resistance	30 mΩ	Battery arrangement	24S
Discharge Cut-Off Voltage	2.5 V	Total Nominal Voltage	76.8 V
Max Charge Voltage	3.65 V

Table 1. Technical parameters of lithium-ion battery [31].

Parameter	Specification	Parameter	Specification
Type	LiFePO4	Max Charge Current	1C
Nominal Capacity	1800 mAh	Weight	46±1 g
Minimum Capacity	1775 mAh	Diameter	18 mm
Nominal Voltage	3.2 V	Height	65 mm
Internal Resistance	30 mΩ	Battery arrangement	24S
Discharge Cut-Off Voltage	2.5 V	Total Nominal Voltage	76.8 V
Max Charge Voltage	3.65 V

Table 2. Thermophysical properties of HFE-7100 [31].

Properties	Values	Unit
Boiling point	334.15	K
Specific heat	1170	J/(kg⋅K)
Latent heat vaporization	112	kJ/kg
Thermal conductivity	0.068	W/(m⋅K)
Liquid density	1418.64	kg/m3
Vapor density	0.98	kg/m3
Kinematic viscosity	3.008×10-7	m2/s

Table 2. Thermophysical properties of HFE-7100 [31].

Properties	Values	Unit
Boiling point	334.15	K
Specific heat	1170	J/(kg⋅K)
Latent heat vaporization	112	kJ/kg
Thermal conductivity	0.068	W/(m⋅K)
Liquid density	1418.64	kg/m3
Vapor density	0.98	kg/m3
Kinematic viscosity	3.008×10-7	m2/s

2.2. Governing Equations

Some assumptions are made here to simplify the problem:

(1)

Constant physical properties are maintained for both liquid coolant and battery components.

(2)

Perfect interfacial contact is established between solid and fluid domains (no-slip boundary condition).

(3)

Incompressible flow is adopted for coolant.

(4)

Buoyancy effects on fluid flow are excluded from the simulation.

(5)

Thermal radiation is omitted in the heat transfer simulation.

For the HFE-7100 liquid, the continuity equation can be written as:

$${\displaystyle \frac{\mathit{\partial }{\rho }_l}{\mathit{\partial }t}}+\nabla ·({\rho }_{l}u)=0$$

(1)

where ρ_l is the density of the fluid, u is velocity component and t is time.

Momentum equation of single-phase HFE-7000 flow:

$${\rho }_l({\displaystyle \frac{\mathit{\partial }u}{\mathit{\partial }t}}\text{}+u\cdot \nabla u)=-\nabla p+\nabla ({\mu }_{eff}\text{}\nabla u)$$

(2)

where p is the pressure, µ_eff is the effective viscosity and has the form of:

$${\mu }_{eff}=\mu +{\mu }_{\mathrm{t}}$$

(3)

where µ is molecular viscosity, µ_t is turbulent viscosity. The turbulent viscosity µ_t is directly related to the turbulent kinetic energy k and the dissipation rate ε:

$${\mu }_{\mathrm{t}}=C_{\mu }{\rho }_l{\displaystyle \frac{k^2}{\epsilon }}$$

(4)

The transport equations in k-ε model are shown as:

$${\displaystyle \frac{\mathit{\partial }({\rho }_{l}k)}{\mathit{\partial }t}}+\nabla ·({\rho }_{l}uk)=\nabla ·\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_k}}\right)\nabla k\right]+G_k-{\rho }_l\epsilon$$

(5)

$${\displaystyle \frac{\mathit{\partial }({\rho }_l\epsilon )}{\mathit{\partial }t}}+\nabla ·({\rho }_{l}u\epsilon )=\nabla ·\left[\left(\mu +{\displaystyle \frac{{\mu }_{\mathrm{t}}}{{\sigma }_{\epsilon }}}\right)\nabla \epsilon \right]+C_1{\displaystyle \frac{\epsilon }{k}}{G}_k-C_2{\rho }_l{\displaystyle \frac{{\epsilon }^2}{k}}$$

(6)

where G_k is the production of turbulence kinetic energy. Values of the constants in the k-e equations are given as follows: C_μ=0.09, C₁=1.44, C₂=1.92, σ_k=1.0, σ_ε=1.3.

The energy conservation equation of the coolant is:

$${\rho }_l{c}_{p,l}\left({\displaystyle \frac{\mathit{\partial }T}{\mathit{\partial }t}}+u\cdot \nabla T\right)=\nabla \cdot ({\lambda }_{\mathrm{eff}}\nabla T)$$

(7)

The energy conservation equation of the battery cell is expressed as:

$${\rho }_s{C}_{p,s}{\displaystyle \frac{\mathit{\partial }T}{\mathit{\partial }t}}=\nabla ·({\lambda }_s\nabla T)+S_b$$

(8)

where ρ_s, C_p,s, and T are the battery material density, specific heat capacity, and temperature, respectively. λ_s is the thermal conductivity of the battery. The heat source term, S_b, is the volumetric heat generation rate in the battery cell which can be divided into reversible heat and irreversible heat according to their origins, and can be presented as:

$$S_b=I(V_t-U)+I\left(T{\displaystyle \frac{\mathit{\partial }U}{\mathit{\partial }T}}\right)$$

(9)

where I represents the current, V_t represents the voltage during operation, U represents the voltage when the circuit is open, T represents the temperature of battery, and ∂U/∂T represents the coefficient of entropy.

Discharge current I at any discharge rate can be determined by multiplying the C_rate by the battery’s nominal capacity:

$$I=C_{\mathrm{rate}}\cdot C_{\mathrm{cap}}$$

(10)

The ideal pump power is calculated as:

$$P_{\mathrm{w}}=\mathrm{\Delta }P\cdot V$$

(11)

where P_w is pump power, ΔP is the pressure drop between the coolant inlet and outlet, V is the volumetric flow rate of the coolant.

2.3. Solution Strategy

2.3.1. Boundary Conditions

In this study, the coolant velocity is specified as an input parameter, with the inlet of the computational domain defined as a velocity inlet boundary. So, the velocity inlet boundary condition is applied, with the flow direction orthogonal to the inlet plane:

$${\mathrm{u}}_{\mathrm{y}}={\mathrm{u}}_{\mathrm{z}}=0，{\mathrm{u}}_x={\mathrm{u}}_{in}，T=T_{in}$$

(12)

On the outer surfaces of the battery module, a no-slip wall boundary condition is applied, and a convective heat transfer coefficient h of 20 W/(m²·K) is assigned to model heat dissipation. T_amb is the ambient temperature, and satisfies the following formula:

$$-{\lambda }_s{\left.{\displaystyle \frac{\mathit{\partial }\mathrm{T}}{\mathit{\partial }n}}\right|}_{wall}=h\left(T_{wall}-T_{amb}\right)$$

(13)

At the outlet, the boundary is defined as an outflow condition to accommodate variable flow rates:

$${\displaystyle \frac{\mathit{\partial }{\mathrm{u}}_{\mathrm{x}}}{\mathit{\partial }\mathrm{x}}}={\displaystyle \frac{\mathit{\partial }{\mathrm{u}}_{\mathrm{y}}}{\mathit{\partial }\mathrm{x}}}={\displaystyle \frac{\mathit{\partial }{\mathrm{u}}_{\mathrm{z}}}{\mathit{\partial }\mathrm{x}}}=0，{\displaystyle \frac{\mathit{\partial }T}{\mathit{\partial }\mathrm{x}}}=0$$

(14)

At the fluid-solid interface, continuity of temperature and heat flux is enforced to ensure energy conservation:

$$T_l=T_s,{\lambda }_l{\displaystyle \frac{\mathit{\partial }T}{\mathit{\partial }n}}{|}_l={\lambda }_s{\displaystyle \frac{\mathit{\partial }T}{\mathit{\partial }n}}{|}_s$$

(15)

2.3.2. Mesh Independence Check

To ensure the reliability of numerical simulations for battery immersion cooling, a comprehensive grid independence study was conducted by systematically refining the mesh resolution from coarse to medium and fine, as depicted in Table 3. The results indicate that the maximum battery temperature under the medium grid density deviates by less than 0.86% compared to the finest mesh, while temperature gradients exhibit a maximum deviation of 1.57%. The detailed mesh structure of the computational domain is illustrated in Figure 2. Given these results, the medium grid configuration was ultimately selected as a balance between computational efficiency and solution accuracy, ensuring y⁺ values below 5 in boundary layers. The simulations were performed using a pressure-based solver with the pressure correction algorithm for pressure-velocity coupling. The Navier-Stokes equations were solved via a SIMPLEC algorithm, while the convective and diffusive terms in the energy equation were discretized using a second-order upwind scheme. Convergence criteria were set to 10⁻⁴ for continuity and momentum equations and 10⁻⁶ for energy equations throughout all simulations.

2.4. Response Surface Design

The Response Surface Method (RSM) is a statistical experimental design and data analysis technique used to systematically study the quantitative relationship between multiple influencing factors and a response variable. By constructing a mathematical model (such as a quadratic polynomial equation), RSM efficiently identifies the effects of key parameters on the output while minimizing the number of required experiments. Its main advantages include significantly reducing experimental costs, avoiding the high resource consumption associated with full factorial designs, and enabling precise determination of optimal operating conditions. Additionally, it allows for accurate prediction of response values under various conditions, thereby enhancing both the efficiency and reliability of the experimental outcomes. In this study, RSM was employed to establish the nonlinear response relationships between operational parameters (e.g., flow rate, coolant temperature), battery arrangement parameters (e.g., partition length), and battery temperature distribution, providing a theoretical foundation and practical guidance for optimizing thermal management strategies.

In this study, the partition length (L_p), battery discharge rate (C_rate), coolant volumetric flow rate (V), coolant inlet temperature (T_in), and ambient temperature (T_amb) were selected as design variables, while the battery maximum temperature (T_max), average temperature (T_ave), and pump power (P_w)were chosen as response values. A response surface methodology with five factors at three levels and two response values was employed. The coding rule was defined as follows: 0 for the center level, -1 for the minimum level, and 1 for the maximum level. The coded matrix based on the Box-Behnken Design (BBD) is presented in Table 4.

Using Design expert 13.0 software and following the design principles of the BBD, a total of 46 experimental runs, including six center points, were designed for the five factors at three levels. The design scheme and results of the BBD experiments are summarized in Table 5.

3. Results and Discussion

3.1. Model Validation

Prior to analyzing the impact of battery design parameters on thermal behavior, validating the accuracy of the numerical model is essential. This study uses experimental data from Febriyanto et al. [31] to verify the numerical model. Figure 3 presents a comparison between the simulated data and the experimental results from Ref. [31]. Figure 3a shows the comparison of temperature values at the observation points across the six channels under 1C and 2C discharge rates. As can be observed from the figure, the simulated temperatures at all observation points are in good agreement with the experimental measurements. Under the 1C discharge rate, the maximum absolute error between the simulated and experimental values is 0.7 °C, with a corresponding relative error of 2.82%. Under the 2C discharge rate, the maximum absolute error is 1.0 °C, and the maximum relative error is 3.71%. Figure 3b compares the average temperatures across all six observation points between the simulation and experiment. It can be seen that, under the same discharge rate, the average temperature decreases with increasing coolant flow rate, and the simulated results align well with the experimental data. The maximum error occurs under the 1C discharge rate with a coolant flow rate of 0.5 LPM, where the absolute error is 0.58 °C, corresponding to a relative error of 2.17%. The results demonstrate a high degree of consistency between the numerical simulation outcomes and the experimental data, not only validating the model’s capability to accurately predict the actual temperature distribution but also confirming its reliability as a foundation for further investigations.

3.2. Analysis of Variance

3.2.1. ANOVA of Battery Temperature

Table 6 is the ANOVA results for T_max after backward elimination. The ANOVA results indicate that the constructed quadratic polynomial model possesses extremely high statistical significance. The F-value is a ratio obtained by dividing the model mean squares by the error of mean squares. It is a statistic used to test the overall significance of the model. A higher F-value indicates that the model explains a substantial portion of the variability in the data, rather than variations caused by random noise. This helps determine whether the independent variables significantly affect the dependent variable. The p-value is the probability used to decide whether to reject the null hypothesis. The null hypothesis typically states that all model coefficients are equal to zero (i.e., the model has no explanatory power). A smaller p-value suggests that the observed data (or more extreme results) would be less likely to occur under the null hypothesis. Generally, if the p-value is less than 0.05, the result is considered statistically significant. In Table 6, the model’s F-value reached 4768.12, significantly greater than 1, with a corresponding p-value less than 0.0001. This demonstrates that the model’s ability to predict T_max is statistically very significant. According to the statistical standards of response surface methodology (RSM), a model p-value < 0.01 indicates good fitting accuracy, allowing this response surface approximation model to be employed in subsequent optimization design. In this study, the p-values for all main effects and interaction effects were significant (all < 0.01 or < 0.05), further confirming the model’s validity.

The ANOVA results reveal the significance and relative importance of each factor’s influence on T_max. Based on the magnitude of the F-values, the factors are ranked in descending order of influence on T_max as follows: T_in > C_rate > V > L_p > T_amb. T_in had the highest F-value (48033.91) and the largest Sum of Squares (336.81), indicating it is the most significant factor affecting T_max, with a contribution far exceeding that of other variables. Regarding interaction effects, the interactions V·T_in and V·T_amb were the most significant (F=26.37), followed by C_rate·V (F=17.47). The existence of these interaction terms indicates that the effects of the factors on T_max are not independent; rather, they exhibit complex synergistic or antagonistic effects. For instance, the significant V·T_in and V·T_amb interactions suggest that optimizing the level of V requires simultaneous consideration of the levels of T_in and T_amb to achieve the optimal T_max value.

The R² value measures the proportion of total variation explained by the model. Ranging from 0 to 1, a value closer to 1 indicates a better fit of the model to the data. However, R² has a limitation: it tends to increase (or at least not decrease) as more independent variables are added to the model, even if these variables are unrelated to the response variable. This may lead to overfitting. To address this limitation of R², the adjusted R² is introduced. It accounts for the number of independent variables in the model by imposing a penalty for each additional variable. Consequently, adjusted R² increases only if a newly added variable genuinely enhances the model’s explanatory power. This makes it a more reliable metric for comparing models with different numbers of independent variables. Predicted R² evaluates the model’s predictive capability. Based on cross-validation principles, it partitions the dataset, trains the model on one subset, and tests its predictive accuracy on the remaining subsets. A higher predicted R² indicates better performance on unseen data. If predicted R2 is significantly lower than R², it suggests potential overfitting. The S-value represents the average deviation between observed and predicted values. As the standard deviation of model residuals, it estimates the uncertainty of predictions. Generally, a smaller S-value implies closer alignment of predictions with actual observations, reflecting higher model precision. For a model to demonstrate high fitting accuracy and robust predictive capability, the following conditions should be met: R² and adjusted R² are both high and relatively close, indicating strong explanatory power without over-reliance on redundant variables. Predicted R² is high and approximates R², confirming strong generalization ability and predictive performance on new data.

The predictive capability of the response surface model was evaluated using the Predicted R² metric. The result shows a Predicted R² of 99.94%, differing only 0.02% from the Adjusted R² of 99.91%. This aligns with the model adequacy criterion stipulating a difference of less than 0.2% between Predicted R² and Adjusted R². This demonstrates that the model not only fits the experimental data exceptionally well but also possesses robust predictive power for new data, establishing it as a reliable basis for subsequent optimization processes.

Table 7 is the ANOVA results for T_ave after backward elimination. The ANOVA table indicates that the model is highly significant as a whole, with an F-value of 4652.34 and a corresponding p-value much smaller than 0.0001, demonstrating that the model’s predictions for T_ave are statistically significant. The model sum of squares is as high as 353.79, while the residual sum of squares is only 0.235. This large discrepancy further confirms the excellent fitting performance of the model. From the ANOVA results, it can be observed that all examined factors and their interactions significantly affect T_ave (p < 0.05). Among them, T_in exhibits the highest F-value (48,760.41), indicating that it has the most significant influence on T_ave. This is followed by L_p (F = 477.83) and C_rate (F = 831.46), while V (F = 62.43) has a relatively smaller but still significant effect. Regarding interaction effects, the combination of V and T_amb (F = 26.75) shows the lowest F-value among interactions; however, its p-value remains below 0.0001, indicating that although its impact is less pronounced than main effects, it still plays a non-negligible role. In terms of quadratic terms, L_p·L_p (F = 16.98), C_rate·C_rate (F = 14.01), and V·V (F = 62.43) are all statistically significant, suggesting clear nonlinear relationships between these factors and T_ave. The model also demonstrates excellent goodness-of-fit performance, with very high R², adjusted R², and predicted R² values. The standard deviation S of the model is only 0.0831, indicating minimal average deviation between predicted and actual values, which reflects the model’s high accuracy. These results confirm that the response surface model accurately characterizes the relationship between T_ave and the design parameters.

3.2.2. ANOVA of Pump Power

Table 8 is the ANOVA results for P_w after backward elimination. The ANOVA results indicate that the response surface model possesses extremely high significance and predictive capability for the P_w response value. The overall model F-value is remarkably high at 62311.83, with a corresponding p-value less than 0.0001, demonstrating the model’s high statistical significance. The model’s multiple coefficient of determination (R²) is 99.99%, the Adjusted R² is also 99.99%, and Predicted R² is 99.97%. All three values are very close to 100%, with negligible differences (a mere 0.02%), indicating an excellent fit to the data, explaining nearly all the variation in the response value. The standard deviation (S-value) is 0.0304, reflecting minimal differences between the model’s predicted values and actual values, confirming the model’s high precision.

The ANOVA results further indicate that all main effects and interaction/squared terms exert an extremely significant influence on the P_w response value (P<0.01), though the magnitude of their effects varies considerably. By comparing the magnitude of the F-values, the hierarchy of influences of the factors on the P_w response value can be determined: The V variable and its squared term have the most significant impact, with a main effect F-value of 275000 and a squared term F-value of 35538.51, substantially higher than those of other factors. This indicates that the V variable not only has a highly significant effect but also exhibits a strong non-linear relationship with the P_w response value. The L_p variable and its squared term also exhibit a highly significant influence, but their F-values are relatively smaller (main effect F-value = 67.45, squared term F-value = 37.33), indicating a lesser impact compared to the V variable. The F-value for the interaction term L_p·V is 10.79. Although significant, it is considerably lower than the F-values of the main effects and squared terms, suggesting that the interaction effect between L_p and V on the P_w response value is relatively weak. This interaction’s significance might be attributable to the characteristics of the experimental design rather than indicating a substantial synergistic or antagonistic effect within the actual process. Consequently, during optimization, primary consideration should be given to the main effects and squared terms, with the interaction effect treated as a secondary factor for adjustment.

3.3. Regression Model of Responses

The regression equation for evaluating the objective functions G including T_max, T_ave and P_w is obtained as follows:

$$\begin{array}{l}G=a_0+a_1\cdot L_{baffle}+a_2\cdot C_{rate}+a_3\cdot V+a_4\cdot T_{in}+a_5\cdot T_{amb}+a_6\cdot L_{baffle}\cdot V+a_7\cdot C_{rate}\cdot V\\ +a_8\cdot T_{in}\cdot V+a_9\cdot T_{amb}\cdot V+a_{10}\cdot L_{baffle}^2+a_{11}\cdot C_{rate}^2+a_{12}\cdot V^2\end{array}$$

(16)

The coefficients of the regression response surface model are listed in Table 9. Figure 4 shows the simulated values and regression fitted values of T_max and P_w in the quadratic response surface regression model. The maximum relative errors between the simulated and fitted values of T_max and P_w are within acceptable ranges, indicating good agreement and suggesting that the established response surface regression model is highly reliable.

3.4. Response Surface Analysis

Figure 5 illustrates response surface contours, demonstrating the interaction effects of different parameters on T_max and P_w. Figure 5a, 5b, and 5c indicate that T_max decreases as the coolant flow rate V increases, but increases with higher C_rate, T_in, and T_amb. Overall, the most significant factors affecting T_max are T_in and C_rate, whereas the effects of V and T_amb are relatively minor. Theoretically, to control T_max, priority should be given to reducing T_in and C_rate, i.e., minimizing the coolant inlet temperature and battery charge/discharge rate. However, lowering T_in means that after the coolant absorbs heat, it needs to be cooled down again through an external cooling cycle, which increases energy consumption. Moreover, the battery charge/discharge rate is influenced by many factors during vehicle operation and cannot be simply determined by the thermal management system alone. Therefore, adjusting the coolant flow rate becomes an effective method for controlling T_max. It should be noted that when V exceeds 1.1 LPM, further increasing the coolant flow rate has little effect on reducing T_max.

Figure 5d illustrates the interaction effect of V and L_p on P_w. As shown in the figure, pump power consumption significantly increases with the coolant flow rate V. Meanwhile, as L_p changes, P_w initially decreases and then increases, though its impact is much less pronounced compared to V. Specifically, when the coolant flow rate is below 0.7 LPM, the pump’s power consumption is lower.

In summary, regulating the coolant flow rate is one effective approach to controlling the T_max of the electric vehicle battery system. However, attention must be paid to the diminishing efficiency gains at high flow rates. Additionally, considering the substantial increase in pump power consumption with increased flow rate, the design of the thermal management system should balance cooling performance and energy consumption.

According to the response surface analysis, considering factors T_max, T_ave and P_w as the objective function, the optimal point of the design factors is shown in Table 10. At this optimal point, the cooling system achieves minimal pump power while ensuring effective heat dissipation, thus balancing coolant flow rate and pump power. The composite desirability value is 1, indicating excellent overall performance, with all design parameters working synergistically as expected.

4. Conclusions

This paper investigates the immersion cooling performance of 18650 lithium-ion power batteries in a serpentine channel. Using a Box-Behnken Design, several key design factors were selected, including coolant flow rate, coolant inlet temperature, battery discharge rate, partition length, and ambient temperature. The maximum battery temperature and pump power consumption were chosen as the response variables. Computational Fluid Dynamics methods were employed to simulate and calculate these responses. The Analysis of Variance (ANOVA) method was applied to systematically evaluate the statistical significance and relative importance of each design variable on the two response variables. Key conclusions are summarized as follows:

(1) For the T_max response, the ANOVA revealed a clear hierarchy of factor influence: the main effects were ranked as T_in > C_rate > V > L_p > T_amb. Among these, the T_in factor demonstrated overwhelming dominance, with its exceptionally high F-value reaching 48033.91 and a sum of squares of 336.81. These values significantly exceeded those of other variables, indicating that T_in is the most critical factor governing T_max. Additionally, two interaction effects, V·T_in and V·T_amb, also showed significant effects (F=26.37), suggesting their synergistic influences cannot be ignored.

(2) For the P_w response, the ANOVA results provided equally important insights: both the main effect of V and its squared term (V²) exhibited extreme significance (F-values of 275000 and 35538.51, respectively). This indicates that the effect of V on P_w is not only highly significant but also exhibits a strong non-linear relationship. In contrast, the main effect of L_p and its squared term (L_p²) had relatively minor effects (F-values of 67.45 and 37.33, respectively). Although the L_p·V interaction effect was statistically significant (F=10.79), its magnitude was far lower than the main effects and squared terms. This finding provides crucial guidance for modeling the P_w response: priority should be given to the main effect and non-linear effect of V, while interaction effects can be considered secondary adjustment factors.

(3) After optimization targeting minimization of T_max, T_ave and P_w, the optimal design values for L_p, C_rate, V, T_in, and T_amb were determined to be 89.5 mm, 1.08 C, 0.51 LPM, 20 °C, and 25.62°C respectively. The corresponding optimized values are: T_max = 22.87°C, T_ave = 21.67°C, and P_w = 0.279 mW.

This study clarifies the significant effects and hierarchical structure of various factors on T_max and P_w, thereby providing a theoretical foundation and practical guidance for parameter optimization in complex systems.