Numerical simulation of plain fin-and-tube heat exchangers:
CFD simulation of a fin-and-tube heat exchanger:
Heat transfer, fluid flow, and turbulence model analysis using 3D open-source CFD code
Master of Science Thesis
Computational Chemical Engineering
Group for Chemical Fluid Flow Processes
Aalborg University Esbjerg
Neils Bohrs Vej 8
DK-6700 Esbjerg
Abstract
Student:
Anna Margrete Hansen
Project Period:
May – November 2008
AAUE Supervisors:
Bjørn H. Hjertager
Tron Solberg
Vestas Aircoil A/S
Supervisor:
Claus Ibsen
Issues: 6
Pages: 89
Enclosed: CD-ROM
Preface
This Master Thesis project was carried out at Aalborg University – Esbjerg (AAUE) to fulfill the 10th semester requirements leading up to the Master’s Degree in Chemical Engineering with a specialisation of Computational Chemical Engineering. The purpose of this report is to document work performed during the project period, and to demonstrate an understanding and competence in using material covered in relevant courses.
This report describes the use of Computational Fluid Dynamics (CFD) to simulate air flow and heat transfer occurring in a two-row tube-and-fin heat exchanger. A finite volume model was created and used for CFD simulations. Experimental work found in the literature validated the CFD results. CFD concepts and simulation methods are briefly described to cover basic theory and practical use as related to this project.
The project report includes three parts: 1) the main report, 2) an appendix, and 3) a CD-ROM. The report consists of:
• Nomenclature and Introduction
• CFD and heat exchangers
• Summary of CFD simulations carried out.
• Description of mathematical theory and numerical solutions.
• Summary of results.
The final summary includes discussion and conclusion sections. References conclude the report, and are listed with square brackets placed after a section (refers to the entire section) or when applicable, directly after a specific section or illustration. The Appendix follows and includes details of CFD computations and the Excel files created for the flow and heat transfer calculations in this project. The CD-ROM includes all CFD simulations, MS Excel calculations, and data files created for this project and a copy of the report and appendix.
I would like to thank Professor Bjørn Hjertager for his enthusiasm in teaching the related courses, as well as Professor Tron Solberg and Ph.D. student Rolf Hansen for their help and guidance with this project. I would also like to thank Professor Birgit Storm for her support and encouragement throughout my time as a master’s degree student at AAUE.
Esbjerg, the 2nd of December, 2008
_________________________________
Anna Margrete Hansen
Table of Contents
Abstract i
Preface iii
Nomenclature vii
1 Introduction 1
1.1 Previous Research and Experiments 2
1.2 Problem Formulation 4
1.3 Project Outline 5
2 Model Description 7
2.1 Governing equations and numerical schemes. 8
2.2 Computational Domain and Boundary Conditions 9
2.3 Numerical mesh 11
2.4 Performance parameters 14
2.4.1 Dimensionless Groups 14
2.4.2 Heat Transfer and Efficiency 15
2.4.3 Pressure Drop 19
3 Computational Fluid Dynamics 23
3.1 CFD Computational Tools 23
3.2 CFD Governing Equations 25
3.3 Turbulence Modelling 26
3.3.1 k-εpsilon Turbulence Model 29
3.3.2 SST k-omega Turbulence Model 33
3.4 Finite Volume Method and Differencing Schemes 35
3.5 Solution Algorithms 37
4 Preliminary Results and Observations 41
4.1 Results of Grid Independence Test 41
4.2 Pressure Tracer in Transient Case 42
4.2 Characteristics of Flow 43
4.3 Characteristics of Heat Transfer 46
4.4 Summary of Preliminary Observations 48
5 Numerical Results and Discussion 49
5.1 Friction Factor 49
5.2 Colburn j-Factor 51
5.3 Summary of Results 52
6 Discussion 53
7 Conclusion 55
8 References 56
Appendix 60
A1 CFD Computations 62
A1.1 CFD Governing Equations 62
A1.2 Reynold’s Averaged Navier Stokes (RANS) Equations 67
A1.3 Finite Volume Method and Finite Differences 70
A2 Pressure drop and friction factor calculations 74
A3 Fin efficiency calculations 80
A4 Inside tube (water flow) heat transfer coefficient 82
A5 Colburn j-factor calculations 84
B Contents of Enclosed CD-ROM 88
Nomenclature
|A |Area |[m2] |Greek Letters |
|Ao |Total surface area |[m2] |δ |Tube wall thickness [m] |
|Ato |External tube surface area |[m2] |ε |HE effectiveness |
|C |Heat capacity rate |[W/K] |η |Fin efficiency |
|Cp |Specific heat |[J/kg K] |ηo |Surface effectiveness |
|Dc |Fin collar outside diameter |[m] |μ |Dynamic viscosity [kg/ms] |
|Di |Inside tube diameter |[m] |ρ |Density |
| | | | |[kg/m3] |
|Do |Tube outside diameter |[m] |σ |Contraction ratio of |
| | | | |cross-sectional area |
|Dh |Hydraulic diameter |[m] | | |
|f |Friction factor | | | |
|Gc |Mass flux of air based on |[kg m2/s] | |
| |minimum flow area | |Subscripts |
|H |Fin spacing |[m] |1 |Air side inlet |
|h |Heat-transfer coefficient |[W/m2 K] |2 |Air side outlet |
|j |Colburn factor: Nu/RePr1/3 | |air |Air side |
|k |Thermal conductivity |[W/m2K] |ave |Average value |
|Kc |Abrupt contraction | |b |Base surface |
| |pressure-loss coefficient | | | |
|Ke |Abrupt expansion | |i |Tube side |
| |pressure-loss coefficient | | | |
|L |Depth of heat exchanger |[m] |in |Inlet |
| |in airflow direction | | | |
|N |Number of tube row | |f |Fin surface |
|[pic] |Mass flow rate (mass*velocity) |[kg/s] |m |Mean value |
|NTU |Number of transfer units: UA/Cmin | |min |Minimum value |
|Nu |Nusselt number: h/(k/Dh ) | |max |Maximum value |
|Δp |Pressure drop |[Pa] |o |Total surface |
|Fp |Fin pitch |[m] |out |Outlet |
|Pl |Longitudinal pitch |[m] |water |Water side |
|Pr |Prandtl number: μ Cp/k | |w |Wall of tube |
|Pt |Transverse pitch |[m] | | |
|[pic] |Heat-transfer rate |[W] | | |
|[pic] |Max. possible [pic]: Cmin(Tw,in-Tair,in) |[W] | | |
|rc |Tube outside radius |[m] | | |
|Re |Reynolds number: (ρ*V*D)/μ | | | |
|Req |Equiv. radius for circular fin |[m] | | |
|r |Tube inside radius |[m] | | |
|t |Fin thickness |[m] | | |
|T |Temperature |[°C] | | |
|U |Overall heat-transfer coefficient |[W/m2 K] | | |
|XL |[pic] |[m] | | |
|XM |Pt /2 |[m] | | |
1 Introduction
Vestas Aircoil A/S produces compact tube-and-fin heat exchangers for ship motors, as well as other types of heat exchangers and cooling towers (Figure 1). The heat exchanger cools heated, compressed air from the motor with cooling water. Fins are used to increase heat transfer area on the air side, since the air has the largest influence on the overall heat transfer resistance.
A test rig has previously utilized Vestas Aircoil A/S in order to empirically determine heat transfer and pressure drop correlations and thereby determine relevant heat transfer parameters for the heat exchangers. Building and testing prototype heat exchangers are processes which are both expensive and time-intensive, and therefore CFD is an attractive way to develop new heat exchangers in the future.
Vestas Aircoil A/S is investigating the possibilities for developing new heat exchangers based on CFD, since it is less expensive than experimental tests and can give better insight into the local flow and heat transfer characteristics occurring within the heat exchanger.
Open-source CFD code OpenFOAM is used for this project, since other commercial codes such as Fluent and Ansys CFX require expensive license fees which are so high as to be prohibitive for most small- and medium-sized companies to justify the cost. With open-source code, the only costs are the computer hardware and the engineer’s time used for setting up the case.
This project involves building a model of a fin-and-tube heat exchanger geometry using open-source software, creating a suitable mesh, setting up the cases (choosing solvers, numerical solution methods, etc.), making the CFD calculations with OpenFOAM, and comparing results to known experimental data. Since the data available from previous Vestas Aircoil A/S testing is confidential and not necessarily comprehensive enough for CFD validation, experiments done on fin-and-tube heat exchangers and reported in the literature are used for validation.
The following subsections describe the project in more detail. First, a summary of other research in the heat-transfer field related to tube-and-fin heat exchangers is presented to put this study into perspective with the other work available in the open literature. Section 1.1 also includes a description of the experimental research used as a basis and validation for this study [Wang et al., 2006]. The problem formulation presents the actual heat exchanger configuration simulated for this project. Finally, the project outline presents specific activities during the project period and an outline of this report.
1 1.1 Previous Research and Experiments
This study is a CFD simulation of the heat transfer and fluid flow of a two-row heat exchanger previously tested experimentally and reported in the literature. The research article chosen for this purpose [Wang et al., 2006] describes experimental results from 15 heat exchangers of different geometrical parameters such as number of tube rows, fin spacing and fin thickness. In the experiments, heat exchangers were tested with an induced flow open wind tunnel, and results presented in graphs of friction factor and Colburn j-factor against Reynolds numbers. The next sections describe other related research done in the past, followed by a more detailed description of the experimental work carried out by Wang et al. (1996).
Previous Research
There has been a variety of work carried out to study tube-and-fin heat exchangers. Previously, much of the research was experimental, as theoretical prediction of heat-transfer coefficients is complex due to the airflow pattern occurring in the exchangers. However, more recently there have also been more numerical studies carried out, as the use of CFD is becoming more widespread.
Previous experimental work has focused on obtaining data for analysis, future design, or to create or verify empirical correlations. Wang et al. (2006) carried out experiments to test the dependence of heat transfer and pressure drop on the geometrical parameters of 15 different tube-and-fin heat exchanger samples, to determine how fin spacing, fin thickness and number of tube rows affect the Colburn j-factor and friction factor. This is the study used for validation in this project (with experimental details described later in this section). Kayansayan (1994) characterized heat-transfer in tube-and-fin heat exchangers for 10 configurations for Reynolds numbers ranging from 100 to 30,000, with the Reynolds number characteristic dimension being the tube collar thickness, and studied in particular the effect of fins on heat transfer. Yan and Sheen (1999) made a study to compare plate, wavy, and louvered fin-and-tube heat transfer and pressure drop characteristics using different evaluation methods for the air side performance. Infrared thermographic experiments have been carried out to characterize the temperature distribution on the fins and calculate fin local convective heat transfer coefficients of staggered and in-line tube-and-fin heat exchanger arrangements [Ay et al, 2002]. Correlations for both staggered and in-line heat exchangers have been developed to predict the friction factor and Colburn j-factors [Gray and Webb, 1986] [McQuiston, 1978]
Numerical studies have included simulations based on finite differences and CFD. A 2D second-order finite differencing analysis on one-and two-row tube-and-fin heat exchangers has been carried out to compare heat transfer and pressure drop between exchangers containing circular and elliptical tubes [Rocha et al.]. Analytical methods for determining fin efficiency have been compared using 2D SimTherm® software for numerically solving the heat conduction equation. This study by Perrotin and Clodic (2003) included comparisons between the commonly reported method for determining fin efficiency (which utilizes modified Bessel functions of the first and second kind) and the more simplified versions for the same calculation: the equivalent circular fin and sector methods. Finite differencing has been used for estimating the heat transfer coefficient on the fins [Chen et al. 2006]. Tao et al. (2007) developed a 3-D code to study fin-and-tube heat exchangers, using a body-fitted coordinate system based on the Poisson equation.
Three studies were found in the literature search which used CFD to simulate flow and heat transfer in tube-and-fin heat exchangers. All of them used the Fluent CFD program and were directed at comparing the heat transfer and pressure drop of heat exchangers with different geometrical characteristics [Erek et al., 2005] [Sahin et al., 2007] [Tutar and Akkoca, 2002]. There were no articles in the literature found regarding the use of open-source CFD software OpenFOAM to simulate tube-and-fin heat exchangers. However, there has been work done by Mangani et al. (2007) to study the development and validation of the CFD computational code used in the OpenFOAM software. It was determined in this study that the OpenFOAM libraries accurately reproduced flow conditions, a conclusion which was verified with both experimental data and commonly used commercial software.
The literature review has shown that virtually no CFD simulations on tube-and-fin heat exchangers using OpenFOAM have been published in the open literature. Furthermore, the CFD studies found all dealt with the effect of geometrical parameters on the heat transfer and pressure drop characteristics. In this study of tube-and-fin heat exchangers, the simulation results from just one heat exchanger geometrical configuration: a two-row, staggered tube-and-fin arrangement, simulating pressure drop and heat flow for a range of Reynolds numbers from approximately 330 to 7000. However, for this study, the CFD simulations are carried out using the open-source CFD software OpenFOAM, and different flow models are used for simulations: a laminar flow model and turbulence models k-epsilon and SST k-omega.
Experimental Work
The experiments carried out by Wang et al. (1996) were conducted using a forced draft wind tunnel (Figure 2). An air straightener was used to keep flow moving in the x-direction, an 8-thermocouple mesh was placed in the inlet and a16-thermocouple mesh in the outlet (locations of which determined by ASHRAE recommendations [ASHRAE, 1993]. All equipment for data acquisition (thermocouples, pressure transducer, airflow measurement station, and flow meter were checked for accuracy prior to running the experiments.
Water at the inlet was held at 60ºC, air flow velocities were tested in the range from 0.3 m/s to 6.2 m/s. Energy balances were monitored during the experiment for both the hot- and cold-side and reported to be within 2 %. The uncertainties for the primary measurements (mass flow rate for air and water, pressure drop, and temperature of the water and air) were very small and therefore these measurements can be assumed to be accurate.
Calculated values for the friction factor f and the Colburn j-factor, however, have higher estimated uncertainties at the lower Reynolds numbers. The calculated friction factor f has an uncertainty of ± 17.7 % at Reynolds number 600 (± 1.3 % at Reynolds number 7000). The Colburn j-factor has an estimated uncertainty of ± 9.4 % at Reynolds number 600 (± 3.9 % at Reynolds number 7000).
[pic]
1) inlet
2) honey cone straightener
3) developing section
4) T/C inlet temperature measuring station
5) pressure tap (inlet)
6) test unit
7) pressure tap (outlet)
8) T/C outlet temperature measuring station
9) code tester for measurement of air flow rate
10) setting means
11) nozzle pressure tap (inlet)
12) nozzle pressure tap (outlet)
13) multiple nozzles plate
14) setting means
15) variable exhaust fan system
16) discharge
17) water pump
Figure 2. Illustration of experimental set-up for heat exchanger testing [Wang et al., 1996].
2 1.2 Problem Formulation
For this project, the geometrical parameters for a two-row heat exchanger based on experimental research [Wang et al., 2006] are used to build a CFD model, and results read from the graphs (friction factor and Colburn j-factor against Reynolds number) in the article are used to validate the results of the CFD simulations. The parameters of interest: friction factor f and Colburn j-factor are widely used in industry to characterize pressure drop and heat transfer, respectively, and thereby determine heat exchanger performance and suitability for specific duties. Determining and using these parameters for performance prediction is part of the heat exchanger design process. More detailed descriptions of the performance parameters are given Section 2.4: Performance Parameters.
The two-row fin-and-tube heat exchanger studied has a staggered tube arrangement, as illustrated in Figure 3. Analyzing flow and heat transfer using CFD can make calculations to predict heat
exchanger performance. However, it is not possible to perform CFD simulation on the entire heat exchanger, due to the large number of volumes and calculations required. Therefore, a small section of a heat exchanger consisting of one channel of air between two fins, with the air flowing by two tubes is modelled for this project (illustrated and described in Sections 2.2 and 2.3). Simulations of the air flow through this passage are carried out, while relevant characteristics of the air flow are sampled and averaged at the inflow, minimum free-flow area(s), and outflow. The characteristics sampled are: flow velocity (in all three directions: x, y, and z), temperature, pressure, and turbulence model parameters k, epsilon, and omega. These measurements are then used for calculating relevant performance parameters such as pressure drop, friction and Colburn factors, heat transfer rate, Reynold’s number, etc., (described in Section 2.4: Performance Parameters).
1.3 Project Outline
For comparison to the graphs in the validation research article [Wang et al., 2006], flows of ten different Reynolds numbers (based on tube collar diameter and minimum free-flow velocity) are simulated ranging from approximately 330 to 7000 with inlet frontal air velocities ranging from 0.3 to 6.2 m/s. Water flows at 60º C through the tubes and cold air through the fins.
For determining which turbulence model most accurately represents heat exchanger flow and heat transfer at the different flow regimes (laminar, transitional, and turbulent), three flow models were chosen for the simulations. They are: laminar, k-epsilon turbulence model, and Menter SST k-omega turbulence model. Two steady-state OpenFOAM solvers were used for the 60 simulations (2 solvers * 10 velocities * 3 turbulence models): ‘simpleFoam’ (flow calculations only), and ‘rhoSimpleFoam’ (both flow and temperature simulations). Finally, to investigate the possibility of transient patterns occurring in the flow, one transient simulation in rhoTurbFoam was carried out.
All computational work is carried out using open source software: pre-processing software Salomé for geometry and meshing, OpenFOAM CFD solver, and ParaView post-processing for visualization. Pressure drop and heat transfer results and comparisons of the different turbulence models and solvers are reported and discussed.
After this introduction section, the usual relevant topics are covered: a model description of the heat exchanger to include governing equations, computational domain, and mesh. Then the performance parameters related to pressure drop and friction factor are presented, which is followed by a section on computational fluid dynamics (CFD), including equations, turbulence models, and solution algorithms. Finally, the results are presented in two sections: Preliminary observations and then final results. The paper concludes with discussion and conclusion sections.
2 Model Description
This section describes the heat exchanger model and performance parameters used in characterizing heat transfer and pressure. The model heat exchanger for this project is presented, and information about the heat exchanger, fin-and-tube efficiency, pressure drop, and the dimensionless groups used in the calculation process are presented.
Heat exchangers are used for transferring thermal energy between fluids, surfaces, or combinations of these, when they are at differing temperatures and in thermal contact. Typical applications include heat recovery, pasteurization, distillation, and heating or cooling of a particular fluid stream. The fluids can be separated by a wall or in direct contact. If there is a wall acting as the heat transfer surface separating fluids, appendages, or fins, can be connected to it in order to increase the heat transfer surface area.
Classification
Heat exchangers can be classified according to construction type, flow arrangement, or surface compactness, among other types possible. If the classification is by construction, the types of heat exchangers are: plate, tubular, extended surface, or regenerative. If classification is by flow arrangement, the types can include single-pass or multi-pass of counter-flow, parallel-flow, cross-flow, or combinations of flow.
Fin-and-Tube Heat Exchangers
The fin-and-tube heat exchanger (HE) studied in this project is classified as extended surface, single-pass with cross-flow (simplifying the header design at the entrance and exit). This type of heat exchanger is widely used in various thermal engineering applications, including chemical plants, food industries, HVAC, automotive, aircraft, and more. They consist of a block of parallel continuous fins with round tubes mechanically or hydraulically expanded into the fins, a popular heat exchanger designed for fluid to flow in the tubes and gas between the fins (see Figure 3).
The advantages to using more compact heat exchangers such as the fin-and-tube are many. The extended surfaces (fins) are designed to increase the heat transfer area per unit volume, resulting in compact units of reduced space and weight (up to 10 times greater surface area per unit volume when compared to shell-and-tube exchangers), with higher heat transfer coefficients than other less compact heat exchanger types. There is also flexibility when designing the surface area distribution between the hot and cold sides. Substantial cost savings are expected. For sensitive materials, tighter temperature control is an advantage, improving product quality. Multiple fluid streams can be accommodated.
There are also limitations to using fin-and-tube heat exchangers. Normally one side must be a gas or liquid with a low coefficient of convection, h. They are difficult to mechanically clean, requiring non-corrosive clean fluids. Temperature and pressure limits are lower than some other types due to brazing or mechanical expansion when joining the fins to the tubes (though pressure can be high on the tube side).
[Rohsenow et. al, 1998]
Typical Materials and Geometry: Fin-and-tube HE
Fins are commonly made of aluminium, while tubes are made of copper. Typical geometry: fin thickness 0.11 to 0.13 mm, tubes with outside diameter 10 mm, transverse pitch 25 mm, longitudinal pitch 22 mm, and fin density of 6 to 16 fins per inch. (However, smaller tube spacing and tube diameters are becoming more widespread). The geometry of the heat exchanger in this project is close to the typical geometrical ranges listed above, with the details found in Section 2.2: Computational Domain and Boundary Conditions.
[Baggio and Fornasieri, 1994]
1 2.1 Governing equations and numerical schemes.
The governing equations for this project are the three-dimensional continuity, Navier-Stokes for momentum, energy, and scalar transport equations for steady-state flow, and can be written (generally) as follows:
Continuity equation:
Equation 1 [pic]
Momentum equation:
Equation 2 [pic]
[pic]
Energy equation:
Equation 3 [pic]
General transport equation (for scalars):
Equation 4 [pic]
The general equations 1-3 are used in the CFD computations to calculate the flow field for both thermal and fluid (air) dynamics, solving for heat transfer and pressure drop. They are discretised and solved by the finite volume method using OpenFOAM, an open-source CFD code. It is solved on a staggered grid using solvers for laminar and turbulent flow, with the latter solution solved using the Reynolds Averaged Navier-Stokes equations (RANS) with both k-epsilon and SST k-omega turbulence models. To ensure coupling between velocity and pressure, the SIMPLE algorithm is used. More regarding the details of the governing equations, numerical schemes, turbulence modelling, and solution algorithms used for this project are reviewed in Section 3: Computational Fluid Dynamics, with a summary of computation procedures in Appendix A1: CFD Computations.
2 2.2 Computational Domain and Boundary Conditions
The pre-processing open-source software Salomé is used to create and mesh the computational model. A diagram of the studied model is shown in Figure 4, and consists of the air flow area between two fins of plain fin geometry and around the surfaces of two rows of tubes, and a schematic of the model with dimensions is shown in Figure 5, with the geometrical values listed in Table 1.
Table 1. Geometric dimensions of heat exchanger model
|Geometric Parameter | |
|Fin thickness t |0.130 mm |
|Fin pitch Fp |2.240 mm |
|Fin collar outside diameter Dc |10.23 mm |
|Transverse pitch Pt |25.40 mm |
|Longitudinal Pitch Pl |22.00 mm |
|Tube wall thickness δ |0.336 mm |
|Number of tube rows |2 |
[pic]
The computational domain is actually 8 times the original heat transfer area (as illustrated in Figure 5), and is defined by 0 < x < 8Pl, 0 < y < Pt/2, and 0 < z < Fp., while the actual modelled heat exchanger length is equal to twice the longitudinal pitch Pl. The volume representing the air which passes through the gap between the two fins is extended upstream from the inlet and downstream from the outlet in order to reduce oscillations and ensure a representative flow in the computational domain of the actual heat exchanger. The heat exchanger model with its extended volume is illustrated in Figure 5, while the actual area of interest for the heat exchanger simulation is shown later in Figure 9 (and the middle section of Figure 5).
Figure 5. Computational domain, including boundary conditions (BC) and extended flow volumes.
The computational domain has contains boundary conditions as shown in Figure 5 with the following conditions:
• Tube surfaces, Dirichlet BC:
T = Tw,
Air velocity: u = v = w = 0.
• Fins, Dirichlet BC:
T = Tfw
Air velocity: u = v = w = 0
• Inlet, Dirichlet BC:
Uniform velocity u = uin,
v = w = 0
T = 5 ºC.
• Outlet, Neuman BC:
Zero gradients, u, v, w, pressure, and temperature. (One-way),
• Free stream planes: (top and bottom planes of the extended surface areas):
‘slip’ conditions?: (∂u/∂z)=0, (∂v/∂z) = 0, w = 0, (∂T/∂z) = 0.
• Side planes: symmetry planes
(∂u/∂y)=0, v = 0, (∂w/∂y) = 0, (∂T/∂y) = 0
The entire computational domain was made up of 50,375 finite volumes, with a structured grid throughout most of the domain, while the areas around the tubes are more unstructured. The cell number was chosen based on the results of a grid independence test (more information about the mesh and grid independence test is found in the next section, Section 2.3: Numerical Mesh, while results of the grid independence test are presented in Section 4.1: Results of Grid Independence Test).
3 2.3 Numerical mesh
The geometry and computational domain were created using the open-source geometry creation and meshing software Salomé. The OpenFOAM version 1.4.1 used for this project did not have sophisticated geometry and meshing capabilities (though the very latest version is improved with the ‘snappy hexmesh’), and therefore another meshing program was required. For this purpose, Salomé was chosen since it is also open-source. Both tetrahedral and hexagonal meshes were created and used for a grid independence study for determining the ideal number of cells to use in the mesh. This section describes the general process of creating the geometry and meshes made for this project.
Tetrahedral meshes First the heat exchanger geometry was created in the geometry module of Salomé. Then the mesh module was used to create either tetrahedral or hexahedral cells. Figure 6 illustrates a mesh with 8465 tetrahedral cells. The program gives choices for various cell properties such as average length, but it was not possible to refine the mesh in a particular manner without making several partitions (which was difficult in Salomé); the program creates the mesh based on vague input parameters. As can be seen in Figure 6, the tetrahedral mesh is unstructured, [pic]
Figure 6. A mesh made up of 8465 tetrahedral cells.
with finer mesh around the curved areas of the tubes (Salomé automatically created a finer mesh in certain areas of the mesh). The extended volume has larger cells, thereby saving computations for the areas of the most change. It was not possible to make more than three layers of cells through the thickness of the geometry (z-direction) or refine the mesh in a specific manner.
An example of one of the problems encountered while making the mesh is illustrated in Figure 7, which is the other side of the mesh in Figure 6 (Fig. 6 is ‘flipped over’). It can be seen that instead of making a complete 3-dimensional mesh, some of it was 2-D, while parts of it were 3-D.
[pic]
Figure 7. Illustration of a meshing problem: the underside of the mesh shown in Figure 6.
Tetrahedral meshes were created with cell numbers ranging from approximately 8,000 to 150,000 cells to be used for the grid independence test (described later in this section). Once the mesh is created in Salomé, it could be converted into the OpenFOAM program using the ‘ideasUnvToFoam’ command. This command for mesh conversion, however, does not work properly in OpenFOAM 1.4.1 due to a bug, and a new C++ file containing the correction had be placed into the source code and compiled. The fix was found on a discussion group in the OpenFOAM website [OpenFOAM, 2008].
When the CFD simulations were first attempted, other problems with the calculations arose due to geometrical parameters. The problems were discovered after importing the mesh into OpenFOAM and trying to run the files. The cases could not make calculations, and it was necessary to investigate the problems. The ‘checkMesh’ command in the OpenFOAM console window is used to investigate information about the geometry and mesh (which should always be run after creating a mesh in OpenFOAM to verify the mesh is okay).
The first problem discovered was due to the scalar units of the mesh created in Salomé. The geometry and mesh were originally created using meters as the unit of length (i.e. fin pitch 2.24 mm was created as 0.00224 m). Although this is the correct scale, there was a problem with the computations in OpenFOAM. To solve this problem, the geometry was re-created and the mesh calculated in Salomé using millimetres as the length unit, and then upon conversion into OpenFOAM, the ‘transformPoints’ command is used to scale the geometry into the correct dimensions in meters.
Other problems arose with the computations due to the non-orthogonality and skewness of the cells in the mesh. Orthogonal meshes contain perpendicular gridlines at intersections, unlike non-orthogonal meshes where grid lines do not intersect one another at 90º angles. Another problem related to the geometry was the skewness between cells. If the non-orthogonality is below 30º or skewness is less than 10, specific inputs into the program are used, which in this case they are necessary, as the ‘checkMesh’ command showed the average non-orthogonality to be 2.4º, and maximum skewness to be 0.74. In the ‘fvSolution’ file, the non-orthogonality correction loops were set to ‘zero’, and the ‘snGradScheme’ (surface normal gradient) in the fvSchemes file set to ‘uncorrected’.
[pic]
[Hjertager, 2008] [Versteeg and Malalasekera, 2007]
Hexagonal Meshes After creating the tetrahedral meshes, it was determined that a more structured grid would be better to use, since that is the type most commonly used and known to provide accurate calculations. It was also determined important to have the capability of making different boundary conditions for the extended volumes, as compared with the actual heat exchanger core to be simulated. In this way, the free-flow air coming into the inlet (and outlet) of the core has the correct boundary conditions (i.e. temperature and x-direction velocity) and are not forced to have the same conditions as the heat exchanger, which would cause the simulations to be less accurate.
The Salomé version used for creating the tetrahedral mesh described previously could not make the more structured orthogonal mesh desired, and it was necessary to download and install the newest version: Salomé 3.2.9. Then, to be capable of creating different boundary conditions for the ‘top’ and ‘bottom’ of the model, it was necessary to divide the geometry up into partitions in Salomé. To do this, a new geometry was created again, and a new mesh with partitions built into the mesh construction. This mesh can be seen in Figure 8, where the partitions are visible. As with the tetrahedral meshes, there was a range of meshes created for use in the grid independence study; hexagonal mesh cell numbers ranged from approximately 3,000 to 117,000 cells. After the grid independence study, the mesh of approximately 50,000 cells was chosen, and is shown in Figure 9.
[pic]
It can be seen from Figure 9, that this mesh is much more structured than the tetrahedral cell mesh. However, there are still some unstructured areas around the tube walls. It was attempted to make a finer mesh in the central heat exchanger core, and coarse mesh in the extended volume areas, but the process was not successful; the calculations could not connect the sections of the geometry through the partitions. Therefore, the final mesh is as shown in Figure 9.
Grid Independence Testing
In all, there were 11 different grid systems investigated to determine how fine the grid must be and to validate the solution independency of the grid. The tetrahedral meshes contained the following number of fluid elements (approximately): 8,000, 25,000, 50,000, 75,000, 100,000, and 150,000, while the hexagonal meshes contained approximately 3,000, 30,000, 50,000, 67,000, and 117,000 cells. Simulations were run on all cases using simpleFoam (more information on this in section 3: CFD and OpenFOAM). The pressure drop between inlet and outlet was found for each simulation, and the results compared to determine when the grid is considered independent. The results of this test are shown and discussed in Section 4.1: Results of Grid Independence Test.
Summary of Meshing Process
Both tetrahedral and hexagonal meshes have been created for this project. The open-source software Salomé was utilized for creating both the geometry and the mesh, which could be converted to an OpenFOAM using the ‘ideasUnvToFoam’ command. Cell numbers for the meshes range from approximately 3,000 to 150,000 cells. A grid independence test was carried out to determine the optimal number of cells to use for the CFD simulations.
2.4 Performance parameters
This section describes how heat transfer and pressure drop are characterized. Included are dimensionless groups, equations for heat transfer and efficiency calculations, and equations for making pressure drop calculations. Following the descriptions of the performance parameters, the values as read from the graphs in the research done by Wang et al. (1996) for friction factor f and Colburn j-factor vs. Reynolds number are tabulated and graphed.
1 2.4.1 Dimensionless Groups
Accurate characterisation of the flow friction and heat transfer is very important in rating and sizing heat exchangers. Dimensional groups are used for this characterisation: heat transfer defined with the Colburn factor j and pressure drop defined by friction factor f. Below is a summary of the dimensionless groups used in this project, the equation to calculate it, and a brief description.
Reynolds number Re
The Reynold’s number Re represents the ratio of flow inertial forces to viscous forces. The Reynold’s number characteristic dimension for this study is the tube collar diameter Dc.
Equation 5 [pic]
where V is the minimum free-flow air velocity (in the minimum flow cross-section of the tube row), and is calculated:
Equation 6 [pic]
Fanning friction factor f
The Fanning friction factor is the ratio of wall shear stress to the flow kinetic energy. It is related to pressure drop in tube-and-fin heat exchangers as:
Colburn j-factor
The Colburn j-factor is the ratio of convection heat transfer (per unit duct surface area) to the amount virtually transferable (per unit of cross-sectional flow area):
Equation 8 [pic]
Nusselt number Nu
The Nusselt number is the ratio of convective conductance h to molecular thermal conductance k/Dh.
Equation 9 [pic]
The Nusselt number is based on the hydraulic diameter Dh. There are different calculations for this available in the literature. The hydraulic diameter in this study is the ratio of the 4 times the minimum free flow air-side area to the wetted perimeter (ratio of air-side surface area to heat exchanger length), and is given by the following expression [Fornasieri and Mattarolo, 1991]:
Equation 10 [pic]
Prandtl number Pr
The Prandtl number Pr is the ratio of a fluid’s momentum diffusivity to thermal diffusivity.
Equation 11 [pic]
[Rohsenow et. al, 1998] [Wang et al., 2006]
2 2.4.2 Heat Transfer and Efficiency
The heat transfer parameter Colburn j-factor can be determined from experimental data and is important in determining heat exchanger performance. To determine the value of j from the simulations, a series of equations related to heat transfer, efficiency, and the Nusselt and Reynolds numbers are worked through. These equations for determining the Colburn j-factor are provided in this section.
The overall heat transfer rate [pic] can be found in terms of the heat capacity rates ([pic]) and temperature differences between the inflow and outflow temperatures of the hot or cold side:
Equation 12 [pic] and
Equation 13 [pic]
where [pic] refers to the true mean temperature difference.
Efficiency determinations are used to find the heat transfer coefficient U for design and analysis of heat exchangers. The three methods most used for determining this are: ε-NTU, P-NTU, and LMTD. They are briefly described below. The ε-NTU method is described in more detail, since this is the method used for calculating efficiency in this project.
ε-NTU Method: The efficiency factor ε is a function of NTU, C* (ratio of minimum heat capacity rate to maximum heat capacity rate) , and flow arrangement.
P-NTU Method: The efficiency factor P is a function of NTU, R (ratio of temperature difference on one side to temperature difference on the other), and flow arrangement. This is similar to the ε-NTU method, but uses R instead of C*.
LMTD Method: The efficiency factor F is a function of P, R, and flow arrangement, and is a ratio of actual mean temperature difference to the log-mean temperature difference (LMTD).
The ε-NTU method is the one used in the research article for validation. The calculations were based on unmixed-unmixed cross-flow configuration, and expressed as [Kays and London, 1998] [Wang et al., 2006]:
Equation 14 [pic]
where
Equation 15 [pic]
Equation 16 [pic]
Equation 17 [pic]
NTU, or “Number of Transfer Units”, is the ratio of UA (overall conductance) to min. heat capacity rate Cmin. NTU determines the ‘thermal size’ of a heat exchanger. The overall heat transfer resistance is determined with the following equation:
Equation 18 [pic]
In the case of fouling, Equation 15 will also include resistance terms relevant to fouling. There may also be some thermal resistance due to the bond between the fin and tube – this resistance would be added to the lower right side of Equation 18.
To find the water-side heat transfer coefficient hi, Equation 19 for fully developed turbulent forced convection through a duct, the Gnielinski correlation is used [Gnielinski, 1976] [Wang et al., 2006]:
Equation 19 [pic]
where
Equation 20 [pic]
The fin surface effectiveness[pic] is the ratio of actual heat transfer to the maximum possible heat transfer which could occur if the base and fin are both at the same temperature. This is described in the following equation set, where A0 = Afin + Abase and the actual fin efficiency[pic] is calculated using the Schmidt approximation for staggered plate-fin geometry [Schmidt, 1949] [Wang et al., 2006]:
Equation 21 [pic]
Equation 22 [pic]
where
Equation 23 [pic]
Equation 24 [pic]
Equation 25 [pic]
When the equations are solved and the air-side heat-transfer coefficient ho has been calculated, the Colburn j-factor can be calculated by first solving for the Nusselt and Reynolds numbers, then using the calculated values in the final calculation for the j-factor. As shown previously in equations 5, 8, and 9, the Nusselt number, Reynolds number, and j-factor are calculated:
[pic] [pic] [pic]
[Rohsenow et. al, 1998] [Wang et al., 2006]
Experimental values for Colburn j-factor
To validate the simulations in this project, the equations previously mentioned (12-25) are used to calculate Colburn j-factors from simulation results and compared with the experimentally determined results of Wang et al. (1996) for the particular heat exchanger geometry are used. The values of the Reynolds number and j-factor as read from the graph in the article for the specific heat exchanger geometry studied in this project are listed in Table 2 and presented in the graph in Figure 10.
Table 2. Colburn j-factor vs. Re,
experimentally determined.
|j |Re |
|0,0420 |330 |
|0,0270 |600 |
|0,0230 |790 |
|0,0170 |1300 |
|0,0140 |1700 |
|0,0120 |2900 |
|0,0094 |4300 |
|0,0090 |5200 |
|0,0084 |6200 |
|0,0081 |7200 |
[Wang et al., 1996]
It can be seen in Table 2 and Figure 10 that the Colburn j-factor decreases with increasing Reynolds number. It ranges from 0.042 at the lowest Reynolds number to 0.0081 at the high Reynolds number. A plot of this type (j vs. Re) for a typical circular tube normally have a more distinct dip in the transition region than in this graph for air flow through a heat exchanger. However, a slight change is seen at Reynolds number 1300, and the graph appears to level off again at Reynolds number 2900. It can therefore be determined from this graph that the laminar flow region goes up to around Reynolds number 1300 (it is difficult to determine exactly without more data points), with a transition region after that, and the turbulent flow begins at a point around Reynolds number 2900. Values for the Colburn j-factor determined from simulations, and followed by the appropriate calculations for this project are compared to these experimental values. During the comparisons, it is kept in mind that the uncertainties in the Colburn j-factor values (experimentally determined values) can be high for the lower Reynolds numbers (uncertainty is ± 9.4 % at Reynolds number 600, and may be even higher at Reynolds number 330, the value of which was not provided in the article.)
[Rohsenow et. al, 1998][Wang et al., 2006]
[pic]
Figure 10. Plot of Colburn j-factor vs. Reynolds number, determined experimentally [Wang et al., 1996]
3 2.4.3 Pressure Drop
The pressure drop determines the amount of pumping power needed to run a heat exchanger. It is therefore important to characterize the pressure drop for design. This section describes how the pressure drop relates to the pumping power, followed by a description of what causes the pressure drop and finally the pressure drop equations for tube-and-fin heat exchangers are presented.
Pumping power P is often seen as an important design constraint because the pressure drop in a heat exchanger (along with associated pressure drops in the inlet/outlet headers, nozzles, ducts, etc.) is proportional to the amount of fluid pumping power needed for the heat exchanger to function, as given by the following expression:
Equation 26 Pumping power [pic]
The overall pressure drop consists of two parts: (1) the pressure drop in the heat exchanger core, and (2) the pressure drop from associated devices the fluid flows through before and after the heat exchanger core (i.e. inlet/outlet manifolds, nozzles, valves, fittings, ducts, etc.).
The core pressure drop is due to the following:
1. Friction from the fluid flowing across the heat transfer surface (i.e. skin fraction, form drag, internal contractions and expansions).
2. Momentum effect (fluid density changes causing a pressure drop).
3. Sudden contraction or expansion at inlet and/or outlet.
4. Gravity effects (if there is a change in elevation between the inlet and outlet of the exchanger – normally negligible with gases) causing a pressure drop (static head) from the change in elevation. This pressure drop is given by the following expression: Δp= ± (ρmgL/gc), with the “+” used in the case of vertical up-flow, while the “-” is used for vertical down-flow.
The actual calculation for pressure drop depends on the specific type of heat exchanger being studied. For fin-and-tube heat exchangers, the pressure drop equation is given in Equation 27 [Wang et al., 2006]:
Equation 27
However, the entrance and exit loss effects Kc and Ke become zero when flow is normal to the tube banks or through wire matrix surfaces, resulting in the following equation:
Equation 28 [pic]
The definitions for the specific variables are expanded from the Nomenclature section:
G: u*ρ. G is the mass velocity entering the core based on minimum free-flow area.
gc: A gravitational constant (equals 1 when working with SI units).
[pic]: This is the specific volume (1/ρ) at inlet temperature.
[pic]: This is the specific volume (1/ρ) at outlet temperature.
[pic]: This is the average specific volume ([pic]+[pic])/2.
σ : Sigma represents the ratio of minimum free-flow area to frontal area.
[pic]: Flow cross-sectional area.
In this project, Equation 28 is used to calculate the friction factor for the CFD simulations, since the pressure measurements are taken at the inflow and outflow of the computational domain, and the entrance and exit loss effects would occur prior to and after the inflow and outflow, respectively.
Experimental values for Fanning friction factor f
From the pressure drop data determined from the simulations, friction factors can be calculated using the pressure drop equation (solving for friction factor f). To validate the pressure-loss simulations in this project, the Fanning friction factor f determined experimentally from Wang et al. (1996) are used. The values for Reynolds number and friction factor f as read from the graph in the article for the specific heat exchanger geometry studied in this project are listed in Table 3 and presented in the graph in Figure 11.
Table 3. Friction factor f vs. Re,
experimentally determined.
|f |Re |
|0,1100 |330 |
|0,0730 |600 |
|0,0630 |790 |
|0,0460 |1300 |
|0,0420 |1700 |
|0,0330 |2900 |
|0,0270 |4300 |
|0,0240 |5200 |
|0,0220 |6200 |
|0,0210 |7200 |
[Wang et al., 1996]
It can be seen in Table 3 and Figure 11 that like the Colburn j-factor, the Fanning friction factor f also decreases as the Reynolds number increases. It ranges from 0.11 at the lowest Reynolds number to 0,021 at the high Reynolds number. A slight change can be seen at Reynolds number 1300, and the graph levels off again between Reynolds number 1300 and 1700. It can therefore be determined from this graph that the laminar flow region goes up around Reynolds number 1300, with a transition region thereafter, and from the slight change in graph again, the turbulent flow regime seems to begin at around Reynolds number 1700 (or possibly it is closer to 2900 since that is what was found on the j-factor graph). However, more data points would be necessary to accurately determine the critical Reynolds values for the different flow regimes. Values for the friction factor f as calculated from the pressure drop values of the simulations are compared to the experimental values from Wang et al. (2006). During the comparisons, it is kept in mind that the uncertainties in the Fanning friction factor f values (experimentally determined values) can be high for the lower Reynolds numbers (the uncertainty is given as ± 17.7 % at Reynolds number 600, and may be even higher at Reynolds number 330, the value of which was not provided in the article.)
[Rohsenow et. al, 1998][Wang et al., 2006]
[pic]
Figure 11. Plot of friction factor f vs. Reynolds number, determined experimentally [Wang et al., 1996]
3 Computational Fluid Dynamics
The purpose of this project is to use open-source CFD software to simulate pressure loss and heat transfer in a heat exchanger and validate the simulation with an actual experimental results from the literature. Different solvers and turbulence models are used to try to determine the most accurate CFD method for predicting pressure loss and heat transfer in this type of compact fin-and-tube heat exchanger.
Computational fluid dynamics (CFD) is a computer-based simulation method for analysing fluid flow, heat transfer, and related phenomena such as chemical reactions. This project uses CFD for analysis of flow and heat transfer (not for analysis of chemical reactions). Some examples of application areas are: aerodynamic lift and drag (i.e. airplanes or windmill wings), power plant combustion, chemical processes, heating/ventilation, and even biomedical engineering (simulating blood flow through arteries and veins). CFD analyses carried out in the various industries are used in R&D and manufacture of aircraft, combustion engines, as well as many other industrial products.
It can be advantageous to use CFD over traditional experimental-based analyses, since experiments have a cost directly proportional to the number of configurations desired for testing, unlike with CFD, where large amounts of results can be produced at practically no added expense. In this way, parametric studies to optimise equipment are very inexpensive with CFD when compared to experiments.
This section briefly describes the general concepts and theory related to using CFD to analyse fluid flow and heat transfer, as relevant to this project. It begins with a review of the tools needed for carrying out the CFD analyses and the processes required, followed by a summary of the governing equations and turbulence models (with details given in Appendix A1: CFD Computations) and finally a discussion of the discretisation schemes (details in Appendix A1: CFD Computations) and solution algorithms is presented.
[Versteeg and Malalasekera, 2007]
1 3.1 CFD Computational Tools
This section describes the CFD tools required for carrying out a simulation and the process one follows in order to solve a problem using CFD. The hardware required and the three main elements of processing CFD simulations: the pre-processor, processor, and post-processor are described.
There is a variety of commercial CFD software available such as Fluent, Ansys CFX, ACE, as well as a wide range of suitable hardware and associated costs, depending on the complexity of the mesh and size of the calculations. Commercial CFD packages can cost up to about $20000 (US Dollars) per year for licenses, maintenance, and support. Complicated transient cases with fine meshes will require more powerful computer processors and RAM than simpler cases with rough meshes. A typical engineering workstation (i.e. 32 GB processing RAM with quad processors) at a cost of approximately $3000-$5000 (US Dollars), or a combination of several processors running in parallel, is probably the minimum investment needed to get started.
The work for this project was carried out on a HP Pavilion laptop with dual processors totalling 2 GHz RAM, running on Linux Operating System downloaded free from Caelinux. The download from Caelinux included open-source software Salomé for geometry construction and meshing, OpenFOAM for the CFD calculations, paraView for visualization of results, along with other useful scientific and mathematics related software. Calculations for this project were carried out for approximately 50,000 cells (CFD calculations are often made for 1-2 million cells – or more). On my system, the steady-state solvers took between 1-3 hours to finish calculations, while the transient simulation took 2-3 days running in parallel on both processors.
One of the purposes of this project is to use all open-source CFD software instead of commercial software for the simulations. This type of software is advantageous for smaller companies to use, as the cost of commercial CFD package licenses can be prohibitive.
To run a simulation, three main elements are needed:
1. Pre-processor: A pre-processor is used to define the geometry for the computational domain of interest and generate the mesh of control volumes (for calculations). Generally, the finer the mesh in the areas of large changes, the more accurate the solution. Fineness of the grid also determines the computer hardware and calculation time needed. The open-source pre-processor used for this project is called Salomé.
2. Solver: The solver makes the calculations using a numerical solution technique, which can use finite difference, finite element, or spectral methods. Most CFD codes use finite volumes, which is a special finite difference method. First the fluid flow equations are integrated over the control volumes (resulting in the exact conservation of relevant properties for each finite volume), then these integral equations are discretised (producing algebraic equations through converting of the integral fluid flow equations), and finally an iterative method is used to solve the algebraic equations. (The finite volume method and discretisation techniques are described more in the next sections. OpenFOAM CFD code is used for solving the simulations in this project.
3. Post-Processor: The post-processor provides for visualisation of the results, and includes the capability to display the geometry/mesh, create vector, contour, and 2D and 3D surface plots. Particles can be tracked throughout a simulation, and the model can be manipulated (i.e. changed by scaling, rotating, etc.), and all in full colour animated graphics. ParaView is the open-source post-processor used for this project.
Problem-Solving with CFD
There are many decisions to be made before setting up the problem in the CFD code. Some of the decisions to be made can include: whether the problem should be 2D or 3D, which type of boundary conditions to use, whether or not to calculate pressure/temperature variations based on the air flow density, which turbulence model to use, etc. The assumptions made should be reduced to a level as simple as possible, yet still retaining the most important features of the problem to be solved in order to reach an accurate solution.
After the above decisions are made, the geometry and mesh can be created. The grid should be made as fine as required to make the simulation ‘grid independent’. To determine the fineness required, a grid dependence study is normally carried out by making a series of refinements on an initially course grid, and carrying out simulations on each to determine when the key results of interest do not change, at which point the grid is considered independent. In this project, a grid with approximately 50,000 cells was chosen after carrying out such a study (described in the next section). To reach a converged solution, relaxation factors and acceleration devices can be chosen. In this project, relaxation factors for all the parameters to be solved and the GAMG smooth-solver for pressure were used to assist in convergence and speed optimisation.
Finally, to ensure accuracy of the simulations, they should be validated against experimental data. This project’s simulation results are compared to an experimental study reported in the literature [Wang et al., 1996]. Details of the heat exchanger geometry, initial boundary conditions related to flow and temperature were followed as closely as possible when building this CFD model.
[Versteeg and Malalasekera, 2007]
OpenFOAM Solvers
Three solvers were used for simulating the ten cases representing each of the Reynolds numbers. The simpleFoam solver solves only for flow fields, while rhoSimpleFoam solves for flow and temperature under steady-state conditions, and rhoTurbFoam also solves flow and temperature, but for transient cases.
2 3.2 CFD Governing Equations
This section is a summary of the governing equations used in CFD to mathematically solve for fluid flow and heat transfer, based on the principles of conservation of mass, momentum, and energy. Details of how they are actually used in the CFD computations are described in Appendix A1: CFD Computations.
Conservation Equations
The conservation laws of physics form the basis for fluid flow governing equations (previously listed as Equations 1-3 in Section 2.1: Governing Equations and Numerical Schemes). The laws are:
• Law of Conservation of Mass: Fluid mass is always conserved. (Equation 1)
[pic]
• Newton’s 2nd Law: The sum of the forces on a fluid particle is equal to the rate of change of momentum. (Equation 2)
[pic]
• First Law of Thermodynamics: The rate of head added to a system plus the rate of work done on a fluid particle equals the total rate of change in energy. (Equation 3)
[pic]
The fluid behaviour can be characterised in terms of the fluid properties velocity vector u (with components u, v, and w in the x, y, and z directions), pressure p, density ρ, viscosity μ, heat conductivity k, and temperature T. The changes in these fluid properties can occur over space and time. Using CFD, these changes are calculated for small elements of the fluid, following the conservation laws of physics listed above. The changes are due to fluid flowing across the boundaries of the fluid element and can also be due to sources within the element producing changes in fluid properties. This is called the Euler method (tracking changes in a stationary mass while particles travel through it) in contrast with the Lagrangian method (which follows the movement of a single particle as it flows through a series of elements).
A review of the governing equations (1)-(3) and how they are used in CFD modelling is given in Appendix A1: CFD equations.
[Hjertager, 2007][Versteeg and Malalasekera, 2007]
3 3.3 Turbulence Modelling
The Reynolds number studied for the heat exchanger flow in this project range from approximately 330 to 7000, which means that the flows can be laminar, transitional, or turbulent at the different flow velocities. Because these flow regimes behave differently, it can be necessary to model the flow in different ways. In this project, two turbulence models (k -epsilon and Menter SST k-omega) are utilized in order to investigate which is best to use for the different types of flow. This section describes some of the commonly used turbulence models, and then presents relevant details regarding the specific turbulence models used in this project. A laminar flow model is also used for comparison with the turbulence models.
Table 4 lists commonly used turbulence models and a short description of the model. Following the table is a brief discussion of the other methods which can be used to predict turbulence (besides using RANS turbulence models). After this, a more detailed description of the turbulence models used in this project is presented.
As can be seen in Table 4, there are a large number of different turbulence models. The turbulence models listed are all based on variances of the Reynolds averaged Navier-Stokes (RANS) equations, in which average values for turbulent fluctuations are used for modelling the turbulence (described in Appendix A1).
There are turbulence models for compressible and incompressible flows, and most of the listed models have been built into OpenFOAM source code. Several variations of the k-epsilon model have been made, as well as low-Reynolds numbers modifications of it. The high Reynold’s number models listed use log-law type wall functions (i.e. Menter SST k-omega). The low Reynolds number models calculate flow to the wall, and with these models, it is important for the y+ value to be approximately 1, whereas with high-Reynolds models, y+ should range from approximately 30-60 to 300-400 in the log layer.
Table 4. Commonly used turbulence models. [Hjertager, 2005 & 2008] [Versteeg and Malalasekera, 2007]
|Turbulence Model |Flow* |Model Type** |Description | |
|Spalart-Allmaras |C, I |One-equation |Mixing-length model for external flow | |
|k- ε |C, I |Two-equation |Standard k- εpsilon model | |
|Two-layer k-ε | |Two-equation |Viscous and turbulent flow regions k-ε model | |
|RNG k-ε |C, I |Two-equation |Renormalisation group k-εpsilon model | |
|Wilcox k-ω |n/a |Two-equation |Turbulence frequency ω = ε /k as 2nd variable | |
|Menter SST k- ω |C, I |Two-equation |High Re shear stress transport k-omega model | |
|Realizable k-ε |C, I |Non-linear |Eddy viscosity model for high Re flows | |
|LRR R- ε |C, I |RSM |R- ε Reynold’s stress (RSM) model | |
|Launder-Gibson R- ε |C, I |RSM |RSM R- εpsilon model + wall reflection terms | |
|Algebraic stress | |Two-equation |Anisotropy of stresses (RSM) in simpler form | |
|Lam Bremhorst k-ε |I | |Low Re k-εpsilon model | |
|Q-z |I | |Gibson Dafa’Alla’s Q-zeta model for low Re | |
|Lien k-ε |I |Cubic |Cubic non-linear k-εpsilon model | |
|Lien low Re k-ε |I |Cubic |Cubic non-linear low Re k-εpsilon model | |
|Lien Leschziner k-ε |I | |Low Re k-εpsilon model | |
|Launder Sharma k-ε |C | |Low Re k-εpsilon model, combusting flows | |
|Shih k-ε |I |Quadratic |Shih’s quadratic non-linear k-εpsilon model | |
*Flow code availability in OpenFOAM: C=compressible, I=incompressible, n/a = not available
**Number of extra transport equations to be solved along with the RANS equations.
There are also other methods that can be used to predict the turbulence (besides the RANS models listed in Table 4). Empirical correlations (i.e. using friction factor, Nusselt number, or Sherwood
number) can predict simple turbulent flow without requiring CFD calculations or large computers. LES (large eddy simulation) solves for motion of the larger eddies while DNS (direct numerical simulation, the most accurate and time-consuming/computationally expensive method of all for calculating turbulence) solves for all motion using the Navier-Stokes equations. LES and
DNS are computationally very expensive due to requirements including very fine 3D mesh, transient transient behaviour. Due to hardware requirements, LES and DNS were not used for this project.
In this project, the k-epsilon and Menter k-omega SST models, both of which utilize RANS equations, are used as turbulence models. The RANS equations are described some detail in Appendix A1: CFD Computations. First, a short explanation for why other turbulence models were ruled out for this project is given, and then the basic turbulence concepts, k-epsilon and k-omega SST models are explained.
To decide which turbulence models to use for this project, considerations were taken regarding relevance to the particular type of flow involved (applicability), complexity of the physics and time for computing, whether the model is for compressible or incompressible flow and availability in OpenFOAM, and how well-known the model is (for accuracy). Keeping these considerations in mind, the 8 turbulence models available for compressible flow (compressible flow modelling is
necessary, since air is flowing through the heat exchanger) in OpenFOAM were considered:
• Spallart-Allmaras,
• standard k-epsilon,
• RNG k-epsilon,
• Menter SST k-omega,
• Realizable k- εpsilon,
• LRR R- εpsilon,
• Launder-Gibson R- εpsilon, and
• Launder-Sharma k- εpsilon.
The Spallart-Allmaras is a fairly new one-equation mixing length model which solves for the turbulent kinematic viscosity using a modelled transport equation. It specifies the eddy length scale with an algebraic formula, is designed for external aerodynamic boundary layer economical computations, but cannot accurately describe flows involving separation and recirculation. This project involves internal flow with both separation and recirculation, making this model irrelevant for this project.
The RNG (renormalization group) k-epsilon model uses ‘highly abtruse’ mathematics to extend the eddy viscosity turbulence model and changes the governing equations by removing smaller scales of motion, replacing them with large motions and modifying the viscosity term. Results have been good for backward-facing steps, while other results have been mixed. The RNG model is not very commonly used, has a high computational overhead, and therefore not used for this project.
The ‘realizable’ k-epsilon model is a non-linear version of the k-epsilon model. It retains the two-equation k-epsilon equations, but expands the model by including additional effects to account for Reynold’s stress anisotropy without actually using the seven extra equations used in the RSM (described in the next paragraph) to exactly model the Reynolds stresses. The turbulent viscosity expression and the rate of dissipation of kinetic energy equation of the standard k-epsilon model are both changed in the realizable k-epsilon model to take into account that turbulence does not always adjust itself instantaneously while moving through the flow domain, meaning that the Reynolds stress is partially dependent on the mean strain rate itself. This means that the non-linear realizable k-epsilon model allows for the phenomena of the state of turbulence lagging behind the changes disturbing the turbulence production and dissipation balance. This turbulence model was not used on the project, but like the RSM, could be considered for future work with this type of heat exchangers to account for anisotropic Reynold’s stresses in the flow.
The Reynold’s Stress Equation Models (RSM): LRR R-epsilon and Launder-Gibson R-epsilon models, calculates the anisotropic Reynolds stresses present in typical flows of complex strain fields or significant body forces. RSM is the most complex of the turbulence models, and was designed to address the problems of the k-e model (which cannot predict flows in long non-circular duct because of the isotropic modelling of the Reynold’s stresses). The RSM can therefore accurately account for the Reynold’s stress field directional effects. Because of the many Reynold’s stresses to model, there are seven extra partial differential equations to solve, making computing costs very high. The RSM is not as widely validated as the k-epsilon model and not used for this project due to the computing requirements, although it could be an interesting continuation to the project in the future, since the long thin channels between the heat exchanger fins probably means the Reynold’s stresses in the heat exchanger are anisotropic.
The Launder-Sharma k-epsilon model is designed for compressible and combusting low-Reynold’s number flow. Not much information regarding this model was found.
After considering the turbulence models available in OpenFOAM, it was decided to use the k-epsilon and Menter SST k-omega turbulence models. These include two equations for expressing turbulence in these models: (1) for the turbulent kinetic energy k (to express the turbulence velocity), and (2) for the rate of dissipation of the turbulent kinetic energy e (to express the turbulence length scale) in the k-epsilon model or of the specific dissipation rate ω in the k-omega model. The k-epsilon model is the simplest and most widely validated turbulence model, with only initial or boundary conditions required to be supplied, and has had good performance in the past with certain types of flows, although performance has not been the best for flows with curved boundary layers, swirling or anisotropic flows. The Menter SST k-omega model improves the k-ε model at the near-wall by using a k-w model at the near-wall region while retaining the k-epsilon model for the free-stream turbulent region far away from the wall. It is also well-known and used in industry. For this project, these two models: results from using the k-ε and k-ω models are compared, along with the laminar model. For future work with this type of flow, the ‘realizable’ k-epsilon and RSM models would be of interest to test in order to study the effect of anisotropic Reynold’s stresses on the flow and simulation results. The RANS equations are described in Appendix A1, while the turbulence models used in this project (based on the RANS equations): the k-epsilon and SST k-omega models, are described in the next sections.
[Hjertager, 2005] [Versteeg and Malalasekera, 2007]
1 3.3.1 k-εpsilon Turbulence Model
Turbulence models are required to predict the Reynolds stresses and scalar transport terms (described mathematically in Appendix A1) in order to close the system of expanded time-averaged RANS equations (continuity, momentum, and scalar transport equations), given in Equations 53, 54, and 56 in Appendix A1. Examples of established turbulence models are given in Table 4. In this section, the k-εpsilon turbulence model is described, starting with a description of the terms used for turbulent viscosity, turbulent velocity, and turbulent length scales, followed by a summary of the equations used in the actual turbulence model flow calculations, and finally a list of the required constants and their values is presented.
Turbulent eddy viscosity, velocity, and length scale
The k-εpsilon model presumes that the average flow is affected by the viscous and Reynolds stresses ([pic]) acting on it. In the time-averaged RANS momentum equations, both stresses are given on the right side to reflect their function as extra turbulent stresses on mean velocity components U, V, and W. As viscosity is the central concept for modelling the stresses, it should be noted that the turbulent viscosity [pic] is expressed as:
Equation 29a [pic] Equation 29b [pic]
From Equation 29a, it can be seen that the turbulent viscosity is a product of density and two new variables representing turbulent velocity [pic] and turbulent length scale[pic]. Turbulent velocity can be described as the typical velocity occurring in the largest eddies and can also be related to the turbulent kinetic energy according to Equation 29b. The turbulent length scale is the average length of the same eddies.
The new variables, [pic] and [pic] form the basis for the ‘two-equation’ k-ε turbulence model, meaning that in addition to the RANS equations, two more equations are required to solve for turbulent velocity and turbulent length using the model.
The length scale [pic] (for large eddies) is used in the k-ε model to define the length scale ε (for small eddies), for which a transport equation is used in the model, and represents the dissipation of the turbulent kinetic energy. The dissipation is expressed as:
Equation 30 [pic] or as related to viscosity [pic]
As seen in the second expression of Equation 30 for ε, dissipation of the turbulent kinetic energy is proportional to the rate of deformation of the eddies.
Reynolds stresses [pic] increase with deformation rate, and are related to viscosity, mean rate of deformation, and turbulent kinetic energy with Boussinesq’s proposal expressed as:
Equation 31 [pic]
Equation 31 can be used to calculate Reynold’s stresses in the final step of turbulence modelling. It is seen from this equation that the Reynold’s stresses are considered proportional to the dissipation rate reduced by the eddy turbulent kinetic energy. (The Kronecker delta [pic] (equal to 0 when i ≠ j, equal to 1 when i = j) ensures that the normal Reynolds stresses are each appropriately accounted for.) It can also be seen from Equation 57 (in the Appendix A1) for kinetic energy, that the kinetic energy allocates an equal third for each normal stress component (isotropic assumption). This is the reason for the inherent inaccuracy of the k-epsilon model, making it incapable of describing anisotropic flow.
Other scalar flow properties such as mass and heat can also be modelled using time-averaged values. Similar to the turbulent momentum transport’s proportionality to average velocity gradients, turbulent scalar transport is proportional to mean scalar value gradients and can be expressed as:
Equation 32 [pic]
where [pic] refers to turbulent (eddy) diffusivity.
As can be seen from Equation 32, the turbulent scalar property transport occurs with the same mechanism as in transport of momentum (mixing of eddies, represented by[pic]). For this reason, it can be assumed by ‘Reynold’s analogy’ that the value of [pic] is similar to[pic], the turbulent viscosity. The ratio of [pic] to [pic] is defined as the Prandtl/Schmidt number σt, and has a value which is normally constant with a value around unity.
The next section presents the k-epsilon model and the two extra k and ε transport equations (PDEs) for closing the system of time-averaged RANS equations. The model is based on the mechanisms causing changes to turbulent kinetic energy (i.e. turbulent viscosity and velocity fluctuations).
Turbulence Equations
The Reynolds Averaged Navier-Stokes equations are used for deducing the transport equation for the turbulent kinetic energy k and turbulent dissipation rate ε. The two new transport equations and dimensionless constant determination for the calculations are presented here. (RANS equations are described in Appendix A1).
Transport Equation for k
To obtain the equation for turbulent kinetic energy k, complicated algebra and rearrangements are made to the time-averaged continuity equation and the time-averaged Navier-Stokes equations for momentum. The mathematical manipulations are very extensive, and therefore only a short description of what is done, followed by the resulting equations is presented.
The 3 continuity equations are each multiplied by the respective fluctuating velocity component and then added together. The same process is carried out for the Reynolds equations for momentum. The two resulting equations are subtracted and rearranged extensively. Terms for viscous dissipation of turbulent kinetic energy (Equation 30) and the Boussinesq assumption (Equation 31) related to Reynold’s stress and turbulent viscosity are plugged in to the system. This results in the turbulent kinetic energy equation, as shown below in Equation 33. [Tennekes and Lumley, 1972]
Equation 33 [pic]
Terms (III) and (V) are replaced using scalar diffusion transport terms for (III) and the time-averaged term for (V) to result in the following:
[pic]
where [pic] is a constant turbulent Schmidt number for k
The terms (I)-(V) in Equation 33, the turbulent kinetic energy k transport equation, can be interpreted as the following:
|(I) |Transient term |Accumulation of k (rate of change of k) |
|(II) |Convective transport |Transport of k by convection |
|(III) |Diffusive transport |Difffusive transport of k by pressure, viscous stresses, and Reynolds |
| | |stresses (must be modelled) |
|(IV) |Production term |Rate of production of k due from the mean flow |
|(V) |Viscous dissipation |Rate of viscous dissipation of k (must be modelled) |
[pic] is a constant (turbulent Schmidt number) of the k equation
The above equation for transport of k, along with the equation for transport of ε, constitute the two additional transport equations to be solved in addition to the RANS equations in the k-epsilon turbulence model. The next section presents the equation for transport of ε.
Transport Equation for ε
Equation 34 [pic]
Similar to the transport equation for k, the transport equation for ε includes the terms I-V:
|(I) |Accumulation, or rate of change, of ε, |
|(II) |Rate of destruction of ε. |
|(III) |Diffusive transport of ε, |
|(IV) |Rate of production of ε, and |
|(V) |Transport of ε by convection, |
[pic],[pic], and[pic] are constants of the ε equation.
Turbulent Viscosity
Another parameter to find before calculating the Reynolds stresses is the turbulent viscosity[pic]. To determine this, Equations 29a, 29b, and 30 (as previously presented) are utilized. They are:
[pic] [pic] [pic], (or [pic] from a dimensional standpoint)
Using dimensional analysis of the above equations, eddy viscosity [pic] is found to be:
Equation 35 [pic] where [pic] is a constant.
Reynold’s Stresses
Finally, the Reynold’s stresses can be calculated using the Boussinesq assumption (Equation 31):
[pic]
The equations for the k-epsilon turbulence model have been presented in this section. They are the equations for turbulent kinetic energy k, dissipation of turbulent kinetic energy e, turbulent viscosity[pic], and Reynold’s stresses. Finally, the equations contain five constants, the values of which are listed below in Table 6.
Table 6. Constant values for k-epsilon turbulence model equations.
|[pic] |[pic] |[pic] |[pic] |[pic] |
|0.09 |1.44 |1.92 |1.0 |1.3 |
The next section explains SST k-omega turbulence model in contrast to the k-epsilon model, lists the equation and values of the constants.
[Hjertager, 2005] [Versteeg and Malalasekera, 2007]
2 3.3.2 SST k-omega Turbulence Model
The shear-stress transport (SST) k-omega turbulence model is a type of hybrid model, combining two models in order to better calculate flow in the near-wall region. It was designed in response to the problem of the k-epsilon model’s unsatisfactory near-wall performance for boundary layers with adverse pressure gradients. It utilizes a standard k-epsilon model to calculate flow properties in the free-stream (turbulent) flow region far from the wall, while using a modified k-ε model near the wall using the turbulence frequency ω as a second variable instead of turbulent kinetic energy dissipation term ε. Since the air in this project’s heat exchanger is flowing between two flat fins very close to each other, it is expected that the boundary layer flow has a strong influence on the results, and properly modelling this near-wall flow could be important for accuracy of the calculations. Therefore the SST k-omega turbulence model has also been chosen for CFD simulations in this project.
This SST k-omega model is similar to the k-epsilon turbulence model, but instead of ε as the second variable, it uses a turbulence frequency variable omega, which is expressed as ω = ε/k [s-1]. The SST k-omega model computes Reynolds stresses in the same way as in the k-epsilon model, using Equation 48.
The transport equation for turbulent kinetic energy k for the k-ω model is:
Equation 36 [pic]
[pic]
The terms (I)-(V) in Equation 36, the turbulent kinetic energy k transport equation for the SST k-omega turbulence model, can be interpreted as the following:
|(I) |Transient term |Accumulation of k (rate of change of k) |
|(II) |Convective transport |Transport of k by convection |
|(III) |Diffusive transport |Turbulent diffusion transport of k |
|(IV) |Production term |Rate of production of k |
|(V) |Dissipation |Rate of dissipation of k |
[pic] and [pic] are equation constants.
The transport equation for turbulent frequency ω for the k-ω model is:
Equation 37
[pic]
The general description for each of the terms in Equation 37 (I) to (V) are the usual terms for accumulation, convection, diffusion, production, and dissipation of ω. The last term (VI) is called a ‘cross-diffusion’ term, an additional source term, and has a role in the transition of the modelling from ε to ω.
The constants for the SST k-omega turbulence model are listed in Table 7.
Table 7. Constant values for SST k-omega turbulence model equations.
|[pic] |[pic] |[pic] |[pic] |[pic] |[pic] |
Additional modifications have been made to the model for performance optimisation. There are blending functions added to improve the numerical stability and make a smoother transition between the two models. There have also been limiting functions made to control the eddy viscosity in wake region and adverse pressure flows.
[Versteeg and Malalasekera, 2007]
Summary of Turbulence Models
There were many turbulence models to choose from to model the heat exchanger. After considering the availability of models in OpenFOAM, the applicability, whether it was computationally economical, and accuracy based on past validation, the k-epsilon and k-omega models were chosen to carry out simulations to comparison with the laminar flow model. The k-epsilon model is a well-known model and the most validated, and is also computationally economical, while the SST k-omega model is more accurate for boundary layer flow. Due to the thin channel and influence of the fin walls on the flow of air in the heat exchanger, it was decided to also carry out simulations using the SST k-omega model. The calculated values for turbulent intensity I, ´turbulent kinetic energy k,
turbulent dissipation ϵ, and omega ω, which are used for the simulations in this project are listed in Table 5.
| | |*I = 0.16 * Re-1/8 |k = 3/2(I*U)² |**ϵ = (C1⁰·⁷⁵*k¹·⁵) |ω = ϵ / k |
| | | |(Eqn. 47) |/(0.1*0.0254) | |
|U |Reynolds number (Re) |Intensity |kinetic energy |dissipation |frequency |
|[m/s] | |( I) |(k) |(ϵ) |(ω) |
|0,3 |330 |0,0775 |0,0008 |0,0015 |1,8424 |
|0,5 |600 |0,0719 |0,0019 |0,0055 |2,8495 |
|0,7 |790 |0,0695 |0,0035 |0,0137 |3,8545 |
|1,1 |1300 |0,0653 |0,0077 |0,0440 |5,6914 |
|1,5 |1700 |0,0631 |0,0135 |0,1010 |7,5051 |
|2,5 |2900 |0,0591 |0,0327 |0,3827 |11,7007 |
|3,7 |4300 |0,0562 |0,0649 |1,0702 |16,4850 |
|4,5 |5200 |0,0549 |0,0916 |1,7928 |19,5786 |
|5,4 |6200 |0,0537 |0,1262 |2,9003 |22,9834 |
|6,2 |7200 |0,0527 |0,1602 |4,1503 |25,8997 |
Table 5. Values used in this project for turbulence intensity I, kinetic energy k, dissipation e, and frequency w.
* [CFD Online, 2008] ** [OpenFoam PG, 2007]
The method of modelling the Reynolds stresses was described, and explains why the k-epsilon and k-omega equations are unable to properly simulate anisotropic flow conditions. The k-epsilon and k-omega model equations (two new equations for k and ε, or ω, are solved in addition to the RANS equations) were described, along with the comparing and contrasting the differences between the two turbulence models. Values of the relevant constants required for turbulence modelling were also provided for both turbulence models.
The final two sections describe the remainder of topics in this section covering computational fluid dynamics, and include descriptions of the finite volume method, numerical schemes, and solution algorithms, as relevant to this project, with details provided in Appendix A1: CFD Computations.
4 3.4 Finite Volume Method and Differencing Schemes
The mass, momentum, and scalar transport equations are integrated over all the fluid elements in a computational domain using CFD. The finite volume method is a particular finite differencing numerical technique, and is the most common method for calculating flow in CFD codes. This section describes the basic procedures involved in finite volume calculations. The computations are described in detail in Appendix A1: CFD Computations.
The finite volume method involves first creating a system of algebraic equations through the process of discretising the governing equations for mass, momentum, and scalar transport. To account for flow fluctuations due to turbulence in this project, the RANS equations (described in Appendix A1) are discretised instead when the cases are run using the k-epsilon or k-omega turbulence models. When the equations have been discretised using the appropriate differencing scheme for expressing the differential expressions in the integral equation (i.e. central, upwind, hybrid, or power-law, or other higher-order differencing schemes), the resulting algebraic equations are solved at each node of each cell.
The discretisation procedure is given in Appendix A1, while differencing schemes, numerical techniques, and the decisions made for numerical schemes used in this project are given in this next section. Solution algorithms and how velocity and pressure are coupled in order to compute the entire flow field are also described.
Differencing Schemes
The accuracy of the solutions depends on the differencing schemes used in the discretisation process. There are different differencing schemes available: central, upwind, hybrid, power law, and various higher-order differencing schemes. The central differencing scheme is described in some detail in Appendix A1. The upwind differencing scheme is the scheme used in the OpenFOAM simulations for this project.
Central Differencing Scheme
The central differencing method can be characterised as conservative, bounded at lower Peclet numbers, but neither bounded nor transportive at higher Peclet numbers. Since a direct linear average from the two faces is used for the property, the result of the differencing calculation will be incorrect since the flow is mainly in a single direction (as is the case with the air flow in the heat exchanger in this project), and therefore highly influenced by the property value just upstream from the node, and the property value should be closer to this value, rather than a simple average of the upstream and downstream value. In the case of central differencing, u must be very low, or the Peclet number must be less than 2 for the scheme to be considered bounded. Keeping this and the other mentioned considerations in mind, the central differencing method is not appropriate for this project.
Upwind Differencing Scheme
The upwind differencing scheme is designed to be capable of accounting for directionality of flow. It is used often in CFD codes for flow problems, which cannot be calculated with the central differencing scheme. If the flow is in one direction, for example the positive x-direction (with P node between West and East nodes), the property value at the xw face is taken to be the same as the value at the West node.
The upwind differencing scheme is conservative, bounded, and transportive. It contains a first-order Taylor series truncation error, which can result in a type of numerical error called false diffusion (over-estimation of diffusion) if the flow is not directly aligned with the grid, causing a diffusion-like transported property distribution, even where no diffusion is actually occurring. Making a finer grid improves the accuracy, but it is still not completely correct. For results of accurate flow calculations, the upwind differencing scheme is not always the best, and a higher order scheme for reducing discretisation errors can be chosen instead.
For this project, the upwind or hybrid scheme (described below, but not available in OpenFOAM) would be preferred to use. In OpenFOAM many schemes are available, and the upwind differencing scheme as described above is used. For continuation of this project, a higher-order upwind or hybrid scheme could be chosen.
Differencing scheme calculations are descibed in Appendix A1.
Solution Technique
When discretised equations are created for each node in the computational domain, using the relevant differencing scheme and modifying the first and last node to incorporate boundary conditions, the problem is ready to be solved using a matrix solution technique. The TDMA (tri-diagonal matrix algorithm) is based on forward eliminations followed by backward substitution and variations of this (i.e. iterative line-by-line methods for more complicated two- and three-dimensional problems) are used together with a SIMPLE (or variation) algorithm for linking pressure to velocity. Gauss-Seidel point iterative techniques are also available. The scalar property of interest is found by solving this system of linear algebraic equations.
The solution algorithm SIMPLE (semi-implicit method for pressure-linked equations), used to solve for the velocity field in all three directions and the pressure, is described in the next section.
[Hjertager, 2007] [Versteeg and Malalasekera, 2007]
5 3.5 Solution Algorithms
The value of the scalar properties of interest (i.e. temperature) at a particular location in the computational domain depends on the flow’s direction and velocity, which must also be solved for in the calculation process. To There are many algorithms available for this purpose, the most popular are the SIMPLE and PISO methods. This section describes the SIMPLE algorithm and compares it to the PISO algorithm. A short overview of the ‘staggered grid’ is also given.
In order to calculate the entire flow field, the momentum equations in all three directions and the continuity equation must be solved, which include terms for each velocity component and the pressure gradient. There is therefore a coupling between pressure and the three velocity components.
In cases where the pressure gradient is known, the usual discretisation techniques can be used to obtain discretised equations for velocity, since pressure can be calculated as any other scalar. However, there is no transport equation for pressure alone, and therefore pressure must be found using another method when the pressure gradient is not known. If the flow is compressible, the transport of density is found using the continuity equation and a scalar property such as temperature is found by the combination of the momentum equations and energy equation. Then the solutions for density and temperature are used for finding the pressure using the equation of state. When fluid is incompressible the density is not linked to pressure, a guess-and-check technique such as SIMPLE is required in order to solve the entire flow field.
The Staggered Grid
The flow and relative transport equations are discretised. Then it must be decided where the velocity component values will be stored. If they are stored at the same nodes as pressure and other scalar variables, a particular problem of a ‘checker-board’ pressure field arises due to linear interpolation of pressure during the process of discretising the equations. It results in the pressure having discretised gradients of zero at all nodal points, and the discretised momentum equations are not properly representing the pressure.
To solve the ‘checkerboard’ problem, the velocity components can be defined at the faces, while scalar variables (i.e. pressure, density, temperature) are stored at the usual nodal points. This results in the scalar variable control volumes (sometimes called the pressure control volume) being different from the velocity control volumes, and realistic pressure gradient results.
The SIMPLE Algorithm
The SIMPLE algorithm is a guess-and-correct technique to determine the values for pressure on a staggered grid. It is iterative and must be done in the specific order when other scalars are also calculated. The general procedure for the technique is shown in Figure 12, which is followed by a description of the steps in the algorithm.
[pic]
The steps in the SIMPLE algorithm, as shown in Figure 12, are summarized as:
START
Estimate a starting guess for the pressure field p*.
STEP 1
Solve the discretised momentum equations for the velocity components, using the line-TDMA method, and based on the pressure guess p*. For a three-dimensional case as the one in this project, 6 discretised momentum balances are solved for each of the six neighbors for node P (W, E, S, N, B, and T). This step results in finding values for u*, v*, w*, based on p*.
STEP 2
Solve for pressure correction equation to find the pressure correction p’. This is done with the discretised continuity equation by finding the mass imbalance bP and between total mass flow inflow and the total mass outflow of the six ‘guessed’ velocities (calculated based on the guessed pressure p*).
STEP 3
Correct the pressure and velocity components using the pressure correction, where the correction (p’) is added to the initial guess (p*), to get the new pressure field p, or:
Equation 38 [pic]
This step of the SIMPLE algorithm results in calculated values for velocity components and pressure: p, u, v, w, φ*, after correcting the guesses, which satisfy the continuity equation.
[pic]
STEP 4
Solve the other discretised transport equations using the line-TDMA method to get calculated values for the remaining scalar variables φ.
CONVERGENCE
After step 4, the outputs are tested for convergence (meaning that the mass imbalance is very close to zero), and if this is not within the value required for convergence, the program loops back to the beginning, using the newly calculated pressure, velocity, and other scalar values as the next starting guess. The process continues until convergence occurs (iteration).
Relaxation Factors
Often the pressure correction is too large, causing unstable calculations and divergence rather than convergence. Due to this problem, the iteration procedure must be slowed down to under-relax the pressure corrections. This is done by utilizing an under-relaxation factor[pic] of between 0 and 1, which is multiplied by the correction factor so that only a fraction of the originally calculated correction factor is actually used in the calculation, (for example with the pressure correction):
Equation 39 [pic]
Under-relaxation factors are also used for the velocity components. Oscillatory or divergent solutions are a result of too large values for [pic], while very slow convergence results from too small of values for the under-relaxation factor. Therefore, the correct under-relaxation factor is important for a converged solution, but cannot be determined in general and must be determined for each specific CFD case. In OpenFOAM, the default values are 0.3 for [pic] (for pressure), and 0.7 for the velocity. These were used for the cases solved using the simpleFoam solver, while for the rhoSimpleFoam cases (and the single rhoTurbFoam), it was necessary to reset these values to 0.2 for pressure [pic] and 0.5 for velocity, and in rhoTurbFoam to 0.1 and 0.3 for pressure and velocity correction factors, respectively.
There are also other versions of the SIMPLE method available: SIMPLER (SIMPLE-Revised) and SIMPLEC (SIMPLE-Consistent). The SIMPLER version involves using a discretised pressure equation rather than the SIMPLE correction equation used for pressure, while the velocity field is still calculated using corrections. The SIMPLEC version omits terms from the velocity correction equation.
[Hjertager, 2007] [Versteeg and Malalasekera, 2007]
The PISO Algorithm
The PISO (Pressure Implicit with Splitting of Operators) algorithm can be considered a continuation of the SIMPLE algorithm. PISO also involves making a guess for pressure, p*, and calculating the velocity components u*, v*, and w* based on this p*. There is a correction step to pressure and velocity fields, just as with the SIMPLE method. However, a second, additional, corrector step using the results from the first correction is carried out just after calculating the first correction. When the values for corrected pressure and velocity fields (now corrected twice) are determined, the other scalar transport equations are solved. The PISO algorithm requires more storage for the extra correction step, resulting in more computational resources needed.
[Hjertager, 2007] [Versteeg and Malalasekera, 2007]
Summary of Solution Algorithms
The SIMPLE algorithm is used for this project, in the solvers simpleFoam and rhoSimpleFoam. It is a known algorithm and used in many CFD codes. SIMPLER can be used to more efficiently use computer processing time (even though there are more calculations) since it uses the correct value for pressure. For certain types of flow, SIMPLEC and PISO can be just as efficient as SIMPLER. Which algorithm to use depends on the specific case being studied: the flow conditions, degree of dependence of momentum and scalars, and the specific numerical schemes used.
It is assumed that the air flow in the heat exchanger reaches a steady-state and does not contain fluctuations in time. Therefore steady-state calculations (rather than transient) are the focus. However, to investigate the possibility of this, one simulation was carried out and results looked as though they were steady-state. To calculate this situation, the solution algorithms can be calculated over time, in which case PISO is probably preferred, while SIMPLE can be converted to run transient calculations, and the continuity imbalance will contain an extra transient term to account for the time integrals of the variables, and an iteration process is used so that the solution converges for each time step (computationally expensive). In the case of this project, since steady-state behaviour is studied, the transient solution procedures are not described in detail.
4 Preliminary Results and Observations
This section presents some preliminary results and observations from the simulations carried out in this project. The simulations included ten different Reynolds numbers (based on ten inlet airflow velocities), all simulated in three flow models: laminar, k-epsilon turbulence model, and SST k-omega turbulence model. First the results of the grid independence test are presented. Then a description of the single transient case which was carried out in order to investigate if the air flow behaviour in this heat exchanger can be considered steady-state. Finally, some initial observations of the flow and temperature fields for cases from low and high Reynolds numbers are compared and contrasted. The next section (Section 5: Results) reviews the final results of the project.
1 4.1 Results of Grid Independence Test
There were a total of 11 different grid systems investigated to determine how fine the grid should be and to validate the solution independency of the grid. The mesh systems were developed using Salomé, as described in Section 2.3: Numerical mesh. First, unstructured tetrahedral grids were made, and then a more structured hexagonal mesh was created containing some unstructured areas around the tube walls (as seen in Figure 3). The meshes contained the following number of fluid elements (approximately): 8000, 25000, 50000, 75000, and 100,000. The predicted averaged pressure drops for the five grids are listed in Table 7 and illustrated in Figure 19.
Figure 13. Grid independence test for both tetrahedral and hexagonal meshes, with polynomial trend line illustrating hexagonal mesh trend; Re = 6200.
The 11 cases were run using the simpleFoam solver using the k-epsilon turbulence model and illustrated in Figure 13, with a trend-line created for the hexagonal mesh points. It can be seen in Figure 13 that the simulated pressure drop starts at a lower value, then stabilizes from total number of finite volumes from 25,000 to about 75,000, after which the pressure drop appears to increase again. The higher pressure drop calculated with the 150,000 cells may not be a stable value, since the
k-epsilon model cannot accurately predict the flow on fine grids near walls, since it is built to model free-stream flow. This may be the reason for a possible unstable pressure drop value from the simulation with a finer mesh. Had the mesh been sufficiently refined in the local areas around the tube walls and against the fins, and the appropriate turbulence model utilized (capable of calculating the wall effects), a grid independence test would not show an appreciable difference in the results (with respect to the total cell number) since the grid sensitive areas would have been taken into account [Song and Nishino, 2008].
As described previously, it proved difficult to use Salomé to refine the mesh only in the heat exchanger core model (locally refined mesh) without having computation problems due to the partitions in OpenFOAM, though this would have been a better use of computational resources (i.e. refine the mesh in the core and around the tubes, while making a coarser mesh for the extended volumes). The computer used for simulations in this study was not capable of efficiently creating in Salomé the meshes containing 150,000 or more cells either, and for these reasons, it was necessary to use a coarser grid without locally refining it), and the mesh of approximately 50,000 cells was chosen. This results of this grid lies in between the lower and upper bounds of the ‘stabilized’ pressure drop data (between grids with 25,000 volumes and 75,000 volumes).
Table 8. Grid independence test results.
|tetrahedral mesh |hexagonal mesh |
|number of cells |P drop |number of cells |P drop |
|8465 |33 |3000 |33 |
|25000 |37 |30000 |37 |
|50000 |40 |50000 |38 |
|75000 |38 |67000 |38 |
|100000 |40 |117000 |42 |
|150000 |49 | | |
2 4.2 Pressure Tracer in Transient Case
To confirm that the flow in the heat exchanger can be considered steady-state and the steady-state solvers can be accurately used for this project, a transient case was run using the ‘rhoTurbFoam’ solver in OpenFOAM (with inlet airflow velocity of 6.2 m/s). This was done in order to find out if there are fluctuation patterns in the airflow, making it an unsteady flow. In this case, a pressure tracer was written into the program to follow the changes occurring in pressure in ‘real-time’, and see if it is possible there is a pattern of fluctuations.
The first second of airflow through the heat exchanger is simulated over very small time-steps. To reach a solution more quickly, the final results for all fields from a steady-state k-epsilon case (with the same inlet airflow velocity) was used as the starting values for this transient case. The simulation is run for 1 second, and pressure tracers were written into the program at the inlet and outlet of the heat exchanger core model in order to track any possible fluctuations occurring in the short period of time. The values for the measurements recorded for pressure during the first 0.1 s. are illustrated in Figure 20. There are some fluctuations in the very beginning, but values quickly stabilize and very little fluctuation, if any, is detected at the 0.1 second mark. An animation was also carried out in ParaFoam, which also showed a steady-state was reached. Therefore it was determined the steady-state solvers could be used for solving this case.
Figure 14. Pressure tracer measurements at inlet and outlet for the first 0.1 s. of a transient case.
3 4.2 Characteristics of Flow
This section describes the initial observations found using paraView after running the CFD simulations in OpenFOAM. The characteristics of low-Reynolds flow and high-Reynolds flow are compared with contour plots of velocity with vectors (Figures 15 and 16), and contour plot of turbulent kinetic energy k (Figures 17 and 18).
[pic]
Velocity Observations
The flow patterns of the two cases at low and high-Reynolds numbers (inlet air velocity 0.3 m/s vs. 6.2 m/s) are similar. The air enters at the inlet on the left and flows in the direction of the arrows, and exits at the outlet on the right-hand side.
In both cases, as the air flows around the first tube, it begins to speed up and then the air velocity increases again as it goes around the second tube. This is verified by the samples taken in the case files for average velocities at the minimum free-flow areas, which showed that the velocity going around the second tube is faster than that going around the first tube. (Data is in Appendix A2). The minimum free-flow area is the area of the heat exchanger between two transverse tubes, so the area just above tube one or just below the second tube are the minimum free-flow areas. The flow is forced to speed up, as the tubes act as a type of pipe contraction in the air flow channel.
The highest velocity areas are just off the streamlines flowing directly around the tubes, and located at the area of minimum free-flow. In the case of the case with inlet velocity of 0.3 m/s, the top velocity at the second tube is 0.88 m/s, nearly 3 times the inlet velocity. For the 6.2 m/s inlet flow case, the top velocity reaches 15 m/s, more than twice the inlet velocity.
It is observed hat the size of the tubes impact the Reynolds number of the air flowing around them, since with larger tubes (at the same distance from each other), there would be an even smaller minimum free-flow area if the transverse pitch remained the same. In this study, the characteristic length for the Reynolds number is the tube collar diameter, and it can be seen here, that increases in this parameter (while keeping transverse pitch the same) can induce higher velocities and with it a higher turbulence and Reynolds number.
In the case of higher air flow, the recirculation zones behind each tube contain small backflow areas, (though difficult to see in Figure 16). The second recirculation zone appears larger. In the case with 0.3 m/s inlet velocity, however, a recirculation zone was not noticeable as it was for the case with inlet velocity of 6.2 m/s.
Kinetic Energy k distribution
The kinetic energy contour plots can be seen to verify previous observations made regarding flow. It is seen in Figure 17, which illustrates the kinetic energy k distribution for the low Reynolds number case. There is no kinetic energy increase in the areas behind the tubes for this case. The kinetic energy increases (slightly) in a different area corresponding to the increase in velocity as the air flows around the second tube. The other area of Figure 17, the plot of lower Reynolds number, which is exhibiting higher kinetic energy is in the area of higher temperatures, which can be seen from the graph in Figure 4.3. However, this illustrates that even at very low flow rates, some turbulent kinetic energy could still be present.
[pic]
For kinetic energy in the higher-Reynolds number case, an increase in kinetic energy is found clearly after the first tube, in the same area as the recirculation zones observed in the higher Reynolds values. According to this plot, then, the second recirculation zone is not as turbulent as the first recirculation zone. This makes sense because the direction of flow has changed as the air moves between the two tubes, and is directed more ‘downward’ at an angle (as shown by the vectors in Figure 15. The flow rounds the tube at an angle making less of an impact with the tube and ‘missing’ the recirculation zone.
Flow and kinetic energy plots have been compared and illustrated the effect flow has on kinetic energy. The next section describes the heat flow characteristics, as visualised in paraFoam for the same cases as for flow and kinetic energy. It is found that the temperature changes, flow, and kinetic energy can be shown to be connected.
4 4.3 Characteristics of Heat Transfer
The contours illustrating the local temperature distributions for the same cases as in the previous section are illustrated in Figures 19, 20, and 21, with the latter being a cross-section of the air flow channel through the thickness (z-direction) in the middle of the channel (it was ‘cut’ at y = 6.35 mm, from z = o to z = 2.24 mm, with the cut running along the x-axis from 0 to 44 mm – the length of the heat exchanger).
The first most noticeable difference between the two Reynolds number heat transfer characteristics is that once steady-state is reached, the slower-moving air (0.3 m/s) is heated up much more in the first two rows than in the case of higher Reynolds number flows. This must be due to the fact that the air flows so slowly, that there is much more time to absorb the heat (longer ‘residence time’). Had the initial inlet conditions been made cyclical, then comparisons could be made deeper into the heat exchanger (for example after 10 or 12 rows) and see how the heat transfer compares.
[pic]
Although streamlines are not physically drawn onto Figures 19 and 20, they can be seen fairly clearly with the color contrast lines. It is seen that the temperature streamlines run practically perpendicular to the velocity streamlines in the beginning of the airflow channel, with the isothermal streamlines running vertical and the velocity of the flow horizontal. This acts as a cross-flow heat exchange, with the flow directly bringing the heat with it. It can then be seen that after the air has flowed through the initial section of the heat exchanger, this ‘synergy’ between flow and heat transfer is no longer as effective. This means that the heat transfer coefficient is changing according to the streamline the flow is in at the time. Tao et al. (2006) has studied this ‘field synergy principle’ effect and reports that the intersection angle between the temperature gradient and the flow’s velocity is the main mechanism by which to increase the heat transfer.
It can be seen in Figure 20 that the higher Reynolds number flow has not only a lower temperature change than the previous example, but also a different pattern (different kinds of isothermal streamlines). The largest temperature changes for this case are occurring in the recirculation and ‘slow velocity’ zones (shown previously in the vector and velocity contour plot) just after each of the tubes. As in the slow-moving flow in the case with 0.3 m/s velocity, the slow-moving areas of the heat exchanger are also better able to absorb heat. The staggered tube arrangement is designed to have these slower-moving and recirculation areas to keep the heat flowing to the air, but at the same time, not allowing recirculation zones to ‘stagnate’ as can occur in inline arrangements where these zones do not keep flowing [Jang et al., 1995].
Figure 20 illustrate the contours across the z-direction in the middle of the airflow channel (which flows directly between the two tubes). It can be seen that the fins on the top and bottom are heated, and only a very thin ‘boundary’ layer of air has time to be heated after first entering the channel. It can be seen where the flow with heat from the first tube warms up the air, but then flows away again as all the air flow has to go around the second tube.
4.4 Summary of Preliminary Observations
From the initial qualitative observations of the graphs for flow velocity and direction, kinetic energy, and temperature fields, it is clear all are interconnected. The direction of flow streamlines compared to heat transfer isothermals determines how well heat is transferred, as illustrated in the grayscales for comparison in Figure 22. Kinetic energy is determined both by the flow velocity and temperature, and pockets of slower-moving and re-circulating air after flowing by tubes transfers the bulk of the heat in faster-moving flow.
5 Numerical Results and Discussion
This section concludes the project with the final results presented for the simulations carried out on the two-row fin-and-tube heat exchanger. Introductory simulations were carried out in the ‘simpleFoam’ solver without temperature changes calculate to get an idea of how the Salomé and OpenFOAM software work together and to see if the calculations were in a ‘ballpark’ range of the experimental values for different Reynold’s numbers before starting a real simulation plan. These initial simulation results were not expected to be accurate since the density differences in the cold and heated air are not accounted for in the ‘simpleFoam’ solver. Then simulations with ‘rhoSimpleFoam’ were carried out using a laminar flow model and two turbulence models: k-epsilon and SST k-omega. The results of these simulations were used to calculate the Fanning friction factor f and Colburn j-factor to characterize the capability of each model to accurately predict pressure drop and heat transfer in the heat exchanger. First the results for friction factor are given (both for the initial studies and those modelling temperature as well), followed by Colburn j-factor results and discussion.
1 5.1 Friction Factor
Initial Simulation with ’simpleFoam’ Solver.
The first simulations were carried out to find out how to get the Salomé meshing program and OpenFOAM CFD program to work and learn to use the files as related to this project. Temperature was not modelled in these, and therefore results were not expected to be very accurate. For these simulations, only the friction factor was calculated and compared with experimental data. However, as can be seen in Figure 23, it was found that the patterns in the graph were somewhat similar for the experimental and simulated values, and that the flow models followed the same pattern, with a gradually decreasing friction factor as Reynolds number increases. Also, the different flow models achieved nearly the same results. The simulations with temperature included in the calculations were then run to see how the simulations actually compared with the experiments.
[pic]
Figure 23. Initial ’trial’ simulations in OpenFOAM using ‘simpleFoam’ solver.
Simulations with ’rhoSimpleFoam’ Solver: Fanning friction factor results
The values for the friction factor f for the 10 test samples are plotted against Reynolds Number (ReDc) in Figure 24.
[pic]Figure 24. Fanning friction factor f against Reynolds number Re for different inlet airflow velocities and
flow models (laminar, and turbulence models k-epsilon and k-omega) with the same geometrical parameters.
It can be seen from Figure 24, that in all cases the friction factor was decreasing with increasing Reynolds number. All of the models underestimated the friction factor, including the transient rhoTurbFoam case. At the lower laminar flows (the first four points at the lower Reynolds number), of Reynolds number from 330 to 1300, all the models found nearly identical results. As the flow moved into transition, it appears the laminar flow model came much closer to the experimental values. At the transition point from transition to turbulent, which appears to have a critical Reynold’s value of between 1700 and 2900 (the exact value is not known, since there weren’t enough data points given to be sure), once again, none of the models were better than another. After moving into turbulent flow, however, the k-Omega SST had the best accuracy, and the laminar flow model able to model the friction factor compared to the k-Epsilon model.
[pic] Figure 24b. Fanning friction factor f: Error by flow model vs. Re.
The results fit reasonably with what how the different models calculate the flow. Both the k-omega and k-epsilon are two-equation models created for calculate turbulence, turbulent kinetic energy, etc., and therefore in laminar flow the additional turbulence terms do not increase the accuracy since there is no turbulence to model in laminar flow anyway. The k-epsilon model was the least accurate flow model of them all, and this is probably due to that the equations only model kinetic energy and dissipation and are accurate for free-flowing fluids, and therefore the friction factor against the wall is not capable of being accurately calculated.
2 5.2 Colburn j-Factor
The heat transfer characterisation parameter Colburn j-factor has been calculated from the simulation results for the different flow models. The error is shown in Figure 25b. The flow models showed very clear differences in abilities to simulate heat transfer at the different Reynol ds numbers. In laminar flow, the laminar flow model was the best for predicting the j-factor, as would be expected. The transition heat flow was best characterized with the k-omega turbulence model, while turbulent heat flow was best calculated using the k-epsilon model, although at the very highest Reynolds number, 7000, none of the models were accurate.
[pic]Figure 25. Colburn j-factor against Reynolds number Re for different inlet airflow velocities and
flow models (laminar, and turbulence models k-epsilon and k-omega) with the same geometrical parameters
[pic]
Figure 25b. Colburn j-factor error: Error by flow model vs. Re.
3 5.3 Summary of Results
Preliminary simulations were run to see how the simulations are run, and if the results are in ballpark range of the experimental values. It was found the friction factor was not extremely far from the experimental values. However, as there are no heat flow in the simulations (while there was heat flow in the experiments), the solutions were not expected to match very closely anyways.
Simulations were then run to determine friction factor f and Colburn j-factor. The results found that different flow models performed best according to Reynolds number and whether accurate solutions were desired for pressure drop or heat transfer.
The optimal conditions for simulating the two-row fin-and-tube heat exchanger are summarized below. First, the optimal simulation conditions for determining the friction factor are shown, and then Colburn j-factor after that. The error as found for flow model vs. Reynolds number is also shown.
Simulation conditions for determining Fanning friction factor f.
|Flow regime |Flow model |Percent Error (vs. Re) |
|Laminar Re 330-1300 |No optimal model found |n/a |
|Transitional Re 1300-2900 |Laminar flow model |-0.5 % |
|Turbulent Re 2900-6200 |SST k-omega turbulence model |4 % - 11.5 % |
|Flow regime |Flow model |Percent Error (vs. Re) |
|Laminar Re 330-1300 |No optimal model found |n/a |
|Transitional Re 1300-2900 |SST k-omega turbulence model |2.4 – 17.8 % |
|Turbulent Re 2900-6200 |k-epsilon turbulence model |1.1 % - 9.4 % |
6 Discussion
Simulations for this project were carried out following as closely as possible the same operating conditions and geometrical configurations of the two-row tube-fin heat exchanger, with tube collar diameter of 10.23 mm and fin pitch 2.23 mm, as presented in the paper by Wang et al. (1996). The Reynolds number ranges from 330 to 7000, which correspond to the frontal air velocity at the inlet ranging from 0.3 to 6.2 m/s.
The work done for this project has shown that it is possible to make practical simulations of heat flow and pressure drop for a tube-and-fin heat exchanger using open source CFD software, and validate the results against experimental data. Data resulting from the simulations should be as accurate as possible, and therefore some considerations can be taken in future work to attempt to further improve the simulation conditions/calculations and the accuracy of the results. These improvements could include changes to the following areas of CFD simulation:
• A more comprehensive grid independence test.
• Changes to the fin temperature based on fin effectiveness calculations. The efficiency equation given in Baggio and Fornasieri (1994) assumes a uniform air and fin temperature, which is not the case in this project. As shown by Ay et al. (2001) and results from infrared thermography measurements, the local convective heat transfer coefficient changes across the fin according to various parameters. It was shown that there is a lower temperature at the leading edge of the plate-fin, and a sharper temperature gradient on the fin surface where the boundary layer increases and destroyed as the fluid flows around the tubes (for the first two rows of tubes). Once the flow has gone around the tubes, the temperature gradient decreases from airflow being swept into the wake. However, by the third row of tubes, the wake pattern changes again, and further variations of the heat transfer coefficient can be seen throughout the heat exchanger, with specific patterns depending on the Reynolds number (among other parameters). The article included studies of staggered and in-line fin-and-tube heat exchangers. The temperature gradients can also be studied in relation to the synergy principle presented by Tao et al. (2007), where the heat transfer coefficient is shown (qualitatively) to change according to the angle of local isothermal streamlines to the temperature field. The general pattern of heat transfer coefficient in this paper is similar to that described in the Baggio and Fornasieri (1994) paper.
• Clarification of the air temperature at inlet and the hot water flow rate through the tubes, as this information was not provided in the Wang et al. (1996) article.
• Improve the final mesh to be used (structured vs. unstructured or hybrid, determine areas of geometry requiring finer mesh, etc.). Versteeg and Malalasekera (2007) suggest that non-structured grids can be more accurate and efficient than structured grids.
• Use a solver for conjugated heat transfer analysis to include the interactions between the air, tube wall, and water, which has recently become available in the OpenFOAM version 1.5.1.
• Use cyclic boundary conditions for inlet flow to investigate pressure drop and heat flow characteristics deeper into the heat exchanger.
• Use the new geometry and mesh creation program snappyHexMesh available now in OpenFOAM (instead of Salomé, which proved difficult to use at times, and no technical support was available).
• Use a low Reynolds number turbulence model to better simulate the turbulence at the lower Reynolds numbers not accurately modelled by any of the flow models.
• Make more use of the openFOAM discussion boards and online information, since no technical support is available. (Learn the OpenFOAM more thoroughly).
• Possible changes to the turbulence model in OpenFOAM, or solving procedures – which is possible since it is open source C++, and changes only require basic programming skills in object-oriented programming, making it relatively simple to implement new turbulence models, solver algorithms, boundary condition types, and physical models. This is an advantage over commercial software, where access to the code is unavailable.
• Finally, as discussed in the book by Kays and London (1998), the properties of air are highly temperature-dependent, and many of the calculations do not account for these changes, but instead use an average value, which can substantially affect the flow at a particular cross section according to the temperature profile (for example in this case, flow characteristics were determined using average temperatures taken at the inlet and outlet).
• Run the steady-state versions of k-omega turbulence models further to see if they can converge better, since the curve of Colburn j-factor vs. Reynolds number seems to be unstable with a noticeable fluctuation at the inlet flow velocity of 3.7 m/s (corresponding to Reynolds number 4300).
• Implementation of anisotropic turbulence models to correct for the differences in flow according to the direction, i.e. use the realizable k-epsilon turbulence model or use of the RSM (Reynolds stress equation model) turbulence models.
As can be seen from the preceding list, which does not consider all the possible aspects, there is much to be considered for ensuring accurate simulations of the fin-and-tube heat exchangers. To summarize, considerations should be taken for: heat exchanger geometry and mesh, fin temperature, boundary conditions (air and water temperatures, and cyclic inlet), turbulence model variations, OpenFOAM use and programming, convergence, and temperature-dependent properties of air. All of these considerations are subjects of interest which can be studied in the open literature as described in this paper and listed in the references.
7 Conclusion
The objective of this project was to make CFD simulations using open source software, and validate the results against experimental data. The system to study was a fin-and-tube heat exchanger. The purpose of the work was to investigate the possibilities of eventually using CFD calculations for design of heat exchangers instead of expensive experimental testing and prototype production.
To analyse the flow and heat transfer characteristics of the heat exchanger, a model of a two-row fin-and-tube heat exchanger was created using open source Salomé software to create the geometry and mesh. The resulting mesh (after a grid independence test was carried out) was used for running a variety of simulations using a laminar flow model and two turbulence models for comparison of results. Ten different inlet flow velocities ranging from 0.3 m/s to 6.2 m/s and corresponding to Reynolds numbers ranging from 330 to 7000 were simulated in the three different flow models (laminar, k-epsilon turbulence model, and SST k-omega turbulence model). A sampling dictionary was written into the CFD model to record pressure and temperature measurements at the inlet and outlet of the heat exchanger model. Using the simulation results and some specific non-dimensional numbers, calculations related to heat flow and pressure loss can be carried out to determine the Fannning friction factor and Colburn j-factor for comparison with the literature values used for the validation.
It was found that the flow model accuracy depended on the flow regime and whether the friction factor f or j-factor was being determined. From the experimental values given in the literature, the laminar flow region for this particular geometry of heat exchanger switched to transitional at around Reynolds number 1300, and moving to transitional around Reynolds number 2900. The Reynolds number has a characteristic dimension of the tube collar outside diameter.
For friction factor determination, little difference is found between the flow models simulating laminar flow, while in transitional flow, the laminar flow model produced the most accurate results (for friction factor) and the SST k-omega turbulence model was more accurate in turbulent flow regimes. For heat transfer, the laminar flow model calculated the most accurate j-factor, while for transitional flow the SST k-omega turbulence model was more accurate and the k-epsilon turbulence model was best for heat transfer simulations of turbulent flow.
The flow model can be chosen based on what is being studied (heat flow or pressure drop) and the flow regime. In conclusion, it is found that the pressure drop and heat transfer characteristics of a fin-and-tube heat exchanger can be determined to within a reasonable accuracy with CFD computations carried out in open source software, and that OpenFOAM can be used to carry out practical work in the design process of heat exchangers.
8 References
ASHRAE. “Handbook Fundamentals”, SI Edition (1993).
Ay, Herchang; Jang, Jiin Yuh; Yeh, Jer-Nan. “Local heat transfer measurements of plate finned-tube heat exchangers by infrared thermography”, International Journal of Heat and Mass Transfer, Vol. 45 (2002), pp. 4069-4078.
Baggio, P.; Fornasieri, E. “Air-side heat transfer and flow friction: theoretical aspects”, in Recent developments in finned tube heat exchangers. Energy Technology (1994) pp. 91-159
CFDonline. wiki/Turbulence_intensity, accessed June 2008.
Chen, Han-Taw; Chou, Juei-Che; Wang, Hung-Chih. “Estimation of heat transfer coefficient on the vertical plate fin of finned-tube heat exchangers for various air speeds and fin spacings”, International Journal of Heat and Mass Transfer, Vol. 50 (2006) pp. 45-57.
Erek, Aytunc; Õzerdem, Baris; Bilir, Levent; Ilken, Zafer. ”Effect of geometrical parameters on heat transfer and pressure drop characteristics of plate fin and tube heat exchangers”, Applied Thermal Engineering, Vol. 25 (2005) pp. 2421-2431.
Fornasieri, E.; Mattarolo, L. “Air-side heat transfer and pressure loss in finned tube heat exchangers: state of art”, Proceedings of the European Conference on Finned Tube Heat Exchangers, Lyon, France, (April 1991).
Gnielinski, V. “New Equation for heat and mass transfer in turbulent pipe and channel flow”, International Chemical Engineering, (1976) 359-368.
Gray, D. L.; Webb, R.L. “Heat transfer and friction correlations for plate fin-and-tube heat exchangers having plain fins”, Proceedings of the Ninth International Heat Transfer Conference, San Francisco (1986).
Hjertager, Bjørn H. “Turbulence Theory and Modelling”, Lecture Notes, Aalborg University Esbjerg, Denmark (2005).
Hjertager, Bjørn H. “Computational Analysis of Fluid Flow Processes”, Lecture Notes, Aalborg University Esbjerg, Denmark (2007).
Hjertager, Bjørn H. “Introduktion til Open Source CFD beregninger I OpenFOAM”, Course Notes, Aalborg University Esbjerg, Denmark (2008).
Jang, Jiin-Yuh; Wu, Mu-Cheng; Chang, Wen-jeng. “Numerical and experimental studies of three-dimensional plate-fin and tube heat exchangers”, International Journal of Heat and Mass Transfer, Vol. 39, No. 14 (1996) pp. 3057-3066.
Kayansayan, N. “Heat transfer characterization of plate fin-tube heat exchangers”, International Journal of Refrigeration, Vol. 17, No. 1 (1994) pp. 49-57.
Kays, W.M.; London, A.L. “Compact Heat Exchangers”, Sub-edition 3, Krieger Publishing Company, New York (1998).
Mangani, L.; Bianchini, C.; Andreini, A.; Vacchini, B. ”Development and validation of a C++ object oriented CFD code for heat transfer analysis”, ASME-JSME 2007 Thermal Engineering and Summer Heat Transfer Converence,Vancouver, Canada (July 2007).
McQuiston, F. C. “Correlation for heat, mass and momentum transport coefficients for plate-fin-tube heat transfer surfaces with staggered tube”, ASHRAE Trans. Vol. 84 (1978), pp. 294-309.
OpenFOAM discussion board. , accessed July 2008.
OpenFOAM documentation. “Programmer’s Guide - PG” and “User’s Guide - UG”, OpenCFD Limited (2007).
Perrotin, Thomas. “Fin efficiency calculation in enhanced fin-and-tube heat exchangers in dry conditions”, International Congress of Refrigeration, ICR0026 (2003).
Rocha, L. A. O.; Saboya, F. E. M.; Vargas, J. V. C. “A comparative study of elliptical and circular sections in one- and two-row tubes and plate fin heat exchangers”, International Jouernal of Heat and Fluid Flow, Vol. 18 (1997), pp. 247-252.
Rohsenow, Warren; Hartnett, James P.; Cho, Young I. “Handbook of Heat Transfer”, 3rd Edition, McGraw-Hill Companies, Inc., New York (1998).
Song, Gil-Dal; Nishino, Koichi. “Numerical Investigation for Net Enhancement in Thermal-Hydraulic Performance of Compact Fin-Tube Heat Exchangers with Vortex Generators”, Journal of Thermal Science and Technology, Vol. 3, No. 2 (2008) pp. 368-380.
Sahin, Haci Mehmet; Dal, Ali Riza; Baysal, Esref. “3-D Numerical study on the correlation between variable inclined fin angles and thermal behaviour in plate fin-tube heat exchanger”, Applied Thermal Engineering, Vol. 27 (2007) pp. 1806-1816.
Schmidt, Th. E. “Heat transfer calculations for extended surfaces”, Refrigeration Engineering, April (1949) pp. 351-357.
Tao, Y. B.; He, Y. L.; Huang, J.; Wu, Z. G.; Tao, W. Q. “Three-dimensional numerical study of wavy fin-and-tube heat exchangers and field synergy principle analysis”, International Journal of Heat and Mass Transfer, Vol. 50 (2007), pp. 1163-1175.
Tennekes, H.; Lumley, J. L. “A First Course in Turbulence”, MIT Press, Cambridge, MA (1972).
Versteeg, H K; Malalasekera, W.”An Introduction to Computational Fluid Dynamics, The Finite Volume Method”, Second edition, Pearson Education Limited, Essex, England (2007).
Tutar, Mustafa; Akkoca, Azize. “A computational study of effects of different geometrical parameters on heat transfer and fluid flow in a wavy and plain fin and tube heat exchanger”, Proceedings of ESDA2002: 6th Biennial Conference on Engineering Systems Design and Analysis, Instanbul, Turkey (July 2002).
Yan, Wei-Mon; Sheen, Pay-Jen. “Heat transfer and friction characteristics of fin-and-tube heat exchangers”, Volume 43 (2000), pp. 1651-1659.
Wang, Chi-Chuan; Chang, Yu-Juei; Hsieh, Yi-Chung; Lin, Yur-Tsai. “Sensible heat and friction characteristics of plate fin-and-tube heat exchangers having plane fins”, International Journal of Refrigeration, Vol. 19, No. 4 (1996) pp. 223-230.
Appendix
1 A1 CFD Computations
1 A1.1 CFD Governing Equations
An illustration of a fluid element for CFD calculations is shown in Figure 22. The element has dimensions δx, δy, and δz, with the center point at (x, y, z) and six faces N, S, E, W, T and B (North, South, East, West, Top, Bottom). Each fluid property (velocity, pressure, density, viscosity, thermal conductivity, and temperature) therefore can be represented as a function of space and time with: [pic](x, y, z, t), p(x, y, z, t), ρ(x, y, z, t), µ(x, y, z, t), k(x, y, z, t), and T(x, y, z, t).
As summarized above, the fluid properties to be calculated using CFD are the 3 velocity components, pressure, density, viscosity, thermal conductivity, and temperature. In all, there are 8 variables. The 8 equations therefore needed for solving these are: (1) mass balance, (2)-(4) momentum balance in 3 directions, (5) energy equation, (6) equation of state, (7)-(8) empirical relations describing viscosity and thermal conductivity. Each equation is presented here.
1) Mass balance (Continuity equation): “Fluid mass is conserved.”
Equation 40 also written as [pic]
or [pic] which can be summarised as:
[mass accumulation over time] = [sum of all inflows] – [sum of all outflows]
where δ represents area, and can be illustrated in the following diagram (Figure 8) of an element representing the one illustrated in Figure 7.
Before presenting the remaining 7 equations, a brief illustration of how equation 40 was derived is given here. The momentum and energy balances are derived similarly for a fluid element. Therefore instead of writing them out in detail, only the summary of equations / expressions will be given, and for details of the derivation, similarities can be assumed with the following mass balance illustration shown here.
Consider an element similar to that shown in Figure 8. The fluid element’s rate of mass increase is:
Equation 41 [pic]
Equation 28 represents the rate of increase of mass over time within the fluid element, where ‘δ’ represents the distance between faces for the specific x, y, or z direction. The mass flow rate across the control volume faces must now also be accounted for.
The net rate of flow is the sum of mass inflow subtracted by the sum of mass outflow. The first two terms of a Taylor series expansion can accurately express fluid properties at the faces. Therefore, the mass flow in the x-direction through the W and E faces (at a distance of (1/2)*(δx) from the center of the element) is expressed as:
[pic] for the west face W, and [pic] for the east face E.
The mass flow in the y-direction through the S and N faces, and in the z-direction through the B and T faces can be similarly expressed. All of these are illustrated in Figure 24 and summarized after the illustration using Equation 42.
[pic]
Figure 28. Fluid element illustrating flows of inflows and outflows of mass on all six faces.
[Versteeg and Malalasekera, 2007]
As can be seen in Figure 24, the overall mass flow rate across the element’s faces is represented by the following expression, in which the entire control volume is taken into account by multiplying the mass rate in a particular direction by the two remaining dimensions:
Net mass flow =
Equation 42 [pic]
Equation 41 (rate of increase of mass inside the element) is equated to Equation 42 (net mass flow rate into the control volume across its faces). The terms are arranged on the left side of the equation, and divided by the control volume[pic] to get the continuity equation for compressible fluids (Equation 27):
If the fluid is incompressible, the density is constant and equation becomes div [pic] = 0.
(2-4) Newton’s 2nd Law: Momentum balances (3 equations for each direction)
“The rate of increase in momentum is equal to the sum of the forces acting on the body.”
There are three equations for the momentum rates of increase in each of the three x, y, and z directions. There are surface forces (pressure and viscous) and body forces (gravity, centrifugal, Coriolis, and electromagnetic) included in the momentum equations. The surface forces are stated explicitly, while relevant body forces are represented in the F term of each equation. The use of the total derivative should be noted, representing convective differentiation on the left side of the equations above, i.e.:[pic].
(2) x-component momentum equation
Equation 43 [pic]
(3) y-component momentum equation
Equation 44 [pic]
(4) z-component momentum equation
Equation 45 [pic]
(5) Energy balance: First Law of Thermodynamics
“The energy change of rate in a fluid element equals the heat addition rate
in the fluid element added with the rate of work done on the element.”
The energy equation for net energy (e) flow includes terms to account for energy flux from conduction, work done by surface forces (i.e. shear stress causing friction and heat loss – irreversible, or change in volume due to pressure), and finally a source term to account for chemical reactions or radiation (either an energy source or sink).
Equation 46 [pic]
6) Equation of state
Fluids are nearly always in thermodynamic equilibrium. When two state variables of a substance are known, an equation of state can be used to describe the other thermodynamic variables (these are pressure p, density ρ, specific internal energy i, and temperature T. The (simplified) perfect gas equation of state is:
Equation 47 [pic]
Density changes as a result of temperature (or pressure) variations in the flow field of compressible fluids. This equation links the mass conservation and momentum equations to the energy equation, making it possible to calculate temperature changes, for example, from changes in density when two of the state variables are known (i.e. density and pressure).
In the case of incompressible fluids, the energy and mass/momentum equations cannot be linked, and the mass and momentum equations must be used to solve for the flow field, and it is not necessary to solve for the energy equation.
In this project, the simpleFoam solver only solves the flow field, while rhoSimpleFoam solves for flow and temperature. The equation of state is used in the calculations to solve for the temperature field for the rhoSimpleFoam cases, but not used in simpleFoam calculations.
7) Newton’s law of viscosity (empirical relation for viscosity µ)
Equation 48 [pic]
Equation 48 can then be inserted into the momentum balance to solve for viscosity.
8) Fourier’s second law of thermodynamics (empirical relation for thermal conductivity k
Equation 49 [pic]
Equation 49 can then be inserted into the energy equation (Eqn. 33) to solve for viscosity.
General Transport Equation
From the insertions into the momentum and energy balances, the general transport equation for property [pic] is:
Equation 50 [pic]
which also can be expressed as:
[pic]
The first term (I) on the left side of Equation 50 represents the rate of change of property [pic] in the control volume over time [accumulation]. The second term (II) on the left side is the convective term [convection], or the net rate of change in [pic] in the control volume due to flow of [pic] across the boundaries of the volume. The first term (III) on the right side represents a diffusive term where [pic] is the diffusion coefficient (i.e. k or μ) [diffusion]. The final term (IV) is the source term for the net rate of production of the property [pic] (i.e. in the case of heat-producing chemical reactions) inside the control volume [internal production source].
The general transport equation is the basis for the CFD finite volume calculations. The property of interest (i.e. u, v, w, or T) is inserted into the equation as[pic], the appropriate diffusion coefficient selected for the variable[pic], and relevant source terms are included. Then each of the mass, momentum, and energy conservation equations are transformed into equations which can be used in the finite volume method to make the CFD computations. The equations are integrated over all the control volumes in the geometry of interest, taking the solution from one fluid element and using it as a start value for the next element.
[Hjertager, 2007][Versteeg and Malalasekera, 2007]
2 A1.2 Reynold’s Averaged Navier Stokes (RANS) Equations
This section describes the RANS equations, which form the basis for turbulence modelling.
Time-Averaged Properties
The variation in the velocity u in the x-direction as a function of time (transient flow) at a particular point in a fluid element in turbulent flow may be illustrated as in Figure 7 below.
[pic]
Figure 29. Transient variations of velocity u in turbulent flow. [Hjertager, 2005]
As seen in Figure 25, there are many small fluctuations occurring in turbulent flow. The fluctuations are due to vertical eddy motions which create strong mixing and a momentum exchange that causes acceleration of slower moving layers and deceleration of faster moving layers. Because of this momentum exchange, there are shear stresses, which are known as Reynolds stresses, from the turbulence. In addition to velocity fluctuations, concentration or heat fluxes can also occur over the fluid element faces.
Most engineering purposes do not require solving for all of the details of the property fluctuations due to eddy motion. Therefore time-averaged properties of the flow are normally used for making calculations. In this case (Figure 25), velocity is characterized using the time-averaged value for velocity u, expressed as ‘[pic]’, while the individual measured variations from the mean are expressed as ‘[pic]’. However, in turbulent flow the other flow properties, designated as φ (i.e. pressure, temperature, density, other velocity components, etc.) also have fluctuations. Therefore, the variables can be generalized as a steady mean component [pic] (i.e., representing[pic]), which can be defined with the following expressions:
Eqn 51 a [pic] Eqn 51b [pic] Eqn 51c [pic] Eqn 51d [pic]
• Equation 51a describes the actual value for property φ to be the sum of the steady mean component of φ ([pic]) and the fluctuating time varying component [pic].
• Equation 51b describes the mean property [pic] as the time average of φ.
• Equations 51c and 51d illustrate that the time average of the fluctuations equals 0, while the variances of the property fluctuations does not equal zero (also second moments made from pairs of different variables (i.e. [pic]) are not equal to zero either.
[Hjertager, 2005] [Versteeg and Malalasekera, 2007]
Time-Averaged Navier-Stokes Equations
The Navier-Stokes equations for mass (Equation 27) and momentum (Equations 30-32) are expanded to include the time-averaged properties affected by the turbulence. As an example, the process for obtaining the time-averaged mass balance is described in some detail. The momentum balances and scalar transport equation follow a similar process (although more complicated), and the resulting time-averaged equation only is presented for these equations.
Time-Averaged Mass Balance
The continuity equation is:
Considering the density ρ and velocity components u, v, and w as a sum of steady mean component [pic] and fluctuating time-varying component [pic], the values can be plugged into the continuity equation:
Equation 52 [pic]
Now replacing with the time-averaged terms, Equation 39 is then changed into the following:
Equation 53 [pic]
which can also be expressed as: [pic]
when the flow variables [pic], (u, v, and w) are replaced with:
[pic]= [pic] + u u = U + u’ v = V + v’ w = W + w’ p = P + p’
to account for the effect of fluctuations on the mean flow.
Time-Averaged Momentum Equations
A process similar to that for the continuity equation above (replacing terms with time-averaged terms) is carried out for the momentum equations (Equations 30-32). The resulting equation for the momentum equation in the x-direction is given as:
Equation 54 [pic]
The equations for the other directions (y and z) are found using the same process, but not listed here. Density is assumed constant, and the new terms ( ) appear as a result of the process of time averaging. The additional terms can be interpreted as extra turbulent stresses, known as the Reynolds stresses, to the mean velocity components (U, V, W), and are due to turbulent eddies causing fluctuating velocities and convective momentum. The Reynolds stresses are defined with the stress tensor (consisting of 3 normal stresses and 3 shear stresses):
Equation 55 [pic] , and in matrix form: = [pic]
The stress tensor [pic] is modelled prior to solving Equation 41, by relating the stresses to the time-averaged velocity using turbulence models (described in Section 3.4.2: k-εpsilon Turbulence Model).
Time-Averaged Scalar Transport Equation
When deriving the time-averaged scalar transport equation (for scalar quantities like temperature), additional turbulent transport terms ( ) also appear.
Equation 56 [pic]
The RANS equations for time-averaged mass, momentum, and scalar transport equations have now been presented. They are used in this project’s turbulence models to calculate the flow properties.
Kinetic Energy
In order to make CFD calculations including effects of turbulence, the total turbulent kinetic energy k per unit mass at a particular point in a fluid element, and the turbulence intensity I, are determined for use in modelling of the turbulence. These are defined with the following expressions:
Equation 57 [pic]=[pic]
Equation 58 [pic]
For the k-epsilon and k-omega turbulence models used in this project, the velocity fluctuations are estimated and k calculated using Equation 57. Turbulence intensity can then be found using the calculated k value and the reference mean flow velocity Uref using Equation 58.
[Versteeg and Malalasekera, 2007] [Hjertager, 2005]
3 A1.3 Finite Volume Method and Finite Differences
Equation Integration and Discretisation
A steady-state one-dimensional heat transfer through diffusion is governed by the general transport equation (Equation 59). When the transient and convective terms are deleted, the equation becomes:
Equation 59 [pic]
Integration of Equation 59 yields:
Equation 60 [pic]
In the one-dimensional steady-state heat transfer diffusion problem, the equation becomes:
Equation 61 [pic]
where enthalpy h = φ =[pic] ,
[pic] , and
[pic]
It can be seen from the expressions above, that the enthalpy is the scalar property to be transported. The transport coefficient for [pic] enthalpy is the thermal conductivity k divided by the specific heat capacity [pic] at the average temperature of the cell. The rate of heat transfer by diffusion (term (III) in the transport equation), is then:
[pic], and the transport equation to use is
Equation 62 [pic]
Equation 62 is the one-dimensional energy equation, to be discretised and solved for diffusion.
Grid Generation
To discretise the energy equation over the fluid elements in a computational domain, the following single fluid element is illustrated:
In the heat diffusion problem, heat flow in three dimensions (N–S, W-E, T-B) is explained. Therefore, one row of fluid elements, each with a central node like the one above, is illustrated.
| | | | | |
| | | | | |
|[pic] |[pic] |[pic] |[pic] |[pic] |
| | | | | |
Figure 31. One-dimensional grid with nodal points.
Discretisation
In three-dimensional calculations, a similar set-up is made for the N, S, T, and B faces and corresponding central nodes. For a one-dimensional calculation, the integration from the west and east faces on either side of node P is carried out according to equation 60, resulting in Equation 63:
Equation 63 [pic]
where A is the cross-sectional area of the fluid element,
ΔV is the volume and
[pic] the average source in the fluid element.
As seen in Equation 63, the governing equation is integrated across the fluid element, with node P described with a discretised equation.
Temperature gradients at xe and xw must be known for Equation 63 to be useful. The temperature T (or other scalar property φ) and thermal conductivity k (or diffusion coefficient Γ) are determined for the nodes, and therefore the property value gradients at the faces xe (halfway between nodes P and E) and xw (between nodes W and P) between must be approximated. There are different types of differencing schemes for this purpose.
Central Differencing Scheme
Direct linear averaging (assuming uniform grid) of the values can be used, which is termed the central differencing method. The other differencing methods are compared in the next section (upwind differencing, hybrid, or power-law, and other higher-order differencing schemes), and the choice made for this project given. When the differencing method is chosen (in this case central differencing), the equation takes the form:
Equation 64 [pic]
The equation is then rearranged with all P scalar variables on the left side, and W and E variables on the other, with specific terms given the names aP, aW, aE, and b, as shown in the following equation:
Equation 65[pic]
[pic]
The discretized equation with the representative terms [pic], [pic], [pic], and [pic] is then:
Equation 66 [pic]
for a three-dimensional problem, a similar
expression is found for the N, S, T, and B faces:
[pic][pic]
If the problem includes convection in addition to diffusion, the coefficients for a will include an additional term to account for the convection. For example, F could represent convective mass flux (ρu), while D represents the diffusion conductance (k/δx) resulting in [pic] now being equated to [Dw+(Fw/2)].
The central differencing scheme does not always reach correct solutions when both convection and diffusion forces are involved. To determine how realistic a differencing scheme is and ability to reach a converged solution, the following terms are used:
• Conservativeness: The flux of a property leaving one face in the computational domain should always be the same as the flux entering the adjacent cell (across the shared face).
• Boundedness: The calculated values for a scalar property should always lie within the boundary values (bounded), and all coefficients a (i.e.[pic],[pic], etc.) should be of the same sign (this way an increase or decrease in the value at a particular node should result in an increase or decrease, respectively, of the property value at its neighbouring nodes).
• Transportiveness: This is described with the Peclet number, a ratio of the relative strength of convection to that of diffusion, as in the following expression:
Equation 67 [pic]
If there is only convective flow involved, then the Peclet number approaches infinity and the values of properties at neighboring nodes are strongly influenced by the direction of flow. However, if only diffusion is occurring, the Peclet number approaches zero, and the values of the property spreads at the same rate in all directions, rather than in a single direction. It can be seen from Equation 67, that the Peclet number depends on fluid properties density and diffusive coefficient, flow, velocity, as well as the grid itself ([pic]).
Hybrid Differencing Scheme
The hybrid differencing scheme utilizes the central differencing scheme for small Peclet numbers while calculating with the upwind differencing schemes for larger Peclet numbers (>2). Since it takes the favourable properties of both differencing schemes, it is conservative, bounded, and transportive. It is stable and used often in CFD codes, with the disadvantage of containing a first-order Taylor series truncation error. It is not available in OpenFOAM.
Power Law Differencing Schemes
The power-law scheme is more accurate than the hybrid. When the Peclet number is greater than 10, diffusion is set to zero, while if the Peclet number is less than 10, a polynomial expression is used to find the flux, resulting in a more exact solution for one-dimensional problems.
Other Differencing Schemes
Since the first-order Taylor series truncation error can be significant and cause false diffusion, higher-order differencing schemes have been created to minimise the error. More neighbouring points are included to widen the influence on the property value of a particular node, while the flow direction is also taken into account. An example is the QUICK (quadratic upwind interpolation for convective kinetics). Another differencing scheme is the TVD (total variation diminishing) scheme, which reduces oscillations in results, which can occur when solving for turbulent energy or rate of dissipation, for example.
2 A2 Pressure drop and friction factor calculations
| | | | | | | |
| | | | | | |
| | | | | |
| | |air frontal U |Reynolds Re |j |f |P drop |
| | | | |
| | |simpleFoam | |
| | |air frontal U |Reynolds Re |
| | |simpleFoam | |
| | |air frontal U |Reynolds Re |P drop |
| | | |simpleFoam | |
| | |air frontal U |Reynolds Re |P drop |f |UminFF1 |UminFF2 |
| | | | | | | |
| | | | | | | |
| | | | | | | |
| | | | | | |
| | | | | |
| | |inlet u |Reynolds Re |j |f |P drop |
| | | | | | | | |
| | | | |density at outflow T =360.77819*T^-1.00336 | |
| | | | |
| | |rhoSmpleFoam | |
| | |P drop |
| | |rhoSimpleFoam |
| | |P drop |
| | | rhoSimpleFoam |
| | |P drop |f |
| | | | |
| | | |
| | |rhoTurbFoam |
| | |P drop |f |UminFF1 |
| | | |
|values for the fin surface effectiveness |are at the different Reynolds numbers. | |
| |
|in the article. | | | | |
| |
|Req/r with Eqn. 23, Φ with Eqn. 22, m with Eqn. 21, η with Eqn. 20and ηowith Eqn. 19. |
| |
| |
|The j-factors are read from the graph in the article. | | | | |
| | |
| | | | | | |
| | | | | | | |
| | | | | | | |
| | | | | | | |
| | | | | | |
| | | | | | |
| | | | | | |
| | | |
| |
|Eqn 8, which requires the heat coefficient h for the air-side. | | | |
| |
|using Eqns 10 and 11, and then h can be found with Eqn 16. Then the j-factor can be | |
|calculated as described above. | | | | | | |
| | | | | | | |
| | |use and from empirical results (Sheet 2) |
| | | | | | | |
| | | | | | | |
| | | | | | | |
| | | | | | |
| | | | | | |
| | | | | | |
| | | |
| | | | | | | | |
| | | | | | | | |
| | | | | | | | |
| | |
| | |
| | | | | | |
| | | | | | |
| | | | | | |
| | |
|Re |u [m/s] |Tout [C] |Tavg [C] |ρ [kg/m3] |Cp [J/kgK] |k [W/mK] |
| | | | | | | | |
| | | | | | | | |
| | | | | | | | |
| | |
|tubeFin.doc |Report with appendix |
|simpleFoam |All OpenFOAM simulation files run with the ’simpleFoam’ named solver and |
|(folder containing the files |named by flow model (laminar, k-epsilon, or k-omega) and inlet airflow velocity |
|listed at right) |laminar |
| |k-epsilon |
| |k-omega |
| | |
| |laminar_0.3 |
| |kEpsilon_0.3 |
| |kOmega_0.3 |
| | |
| |laminar_0.5 |
| |kEpsilon_0.5 |
| |kOmega_0.5 |
| | |
| |laminar_0.7 |
| |kEpsilon_0.7 |
| |kOmega_0.7 |
| | |
| |laminar_1.1 |
| |kEpsilon_1.1 |
| |kOmega_1.1 |
| | |
| |laminar_1.5 |
| |kEpsilon_1.5 |
| |kOmega_1.5 |
| | |
| |laminar_2.5 |
| |kEpsilon_2.5 |
| |kOmega_2.5 |
| | |
| |laminar_2.5 |
| |kEpsilon_3.7 |
| |kOmega_3.7 |
| | |
| |laminar_2.5 |
| |kEpsilon_4.5 |
| |kOmega_4.5 |
| | |
| |laminar_2.5 |
| |kEpsilon_5.4 |
| |kOmega_5.4 |
| | |
| |laminar_2.5 |
| |kEpsilon_6.2 |
| |kOmega_6.2 |
| | |
|rhoSimpleFoam |All OpenFOAM simulation files run with the ’rhoSimpleFoam’ solver and |
|(folder containing the files |named by flow model (laminar, k-epsilon, or k-omega) and inlet airflow velocity |
|listed at right) |laminar |
| |k-epsilon |
| |k-omega |
| | |
| |laminar_0.3 |
| |kEpsilon_0.3 |
| |kOmega_0.3 |
| | |
| |laminar_0.5 |
| |kEpsilon_0.5 |
| |kOmega_0.5 |
| | |
| |laminar_0.7 |
| |kEpsilon_0.7 |
| |kOmega_0.7 |
| | |
| |laminar_1.1 |
| |kEpsilon_1.1 |
| |kOmega_1.1 |
| | |
| |laminar_1.5 |
| |kEpsilon_1.5 |
| |kOmega_1.5 |
| | |
| |laminar_2.5 |
| |kEpsilon_2.5 |
| |kOmega_2.5 |
| | |
| |laminar_2.5 |
| |kEpsilon_3.7 |
| |kOmega_3.7 |
| | |
| |laminar_2.5 |
| |kEpsilon_4.5 |
| |kOmega_4.5 |
| | |
| |laminar_2.5 |
| |kEpsilon_5.4 |
| |kOmega_5.4 |
| | |
| |laminar_2.5 |
| |kEpsilon_6.2 |
| |kOmega_6.2 |
| | |
|rhoTurbFoam |OpenFOAM simulation file run with ‘rhoTurbFoam’ solver and kOmegaSST turbulence model with inlet airflow |
| |velocity 6.2 m/s |
|pDrop_f__j.xls |Excel file with calculations for determining the friction factor f (Appendix A2) |
|Heat_tx_calcs.xls |Excel file with calculations for determining Colburn j-factor (Appendix A3, A4, and A5) |
|transPressure.xls |Excel file for the transient pressure trace |
|Samples |File folder containing the text files (.xy) for the pressure and velocity meaurements recorded from inflow and |
| |outflow during the simulations. Arranged by file type as in the above files. |
-----------------------
Hybrid recorder
[pic]
hot water
thermostat
reservoir
x
Free stream plane
(slip BC)
fin
(wall BC)
Tube
(wall BC)
Outlet
(Outlet BC)
Dc=10.23 mm
N
[pic]
Figure 30. CFD Fluid element for calculating changes in fluid property, with two coordinate systems in three dimensions x, y, z, and according to faces N, S, E, W, T and B.
[Hj ertager, 2007] [Versteeg and Malalasekera, 2007]
Free stream plane
(slip BC)
Side plane
(symmetric plane BC)
Inlet
(inlet BC)
Three-dimensional CFD simulations are carried out to investigate heat transfer and fluid flow characteristics of a two-row plain fin-and-tube heat exchanger using Open FOAM, an open-source CFD code. Heat transfer and pressure drop characteristics of the heat exchanger are investigated for Reynolds numbers ranging from 330 to 7000. Model geometry is created, meshed, calculated, and post-processed using open source software. Fluid flow and heat transfer are simulated and results compared using both laminar and turbulent flow models (k-epsilon, and Menter SST k-omega), with steady-state solvers to calculate pressure drop, flow, and temperature fields. Model validation is carried out by comparing the simulated case friction factor f and Colburn factor j to experimental results from the literature. For friction factor determination, little difference is found between the flow models simulating laminar flow, while in transitional flow, the laminar flow model produced the most accurate results and the k-omega SST turbulence model was more accurate in turbulent flow regimes. The most accurate simulations for heat transfer in laminar flow are found using the laminar flow model, while heat transfer in transitional flow is best represented with the SST k-omega turbulence model, and heat transfer in turbulent flow is more accurately simulated with the k-epsilon turbulence model. Reasonable agreement is found between the simulations and experimental data, and the open-source software has been sufficient for simulating the flow fields in tube-fin heat exchangers.
[pic]
[pic] [pic]
Figure 1. Vestas Aircoil A/S heat exchanger and close-up of fin-and-tube arrangement.
W
[pic]
Figure 9: CFD computational domain for heat exchanger simulation in this project.
[pic]
Figure 3: Typical fin-and-tube heat exchanger section with
staggered tube arrangement. [Song and Nishino, 2008]
0.00597 0.881 1.76 2.68 3.50
[pic]
[pic]
Figure 18. Contours of turbulent kinetic energy k distribution, SST k-omega model, inlet air velocity 6.2 m/s.
278 292 306 319 333
[pic]
[pic]
Figure 21. Contours of temperature Field SST k-omega flow model, z-direction, 6.2 m/s inlet air velocity.
Equation 7
[pic]
-for tube-and-fin heat exchangers
y
PC
278 292 306 319 333
[pic]
[pic]
Figure 20. Contours of temperature field, SST k-omega flow model, 6.2 m/s inlet air velocity.
1
B
[pic]
[pic]
Figure 22. Comparison of Figures 15 and 19 to illustrate synergy between flow direction and temperature streamlines. (Flow velocity streamlines shown in top picture, and isothermal pattern shown in the lower picture).
[pic]
entrance effect
Flow acceleration
due to density
differences
core
friction
exit effect
z
(I) (II) (III) (IV)
or
[pic]
[pic]
Figure 27. Change of mass in fluid element.
[Hjertager, 2007]
sum of all outflows
Sum of all inflows
mass accumulation
over time
[pic]
[pic]
(x, y, z)
S
E
B
T
W
N
Figure 4. Illustration of the main computational domain and geometric parameters of the heat exchanger model studied (z-direction not shown).
278 292 306 319 333
[pic]
[pic]
Figure 19. Contours of temperature field, SST k-omega flow model, 0.3 m/s inlet air velocity.
*
*
*
*
*
*
[pic]
Figure 26. CFD Fluid element for calculating changes in fluid property, coordinate systems in three dimensions x, y, z, and according to faces N, S, E, W, T and B. [Hjertager, 2007] [Versteeg and Malalasekera, 2007]
[pic]
(I) (II) (III) (IV) (V)
(I) (II) (III) (IV) (V)
(I) (II) (III) (IV) (V)
(I) (II) (III) (IV) (V) (VI)
where [pic], [pic], [pic], and [pic] are constants.)
[pic]
START
([pic]) ([pic]) ([pic]) ([pic])
STOP
Initial guesses for p*, velocity components u*, v*, w*, and other scalar properties φ* (i.e. T).
STEP 1: Solve discretised momentum equations
STEP 2: Solve for pressure correction equation
STEP 3: Correct pressure and velocities
STEP 4: Solve the other discretised transport equations
Convergence?
Yes
No
Set solved values equal to new initial guesses.
p*=p, u*=u, v*=v, w*=w, φ= φ*
Figure 12. The SIMPLE Algorithm [Versteeg and Malalasekera, 2007]
p’
p, u, v, w, φ*
φ
u*, v*, w*
[Hjertager, 2007]
44 mm (2 x Pl)
11 mm (Pl/2)
Pl = 22 mm
12.7 mm (Pt/2)
2
17
3
4
5
6
7
8
9
10
11
12
13
14
15
16
[pic]
Figure 8. CFD computational domain for heat exchanger illustrating partitions.
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
[pic]
time (s)
Pressure [Pa]
0.00 0.221 0.442 0.668 0.884
[pic]
[pic]
Figure 15. Contours and vector plot for U velocity field, SST k-omega flow model, inlet air flow 0.3 m/s.
T
0.00 3.75 7.50 11.2 15.0
[pic]
[pic]
Figure 16. Contours and vector plot for U velocity field, SST k-omega flow model, inlet air flow 6.2 m/s.
8.25e-05 0.00109 0.00217 0.00325 0.00433
[pic]
[pic]
Figure 17. Contours of turbulent kinetic energy k distribution, SST k-omega model, inlet air velocity 0.3 m/s.
E
S
(x, y, z)
δxe
δxw
BC=constant
xw
W
Δx
xe
P
(x, y, z)
E
BC=constant
................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related download
- numerical simulation of plain fin and tube heat exchangers
- theory and applications of boundary element
- california state university northridge
- free surface cfd prediction of components of ship
- introduction of computational fluid dynamics
- paper no 200000
- numerical analysis of auto ignition of ethanol
- kursus mekanik bendalir
- write here the title of your paper researchgate
Related searches
- list of autoimmune diseases and symptoms
- types of home loans and qualifications
- list of all stocks and prices
- list of autoimmune diseases and sympto
- office of state lands and investments
- university of scranton tuition and fees
- department of labor wage and hour
- side effects of baking soda and water
- latent and sensible heat flux
- solidworks simulation of assembly
- neck and shoulder heat massager
- sensible and latent heat fluxes