INS2D Version 2.2 Incompressible Navier-stokes flow solver Author: Stuart Rogers NASA Ames Research Center Description: This code solves the incompressible Navier-Stokes equations in two-dimensional generalized coordinates for both steady-state and time varying flow. The equations are formulated into a hyperbolic set of PDE using the method of artificial compressibility. The convective terms are differenced using an upwind biased flux-difference splitting. The equations are solved using an implicit line-relaxation scheme. The code is written for single or multiple-zone calculations. It can utilize either pointwise continious zonal interfaces, or random overlapped zonal interfaces if a PEGSUS interpolation database is supplied. The flow solver contains pre-coded boundary conditions for slip and noslip walls, symmetry planes, inflow and outflow boundaries, and far-field boundaries. For more information on the algorithm, see the following technical papers, which can be accessed through the WWW URL: http://www.nas.nasa.gov/~rogers Rogers, S. E. and Kwak, D., "An Upwind Differencing Scheme for the Time Accurate Incompressible Navier-Stokes Equations," AIAA Paper 88-2583, June 1988. Published in AIAA J., 28(2), February, 1990, pp. 253--262. Rogers, S. E. and Kwak, D., "An Upwind Differencing Scheme for the Steady-state Incompressible Navier-Stokes Equations", NASA TM 101051, November 1988. Published in Journal of Applied Numerical Mathematics, 8(1), 1991, pp. 43--64. Rogers, S. E., Kwak, D., and Kiris, C., "Numerical Solution of the Incompressible Navier-Stokes Equations for Steady-State and Time-Dependent Problems, AIAA Paper 89-0463, January 1989. Also in Proceedings of the Tenth Australiasian Fluid Mechanics Conference, Melbourne, Australia, Dec. 1989. Published in AIAA J., 29(4), April 1991, pp. 603--610. Summary of recent changes to the code: -------------------------------------- - addition of curvature terms to Spalart-Allmaras turbulence model - addition of chimera ORPHAN point averaging - changes for version 2.1-27, which involve porting the code to several new hardware platforms, require the following changes to the geom.f if you are upgrading from a previous version: - change amax1 to max - change IRIS to sgi - change real to __REAL - change integer to __INTEGER - #include precis.h to any subroutine which does not include common.f Summary of Input files: ----------------------- To run the code for a specific geometry, the following files are used to uniquely define the problem: standard input: File read by code via standard input which contains namelists of input variables which control the code. See below for a complete description of these variables. Required file. bcmain.dat: defines surfaces using computational indices and specifies type of boundary conditions to be implemented on this surface. Required file. grid file: Grid file (plot3d format) containing computational grid. Required file. The first occurance of a file named x.fmt, x.dat, xy.fmt, xy.dat, xyz.fmt, xyz.dat, grid.fmt, grid.dat, g.fmt, g.dat, or g000.dat is automatically read by flow code. This can be any type of whole, 2D plot3d type grid file, either formatted or unformatted, single or multiple zone, with or without iblank array. bcpzn.dat: Boundary condition file used to define patched (pointwise continous) zonal interfaces. Optional file. bcozn.dat: Boundary conditions file used to define overlap chimera type zonal interfaces. Optional file. bcintrp.dat: Boundary condition file used to define interpolation stencil for updating (singular) boundaries by interpolating values from surrounding points in the same grid. For example, the wake boundary of a c-grid around an airfoil could be updated this way. Optional file. geom.f: Many geometries can be run without modifying any source code. The subroutines which are most commonly modified are all located in the source file geom.f. This includes user coded boundary conditions, user coded initial conditions, and output routines. File Output: The following fortran units are used for output: ------------ standard output: Writes input variable namelists and convergence history 20: Restart file output 21: Restart file input 30: Plot3d grid file output 31: Plot3d solution file output 32: Plot3d solution file containing Reynolds stress terms for turbulent calculations 33: Plot3d solution file containing functions from a turbulence model; varys with each turbulence model. 34: Plot3d solution file containing passive scalar field in the density position (fun 100). 40: filename = insaux.out, for auxiluary/debug output 50: grid file for swept infinite wing option(3D) 51: solution file for swept infinite wing option(3D) 99: scratch file for I/O processing The plot3d solution file contains fields of density, 2 components of momentum, and the internal energy per unit volume. The density field is set to unity, thus the momentum components are equal to the velocity components. The last field is computed using the formula: qp3d(4) = pressure/0.4 + 0.5*(u*u + v*v) This means that when you use function 110 (pressure) in plot3d, you get the pressure field as computed by ins2d. Beware, however, that several plot3d functions which use definitions based on a compressible equation of state and/or isentropic relation do not work: Cp, Mach number, etc. Note on file input/output for double-precision version on workstation: no change is made to reading/writing of unformatted data, so this is done using 64-bit words, with the exception of the plot3d output routines. In the double precision version, these write out 32 bit words to be compatable with 32-bit versions of post processing programs. Memory Usage: ------------- The total memory used by the program (in words) is approximatly: 75*imaxtot dynamically allocated + 3*imaxtot*(igmr+4) dynamically allocated + iunsdim*imaxtot: dynamically allocated + iswpvel*imaxtot: dynamically allocated + iperilrelax*imaxtot: dynamically allocated + nturb*imaxtot: dynamically allocated + 7*ioptot dynamically allocated + 6*idctot dynamically allocated + 11*inttot dynamically allocated + 300*jkmax: hardcoded + 34*nzne: hardcoded + 32*ibcmax: hardcoded + 180,000: overhead where imaxtot = total number of grid points ioptot = total number of chimera overlap (bcozn.dat) points idctot = total number of defect correction (bcdczn.dat) points inttot = total number of interpolation (bcintrp.dat) points nzne = parameter >= number of zones (hardcoded) jkmax = parameter >= max(jmax(nz),kmax(nz)),nz=1,nzone (hardcoded) ibcmax = parameter >= no. of bc surfaces in bcmain.dat (hardcoded) igmr = ngmrres when impsch=5, 0 otherwise iunsdim = 13 for unsteady calculations; 0 for steady-state calcs. iswpvel = 1 when swpang > 0; 0 otherwise iperilrelax = 9 when using line-relaxation on a periodic grid;0 otherwise nturb = 1 for Baldwin-Barth turbulence model = 1 for Spalart-Allmaras turbulence model = 3 for K-Omega turbulence model Description of arrays: ---------------------- Array indices: j is the index in the xi-direction k is the index in the eta-direction q(j,k,1) pressure q(j,k,2) u-component of velocity q(j,k,3) v-component of velocity s(j,k,n) storage of right-hand-side of the equations, n=1,3 dq(j,k,n) storage of delta-q terms x(j,k) x-coordinate of grid point y(j,k) y-coordinate of grid point qn(j,k,n) solution at previous time step, n=1,3 qnn(j,k,n) solution at second previous time step, n=1,3 ib(j,k) iblank array for denoting boundaries and holes rtxy(j,k,m,n) metric terms at each point, m=1,2; n=1,3 rtxy = d(xi_m)/d(x_n) / jacobian, where xi_m = xi, eta for m=1,2 x_n = t, x, y for n=1,3 dj(j,k) Jacobian of the metric transformation for current grid djn(j,k) Jacobian of the metric transformation for previous time level grid djnn(j,k) Jacobian of the metric transformation for second previous time level grid vnut(j,k) turbulent eddy viscosity turvar(j,k,2) turbulent variable used for 1 or 2 eq. turbulence models turvarn,turvarnn values of turvar from previous two time steps, used only in time-dependent integration mode. bjm, btc, bjp, bkm, bkp, bjmm, bjpp, bkmm, bkpp: These arrays hold the 3x3 blocks for the implicit side of the equations at each point in the entire flow field. For a point (j,k), the storage of the pentadiagonal system goes as follows: bjmm(j,k,m,n) stores (j-2,k) terms bjm (j,k,m,n) stores (j-1,k) terms btc (j,k,m,n) stores (j,k) terms bjp (j,k,m,n) stores (j+1,k) terms bjpp(j,k,m,n) stores (j+2,k) terms bkm (j,k,m,n) stores (j,k-1) terms bkmm(j,k,m,n) stores (j,k-2) terms bkp (j,k,m,n) stores (j,k+1) terms bkpp(j,k,m,n) stores (j,k+2) terms for m=1,3; n=1,3 scr1, scr2: Scratch memory space used for several things including storage arrays used for lu decompostion of periodic tridiagonal system. Variables in namelist DATAIN: ----------------------------- This namelist must occur once in input file. nbcimp The frequency with which the line-relaxation sweep routines will implicitly update boundary conditions. Default = 1. ndtslow Used to systematically reduce time-step size during run. Every ndtslow iterations, dtau and dt are reduced by a factor of 0.5. This should be used ONLY when running steady-state cases which have trouble converging and show oscillatory behavior usually associated with a physical (not numerical) unsteadiness. Use this with CAUTION, as small values of ndtslow will cause the time-step to become very small very quickly, making it appear that you have a steady-solution, when all you have is a time-varying solution to the INS equations "frozen" in time. Recommend values > 50. Default = 0. ngmrmax Maximum number of iterations used in gmres solver. Default = 10. ngmrres Maximum number of iterations used in gmres solver before performing a restart. Controls the amount of memory needed by the gmres solver. Default = 10. niter The frequency with which the convergence history will be written by iterout. Default = 1. nlnsrch Non-zero values of nlnsrch will cause the code to utilize a line-searching algorithm to scale the change in the mean-flow quantities. Results in between 1 and nlnsrch extra evaluations of the right-hand-side at each (sub)iteration. Default = 0. np3d The frequency with which plot3d files will be written by p3dout. If np3d=0, 2D grid and solution files are written to files fort.30 and fort.31 at the end of the calculation. For np3d>0, the grid and solution files are written to to files fortxxxx.30 and fortxxxx.31, respectively, every np3d iterations or time-steps; where xxxx=0001, 0002, etc, as determined by the current time-step. Default = 0. nqfile The frequency with which restart files will be written by diskio. nsub The maximum number of sub-iterations allowed during time-accurate calculations. ntmax The number of time steps to be run. iaxisym Turns on axisymmetric option; the y-direction is assumed to be the radial coordinate. idefect Controls input and use of defect correction for overset grids. idefect>0 causes program to read in bcdczn.dat and implement defect correction algorithm based on the contents of this file. For steady-state calculations, defect correction terms are only added after nt>idefect, where nt = iteration number. For the moving grid case (igrid > 0 and iss=0), a new defect correction file named bcdczn***.dat will be read in for each time step; see section below on Moving Grids. Default = 0. igrid Switch for an unsteady grid motion. If igrid=0, the grid is only read in at start of exectution, and the grid is stationary. If igrid > 0 and iss = 0, a time-varying sequence of grids name g***.dat will be read in, one before each time step. See section below on moving grids. Default=0. iimpio Controls output of implicit scheme. For line-relaxation, non-zero value causes the error to be evaluated after each sweep and printed. Default = 0. iintrp Controls input of interpolation boundary condition info. If 0, no single zone interpolation is used. If non-zero, code will attempt to read in file named bcintrp.dat. boundary. Default = 0. imatvec Switch for gmres process which controls how matrix-vector products will be evaluated. imatvec=0 uses analytic approximate Jacobian LHS matrix. imatvec=1 uses frechet numerical derivative. Default = 0. impsch Switch for implicit scheme: 1: line relaxation 2: point relaxation 3: ilu factorization 4: lusgs factorization 5: gmres solver: generalized minimum residual Default = 5. iover Controls input of overlap boundary condition info. If iover=0, there are no overlap (Chimera) boundaries. If iover=1, the code will attempt to read in a file named bcozn.dat containing overlap interpolation info. If iover>0 and iss=0, a new overlap boundary condition file named bcozn***.dat will be read in for each time step; see section below on moving grids. See section below for file format of bcozn.dat file. Default=0. ipatch Controls input of patch boundary condition info. If 0, no patched boundaries. If ipatch=1, code will attempt to read in file named bcpzn.dat containing patch boundary info. If ipatch>0 and iss=0, a new patch boundary condition file named bcpzn***.dat will be read in for each time step; see section below on Moving Grids. See section below for file format bcpzn.dat files. Default=0. iprec Selects preconditioner gmres implicit solver. 1: line relaxation 2: point relaxation 3: ilu factorization 4: lusgs factorization Default = 0. iplot Controls the frequency of the calling of the plot routine, used for graphical interface to iris workstation. Default=0 iscalar Switch to control calling of passive scalar solver at end of computation. Call solver if iscalar.ne.0. Default=0 iss Switch for steady-state (iss=1), or time-accurate (iss=0) calculations. istart Switch which controls whether the code is started from initial conditions (istart=0), or whether from a restart file to be read in by diskio (istart=1 or 2). If istart=2, time and nt are reset to 0 at the restart. istart=3 denotes a run at a different angle of attack than the previous run, and the velocity components from the previous run are rotated to the new angle of attack as specified in the input. isubio Switch which controls the output of convergence history during the subiterations. If isubio=0, output is written only at end of subiterations. Else it is written for every subiteration. iturb Switch to choose the turbulence model used. Baldwin-Barth(2), Spalart-Allmaras(3), K-Omega(4), Durbin-Mansour(6) Default=3 iupsch Switch to choose the type of differencing, flux-difference splitting(1), or flux-vector splitting(2). Default = 1. ivis Switch which controls the type of viscous fluxes used. ivis = 1: constant viscosity and orthogonal grid 2: constant viscosity and non-orthogonal grid 3: variable viscosity and orthogonal grid 4: variable viscosity and non-orthogonal grid beta The artificial compressibility parameter. See notes at the end of this document on selecting beta and dtau. dtau The time step in pseudo-time. Default = 1.e12. See notes at the end of this document on selecting beta and dtau. dt The time step in real-time, used during time- accurate calculations. For steady-state turbulent calculations, most turbulence models use this as their time step. eigeps The epsilon used in forming the plus/minus eigenvalues. This varible will have the effect of smoothing out problems which might occur where the contravarient velocity changes sign. Default = 0.1. epscon Used as a convergence critera. For time-dependent cases (iss=0) this is the maximum allowable value of the maximum divergence of velocity during the sub-iterations; when this is reached, the subiterations are terminated and the next time-step is begun. For steady-state cases (iss=1) this determines the maximum divergence of velocity allowed for a converged steady-state solution; when the maximum divergence of velocity is less than epscon, the calculations are stopped. Default = 1.e-4. flxlim Controls the order of the explicit numerical fluxes, 0.0 gives uniform third-order fluxes, 1.0 gives first-order. Values inbetween give a linear combination of the two. Using small values (0.01) can add a little dissipation to stabilize the equations, using larger values (0.1 - 1.0) for the first 50 or so iterations provides very good initial conditions, and extremely fast convergence. The accuracy of the final solution, however, might be seriously degraded if any value greater than 0.0 is used at the end of the computation. Default = 0.0 gmrtol Tolerance for residual ||A*dq - b|| when solving the system of equations using gmres. Default = 1.e-6 pin Total pressure used at inflow boundaries. Default = 1.5 pout Static pressure used at outflow boundaries. Default = 1.0 Both pin and pout are the dimensionless values of pressure as used by the code: pressure is nondimensionalized by (p - p_ref)/(2*dynamic pressure); dynamic pressure = 1/2*rho*U_ref*U_ref. p_ref is arbitrary, but by default is selected such that the inflow/far-field pressure is 1.0. preleft Logical variable which selects left or right preconditioning for the gmres algorithm. Default = .true. reynum The Reynolds number based on the units of length and velocity input to the code: If you run the code with the inflow/farfield velocity scaled to unity, give the Reynolds number based on your inflow/farfield velocity. If you set the inflow velocity in the code in some dimensional units (i.e. m/sec), give the Reynolds number based on that unit (i.e. 1 m/sec.) If you non- dimensionalize the coordinates of the grid such that the characteristic length is unity, give the Reynolds number based on that length. If the grid coordinates are given in dimensional units (i.e. cm), give the Reynolds number based on that unit (i.e. Reynolds number per cm). swpang Sets sweep angle for swept infinite wing option. Non-zero values causes the code to solve the uncoupled z-momentum equation assuming zero gradients in the z-direction, where the z-axis is aligned with the leading-edge of an infinite wing. and the xy-plane (containing the input grid) is aligned perpendicular to the leading-edge. The input Reynolds number should be that based upon the chord in this xy-plane, and the actual freestream velocity. This way, the grid and the Reynolds number do not change with swpang for a given wing. The calculation produces a solution to the w-component of velocity which couples into the rest of the flow variables only through the eddy-viscosity. Currently only coded for the Baldwin-Barth turbulence model. Default = 0. underr Variable which enables under-relaxation to be used in the line-relaxation updates. Set to a value 0 < underr < 1.0 for under-relaxation. Helps reduce instabilities. Values greater than 1.0 (over-relaxation) almost always cause the code to become unstable. Default=1.0 Variables in namelist ZONEIN: ----------------------------- This namelist must occur once for every zone. kpr Logical variable which determines if the grid is periodic in the k-direction for this zone. The k=1 and k=kmax line should be neighbors and not be coincident. njswp,nkswp Number of sweeps used by the line- and point-relaxation algorithms. For line-relaxation, 1 j-sweep is either a forward or a backward sweep through the entire domain solving the equations along a line of constant k; a k-sweep solves lines of contant j. Line-relaxation works best when using constant lines which are normal to a no-slip surface; thus for body normal grid lines in the eta-direction, use nkswp>0, njswp=0. For point-relaxation, the number of sweeps used is equal to max(njswp, nkswp); where a point-relaxation sweep is a forward+backward sweep through the entire domain. ntjswp Similar to njswp,nkswp except these control the relaxation schemes when solving the turbulence model. iflxodr Selects first-order (iflxodr=1), third-order (iflxodr=3), or fifth-order (iflxodr=5) stencil for flux-differencing splitting of convective terms. Default = 3. itran Controls handline of transition in turbulence model. Effect is dependent on which turbulence model is used. cdiss Coeffient of dissipation used to multiply the dissipative terms in the flux-difference splitting scheme. Default = 1.0 xtr1,xtr2 Given to turbulence model for specifing transition for airfoils with c-grids. Grid file --------- The code accepts PLOT3D 2D whole files. These can be single- or multi-zone, formatted or unformatted, with or without iblank arrays. The flow solver will automatically determine which type of file you have. The following are the appropriate read statments for the different types of files, for a Fortran unformatted file. The Fortran formatted file type uses identical records except with a wildcard format specification. Single zone, no iblank array: read(unit) jmax,kmax read(unit) ((x(j,k),j=1,jmax),k=1,kmax), & ((y(j,k),j=1,jmax),k=1,kmax) Single zone, with iblank array: read(unit) jmax,kmax read(unit) ((x(j,k),j=1,jmax),k=1,kmax), & ((y(j,k),j=1,jmax),k=1,kmax), & ((ib(j,k),j=1,jmax),k=1,kmax) Multiple zone, no iblank array: read(unit) nzone read(unit) (jmax(nz),kmax(nz),nz=1,nzone) do 10 nz=1,nzone read(unit) ((x(j,k),j=1,jmax(nz)),k=1,kmax(nz)), & ((y(j,k),j=1,jmax(nz)),k=1,kmax(nz)) 10 continue Multiple zone, with iblank array: read(unit) nzone read(unit) (jmax(nz),kmax(nz),nz=1,nzone) do 10 nz=1,nzone read(unit) ((x(j,k),j=1,jmax(nz)),k=1,kmax(nz)), & ((y(j,k),j=1,jmax(nz)),k=1,kmax(nz)), & ((ib(j,k),j=1,jmax(nz)),k=1,kmax(nz)) 10 continue bcmain.dat: boundary condition input file ----------------------------------------- Each line of input in this file should include values for ibcval, izone, jbeg,jend, kbeg,kend. The last 5 integers are a set of indices which describe a computational surface. Negative values of jbeg, jend, kbeg, or kend are interpreted as: jbeg = jmax + 1 - abs(jbeg). Thus jbeg=-1 is equivalent to jbeg=jmax; jbeg=-2 is equivalent to jbeg=jmax-1. The value of ibcval specifies which boundary condition will be applied to that surface. Any surfaces on the computational boundaries not specified here will by default remain at their initial condition values. Only the wall boundary conditions can be implemented on interior points. Comments are allowed anywhere in the file using a "c" in the first column. Tab characters are allowed in the input records. NOTE: The no-slip wall and slip wall boundary conditions can be applied at any computational surface, interior or boundary. If a wall boundary condition is to be specified in an interior plane, it is important that the wall be at least 3 grid-planes thick. Two of these planes would be needed to specify the wall faces, and a third plane sandwiched inbetween the walls which must be blanked out with a negative value for ibcval. This is necessary to prevent the differencing stencil from performing a differencing "through" the wall. ibcval Description ------ ----------- 0 No-slip wall, with wall normal vector pointing in the positive computational direction 1 No-slip wall, with wall normal vector pointing in the negative computational direction 4 Moving no-slip wall, with wall normal vector pointing in the positive computational direction. Pressure is updated by specifying zero pressure gradient in the wall normal direction. If igrid=0, no change in wall velocity is specified (maintains initial conditions). If igrid>0, wall velocity is set equal to grid velocity. 5 Moving no-slip wall, with wall normal vector pointing in the negative computational direction. Pressure is updated by specifying zero pressure gradient in the wall normal direction. If igrid=0, no change in wall velocity is specified (maintains initial conditions). If igrid>0, wall velocity is set equal to grid velocity. 10 Slip wall with the wall normal pointing in the positive computational direction 11 Slip wall with the wall normal pointing in the negative computational direction 20 Inflow boundary with constant velocity and extrapolation of pressure 21 Inflow boundary with total pressure held constant and everywhere equal to pin (namelist input variable); and flow normal to the surface, extrapolate pressure 22 Inflow boundary with total pressure locally held constant as specified in initial conditions; flow held normal to surface, and extrapolate pressure. 23 Inflow boundary with total pressure locally held constant as specified in initial conditions; flow held normal to surface, and extrapolate normal velocity component. 25 C-grid outer boundary for airfoils which imposes velocity field based on freestream velocity at current angle of attack superimposed with a velocity field from a point vortex centered at (0.25,0) (quarter chord). Pressure is computed using a characteristic relation. The strength of the vortex is based on the current value of lift coefficient stored in the global variable clift, which should be computed and updated after every iteration. This is usually done in subroutine iterout. The routines also use the variable alpha to determine the incidence of the free-stream. 26 Same as ibcval=25, except that pressure is extrapolated. 30 Outflow boundary using characteristic relations for velocity and constant static pressure (pout). 31 Outflow boundary using extrapolated velocity and constant static pressure (pout). 32 Outflow boundary using characteristic relations for velocity and constant (initial condition) pressure. 33 Outflow boundary using extrapolated velocity and constant (initial condition) static pressure. 40 Symmetry boundary conditions: currently only works for a j-surface if the xi-computational direction is aligned with the x-axis. Similar for k surfaces. 50 Far-field boundaries which will have inflow in one area and outflow in another, such as in an o-grid. For points with inflow, it uses ibcval = 20 For points with outflow, it uses ibcval = 31 60 c-grid wake cut boundary: these points will be updated by averaging values from surrounding points. Normal surface vector points in positive computational direction. 61 c-grid wake cut boundary: these points will be updated by averaging values from surrounding points. Normal surface vector points in negative computational direction. 100 User coded boundary condition: code calls subroutine bcuser before beginning of implicit solution. In this routine code can be written to specify values in the dq array. Useful for implicitly specifying analytic or known conditions. < 0 Any value less than zero indicates that this region is to be blanked out, such that the code will never use any values at those points. This region does not have to be a surface, it can be a volume. bcpzn.dat: Patched grid interface input file -------------------------------------------- The bcpzn.dat file is used by INS2D to update grid interface boundaries on patched grids. The interface requires pointwise continuity at the boundaries between the grids, and an overlap of points such that the boundary points in one grid lie on top of interior points of another grid. Each interface consists of two input lines: one containing indices describing a computational line of points in the target zone, another describing a computational line of points in the base zone. The values of the dependent variables from the base zone points are used to update the dependent variables in the target zone points. The target line is listed first, the base line is listed second. Each line lists the zone number, the beginning, ending, and increment for the j-index, and the beginning, ending, and increment for the k-index. nz jbeg jend jinc kbeg kend kinc 1 1 11 1 20 20 1 2 2 2 1 31 21 -1 In this example, the k=20 line in zone 1 from j=1,11 will be updated with values from a j=2 line in zone 2. Note that the zone 2 k-direction aligns with negative j-direction in zone 1. Also note that for each set of indices, you must have either jbeg=jend, or kbeg=kend. bcozn.dat: Overlaid grid interface input file --------------------------------------------- The bcozn.dat file is used by INS2D to update grid interface boundaries on overlaid grids, also known as chimera grids. The bcozn.dat file contains interpolation weights and grid indices. Typically the interpolation information is generated by the PEGSUS software, or a similar program. The interpolation output from PEGSUS needs to be converted to the bcozn.dat format. This can be done using the utility program pegconv.f. The interpolation is used to update points on either an outer or fringe boundaries. These points are refered to as the TARGET points. Each target point has a corresponding BASE grid cell, usually in a different grid and which surrounds the target point. In two-dimensional grids, the four grid points which comprise the base grid cell are the base points. Values at the base points are interpolated and injected into the target point. The bcozn.dat file is in ascii format (Fortran formatted). The data is organized as a list of the base points found in each zone, and their corresponding target point indices. The file is composed of the following records: do 10 nz=1,nzone read(*,*) ibpnts(nz) read(*,*) (jbase(nz,i),kbase(nz,i), & dxb(nz,i),dyb(nz,i), & nztar(nz,i),jtar(nz,i),ktar(nz,i),i=1,ibpnts(nz)) 10 continue where ibpnts(nz) is the number of base points in zone nz; jbase and kbase are the indices of the lower left corner (in computational space) of the base grid cell in zone nz; dxb and dyb are the interpolation weights; jtar and ktar are the indices of the target point in zone nztar. The interpolation formula used to update points is best given through the following example. Also, see subroutine bcozndq and subroutine bcoznds for other examples. In this example, b denotes the zone containing the base point and t denotes the zone containing the target point. real qb(jmaxb,kmaxb,3), qt(jmaxt,kmaxt,3) j = jbase(i,nzb) k = kbase(i,nzb) a1 = qb(j,k,n) a2 =-qb(j,k,n) + qb(j+1,k,n) a3 =-qb(j,k,n) + qb(j,k+1,n) a4 = qb(j,k,n) - qb(j+1,k,n) - qb(j,k+1,n) + qb(j+1,k+1,n) jt = jtar(nzb,i) kt = ktar(nzb,i) qt(jt,kt,n) = a1 +a2*dxb(i) +a3*dyb(i) +a4*dxb(i)*dyb(i) bcdczn.dat: Defect correction on overset grids, interface input file --------------------------------------------------------------------- bcintrp.dat: Interpolation input file ------------------------------------- Moving Grids ------------ The code will handle moving grids during an unsteady calculation. In this case, the value of igrid is to be set greater than zero. This will cause the code to read in a new grid file at the beginning of each time-step. The grid files must be named gxxx.dat, where xxx is 000, 001, 002, 003, ... , nnn. Set igrid equal to the number nnn, which is the last number in this sequence of grid file names. If igrid < ntmax, then the code assumes that the motion is periodic in time. In this case, after using the grid file gnnn.dat, the next time step uses the grid file g001.dat. For a periodic motion, the grid g000.dat should be an exact copy of gnnn.dat, where nnn=igrid. Restarts should work OK at any point in this sequence. Any no-slip surfaces which are moving in their grids should use boundary conditions of ibcval=4 or 5. The velocity boundary condition for these walls is specified to match the velocity of the sequential motion of the grid. For example, the velocity on the moving wall will be set to the distance a boundary point has moved from one time-step to another divided by the physical time step dt. The zones in a sequence of grid files may move relative to one-another. In this case a different set of zonal boundary conditions must be applied for each time step. For overset zones, set iover>1 and provide files named bcozn000.dat, bcozn001.dat, etc. For patched zones, set ipatch>1 and provide files named bcpzn000.dat, bcpzn001.dat, etc. For defect correction overlaid zones, set idefect>1, and provide files named bcdczn000.dat, bcdczn001.dat, etc. Typically one would have the same number of bc files as grid files in the sequence (igrid=iover or igrid=ipatch). Additional Notes: (Or, what to do to make the code run.) ---------------------------------------------------------- Steady-state convergence: One of the most important aspects in make the code converge is to select appropriate values for beta and dtau. If stability is not a problem, keep dtau very large. Otherwise, reduce dtau to a value on the order of 1.0. This assumes that the velocity and coordintes have been scaled by the characteristic velocity and lengths, and they are unity. Thus, dtau=1.0 will convect along a characteristic length in one time-step. If dtau is reduced a lot less than this, it will really slow down the convergence, as it will take many time-steps to propagate the disturbances away from the body or out of the flow domain. So keep dtau near 1.0, and try to stablize by selecting beta: Typical values for beta are 1 to 100 depending on geometry, grid, and reynolds number (assuming you non-dimensionalize your grid with characteristic length). Internal flows tend to do better in the lower range, but turbulent and fine-grid cases tend to require larger beta. Try running a couple of short runs with beta=1,10,100, to see what order of magnitude it likes for your case.
Updated: Wednesday, 18-Oct-2000 06:52:04 PDT
WebWork: Stuart E. Rogers <rogers@nas.nasa.gov> NASA Responsible Official: Stuart Rogers |