Influence of Reynolds Number and Heat Flux when Cooling an Array of Cylinders
Erik Fagerström
Space Engineering, master's level 2020
Influence of Reynolds Number and Heat Flux when Cooling an Array of Cylinders
Erik Fagerstr¨om
Master of Science Program, Space Engineering
Lule˚a University of Technology
Division for Fluid Dynamics and Experimental Mechanics
February 2020
Preface
This report is the final result of my master thesis for the Master of Science program in Space Engineering, specialisation Aerospace Engineering. The work was done at the division for Fluid Mechanics at Lule˚a University of Technology.
I would want to express my sincere gratitude to my examiner and supervisor Anna-Lena Ljung.
Her guidance, feedback and help has been very appreciated. I would like to thank Staffan Lund- str¨om for the opportunity to do the thesis work. Also, I would thank the Division for Fluid Mechanics for welcoming me during my thesis work.
Last I would like to say thanks to my family Jan, Ulla-Britt, Hugo, Maximilian and towards all my friends who have supported me all the way and never stopped believing in me.
Erik Fagerstr¨om Lule˚a, February 2020
Abstract
When batteries are used, it is of interest to know how high the temperature will become and if the cooling is efficient enough. In this thesis project, the geometry used are batteries placed in an array shape, inside a box. This gives the possibility of assuming porosity conditions when looking at the box from the outside. With the use of the commercial program ANSYS CFX, simulations of convective cooling of cylinders was performed. The goal was to find how large influence the Reynolds number and heat flux from the batteries have on the fluid flow and if there is a possibility to find a prediction equation to determine certain temperatures even before simulating. To determine the influence, the variables Reynolds number, heat flux and type of fluid was varied one at a time. The Reynolds number change the inlet velocity of the fluid and a higher Reynolds number will transport the heat away faster. The heat flux was changed to see the influence on how the temperature of the fluid changed and how the velocity profiles varied due to the amount of power generated by the batteries. Increasing the power output would increase the temperature of the fluid and change the velocity profile difference. Varying the fluid was done to see how the thermal conductivity influenced the fluids behaviour while heated up. The three fluids used in this thesis was water, engine oil and air. There was a clear difference between the three fluids with engine oil being the most efficient while air was considerably worse.
Contents
1 Introduction 1
1.1 Previous work . . . . 1
1.2 Purpose . . . . 1
1.3 Project limitations . . . . 1
2 Theory 2 2.1 Porous media . . . . 2
2.2 Fluid mechanics . . . . 2
2.2.1 Conservation of mass . . . . 2
2.2.2 Conservation of momentum . . . . 3
2.2.3 Conservation of energy . . . . 4
2.3 Computational fluid dynamics (CFD) . . . . 4
2.3.1 Finite volume method (FVM) . . . . 4
2.3.2 Mesh . . . . 4
2.3.3 Wall functions . . . . 5
2.4 Turbulence . . . . 5
2.4.1 k-ε . . . . 6
2.4.2 k-ω . . . . 6
2.4.3 SST - shear stress transport . . . . 7
2.5 Heat transfer . . . . 7
3 Method 9 3.1 Step of work . . . . 9
3.1.1 Reproduce the earlier result . . . . 9
3.1.2 Introducing heat generation . . . . 9
3.1.3 Different fluids . . . . 10
3.1.4 Introducing cylinders . . . . 10
3.2 Simulation setup . . . . 10
3.2.1 Mesh . . . . 10
3.2.2 Boundary condition (BC) . . . . 10
3.2.3 Turbulence model . . . . 10
3.2.4 Fluid setup . . . . 10
3.2.5 Cylinders defined as battery’s . . . . 11
3.2.6 Solver . . . . 11
3.3 Post processing . . . . 11
3.3.1 Temperature relation . . . . 12
4 Results 13 4.1 Mesh study and setup . . . . 13
4.2 Replication of earlier results . . . . 15
4.3 Heat transfer without cylinders . . . . 17
4.3.1 Difference between velocity profile due to temperature difference . . . . 19
4.3.2 Wall heat transfer coefficient . . . . 22
4.3.3 Maximum temperature of the outlet . . . . 23
4.4 Heat transfer with cylinders . . . . 25
4.4.1 Side view of the flow . . . . 27
5 Discussion and Conclusion 29
5.1 Validation . . . . 29
5.2 Mesh study . . . . 29
5.3 Heat transfer without cylinders . . . . 29
5.3.1 Velocity profiles . . . . 29
5.3.2 Wall heat transfer coefficient . . . . 30
5.3.3 Temperature relations . . . . 30
5.4 Heat transfer with cylinders . . . . 30
5.5 Conclusion . . . . 30
5.6 Future work . . . . 31
6 APPENDIX 33 6.1 Mesh studies . . . . 33
6.2 Heat transfer without cylinders . . . . 37
6.2.1 Maximum temperature of the fluid . . . . 37
6.2.2 Average temperature of the outlet . . . . 40
1 Introduction
This thesis work has investigated how flow through a porous media can be used to cool batteries and how the fluid is affected by the heat generated. This is quite interesting at this time period where electronic components are getting smaller. If it is possible to approximate the electronic box as a porous media this could be a way to get a fast idea of the temperature development in the electronic box. In the vehicle industry, electronic cars are getting more popular and improved upon.
It is therefore of interest to know how to most efficient cool down the batteries which are one of the key component in the cars. In this work, the porous media will be an array of cylinders where the cylinders will represent batteries generating heat. With the help of the commercial simulation program ANSYS, Fluid Flow (CFX) simulations will be analysed.
1.1 Previous work
In previous work by Larsson et al. (2018) [9], a study of the turbulence through an array of 60 cylinders (6x10) has been done both with simulations and experiments. In Hellstr¨om et al. (2006) [11], flow through porous media at different Reynolds numbers was studied. Both Larsson et al.
and Hellstr¨om er al., used the definition of particle Reynolds number that is common to use in flow through porous media can be seen as
Rep=
Q ADp
νφ . (1)
Here QA is the volumetric flow rate over the cross-sectional area, Dp is the particle diameter (the diameter of the cylinders), ν is the kinematic viscosity and φ is the bed porosity (φ = 1 − n where n is the the porosity (2)). This equation is a modified version of Erguns equation that originates from Darcys law.
1.2 Purpose
The purpose of this thesis is to study how the fluid velocity profiles change when the fluid is heated up from the batteries and also how cooling works for battery’s with three different fluids when assuming the geometry is a porous media.
1.3 Project limitations
One of the largest limitation is to choose the proper limits of Reynolds numbers and heat generation to look at. There is also a limiting factor of using only three different fluids. The project is furthermore based on simulations, and no experiments are carried out in this work.
2 Theory
2.1 Porous media
There are a lot of different examples of a porous media. Some examples are; soil, porous or fissured rocks, ceramics, filter paper and a loaf of bread. Some less obvious examples are large geologic for- mations of karstic limestone, where the passages in the stone can be far apart and of different sizes [2].
The porosity of a porous media is defined by
n =Uv
Ub
(2) where Uv is the volume of the void space inside the porous medium and Ub is the bulk volume of the whole geometry. It gives a ratio on how much of the geometry that is void (fluid) and how much that is solid.
A term used while working with porous media is permeability, K, [m2] [4], it comes from Darcy’s Law and is used in hydrology. It is used in uniform one directional flow and is a proportional between flow rate and pressure difference in the direction of flow stated as
v = −K µ
∂p
∂x (3)
where ∂P∂x is the pressure difference between the inlet and outlet.
2.2 Fluid mechanics
In fluid mechanics there are a few important physical principles; mass is conserved, Newtons’s second law and the energy is conserved. These principles give the governing equations; continuity equation, momentum equation and energy equation. They together describe how the fluid will behave. All the velocities used in the thesis come from pre-determined Reynolds numbers to be able to study the geometry with different fluids. Using equation (1) that is the particle Reynolds number, the inlet velocity for the system can be determined by v = QA[mm3/s2 = ms].
2.2.1 Conservation of mass
Due to the fact that mass cannot be created or destroyed, the fluid particles that occupy a certain volume can be called material volume V(t). Even when moving and deforming the material volume, it will contain the same fluid particles and nothing comes in and nothing leaves the volume. It gives then that the material volumes surface A(t) moves with the fluid velocity v because the fluid particles will remain inside the material volume V(t) [12]. The conservation of mass for the material volume in the fluid can then be described as
d dt
Z
V (t)
ρ(x, t)dV = 0. (4)
This can better be shown with Reynolds transport theorem d
dt = Z
V ∗(t)
F (x, t)dV = Z
V ∗(t)
∂F (x, t)
∂t dV + Z
A∗(t)
F (x, t)b · ndA (5)
where V*(t) and A*(t) is the control volume respective control area, n is the normal of A*(t) and b is the local velocity of A*(t). To display the mass conservation with Reynolds transport theorem put F = ρ and b = u, this gives
Z
V (t)
∂ρ(x, t)
∂t dV + Z
A(t)
ρ(x, t)u(x, t) · ndA = 0. (6)
By substituting Gauss’ divergence theorem Z Z Z
V
∂Q
∂xi
dV = Z Z
A
niQdA (7)
with the surface integral in equation (5) this leads to the continuity equation, [12]
Dρ Dt =∂ρ
∂t + u · ∇ρ = 0 (8)
for incompressible flow.
If the density is constant the derivative DρDt will be zero due to no change throughout the flow field and leads to the following equation for incompressible flow
∇ · u = 0. (9)
2.2.2 Conservation of momentum
Conservation of momentum is developed from Newton’s second law, with material volume V(t) and surface area A(t) it can be directly stated as follows [12]
d dt
Z
V (t)
ρ(x, t)u(x, t)dV = Z
V (t)
ρ(x, t)gdV + Z
A(t)
f(n, x, t)dA. (10)
Here ρgdV is the body force acting on the fluid element dV. Using Reynolds transport theorem and Gauss’ theorem equation (10) becomes
ρDuj
Dt = ρgj+ ∂
∂xi
(Tij) (11)
that sometimes is called Cauchy’s equation of motion. This relates to fluid-particle acceleration towards the net body ρgiand ∂T∂xij
j is the surface force on the particle, the stress tensor Tij is related to the velocity field [12].
The stress tensor for an incompressible fluid is stated as follows Tij = −pδij+ 2µ(Sij−1
3Smmδij) + µvSmmδij. (12) Substituting equation (12) into equation (11) gives the momentum conservation equation for a Newtonian fluid. For an incompressible fluid it can be developed into
ρDu
Dt = −∇p + ρg + µ∇2u. (13)
Equation (13) is the Navier-Stokes momentum equation for an incompressible fluid [12].
2.2.3 Conservation of energy
The conservation of energy principle is represented by the first law of thermodynamics, it states that energy can transfer from one to another but not be created or destroyed [15]
d dt
Z
V (t)
ρ
e +1
2|u|2
dV =
Z
V (t)
ρg · udV + Z
A(t)
f · udA − Z
A(t)
q · ndA. (14)
In the equation above q is the heat flux and is negative due to the heat being transferred out. Using Reynolds transport theorem and then Gauss’ theorem the conservation of energy equation ends up
as De
Dt = −pDv Dt +1
ρτijSij−1 ρ
∂qi
∂xi
. (15)
2.3 Computational fluid dynamics (CFD)
In fluid mechanics there exists a lot of uncertainty, i.e. it is possible to simulate with DNS (direct numerical solution) but it demands an incredible fine mesh and a large computing capacity. That is why in CFD the solution usally iterated over and over again with certain assumptions on how the turbulence will behave.
2.3.1 Finite volume method (FVM)
For the FVM the conservation equations in integral form are used as a starting point. The solution domain is split up into a finite amount of contiguous control volumes (CV) and the conservation equations is then applied to all CV. All of the CV have a computational node at the centroid where the variable values are calculated. With interpolation the CV values are projected to the surface of the CV [5].
2.3.2 Mesh
The geometry is divided into smaller elements and the general name for it is mesh. When running a simulation the result will be an approximation of reality. How good and accurate the result is defined by how good the mesh is, the more elements in the mesh will generate a finer mesh, the less amount of elements gives a courser mesh. The back side of a finer mesh is that it takes a longer time to simulate due to more elements. In CFD it is needed to choose a balance between accuracy and time. To know how good a mesh is and how approximately wrong it is from reality with regard to numerical error the Richardson extrapolation (RE) method can be used with three different meshes, (course, medium and fine) [3].
The first step is to define a representative cell, mesh or grid size h. It can be found by using equation (16 & 17). N is the number of elements, ∆Viis the volume of cell i and ∆Ai is the area of cell i.
Calculating for two dimensions
h =
"
1 N
N
X
i=1
(∆Ai)
#1/2
(16) and calculating for three dimensions
h =
"
1 N
N
X
i=1
(∆Vi)
#1/3
. (17)
The second step is to determine values for a key variable to be studied and in the equations to follow will be called φ. Next step is to make a grid refinement factor r = hhcourse
f ine , an experimentally determined value for it states that it should be 1.3 [3]. For this step it requires that h1< h2 < h3
and r21 = hh2
1, r32 = hh3
2. It is then required to find the variable p through iterations of equations (18-20) where ε32= φ3− φ2, ε21= φ2− φ1. The iteration starts with assuming q(p)=0.
p = 1
ln(r21)
ln
ε32
ε21
+ q(p)
(18)
q(p) = ln r21p − s r32p − s
(19)
s = 1 · sign ε32 ε21
(20) After p has been found the next step is to determine an extrapolated value for the key variable φ with
φ21ext=rp21φ1− φ2
rp21− 1 . (21)
After the extrapolated value has been found the approximately and extrapolated relative error can be found with
e21approx =
φ1− φ2
φ1
(22) and
e21ext=
φ21ext− φ1
φ21ext
. (23)
2.3.3 Wall functions
Close to the wall, the flow is divided into three parts: Viscous sublayer (y+ < 5), Buffer layer (5 < y+< 30) and Log-law region (30 < y+). In these regions the velocity is scaled with uτ that is the shear velocity. It is related to the y+ value which is the dimensionless wall distance, which can be found by
y+=ρuτy
µ . (24)
For the viscous sublayer u+ and y+follow each other linearly and for the log-law region it scales by u+= 1
κln(y+) + B (25)
where κ is called von Karman constant and B is a constant, usually κ = 0.41 and B=5.2 [5]. This is used to determine the flow behaviour for CFD close to the wall and y+ should be in the viscous sublayer region to effectively capture the small effects of the flow.
2.4 Turbulence
The fluid flow can be laminar, turbulent or in a transition region. In the late 19th century Osborne Reynolds did this for a flow in a pipe and came up with the dimensionless variable Reynolds number (equation (26)) [6]. In the equation v is the velocity, d is the diameter of the pipe and ν is the
kinematic viscosity. For low Reynolds number the flow is laminar. When it reaches a critical value it will turn turbulent
Re = vd
ν . (26)
In ANSYS, different kind of turbulence models can be chosen and the most common will be described in this section [7].
In turbulence, the Reynolds stresses are used and in terms of eddy viscosity it is modeled as in τtij = 2µt
Sij−Snnδij
3
−2ρδij
3 (27)
where µtis the eddy viscosity, ρ is the fluid density, Sij is the mean-velocity strain-rate tensor, k is the turbulent kinetic energy and δij is the Kronecker delta [8]. The eddy viscosity is then defined by a function of turbulent kinetic energy (k ) and turbulent dissipation rate (ε or ω).
2.4.1 k-ε
The k-ε model has given reasonably good results for free-share-layer flows with small pressure gra- dients [8]. The k stand for turbulence kinetic energy and ε rate of dissipation. For the k-ε model the eddy viscosity is defined as
µt=cµfµρk2
ε (28)
where cµ is the model coefficient and fµ is the damping function. This leads to the turbulence transport equations for the k-ε with
∂ρk
∂t + ∂
∂xj
ρuj ∂k
∂xj −
µ +µτ
σk
∂k
∂xj
= τijSij− ρε + φk (29) being the turbulence transport equation and
∂ρε
∂t + ∂
∂xj
ρujε −
µ +µτ
σk
∂k
∂xj
= cε1
ε
kτijSij− cε2f2ρε2
k + φε (30)
being the Energy dissipation transport theorem [8].
2.4.2 k-ω
The k-ω model works well close to boundaries and attached flow and less well in the free stream due to difficult to exercise control over turbulence in the free stream. For the k-ω model the eddy viscosity is defined as follows
µt=ρk
ω. (31)
This leads to the turbulence transport equations for the k-ω with
∂ρk
∂t + ∂
∂xj
ρujk − (µ + σ∗µt)∂k
∂xj
= τtijSij− β∗ρωk (32) being the turbulence transport equation and
∂ρω
∂t + ∂
∂xj
ρujω − (µ + σµt)∂ω
∂xj
= αω
kτtijSij− βρω2 (33) being the Energy dissipation transport theorem [8].
2.4.3 SST - shear stress transport
The SST model is a combination of both the k-ε and k-ω. It uses k-ω close to the wall and k-ε at the free stream taking the advantages of both. For the SST model the eddy viscosity is defined as
µt=
ρk ω
max[1;ΩFa 2
1ω] a1= 0.31 (34)
where Ω is the absolute value of vorticity and F2 is and auxiliary function seen as
F2= tanh
max
"
2
√ k
0.09ωy;500µ ρy2ω
#!2
. (35)
For the two SST transport equations, a blending function F1 is used to model coefficients of the original ω and ε equations. This leads to the turbulence transport equations stated as
∂ρk
∂t + ∂
∂xj
ρujk − (µ + σkµt)∂k
∂xj
= τtijSij− β∗ρωk (36) being the turbulence transport equation and
∂ρω
∂t + ∂
∂xj
ρujω − (µ + σωµt)∂ω
∂xj
= Pω− βρω2+ 2(1 − F1)ρσω2
ω
∂k
∂xj
∂ω
∂xj (37)
being the Energy dissipation transport theorem [8].
In equation (37) the last term is the cross-diffusion term from the ε equation in the transformed ω equation and Pω≈ γρΩ2. The cross diffusion term F1 is stated as
F1= tanh
max
"
max
" √ k
0.09ωy;500µ ρy2ω
#
; 4ρωω2k CDkωy2
#!4
(38)
where CDkω is the cross-diffusion of the k-ω model as follows CDkω= max 2ρσω2
ω
∂k
∂xj
∂ω
∂xj
; 10−20
. (39)
2.5 Heat transfer
Two main phenomena of heat transfer will apply in this thesis work. The first is conduction from the cylinders and second will be the fluid with forced convection. This will work in relation of the cylinders heating up the fluid and the fluid cooling down the cylinders creating a fluid solid interface with a thermal boundary layer.
The conduction is expressed by Fourier’s law of heat conduction, Q˙cond= −kAdT
dx (40)
where k is the thermal conductivity, dTdx is the temperature gradient and A is the Area.
Fourier’s law of heat conduction can also be expressed with heat flux ˙q [W m−2] [14] stated as
˙
q = −k∂T
∂x. (41)
The convection is expressed by Newton’s law of cooling,
Q˙conv= hAs(Ts− T∞) (42)
where h is the convective heat transfer coefficient [W m−2K−1], Tsis the surface temperature, T∞ is the temperature of the fluid far away and Asis the surface area affected by convection [15].
For the case of no cylinders the convection boundary condition is used, i.e. ”the heat conduc- tion at the surface in a selected direction = heat convection at the surface in the same direction”
[14]. For the case with cylinders under steady conditions the energy balance equation ”Rate of heat transfer from the solid = Rate of energy generation within the solid” and it can be viewed with
Q = ˙gV˙ (43)
where ˙g is the heat generation per volume [W m−3].
Equation (43) and (42) can be combined into
Ts= T∞+ ˙gV
hAs (44)
for a cylinder it becomes
Ts= T∞+ ˙gr0
2h (45)
where r0 is the radious of the cylinder.
By knowing the surface temperature, an equation that states the temperature of the inside of a cylinder can be found and is as seen by [14]
T (r) = Ts+ ˙g
4k(r02− r2). (46)
3 Method
The thesis work was done in the commercial simulation program ANSYS CFX that use a hybrid finite-element/finite-volume approach [1]. The work will revolve around simulating the flow in a setup as seen in Figure 1. The first step is to start with achieving similar results as obtained in the experiment work by Larsson et al. (2018) [9]. Acquiring similar results will validate that the geometry and mesh are setup correctly. Next the heat flux is introduced to the cylinder walls to investigate how the flow change when the temperature of the fluid increase along the flow path.
To investigate the influence of material parameters a test with fluids that has a different thermal conductivity and viscosity is tested. The last part is to introduce cylinders that generate heat from a point source.
(a) With the cylinders. (b) Without the cylinders.
Figure 1: The geometry’s used in ANSYS.
3.1 Step of work
3.1.1 Reproduce the earlier result
Validating with previous experiments helps to know if the setup for simulations are correct. Impor- tant is to ensure that the y+ value is below 1 and to make the simulations as computational efficient as possible. The comparison can be done by looking at the velocity profiles to see if the flow behaves the same. The earlier experiments has only been done for a cold flow with no heat generation so the validation will be performed under the same conditions. For reproducing the result water was chosen as the fluid due to lack of thermal knowledge of the working fluid used in earlier experiments [9].
3.1.2 Introducing heat generation
After the initial settings have been figured out to make the simulations work like the experiment, the next step was to introduce heat generation from the cylinders. Initially the right image of Figure 1 was used. It is holes instead of cylinders so a heat flux was introduced along the walls, in the direction out towards the fluid. This was done to be able to study how the flow changed when heat was introduced, before adding the cylinders with an interface. The step is necessary because it can show how the fluid will behave with heat. It also helps to verify if the interface between solid cylinders and fluid is added correct by comparing the velocity profiles with and without cylinders.
3.1.3 Different fluids
The next step is to check how the heat transport changes the behaviour of the fluid depending on the different fluid properties. In this thesis the different fluids used are water, air and engine oil. The point of interest is to investigate how the power generated by the cylinders affect the temperature and velocity profiles in the fluid.
3.1.4 Introducing cylinders
The last step is to introduce cylinders to the geometry. The heat generation will be applied as a point source in the middle of the cylinders. Running simulations will then show how the cylinders are cooled by the flow. It will also show if and how much the fluid velocity differs when using heat flux with or without cylinders.
3.2 Simulation setup
3.2.1 Mesh
The mesh was set up using ANSYS workbench. A mesh study was done for a few different meshes ranging from 2000000 to 7000000 elements. Richardson extrapolation was then used to determine the error in simulated result.
3.2.2 Boundary condition (BC)
• Inlet - Fixed inlet velocity coming from predetermined Reynolds numbers and a fixed temper- ature of the fluid.
• Outlet - Relative Pressure 0 [P a].
• Top, bottom and side walls are put to no slip and smooth wall.
• Cylinders - For non heating it was set to the same BC as above with smooth walls. For heat generation a certain heat flux was set and for cooling of cylinders a point source was set in the middle of the cylinder with a conservative interface flux between the fluid and solid consisting of a thin material representing a shell of a battery.
3.2.3 Turbulence model
To determine the turbulence model for the simulations, test cases were done with different turbulence models. Reynolds number, Re, of Re = 950 with water as a fluid was chosen for the test case, and the simulations was performed and then compared at a cross-section of the geometry. The criteria for the needed model is to be able to solve heat and velocity close to the cylinders in the boundary layer and also give to an accurate representation of the free stream velocity between the cylinders.
3.2.4 Fluid setup
• Water - Used the predefined material ”Water” from ANSYS Engineering data sources.
• Engine oil - Used the predefined material ”Engine Oil” from ANSYS Engineering data sources.
• Air - Used the predefined material ”Air Ideal Gas”, in the material the transport properties were changed to Sutherland’s Formula for dynamic viscosity and thermal conductivity. They are found in Table 1 [10].
Table 1: Sutherland’s constants for air.
Dynamic viscosity:
Thermal Conductivity:
Tref 273.1K Tref 25C
µref 1.716 · 10−5N s m−2 kref 0.0261 W m−1 K−1
S 111K S 111K
3.2.5 Cylinders defined as battery’s
When the cylinders were implemented in the geometry they were set to Manganeese and the thin material for the conservative interface flux was put to be Zinc. It would represent a Zinc Carbon battery (AA battery) [13]. The point source would represent a carbon rod.
3.2.6 Solver
In the solver, an expert parameter was inserted in the physical model. Tbulk for htc was activated and set to the inlet temperature of the fluid, 298.15K (25◦C). This was done so that even if the fluid is increasing in temperature the wall heat transfer coefficient would be calculated based on the inlet temperature.
3.3 Post processing
After a simulation is done, six outputs are investigated such as
• A normalized velocity profile from the bottom to the top of the geometry at both the outlet and between the two last cylinders for the middle flow path. The lines can be seen in Figure 2 where the inlet is at the left side. This is used to make a graph showing how heating changes the velocity. By using a normalized velocity (velocity divided by the inlet velocity) comparison between simulations with different initial conditions is valid.
• The maximum temperature of the fluid. This is investigated to determine the temperature closest to the cylinder walls and also an indication if the fluid will start to boil.
• The heat transfer coefficient. This is investigated to study how the wall heat transfer coefficient will change depending on different Reynolds numbers or amount of power generated by the batteries.
• Average outlet temperature. This is investigated to get an idea of how much warmer the fluid has become while passing the cylinders.
• The maximum outlet temperature. This is investigated to know if the fluid have passed the boiling point and is done to know if it that specific scenario is possible or not. If it is above the boiling point it shows that the fluid would change aggregation form and complications could arise in the system with pressurisation.
• The maximum and minimum temperature of the cylinders. This is investigated to see how the cylinders in the array will affect each other while cooling.
Figure 2: The two yellow lines at the right side is the measuring lines.
3.3.1 Temperature relation
By knowing the temperature for different Reynolds numbers and heat generation, an equation to calculate the temperature depending on both flow and heat could be found. The first step is to make a graph of temperature and power, P. From the graph a curve fitting could be done, in this thesis a linear fitting is used. A relation can be seen as
T emperaturef luid= θ · P + m (47)
where m is the inlet temperature.
By doing it for different Reynolds numbers it gives a different constant θ. Next is to create a new graph with Reynolds number and θ. Using curve fitting, an equation can be found where θ depends on the Reynolds number, θ = f (Re). This gives a final equation
T emperaturef luid= f (Re) · P + m (48)
where the temperature depends on Reynolds number and the power produced by the cylinders.
The temperature equations could experimentally be found this way and be done for the differ- ent places of interest, in this case maximum/average outlet, maximum of fluid and temperature of the cylinders.
4 Results
The result consist of three parts, starting with validation of the setup. Secondly comes the part of heat transfer without cylinders to investigate what influence the Reynolds number, heat flux and fluid have on the flow. In the last part the cylinders are included to show the difference in cooling.
4.1 Mesh study and setup
For the mesh study permeability is chosen as key variable, and the difference between the meshes can be seen in Table 2.
Table 2: Different meshes tested for Re = 465.
Number of elements Max face size [10−3 m] Permeability [10−5 m2] Max y+
2237825 1.6 0.20704732337894 1.61575
3199221 1.35 0.21353114569112 1.69817
3921516 1.2 0.21624297881125 1.62936
5162217 1.05 0.220748712245734 1.66983
6018840 0.98 0.223152962113225 1.61276
6910703 0.92 0.22494724022041 1.68441
From the data Figure 3 is constructed where an extrapolated value for the permeability is deter- mined along with an approximate error for the mesh and the extrapolated permeability. Richardson’s extrapolation were used on the three finest meshes. The values of r21and r32were 1.1414 respectively 1.1569 and determined to be acceptable close to the experimental value of 1.3. Last the approximate relative error (equation (22)) is 1.089% and the extrapolated error (equation (23)) is 2.479%.
Figure 3: Permeability plotted against max face size from Table 2.
The mesh chosen for the simulations was the mesh with approximately 5 000 000 elements due to high accuracy and shorter simulation time and the mesh can be seen in Figure 4. All the other meshes from Table 2 can be seen in the Appendix section 6.1.
Figure 4: Re = 465, mesh with approximately 5 000 000 elements.
4.2 Replication of earlier results
In Figure 5, the flow can be seen for Re = 45, 465 and 950. The case of Re = 45 has a oval shape in the y-z plane where it is faster in the middle of the free stream flow. For Re = 465 the flow has developed into a more of bone shaped velocity profile in the y-z plane. In the narrow area it is clearly visible that the velocity is fastest at the center of the free stream flow. Closer to the top and bottom the fast velocity area expands towards the cylinder walls as seen in image (b) Figure 5. The maximum normalized velocity is furthermore lower then the case of Re = 45. For the case of Re = 950 it is resembling Re = 465 with the bone shape but the maximum normalized velocity is lower for Re = 950. Comparing this with the results from earlier experiments [9] Re = 45 have the same kind of round shape for both the wide and small area. For Re = 465 have the same rectangular shape with a small tendency of narrowing at the middle and the small area shows the resemblance of the bone shape. For Re = 950 it was less resemblance since the red areas went all the way to the top and bottom edges in experiments by Larsson et al. (2018), but both still has the narrowing in the middle of the flow.
(a) Re = 45 (b) Re = 465
(c) Re = 950
Figure 5: Replication of Larsson et al. (2018) Figure 5.
In Figure 6, the velocity profiles between cylinders in the middle of the flow can be seen. The smallest passage is between the two last cylinders and the widest passage is at the outlet. The shapes resembles the experiments but a bit higher values for the normalized velocity is observed.
Figure 6: Replication velocity profile of Larsson et al. (2018).
4.3 Heat transfer without cylinders
This section will start with a general view of the flow when the heat flux is set according to power = 50W. Then the different subsections will be presented. The first subsection will show a comparison between the normalized velocity profiles for the flow when there is a heat flux subtracted to when there is no heat flux. This shows how the velocity will change at the end of the geometry, between the top and bottom for different fluids. Next subsection shows how the wall heat transfer coefficient will change for different fluids. The last subsections describes the path to find a prediction equation for maximum outlet temperature. Maximum temperature of fluid and average outlet temperature can be found in the Appendix 6.2.
In Figure 7-9 temperature, velocity and turbulent kinetic energy is shown for water with Re = 45 (left image) and Re = 950 (right image). From Figure 7 it can clearly be seen that for a higher Reynolds number the heat is transported away faster. With a lower Reynolds number the heat is starting to accumulate for each cylinder. Figure 8 shows the normalized velocity of the fluid.
It can be seen that for a lower Reynolds number the velocity increase along the flow, compared to the inlet velocity is higher than for higher Reynolds numbers. This shows that the higher the inlet velocity the smaller the normalized maximum velocity will appear in the flow. Figure 9 shows
how the turbulence kinetic energy increase with Reynolds number. The inlet velocity is set with a medium turbulence (Intensity=5%) setting in ANSYS. For a low Re the flow turbulent kinetic energy is more negligible for the free stream indicating laminar flow at Re = 45. When the Re is higher the turbulent kinetic energy has developed into the free stream and at the cylinders.
(a) Re = 45. (b) Re = 950.
Figure 7: Temperature plane of the geometry at half the height of the geometry with water as the fluid and the inlet on the right side.
(a) Re = 45. (b) Re = 950.
Figure 8: Normalized velocity plane of the geometry at half the height of the geometry with water as the fluid and the inlet on the right side.
(a) Re = 45. (b) Re = 950.
Figure 9: Turbulent kinetic energy plane of the geometry at half the height of the geometry with water as the fluid and the inlet on the right side.
4.3.1 Difference between velocity profile due to temperature difference
Next the velocity profiles from the top to the bottom of the geometry is compared from the center of the flow at two places. The ”Narrow” area is between the two last cylinders and the ”Wide” area is at the Outlet of the geometry as seen in Figure 2. To display how addition of heat influence the fluid flow, Figure 10-12 shows the velocity profile without heat flux subtracted to the velocity profile with heat flux for Re = 45, 465 and 950, respectively. Figure 10 shows the difference between Re for water. For Re = 45 there is a small difference (10−4) where the fluid in the wide part becomes slower and the narrow part becomes faster. Figure 11 shows the difference for engine oil, where the velocity only differs for 50 watt for the three different Reynolds number. When the Reynolds number increase, the velocity is increasing at the wide area and the largest difference between when heat is included or not is observed at the top and bottom of the geometry. The narrow area remains in the difference span of -0.05 to 0 when heat is included or not. For the wide area the difference between when heat is included or not is increasing up to approximately 0.35 close to the top and bottom.
Figure 12 shows the difference in air. When the power is increased it generates a larger normalized velocity difference between the two case. The difference for air compared to water and engine oil is that the air is decreasing in difference while water and engine oil shows a larger difference for higher Reynolds number. This could be from the fact that for higher Re the maximum temperature for the different power cases is lower.
(a) Re = 45. (b) Re = 465.
(c) Re = 950.
Figure 10: Velocity profile illustrating the difference between simulations of water (heat-no heat).
(a) Re = 45. (b) Re = 465.
(c) Re = 950.
Figure 11: Velocity profile illustrating the difference between simulations of engine oil (heat-no heat).
(a) Re = 45. (b) Re = 465.
(c) Re = 950.
Figure 12: Velocity profile illustrating the difference between simulations of air (heat-no heat).
4.3.2 Wall heat transfer coefficient
In Figure 13 the wall heat transfer coefficient, h can be seen for different Re and fluids depending on the power. In the left image it shows h for water and engine oil. The the value of h does not show any change depending on the amount of power generated for the different Re. This could be due to the temperature does not change that much and thermal conductivity, k, is set to a constant value in the simulations. The right image shows that there is a change in the value for air depending on the amount of power generated. For air the temperature change a lot and the k is not constant due to the use of Sutherland’s formula. It can also be seen that the scale of h for the two different images is quite different. This could be from h depending on the temperature difference in the denominator.
Air has a temperature that is a lot higher then water and engine oil resulting in h have a lower value.
(a) Wall heat transfer coefficient for engine oil and water. (b) Wall heat transfer coefficient for air.
Figure 13: Wall heat transfer coefficient for different Reynolds depending on power.
4.3.3 Maximum temperature of the outlet
In this part, the maximum outlet temperature of the fluid is investigated. Figure 14 shows the maximum outlet temperature for different amount of power. The figure is split up into two graphs due to the fact that air has a lot higher temperature while engine oil and water are closer together in temperature range.
(a) Water and engine oil, the upper green line represents boiling temperature for water and the lower green line il- lustrates the inlet temperature.
(b) Air.
Figure 14: Maximum outlet temperature of fluid.
From Figure 14, it can be seen that temperature depending on watt is a linear relation. For engine oil and water the fluids properties are constant but for air the thermal conductivity and viscosity are varied with Sutherland’s formula. In Table 3 the linear equations for temperature can
be seen for engine oil of different Reynolds number. The constant for the different Reynolds are all 25 that is the inlet temperature, this confirms that linearization in Figure 14 is a valid approximation for engine oil.
Table 3: The linear equations for temperature, depending on power for different Reynolds numbers of engine oil, equation (47).
Reynolds number Temperature (Tmax=θ· P+25)
45 Tmax=0.02435·P +25
350 Tmax=0.002476·P +25
465 Tmax=0.001302·P +25
950 Tmax=0.0006755·P +25
1500 Tmax=0.000302·P +25
2000 Tmax=0.0002·P +25
In Figure 15, the left image shows the curve fit for θ (from equation 47) plotted against the Reynolds number for engine oil. The curve fit for all the data points shows a good correlation in the beginning but after Re = 465 the fitting starts to oscillate. A shape-preserving interpolant shows that for the data points a linear fitting is a better option for Reynolds above Re = 465. This could relate to the flow going from a laminar flow to a turbulent flow. The right image shows the predicted temperature for different watt.
(a) Reynolds numbers relation to the constant θ. (b) With the predicted average outlet temperature of the fluid at the cylinder wall.
Figure 15: Validation of prediction of temperature depending on power and Reynolds number.
The equation found that gives a relation for the average outlet temperature of engine oil can be seen as
Tmax outlet= P · (9.14 · 10−14· Re4− 3.108 · 10−10· Re3+ 3.644 · 10−7· Re2− 0.0001765 · Re + 0.03158) + 25 Re ≤ 465 Tmax outlet= P · (−4.566 · 10−7· Re + 0.00107) + 25 Re > 465.
(49)
Equation (49) had to be split up into two and started to oscillate for higher Reynolds number. A forth degree polynomial fitted best for Reynolds below Re = 465 and a linear fit for Reynolds above Re = 465.
From left image of Figure 14 the water has a linear increase depending on watt and Table 4 shows that it could be a valid approximation because the constant is approximately 25. For air it can be seen in the right image of Figure 14 that it does not increase linearly and Table 4 shows that it is not linearly dependent since to the constant is far from 25.
Table 4: The linear equations for temperature depending on power for different Reynolds numbers of water and air, equation (47).
Water
Reynolds number Temperature (Tmax=θ· P+25)
45 Tmax=4.58· P +24.93
465 Tmax=0.5332· P +24.99
950 Tmax=0.2133· P +25.03
Air
Reynolds number Temperature (Tmax=θ· P+25)
45 Tmax=429.5· P -0.198
465 Tmax=80.38· P +132.5
950 Tmax=49.53· P +114.2
4.4 Heat transfer with cylinders
In this section the cylinders are added to the geometry, a point source of heat is added in the middle of each cylinder. In Figure 16-18 the inlet is at the right side and the point heat source is generating 50W each and the fluid is water. There is two temperatures showing in the images. One is the temperature inside the cylinders and one is the temperature of the fluid.
It can be seen in Figure 16 that temperature is starting to increase behind the cylinders in the same way as in the left image of Figure 7. Looking at the inside of the cylinder it can clearly be seen that the first row is colder at the edges compared to the last row. In Figure 17, with a faster flow the heat is not influencing the cylinders behind as much. The cylinders are getting cooled and that can be seen with all the rows of cylinders have around the same temperature in the middle. In Figure 18, the water is not significantly heated up and that gives an even cooling to the cylinders.
All the cylinders are evenly cooled and it can be seen that they are warmest in the middle.
Figure 16: Temperature of the water and the cylinders for Re = 45 and 50 watt.
Figure 17: Temperature of the water and the cylinders for Re = 465 and 50 watt.
Figure 18: Temperature of the water and the cylinders for Re = 950 and 50 watt.
4.4.1 Side view of the flow
The side view at the middle of the flow where the inlet is at the left side in Figure 19. It can be seen that the point source of heat is affecting the velocity of the fluid,displayed by an increase of velocity in the middle of the flow.
Figure 19: Side view of the flow for water at Re = 45 and 50 watt with cylinders.
