An Introduction to Computational Fluid Dynamics Chapter 20 in Fluid Flow Handbook By Nasser Ashgriz & Javad Mostaghimi Department of Mechanical & Industrial Eng. University of Toronto Toronto, Ontario

1 1 Introduction:................................................................................................................ 2 2 Mathematical Formulation.......................................................................................... 3 2.1 Governing equations.............................................................................................. 3 2.2 Boundary Conditions............................................................................................. 5 2.2.1 Example....................................................................................................... 7 3 Techniques for Numerical Discretization................................................................... 9 3.1 The Finite Difference Method............................................................................... 9 3.2 The Finite Element Method................................................................................. 11 3.3 The Finite Volume Method ................................................................................. 14 3.4 Spectral Methods................................................................................................. 15 3.5 Comparison of the Discretization Techniques..................................................... 16 4 Solving The Fluid Dynamic Equations..................................................................... 17 4.1 Transient-Diffusive Terms .................................................................................. 17 4.1.1 Finite Difference Approach....................................................................... 17 4.1.2 Finite Element Approach........................................................................... 21 4.2 Transient-Convective Terms ............................................................................... 24 4.3 Shock Capturing Methods ................................................................................... 26 4.4 Convective-Diffusive Terms ............................................................................... 27 4.5 Incompressible Navier-Stokes Equations............................................................ 30 4.5.1 Pressure-Based Methods............................................................................ 30 5 Basic Solution Techniques........................................................................................ 34 5.1 Direct Method...................................................................................................... 34 5.2 Iterative Methods................................................................................................. 34 5.2.1 Jacobi and Gauss-Seidel methods.............................................................. 35 5.2.2 Relaxation methods. .................................................................................. 37 5.2.3 ADI Method: ............................................................................................. 38 5.3 Convergence and Stability................................................................................... 39 5.4 Von Neuman Stability Analysis .......................................................................... 39 5.5 Convergence of Jacobi and Gauss-Seidel Methods (iterative methods): ............ 41 6 Building a Mesh........................................................................................................ 44 6.1 Element Form ...................................................................................................... 44 6.2 Structured Grid .................................................................................................... 46 6.2.1 Conformal mapping method...................................................................... 47 6.2.2 Algebraic method ...................................................................................... 47 6.2.3 Differential equation method..................................................................... 47 6.2.4 Block-structured method ........................................................................... 47 6.3 Unstructured grid................................................................................................. 47 7 References................................................................................................................. 49

2 1 Introduction: This chapter is intended as an introductory guide for Computational Fluid Dynamics CFD. Due to its introductory nature, only the basic principals of CFD are introduced here. For more detailed description, readers are referred to other textbooks, which are devoted to this topic.1,2,3,4,5 CFD provides numerical approximation to the equations that govern fluid motion. Application of the CFD to analyze a fluid problem requires the following steps. First, the mathematical equations describing the fluid flow are written. These are usually a set of partial differential equations. These equations are then discretized to produce a numerical analogue of the equations. The domain is then divided into small grids or elements. Finally, the initial conditions and the boundary conditions of the specific problem are used to solve these equations. The solution method can be direct or iterative. In addition, certain control parameters are used to control the convergence, stability, and accuracy of the method. All CFD codes contain three main elements: (1) A pre-processor, which is used to input the problem geometry, generate the grid, define the flow parameter and the boundary conditions to the code. (2) A flow solver, which is used to solve the governing equations of the flow subject to the conditions provided. There are four different methods used as a flow solver: (i) finite difference method (ii) finite element method, (iii) finite volume method, and (iv) spectral method. (3) A post-processor, which is used to massage the data and show the results in graphical and easy to read format. In this chapter we are mainly concerned with the flow solver part of CFD. This chapter is divided into five sections. In section two of this chapter we review the general governing equations of the flow. In section three we discuss three standard numerical solutions to the partial differential equations describing the flow. In section four we introduce the methods for solving the discrete equations, however, this section is mainly on the finite difference method. And in section five we discuss various grid generation methods and mesh structures. Special problems arising due to the numerical approximation of the flow equations are also discussed and methods to resolve them are introduced in the following sections.

