Numerical methods in heat transfer and fluid dynamics

MSc Energy Enginering Numerical methods in heat transfer and fluid dynamics MASTER THESIS Author: Sergio Zavaleta Camacho Director: Carlos David Pérez Segarra Announcement: September 2017 Escola Tècnica Superior d’Enginyeria Industrial de Barcelona Numerical methods in heat transfer and fluid dynamics Page 1 Summary Numerical methods in fluid dynamics and heat transfer are experiencing a remarkable growth in terms of the number of both courses offered at universities and active researches in the field. There are some software packages available that solve fluid flow problems. Nevertheless, Computational Fluid Dynamics (CFD) codes are progressively being accepted as design tools by the industry. Nowadays users of CFD need to be fairly knowledgeable, which requires instruction of both students and working engineers. The present text is a starting point to immerse the student in the tissues of the field. The two main objectives of this project are: to acquire a basic training in the numerical resolution of the governing equations in the heat transfer and fluid dynamics, and to get used to CFD and Heat Transfer (HT) codes and acquire the skills to critically judge their quality, this is, apply code verification techniques, validation of the used mathematical formulations and verification of numerical solutions. In the present text, fundamental methods for solving the transport phenomena are covered. Chapter 1. ‘Discretization and solvers’ contains the fundamental numerical method since the physical phenomena must be described through appropriate differential equations. Chapter 2. ‘Heat conduction methods’ is the construction base of the numerical method, therefore emphasis on concepts and calculation details are given here. Chapter 3. ‘Analysis of the general convection-diffusion equation’ is focused on the interaction of convection and diffusion, with the flow field known in advance. Finally, the calculation of the velocity field itself is treated in Chapter 4. ‘Incompressible flow method using the Navier-Stokes equations’. This chapter represents an effort to employ the Fractional Step Method (FSM) in the solution of the Navier-Stokes equations with the aim to obtain solutions for diverse Reynolds numbers and mesh refinements. The problems presented and solved are intended to be a material base over which analysis, discussion and conclusions are developed. The Smith-Hutton problem is addressed since many of the features commonly encountered in practical convection-diffusion problems are here present. Different numerical schemes are submitted and their pros and cons are described. Moreover, the robustness of the Fractional-Step Method (FSM) has been demonstrated using the Driven cavity flow benchmark problem. Detailed accurate results have been presented for this problem. Up to 128x128 computational points and Reynolds as high as 3200 have been considered. Keywords – numerical methods, fluid dynamics, heat and mass transfer, convectiondiffusion, convective schemes, Smith-Hutton, incompressible flow, Navier-Stokes, fractionalstep method, staggered meshes, Driven cavity flow. Page 2 Master Thesis Index SUMMARY ___________________________________________________ 1 INDEX _______________________________________________________ 2 GLOSSARY __________________________________________________ 5 PREFACE ____________________________________________________ 6 INTRODUCTION _______________________________________________ 7 Objectives ................................................................................................................ 8 Scope ...................................................................................................................... 8 1. DISCRETIZATION AND SOLVERS ___________________________ 10 1.1. Finite volume method ................................................................................... 10 1.1.1. Approximation to surface integrals ................................................................... 11 1.1.2. Approximation to volume integrals ................................................................... 12 1.1.3. Implementation of boundary conditions ............................................................ 12 1.1.4. The algebraic equation system ........................................................................ 13 1.2. Solution of Linear Equation Systems ........................................................... 13 1.2.1. Gauss-Seidel method ...................................................................................... 13 1.2.2. Tri-diagonal matrix algorithm ............................................................................ 14 1.2.3. Line-by-line method ......................................................................................... 15 1.2.4. The relaxation factor ........................................................................................ 16 2. HEAT CONDUCTION METHODS ____________________________ 18 2.1. Heat conduction discretized equation .......................................................... 18 2.2. The Four materials problem ......................................................................... 22 2.2.1. Problem definition ............................................................................................ 22 2.2.2. Code development ........................................................................................... 23 2.2.3. Algorithm .......................................................................................................... 28 2.2.4. Verification ....................................................................................................... 29 2.2.5. Validation ......................................................................................................... 30 2.2.6. Results ............................................................................................................. 33 3. ANALYSIS OF THE GENERAL CONVECTION-DIFFUSION EQUATION ______________________________________________ 37 3.1. The convection-diffusion equation ............................................................... 37 3.1.1. Domain and temporal discretization ................................................................. 38 3.1.2. Mass conservation equation ............................................................................ 38 Numerical methods in heat transfer and fluid dynamics Page 3 3.1.3. Discretization of the convection-diffusion equation .......................................... 39 3.1.4. Evaluation of the convective terms .................................................................. 40 3.1.5. Normalized variables ....................................................................................... 42 3.1.6. Final form of the discrete convection-diffusion equation .................................. 43 3.1.7. Boundary conditions........................................................................................ 44 3.2. The Smith-Hutton problem ........................................................................... 45 3.2.1. Problem definition ........................................................................................... 45 3.2.2. Code development .......................................................................................... 46 3.2.3. Algorithm ......................................................................................................... 51 3.2.4. Verification ...................................................................................................... 52 3.2.5. Results ............................................................................................................ 55 3.2.6. Discussion ....................................................................................................... 57 4. INCOMPRESSIBLE FLOW METHOD USING THE NAVIER-STOKES EQUATIONS _____________________________________________ 65 4.1. Introduction to the Fractional Step Method ................................................... 65 4.1.1. Application of the Helmholtz-Hodge decomposition theorem .......................... 66 4.1.2. Time integration method ................................................................................. 67 4.1.3. Determination of ......................................................................................... 69 4.1.4. The Poisson equation ..................................................................................... 69 4.2. The Driven cavity flow problem .................................................................... 72 4.2.1. Problem definition ........................................................................................... 72 4.2.2. Code development .......................................................................................... 73 4.2.3. Algorithm ......................................................................................................... 80 4.2.4. Verification ...................................................................................................... 82 4.2.5. Results ............................................................................................................ 95 4.2.6. Discussion ....................................................................................................... 99 CONCLUSIONS _____________________________________________ 101 FUTURE ACTIONS ___________________________________________ 103 BIBLIOGRAPHY _____________________________________________ 104 APPENDIX _________________________________________________ 106 A. Previous calculations .................................................................................. 106 B. Forced convection inside an isothermal tube ............................................. 112 Numerical methods in heat transfer and fluid dynamics Page 5 Glossary 1D one-dimensional 2D two-dimensional 3D three-dimensional CDS central difference scheme CFD computational fluid dynamics CV control volume CTTC Heat and Mass Transfer Technological Center EDS exponential difference scheme FD finite difference FE finite elements FSM fractional-step method FV finite volume GS Gauss-Seidel method HT heat transfer NS Navier-Stokes QUICK quadratic interpolation for convective kinematics SMART sharp and monotonic algorithm for realistic transport SUDS second-order upwind difference scheme TDMA tri-diagonal matrix algorithm UDS upwind difference scheme Pág. 6 Memoria Preface This project is worked out as a result of my attendance to the courses Heat and Mass Transfer, Computational Engineering and Numerical Methods in Fluid Dynamics and Heat Transfer given by Professor Carlos-David Pérez-Segarra in the Tarrasa School of Industrial, Aerospace and Audiovisual Engineering (ESEIAAT), and the development of the exercises proposed by the professor. Numerical methods in fluid dynamics and heat transfer are experiencing a remarkable growth in terms of the number of both courses offered at universities and active researches in the field. There are some software packages available that solve fluid flow problems. Nevertheless, Computational Fluid Dynamics (CFD) codes are progressively being accepted as design tools by the industry. Nowadays users of CFD need to be fairly knowledgeable, which requires instruction of both students and working engineers. The present text is a starting point to immerse the student in the tissues of the field. It is designed to deploy the topic in a certain sequence, and the reader is preferably to follow the same sequence. It will not be profitable to jump to a later chapter, as all chapters build upon the material covered in the previous ones. Whilst I have invested a lot of care to avoid typing or spelling mistakes, no doubt some remain to be found by readers. I will appreciate your notifying me of any mistakes you might come across. For that purpose, the author’s electronic mail address is given: [email protected] I wish to record my sincere thanks to Professor C.D. Pérez-Segarra for inviting me to his lessons and his incalculable value in my learning of numerical methods in fluid dynamics and heat transfer, for his tutorship and support since the basic heat transfer course until the end of my Master Thesis. Also, my gratitude towards the Heat and Mass Transfer Technological Center (CTTC-UPC) staff for the study material provided and their quick attendance at all times to my questions and doubts. I am very grateful to my parents Victor and Coti for their affection, to my sister Gabriela for the spelling and grammatical corrections, to my aunt Patricia, to Laurent, Marco, Mariuxi and Garbiñe for their solidarity with me throughout my Master’s studies. S. Z. Barcelona, Spain September 2017 Numerical methods in heat transfer and fluid dynamics Page 7 Introduction Heat transfer (HT) phenomenon plays an important role in many industrial and environmental problems. First and foremost, in the applications of energy production and transformation, there is not a single application in this area that does not involve HT effects in some way or another. It is widely involved in the generation of power from conventional fossil fuels, nuclear sources or the use of geothermal energy sources. Heat transfer processes determine the design of industrial equipment such as boilers, condensers and turbines. Quite often the challenge is to maximize heat transfer rates such as in heat exchangers, or to minimize it like in insulations. Likewise, HT arises in the design of solar energy conversion systems for water and space heating, cooling of electronic equipment, or refrigeration and air conditioning systems. HT issues also occur in air and water pollution problems and strongly influence climate at local and global scale. CFD involves the analysis of fluid flow and related phenomena using numerical solution methods. Due to the rapid advances in CFD and the potential it provides to analyze, on a fundamental basis, systems of considerable interest to the energy engineer. Nowadays users of CFD need to be fairly knowledgeable, which requires instruction of both students and working engineers. It can be anticipated that the importance of CFD as a “workhorse” for the energy engineering community will rapidly increase soon. In the present text, fundamental methods for solving the transport phenomena are covered. Chapter 1. ‘Discretization and solvers’ contains the fundamental numerical method since the physical phenomena must be described through appropriate differential equations. Chapter 2. ‘Heat conduction methods’ is the construction base of the numerical method, therefore emphasis on concepts and calculation details are given here. Chapter 3. ‘Analysis of the general convection-diffusion equation’ is focused on the interaction of convection and diffusion, with the flow field known in advance. Finally, the calculation of the velocity field itself is treated in Chapter 4. ‘Incompressible flow method using the Navier-Stokes equations’. This chapter represents an effort to employ the Fractional Step Method (FSM) in the solution of the Navier-Stokes equations (NSE) with the aim to obtain solutions for diverse Reynolds numbers and mesh refinements. The problems presented and solved are intended to be a material base over which analysis, discussion and conclusions are developed. The Smith-Hutton problem is addressed since many of the features commonly encountered in practical convection-diffusion problems are here present. Different numerical schemes are submitted and their pros and cons are described. Moreover, the robustness of the Fractional-Step Method (FSM) has been demonstrated using the Driven cavity flow benchmark problem. Detailed accurate results Page 8 Master Thesis have been presented for this problem. Up to 128x128 computational points and Reynolds as high as 3200 have been considered. Objectives The objectives of the project are:  Acquire a basic training in the numerical resolution of the governing equations in the fluid dynamics and heat and mass transfer.  Acquire practical experience in the programming, validation and verification of CFD and HT codes.  Get familiarized with the use of CFD & HT codes and acquire the ability to critically judge their quality, this is, apply code verification techniques, validation of the used mathematical formulations and verification of numerical solutions.  Consolidation of the mathematical formulations of phenomena of fluid dynamics and heat and mass transfer.  Knowledge of numerical integration methodology of the NSE. Scope  Discretization and solvers General approach of the problems involved in the integration of the equations of fluid dynamics and heat and mass transfer. The finite differences methods, the finite volume method, and the standard iterative solvers are presented. The finite volume method is favored in this text.  Conduction heat transfer Development of the methodology explained in basic courses of heat and mass transfer, based on finite volume techniques and structured, orthogonal and domain adaptable discretization meshes. At this stage, the resolution of the systems of the discretization equations is done with direct and iterative methods i.e. Gauss-Seidel method (GS), Tri- Diagonal Matrix Algorithm (TDMA), line by line, and relaxation factors. Moreover, a two-dimensional (2D) conduction HT benchmark problem in transient regime is presented, solved, validated and verified. Numerical methods in heat transfer and fluid dynamics Page 15 (11) Which can be rearranged to look like Eq (9). In other words, the coefficients and then stand for (12) (13) These are recurrence relations, since they give and in terms of and . To start the recurrence process, we note that Eq. (7) for is almost of the form Eq. (9). Thus, the values of and are given by and (14) At the other end of the , sequence, we note that . This leads to , and hence from Eq. (9) it is obtained (15) Now we start the back substituion via Eq. (9). The algorithm results 1. Calculate and from Eq. (14) 2. Use the recurrence relations (13) to obtain and for i = 2,3, …, N. 3. Set 4. Use Eq. (9) for , to obtain , , …, , , . 1.2.3. Line-by-line method A suitable combination of the TDMA for one-dimensional situations and the GS can now be formed. We shall choose a grid line (say, in the y direction), assume that the temperatures (or the property of interest) along the neighboring lines (i.e., the x- and z-direction neighbors of the points on the chosen line) are known from their most-recent values, and solve for the temperatures (or the property of interest) along the chosen line by the TDMA. We shall follow this procedure for all the lines in one direction and repeat the procedure, if desired, for the Page 16 Master Thesis lines in the other direction(s). This method is equally applicable to 2D or 3D. Figure 4 shows the line-by-line scheme. Figure 4. Representation of the line-by-line method. [Patankar, 1980] The discretization equations for the grid points along a chosen line are considered. They contain the temperatures at the grid point (shown by crosses) along the two neighboring lines. If these temperatures are substituted from their latest values, the equations for the grid points (shown by dots) along the chosen line would look like one-dimensional equations and could be solved by the TDMA. This procedure is carried out for all the lines in the y direction and may be followed by a similar treatment for the x direction. The convergence of the line-by-line method is faster, because the boundary-condition information from the ends of the line is transmitted at once to the interior of the domain, no matter how many grid points lie along the line [Patankar, 1980]. The rate of transmission of information in the other direction is like that of the GS. Figure 5. Boundary conditions that make a left-to-right sweep more adventageous. [Patankar, 1980] 1.2.4. The relaxation factor In the iterative solution of the algebraic equations or in the overall iterative scheme employed Numerical methods in heat transfer and fluid dynamics Page 17 for handling nonlinearity, it is often desirable to speed up or to slow down the changes, from iteration to iteration, in the values of the dependent variable. This process is called overrelaxation or underrelaxation. With the line-by-line method, underrelaxation is a very useful device for nonlinear problems. It is often employed to avoid divergence in the iterative solution of strongly nonlinear equations. There are many ways of introducing overrelaxation or underrelaxation. In this text, the preferred way is (16) When the relaxation factor is between 0 and 1, its effect is underrelaxation; that is, the values of stay close to . For a very small value of , the change in becomes very slow. When f is greater than , overrelaxation is produced. Page 18 Master Thesis 2. Heat conduction methods In this chapter, we shall begin the task of constructing a numerical method for solving the general differential equation (17), which governs the physical process of interest here. (17) This equation contains four basic terms. Here we shall omit the convection term (the second term) and concentrate on the remaining three terms. The construction of the method will be completed in Chapter 3, where the treatment of the convection term will be discussed. Omission of the convection term reduces the situation to a conduction-type problem. Heat conduction provides a convenient starting point for our formulation, because the physical processes are easy to understand and the mathematical complication is minimal. The objectives of this chapter, however, go far beyond presenting a numerical method for heat conduction alone. First, other physical processes are governed by very similar mathematical equations. Among these are potential flow, mass diffusion, flow through porous media, and some fully developed duct flows. The numerical techniques described in this chapter are directly applicable to all these processes. Electromagnetic field theory, diffusion models of thermal radiation, and lubrication flows are further examples of phenomena governed by conduction-type equations. Second, this chapter accomplishes much of the preparatory work needed for later chapters. The procedure for the solution of the algebraic equations is presented here in a once-and- for-all manner. Later chapters modify the content of the algebraic equations, but the same solution technique continues to be applicable. 2.1. Heat conduction discretized equation A portion of a 2D grid is shown in Figure 6. For the grid point P (denoting principal), points E and W (denoting east and west) are its x-direction neighbors, while N and S (denoting north and south) are the y-direction neighbors. The CV around P is shown dashed by lines. Its thickness in the z-direction is assumed to be unity. The HT surfaces are indicated with a lowercase letter. Numerical methods in heat transfer and fluid dynamics Page 19 Figure 6. Control volume for the 2D situation. [Patankar, 1980] The differential heat conduction equation is, (18) this equation is derived from the energy conservation equation (19). Figure 7 shows the energy balance flows of a CV. Figure 7. Energy balance of an internal node In general, the energy balance of a 2D CV is given by: (19) Where, Page 20 Master Thesis : Conduction heat flux entering through the west face : Conduction heat flux leaving through the east face : Conduction heat flux entering through the south face : Conduction heat flux leaving through the north face : Internal heat flux generation in the CV due to external sources : Accumulated heat flux in the CV over the time-step Depending on the scheme employed -related to the selection of the parameter - the equation can take three different forms: (20) In the implicit scheme, the heat fluxes are evaluated at the instant , in the explicit scheme, at the instant ; meanwhile, in the Cranck-Nicholson scheme the heat fluxes are averaged over both instants. Each scheme has its advantages and disadvantages, but that discussion is beyond the scope of this text. An implicit scheme is selected to ensure the stability. For simplicity, the superindex will not be explicitly shown in the equations. Then the previous equation using an implicit scheme turns into (21) The accumulated energy in the control volume is approximated to (22) Using the Fourier’s law approximated to a 1st order Taylor series; and using Eq. (22), Eq. (22) turns into (23) : Temperature of the west node; : Temperature of the east node Numerical methods in heat transfer and fluid dynamics Page 21 : Temperature of the south node; : Temperature of the north node : Internal volumetric heat generation flux As the specific capacity in all the nodes is not time-dependent the result of the time-average is the same value (24) Eq. (23) can be turned into the general 2D discretization equation (see Eq. (5)) (25) Adopting the solver form, the discretization coefficients of the internal nodes are (26) (27) (28) (29) (30) (31) Page 22 Master Thesis 2.2. The Four materials problem The Four materials problem is a benchmark heat conduction problem in transient regime problem proposed by the CTTC (2016). The name ‘four materials’ is because the problem is modeled as a rod composed by four different materials, each one with its own thermophysical properties and size. The problem is defined, solved, verified and validated next. 2.2.1. Problem definition A very long rod is composed of four different materials (M1 to M4), represented with different colors in Figure 8. All the lines are parallel to the coordinate axis. The coordinates of the points p1 to p3 are given in Table 1. The properties of the materials are given in Table 2. Each of the four sides of the rod interacts with the surrounding in a different manner, as described in Table 3. The initial temperature field is It is asked to assess the transient response of single nodes of each material. Figure 8. Scheme of the conduction problem2. [CTTC, 2016] 2 The drawing is not scaled to the true dimensions of the materials. Numerical methods in heat transfer and fluid dynamics Page 23 Table 1. Problem coordinates x [m] y [m] p1 0.50 0.40 p2 0.50 0.70 p3 1.10 0.80 Table 2. Thermophysical properties [kg/m3] [J/kgK] [W/mK] M1 1500 750 170 M2 1600 770 140 M3 1900 810 200 M4 2500 930 140 Table 3. Boundary conditions Cavity wall Boundary condition Bottom Isotherm at Top Uniform length Left In contact with a fluid at and heat transfer coefficient Right Uniform temperature (where is the time in seconds) 2.2.2. Code development As this is a 2D heat conduction HT problem in transient regime, the domain is discretized in the x-direction and y-direction, obtaining a finite number of control volumes and nodes. The map temperature of all the nodes for a given instant ‘ ’ can be calculated using the the solver for a given time instant. Hence energy balances shall be applied over all the nodes of the domain to find their discretization coefficients. In this way, the calculated coefficients are the input parameters of the solver. The energy balances consider the boundary conditions and input data -geometric data, physical data, numerical inputs-, previous geometrical calculations and the map temperature of the previous instant ‘ ’. To calculate the temperature map of the next time instant, the time instant is updated and the process is repeated. The detailed solution of the problem is developed next. Page 24 Master Thesis Assumptions  Two-dimensional flow – x-axis and y-axis–, z-axis flow is neglected.  Constant thermophysical properties –density, thermal conductivity, specific capacity– for each material.  HT phenomena: only conduction within the solid, conduction and convection in the left-side and top-side.  Cell-centered discretization: average values over each cell –temperature, density, specific capacity, thermal conductivity–. Domain discretization The nodes are placed in the center of each control volume. Figure 9 shows the control volume discretization represented by small rectangles. The nodes are supposed to be placed in the center of each CV -not shown-. Figure 9. Discretization of the wall Previous calculations, i.e. geometric relationships and harmonic mean, are solved in the Appendix. Numerical methods in heat transfer and fluid dynamics Page 31 (38) Convergence criterion sensitivity Table 4 shows the variation of and the relative error regarding the convergence criterion . For the results are equal, this mean those values are adequate for the simulation. Table 4. Convergence criterion sensitivity Heat flux [W] Relative error e-9 72.43 0.00% e-8 72.43 Base case e-7 72.43 0.00% e-6 72.44 0.00% e-5 72.46 0.03% e-4 72.59 0.21% e-3 73.53 1.51% e-2 83.94 15.89% e-1 84.35 16.46% Time-step sensitivity To assess the time-step sensitivity a convergence criterion and a node density are used. Table 5 shows the variation of regarding the timestep, , it shows the heat flux is stable and the relative error is low even at relatively large time-steps such as 50 seconds or 100 seconds. This is due to the implicit scheme employed. Page 32 Master Thesis Table 5. Time-step sensitivity ∆t Heat flux [W] Relative error 1 73.05 Base case 2 73.05 0.00% 5 73.06 0.01% 10 73.07 0.03% 20 73.09 0.06% 25 73.10 0.07% 50 73.15 0.14% 100 73.26 0.29% 200 73.47 0.58% 250 73.58 0.72% 500 74.12 1.47% Grid sensitivity To assess the grid sensitivity a convergence criterion and a time-step of =1 are used. Figure 14 shows the variation of regarding the node density, . Figure 14. Grid sensitivity Table 6 shows the calculated heat flux values and the related relative error regarding and the related error regarding There is a relative error of 1.71% if the node density is nodes/m. More accurate results are obtained using is nodes/m, although the drawback is that the simulation lasts longer. Then a compromise between number of nodes and relative fast time response would be appropriate e.g. (relative error < 1%). 70 71 72 73 74 75 0 20 40 60 80 100 120 140 Left-side heat flux [W] Node density: n [nodes/m] Numerical methods in heat transfer and fluid dynamics Page 33 Table 6. Grid sensitivity Node density Heat flux [W] Relative error Relative error 120 73.70 1.74% Base case 100 73.51 1.49% 0.25% 80 73.24 1.11% 0.62% 75 73.15 0.98% 0.75% 70 73.05 0.85% 0.88% 65 72.92 0.67% 1.05% 60 72.79 0.49% 1.23% 55 72.62 0.26% 1.46% 50 72.44 Base case 1.71% 45 72.19 0.34% 2.05% 40 71.90 0.74% 2.44% 35 71.50 1.29% 2.98% In Table 7 is proposed a set of numerical parameters based on the previous numerical data sensitivity studies to have accurate results and fast computing response. Table 7. Optimized numerical inputs δ ∆t (s) n (nodes/m) e-4 200 70 The resulting heat flux has a relative error of 0.85% regarding one of the most accurate heat fluxes calculated ( ) and the computational response is around 2 minutes (in a regular computer) for an end time of 10000 s. For further validation, the Crank-Nicholson scheme and the explicit scheme can be employed in the same way to compare the results with the ones of the implicit scheme. The algorithm would remain the same; the only change to consider is modifying the parameter in Equation 1 for the energy balance of the internal nodes, and analogously for the energy balance equations of the left-side nodes and top-side nodes. 2.2.6. Results Figure 15 shows the evolution of four nodes of the wall from until . These four points are the gravity center of each one of the materials. For practicality, each node Page 34 Master Thesis can be thought as a representative temperature of each material8. The temperatures’ evolution can be thought as a superposition of a transient regime from until reaching the steady-state given by left-side, bottom-side and top-side boundary conditions, plus a constant perturbation due to the right-side boundary condition. Thermal inertia is defined as . Material 1 and Material 2 have the lower thermal inertia; this explains why they would last less to reach the fictional steady-state, while Material 3 and Material 4 would last longer. Figure 15. Temporal evolution of the representative nodes’ temperature. Figure 16 shows the temperature maps for the instants , and Figure 17 shows the temperature maps for the instants .9 In Figure 16 the temperature contours significantly change between , and . The wall is heated from the top-side due to and from the right-side due to the temperature increase of the right-side nodes as a function of the time. We may think that the whole wall is going through a transition, from an initial map of 8°C until it reaches its steady state if the right-side boundary condition would not exist. In Figure 17 the temperature contours do not present significant change, but the temperature from the right half of the wall increases, being the maximum temperature 32°C at and 55°C at . This behavior can be thought as coexistence 8 Although this is not technically correct, it is used only for interpretation purposes. 9 Note that the temperature scale changes from map to map and a specific color corresponding to a specific temperature in one map may correspond to a different temperature other maps. 02000 4000 6000 8000 10000 12000 5 10 15 20 25 30 35 40 time (s) Temperature (°C) Material 1 Material 2 Material 3 Material 4 Numerical methods in heat transfer and fluid dynamics Page 35 of a steady regime and a transient regime. The steady state is ‘reached’ by the left half of the wall10, which ‘reaches’ a temperature around 25°C since the left, bottom and top boundary conditions are not time-dependent, while the transient regime occurs due to the right boundary condition and influences the temperature of the right half of the wall. It was verified that the convergence is reached slower with the GS than with the line-by- line method, especially when many grid points are involved. The reason for the slowness is easy to understand, the GS transmits the boundary-condition information at a rate of one grid interval per iteration. Figure 16. Temperature maps at 1000 s, 2000 s, 3000 s and 4000s. 10 This shall be taken as an ‘illustrative explanation” although it is not strictly correct because the rightside boundary condition has an effect over the left half of the wall. 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 1000s (°C) 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 2000s(°C) 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 3000s(°C) 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 4000s(°C) 12 14 16 18 20 22 16 17 18 19 20 21 22 23 19 19.5 20 20.5 21 21.5 22 22.5 23 21 22 23 24 25 26 27 28 Page 36 Master Thesis Figure 17. Temperature maps at 5000 s, 6000 s, 8000 s and 10000 s. 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 5000s(°C) 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 6000s(°C) 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 8000s(°C) 10 20 30 40 50 5 10 15 20 25 30 35 40 Nodes [j] Nodes [i] Temperature map 10000s(°C) 25 30 35 40 45 50 55 24 26 28 30 32 24 26 28 30 32 34 36 38 25 30 35 40 45 Numerical methods in heat transfer and fluid dynamics Page 37 3. Analysis of the general convection-diffusion equation This chapter shows how the FVM is applied to a model of convective transport: the 2D convection-diffusion equation. There are two main goals of this chapter; the first is to expose the FVM, and the second goal is to introduce and compare the numerical schemes meddling the convective term in transport equations. The 2D convection-diffusion equation is a compact, though somewhat non-physical, model of transport of heat, mass and other passive scalars. Applying the FVM to this equation allows different schemes for approximating the convection term to be compared. This chapter should be considered a brief introduction to the topic of convection modeling schemes. The schemes considered are the upwind difference scheme (UDS), the central difference scheme (CDS), the exponential difference scheme (EDS), the second-order upwind scheme (SUDS) and the quadratic interpolation for convective kinematics (QUICK). 3.1. The convection-diffusion equation The general differential equation is (39) where is the diffusion coefficient and the source term. The quantities and are specific to a particular meaning of . The four terms in the general differential equation are the unsteady term, the convection term, the diffusion term, and the source term. The mass conservation equation can also be written as Eq. (39) where , , . Table 8 shows that momentum, energy and species conservation equations can also be written using the convection-diffusion equation. Page 38 Master Thesis Table 8. Conservation equations written with the convection-diffusion form . 3.1.1. Domain and temporal discretization The FVM approach is employed to discretize the general convection-diffusion equation. Structured meshes (centered nodes or centered faces) and unstructured meshes (usually centered nodes) are employed to model the spatial domain [CTTC, 2017]. Time is discretized (see Figure 18) and the time-step has an essential influence over the stability of the solution. Figure 18. Scheme of the time discretization 3.1.2. Mass conservation equation The mass conservation equation on its discretized form shall be obtained as a previous step for discretizing the convection-diffusion equation, in its differential form is . . . t t=0 t1t2t3. . . tntn+1 Numerical methods in heat transfer and fluid dynamics Page 39 (40) Integrating over the volume and time (41) For an easier notation: and . Assuming second-order approximations, implicit scheme , mass flow at the faces: in the positive coordinate direction . (42) 3.1.3. Discretization of the convection-diffusion equation A 2D discretization is carried out but it is easily extended to fully 3D. The Eq. (39) using the methodology explained in 3.1.2. results (43) where , , , Equation (43) is re-written using the mass conservation equation and introducing a linearization of the source term, . . (44) Page 40 Master Thesis 3.1.4. Evaluation of the convective terms Convective terms can adapt different approximation profiles, which are defined next.  CDS This scheme generates a linear interpolation assuming that the cell face is placed between the nodes. (45) Figure 19. CDS scheme  UDS The UDS assumes that the convected variable at the cell face is the same as the upwind cell-center value. (46) Figure 20. UDS scheme Numerical methods in heat transfer and fluid dynamics Page 47 Domain discretization The nodes are placed in the center of each control volume. In Figure 27 is shown the control volume discretization represented by small rectangles is depicted. The nodes are placed in the center of each control volume. The nodes located at the sides have no volume and they are useful just to apply the boundary conditions. Figure 27. Discretization of the wall The nodes are classified as follows:  Internal nodes: ; and .  Left-side nodes: ;  Bottom-side nodes: ;  Right-side nodes: ;  Top-side nodes ; Previous calculations, i.e geometric relationships among nodes, are specified in the Annex A.2 Discretization coefficients calculation Internal nodes: The unsteady term of Eq. (60) vanishes as this is a steady state case. In this case . Employing a UDS scheme, from Eq. (59), this results into Page 48 Master Thesis (66) Modeling Equation (66) like Eq. (25), the discretization coefficients result (67) (68) (69) (70) (71) (72) As then Moreover, using a CDS scheme: Eq. (49) results into the final CDS equation Numerical methods in heat transfer and fluid dynamics Page 49 (73) Then the discretization coefficients result (74) (75) (76) (77) (78) (79) For left-, right- and top-side nodes, . Then (80) (81) In the case of bottom-side nodes, . Then (82) (83) And Then (84) Page 50 Master Thesis (85) Then the discretization coefficients result (86) (87) (88) QUICK and SUDS schemes are better introduced with the previously explained ‘deferred convection approach’11. Eq. (59) is modeled like Eq. (25). Then the discretization coefficients result (89) (90) (91) (92) (93) (94) can be computed using the normalized variables Eq. (50), Eq. (55) and Eq.(56). Solver 11 CDS or UDS schemes may be used in the CVs next to the sides since nodes WW, EE, SS and NN may not exist. Numerical methods in heat transfer and fluid dynamics Page 51 After obtaining all the discretization coefficients, the property map is obtained employing the line-by-line methodology (see Solver in 2.2.2). This time is not the temperature, like in heat conduction, but the property which is solved with the line-by-line method. Then for  (95)  (96) 3.2.3. Algorithm The algorithm proposal for the resolution of the Smith-Hutton problem is presented below. 1. Input data  Physical data: ; numerical data: , convective scheme  Geometrical data: ; boundary conditions. 2. Previous calculations12. (See Annex A2) , ; . 3. Guess the property field 4. Calculate the discretization coefficients of all the nodes according the convective scheme employed. 5. Calculate the property map using the line-by-line or GS solver. 6. Apply the convergence criteria: Is a. If no, . (Refresh the guessed property map to the last calculated map); and go to 4. b. If yes, convergence reached. Go to 7. 7. Final calculations and print results. 12 These variables are calculated using only the input data. Page 52 Master Thesis 3.2.4. Verification Table 27 shows the simulation parameters used for the simulation of the Smith-Hutton problem. Table 9. Simulation parameters of Smith-Hutton problem. Convergence criterion 1E-08 Grid 100x50 Relaxation factor Variable Solver TDMA + GS Simulations were performed for Peclet numbers: , using different convective schemes. Figure 28 shows a comparison between the simulation results represented by continuous lines and the benchmark solution data [Smith & Hutton, 1982], represented by dots. Numerical methods in heat transfer and fluid dynamics Page 53 Figure 28. Smith-Hutton problem. Benchmark comparison. 0.0 0.5 1.0 1.5 2.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ф x [m] ρ/Г = 101 Reference UDS 100x50 SUDS 100x50 QUICK 100x50 CDS 100x50 0.0 0.5 1.0 1.5 2.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ф x [m] ρ/Г = 103 Reference UDS 100x50 SUDS 100x50 QUICK 100x50 0.0 0.5 1.0 1.5 2.0 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ф x [m] ρ/Г = 106 Reference UDS 100x50 QUICK 80x40 Page 54 Master Thesis Table 10 represents the average relative error of the simulations for each scheme and Pe. The same results in a tabular format are presented from Table 11 to Table 13. The relative error of each simulation regarding the benchmark solution data is calculated aside. Table 10. Average relative error. UDS13 CDS SUDS QUICK ρ/Г = 10 2.4% 1.1% 0.7% 0.9% ρ/Г = 103 - X 1.2% 0.7% ρ/Г = 106 - X X 10% Table 11. Simulation results for . ρ/Г = 10 Pos Reference UDS 100x50 CDS 100x50 SUDS 100x50 QUICK 100x50 0 1.989 2.000 0.6% 1.871 5.9% 2.000 0.6% 1.904 4.3% 0.1 1.402 1.406 0.3% 1.383 1.4% 1.379 1.6% 1.383 1.4% 0.2 1.146 1.150 0.3% 1.136 0.9% 1.134 1.0% 1.136 0.9% 0.3 0.946 0.950 0.4% 0.940 0.6% 0.939 0.7% 0.940 0.6% 0.4 0.775 0.780 0.6% 0.771 0.5% 0.770 0.6% 0.771 0.5% 0.5 0.621 0.628 1.1% 0.618 0.5% 0.617 0.6% 0.618 0.5% 0.6 0.480 0.489 1.9% 0.478 0.4% 0.478 0.4% 0.478 0.4% 0.7 0.349 0.360 3.2% 0.348 0.3% 0.348 0.3% 0.348 0.3% 0.8 0.227 0.238 4.8% 0.227 - 0.226 0.4% 0.227 - 0.9 0.111 0.123 10.8% 0.111 - 0.111 0.0% 0.111 - 1 0.000 0.000 - 0.000 - 0.000 - 0.000 - 13 It is difficult to assign an average relative error to UDS when ρ/Г = 103, ρ/Г = 106 because it has large relative errors when the property is close to 0. Numerical methods in heat transfer and fluid dynamics Page 55 Table 12. Simulation results for . ρ/Г = 103 Pos Reference UDS 100x50 SUDS 100x50 QUICK 100x50 0 2.000 2.000 - 2.000 - 2.000 - 0.1 1.999 2.000 0.1% 2.000 0.1% 2.000 0.1% 0.2 1.9997 1.9884 0.6% 2.0000 - 1.9999 - 0.3 1.985 1.891 4.7% 1.997 0.6% 1.993 0.4% 0.4 1.841 1.548 15.9% 1.825 0.9% 1.825 0.9% 0.5 0.951 0.960 0.9% 0.956 0.5% 0.967 1.7% 0.6 0.154 0.409 166% 0.141 8.4% 0.149 3.2% 0.7 0.001 0.111 11000% 0.001 - 0.005 400.0% 0.8 0.000 0.017 inf 0.000 - 0.000 - 0.9 0.000 0.001 inf 0.000 - 0.000 - 1 0.000 0.000 - 0.000 - 0.000 - Table 13. Simulation results for. . ρ/Г = 106 Pos Reference UDS 100x50 QUICK 80x40 0 2.000 2.000 2.000 - 0.1 2.000 2.000 - 2.000 - 0.2 2.000 1.997 0.1% 2.000 - 0.3 1.999 1.937 3.1% 1.999 - 0.4 1.964 1.618 17.6% 1.983 1.0% 0.5 1.000 0.975 2.5% 0.996 0.4% 0.6 0.036 0.372 933.3% 0.021 41.7% 0.7 0.001 0.081 8000.0% 0.001 - 0.8 0.000 0.009 inf 0.000 - 0.9 0.000 0.000 - 0.000 - 1 0.000 0.000 - 0.000 - 3.2.5. Results In order to analyze the results and establish a reasonable comparison between the different schemes it is necessary to fix the convergence criterion and the grid size. The purpose is to get an idea of the orders of magnitude of the number of iterations needed by the different the convective schemes to reach the convergence. Page 56 Master Thesis As the solution is given with four decimal positions, a convergence criterion of is expected to be truthful. As the grid is regular and the nodes position are multiples of 0.1, the number of nodes shall be multiple of 10 to have results at the exact positions; then 100 length-nodes are suitable. As the height of the domain is half the length, the mesh is discretized to 100x50 CVs. Figure 29, Figure 30 and Figure 31, shows the property map for the three cases of the problem. These maps were obtained using a UDS scheme. Note that while increases the map turns more symmetric. Figure 29. Smith-Hutton. Property ϕ map for Pe = 10. Figure 30. Smith-Hutton. Property ϕ map for Pe = 103. 10 20 30 40 50 60 70 80 90 100 5 10 15 20 25 30 35 40 45 50 Nodes [j] Nodes [i] Phi 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 10 20 30 40 50 60 70 80 90 100 5 10 15 20 25 30 35 40 45 50 Nodes [j] Nodes [i] Phi 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 Numerical methods in heat transfer and fluid dynamics Page 63 Table 23. QUICK simulation N°2: Parameters Scheme QUICK 10^3 Grid 50x50 Convergence criterion 1E-04 Table 24 shows that the best relaxation factor found is 0.52. Relaxation factor values lower than this factor will increase the number of iterations considerably. As the grid size has been doubled it is expected that the number of iterations increases. Table 24. QUICK simulation N°2: Influence of the relaxation factor. fr Iterations 0.2 1468 0.3 972 0.4 730 0.5 586 0.51 570 0.52 560 0.53 X 0.6 X Note that if a relaxation factor of 0.53 is used the solution diverges (just by increasing the best relaxation factor by 0.01!). Figure 36.QUICK simulation N°2. 0 500 1000 1500 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 Ierations Relaxation factor Page 64 Master Thesis QUICK Simulation N°3: In this simulation is increased to . Table 25. QUICK simulation N°3: Parameters Scheme QUICK 10^6 Grid 50x50 Convergence criterion 1E-04 Table 26 shows that the best relaxation factor found is 0.42. Relaxation factors greater than this will increase the number of iterations considerably. The number of simulations considerably increases regarding simulation N°2 because of increasing . Table 26. QUICK simulation N°3: Influence of the relaxation factor. fr Iterations 0.35 X 0.4 962 0.42 926 0.43 956 0.44 967 0.5 7356 Note how drastically the iterations increase when the relaxation factor is 0.5 Figure 37. QUICK simulation N°3. 0 500 1000 1500 2000 0.35 0.37 0.39 0.41 0.43 0.45 0.47 0.49 Ierations Relaxation factor Numerical methods in heat transfer and fluid dynamics Page 65 4. Incompressible flow method using the Navier- Stokes equations This chapter shows how the classical fractional step projection method [Chorin, 1968; Temam, 1969; Yanenko, 1971] is applied to incompressible flows of Newtonian fluids to solve the velocity-pressure coupling in the NSE14. The fractional step method (FSM) provides an approach that does not use pressure in a predictor step to estimate the projected velocity, an approximate solution of the momentum equation. Consequently, it is required to couple a pressure equation which determines the minimum perturbation that will make the predictor velocity incompressible. It is also important to recall that the role of the pressure in an incompressible flow is to enforce continuity; in some sense, it is more a mathematical variable than a physical one. [Ferziger, 2002]. There are two main goals of this chapter; the first is to depict the FSM employing staggered meshes for its solution, and the second is to assess the performance the numerical schemes presented in Section 3.1.4 in the solution of a model problem, the driven cavity flow. [Ghia, 1982]. 4.1. Introduction to the Fractional Step Method The simulation of incompressible flows of Newtonian fluids is considered. Under these assumptions, the dimensionless governing NS equations in primitive variables are: (97) (98) where Re is the dimensionless Reynolds number defined as (99) 14 High accurate numerical solutions have evidenced that these equations form an excellent mathematical model for turbulent flow. However, more than one and a half century later the understanding of them remains minimal. For instance, the existence and uniqueness of general solutions remains an open question. [CTTC, 2017] Page 66 Master Thesis And and are the density and the dynamic viscosity of the working fluid; and are the characteristic length and velocity, respectively. 4.1.1. Application of the Helmholtz-Hodge decomposition theorem In this section, the Helmholtz-Hodge decomposition theorem is applied to the incompressible NSE. The results here obtained will help us to better elucidate the role of the pressure field. According to the Helmholtz-Hodge theorem, a given vector field ω, defined in a bounded domain with smooth boundary , is uniquely decomposed in a pure gradient field and a divergence-free vector parallel to . (100) where, (101) The theorem also applies for periodic inflow/outflow conditions. Let be a projector operator. It projects any vector field onto a divergence-free space15 (102) Now, taking the NS equations and applying the projector operator (103) Since the velocity field is incompressible the transient term remains unchanged (see Eq. I) when projected 15 Note that any incompressible field remains unchanged by the action of the projector operator , (I) moreover, the projector operator , has the following property (II) Numerical methods in heat transfer and fluid dynamics Page 67 (104) Whereas the projection of the pressure gradient vanishes (see Eq. II) (105) Therefore, the NS equations can be split in two parts: a divergence free-vector and a gradient of a scalar field. (106) (107) The Helmholtz-Hodge decomposition theorem ensures that this decomposition is unique. Finally applying the divergence operator to Eq. (107) and using projector definition, see Eq. (102), leads to a Poisson equation for pressure. (108) Therefore, for incompressible flows the role of the pressure gradient is to project the vector field into a divergence-free space. Figure 38. Convective and viscous term vector field unique decomposition. [CTTC, 2017] 4.1.2. Time integration method The final form of the fractional step method would depend on the time-integration method chosen. Here, for the sake of clarity, it is proposed a fully explicit time integration scheme. To simplify the notation, momentum equation can be rewritten as Page 68 Master Thesis (109) where stands for the convective and diffusive terms (110) For the temporal discretization, a central difference scheme is used for the time derivative term, (111) a fully explicit second-order Adams-Bashforth scheme for (112) and a first order backward Euler scheme for the pressure-gradient term. Incompressibility constraint is treated implicitly. Thus, we obtain the semi-discretized NSE. (113) (114) Then, the pressure gradient forces (projects) the predictor velocity field to be incompressible . This projection is derived from the well-known Helmholtz- Hodge vector decomposition theorem, whereby the predictor velocity , can be uniquely decomposed into a divergence-free vector, , and the gradient of a scalar field, . This decomposition is written as (115) where the predictor velocity is given by (116) and the pseudo-pressure is Taking the divergence of Eq. (115) yields a Poisson equation for . Numerical methods in heat transfer and fluid dynamics Page 69 (117) Once the solution is obtained, results from the correction (118) Therefore, the algorithm for the integration of each time step is 1. Evaluate . 2. Evaluate from Eq. (116). 3. Evaluate and solve the discrete Poisson equation, Eq. (117). 4. Obtain the new velocity field with Eq.(118). 4.1.3. Determination of Due to stability reasons, explicit temporal schemes introduce severe restrictions on the time step, while implicit discretization would improve the overall stability. Here, we have only considered explicit method in the view of formulation simplicity. Therefore, in our case, the time-step, , must be bounded by the CFL condition [Courant, 1928] given by (119) (120) where the bounding values are must be smaller than unity. In our case, we will follow the recommendations given by Simons (2000) using values and , respectively.16 4.1.4. The Poisson equation Since the formulation is fully explicit, the only system to be solved is the pressure Poisson Eq. (115). An efficient solution of the Poisson equation is critical since it may consume a 16 Note that, since the mesh would remain constant during all the simulation, only the CFL-condition for the convective term, see Eq. (119), has to be recomputed each time-step. Page 70 Master Thesis substantial portion of the total computing time. If we focus in the 1D spatial discretization of the step 3 and after applying finite differences at node P: (121) For the x-component of the velocity: Figure 39.1D spatial discretization. [CTTC, 2017] (122) The discrete approximation of at the node P is independent of . Then velocity fields for unphysical pressure distributions can be obtained. For example, Since at node P is independent of , the final velocity field will verify This is known as ‘the checkerboard problem’. Then it is needed a smarter strategy to couple with the velocity field . The staggered meshes solve the checkerboard problem. Typical staggered control volumes are shown in Figure 41. Figure 40.The checkerboard problem. [CTTC, 2017] Numerical methods in heat transfer and fluid dynamics Page 71 Figure 41. Staggered meshes. [CTTC, 2107] Figure 42 shows that the control volumes for u and v, are displaced with respect to the control volume for the continuity equation. The solution of the continuity equation is decoupled from the momentum equations by a fractional step method which enforces mass conservation by solving a Poisson equation. This procedure, combined with the consistent approximations of the geometric quantities, is done to satisfy the discretized mass conservation equation to machine accuracy, as well as to gain the favorable convergence properties of the Poisson solver [Rosenfeid, 1992]. The unsteady resolution advances with adaptive time steps until a specified condition is reached, e.g. steady state. Figure 42. Control volumes for a staggered grid: for mass conservation and scalar quantities (left), for x-momentum (center) and for y-momentum (right). [Ferziger, 2002] The simultaneous solution of the large number of resulting discrete equations is very costly, especially for three-dimensional cases. An efficient approximate solution can be obtained by Page 72 Master Thesis decoupling the solution of the momentum equations from the solution of the continuity equation by the fractional step method. The basic fractional step (or projection) method was proposed by Chorin (1968). Despite all the inherent difficulties17, the fractional step method has become a very popular technique for solving the incompressible NS equations due to its good performance and code simplicity. Fractional step methods are also referred to as projection methods because the system of equations given can be interpreted as a projection into a divergence-free velocity space. The predictor velocity, is an approximate solution of the momentum equations, but because the predictor velocity is obtained with no pressure gradient contribution it cannot satisfy the incompressibility constraint at the next time level. Then, the Poisson equation determines the minimum perturbation that will make the predictor velocity incompressible [CTTC, 2017]. 4.2. The Driven cavity flow problem The Driven Cavity Flow problem is a model problem that has served repeatedly for testing and evaluating numerical techniques, despite the singularities at the two of its corners. For moderately high values of the Reynolds number, published results are available for this flow problem from many sources [Ghia, 1982]. Some results are also available for high Re [Nallasamy, 1977], but the accuracy of most of these high-Re solutions has generally been viewed with some skepticism because of the size of the computational mesh employed and the difficulties experienced with the convergence. 4.2.1. Problem definition A laminar incompressible flow in a square cavity whose top wall moves with a uniform velocity in its own plane, see Figure 43. Velocities and are zero in the right, left and bottom walls; in the top side, the velocity is horizontal and equals to . The boundary conditions are that the pressure gradient normal to the walls is zero. The cavity is modeled as a 1x1 m2 square. It is required to solve the velocity field (x and y directions), in particular, numerical values corresponding to the velocity profiles for horizontal and vertical lines passing through the geometric center of the cavity. Cases to be assessed are . 17 There are still some controversy related with the temporal accuracy of the method and some debate about how the boundary conditions should be specified. Numerical methods in heat transfer and fluid dynamics Page 79 Figure 47. Boundary condition: prescribed velocity. [CTTC, 2017] Step 3 - Staggered-x mesh Figure 48. Geometric relation of the pressure mesh and the x-staggered mesh. [CTTC, 2017] Note that for the velocity of the x-staggered mesh, and are the correspondent pressures for and respectively. Likewise, corresponds to . Step 3 - Staggered-y mesh Figure 49. Geometric relation of the pressure mesh and the y-staggered mesh. [CTTC, 2017] Page 80 Master Thesis Note that for a velocity at the x-staggered mesh, and are the correspondent pressures for and respectively. Likewise, corresponds to . Step 4 - Choice of the time step. The CFL (Courant-Friedrich-Levy) condition (142) (143) (144) 4.2.3. Algorithm The algorithm proposal for the resolution of the Driven cavity flow problem is presented below. [1] Input data  Physical data:  Numerical data: , , convective scheme.  Geometrical data:  Boundary conditions. at boundaries, except at the upper boundary. [2] Previous calculations  Creation of the regular mesh for the pressure.  Creation of the staggered-x mesh for x-velocity.  Creation of the staggered-y mesh for y-velocity.  Velocity fields creation:  Calculation of [3] Calculation of the predictor velocities: 3.1 Numerical methods in heat transfer and fluid dynamics Page 81 3.2 [4] Guess [5] Calculate the discretization coefficients [6] Calculate [7] Calculate [8] Calculate the next time-step pressure field. [9] Check the pressure field convergence:  If No: Return to [3].  If Yes: Pressure field convergence, go to 10. [10] Calculation of the next time-step velocity fields, [11] Velocities update for the next time-step , , [12] Calculation of the new time-step [13] Steady state condition verification  If No: Return to [3] Page 82 Master Thesis  If Yes: Convergence reached. Go to [14]. [14] Final calculations and print results.  Calculate the velocities at the CVs’ faces  Plot 4.2.4. Verification Table 27 shows the simulation parameters used for the simulation of the Driven cavity flow problem. Table 27. Simulation parameters of the Driven Cavity Flow problem. Grid 32x32,64x64, 128x128 Pressure relaxation factor 0.9 Pressure convergence criterion 1E-06 Velocities convergence criterion 1E-06 Reynolds number 100, 400, 1000, 3200 Convective schemes UDS, CDS, SUDS, QUICK, SMART Solver TDMA + GS The results obtained are compared with the benchmark data. The results obtained are compared with the benchmark data. Table 28 and Table 29 list the numerical values corresponding to the velocity profiles for lines passing through the geometric center of the cavity. Only typical points, rather than the entire large set of computational points, along these profiles have been listed. The points of local maxima and minima for all values of Re are underscored. Numerical methods in heat transfer and fluid dynamics Page 83 Table 28. Benchmark solution: u Table 29. Benchmark solution: v From Figure 50 to Figure 57, verification graphics are depicted comparing x- and yvelocities of the simulations with the benchmark solution data. For Reynolds 100, 400 and 1000 satisfactory results have been obtained adjusting to a proper mesh size. Nevertheless the results for Re 3200 has not been that much accurate. The simulation Page 84 Master Thesis code could not evolve favorably for Re 5000, 7500 and 10000. In Table 30, the average relative error of the simulations for each scheme, mesh and Re is presented. Figure 50. Verification Re 100, x-velocities. Figure 51. Verification Re 400, x-velocities. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 y [m] u [m/s] Verfication - Re 100 Benchmark solution CDS 64x64 QUICK 64x64 SMART 64x64 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 y [m] u [m/s] Verification - Re 400 Benchmark solution CDS 64x64 QUICK 64x64 Numerical methods in heat transfer and fluid dynamics Page 85 Figure 52. Verification Re 1000, x-velocities. Figure 53. Verification Re 3200, x-velocities. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 y [m] u [m/s] Verification - Re 1000 Benchmark solution QUICK 128x128 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 y [m] u [m/s] Verification - Re 3200 Benchmark solution QUICK 128x128 Page 86 Master Thesis Figure 54. Verification Re 100, y-velocities. Figure 55. Verification Re 400, y-velocities. -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 v [m/s] x [m] Verification - Re 100 Benchmark solution CDS 64x64 QUICK 64x64 SMART 64x64 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 v [m/s] x [m] Verification - Re 400 Benchmark solution CDS 64x64 QUICK 64x64 Numerical methods in heat transfer and fluid dynamics Page 87 Figure 56. Verification Re 1000, y-velocities. Figure 57. Verification Re 3200, y-velocities. -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 v [m/s] x [m] Verification - Re 1000 Benchmark solution QUICK 128x128 -0.60 -0.55 -0.50 -0.45 -0.40 -0.35 -0.30 -0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 v [m/s] x [m] Verification - Re 3200 Benchmark solution QUICK 128x128 Page 88 Master Thesis Table 30. Average relative error. Mesh Reynolds UDS CDS QUICK SMART 32x32 Re 100 35.5% 4.8% 2.7% 2.7% Re 400 X 9.8% 8.4% X Re 1000 X 12.2% 7.2% X 64x64 Re 100 10.0% 2.3% 2.2% 2.2% Re 400 X 2.1% 1.6% X Re 1000 X 3.5% 2.6% X 128x128 Re 1000 - - 1.1% X Re 3200 - - 16.1% X Numerical methods in heat transfer and fluid dynamics Page 95 4.2.5. Results Figure 58 shows the pressure map of the domain. Cases Re = {100, 400, 1000} are shown. In fact, pressure field as an absolute value is not relevant since there is no explicit boundary pressure value in the problem. What is relevant is the pressure gradients within the mesh19. Figure 58. Results p [Pa] for Re 100, 400, 1000. 19 Consequently, is not surprising to find negative pressures, this apparent error may be prevented just by fixing a pressure value for one single node, the velocity field is not affected. 0 0 0 0 50 Pressure [Pa] Re 100 10 20 30 40 50 60 10 20 30 40 50 60 -100 -50 0 50 100 150 0 0 0 0 0 0 0 0 50 100 150 200 Pressure[Pa] Re 400 10 20 30 40 50 60 10 20 30 40 50 60 -50 0 50 100 150 200 250 300 350 0 0 0 0 0 0 0 0.1 0.2 0.3 Pressure [Pa] Re 1000 20 40 60 80 100 120 20 40 60 80 100 120 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Page 96 Master Thesis The x-velocity of the central vertical line and the y-velocity of the central horizontal line are shown in Figure 59 and Figure 60, respectively. Cases Re = {100, 400, 1000, 3200} are plotted. There are seen kinks in the u-profiles near y = 1 and in the v-profiles near x = 1. This behavior persists for greater Re cases. Such behavior has been reported by some previous investigators. This would imply that the velocity distributions near these walls are not too much sensitive to mesh size. Figure 59. Results u [m/s] for Re 100, 400, 1000. 0 0 0 0 0 0.2 0.2 0.2 0.2 0.4 0.4 0.4 -0.2 -0.2 -0.2 0.6 0.6 0.6 0.8 0.8 0.8 u [m/s] Re 100 10 20 30 40 50 60 10 20 30 40 50 60 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 u [m/s] Re 400 10 20 30 40 50 60 10 20 30 40 50 60 -0.2 0 0.2 0.4 0.6 0.8 0 0 0 0 0 0 -0.2 -0.2 -0.2 -0.2 -0.2 0.2 0.2 0.2 0.2 0.4 0.4 0.4 0.6 0.6 0.6 0 0 0.8 0.8 0.8 0 u [m/s] Re 1000 20 40 60 80 100 120 20 40 60 80 100 120 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Numerical methods in heat transfer and fluid dynamics Page 97 Figure 60. Results v [m/s] for Re 100, 400, 1000. Figure 61 shows the velocity vectors for the flow configurations for Re = {100, 400, 1000} ; and Figure 62, for Re = 3200. Figure 61. Velocity vectors for Re = 100, 400, 1000. 0 0 0 0.1 0.1 0.1 0.1 0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.2 -0.2 -0.2 -0.2 0.2 0.2 0.2 0.2 -0.3 -0.3 -0.3 -0.4 -0.4 -0.5 0 v [m/s] Re 100 10 20 30 40 50 60 10 20 30 40 50 60 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0 0 0 0 0.1 0.1 0.1 0.1 0.1 0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.2 0.2 0.2 0.2 -0.2 -0.2 -0.2 -0.2 -0.2 -0.3 -0.3 -0.3 -0.3 -0.4 -0.4 -0.4 -0.5 -0.5 0.3 0.3 0 -0.6 v [m/s] Re 400 10 20 30 40 50 60 10 20 30 40 50 60 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0 0 0 0 0.1 0.1 0.1 0.1 0.1 0.1 -0.1 -0.1 -0.1 -0.1 -0.1 -0.1 0.2 0.2 0.2 0.2 -0.2 -0.2 -0.2 -0.2 -0.2 -0.3 -0.3 -0.3 -0.3 0.3 0.3 0.3 -0.4 -0.4 -0.4 -0.4 -0.5 -0.5 -0.5 0 -0.6 -0.6 0.1 0 v [m/s] Re 1000 20 40 60 80 100 120 20 40 60 80 100 120 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Page 98 Master Thesis Figure 62. Velocity vectors for Re = 3200. 020 40 60 80 100 120 140 0 20 40 60 80 100 120 140 Numerical methods in heat transfer and fluid dynamics Page 99 4.2.6. Discussion  CDS and QUICK reached the convergence for Re 100. 400 and 1000, with a comparable accuracy for the lowest Re. CDS is the simplest scheme of secondorder accuracy and offers a good compromise among accuracy, simplicity and efficiency. On the other hand, QUICK performance is much better as Re is increased. QUICK scheme seems to ensure diagonal dominance for the resulting algebraic equations, thus lending the necessary stability property to the evolving solutions while restoring second-order accuracy at convergence.  Among the different schemes UDS performs the worst accuracy. SMART scheme showed an accuracy as good as the QUICK scheme, but presents stability issues.  From simulations carried out with meshes 32x32 and 64x64, it is inferred that 128x128 meshes are adequate for moderate values of Re.  As Re increased, very coarse grids could not be included in the procedure, and mesh refinement is mandatory; nevertheless, the corresponding computer time is quite large, i.e., convergence for with coarse meshes is not possible using any of the schemes presented and finer meshes, from about 256x256, are necessary.  The time step calculated with the CFL condition does not always guarantee the convergence, therefore a timestep corrective factor, , is introduced. (145) This corrective factor is seen to be influenced by the physical problem parameters, namely, Re and mesh discretization. Table 37 shows the time-step corrective factor needed to reach the convergence using a QUICK scheme. Table 37. Time-step corrective factor20 Mesh vs. Re Re 100 Re 400 Re 1000 Re 3200 Re 5000 8x8 100% 100% 90% X - 16x16 100% 100% 80% X - 32x32 100% 100% 75% X - 64x64 40% 80% 50%* X - 128x128 X X 100% 50%* X 20 ‘X’ stands for divergence, and (*) means the simulation did not converge but results were taken as they might be orientative. Page 100 Master Thesis  A pressure relaxation factor, , was used in all the cases to achieve the pressure field convergence.  The relative error is significant near the walls, this may be solved refining the mesh or using an irregular mesh with higher density of CVs near the walls (as seen in Chapter 5).  It is interesting to study the influence of parameters and for a given Re. For example, for a , the computing time is much lower for a set of and than for a set of and .  Difficulties have been experienced with the convergence for high Re, mainly due to the limits of the computer processor used. The size of the computational mesh employed must be increased and the use of a workstation would be much helpful. Numerical methods in heat transfer and fluid dynamics Page 101 Conclusions Heat conduction methods  An understanding of the heat conduction solution methodology is essential. To be able to see the similarities between transfer of momentum and transfer of heat and to regard velocity as, in some ways, analogous to temperature, is a great conceptual help. The use of heat conduction as a building block in the fluid-flow calculation scheme reinforces this conceptual unity.  The TDMA is a very powerful and convenient equation solver whenever the algebraic equations can be represented with the generic conservation equation. In the ‘Four materials problem’ it was noticed that the TDMA requires time proportional to N, the number of CVs, rather than to N2 or N3. The convection-diffusion equation  Many of the features commonly encountered in practical convection-diffusion problems were present in the Smith-Hutton problem. The streamline pattern is fairly complicated, exhibiting a large degree of curvature (a general feature of recirculating flows), and the transported variable displays a rapid change over a small distance within the flow (possibly as a consequence of a source, the mixing of two streams of different temperature, etc.).  According to the relaxation factor analysis applied to the Smith-Hutton problem, it has been determined that there are no general rules for choosing the best value of a relaxation factor. The optimum value depends upon many factors, such as the nature of the problem, the number of grid points, the grid spacing, and the iterative process used. Suitable values of were found from exploratory computations for the given problem. Finding an appropriate relaxation factor for high-order schemes is in some cases a necessary requirement and shall be faced by first exploring the behavior of the system in coarse meshes and progressively refining them. The optimal relaxation factor will depend on the numerical, physical and geometrical factors of the problem.  The study of different schemes for the discretization of the convective term in the Page 102 Master Thesis convection-diffusion equation evidenced that: o There does not appear to be a perfect scheme. On regular coarse meshes, all the methods submitted showed evidence of false diffusion or spatial oscillations or some combination of the above. o UDS is characterized by false diffusion, which completely swamps the real diffusion at high Peclet numbers. The scheme gives physically plausible results, but on practical meshes it appears to be impossible to make accurate predictions at high Pe. o CDS tend to give false spatial oscillations when the mesh refinement is insufficient to resolve variations of the dependent variable that occur in the problem. If the mesh is refined an accurate solution will result. This behavior could be considered an advantage, given that the oscillations are a warning that the grid is insufficient to give an accurate solution. o High-order schemes yield performances on coarse meshes that represent a compromise between diffusive and spatial oscillation errors. It is noteworthy, however, that these schemes do not always give good results at low Pe, a matter of some concern if calculating a complex flow containing regimes of differing Peclet numbers. Incompressible flow method  The robustness of the FSM has been demonstrated using the model problem ‘Driven cavity flow’. Detailed accurate results have been presented for this problem. Up to 128 x 128, i.e. 16384 computational points and Re as high as 3200 have been considered.  Thanks to my simulation experience acquired throughout this project it is important to emphasize that the solution convergence rate can be seriously affected if the coupling among the many governing differential equations is not properly respected either in the interior of the solution domain or at its boundaries. The rate of convergence is also generally strongly dependent on such problem parameters as the Reynolds number, the mesh size, and the total number of computational points. Numerical methods in heat transfer and fluid dynamics Page 103 Future actions  Irregular meshes vs. regular meshes performance comparison. Look for improvements for those nodes near to boundaries.  Extension of the convection-diffusion equation and NS equations for incompressible flow to 3D meshes, this is, including the top and bottom fluxes and nodes.  Use of the blocking-off technique for complex geometries in heat conduction problems.  Introduction of the energy conservation equation to the NS equations. Solution of the temperature field considering the buoyancy and temperature gradient effects. Currently I am working on the benchmark problem called “Differentially-heated cavity”.  Turbulent flow modeling, using the Burgers’ equation in Fourier space. Page 104 Master Thesis Bibliography 1. Suhas V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, McGraw-Hill Book Company, 1980. 2. Ferziger, Joel H. Computational Methods for Fluid. Berlin; Heidelberg; New York; Barcelona; Hong Kong; London; Milan; Paris; Tokyo: Springer, 2002. ISBN 3- 540-42074-6 3. CTTC, Centre Tecnològic de Transferència de Calor, A Two-dimensional Transient Conduction Problem, Universitat Politècnica de Catalunya, 2016. 4. CTTC, Centre Tecnològic de Transferència de Calor, Convection-diffusion equations, Universitat Politècnica de Catalunya, 2017. 5. Darwish, M. S., Moukalled, F. H. Normalized variable and space formulation methodology for high-resolution schemes. Numerical Heat Transfer, Part B: Fundamentals, 26:1,79 — 96; 1994 6. R. M. Smith & A. G. Hutton. The numerical treatment of advection: a performance comparison of current methods, Numerical Heat Transfer, 5:4, 439-461, 1982. 7. CTTC, Centre Tecnològic de Transferència de Calor, Introduction to the Fractional-Step Method, Universitat Politècnica de Catalunya, 2017. 8. Ghia et al., High-Re Solutions for Incompressible Flow Using Navier-Stokes Equations and a Multigrid Method, Journal of Computational Physics 48, 387- 411, 1982. 9. A. J. Chorin. Numerical Solution of the Navier-Stokes Equations. Journal of Computational Physics, 22:745–762, 1968. 10. R. Temam. Sur l’approximation de la solution des ´equations de Navier-Stokes par la methode des pas fractionnaires (II). Archive for Rational Mechanics and Analysis, 33:377–385, 1969. 11. N. N. Yanenko. The Method of Fractional Steps. Springer-Verlag, 1971. 12. M. Rosenfeid, D. Kwak, and M. Vinokur. Development of a Fractional-Step Method for the Unsteady Incompressible Navier-Stokes Equations in Generalized Numerical methods in heat transfer and fluid dynamics Page 111 (174) The mass flow at face ‘f’ is given by (175) The CVs’ surfaces and volume are constant throughout the domain and do not depend on their position. (176) (177) (178) Page 112 Master Thesis B. Forced convection inside an isothermal tube Problem definition  Fluid flow inside an isothermal tube of internal diameter Di and length L.  Inlet flow conditions are known: . Figure 63. Scheme of the isothermal tube problem. The objective of the exercise is the evaluation of the average exit velocity, pressure and temperature: and . The total heat flux exchanged between the wall and the fluid, must also be calculated. Input data. Three different cases:  Fluid: water (case 1), thermal oil Therminol 66 (case 2), air (case 3)  Geometry: (smooth tube)  Fluid conditions at the entrance. (case 1 and 2), (case 3), ), .  Temperature of the tube: Code development Assumptions  One dimensional flow, this implies using mean values at each cross-sectional area and local heat transfer coefficients. Numerical methods in heat transfer and fluid dynamics Page 113  Steady-state  Constant temperature of the wall.  Empirical inputs needed (local heat transfer coefficient, friction factor and thermophysical properties) for the different fluids are evaluated from the expressions given in the Formulae (sections C2, C7 and D0).  Heat transfer phenomena: only convection. Domain discretization Figure 64. Discretization of the isothermal tube. N is the number of control volumes (CV) and N+1 is the number of nodes. The solution approach is to calculate the fluid temperature, pressure and velocity, from node to node and use the conservative equations for each CV. This is, formulate the conservative equations for CV 1 and calculate the fluid properties of node 2, formulate the conservative equations for CV 2 and calculate fluid properties of node 3, and so on until arriving to CV N and calculate the fluid properties of node N+1. The length of each control volume in the axial axis, , is given by: Note that the cross-sectional surface, , is constant and it is given by: Page 114 Master Thesis Mass conservation equation The unsteady term can be neglected since steady state is assumed. The transport term represents the net mass flow that exits the control surfaces. The mass flow is constant in all the control volumes. The density of the fluid can be evaluated using algebraic correlations given in the Formulae Section D. This equation is known as the state equation. The inlet density, , must be evaluated to obtain the mass flow . The fluid density depends on the temperature and pressure of the fluid, this is the case for the air, but it is approximated as a function only of the temperature in the case of water and Therminol 66. From Equation 4 and considering the control surface constant across the tube, the mass conservation equation turns into: Numerical methods in heat transfer and fluid dynamics Page 115 Momentum conservation equation The unsteady term can be neglected since steady state is assumed. The transport term represents the net linear momentum that exits the control surfaces. The third term represents the forces acting over the control surfaces and it is positive in the coordinate positive direction, this implies the fluid forces crossing the inlet and outlet of the control volume and the viscous forces fluid-tube ( ). The viscous tensions are neglected. The fourth term is equal to 0, since the gravity in the axial coordinate does not exist. It is assumed that is constant across the control surface, then it can be placed out of the integral. is given by: is the skin friction coefficient and depends on the cross-sectional shape, the hydraulic diameter (in this exercise , the Reynolds number (Re) and the relative roughness ( ). In this exercise the cross-sectional shape is a circular tube and the relative roughness is (smooth tube). Considering this the friction factor results: Page 116 Master Thesis is the average density in the control volume ‘i’ and it is evaluated at the mean temperature (and the mean pressure in the case of air) of the control volume. Where is the average velocity in the control volume ‘i’ and it is the arithmetic mean of and . Then the Equation 10 can be expressed as: Energy conservation equation The unsteady term can be neglected since steady state is assumed. The transport term represents the net energy that exits the control surfaces, here, the potential energy variation is neglected, and the kinetic energy variation can be neglected for liquids if wanted. The third term represents the net heat flux that enters through the control surfaces. The fourth term Numerical methods in heat transfer and fluid dynamics Page 117 represents the work developed and in this exercise is equal to 0. The third term in this exercise is given by heat exchanged between the tube and the fluid through the cylindrical surface. Considering the heat transfer coefficient constant across a CV, and taking the mean temperature of the fluid in the CV, we obtain If the fluid is air, is treated as a semi-perfect gas, then the enthalpy differential is If the fluid is a liquid (in the case of water and Therminol 66), the enthalpy differential is: The component can be neglected as it is some orders of magnitude smaller than . Then it is obtained In forced convection inside ducts, to obtain the heat transfer coefficient, , the thermophysical properties of the fluid (dynamic viscosity, thermal conductivity and specific heat) must be obtained first, and these are evaluated at the mean fluid bulk temperature, . , , , . refers to the fluid viscosity if it were at the temperature of the tube. Next, the non-dimensional numbers are obtained: Re, Pr, Gz and Nu. Page 118 Master Thesis For laminar flow ( < 2000): For turbulent flow ( > 2000): Finally, the heat transfer coefficient of the CV ‘i', , is given by Finally, this heat transfer coefficient is used in the energy conservation equation. Numerical methods in heat transfer and fluid dynamics Page 119 Algorithm 1. Input data: Physical data: . Numerical data: . 2. Previous calculations: . 3. Initial map: For node i = 1. where . The fluid properties at the inlet are already known. 4. Evaluation of node ‘i+1’.  Guess the fluid properties at node ‘i+1’: , where  Evaluate the thermophysical properties of the fluid: at the mean bulk temperature and pressure, and .  Calculate the non-dimensional numbers (Re, Pr, Gz, Nu), the heat transfer coefficient ( ), and the friction factor .  Calculate the fluid properties from the conservative equations  From the linear momentum conservation equation, calculate .  From the energy conservation equation, calculate .  From the state equation, calculate .  From the mass conservation equation, calculate . 5. Convergence criteria: Is , where If no: . (Refresh the guessed fluid properties to the last calculated values); and go to 4.2. If yes: Go to 6. 6. Is i = N+1? (Evaluation of the last node). If NO: Update the node: i i+1. Go to 4.1 If YES: Go to 7. 7. Final calculations: and . 8. Print results: vs. x, where and x is the axial coordinate. Page 120 Master Thesis Results Figure 65 shows the temperature distribution for the three fluids. The air is the fluid which more effectively reaches the closest outlet temperature to the tube temperature, followed by the water which considerably rises its temperature, and lastly the Therminol 66 which hardly rises its temperature. The nature of the curves shows that air increases its temperature more rapidly at lower temperatures than at higher ones. Water and Therminol 66 increase their temperature with less dependence of the temperature, close to the linearity. Figure 66 shows pressure drop in all the fluids, but the decrease is very small in the case of water and relatively small in the case of air, on the other hand, the pressure decrease in Therminol 66 is considerable. Figure 67shows that the velocity of the liquids hardly increases, this validates the simplification of neglecting the kinetic energy variation in the energy conservation equation for liquids. On the other hand, the velocity increase of the air is more notorious, this explains its pressure drop. Figure 68 shows that the water is by far the fluid which exchanges more heat flux with the hot tube (two orders of magnitude larger than the other fluids), the second is the Therminol 66 and the third is the air. As expected, the lower the fluid temperature the larger is the heat exchanged. Figure 71 shows that the water is by far the fluid which has the largest heat transfer coefficients along the tube (one order of magnitude larger than the other fluids), the air and Therminol 66 come then. In the case of water and air, the heat transfer coefficients increase downstream, but in the case of Therminol 66, the heat transfer coefficient decreases downstream. Moreover, the variation of the heat transfer coefficient is relatively low in the case of air and Therminol 66, and large in the case of water. Figure 73 shows that the viscous shear stress is more significant for Therminol 66 than for water and air, this explains why the pressure drop is significant in this fluid. In the case of Therminol 66 the shear stress decreases downstream. Figure 75, Figure 76 and Figure 77 show the intensive energy distribution (enthalpy and kinetic energy) in water, Therminol 66 and air, respectively. Second order approximations are used for the calculation. The differential forms of intensive energy in a control volume are given by: (1) (2) Numerical methods in heat transfer and fluid dynamics Page 127 Grid sensitivity Table 41 shows the total heat exchanged between the fluids and the tube Table 41. Total heat exchanged between the fluids and the tube. Total heat exchanged Water Therminol 66 Air Q_total [W] 78167 3864 1616 The following figures show the grid sensitivity in function of the fluid outlet temperature, , the fluid outlet pressure, , and the fluid outlet velocity , , and the total heat flux exchanged, . In all the cases, there is a significant relative error for values of N below 10. From onwards, the relative error is practically zero. 010 20 30 40 50 60 70 80 90 100 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 Number of control volumes (N) Relative error [%] Grid sensitivity test in function of Tout Water Therminol 66 Air Page 128 Master Thesis 010 20 30 40 50 60 70 80 90 100 0 0.05 0.1 0.15 0.2 0.25 Number of control volumes (N) Relative error [%] Grid sensitivity test in function of Pout Water Therminol 66 Air 010 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Number of control volumes (N) Relative error [%] Grid sensitivity test in function of vout Water Therminol 66 Air 010 20 30 40 50 60 70 80 90 100 0 2 4 6 8 10 12 Number of control volumes (N) Relative error [%] Grid sensitivity test in function of Qtotal Water Therminol 66 Air Numerical methods in heat transfer and fluid dynamics Page 129 Convergence criterion sensitivity From Figure 78 to Figure 81 the convergence criterion sensitivity is plotted in function of the fluid properties , , and . The x-axis represents the values taken by , and the y-axis represents the relative error regarding the results obtained when . increases its value in steps of 0.001 until arriving to , this is from until . The figures show that from to at least the relative error is practically 0.00%. Hence is a good enough value for validating the code. Consequently, the value of used to obtain the previous results, , is good enough as well. Figure 78. Convergence criterion sensitivity in function of Figure 79. Convergence criterion sensitivity in function of 0100 200 300 400 500 600 700 800 900 1000 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 Delta (x1e-3) Relative error [%] Convergence criteria sensitivity in function of Tout Water Therminol 66 Air 0100 200 300 400 500 600 700 800 900 1000 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 Delta (x1e-3) Relative error [%] Convergence criteria sensitivity in function of pout Water Therminol 66 Air Page 130 Master Thesis Figure 80. Convergence criterion sensitivity in function of Figure 81. Convergence criterion sensitivity in function of Discussion If the objective were reaching the temperature of the tube:  A double-length tube is needed in the case of air for reaching the temperature of the tube in 99.8%  A double-length tube is needed in the case of water for reaching the 0100 200 300 400 500 600 700 800 900 1000 0 0.02 0.04 0.06 0.08 0.1 0.12 Delta (x1e-3) Relative error [%] Convergence criteria sensitivity in function of vout Water Therminol 66 Air 0100 200 300 400 500 600 700 800 900 1000 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Delta (x1e-3) Relative error [%] Convergence criteria sensitivity in function of Qtotal Water Therminol 66 Air Numerical methods in heat transfer and fluid dynamics Page 131 temperature of the tube in 97.4% Figure 82. Temperature distribution for a 10 m. long tube.  A half-diameter tube is needed in the case of water for reaching the temperature of the tube in 99.9%  A half-diameter tube is needed in the case of water for reaching the temperature of the tube in 98.5% 0 5 10 15 20 30 40 50 60 70 80 90 100 X: 10 Y: 94.81 Axial position [m] Temperature [°C] Temperature distribution Tube Water Therminol 66 Air 8 9 10 11 12 13 14 15 88 90 92 94 96 98 100 X: 10 Y: 92.55 Axial position [m] Temperature [°C] Temperature distribution Tube Water Therminol 66 Air Page 132 Master Thesis Figure 83. Temperature distribution for a Di = 10 mm. tube 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 20 30 40 50 60 70 80 90 100 X: 5 Y: 93.6 Axial position [m] Temperature [°C] Temperature distribution Tube Water Therminol 66 Air