Mechanical Engineering Department
ME 692 – Computational Fluid Dynamics
Spring 2002 Ticket: 57541 Instructor: Larry Caretto
Engineering Building Room 2303Mail CodePhone: 818-677-6448
E-mail: 8348Fax: 818-677-7062
Basics of CFDME 692, L. S. Caretto, Spring 2002Page 1
Three – Basics of Computational Fluid Dynamics
Using the Basic Equations
In the first part of these notes we showed how the basic equations of fluid dynamics were derived from balances of mass, species, momentum, and energy. The derivations involve two parts. The first is a balance that involves terms such as the fluid stresses, σyx or σji, the heat flux vector, qx or qi, and the diffusion flux of species K, ji(K). We later used various relations between these terms, which will call “flux terms” as a general reference in these notes, and other flow properties: Newton’s law for fluid stress, Fourier’s law for heat conduction, and Fick’s law for diffusion. Relations of this kind are called constitutive relations that relate flux quantities to other fluid properties, usually by gradient terms.
With this system in place, we showed that the basic balance equations for computational fluid dynamics could be written in the form given by equation [1-81], which is copied below.
[1-81]
In this equation c = 1 in all equations except for the temperature-based energy equations where c may represent cP or cv; Γ(φ) represents some transport property such as the viscosity or thermal conductivity. The “Source” term is a combination of true source terms, such as the production of species by chemical reactions and other terms that do not fit into the classification of transient, convection and diffusive. Table 1-2 shows how this general representation is used for the different balance equations that we will be using.
Equation [1-81] has both first and second order derivatives. In a simple case of a one-dimensional flow, equation [1-81] becomes
[3-1]
For constant property flow, with no body forces, the only part of the source term is the pressure gradient. We can move the properties outside the derivatives, divide by the density, and use the definition of the kinematic viscosity, = μ/ρ, to rewrite equation [3-1] as follows.
[3-2]
Some approaches start with equations like [3-2] to obtain finite-difference or finite-element forms. Other approaches, particularly in compressible flows and some finite-element approaches prefer to work with first order equations. In this approach the various flux terms and their constitutive relations are treated as separate equations in the numerical algorithms. In these approaches the total stress term, σij, is written as the sum of the pressure, p δij, and the viscous stress, ij.
[3-3]
Comparing this equation with equation [1-29] for σij shows that the viscous stress term is given by the following equation.
[3-4]
Over overall momentum balance, in of the stress, was given by equation [1-28], which is copied below.
[1-28]
Using equation [3-3] we can rewrite this general momentum balance equation as follows
[3-5]
Recall that the summation convention over repeated indices means that the term is actually the sum of three different terms. (When equation [3-3] was substituted into equation [1-28] there was a term . However, this term may be written as, since δij = 0 if i is not equal to j.)
Equation [3-5] uses only first order derivatives, but it is necessary to use equation [3-4] to compute the viscous stress terms, ij, in this equation. (Because the viscous stress terms are symmetric – i.e., ij = ji – it is in only necessary to compute six of the nine different ij terms.)
The solution of CFD equations, using differential equations that have only first derivatives, also uses the following equations for mass, energy, and species balance, copied from the previous notes on deriving equations.
[1-18]
[1-41]
[1-75]
We can use equation [3-3] to replace σij, in the energy balance equation by the pressure and the viscous stress, ij. Doing this and using the fact that δij = 0 if i is not equal to j, allows us to rewrite the energy balance equation as follows:
[3-6]
Since the subscripts are used as summation indices, we can replace i by j on the convection term and the heat flux gradient term without changing the result. Doing this allows us to rewrite equation [3-6] as follows
[3-7]
CFD practitioners who use the first-derivative approach derive a general algorithm from writing a vector form of the equation. This form summarizes the continuity equation [1-18], the three momentum equations implied in equation [3-5], the energy equation of equation [3-7] and the species balance of equation [1-75]. The equation they use is written as follows.
[3-8]
Where the terms U, E, F, G, and H are vectors (one-dimensional column matrices) defined by the following equations:
[3-9]
[3-10]
[3-11]
[3-12]
The algorithms that use this format typically solve for the U matrix components. (In equations [3-9] and [3-13], the ui terms refer to the components of the U vector, not to velocity components.) Then the various flow properties are found from the following equations:
[3-13]
The use of these equations in CFD algorithms, which we will discuss later, is similar to the process outlined in the first section of these course notes. The equations for all the variables that we have to use have a similar form. This means that a numerical algorithm that we can apply a numerical algorithm, devised for one physical equation, all the phenomena that we seek to model. In the equations above we have talked about six simultaneous equations; in a multicomponent system we would have one equation for each chemical species in the flow. Thus if we had n chemical species, we would have 5+n equations to solve. The unknowns, as shown above, are the density, three velocity components, the internal energy and the mass fractions of the n chemical species present. We also have the pressure as an unknown. In a compressible flow problem we would have an equation of state, P = P(T,ρ) that could be used to solve for the pressure.
Finite-difference methods and stability
In section two of the notes we applied finite-difference and finite-element methods to a simple ordinary differential equation. The extension of the finite-difference approach used there to partial differential equations is fairly straightforward. As an example of this, consider the one-dimensional equation for heat conduction in a solid with constant properties and no source term.
[3-14]
If you are not familiar with this equation, you can derive it from equation [1-71], shown below, by setting all velocities, dissipation and dilatation to zero. The constant property assumption allows the properties to be brought outside the differential operators. For the one-dimensional case, we have spatial derivatives only with respect to x1 or, more simply, x.
[1-71]
In the conduction equation, we define the thermal diffusivity, = k/(ρ cv) and write equation [3-14] with the single property value, .
[3-15]
Equation [3-15] has an open boundary in time; we do not have to specify the conditions at a future time. We need an initial condition for the temperature at all points in the region. We can write a general initial condition as T(x,0) = T0(x). Similarly we can write arbitrary boundary conditions for the boundaries at xmin and xmax: T(xmin,t) = TL(t) and T(xmax,t) = TR(t). These conditions, as well as values for the geometry and the thermal diffusivity must be specified before we can solve the differential equation.
We can construct finite-difference grids in space and time. For simplicity, we will assume constant step sizes, Δt and Δx. We define our time and space coordinates on this grid by the integers n and i, so that we can write
[3-16]
We can obtain a solution process, known as the explicit method, if we use a forward difference for the time derivative and a central difference for the space derivatives. These are given by the following equations.
[3-17]
These equations are modifications of equations [2-19] and [2-29] from the introductory notes on numerical methods. In those notes we dealt with ordinary derivatives and needed only one subscript for the dependent grid variable. Here, in a multidimensional case, we have temperature as a function of time and distance, T(x,t). Thus we define Tin = T(xi, tn). We use differences between variables like Tin, T n i+1, and T n i-1 to give differences in the x-direction that are used to compute spatial derivatives. The constant n superscript in the space derivative is an indication that we are holding time constant and varying distance. Similarly, the time derivatives are reflected in the finite-difference form by using variations in the superscript n while holding the subscript i constant.
The forward difference for the time derivative is chosen over the more accurate central difference to provide a simple algorithm for solving this problem. There are other problems with the use of the central difference expression for the time derivative that we will discuss later.
If we substitute the finite-difference expressions in [3-17] into equation [3-15], we get the following result.
[3-18]
We have written the combined error term to show the order of the error for both Δx and Δt. We have not written these as separate terms since O() notation simply shows how the error depends on the power of the step size. If we neglect the truncation error, we can obtain a simple equation for the value of the temperature at the new time step.
[3-19]
We see that the grid spacings and the thermal diffusivity appear in a single term. We can define this as a single term to save writing.
[3-20]
With this definition, equation [3-19] may be written as follows.
[3-21]
With this finite difference equation, each temperature at the new time step is related only to temperatures at the previous time step.[*] These are the temperature at the same spatial location and the two neighboring temperatures in the grid. Equation [3-21] is illustrated below.
times f / + / times (1 – 2f) / + / times f=
From the initial conditions, we know all the temperatures at the first time step. That is, Ti0 equals the given initial temperature at x, T0(xi). Similarly, the temperatures at i = 0, and i = imax are known from the boundary conditions. We can illustrate this problem for a simple case where the initial conditions and the boundary conditions are constants. We define the boundaries of our one-dimensional region as x = 0 and x = L. The initial temperature, denoted as Tinit, is the same for all values of x. The constant boundary temperature at x = 0, for all times greater than zero, will be called T0, and the boundary temperature at x = L will be denoted as TL. We can compare our numerical results for this problem with the analytical solution for both the temperature and the temperature gradient. The analyticalsolution to equation [3-15] for the constant boundary conditions used here, is given by the following equations.[*]
[3-22]
This analytical solution can be used to obtain the temperature and the temperature gradient, which is proportional to the heat flux, at any time and space coordinate, for any given values of , L, Tinit, T0, and TL.
In order to obtain a numerical solution, which is only obtained at discrete points, we must specify numerical values of , L, Tinit, T0, and TL. For simplicity let us choose = 1, L = 1, Tinit = 1000 and T0 = TL = 0. We will start with a small computational grid using only five points in the x direction. This gives Nx = 4 and Δx = 0.25. We also pick a time step, Δt = 0.01 and calculate to a maximum time of 0.2. This requires 20 time steps.
Table 3-1 below shows the results of the numerical calculations with the spatial variation in columns and the time variation in rows. There are two initial condition rows. The “t = 0” row shows the initial conditions. The “t = 0+” row shows the temperatures immediately after t = 0 when the boundary conditions are applied. This second row is the starting time point, n = 0, for the numerical calculations. These two initial condition rows illustrate the discontinuity in the temperature at the two points (0,0) and (L,0) in the solution domain. This discontinuity does not affect our ability to obtain an analytical solution. We must treat it explicitly in the numerical solution.
In addition to the modified initial condition row, where we know all the temperatures, we have the boundary condition columns. For all time points at x = 0 and x = L = 1, we set the table entries to the known boundary temperatures. Although the initial and boundary conditions are particularly simple in this example, we could handle initial conditions that were arbitrary functions of x and boundary conditions that were arbitrary conditions of time by placing the appropriate results for the initial and boundary temperatures at the correct places in the table.
Once the modified initial condition row and the boundary condition columns are set, the solution for the remaining temperatures, using equation [3-21] is straightforward. We first use equation [3-20] to compute the value of f.
[3-23]
Once the value of f is known, we can apply equation [3-21] to each node, in successive rows. The calculations below are for the first time step.
[3-24]
These are the results shown in the “n = 1” row of the table. The numerical results are symmetric in x as would be expected from the uniform initial condition and the symmetric boundary condition. The process illustrated in equation [3-24] is used to calculate the remaining results in the table. The final two entries in the table are the exact results, calculated for the different values of x at the final time t = 0.2, and the error. The latter is the absolute value of the difference between the numerical and analytical solution.
Table 3-1
Numerical solution of the Conduction Equation
= 1, Δx = 0.25, Δt = 0.01, tmax = 0.20i = 0 / i = 1 / i = 2 / i = 3 / i = 4
x = 0.00 / x = 0.25 / x = 0.50 / x = 0.75 / x = 1.00
t = 0 / 1000 / 1000 / 1000 / 1000 / 1000
n = 0 / t = 0+ / 0 / 1000 / 1000 / 1000 / 0
n = 1 / t = 0.01 / 0 / 840 / 1000 / 840 / 0
n = 2 / t = 0.02 / 0 / 731.2 / 948.8 / 731.2 / 0
n = 3 / t = 0.03 / 0 / 649.0 / 879.2 / 649.0 / 0
n = 4 / t = 0.04 / 0 / 582.0 / 805.5 / 582.0 / 0
n = 5 / t = 0.05 / 0 / 524.6 / 734.0 / 524.6 / 0
n = 6 / t = 0.06 / 0 / 474.2 / 667.0 / 474.2 / 0
n = 7 / t = 0.07 / 0 / 429.2 / 605.3 / 429.2 / 0
n = 8 / t = 0.08 / 0 / 388.7 / 548.9 / 388.7 / 0
n = 9 / t = 0.09 / 0 / 352.1 / 497.7 / 352.1 / 0
n = 10 / t = 0.10 / 0 / 319.1 / 451.1 / 319.1 / 0
n = 11 / t = 0.11 / 0 / 289.2 / 408.9 / 289.2 / 0
n = 12 / t = 0.12 / 0 / 262.0 / 370.5 / 262.0 / 0
n = 13 / t = 0.13 / 0 / 237.5 / 335.8 / 237.5 / 0
n = 14 / t = 0.14 / 0 / 215.2 / 304.4 / 215.2 / 0
n = 15 / t = 0.15 / 0 / 195.0 / 275.8 / 195.0 / 0
n = 16 / t = 0.16 / 0 / 176.8 / 250.0 / 176.8 / 0
n = 17 / t = 0.17 / 0 / 160.2 / 226.5 / 160.2 / 0
n = 18 / t = 0.18 / 0 / 145.2 / 205.3 / 145.2 / 0
n = 19 / t = 0.19 / 0 / 131.6 / 186.1 / 131.6 / 0
n = 20 / t = 0.20 / 0 / 119.2 / 168.6 / 119.2 / 0
Exact / t = 0.20 / 0 / 125.1 / 176.9 / 125.1 / 0
Error / t = 0.20 / 0 / 5.8 / 8.2 / 5.8 / 0
The numerical calculation of the final temperatures for an end time of = 0.2 can be repeated for time steps of 0.02 and 0.04. These results are shown in the Tables 3-2 and 3-3 below. When the time step is doubled from 0.01 to 0.02 the solution results are similar. However, the error in the temperatures at the final time step has increased by about a factor of three.
When the time step is increases to 0.04, the results are completely different, and these results do not match our ideas of physical reality. We expect that the region will eventually reach a constant temperature equal to the two boundary temperatures of zero. As this cooling takes place, we do not expect to see any temperatures outside the range bounded by the initial temperature of 1000 and the boundary temperature of 0. However, the solutions with a time step of 0.04 have some negative temperatures in the solution. You can verify that these unexpected negative temperatures are correctly computed from the finite difference equation.
Table 3-2
Numerical solution of the Conduction Equation
= 1, Δx = 0.25, Δt = 0.02, tmax = 0.20i = 0 / i = 1 / i = 2 / i = 3 / i = 4
x = 0.00 / x = 0.25 / x = 0.50 / x = 0.75 / x = 1.00
t = 0 / 1000 / 1000 / 1000 / 1000 / 1000
n = 0 / t = 0+ / 0 / 1000 / 1000 / 1000 / 0
n = 1 / t = 0.02 / 0 / 680 / 1000 / 680 / 0
n = 2 / t = 0.04 / 0 / 564.8 / 795.2 / 564.8 / 0
n = 3 / t = 0.06 / 0 / 457.9 / 647.7 / 457.9 / 0
n = 4 / t = 0.08 / 0 / 372.1 / 526.2 / 372.1 / 0
n = 5 / t = 0.10 / 0 / 302.3 / 427.6 / 302.3 / 0
n = 6 / t = 0.12 / 0 / 245.7 / 347.4 / 245.7 / 0
n = 7 / t = 0.14 / 0 / 199.6 / 282.3 / 199.6 / 0
n = 8 / t = 0.16 / 0 / 162.2 / 229.4 / 162.2 / 0
n = 9 / t = 0.18 / 0 / 131.8 / 186.4 / 131.8 / 0
n = 10 / t = 0.20 / 0 / 107.1 / 151.4 / 107.1 / 0
Exact / t = 0.20 / 0 / 125.1 / 176.9 / 125.1 / 0
Error / t = 0.20 / 0 / 18.0 / 25.4 / 18.0 / 0
Table 3-3
Numerical solution of the Conduction Equation
= 1, Δx = 0.25, Δt = 0.04, tmax = 0.20i = 0 / i = 1 / i = 2 / i = 3 / i = 4
x = 0.00 / x = 0.25 / x = 0.50 / x = 0.75 / x = 1.00
t = 0 / 1000 / 1000 / 1000 / 1000 / 1000
n = 0 / t = 0+ / 0 / 1000 / 1000 / 1000 / 0
n = 1 / t = 0.04 / 0 / 360 / 1000 / 360 / 0
n = 2 / t = 0.08 / 0 / 539.2 / 180.8 / 539.2 / 0
n = 3 / t = 0.12 / 0 / -35.3 / 639.6 / -35.3 / 0
n = 4 / t = 0.16 / 0 / 419.2 / -224.2 / 419.2 / 0
n = 5 / t = 0.20 / 0 / -260.9 / 599.3 / -260.9 / 0
Exact / t = 0.20 / 0 / 125.1 / 176.9 / 125.1 / 0
Error / t = 0.20 / 0 / 385.9 / 422.5 / 385.9 / 0
There is something obviously wrong with these results, but there is no obvious problem. We have used the same finite difference equation that had previously produced correct results. The problem that we are facing here results from the translation of the differential equation into a finite difference equation. It is an issue known as stability. Numerical algorithms that are not stable can produce erroneous results whose errors grow without bound.
We can do a heuristic analysis of stability for the algorithm considered here. This is a heuristic analysis because it uses physical reasoning rather than examining the error growth of the finite difference equations. This argument proceeds by considering the physical implications of the finite-difference equation, [3-21].
[3-21]
This equation shows that the temperature at the new time step, for any given value of xi depends on the value of the temperature at the previous time step, at the same value of xi, multiplied by the factor (1 – 2f). Depending on the value of f, the (1 – 2f) factor may be positive, negative or zero. If the factor is negative, we have a physically unrealistic result. For given values of the neighboring temperatures, the new temperature at xi will be smaller, if the temperature value at the previous time step is larger. This is not correct. We expect that larger values of temperature at xi at one time step should lead to larger values of temperature at xi at the next time step. Thus, we do not want to use negative values for the factor (1 – 2f). This sets the following limit on our numerical calculations.
[3-25]
In the three sets of calculations above the value of f was 0.16, 0.32, and 0.64. Thus, the third set of calculations violated the stability limit. Because we used an inadmissible value of f, we obtained incorrect results from our calculation. The stability limit becomes a more serious concern when we are dealing with very small values of Δx. We would have to choose a very small value of Δt to satisfy the stability condition.
The stability limit on Δt differs from the need to choose a small Δt for accuracy. The stability limit places a significant restriction on Δt. If we were to half the step size, Δx, to obtain a more accurate solution, we would have to cut Δt by ¼. Thus, doubling the number of space steps requires a factor of four increase in the time steps or an overall factor of eight in the total nodes in the calculation. This requirement for stability is separate from the requirement for accuracy. Accuracy depends on small sizes for Δt and Δx. We may want to use a larger step size for Δt to obtain a desired accuracy, however we cannot use a step size that violates the stability condition in equation [3-25] if use the explicit method.
The appendix to this section of notes contains a derivation for the truncation error for this approach. There are two equivalent equations for this truncation error. Both of these are written as an infinite series in equation [3-26].
[3-26]
We see that the truncation error depends in a complex way on the two finite-difference grid spacings, Δt and Δx. Choosing f = 1/6 eliminates the first (k = 2) term in the truncation error and gives a higher order truncation error.
The method just described is called the explicit method because all temperatures at the new time step can be written as an explicit function of temperatures at the old time step. Alternative, implicit approaches use the value of the temperature at the new time step. The Crank-Nicholson method is the most popular implicit method. This method uses a finite-difference representation of the conduction equation at a time point midway between the two specified time grid lines. Although there is no computational point there, we can write this time derivative as follows, using a central-difference approximation. In this expression the step size is Δt/2, the distance from the location of the derivative to the location of each of the two points used in its computation.