Project 2 - 3D CFD modelling of Air cooling system and liquid cooling system for battery thermal management

1. Comparative study of thermal performance of air-cooled and liquid-cooled battery modules:-

Air-cooled module:-

The temperature distribution over the module surface with the air-cooling system at the end of the discharge process. The flow rate and temperature of the air at the inlet of the cooling system are 3 L/s and 25 °C, respectively. As the hottest area on the module surface is located near the rear side of the module. Moreover, the hottest area is also concentrated on the left side of the module, where the Outlet manifold is located.

The average temperature distribution of the cells within the air-cooled module at the end of the discharge process for the different flow rates and the inlet temperature of 25 °C is depicted. A significant variation in the average temperature is observed between the cells. As it is seen, for all flow rates, the highest and lowest average temperatures in the air-type BTMS correspond to cell 11 and cell 1, respectively. Moreover, a significant temperature variation is observed between the cells.

The thermal behaviour of the air-cooled module is strongly dependent on the airflow characteristics in the cooling channels between the cells. In this regard, the flow rate, and the average temperature of the air at the inlet of the cooling channels between the cells. Moreover, the temperature distribution of the air within the cooling system is depicted. These results are related to the end of the discharge process, and the flow rate and temperature of the air at the inlet of the cooling system are 3 L/s</math>">L/s and 25 °C. The horizontal axis displays the index of cooling channels. There are 13 cooling channels between the cells, and the cooling channels are numbered the same as battery cells, from the front side of the module (where the inlet and outlet are located) toward the rear side. It is seen that the flow rate inside the cooling channels continually decreases as the channel index number increases from 1 to 11, and after that, it slightly increases. The average temperature of the air at the inlet of the cooling channels rises from channel 1 to channel 7. After that, a small variation in air temperature is observed from channel 7 to channel 11, and then it reduces. The air temperature increases as it passes through the cooling channels. Hence, the air temperature in the outlet manifold is much higher than the inlet manifold. These hydrodynamic and thermal characteristics of the air flow confirm the local and average temperature distribution of the cells within the module. Since the heat transfer from the cells to the air is proportional to the temperature and the flow rate of the air inside the cooling channels, the reduction of air flow and increment of air temperature in the cooling channels from the front side toward the rear side, causes the hottest area of the module to be revealed near the rear side. As the flow rate in the cooling channels 11 and 12 are the minimum ones, that is why cell 11 has the highest average temperature in the module. Moreover, the higher temperature of the air in the outlet manifold results in lower heat dissipation from the left side of the module, and consequently, the hottest area of the module is concentrated on the left side.

Liquid Cooled Module:-

The temperature distribution over the module surface with the liquid cooling system at the end of the discharge process is revealed.  The inlet temperature and flow rate of the water are 25 °C and 0.5 L/min. The legend of the contour is the same as one used for the air cooling BTMS. It is observed that both sides of the module with the liquid cooling system are significantly colder than the other areas, and the hottest area is concentrated in the middle of the module. It is worth mentioning that the temperature distribution on each battery cell in the liquid-cooled module is almost the same. Comparing, it is clear that the liquid cooling system leads to a lower temperature of the module than air cooling due to the larger cooling capacity of water.

It depicts the average temperature distribution of the cells within the liquid-cooled module with at different flow rates and the inlet temperature of 25 °C. To remind the readers, the cells are numbered from the front side of the module toward the rear side. It is seen in this figure that for all flow rates, the maximum and minimum average temperatures correspond to cells 6 and 12, respectively. Moreover, a slight variation is observed in the average temperature between the cells.

In order to pursue the understanding of the temperature distribution within the liquid-cooled module, the temperature of water at the midplane of the upper cooling channel in the cooling plate is displayed. This contour is obtained at the end of the discharge rate for the flow rate of 3 L/min</math>">L/min and inlet water temperature of 25 °C. As it is observed, the water temperature is gradually rising along the flow path. The increase in water temperature leads to a decrease in heat transfer from the cells to water. Accordingly, the heat dissipation from cell 1 to cell 12 gradually decreases on the inlet side, whereas, on the outlet side, the heat transfer reduces from cell 12 toward cell 1. Since there is a very narrow gap between the cells in which the heat transfer is occurring mainly by conduction, the cells located at the ends are colder, and the cells in the middle of the module are hotter due to the accumulation of heat. Among cells 12 and 1, cell 12 is colder because the temperature of the water on both sides of this cell is approximately the same, which is around 26.3 °C. However, there is around 2.5 °C difference in water temperature on different sides of cell 1, which leads to unequal heat dissipation from the sides of this cell. It is worth mentioning again that due to the high cooling efficiency of water, the discrepancy of the temperature of the cells is tiny.

