Computational s for Fluid Dynamics-Fortran Source CodeI-Joel HFerziger and Milovan Peric.doc

Joel H. Ferziger / Milovan PericComputational Methods for Fluid Dynamicsthird, rev. editionFORTRAN Source Code-IDepartment of Disaster Prevention and Mitigation,Beijing Jiaotong UniversityMay, 2010目录一、源代码说明2二、勘误4三、源代码及相关说明63.1 一维对流扩散方程有限差分Fortran源代码63.1.1 P1DS.F程序说明63.1.2 P1DUS.F程序说明73.1.3 P1DS.F源程序83.1.4 P1DUS.F源程序143.2 笛卡尔坐标系下二维热及流体流动Fortran源代码193.2.1 GRID.F程序说明193.2.2 PSC.F程序说明203.2.3 PSCUS.F程序说明213.2.4 PCOL.F程序说明213.2.5 PLOT.F程序说明223.2.6 PCOL1.F程序说明223.2.7 GRID.F源程序233.2.8 PSC.F源程序343.2.9 PSCUS.F源程序443.2.10 PCOL.F源程序543.2.11 PLOT.F源程序883.2.12 PCOL1.F源程序122一、源代码说明Latest update: 1999-04-08-Dear reader,in this directory you will find the following files and sub-directories(directory and file names are written using capital letters to emphasize them; on the computer, all directory and file names are in lower case letters).README - This file.ERRATA - File containing information about typing errors found so far in the book; the latest update is indicated at the beginning of the file.NEWPRINT - Sub-directory containing postscript files in which the most important changes made in the 2nd corrected printing of the first edition and in the second edition are documented. These files have been compressed using gzip. To recover the postscript file after you have transfered the files name.ps.gz, do the following: gunzip name.ps.gzname is the file name before extension ps.gz as you see it in the sub-directory; substitute it accordingly for the various files.The following directories contain computer codes that have been updated to comply with the notation used in the book. They also include many lines of comments to make them easier to understand.1D - Sub-directory, containing codes for solving the one-dimensional steady and unsteady convection/diffusion equation, using Finite Difference method and different schemes for convective terms and time integration. Change to this directory and read the local readme file for further information on files in that directory.2DC - Sub-directory, containing codes for solving the two-dimensional convection/diffusion equation (scalar transport in a given velocity field, steady and unsteady) using Finite-Volume method and Cartesian grid, and for solving the steady and unsteady Navier-Stokes equations using Cartesian grids and a collocated arrangement of variables. It also contains a code for generating Cartesian grid (single or multigrid), and for plotting the results of computations. The latter code produces postscript files of each plot page: grid, velocity vectors, profiles of velocity or temperature, isobars, isotherms or streamlines, and color-fill plots of pressure, temperature or streamlines. Change to this directory and read the local readme file for further information on files in that directory. The subdirectory PIPE contains a version of the code that is set for pipe flow; it includes inlet and outlet boundary conditions, and the global mass flux correction in the SIMPLE-algorithm.SOLVERS - Sub-directory, containing codes for solving the two-dimensional and three-dimensional Laplace equation using a variety of linear equation solvers. The 2D code contains seven solvers to choose from, including ILU (SIP), ICCG, and two multi-grid versions. The 3D code contains ILU (SIP), ICCG and CGSTAB. These solvers can easily be adapted and implemented in other codes solving different equations. Change to this directory and read the local readme file for further information on files in that directory.STAG - Sub-directory, containing a code for solving the two-dimensional, steady Navier-Stokes equations using Finite-Volume method and Cartesian structured grids with a staggered arrangement of variables. Change to this directory and read the local readme file for further information on files in that directory.2DGL - This subdirectory contains the new codes for solving the laminar two-dimensional Navier-Stokes equations (steady or unsteady) using Finite-Volume method and non-orthogonal structured grids with a co-located arrangement of variables. H-, O- and C-grids can be used, with multiple grid levels. Two subdirectories exist: SG - This sub-directory contains two single-grid versions of the code: in one, the grid-non-orthogonality and non-smoothness are specially accounted for when solving the pressure-correction equation and when computing the gradients, and in the other version not. The solution from one grid level serves to initialize the iterations on the finer grid.MG - This sub-directory contains the version of the above code with multi-grid acceleration of outer iterations. More information is available in the README file in this directory. 2DGT - This subdirectory contains the new codes for solving the RANS-equations using Finite-Volume method and non-orthogonal structured grids with a co-located arrangement of variables. H-, O- and C-grids can be used, with multiple grid levels. Two subdirectories exist:KEPS - This sub-directory contains the version of the above mentioned code with the k-epsilon model of turbulence(for high Reynolds numbers) with wall functions.KOMEGA - This sub-directory contains the version of the above mentioned code with the k-omega model of turbulence(for high Reynolds numbers) with wall functions.More information is available in the README file in this directory. The files in the following directories contain codes that have been written and used in the past; they do not comply with the notation and algorithms described in the book completely and will be edited in the near future. See README files in these sub-directories for information on latest updates.3DC - Sub-directory, containing codes for solving the three-dimensional Navier-Stokes equations (steady) using Finite-Volume method and Cartesian grids with a co-located arrangement of variables. The codes currently available will be re-written to comply with the notation used in the book; check this file for information about the update. Change to this directory and read the local readme file for further information on files in that directory.PARALLEL - Sub-directory, containing a parallelized version of the code for solving 2D steady Navier-Stokes equations on Cartesian grids with a co-located arrangement of variables using PVM message-passing software. Change to this directory and read the local readme file for further information on files in that directory.As already indicated above, each sub-directory contains a README file with detailed information about files in that directory. Where indicated, new versions of files will be included in the near future.The codes provided here are intended to demonstrate the implementation of the algorithms described in the book. They can be useful to learn and demonstrate the use of numerical solution techniques; students may be given tasks to modify these codes by implementing different discretization or solution methods, or different boundary conditions. The updated codes are written so that the algorithm is implemented in a most straight-forward way, and that it is easy to recognize the equations from the book; they can be re-arranged to be more efficient or shorter. Also, these codes might serve as starting point for developers of new solution methods, or of solution methods for different problems. For example, although written for a single-block structured grid, the codes in the directories 2DGL and 2DGT can easily be extended to block-structured or unstructured grids, both two- and three-dimensional. This is not attempted here since the codes needed to generate such grids are far more difficult to write than the extension of the flow solver; if an appropriate grid generator is available, the extension of the flow solver is relatively easy (some instructions are included as comments in one of the codes in directory 2DGL).For more information on these and other codes, interested readers maycontact: pericschiffbau.uni-hamburg.de二、勘误ERRATA FILE Last updated on March 21, 2001-First Edition-All of the typing errors mentioned below have been corrected in the 2nd printing of the book; this does not guarantee that there are no typing errors that remained undetected so far.The following major typing errors have been found so far in the first printing of the book:1. On page 105, in the CGSTAB algorithm, last equation: the parameter gamma should be replaced by alpha, i.e. the expression should read: rhok = w - alphak v The solver included in the file LAPL3D.F in the directory SOLVERS is correct.2. On page 170, Eq. (7.51), a minus sign is missing in front of u_y.3. On page 183, Eq. (7.101), a minus sign is missing after equal sign in expressions for the coefficients A_E, A_W, A_N and A_S. The expressions used in the code PSTAG.F (sub-directory STAG) are correct.4. On page 328, after end of paragraph after Eq. (12.24), there is a reference to Eq. (12.24); the equation meant is Eq. (12.23). The same correction applies to the reference to Eq. (12.24) on page 330, at the end of paragraph after Eq. (12.28).5. On page 231, Eq. (8.62), the derivative of p should be with respect to xi, not n, as in Eq. (8.56).6. On page 271, Eq. (9.37): kappa should be above and not below the fraction line.One major error that slipped through to the second printing affects Eqs. (5.44) and (5.45): in Eq. (5.44), on the right-hand side, rho should be replaced by R in the second and third term (superscripts l-1 and l-N_j), and similarly, in Eq. (5.45), R should be replaced by delta for superscripts larger than l. The SIP-solvers in all codes are OK; actually, the error came from looking at the coded version, since in it, the same variable is used for the residual rho, the auxiliary vector R, and the increment delta, as the order of execution of the forward and backward loops allows overwriting of one vector by another.Some spelling errors have also been found but are not listed here, since the meaning of text was not changed - the misspelling is obvious.Second Edition-Unfortunately, some errors remained undetected until the second edition was printed. Those that also existed in the second printing of the first edition are in Sect. 7.3.3 (NOTE: this section does not exist in the first printing of the first edition - the Sect. 7.3.3 there corresponds to Sect. 7.3.4 of the two latest printings; the equations with numbers given below are all OK in the first printing). The equation numbers given below refer only to the second printing of the first edition and to thesecond edition: - in Eqs. (7.22), (7.23) and (7.26), the convective terms on the RHS of the equation should bare a minus sign (see Eq. (7.17);- in Eq. (7.27), the RHS should be divided by Delta t;- in Eq. (7.28), the pressure term should be multiplied by Delta t and divided by rho.- in Eq. (9.27), the sign in front of pu_j should be + instead of -in the second edition and instead of + in the first edition (this term should be identical as in Eq. (9.26)We apologize for these typing errors (and others that might still be found) and ask kindly readers to apply corrections by hand. Thank you.The 2nd corrected printing and the second edition contain - in addition to corrected typing errors and some minor corrections to wording - also few new sections compared to the first edition. The new pages from the 2nd printing and 2nd edition are available as postscript files in the directory NEWPRINT; see readme file in this directory for further information.-Latest modification on March 21, 2001:Only today we learned that Index of the second edition is not correct -it is simply a copy of the index of the first edition! We apologize once more for this error in the preparation of the manuscript and announce that a postscript file containing the new index is included in this directory (index.ps).三、源代码及相关说明3.1 一维对流扩散方程有限差分Fortran源代码3.1.1 P1DS.F程序说明 Program for solving the steady one-dimensional convection/diffusion equation. It has been used to produce results presented in Sect. 3.11. The parameters which can be varied (by typing the values upon request during execution, or by modifying the data in an input file and re-directing input to that file) are:Density (DEN)Velocity (VEL)Diffusion coefficient (GAM)Boundary values (FI at X=X_min and FI at X=X_max)Differencing scheme for convection (UDS or CDS)Limits of solution domain (XMIN and XMAX)Grid expansion factor (EXP)Number of nodes (N)If you want to re-direct input to a file, type file after the name of the executable file. For example, if the executable file is named P1DS and you want to run the program using data from the example input file P1DS.INP, type:P1DS CDS, 2 - UDS)2 IT (1 - EE, 2 - IE, 3 - CN, 4 - I3L)0. 1. 1. 11 (X_MIN, X_MAX, EXP, N)3.1.3 P1DS.F源程序C# PROGRAM FDC#C This program solves the one-dimensional convection-C diffusion equation with Dirichlet boundary conditionsC on both ends. The exact solution is compared withC solutions obtained using FD method and either UDSC or CDS for convection term, CDS for diffusion term.C C 1995 M. Peric, Institut fuer SchiffbauC# PARAMETER (NX=81) COMMON N,NM,FI(NX),AE(NX),AW(NX),AP(NX),Q(NX),X(NX) DIMENSION FIEX(NX) CHARACTER FILOUT*10CC.OPEN FILESC PRINT *, ENTER OUTPUT FILE NAME: READ(*,1) FILOUT 1 FORMAT(A10) OPEN (UNIT=8,FILE=FILOUT)CC.READ INPUT DATAC PRINT *, ENTER: DEN(SITY), VEL(OCITY), DIF(FUSION COEFF.) READ(*,*) DEN,VEL,DIF PRINT *, ENTER BOUNDARY VALUES: FI0, FIN READ(*,*) FI0,FIN PRINT *, CHOOSE CONVECTION SCHEME: 1 - CDS, 2 - UDS READ(*,*) ICCC.DEFINE THE GRID: EX - EXPANSION FACTOR; C N - NUMBER OF NODES INCL. BOUNDARY ONESC PRINT *, ENTER: XMIN, XMAX, EX, N READ(*,*) XMIN,XMAX,EX,N NM=N-1 IF(EX.EQ.1.) THEN DX=(XMAX-XMIN)/REAL(N-1) ELSE DX=(XMAX-XMIN)*(1.-EX)/(1.-EX*(N-1) ENDIF X(1)=XMIN DO I=2,N X(I)=X(I-1)+DX DX=DX*EX END DOCC.INITIALIZE FIELDSC DO I=1,N FI(I)=0. END DO FI(1)=FI0 FI(N)=FIN DENVEL=DEN*VEL ZERO=0.C DO I=2,NMCC.CENTRAL DIFFERENCE CONVECTION APPROX. (CDS)C IF(IC.EQ.1) THEN AEC= DENVEL/(X(I+1)-X(I-1) AWC=-AECCC.UPWIND CONVECTION APPROX. (UDS)C ELSEIF(IC.EQ.2) THEN AEC= MIN(DENVEL,ZERO)/(X(I+1)-X(I) AWC=-MAX(DENVEL,ZE