It is common practice to highlight the contributions due to the surface forces as separate terms in the momentum equation and to include the effects of body forces as source terms. The state of stress of a fluid element is defined in terms of the pressure and the nine viscous stress components shown in adjacent figure. The pressure, a normal stress, is denoted byp. Viscous stresses are denoted by t. The usual suffix notation t, is applied to indicate the direction of the viscous stresses. The suffices i and j in t, indicate that the stress component acts in the j direction ‘on @ surface normal to the i-direction.It ls common practice to highlight the contributions due to the surface forces @5 separate terms in the momentum equation and to include the effects of body forces as source terms. The state of stress of a fluid element is defined in terms of the pressure and the Viscous stresses are denoted by t. The usual suffix notation ©, is applied to indicate the direction of the viscous stresses. The suffices | and j in t, indicate that the stress component acts in the j direction on a surface normal to the i-direction.First we consider the x-components of the forces due to pressure p and stress components t,,, t,, and t,, shown in the figure below. The magnitude of a force resulting from a surface stress is the product of stress and area. Forces aligned with the direction of a co-ordinate axis get a positive sign and those in the opposite direction a negative sign. The net force in the x-direction is the sum of the force components acting in that direction on the fluid element.On the pair of faces (E, W) we have [(-33)- (.-% Se lfeee [di i*|* {«. “Sye)/oe- (-E +3 Jaws: The net force in the x-direction on the pair of faces (N, S) is Senne fede entra Finally the net force in the x-direction on faces T and B is given by -{ t= Seon Sete a ate The total force per unit volume on the fluid due to these surface stresses is equal to the sum of above three equations divided by the volume 6x5y5z:On the pair of faces (E, W) we have [(--33+)- (.-= Soe [fe i] {(x.+% Mal Frs}oe- [-2 + |aws: The net force in the x-direction on the pair of faces (N, S) is ~{ santos.» eto ate ese Finally the net force in the x-direction on faces T and B is given by {= Se Jae (+ Sela Sone The total force per unit volume on the fluid due to these surface stresses is equal to the sum of above three equations divided by the volume 5x5y5z: EP + te) , He , Me wyWithout considering the body forces in further detail their overall effect can be included by defining a source S,,, of x-momentum per unit volume per unit time. The x-component of the momentum equation is found by setting the rate of change of x-momentum of the fluid particle (equation 2 slide 16) equal to the total force in the x-direction on the element due to surface stresses (equation 4 slide 19) plus the rate of increase of x-momentum due to sources: Du _ H-P+ Ta) Nye, Mee "Oo eee ae It is not too difficult to verify that the y-component of the momentum equation is given by De ak + Se and the z-component of the momentum equation by ple Ata , Pn Dt o& a De _ Hy , HP+%) a. a Op + te) Se Sg +e 7Without considering the body forces in further detail their overall effect can be included by defining a source S,,, of x-momentum per unit volume per unit time. The x-component of the momentum equation is found by setting the rate of change of x-momentum of the fluid particle (equation 2 slide 16) equal to the total force in the x-direction on the element due to surface stresses (equation 4 slide 19) plus the rate of increase of x-momentum due to sources: Du _ HP + Tae), Hye , Mer, fate wk ae It is not too difficult to verify that the y-component of the momentum equation is given by De My , APT) , Me Dk a dz and the z-component of the momentum equation by pee a, Hee, HD+ Ea) =—+ Dh wk wy oe* The sign associated with the pressure is opposite to that associated with the normal viscous stress, because the usual sign convention takes a tensile stress to be the positive normal stress so that the pressure, which is by definition a compressive normal stress, has a minus sign. * The effects of surface stresses are accounted for explicitly; the source terms Sia, Say and Sy, im (equation 1 - 3 of slide 20) include contributions due to body forces only. * For example, the body force due to gravity would be modelled by 5.., = 0, S.4,=0 and Sin, = -pg. Energy equation in three dimensions * The energy equation is derived from the first law of thermodynamics, which States that the rate of change of energy of a fluid particle is equal to the rate of heat addition to the fluid particle plus the rate of work done on the particle: Rate of increase Net rate of Net rate of work of energy of = heat addedto + done on. As before, we will be deriving an equation for the rate of increase of energy of a fluid particle per unit volume, which is given by DE De Work done by surface forces * The rate of work done on the fluid particle in the element by a surface force is ‘equal to the product of the force and velocity component in the direction of the force. * For example, the forces given by (Eq 1—300n slide 19) all act in the x-direction. * The work done by these forces is given by [[~-2grys}- (niet ~ [me Me) (nA Jom [(o=-A5P) [soe] ae | feo Aaa) (ASThe net rate of work done by these surface forces acting in the x-direction is given by Auk —p + T,.)) Hut) , Aut.) el im én Pa Slane Surface stress components in the y- and z-direction also do work on the fluid Particle. A repetition of the above process gives the additional rates of work done on the fluid particle due to the work done by these surface forces: [geste Hvty) et ox % é The total rate of work done per unit volume on the fluid particle by all the surface forces is given by adding above equations divided by the volume 5x5y6z. The terms containing pressure can be collected together and written more compactly in vector form -Xeth _ Hep) _ He) _sivigay x: ey & ie (22 Hut) , Aalap + paEnergy flux due to heat conduction * The heat flux vector q has three components: g,, q, and q, The net rate of heat transfer to the fluid particle due to heat flow in e the x-direction is Biven by the ‘-F difference between the rate of heat input across face W and the rate of heat loss across face E: a | dy6 M Gu dyée | ax HA {, , MlFourier’s law of hi Bradient. So ay This can be written in vector form as follows: q= hk grad T Combining equation 2 and 4 of this slide yields the final form of the rate of heat addition to the fluid particle due to heat conduction across element boundaries: ~div q= div( grad 7)Thus far we have not defined the specific energy E of a fluid. Often the energy of a fluid is defined as the sum of internal (thermal) energy i, kinetic energy 1/2 (u + v? + w*) and gravitational potential energy. This definition takes the view that the fluid element is storing gravitational potential energy. It is also possible to regard the gravitational force as a body force, which does work on the fluid element as it moves through the gravity field. Here we will take the latter view and include the effects of potential energy changes as a source term. As before, we define a source of energy S; per unit volume per unit time. Conservation of energy of the fluid particle is ensured by equating the rate of change of energy of the fluid particle (eq 1 slide 22) to the sum of the net rate of work done on the fluid particle (eq 1 slide 24), the net rate of heat addition to the fluid (last eq of slide 25) and the rate of increase of energy due to sources.the above equation is a per practice to extract the chang quation for int of th rnal energy | or e pi energy equation attributable to th by multiplying the x-momentum equation by velocit momentum equation by v gether. e kine emperature T. nponent u ind the z-momentum equation bi y* Subtracting this equation from previous equation and definin term as S,=Se-u.Siy Ids the internal energy equation as liv u + div(é grad T) + t r :bi Ou Ou du ae ce T) +t, a +h 5 “ay? me at "5, Cis Ow Ow a 3 a ae + For compressible flows (eq 1 on slide 27) is often rearranged to give an equation for the enthalpy. * The specific enthalpy h and the specific total enthalpy ho of a fluid are defined as: h=itp/p and hy=h+toe+e+ wy‘specitic heat and divu =U. This allows us to recast the previous equation into a temperature equation as: DT Ou Ou du ad — = div(k grad 7) + t,,— + T,— + F. Om i ox oy dz ax dv For compressible flows (eq 1 on slide 7) is often rearranged to give an equation for the enthalpy. The specific enthalpy h and the specific total enthalpy hy of a fluid are defined as: h=it+p/p and hg=hot(u'+ +o) Combining these two definitions with the one for specific energy E we get hy =i + p/p +4(e+ 2+ w)= E+ p/pEquations of state + The motion of a fluid in three dimensions is described by a system of five partial differential equations: mass conservation, x-, y- and z-momentum equations and energy equation. * Among the unknowns are four thermodynamic variables: p, p, | and T. + Relationships between the thermodynamic variables can be obtained through the assumption of thermodynamic equilibrium. + The fluid velocities may be large, but they are usually small enough that, even though properties of a fluid particle change rapidly from place to place. the fluid can thermodynamically adjust itself to new conditions so quickly that the changes are effectively Instantaneous. * Thus the fluid always remains in thermodynamic equilibrium. * The only exceptions are certain flows with strong shockwaves, but even some of those are often well enough approximated by equilibrium assumptions. + We can describe the state of a substance In thermodynamic equilibrium by means of just two state variables. * Equations of state relate the other variables to the two state variables. © Ifwe use p and T as state variables we have state equations for pressure p and specific internal energy I: p=pp,T) and t= Kp,T)For a perfect gas the following, well-known, equations of state are useful: p=pRT and 1=CT * The assumption of thermodynamic equilibrium eliminates all but the two thermodynamic state variables, + In the flow of compressible fluids the equations of state provide the linkage between the energy equation on the one hand and mass conservation and momentum equations on the other. * This linkage arises through the possibility of density variations as a result of pressure and temperature variations in the flow field. * Uquids and gases flowing at low speeds behave as incompressible fluids. * Without density variations there is no linkage between the energy equation and the mass conservation and momentum equations. * The flow field can often be solved by considering mass conservation and momentum equations only. * The energy equation only needs to be solved alongside the others if the problem involves heat transfer.Sey = Sy = There are also six shearing linear deformation components: The volumetric deformation is given by du de Ow. t+ eadivu ax dy de * Ina Newtonian fluid the viscous stresses are proportional to the rates of deformation. * The three-dimensional form of Newton's law of viscosity for compressible flows involves two constants of proportionality: the first (dynamic) viscosity, 1, to relate stresses to linear deformations, and the second viscosity, A, to relate stresses to the volumetric deformation. * The nine viscous stress components, of which six are independent, are: Ms Wa. a ty =u + Adiv u t= + Adiv u : es Adv u ax ay LsThere are also six shearing linear deformation components: tee oe, 2 a(n oe ” aay Or "Ol ds The volumetric deformation is given by Ou do, Ow oer ox dy a + Ina Newtonian fluid the viscous stresses are proportional to the rates of deformation. * The three-dimensional form of Newton's law of viscosity for compressible flows involves two constants of proportionality: the first (dynamic) viscosity, U1, to relate stresses to linear deformations, and the second viscosity, A, to relate stresses to the volumetric deformation. * The nine viscous stress components, of which six are independent, are: =divu t= ts Adiv u Ty = 112 + A div u ax a ou Ow ul—+ Os ox Maa +Adivu aNot much is known about the second viscosity A, because its effect is small in practice. For gases a good working approximation can be obtained by taking the value of A= (-2/3) w Liquids are incompressible so the mass conservation equation is div u =0 and the viscous stresses are just twice the local rate of linear deformation times the dynamic viscosity. Substitution of these shear stresses into momentum equation (eq 1—3 of slide 20) yields the so-called Navier-Stokes equations De) a alias See allt (ovice\l tol (ex on =- dh + + Suv Or are peda ‘| aye llaywaril| wosl(tlescor|| " Dv ap al (du. av m al (dv. aw S ad s Dy ar oeaeslb |" eee “+3 Elias aay || bv, au, w)] , af (av, au)], af, =-— =—+— =| a] —+— U— + Adi +58, La we BoE) Sloe ay || aa ue ce al aeOften it is useful to rearrange the viscous stress terms as follows: af, au af (au ae)), af (au, dm PY py a raivale 2 (4% ae alu sls |: ls = (a), afm) a( a) [af a) a( a0) a( a) a. Fes} sles} zl) |B) aes) a } Aes »| = div(fe grad u) + [ssid * The viscous stresses in the y- and z-component equations can be recast in a similar manner. * We clearly intend to simplify the momentum equations by ‘hiding’ the bracketed smaller contributions to the viscous stress terms in the momentum source. * Defining a new source by: Su=Su+lsul The Navier-Stokes equations can be written in the most useful form for the development of the finite volume method:Du a =~ em + div(u grad u) + Sy, Dv 9 pe 7 + div(j1 prad v) + Sy, Dw op pP— = —-— + div(gi grad w) + § ee “ii : 'f we use the Newtonian model for viscous stresses in the internal energy equation we obtain after some rearrangement: Di ae ~p div u + div(k grad T) + +, All the effects due to viscous stresses in this internal energy equation are described by the dissipation function ©, which, after considerable algebra, can be shown to be equal to (3) (2) +(2} au =| [> dw)’ ede} | +|—+—1 + tml dz de ae %)| on >*s + Adis u)Conservative form of the governing equations of fluid flow * To summarize the findings thus far, following table gives the conservative or divergence form of the system of equations which governs the time dependent three-dimensional fluid flow and heat transfer of a compressible Newtonian fluid. ® dinamo Ae sinipua) = di grads) +S, 2 aig) = dna rad 9) +S 2 snipouy= 2 sda grad n+ | 2 ipa = = —pdivu+div(k grad 7) ++ S, P= Kp, T and += ip, T) ©, perfect gas p= pRT and i= CT ml & EDifferential and integral forms of the general transport equations * Itis clear from the previous table that there are significant commonalities between the various equations. * If we introduce a general variable o the conservative form of all fluid flow equations, including equations for scalar quantities such as temperature and pollutant concentration etc., can usefully be written in the following form eum + div(pou) = div(T grad a) + Sy Mt Rate of increase Net rate of flow Rate of increase Rate of increase of offlud — + of Oout of = of odue to tot odue to element fluid element diffusion: sources. * The above equation is also called transport equation for property >. * Itclearly highlights the various transport processes: the rate of change term and the convective term on the left hand side and the diffusive term (I = diffusion coefficient) and the source term respectively on the right hand side. * Above equation is also used as the starting point for computational procedures in the finite volume method. * By setting equal to 1, u, v, w and i (or T or hy) and selecting appropriate values for diffusion coefficient [ and sourte terms, we obtain special forms of table shown in previous slide for each of the five PDEs for mass, momentum and energy conservation.‘The key step of the finite volume method, i: the mtegrabon of equaton 1 of Previous slide over 2 three-dimensional contro! volume (CV) Age ¥+ | sow = | div grad ond | Sat iaeeteantuiantisin hal eebd tiedies: Sled Nid te, on onsets term, and in the first term on the right hand side, the diffusive term, are Fewritten a5 integrals over the entire bounding surface of the control volume by weing Gauss's divergence theorem. For a vector @ this theorem states [asa oats ‘The physical interpretation of n.a is the component of vector a in the direction of the vector n normal to surtace element dA. ‘Thus the integral of the divergence of a vector @ over @ volume Is equal to the Component of a in the direction normal to the surface which bounds the vokime ‘summed (integrated) over the entire bounding wrface A.Similarly, the product n . (T grad $), which is also equal to [(-n . (-grad $)), can be interpreted as a positive diffusion flux in the direction of the inward normal vector -n, ie. into the fluid element. The first term on the right hand side of equation 1 of previous slide, the diffusive term, is thus associated with a flux into the element and represents the net rate of increase of fluid property > of the fluid element due to diffusion. The final term on the right hand side of this equation gives the rate of increase of property @ as a result of sources inside the fluid element. In words, this relationship can be expressed as follows: Net rate of decrease Net rate of 7 Rate of increase of 6 due to increase uf @ ee of @inside the + convection acrass = duc to diffusion on control volume the control volume across the control toundecies aeApplying Gauss’s divergence theorem, equation 1 of previous slide can be written as follows: 5 | e + [a (peu d= | in (msgs | sar 4 a a The order of integration and differentiation has been changed in the first term On the left hand side of above equation to illustrate its physical meaning. This term signifies the rate of change of the total amount of fluid property $ in the control volume, The product n.péu expresses the flux component of property $ due to fluid flow along the outward normal vector n, so the second term on the left hand side of above equation, the convective term, therefore is the net rate of decrease of fluid property ¢ of the fluid element due to convection. A diffusive flux is positive in the direction of a negative gradient of the fluid property &, Le. along direction ~ grad . For instance, heat is conducted in the direction of negative temperature gradients, Thus, the product n . (-f grad ) is the component of diffusion flux along the outward normal vector, so out of the fluid element.Similarly, the product n .(F grad ¢), which is also equal to f(-n . (-grad )), can be interpreted as a positive diffusion flux in the direction of the inward normal vector -n, Le. into the fluid element. The first term on the right hand side of equation 1 of previous slide, the diffusive term, is thus associated with a flux into the element and represents the net rate of increase of fluid property of the fluid element due to diffusion. The final term on the right hand side of this equation gives the rate of increase of property as a result of sources inside the fluid element. In words, this relationship can be expressed as follows: Net rate of decrease Net rate of . Rate of increase of @ duc to increase of @ ant of inside the + convection across = duc to diffusion ae. This clarifies that integration of the PDE generates a statement of the conservation of a fluid property for a finite size (macroscopic) control volume. In steady state problems the rate of change term in equation 1 of slide 42 is ‘equal to zero.to the integrated form of the steady transport equation: {n.coounst= [ace gad ges [sar ‘ 1 ew * Intime-dependent problems it is also necessary to integrate with respect to tit tover a small interval At from, say, t until t + At. * This yields the most general integrated form of the transport equation: [elf fons fare [se a aa wa we Classification of physical behaviours of problems After deriving the conservation equations of fluid flows, we focus on the issue of the initial and boundary conditions that are needed in conjunction with the equations to construct a well-posed mathematical model of a fluid flow. First we distinguish two principal categories of physical behavior: > Equilibrium problems > Marching problems “This problem is one-dimensional and governed by the equation kd°T/dx? = 0. Under the given boundary conditions the temperature distribution in the x- direction will, of course, be a straight line. A unique solution to this and all elliptic problems can be obtained by specifying Conditions on the dependent variable (here the temperature or its normal derivative the heat flux) on all the boundaries of the solution domain. Problems requiring data over the entire boundary are called boundary-value problems. & An important feature of elliptic problems ts that a disturbance in the interior of the solution, e.g. a change in temperature due to the sudden appearance of a ‘small local heat source, changes the solution everywhere else. Disturbance signals travel in all directions through the interior solution.Marching problems * Transient heat transfer, all unsteady flows and wave phenomena are examples of problems in the second category, the marching or propagation problems. * These problems are governed by parabolic or hyperbolic equations. * However, not all marching problems are unsteady. * Parabolic equations describe time-dependent problems, which involve ‘Significant amounts of diffusion. * Examples are unsteady viscous flows and unsteady heat conduction. * The prototype parabolic equation is the diffusion equation which is written as: Fe oe + The transient distribution of temperature (again > = T) in an insulated rod of metal whose ends at x = 0 and x = L are kept at constant and equal temperature Ty is governed by the diffusion equation. * This problem arises when the rod cools down after an initially uniform source is switched off at time t = 0. . a ee es a ARMS Ps teins wksSeen Zoo The steady state consists of a uniform distribution of temperature T = T, throughout the rod. The solution of the diffusion equation (on previous slide) yields the exponential decay of the initial quadratic temperature distribution. Initial conditions are needed in the entire rod and conditions on all its boundaries are required for all times t > 0. This type of problem is termed an initial-boundary-value problem. A disturbance at a point in the interior of the solution region (ie. 0 0) can only influence events at later times t >t). The solutions move forward in time and diffuse in space.The occurrence of diffusive effects ensures that the solutions are always smooth In the interior at times t » 0 even if the initial conditions contain discontinuities. The steady state is reached as time t > a and is elliptic. This change of character can be easily seen by setting d¢/dt = 0 in the parabolic equation, The governing equation is now equal to the one governing the steady temperature distribution in the rod. Hyperbolic equations dominate the analysis of vibration problems. in general they appear in time-dependent processes with negligible amounts of energy dissipation, The prototype hyperbolic equation Is the et ? be * 4 ow Ff oe « The above form of the equation governs the transverse displacement (o=viota ‘string under tension during small-amplitude vibrations and also acoustic oscillations, The constant c is the wave speed.* Itis relatively straightforward to compute the fundamental mode of vibration of a string of length L using hyperbolic equation. * Solutions to wave equation and other hyperbolic equations can be obtained by specifying two initial conditions on the displacement y of the string and one condition on all boundaries for times t > 0. * Thus hyperbolic problems are also initial-haundary-value problems. Sobsvon prouder, soect*cabor He tah Reg ord y/ate, £020 for first cycle 0 Quis + wis ceatral — Me i vi soe (3), ~ (ax?¢ & ey wepyaicee SS | YIN OP rLel) Geek a*fegey £ £ What heppens at a boundary? . Coisilr te gure eins wikca nines a portion of a boundary to a flow field, with the y axis perpendicular to the boundary. * Letgrid paint 1be on the boundary, with points 2 and 3 a distance Ay and 24y above the boundary, respectively. * We wish to construct a finite-dilference approximation for bu/6y at the boundary.Assume at the boundary shown in the figure that ucan be expressed by the polynomial u=atby+or Applied successively to the grid points in this figure, the above equation yields at grid point 1 where y = 0: uj=a and at grid point 2 where y = Ay, uw = a+b Ay+c(Ay)? and at grid point 3 where y = 2Ay, uy = a + b(2Ay) +.¢(2Ay)? Solving the above 3 equations for b, we obtain: _ ~3uy + 4u2 = uy 2hyReturning to the polynomial equation and differentiating with respect to y, Ou ==b a +2cy Above equation, evaluated at the boundary where y=0, yields (5) Combining the above 2 equations, we obtain (F) _ ~3uy + 4uy = uy dy, 1 ey This equation is a one-sided finite-difference expression for the derivative at the boundary also called as one-sided because it uses information only on one side of the grid point at the boundary, namely, information only above grid point 1. Also, this equation was derived using a polynomial expression rather than a Taylor series representation.In order to check the order of accuracy of this equation wrt Taylors series, consider a Taylor series expansion about the point 1. worom (SoCs) a Go)et™ On comparing the above equation with the polynomial equation (eq 1 slide 72), we find that the polynomial expression is the same as using the first three terms in the Taylor series. Hence, the polynomial expression is of O(Ay)*. If we compare the numerator of last equation of previous slide, we find that u,, u,, and u; can all be expressed in terms of the polynomial expression and since that equation is of O(Ay)’, then the numerator of last equation of previous slide is also of O(Ay)?. However, in forming the derivative in last equation of previous slide, we divided by Ay, which then makes that equation of O(Ay)’. Thus, we can write from last equation of previous slide ‘Ow 1 + 4ua — uy 2 =) = — > OH) G)-—" adThis is our desired second-order-accurate difference quotient at the boundary. This equation along with equation 1 of slide 71 are called one-sided differences, because they express a derivative at a point in terms of dependent variables on ‘only one side of that point. Moreover, these equations are general; i.e., they are not in any way limited to application just at a boundary; they can be applied at internal grid points as well. DIFFERENCE EQUATIONS Till now, we discussed the representation of a partial derivative by means of an algebraic finite-difference quotient. Most partial differential equations involve a number of Partial derivative terms. When all the partial derivatives in a given partial differential equation are replaced by finite-difference quotients, the resulting algebraic equation is called a difference equation, which is an algebraic representation of the partial differential equation. The essence of finite-difference solutions in CFD is to use the difference quotients to replace the partial derivatives in the resulting in a system of algebraic difference e variables at each grid point. Governing flow equations, quations for the dependentIn the present section, we examine some of the basic aspects of a difference equation. For simplicity, we choose to examine a partial differential equation which is less elaborate than the governing flow equations. For example, let us consider unsteady, one-dimensional heat conduction equation with constant thermal diffusivity, given below. oT a PT oH Ox? We know that this is a parabolic partial differential equation and lends itself to a marching solution with respect to time t. Let us replace the partial derivatives in first term on LHS with finite-difference quotients. The above equation has two independent variables, x and t. If we consider the grid where iis the running index in the x direction and nis the running index in the t direction.When one of the independent variables in a partial differential equation is a marching variable, it is conventional to denote the running index for this marching variable by n and to display this index as a superscript in the finite- difference quotient. Hence, we can write, arty _T*'-T (3) ola A -(s3) s+ & ) ar}, 2 Where, the truncation error is the same as that displayed in previous equations. Also, let us replace the x derivative with a central difference differences, we get (2y- May Th - (aa Where, the truncation error is the same as that displayed in previous slides. Let us write the base equation as: ar PT Be coar* Inserting values of equation 1 and 2 into 3 we have: ecaomena. | Gree | Difference equation “Truncation error * Inthe above equation, the left-hand side is the original partial differential equation, in which the first two terms on the right-hand side are the finite- difference representation of this equation, and the terms in the vertical brackets give the truncation error for the difference equation. * Hence, just the difference equation from above equation can be written as Te! = Ty _ a(t, - 2+ T,) At (ax)?ERRORS AND AN ANALYSIS OF STABILITY ay Clearly, the new difference equation involves products of the depe nden variables, such as [a(T/°*)}.7,""2, [a(Tj"2)).Ti"? and (a(T").7,.°"2 In other words, the resulting difference equation is a nonlinear algebraic equation. An implicit solution would therefore demand a simultaneous solution of a larg system of nonlinear equations which is an exceptionally difficult task. This is a tremendous disadvantage of an implicit approach.The numerical solution of model equation is influenced by two sources of er! Discretization error: The difference between the exact analytical solution of thi partial differential equation and the exact (round-off-free) solution of the corresponding difference equation. The discretization error is simply the truncation error for the difference equation plus any errors introduced by the numerical treatment of the boundary conditions. Round-off error, the numerical error: It is introduced after a repetitive number of calculations in which the computer is constantly rounding the numbers to some significant figure. A = analytical solution of partial differential equation D = exact solution of difference equation N = numerical solution from a real computer with finite accuracy Discretization error = A - D Round-off error =« = ~D From Eq. (4.54), we can write N=D+ex For the remainder of our discussion in this section we will simply call rou error (€) as ‘error’. The amount of error added in each step can be written e prt+ey! ree) Ce tye Gp: — 27 G+) © ’ oe fi ). - a Ar (Ax)? — If Dis thee i n, Le,, it exactly satisfies difference equation, we can write ad Y Dit! Di _ Di. -2D1+ Dis ® Qt et ty adt (Ax) Subtracting the second equation from first, we have © - ®@ we [eee eS e A aAt (Ar) Hence, we see that the error € also satisfies the difference equation. If errors €, are already present at some stage of the solution of this equation, then the solution will be stable if the €'s shrink, or at best stay the same, as the solution progresses from step n to n+1. On the other hand, if the ¢;'s grow larger during the progression of the solution from steps n to n+1, then the solution is unstable or |* Thatis, for a solution to be stable o = For the unsteady, one-dimensional problem of model equation, the round-ot error can be plotted versus x at any given time step (given in figure below). z 4 « Here, we assume that the length of t solved is denoted by L. — S * For convenience we place the origin at the m left boundary is located at -L/2 point of the domain; hence the yis at +1/2. aaThe distribution of € along the x axis is represented by the rather random variation as shown in previous figure. Note that «= =Qatx= ~UZINAT/Z, because there are specifi ed boundary values at both ends of the domain, and hence no errors introduced, Le., the jounda values are always fed in as exact, known numbers. At any given time, the random variation of € with x can be expressed analytically by a Fourier series as follows: dx) =o m ; jk x ini Above equation represents both a sine and a cosine series, since e,," = cos k,,X +isin ie called the wave number. The real part of this equation represents the error. To understand wavenumber, consid st the sine function plotted as a function of x as shown in next figure. ot 2”2 ‘ By definition, the wavelength Alsthe is the = ‘over x encompassing one complete wavelength. Therefore, a form of this sine function is yosi ry kn or y=sink,x Here, the subscript m represents the number of waves fitted into the given 7 interval L (in this case m = 2). Hence we can say that the summation over m denotes the sequential addition of sine and cosine functions with sequentially increasing wave numbers. That is, the original Fourier transform equation is a sum of terms, each representing a higher harmonic. ane(di By definition, the wavelength A is the interval over x encompassing one complete wavelength. Therefore, a form of this sine function is nx . y= sin or ys sink,x Here, the subscript m represents the number of waves fitted into the given interval L (in this case m = 2). Hence we can say that the summation over m denotes the sequential addition of sine and cosine functions with sequentially increasing wave numbers.When taken out for an infinite number of terms, Fourier transform equation car represent a continuous variation of € as a function of x. However, in regard to a practical numerical solution which involves only a finite number of grid points, there is a constraint imposed on the number of terms in Fourier transform equation. To see this more clearly, consider the figure given below, which shows the interval L over which the numerical calculations are being made. The largest allowable wavelength is A,,,,, = L; this is the wavelength for the first term in Fourier transform equation and corresponds to m = 1. In turn, the smallest possible wavelength is that having all three zeros of the sine (or Aeg = 2h cosine) function going through < duel three adjacent grid points. , Hence, the smallest allowable wavelength is Aji, = 24x. 9sIf there are N + 1 grid points distributed over the interval L, tl 5 intervals between these grid points, and hence Ox = L/N. 7 Thus, Amin= 2L/N and hence, we can write, x Raat a2 sv SQN Comparing this equation with previous one, we see that min equatio Ty to N/2. for N+ 1 grid points we can write, N?2, v (@} a A er ae This equation gives the spatial variation at a given time level n. 7 However, we are interested in the variation of e with time. Therefore, we extend the above by assuming the amplitude A, isa function of time. 3 . N/z (x,t) = D> An(tde* - 4 = 42 mal ;* Moreover, if we assume an exponential variation of error with time, tend to grow or diminish exponentially with time. we write [2 Ni a r—_ ex.) = Dots ge * Where, ais constant (which may take on different values for different m's). + Hence, the previous equation represents a final, reasonable form for the variation of round-offerrarin both ‘Spaceandtime. * Since the original difference equation (equation 1 of slide 80) is linear and since the round-off error satisfies the same difference equation as written in last quation of slide 91, then we can substitute the values of its value in previous equation. This can be done as the behavior of each term of the series is the same as the series itself. lence, let us deal with just one term of the series and write €m(% 0) = nlf) = ee * This general equation can be used to understand the stability Sharacteristics. + Let us now proceed to find out how € varies in steps of time. Alt ON thax — Picker Bighalxt 4x) _ 2etettar 4. Alethnlx Ax) a at (Ax)?Divide the previous equation by erent ona un vey Ar ik, Ox te eat ett 2b entlnte >») aA (Ax)* or eo 8 Ee ern 2) —® Recalling the identity in Ox thm Ax cos (EAs) ee + S — AW Equation 2 can be written as 2a At Co aot aye te) - Recalling another trigonometric identity sin? #2 BE _ 1 = 608 (ln Ax) 2 2 ieHence, equation 4 of previous slide becomes finally eu _ 4a ar SOA sin? Kon Ax ~ (xy 2 From the above equation a AM AS a | tet Combining equation 1 of slide 92 with the above both equations, we have \e“|= 1A ge <1 (ay sat Above equation must be satisfied to have a stable solution, wherein the factor nee is called the amplification factor and is denoted by G. 9= Oat the steady state and t~ is a fictitious time. n With these definitions and combining the equation of continuity and artificial density, we may write the incompressible Navier-Stokes system of equations in the form aw aw la aw art May” Tei ™aR,)Where, ~ of ee av” By 0 Sa Li A= ow" ns we rel 6, =[ 8] 0p 00 00 BO : adi _|t 2u 0 0 a _|o vu 0 M=GWwolo vu o| 2taw=l1 o 2 of 43-5 Ow w 00 w vy. Let us now investigate the eigenvalues of Ai, IA; - AT =O where the eigenvalues of Ai (i = 1, 2, 3) are, respectively, (au, we Vie +B), (wy, vi V7 +8), (w, w, wt Jw? +B) ; + Inwhich VB is sthe artificial speed of ean (often called the artificial compressibility factor) with B being chosen adequately (between 0.1 and 10) i ta* The value of B is maintained low enough (close to the convective velocity] ‘overcome stiffness associated with a disparity in the magnitudes of the eigenvalues, but high enough such that pressure waves (moving with infinite speed at incompressible limit) be allowed to travel far enough to balanc viscous effects. 7 * Asa result, the conservation of mass or incompressibility condition is assure means of an artificial compressibility. + In this process, it is possible to obtain the correct pressure distributions. SEMI-IMPLICIT METHOD FOR PRESSURE-LINKED EQUATIONS (SIMPLE) * [tis well known that, if the finite difference equation is written in control volume grids for continuity v;; = 0, this will lead to nonphysical, checkerboard-type oscillations of velocity in each one-dimensional direction (same values repeated at every other node, assuming that the velocity distribution between the adjacent nodes is linear). 107Hence, equation 4 of previous slide becomes finally 4a 2 f*al- At be be ae _ From the above equation tan All pas ‘ “q” eda =¢ ‘Combining equation 1 of slide 92 with the above both equations, we have 6 |. ih 40 At) by Ax [Se l-tenie ae So Abave equation must be satisfied to have a stable solution, wherein the factor e ait “a 2G Is called the amplification factor and Is denoted by G,* Evaluating the inequality in previous equation namely, G < 1, we have two possible situations which mus} wld simultanedualy: 4a At sin 2 km Ar 2 21 I-55 — OH) G)-—" adThis is our desired second-order-accurate difference quotient at the boundary. This equation along with equation 1 of slide 71 are called one-sided differences, because they express a derivative at a point in terms of dependent variables on ‘only one side of that point. Moreover, these equations are general; i.e., they are not in any way limited to application just at a boundary; they can be applied at internal grid points as well. DIFFERENCE EQUATIONS Till now, we discussed the representation of a partial derivative by means of an algebraic finite-difference quotient. Most partial differential equations involve a number of Partial derivative terms. When all the partial derivatives in a given partial differential equation are replaced by finite-difference quotients, the resulting algebraic equation is called a difference equation, which is an algebraic representation of the partial differential equation. The essence of finite-difference solutions in CFD is to use the difference quotients to replace the partial derivatives in the resulting in a system of algebraic difference e variables at each grid point. Governing flow equations, quations for the dependentIn the present section, we examine some of the basic aspects of a difference equation. For simplicity, we choose to examine a partial differential equation which is less elaborate than the governing flow equations. For example, let us consider unsteady, one-dimensional heat conduction equation with constant thermal diffusivity, given below. oT a PT oH Ox? We know that this is a parabolic partial differential equation and lends itself to a marching solution with respect to time t. Let us replace the partial derivatives in first term on LHS with finite-difference quotients. The above equation has two independent variables, x and t. If we consider the grid where iis the running index in the x direction and nis the running index in the t direction.When one of the independent variables in a partial differential equation is a marching variable, it is conventional to denote the running index for this marching variable by n and to display this index as a superscript in the finite- difference quotient. Hence, we can write, arty _T*'-T (3) ola A -(s3) s+ & ) ar}, 2 Where, the truncation error is the same as that displayed in previous equations. Also, let us replace the x derivative with a central difference differences, we get (2y- May Th - (aa Where, the truncation error is the same as that displayed in previous slides. Let us write the base equation as: ar PT Be coar* Inserting values of equation 1 and 2 into 3 we have: ecaomena. | Gree | Difference equation “Truncation error * Inthe above equation, the left-hand side is the original partial differential equation, in which the first two terms on the right-hand side are the finite- difference representation of this equation, and the terms in the vertical brackets give the truncation error for the difference equation. * Hence, just the difference equation from above equation can be written as Te! = Ty _ a(t, - 2+ T,) At (ax)?