2.Module consists of 5 prismatic cells enclosed in an aluminum enclosure

3.Worst case scenario- cell soaked at 45⁰C ambient with heat generation corresponding to 0.5 C discharge

4.Inlet temperature of 20 ⁰C for both air and liquid.

 Influence of the Discharge Rate:-

The heat generation rate of the battery cell depends on the discharge rate. This section considers various discharge rates as 3-current (3C), 4-current (4C), 5-current (5C) and 6-current (6C) and conducts the CFD calculation respectively. It can be observed that Tmax</math>">Tmax and ΔTmax</math>">ΔTmax both increase as the discharge rate increases. In the BTMS, the temperature difference of the air is not large, so the properties of the air can be treated to be independent of temperature. That means the temperature calculation is decoupled from the velocity calculation. Therefore, the change of discharge rate does not affect the flow field in the BTMS. When the discharge rate increases, the discharge current increases and the heat generation rate increases, which will increase the temperature and the temperature difference of the battery pack.

 

 

Therefore, reducing the discharge rate can help to reduce Tmax</math>">Tmax and ΔTmax</math>">ΔTmax of the battery pack. In reality, when the EV starts or accelerates, more power is needed from the battery pack, and the discharge rate is larger. In these situations, the temperature rise and the temperature difference of the battery pack are large, and thermal runaway is likely to be induced. However, the power consumption of the battery pack usually depends on the properties of the battery cell and the actual power requirement of the EV. The discharge rate cannot be controlled effectively. Therefore, reducing the discharge rate is not a realizable way to reduce the temperature of the battery pack.

 

 

 