3 2 Mathematical Formulation 2.1 Governing equations The equations governing the fluid motion are the three fundamental principles of mass, momentum, and energy conservation. Continuity 0 ) = + ���t ����� ���.(��V (1) Momentum ��F �� + ���p ��� ���.��ij = Dt DV (2) Energy �� + ���.q ��� ���t ���Q = + p(���.V) Dt De �� (3) where �� is the fluid density, V is the fluid velocity vector, ��ij is the viscous stress tensor, p is pressure, F is the body forces, e is the internal energy, Q is the heat source term, t is time, �� is the dissipation term, and ���.q is the heat loss by conduction. Fourier���s law for heat transfer by conduction can be used to describe q as: ���k���T = q (4) where k is the coefficient of thermal conductivity, and T is the temperature. Depending on the nature of physics governing the fluid motion one or more terms might be negligible. For example, if the fluid is incompressible and the coefficient of viscosity of the fluid, ��, as well as, coefficient of thermal conductivity are constant, the continuity, momentum, and energy equations reduce to the following equations: 0 .V = ��� (5) V ��F �����2 �� + ���p ��� = Dt DV (6) �� + + ���t ���Q = k���2T Dt De �� (7) Presence of each term and their combinations determines the appropriate solution algorithm and the numerical procedure. There are three classifications of partial differential equations6 elliptic, parabolic and hyperbolic. Equations belonging to each of

4 these classifications behave in different ways both physically and numerically. In particular, the direction along which any changes are transmitted is different for the three types. Here we describe each class of partial differential equations through simple examples: Elliptic: Laplace equation is a familiar example of an elliptic type equation. 0 2u = ��� (8) where u is the fluid velocity. The incompressible irrotational flow (potential flow) of a fluid is represented by this type of equation. Another practical example of this equation is the steady state pure heat conduction in a solid, i.e., ���2T=0, as can be readily obtained from equation (7). Parabolic: The unsteady motion of the fluid due to an impulsive acceleration of an infinite flat plate in a viscous incompressible fluid exemplifies a parabolic equation: ���2u = ���t ���u �� (9) where �� is the kinematic viscosity. Transient diffusion equation is represented with a similar equation. In this type of equations, events propagate into the future, and a monotone convergence to steady state is expected. Hyperbolic: Qualitative properties of hyperbolic equations can be explained by a wave equation. 2 2 2 2 u c2 u ���x ��� = ���t ��� (10) where c is the wave speed. In this case, values of solution depend locally on the initial conditions. The propagation signal speed is finite. Continuous boundary and initial values can give rise to discontinuity. Solution is no more continuous and therefore shocks can be observed and captured in this class of equations. Depending on the flow, the governing equations of fluid motion can exhibit all three classifications.

5 2.2 Boundary Conditions The governing equation of fluid motion may result in a solution when the boundary conditions and the initial conditions are specified. The form of the boundary conditions that is required by any partial differential equation depends on the equation itself and the way that it has been discretized. Common boundary conditions are classified either in terms of the numerical values that have to be set or in terms of the physical type of the boundary condition. For steady state problems there are three types of spatial boundary conditions that can be specified: I. Dirichlet boundary condition: (x, z) y, f1 = �� (11) Here the values of the variable �� on the boundary are known constants f1. This allows a simple substitution to be made to fix the boundary value. For example, if u is the flow velocity, its value may be fixed at the boundary of the domain. For instance, for no-slip and no-penetration conditions on the solid walls, the fluid velocity is the same as the velocity of the wall. II. Neuman boundary condition: ( z) y, x, f 2 = ���n ����� (12) Here the derivatives of the variable �� on the boundary are known f2, and this gives an extra equation, which can be used to find the value at the boundary. For example, if the velocity does not change downstream of the flow, we can assume that the derivative of u is zero at that boundary. III. Mixed type boundary condition: ( z) y, x, f3 b a�� = ���n ����� + (13) The physical boundary conditions that are commonly observed in the fluid problems are as follows: (A) Solid walls: Many boundaries within a fluid flow domain will be solid walls, and these can be either stationary or moving walls. If the flow is laminar then the velocity components can be set to be the velocity of the wall. When the flow is turbulent, however, the situation is more complex. (B) Inlets: At an inlet, fluid enters the domain and, therefore, its fluid velocity or pressure, or the mass flow rate may be known. Also, the fluid may have certain characteristics, such as the turbulence characterizes which needs to be specified. (C) Symmetry boundaries: When the flow is symmetrical about some plane there is no flow through the boundary and the derivatives of the variables normal to the boundary are zero.

6 (D) Cyclic or periodic boundaries: These boundaries come in pairs and are used to specify that the flow has the same values of the variables at equivalent positions on both of the boundaries. (E) Pressure Boundary Conditions: The ability to specify a pressure condition at one or more boundaries of a computational region is an important and useful computational tool. Pressure boundaries represent such things as confined reservoirs of fluid, ambient laboratory conditions and applied pressures arising from mechanical devices. Generally, a pressure condition cannot be used at a boundary where velocities are also specified, because velocities are influenced by pressure gradients. The only exception is when pressures are necessary to specify the fluid properties, e.g., density crossing a boundary through an equation of state. There are typically two types of pressure boundary conditions, referred to as static or stagnation pressure conditions. In a static condition the pressure is more or less continuous across the boundary and the velocity at the boundary is assigned a value based on a zero normal-derivative condition across the boundary. In contrast, a stagnation pressure condition assumes stagnation conditions outside the boundary so that the velocity at the boundary is zero. This assumption requires a pressure drop across the boundary for flow to enter the computational region. Since the static pressure condition says nothing about fluid velocities outside the boundary (i.e., other than it is supposed to be the same as the velocity inside the boundary) it is less specific than the stagnation pressure condition. In this sense the stagnation pressure condition is generally more physical and is recommended for most applications. As an example, consider the problem of flow in a section of pipe. If the upstream end of the computational region coincides with the physical entrance to the pipe then a stagnation condition should be used to represent the external ambient conditions as a large reservoir of stationary fluid. On the other hand, if the upstream boundary of the computing region is inside the pipe, and many diameters away from the entrance, then the static pressure condition would be a more reasonable approximation to flow conditions at that location. (F) Outflow Boundary Conditions: In many simulations there is a need to have fluid flow out of one or more boundaries of the computational region. At such "outflow" boundaries there arises the question of what constitutes a good boundary condition. In compressible flows, when the flow speed at the outflow boundary is supersonic, it makes little difference how the boundary conditions are specified since flow disturbances cannot propagate upstream. In low speed and incompressible flows, however, disturbances introduced at an outflow boundary can have an affect on the entire computational region. It is this possibility that is discussed in this article. The simplest and most commonly used outflow condition is that of a ���continuative" boundary. Continuative boundary conditions consist of zero normal derivatives at the

7 boundary for all quantities. The zero-derivative condition is intended to represent a smooth continuation of the flow through the boundary. It must be stressed that the continuative boundary condition has no physical basis, rather it is a mathematical statement that may or may not provide the desired flow behavior. In particular, if flow is observed to enter the computational region across such a boundary, then the computations may be wrong because nothing has been specified about flow conditions existing outside the boundary. As a general rule, a physically meaningful boundary condition, such as a specified pressure condition, should be used at out flow boundaries whenever possible. When a continuative condition is used it should be placed as far from the main flow region as is practical so that any adverse influence on the main flow will be minimal. (G) Opening Boundary Conditions: If the fluid flow crosses the boundary surface in either directions an opening boundary condition needs to be utilized. All of the fluid might flow out of the domain, or into the domain, or a combination of the two might happen. (H) Free Surfaces and Interfaces: If the fluid has a free surface, then the surface tension forces need to be considered. This requires utilization of the Laplace's equation which specifies the surface tension-induced jump in the normal stress p s across the interface: �� �� = s p (14) where �� represents the liquid-air surface tension and �� the total curvature of the interface7,8,9,10. A boundary condition is required at the contact line, the line at which the solid, liquid and gas phases meet. It is this boundary condition which introduces into the model information regarding the wettability of the solid surface. 2.2.1 Example In this example a converging-diverging nozzle with a distributed inlet ports is considered. Inlet mass flow rate is known and flow exits to the ambient air with atmospheric pressure. Figure 1. Schematic of the flow inside and outside of a converging-diverging nozzle Inlet Outlet

8 Choosing the appropriate boundary conditions can reduce the computer effort. In this example the slice shown in Figure 1 is repeated to produce the whole physical domain. Using the periodic boundary condition at the imaginary planes shown in Figure 2 can reduce the computational domain to a much smaller area. Figure 3 shows the other boundary conditions applied to the problem. Figure 2. Minimizing the computational domain using periodic boundary condition Figure 3. Various Boundary Conditions Periodic Boundaries Inflow Symmetric Boundary Solid Wall Open Boundary Outflow

9 3 Techniques for Numerical Discretization In order to solve the governing equations of the fluid motion, first their numerical analogue must be generated. This is done by a process referred to as discretization. In the discretization process, each term within the partial differential equation describing the flow is written in such a manner that the computer can be programmed to calculate. There are various techniques for numerical discretization. Here we will introduce three of the most commonly used techniques, namely: (1) the finite difference method, (2) the finite element method and (3) the finite volume method. Spectral methods are also used in CFD, which will be briefly discussed. 3.1 The Finite Difference Method Finite difference method utilizes the Taylor series expansion to write the derivatives of a variable as the differences between values of the variable at various points in space or time. Utilization of the Taylor series to discretize the derivative of dependent variable, e.g., velocity u, with respect to the independent variable, e.g., special coordinated x, is shown in Figure 4. Consider the curve in Fig. 4 which represent the variation of u with x, i.e., u(x). After discretization, the curve u(x) can be represented by a set of discrete points, ui���s. These discrete points can be related to each other using a Taylor series expansion. Consider two points, (i+1) and (i-1), a small distance ���x from the central point, (i). Thus velocity ui can be expressed in terms of Taylor series expansion about point (i) as: (���x)2 (���x ) ... 6 2 3 3 3u 2 2u + ��� ��� ���i ��� ��� ��� ��� ��� ���x ��� + ��� ��� ��� ��� ��� ��� ��� ��� ���x ��� + ���x ��� ��� ��� ��� ��� ��� ���x ���u + = ui ui+1 (15) and (���x)2 (���x ) ... 6 2 3 3 3 2 2 + ��� ��� ���i ��� ��� ��� ��� ��� ���x ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� ���x ��� + ���x ��� ��� ��� ��� ��� ��� ���x ���u ��� = u u ui ui���1 (16) These equations are mathematically exact if number of terms are infinite and ���x is small. Note that ignoring these terms leads to a source of error in the numerical calculations as the equation for the derivatives is truncated. This error is referred to as the truncation error. For the second order accurate expression, the truncation error is: ( n! u n=3 n n ��� ���x)n���1 ��� ��� ���i ��� ��� ��� ��� ��� ���x ��� ���

10 Figure 4. Location of points for Taylor series By subtracting or adding these two equations, new equations can be found for the first and second derivatives at the central position i. These derivatives are (���x)2 6 2���x 3 3u ui���1 ui+1 ��� ��� ���i ��� ��� ��� ��� ��� ���x ��� ��� ��� = ��� ���i ��� ��� ��� ��� ���x ���u (17) and 2 2u (���x)2 2ui O(���x)2 ui+1 ui���1 + + ��� = ��� ��� ���i ��� ��� ��� ��� ��� ���x ��� (18) Equations (17) and (18) are referred to as the central difference equations for the first and the second derivatives, respectively. Further derivatives can also be formed by considering equations (15) and (16) in isolation. Looking at equation (15), the first-order derivative can be formed as (���x) 2 2 2 u ui ui+1 ��� ��� ���i ��� ��� ��� ��� ��� ���x ��� ��� ���x ��� = ��� ���i ��� ��� ��� ��� ���x ���u (19) This is referred to as the Forward difference. Similarly, from equation (16) another first- order derivative can be formed, i.e., (���x) 2 2 2u ui���1 ui ��� ��� ���i ��� ��� ��� ��� ��� ���x ��� ��� ���x ��� = ��� ���i ��� ��� ��� ��� ���x ���u (20) u ui+1 ui-1 ui xi-���x xi+���x xi x

11 This is referred to as the Backward difference. As noted by the expressions, difference formulae are classified in two ways: (1) by the geometrical relationship of the points, namely, central, forward, and backward differencing or (2) by the accuracy of the expressions, for instance, central difference is second-order accurate, whereas, both forward and backward differences are first-order accurate, as the higher order terms are neglected. We can obtain higher order approximations by applying the Taylor series expansion for more points. For instance, a 3-point cluster would result in a second order approximation for the forward and backward differencing, rather than a first order approximation: Forward difference: ( ) 4ui+1 3ui 2���x 1 O( ui+2 ���x)2 + + + ��� = ��� ���i ��� ��� ��� ��� ���x ���u (21) Backward difference: ( ) 3ui 4ui���1 2���x 1 O( ui���2 ���x)2 + + ��� = ��� ���i ��� ��� ��� ��� ���x ���u (22) Similarly a 4-point cluster results in a third order approximation for the forward and backward differencing: Forward difference: ( ) 6ui+1 3ui 2ui���1 6���x 1 O(3 ui+2 ���x) + ��� ��� ��� ��� = ��� ���i ��� ��� ��� ��� ���x ���u (23) Backward difference: ( ) 2ui+1 3ui 6ui���1 6���x 1 O( ui���2 ���x)3 + + + + = ��� ���i ��� ��� ��� ��� ���x ���u (24) The above difference equations are used to produce the numerical analogue of the partial differential equations describing the flow. In order to apply this discretization method to the whole flow field, many points are placed in the domain to be simulated. Then, at each of these points the derivatives of the flow variables are written in the difference form, relating the values of the variable at each point to its neighboring points. Once this process is applied to all the points in the domain, a set of equations are obtained which are solved numerically. For more discussion on this topic refer to text books on numerical analysis such as Hildebrand11, and Chapra and Canale12. 3.2 The Finite Element Method In the finite element method, the fluid domain under consideration is divided into finite number of sub-domains, known as elements. A simple function is assumed for the variation of each variable inside each element. The summation of variation of the

12 variable in each element is used to describe the whole flow field. Consider the two- nodded element shown in Figure 5, in which variable u varies linearly inside the element. The end points of the element are called the nodes of the element. For a linear variation of u, the first derivative of u with respect to x is simply a constant. If u is assumed to vary linearly inside an element, we cannot define a second derivative for it. Since most fluid problems include second derivative, the following technique is designed to overcome this problem. First, the partial differential equation is multiplied by an unknown function, and then the whole equation can be integrated over the domain in which it applies. Finally the terms that need to have the order of their derivatives reduced are integrated by parts. This is known as producing a variational formulation. Figure 5. A two-noded linear element As an example, we will develop the finite element formulation of the Laplace's Equation in one dimensions: 0 2 2u = dx d (25) where velocity u is a function of the spatial coordinates x. We multiply equation (25) by some function W and integrate it over the domain of interest denoted by ���: 0 2 2 = ��� ��� ��� ��� ��� ��� ���W d��� dx u d (26) Equation (26) can be integrated by parts to result in: xi+1 xi-1 ui ui+1 ui-1 u x

13 0 = ��� ��� ��� ���W ��� ��� + ��� ��� ��� ������ ��� ��� ��� ��� d�� nx dx du d��� dx du dx dW (27) where �� denotes the boundary of the domain ��� and nx is the unit outward normal vector to the boundary ��. The second order derivative in equation (26) is now transformed into products of first order derivatives. Equation (27) is used to produce the discrete form of the partial differential equation for the elements in the domain. Equation (27) is known as the variational form of the partial differential equation (25). Although this technique reduces the order of the derivatives, it introduces the terms corresponding to the boundary of the domain into the governing equation (27). We will now divide the domain into several elements and assume a function for the variation of the variable u in each element. If a two-noded linear element is assumes (see Fig. 5), the variation of u in each element can be represented by ��� ��� ��� ��� ��� ��� ��� ��� + = )��� (ui+1 i xi���1 xi+1 xi���1 xi ui���1 ui���1 u (28) or ��� ��� ��� ��� ��� ��� ��� ��� + ��� ��� ��� ��� ��� ��� ��� ��� = i xi���1 xi+1 xi���1 xi ui+1 xi���1 xi+1 xi xi+1 ui���1 u (29) The terms in the brackets are called the shape functions and are denoted as Ni���s. ui-1 and ui+1 are the nodal values of the variable u and are denoted as ui���s. Therefore, the variable u can be written in the following form + = i+1 i ui+1 N Ni���1ui���1 u (30) Thus, the shape functions corresponding to the two-nodal linear element, represented by equation (28) are ��� ��� = i���1 xi���1 xi+1 xi xi+1 N (31) and ��� ��� = i+1 xi���1 xi+1 xi���1 xi N (32) We can now determine the derivatives of the variable u, using equation (30): m i=1 ui dx dNi dx du ��� = (33)

14 where m is the number of nodes on the element. Note that ui���s are nodal values of u and they are not variables, therefore, they are not differentiated. In order to solve equation (27) we still need to describe the function W. There are several methods, which are used for the specification of the variable W. However, the most common method is the Galerkin method in which W is assumed to be the same as the shape function for each element. Therefore, equation (27) is discretized by using equations similar to equation (33) for the derivatives of the variable, and equations similar to equations (31) and (32) for W. For every element there can be several equations depending on the number of the nodes in that element. The set of equations generated in this form are then solved together to find the solution. The above formulation was based on a linear variation of the variable in each element. If higher order variations are used, e.g., quadratic or cubic, second derivatives will appear which require more points to describer them. This makes the computation more cumbersome. References [3] and [4] are recommended for more detailed discussion of FEM. 3.3 The Finite Volume Method The finite volume method is currently the most popular method in CFD. The main reason is that it can resolve some of the difficulties that the other two methods have. Generally, the finite volume method is a special case of finite element, when the function W is equal to 1 everywhere in the domain. This technique is discussed in detail by Patankar.34 A typical finite volume, or cell, is shown in Fig. 6. In this figure the centroid of the volume, point P, is the reference point at which we want to discretize the partial differential equation. Figure 6. A finite volume in one dimension. The neighboring volumes are denoted as, W, volume to the west side, and E, the volume to the east side of the volume P. For the one-dimensional finite volume shown in Fig. 6, the volume with centroid P, has two boundary faces at w and e. W P E w e