The solution of differential equations of fluid flow by numerical program Levent Yilmaz1 1Civil Engineering Department, Technical University oflstanbul, Maslak, 80626 Istanbul, Turkey; e-mail: lyilmaz@itu.edu.tr Abstract: The finite difference methods are used in the solution of engineering problems and in fundamental studies of the behavior of fluids. Considering the basic equations of energy, momentum, and mass conservation in a form in which the fluid properties are constant. In the differential equations the limiting assumptions of linear behavior, incompressibility, continuum properties, and simplifications of boundary properties are given for the solution. Keywords: fluid flow, finite difference method, incompressibility Introduction The basic equations of energy, momentum, and mass conservation are given in the form in which the fluid properties are constant. These may be written as follows: Received: June 21, 2007 Accepted: July 10, 2007 Energy: f ? \ p — + uN I = kV2T- PV.u + X$.u)2 + [M> dt (1) V Momentum: ft) s P — + uN M = pg~VP + (Â + jU)V(V.w)+jUV2; ldt J Mass: id ') — + u.V p = -p(V.u) (!) Preliminary note where I is the internal energy per unit mass, T is the temperature, P is the pressure, u is the velocity,p the density, andcP the fTictional heating. k,fl, and/l are the thermal conductivity and the coefficients of dynamic and kinematic viscosity and in the Equation of (2) g is given as the body force. The equations considered here are of flows in which the density variations are small and can be neglected in the sense of the Boussinesq approximation. The basis of this approximation is that there are flows in which the temperature varies little, and therefore the density varies little, yet in which the buoyancy drives the motion. Thus the variation in density is neglected everywhere except in the buoyancy term. Let />/, denotes the density at the bottom where temperature is Tb. For small temperature difference between the top and bottom layer we can write p^pb [l-a(T-T)] where a is a coefficient of volume expansion. For liquid a is in the range lO"3 to lO"4. For a temperature variation of moderate amount we have M=kz^=„T-Tb i« i Pb Pb but the buoyancy term g (p — pb) is the same order of magnitude as the inertia, acceleration or the viscous stress so is not negligible. The differential in density in the continuity equation are of the order a and hence neglected to give = o (4) dx,. as for an incompressible fluid. Then the stress tensor becomes / r=-PS,+M ^ du. du,^ —L +—. dx. dx, \ J 'J Again the momentum equation becomes -gPb [l-«(T-Tb )]5,3- ax. du . du — + u. —1 dt ' dx,. V 1 p V which after using the Boussinesq approximation becomes (after a little manipulation) whereA = d2 /9*2 is the Laplacian operator. 1/2 a and thus the ratio of(P term to the term due to conduction ofheat is of the order of jiagd k and this is in the range from 10 "7 to 10"8 for gases and liquids and the viscous dissipation is neglected. The term Thus though8wk ldxk = 0 holds in the equation of continuity, we should not use this relation in the heating due to compression term for gases. For liquids however, this heating is negligible. Thus the final form of the energy equation is whereK = k lpbcp for a perfect gas and k lpbc for liquids. Equations (4), (5) and (6) are called Boussinesq equations and describe the motion of a Boussinesq fluid. The density variation is neglected in all terms except that involving the body force g (3). In this case we take Using perfect gas relations p = pRT, R = (cp — p ) anda = 1l T weget (6) p = po+8p = po-poa{T-T0) where a is the volume expansion coefficient of the fluid. Neglecting viscous heating 0 we shall take I = cvT, where cv is the heat capacity at constant volume. Grouping the constant fluid properties in the form X = and v = — (8) Pcv P we may rewrite (1) through (3) as — + u.V [dt , T = -x (V2 T (9) 9 „ — + uN dt i ^ ^ VP g-a(T-T0)g- —+vV2u (10) p Vm = 0 (11) The vorticity is defined by (0 = Vxu (12) Equation (10) may be expressed in terms of vorticity in the form jjUM.V^ = -aV(T)xg+vV2ffl (13) Equations (9), (11), (12) and (13) can be used to establish a set of two-dimensional equations that then can be expressed in finite difference form. In two dimensions the zero divergence of velocity (11) permits us to make use of a stream function defined by the equations m = and v = (1#) dy dx in a rectangular coordinate system. Through (12) the mass conservation equation leads to the relation ^V+^V^ (15) a*2 ¿y2 The vorticity has only one component in the two-dimensional case and is hence treated as a scalar. If g is taken as a constant body force (gravity) directed in the negative y direction, the momentum equation (13) becomes d IT = u V= v Figure 1. A representative section of the lattice of points used in the finite calculation (Abrahamand Tiller, 1972) We consider the cell as a volume element of fluid and prescribe mass conservation by Mi+1J+1/2 Ui,j+1/2 Vi+1/2,j+1 Vi+1/2,j _ q Ax Ay When the velocities are given by (18) u i, j+1/2 y^i, i+1 y^i, i i ■ —--- and v Ay i 1 2, i . Vi, Ax (19) the continuity equation is automatically satisfied in difference form. Further, the continuity condition is satisfied for groups of four cells in the form ui+1,j~ui-hJ | vi,j+i~ vi,-M _Q /20) 2Ax 2Ay if we prescribe the last equation u — Wi+i, j+i~ Wi+i, y-1 ui+1, /+1/2 ui+1, y—1/2 Mi+1 i -_- (21) w 2Ay 2 Also, if the convection of energy is given in the form (»tLj-(»tLJ | (vT)i,;+x~ ('vT\(22) 2&x 2Ay no extraneous source of sinks of heat will result (ABRAHAM AND TILLER, 1972). This is true both because of the manner of centering the velocities and because of the use of differences of products of the variables. An analogous form to Equation (22) above also may be used for the transport of vorticity. The vorticity at the point i,j may be defined in terms of nearest neighbor velocities as Ui, j-1/2 _ Ui, j+1 /2 Vi'+1/2, j+1~ Vi-1 /2, j <°U j=—-A---+-—K---(2") Ay Ax Using Equation (19) above, a = (Vi,w -Vi,w-^i,1 - Vi',w) , ^i,w -WHUw)- (V-U,w -Vi,w (24) iJ Ay2 Ax2 An explicit differencing scheme with centered time differences is used in this formulation (ABRAHAM and Tiller, 1972). A variable value at the current time is designated by a superscript n. Backward and forward time values are labeled with superscripts n-1 and n+1, respectively. Tn+1 _1 tJ ~ ^ J = f(un,vn,Tn) (25) 2At J v ) For numerical stability it is necessary to use a special form for the conduction terms , also involving backward and forward time values (RICHTMYER, 1957). In these terms the Dufort-Frankel method is used as , (DUFORT AND FRANKEL , 1953) , rpn+1 ÎTTH-1 pn _ i,j_i ,j ''J ~ o (26) The two-dimensional conduction terms are fc+i- u -'-T: .v J, fc, ^-^fe-f:, J at 2 A^2 (27) It is convenient for variable mesh specification to introduce a factor f = Ax /Ay . Then if Ay = a, Ax = yh, the complete set of difference equations becomes (Abraham and Tiller, 1972) (a) Energy conservation (Abraham and Tiller, 1972) T i_ 1/ 1 , 2Xà f 2+1 a2 f2 ^ 8t iuTYi+1, w - (mU ) i+1, j f (vr)n^1 2* & T +Tn ^ f2 un--^ ^ ur)} (28) (b) Vorticity (Abraham and Tiller, 1972) )i i f +(v®L+1- (v©L- 1 (Tn+1 _ Tn+1 fa 1 j ^ j) 2v St nn ®n+ 1,w + ft>i 1, , f2 nn ■+ dite) i " FJUHTMAP (Mç pmi) PRIHTPAÄA (Turunda: pint) PRINTDATA. (Diiipimi) -4- CHECKSTA (OietJ: StsMHy) —1 VtONVECCAL ifDr-* (awHtion) HCONVECCAL DUFUCAL ^fflfusillj Figure 2. Block Diagram of the numerical program For understanding of the numerical methods, it is included a skeleton program for computing incompressible flows. Only the main equation of the incompressible flows are used, giving just the basic flow dynamics in the absence of all but the convective or inertial and viscous forces. There is a certain completeness in the given program in that additions may be made as insertions without involving major reprogramming. We begin the block program by entering INPUT DATA which reads input data. The input data include a minimum of input parameters such as grid size, grid aspect ratio, print interval, and Reynolds number. The INPUT DATA program also is used to consolidate parameters to provide for simplified constants that will later speed the computation. It also determines the size of the initial time step. The program INPUT GEOM, like INPUT DATA, reads input data but is specialized to give a flexibility of geometry with a minimum number of data words. The boundary must consist of straight line segments passing through mesh points, but this is the only restriction on the geometry. INPUT GEOM also establishes the initial values for boundary points. The program PRINT GEOM gives a symbolic printout of input geometry for checkout purposes. The symbols used characterize the type of bundary condition that is to be imposed. Before starting the computation it is convenient to introduce an additional program GUESS SOL to establish initial guess values based on already prescribed boundary values. An intelligent guess enables us to speed up the convergence of our initial yf solution. This initial guess varies from a crude layout of values given by a linear interpolation between boundary values all the way to known exact solutions. If some initially specified distribution is to be perturbed, this is also programmed in the subroutine (GUESS SOL). The program enters the subroutine CALCU PSI for actual solution of Poisson's equation. Here, through an iteration averaging process, it is obtained the stream function distribution. Most often the stream function is the potential flow solution for the given geometry. After obtaining convergence to some specified accuracy it is established boundary values of vorticity at no-slip surfaces. This is done in the subroutine WATT. At this point in the calculation we have a complete solution with a consistent set of numbers for yf andft). It also may be desirable to include input parameters or functions evaluated from the data for identification relative to the printouts. The program PRINTMAP is used to print a sheet of normalized and highly rounded values to give a picture that bears some resemblance to a contour map. A good printer map in some cases may provide a substitute for a contour map (ABRAHAM AND Tiller, 1972). The PRINTDATA program simply prints out the data values at all grid points in a layout that is geometrically similar to the given maps but usually requires several sheets to contain the information. Finally PRINTPARA serves as a labelling program, providing orientation information for the given printouts. This is a logical point to terminate the program. However, we have not yet advanced the solutions in time, but we have only obtained a consistent initial solution. From this initial solution we now may reevaluate the size of time step required to maintain stability. This is done in the program CHECKSTA which tells us to reduce the time step that was initially furnished, we can use CHECKSTA to reevaluate any constants that are affected by this change. The arrow in the block diagram symbolizes a repeated reentry into the program until stability conditions are appropriate for proceeding with the computation. The next set of subroutines allows for the advancement of the vorticity distribution to correspond to a new time, one time increment as At established by CHECKSTA beyond the previous solution. VCONVECCAL computes the convection of the vorticity distribution in the coordinate direction of the u velocity in the x-coordinate. Similarly, HCONVECCAL gives the convection associated with the v velocity in the y-direction. Finally, DIFFUCAL incorporates the diffusion of vorticity that will occur over the time interval. The time-advanced distribution of vorticity will induce a change in the flow so we now must evaluate a new stream function distribution by reentering CALCU PSI. Each passage through this loop of subroutines amounts to the advancement in time of the solution by a small discrete time interval. We may or may not record data at each time step, but the stability is tested every time step. One important point about the program listing is that there are no details missing so far as the physics of the calculation is concerned. Because the program has been References Abraham, F. F., and Tiller, W. A. (1972): An introduction to computer simulation in applied science. Plenum Press, New York - London. Dufort, E. and Frankel, S. (1953): Stability conditions in the numerical treatment of parabolic differential equations. Math. Tables and Other Aids to Computation:, Vol. 7, pp. 135-152. Fromm, J., (1963): A method for computing nonsteady, incompressible, viscous fluid flows. Los Alamos Scientific extracted from a more extensive set of programs, there are more parameters defined in the data lists than are here applicable. The subroutine calls sometimes contain redundant arguments. Generally, the best procedure for correction a program from listings is to draw flow charts, including all detail, for each subroutine. INPUT data involves array size information and data as read by INPUT DATA and INPUT GEOM. The array size data is given in the first six lines of the mainline program, while the remaining data is included at the end of the program deck. Logical arrays are used to establish boundary condition for both yr andft) fields. Laboratory Report:, LA 2910, Los Alamos. RiCHTMYER, R., (1957): Difference Methods for Initial Value Problems. Interscience Publishers, Inc., New York. Rushton, K. R., and. Redshaw, S. C. (1979): Seepage and Groundwater Flow: Numerical Analysis Description. Chichester John Wiley & Sons, pp. 339.