It can be observed that the maximum temperature rise (Tmax−T0</math>">Tmax−T0) and the maximum temperature difference ΔTmax</math>">ΔTmax are independent of T0</math>">T0. Note that T0</math>">T0 does not affect the flow field in the BTMS. When T0</math>">T0 is increased or decreased by a certain value, the whole temperature field is increased or decreased by the same value. In reality, the inlet air temperature can be reduced through refrigeration equipment, such as the air conditioner. However, reducing the inlet air temperature can only reduce the absolute temperature, but cannot reduce the temperature difference of the battery pack. Moreover, the introduction of the refrigeration equipment will require additional power consumption. Therefore, it can be concluded that reducing the inlet air temperature is not an effective way to reduce the temperature difference of the battery pack.

 

 

 

 

 

Meanwhile, the velocities in the cooling channels increase as d0</math>">d0 decreases, which will increase the convective heat transfer coefficient on the surface of the battery cells and reduce the cell temperatures. In the parallel air-cooled BTMS, the pressure drop of the airflow is caused by the drag force along the inlet duct and the outlet duct, the drag force along the cooling channels and the local pressure loss due to the swerve of the flow. The reduction of d0</math>">d0 will remarkably increase the drag force along the cooling channel, making the drag force along the cooling channel dominate compared to other pressure losses. Therefore, the discrepancy of the airflow rates among the cooling channels is reduced when d0</math>">d0 is reduced. Meanwhile, the increase of the drag force along the cooling channel increases the power consumption to drive the airflow. T hat compared to the system with d0=3</math>">d0=3 mm, Tmax</math>">Tmax and ΔTmax</math>">ΔTmax of the system with d0=1</math>">d0=1 mm are reduced by 3.8 K and 2.8 K respectively, but Wp</math>">Wp are increased by more than 10-times. Therefore, it is not an effective way to improve the cooling performance of the BTMS through reducing the cell spacing.

 

 

 

In this section, an optimization strategy is used to optimize the plenum widths (w1</math>">w1 and w2</math>">w2). First, the value of w2</math>">w2 is fixed, and the value of w1</math>">w1 is optimized. Then, the value of w1</math>">w1 is fixed and the value of w2</math>">w2 is optimized. Take the optimization of w1</math>">w1 as an example, the detailed steps of the optimization strategy are shown as follows.

• Set w2</math>">w2 as the fixed value and set the initial range of w1</math>">w1 as w11,w12</math>">[w11,w12]. Calculate the middle point of the range w11,w12</math>">[w11,w12] as w13=w11+w12/2</math>">w13=(w11+w12)/2.
• Let w1</math>">w1 equal w11</math>">w11, w12</math>">w12 and w13</math>">w13, respectively. Evaluate the values of ΔTmax</math>">ΔTmax of these three BTMSs using the CFD calculation, respectively, denoting them as ΔTmax,1</math>">ΔTmax,1, ΔTmax,2</math>">ΔTmax,2 and ΔTmax,3</math>">ΔTmax,3.
  If ΔTmax,1<ΔTmax,3</math>">ΔTmax,1<ΔTmax,3,
   let w11=w11,w12=w13,w13=w11+w12/2</math>">w11=w11,w12=w13,w13=(w11+w12)/2, then return to Step 2 and continue the process.
  else if ΔTmax,2<ΔTmax,3</math>">ΔTmax,2<ΔTmax,3,
   let w11=w13,w12=w12,w13=w11+w12/2</math>">w11=w13,w12=w12,w13=(w11+w12)/2, then return to Step 2 and continue the process.
  else
   go to Step 3.
• The range for w1</math>">w1 is reduced through using the following strategy. Let w14=w11+w13/2</math>">w14=(w11+w13)/2 and w15=w13+w12/2</math>">w15=(w13+w12)/2. Evaluate the values of ΔTmax</math>">ΔTmax of the two BTMSs with w1=w14</math>">w1=w14 and w1=w15</math>">w1=w15 using CFD calculation, respectively, denoting them as ΔTmax,4</math>">ΔTmax,4 and ΔTmax,5</math>">ΔTmax,5.
  If ΔTmax,4<ΔTmax,3</math>">ΔTmax,4<ΔTmax,3,
   let w11=w11,w12=W13,w13=w14,ΔTmax,1=ΔTmax,1,ΔTmax,2=ΔTmax,3,ΔTmax,3=ΔTmax,4</math>">w11=w11,w12=W13,w13=w14,ΔTmax,1=ΔTmax,1,ΔTmax,2=ΔTmax,3,ΔTmax,3=ΔTmax,4.
  else if ΔTmax,5<ΔTmax,3</math>">ΔTmax,5<ΔTmax,3,
   let w11=w13,w12=w12,w13=w15,ΔTmax,1=ΔTmax,3,ΔTmax,2=ΔTmax,2,ΔTmax,3=ΔTmax,5</math>">w11=w13,w12=w12,w13=w15,ΔTmax,1=ΔTmax,3,ΔTmax,2=ΔTmax,2,ΔTmax,3=ΔTmax,5.
  else
   let w11=w14,w12=w15,w13=w13,ΔTmax,1=ΔTmax,4,ΔTmax,2=ΔTmax,5,ΔTmax,3=ΔTmax,3</math>">w11=w14,w12=w15,w13=w13,ΔTmax,1=ΔTmax,4,ΔTmax,2=ΔTmax,5,ΔTmax,3=ΔTmax,3.
• If w13−w11</math>">(w13−w11) and w12−w13</math>">(w12−w13) are both smaller than a specified threshold, the optimization process is stopped; otherwise, return to Step 3, and continue the process.

• The temperature rise and the temperature difference of the battery pack increase as the discharge rate increases, but the discharge rate cannot be controlled effectively.
• Reducing the inlet air temperature can reduce the absolute temperature of the battery pack, but cannot effectively reduce the temperature rise and the temperature difference of the system.
• As the inlet airflow rate increases, the maximum temperature of the battery pack is reduced, but the temperature difference of the battery pack and the power consumption to maintain the flow rate are both increased.
• As the cell spacings decrease, the temperature and the temperature difference of the battery pack are both reduced effectively. However, the decrease of the cell spacing increases the power consumption of the system significantly.
• The angles of the plenums remarkably influence the performance of the BTMS. The maximum temperature and the maximum temperature difference of the battery pack can be reduced effectively through optimizing the angles of the plenums, without increasing the total volume and the power consumption of the BTMS.

Code:

clc

clear all

close all

<!-- [if mso]>
<![endif]-->

% Given input values

Qo = [0.005, 0.010, 0.012, 0.015, 0.020];

T_max_air = [329.5, 327.2, 326.5, 325.6, 324.3];

T_max_liquid = [331.2, 328.7, 327.9, 326.8, 325.2];

<!-- [if mso]>
<![endif]-->

% Plot comparative study

figure;

subplot(2, 1, 1);

plot(Qo, T_max_air, '-o', 'LineWidth', 2, 'DisplayName', 'Air-Cooled');

hold on;

plot(Qo, T_max_liquid, '-s', 'LineWidth', 2, 'DisplayName', 'Liquid-Cooled');

xlabel('Discharge Rate (Qo)');

ylabel('Maximum Temperature (Tmax) [K]');

title('Comparative Study of Thermal Performance');

legend('Location', 'Northwest');

grid on;

<!-- [if mso]>
<![endif]-->

% Initialize optimization parameters

w2_fixed = 0.02; % Fixed value of w2

w11_initial = 0.01; % Initial range of w1 start

w12_initial = 0.04; % Initial range of w1 end

threshold = 1e-5; % Specified threshold for stopping optimization

<!-- [if mso]>
<![endif]-->

% Initialize arrays to store optimization results

optimized_w1_values = zeros(length(Qo), 5);

optimized_w2_values = zeros(length(Qo), 5);

<!-- [if mso]>
<![endif]-->

% Loop over different discharge rates (Qo values)

for k = 1:length(Qo)

w11 = w11_initial;

w12 = w12_initial;

 

for i = 1:5

w13 = (w11 + w12) / 2;

% Evaluate ΔTmax for w1 = w11, w12, w13 using CFD calculation

DeltaTmax_1 = cfdCalculation(w11, w2_fixed);

DeltaTmax_2 = cfdCalculation(w12, w2_fixed);

DeltaTmax_3 = cfdCalculation(w13, w2_fixed);

 

if DeltaTmax_1 < DeltaTmax>

w12 = w13;

w13 = (w11 + w12) / 2;

elseif DeltaTmax_2 < DeltaTmax>

w11 = w13;

w13 = (w11 + w12) / 2;

else

w11 = w11;

w12 = w12;

w13 = w13;

end

 

w14 = (w11 + w13) / 2;

w15 = (w13 + w12) / 2;

 

DeltaTmax_4 = cfdCalculation(w14, w2_fixed);

DeltaTmax_5 = cfdCalculation(w15, w2_fixed);

 

if DeltaTmax_4 < DeltaTmax>

w12 = w13;

w13 = w14;

DeltaTmax_1 = DeltaTmax_1;

DeltaTmax_2 = DeltaTmax_3;

DeltaTmax_3 = DeltaTmax_4;

elseif DeltaTmax_5 < DeltaTmax>

w11 = w13;

w13 = w15;

DeltaTmax_1 = DeltaTmax_3;

DeltaTmax_2 = DeltaTmax_2;

DeltaTmax_3 = DeltaTmax_5;

else

w11 = w14;

w12 = w15;

w13 = w13;

DeltaTmax_1 = DeltaTmax_4;

DeltaTmax_2 = DeltaTmax_5;

DeltaTmax_3 = DeltaTmax_3;

end

end

 

% Store optimized values

optimized_w1_values(k, i) = w13;

optimized_w2_values(k, i) = w2_fixed;

end

<!-- [if mso]>
<![endif]-->

<!-- [if mso]>
<![endif]-->

% Plot optimization results

subplot(2, 1, 2);

for i = 1:5

plot(Qo, optimized_w1_values(:, i), '-o', 'LineWidth', 2, 'DisplayName', ['w1 Iteration ', num2str(i)]);

hold on;

end

xlabel('Discharge Rate (Qo)');

ylabel('Optimized w1');

title('Optimization of w1');

legend('Location', 'Northwest');

grid on;

<!-- [if mso]>
<![endif]-->

figure; % Create a new figure for w2

for i = 1:5

plot(Qo, optimized_w2_values(:, i), '-s', 'LineWidth', 2, 'DisplayName', ['w2 Iteration ', num2str(i)]);

hold on;

end

xlabel('Discharge Rate (Qo)');

ylabel('Optimized w2');

title('Optimization of w2');

legend('Location', 'Northwest');

grid on;

<!-- [if mso]>
<![endif]-->

% Function to perform CFD calculation and return ΔTmax

function DeltaTmax = cfdCalculation(w1, w2)

% Perform CFD calculation and return ΔTmax

% Replace with actual CFD calculation based on your study

DeltaTmax = w1 * w2; % Placeholder calculation

end

 

RESULT: