INS3D-UP Version 2.1 Incompressible Navier-stokes flow solver Author: Stuart E. Rogers NASA Ames Research Center Description: This code solves the incompressible Naiver-Stokes equations in three-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, or with a gmres (generalized minimum residual) iterative solver. 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: -------------------------------------- - changed handling of orphan points, no longer need file ORPHAN; instead check for iblank=101 values in input grid file if iorphan=1 (the default value) and the input grid has an iblank array. - added periodic o-grid averaging boundary conditions. - out-of-core run-time option with impsch=1 (line-relaxation): incore=0 causes code to buffer zones out of memory, resulting in possibly significantly less memory usage with about a 5% CPU increase. incore=0 should produce exact same results and convergence as incore=1. To upgrade to this version from previous versions, you will have to add some code to subroutine iterout in your geom.f file; see the geom.f file that came with the release. - For overset grids, now use INTOUT file instead of bcozn.dat file; no longer use bcpzn.dat file. - Checks for existence of XINTOUT as well as INTOUT to read in overset grid boundary conditions. - 2nd-order accurate time-integration in both turbulence models for unsteady flows. - Various bug fixes, mostly for Spalart-Almaras turbulence model. - Addition of GMRES implicit solver, impsch=5. This is faster per iteration, faster convergence, more robust, but requires a lot more memory. - Changed inflow boundary conditions for pressure: replace characteristic relation with extrapolation. Found that the characteristic relation had a fairly strong dependence on beta. - Add impsch=6: GMRES version which solves only 1 zone at a time, as opposed to impsch=5 which solves entire domain simultaneously. This will converge as well as impsch=5 for some cases, and for others it converge quite a bit slower. It is still better than impsch=1. With this option only is it feasible to run with an out-of-core option. - out-of-core run-time option with impsch=6 (gmres): incore=0 causes code to buffer zones out of memory, resulting in possibly significantly less memory usage with about a 20-40% CPU increase. incore=0 should produce exact same results and convergence as incore=1. To upgrade to this version from previous versions, you will have to add some code to subroutine iterout in your geom.f file; see the geom.f file that came with the release. - Added periodic O-grid averaging boundary condition Summary of Input files: ----------------------- The following files are used to uniquely define the flow problem: Standard Input: File read by code via standard input which contains namelists of input variables which control the code. See below for description of the namelists and their 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 grid.in is auto- matically read by flow code. This can be any type of whole, 3D plot3d type grid file, either formatted or unformatted, single or multiple zone, with or without iblank array. INTOUT or XINTOUT: Boundary conditions file used to define overlap chimera type zonal interfaces. 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 The plot3d solution file contains fields of density, 3 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(5) = pressure/0.4 + 0.5*(u*u + v*v + w*w) This means that when you use function 110 (pressure) in plot3d, you get the pressure field as computed by ins3d. 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: ------------- For incore, steady problems, the memory usage in the code is approximately 40*imaxtot words for impsch=1, or 220*imaxtot words for impsch=5 or 6, where imaxtot = total number of grid points. For incore=0, replace imaxtot with imaxnz, where imaxnz = number of grid points in the largest zone. For unsteady problems, add 12*imaxtot (incore=1) or 12*imaxnz (incore=0). More precisely, the total memory used by the program (in words) is: with impsch=1, incore=1: 34*imaxtot dynamically allocated + 12*imaxtot*iuns dynamically allocated + imaxtot*nturb dynamically allocated + 10*ioptot dynamically allocated + (200 + 52*lenvec)*jklmax hardcoded + 47*nzne: hardcoded + 22*ibcmax: hardcoded + 210,000: overhead with impsch=1, incore=0: 34*imaxnz dynamically allocated + 10*imaxnz*iuns dynamically allocated + imaxnz*nturb dynamically allocated + 10*ioptot dynamically allocated + (200 + 52*lenvec)*jklmax hardcoded + 47*nzne: hardcoded + 22*ibcmax: hardcoded + 210,000: overhead with impsch=5 or 6, incore=1: 162*imaxtot dynamically allocated + 4*(ngmrres+4)*imaxtot dynamically allocated + 12*imaxtot*iuns dynamically allocated + imaxtot*nturb dynamically allocated + 10*ioptot dynamically allocated + (200 + 52*lenvec)*jklmax hardcoded + 47*nzne: hardcoded + 22*ibcmax: hardcoded + 210,000: overhead with impsch=6, incore=0: 162*imaxnz dynamically allocated + 4*(ngmrres+4)*imaxnz dynamically allocated + 12*imaxnz*iuns dynamically allocated + imaxnz*nturb dynamically allocated + 10*ioptot dynamically allocated + (200 + 52*lenvec)*jklmax hardcoded + 47*nzne: hardcoded + 22*ibcmax: hardcoded + 210,000: overhead where imaxtot = total number of grid points imaxnz = number of grid points in largest zone iuns = 0 for steady-state calculations = 1 for unsteady calculations nturb = 4 for Baldwin-Barth turbulence model = 1 for Spalart-Almaras turbulence model ioptot = total number of chimera overlap (INTOUT) points jklmax = parameter >= max(jmax(nz),kmax(nz),lmax(nz)),nz=1,nzone nzne = parameter >= number of zones ibcmax = parameter >= no. of bc surfaces in bcmain.dat lenvec = parameter: max. length of vector used during vectorized block-tridiagonal solution process Description of arrays: ---------------------- Array indices: j is the index in the xi-direction. k is the index in the eta-direction. l is the index in the zeta-direction. q(j,k,l,1) pressure q(j,k,l,2) u-component of velocity q(j,k,l,3) v-component of velocity q(j,k,l,4) w-component of velocity s(j,k,l,n) storage of right-hand-side of the equations, n=1,4 dq(j,k,l,n) storage of delta-q terms x(j,k,l) x-coordinate of grid point y(j,k,l) y-coordinate of grid point qn(j,k,l,n) solution at previous time step, n=1,4 qnn(j,k,l,n) solution at second previous time step, n=1,4 rtxyz(j,k,l,m,n)metric terms at each point, m=1,3; n=1,4 rtxyz = d(xi_m)/d(x_n) / jacobian, where xi_m = xi, eta, zeta for m=1,2,3 x_n = t, x, y, z, for n=1,4 dj(j,k,l) jacobian of the metric transformation vnut(j,k,l) turbulent eddy viscosity turvar(j,k,l) field variables used for various turbulence models bm,bd,bp,ff: These arrays hold the 4x4 blocks for the implicit side of the equations at each point in a line. Variables in namelist DATAIN: ----------------------------- This namelist must occur once in standard input file. 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. ntmax The maximum number of physical time steps or iterations to be performed. Setting ntmax=0 causes the code to stop after reading in grid, setting initial conditions or reading the restart, and writing out plotting files. Setting ntmax=-1 causes the code to read in the grid and compute and save the distance function used by the turbulence model. np3d The frequency with which plot3d files will be written by p3dout. Default = 100. nqfile The frequency with which a restart file will be written by diskio. Default = 100. niter The frequency with which the convergence history will be written by iterout. Default = 1. nsub The maximum number of sub-iterations that can be used during one physical time step. This should be set to 1 for steady-state computations, and should be larger than 1 for time-dependent calculations. Default = 1. igrid Switch which controls if the grid is regenerated at each time step(igrid = 1) or if only at begining of run(=0). Default = 0. impsch Switch for implicit scheme: 1: line relaxation 3: ilu factorization 4: lusgs factorization 5: gmres solver: generalized minimum residual solving the entire solution at once. 6: gmres solver: generalized minimum residual solving one zone at a time. The only reason to run this instead of impsch=5 is when you need to run with incore=0. Default = 1. incore Out-of-core option for multiple zone grids: setting incore=0 causes the code to keep only 1 zone in core memory at a time; the other zones are written to disk. Otherwise, all zones are kept in core memory at all times. Can signficantly reduce memory requirements while increasing CPU and turnaround time due to time spent doing disk I/O. Using incore=0 requires about 60*imaxtot words in temporary disk space, where imaxtot is the total number of grid points. This option is not available with impsch=5. Default = 1. iorphan Non-zero values of iorphan cause the code look for orphan points which have been designated in the input grid file with an iblank value of 101. This iblank value is the default behavior for orphan points in the PEGSUS program, as written to the file IBPLOT, and processed by the utility program mergepeg, which creates the composite grid file. The orphan points are updated in the flow code by averaging all surrounding points. Default = 1. iover Controls input of overlap boundary condition info. If 0, no overlap (Chimera) boundaries. If non-zero, code will attempt to read in file named either XINTOUT or INTOUT, containing overlap interpolation info. Default = 0. iprec Selects preconditioner for the gmres implicit solver. 0: no preconditioning 3: ilu factorization 4: lusgs factorization Default = 3. iss Switch used to control steady-state runs (iss=1) or for time-accurate calculations (iss=0). Default = 1. istart Switch which determines if initial conditions are to be used (istart=0), or if a restart file is to be read in by diskio (istart=1 or 2). If istart=2, time and nt are reset to zero. Default = 0. 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. Default = 1. iturb Switch to choose the turbulence model used: 2 = Baldwin-Barth one-equation model 3 = Spalart-Allmaras one-equation model Default=3. ivis Determines what assumptions are to be used in calculating the viscous fluxes. ivis = 1: Constant viscosity and orthogonal grid. 2: Constant viscosity and non-orthogonal grid. 3: Variable viscosity and non-orthogonal grid. ivis needs to be set to 3 for turbulent calculations. beta The artificial compressibility parameter. See discussion near end of document on choosing this parameter. dtau The time step in pseudo-time. dt The time step in real-time, used only in time-accurate calculations. eigeps The epsilon used in forming the plus/minus eigenvalues. This value should not be changed. 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. gmrtol Tolerance for residual ||A*dq - b|| when solving the system of equations using gmres. Default = 1.e-2 omegax The non-dimensional rate of rotation about the x-axis for steady-rotating reference frame geometries. For omegax .ne. 0., the code includes source terms which account for the centrifugal and coliolis forces. In this case, the velocity component variables in the code are the relative velocity components in the rotating reference frame; the velocities are transformed into the absolute frame of reference before being used by the turbulence models. The non-dimensional value of omega is obtained from the dimensional value (in radians/[sec,min,etc]) and multiplying by the characteristic length and dividing by the characteristic velocity. Example: rotation = 10,000 RPM = 1047 rad/sec; U_infinity = 150 ft/sec = 1800 in/sec; length=1 in; omegax = (1047 rad/sec)*(1 in)/(1800 in/sec) = 0.582 Default = 0. omegay The non-dimensional rate of rotation about the y-axis for steady-rotating reference frame geometries. Default = 0. omegaz The non-dimensional rate of rotation about the z-axis for steady-rotating reference frame geometries. Default = 0. pin Total pressure used at inflow boundaries. Default = 1.5 pout Static pressure used at outflow boundaries. Default = 1.0 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). 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 Number of times the line-relaxation algorithm sweeps the entire domain solving the equations along a line of constant k and l. nkswp Number of times the line-relaxation algorithm sweeps the entire domain solving the equations along a line of constant j and l. nlswp Number of times the line-relaxation algorithm sweeps the entire domain solving the equations along a line of constant j and k. ntjswp, ntkswp, ntlswp Number of line-relaxation sweeps used by the turbulence model. iaxsym Specifies that a singular axis boundary lies in a symmetry plane. A non-zero value of iaxsym causes the boundary conditions for the singular axis in this zone to zero out one or more velocity components, to enforce the symmetry condition. The value of iaxsym should be a 3 digit binary number, where the left digit corresponds to the u-component, the middle digit to the v-velocity component, and the right to the w-component. If the digit is zero, that velocity component is zeroed out, otherwise it retains the value computed with the axial singularity routines. The axial boundary condition must be specified for this parameter to work: see ibcval=70,71,72. Default=111. iflxodr Selects third-order (iflxodr=3) or fifth-order (iflxodr=5) stencil for flux-differencing splitting of convective terms. Default = 3. cdiss Coeffient of dissipation used to multiply the dissipative terms in the flux-difference splitting scheme. Default = 1.0 It is strongly recommended that this is never changed. Grid file --------- The code accepts PLOT3D 3D 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,lmax read(unit) (((x(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax), & (((y(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax), & (((z(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax) Single zone, with iblank array: read(unit) jmax,kmax,lmax read(unit) (((x(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax), & (((y(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax), & (((z(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax), & (((ib(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax) Multiple zone, no iblank array: read(unit) nzone read(unit) (jmax(nz),kmax(nz),lmax(nz),nz=1,nzone) do 10 nz=1,nzone read(unit) (((x(j,k,l),j=1,jmax(nz)),k=1,kmax(nz)),l=1,lmax(nz)), & (((y(j,k,l),j=1,jmax(nz)),k=1,kmax(nz)),l=1,lmax(nz)), & (((z(j,k,l),j=1,jmax(nz)),k=1,kmax(nz)),l=1,lmax(nz)) 10 continue Multiple zone, with iblank array: read(unit) nzone read(unit) (jmax(nz),kmax(nz),lmax(nz),nz=1,nzone) do 10 nz=1,nzone read(unit) (((x(j,k,l),j=1,jmax(nz)),k=1,kmax(nz)),l=1,lmax(nz)), & (((y(j,k,l),j=1,jmax(nz)),k=1,kmax(nz)),l=1,lmax(nz)), & (((z(j,k,l),j=1,jmax(nz)),k=1,kmax(nz)),l=1,lmax(nz)), & (((ib(j,k,l),j=1,jmax(nz)),k=1,kmax(nz)),l=1,lmax(nz)) 10 continue bcmain.dat: boundary condition input file ----------------------------------------- Each line of input in this file should include values for ibcval, izone, beg,jend, kbeg,kend, lbeg,lend. The value of ibcval specifies which boundary condition will be applied, the remaining indices describe the zone number and j,k,l range which describes a computational surface where the boundary condition will be applied. Any surfaces on the computational boundaries not specified here will by default remain at their initial condition values. 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. 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, slip wall, inflow, and outflow 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. The initial condition velocity is maintained. 5 Moving no-slip wall, with wall normal vector pointing in the negative computational direction. The initial condition velocity is maintained. 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 the pressure 21 Inflow boundary with constant total pressure (pin) and flow normal to the surface, and extrapolation of the pressure 22 Inflow boundary with constant total pressure (initial conditions) and flow normal to the surface, extrapolation of the pressure 30 Outflow boundary using characteristic relations for velocity and uniform constant static pressure=pout. 31 Outflow boundary using extrapolated velocity and uniform constant static pressure=pout. 32 Outflow boundary using characteristic relations for velocity and constant static pressure as set by the initial conditions. 33 Outflow boundary using extrapolated velocity and constant static pressure as set by the initial conditions. 40 Symmetry boundary conditions for x-constant symmetry plane. Cannot be applied at interior planes. 41 Symmetry boundary conditions for y-constant symmetry plane Cannot be applied at interior planes. 42 Symmetry boundary conditions for z-constant symmetry plane Cannot be applied at interior planes. 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, this uses ibcval = 20 For points with outflow, this uses ibcval = 31 60-63 C-grid flow-through conditions: for a boundary grid plane with a c-topology in which the two halves of the c collapse onto the same plane. These boundary points are updated by averaging values above and below the plane. Only one half of the c-grid plane needs to be given as a boundary condition, both sides will be updated. ibcval=60,61 are used when the vector normal to the plane points in the positive computational direction (ie, j=1 plane). Use ibcval=62 or 63 when this vector points in the negative computational direction (ie, j=jmax plane). Using indicinal notation to represent j,k,l: when the boundary plane is in the i-plane and the c-direction is in the i+1 computational direction, use ibcval=60 or 62; when the c-direction is in the i-1 computational direction, use ibcval=61 or 63. To summarize: Boundary c-direction ibcval plane j=1 k 60 j=1 l 61 j=jmax k 62 j=jmax l 63 k=1 l 60 k=1 j 61 k=kmax l 62 k=kmax j 63 l=1 j 60 l=1 k 61 l=lmax j 62 l=lmax k 63 70,71,72 Singular axis boundary condition: where a computational boundary plane collapses on to a singular line. The value of ibcval identifies which computational direction wraps circumferentially around the singular line: 70=j-index, 71=k-index, 72=l-index. If the axis lies in a symmetry plane, must change the input value of iaxsym in namelist ZONEIN for this zone. 80 Periodic O-grid averaging (nonconservative). With an o-grid periodicity in the j-direction, the j=1 and j=jmax planes should be pointwise coincident. Specifiy on one of the two boundary planes. Updates both planes by averaging the values on either side of the periodic plane. This boundary condition will lead to inaccuracies in the solution, and should only be used as a last resort. 81-89 Reserved for spatially periodic BC (defered implementation) 90 Extrapolate all four variables from the neighboring surface in positive computational direciton. 91 Extrapolate all four variables from the neighboring surface in negative computational direciton. 100 User coded boundary condition: code calls subroutine bcuser during 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. INTOUT and XINTOUT: Overlaid grid interface input file ------------------------------------------------------ The INTOUT and XINTOUT files are used by INS3D to update grid interface boundaries on overlaid grids, also known as chimera grids. This file contains interpolation weights and grid indices. It is the standard interpolation file output from both the PEGSUS and DCF software packages. The interpolation is used to update points on either an outer boundary or a fringe boundary. The following shows the contents of the INTOUT file. This set of 4 records are repeated once for each zone in the grid: read(2) ibpnts(nz),iipnts(nz),iieptr(nz),iisptr(nz) read(2) (ji(i),ki(i),li(i),dj(i),dk(i),dl(i),i=1,iipnts(nz)) read(2) (jb(i),kb(i),lb(i),ibc(i),i=1,ibpnts(nz)) read(2) (iblank(i),i=1,jmax(nz)*kmax(nz)*lmax(nz) The XINTOUT file is a variation of the INTOUT file but with a different ordering of the data which results in significantly faster I/O. The code first checks for the existence of XINTOUT; if it exists, it will not use an INTOUT file. For each zone, the XINTOUT file should have the following four records: read(2) ibpnts(nz),iipnts(nz),iieptr(nz),iisptr(nz) read(2) (ji(i),i=1,iipnts), & (ki(i),i=1,iipnts), & (li(i),i=1,iipnts), & (dj(i),i=1,iipnts), & (dk(i),i=1,iipnts), & (dl(i),i=1,iipnts) read(2) (jb(i),i=1,ibpnts), & (kb(i),i=1,ibpnts), & (lb(i),i=1,ibpnts), & (ibc(i),i=1,ibpnts), read(2) (iblank(i),i=1,jmax(nz)*kmax(nz)*lmax(nz) where ji,ki,li are the indices of the interpolation cell which will be used to supply information to the boundary points; dj,dk,dl are the interpolation weights; jb,kb,lb are the indices of the boundary points which receive the new information, and iblank is an integer array which has a value of 1 at active field points, and a value of zero at hole and boundary points. iisptr and iieptr are start and end pointers for interpolated data from this grid into a global array chimera-boundary storage array, known as QBC. These can be computed as: iisptr(1) = 1 iisptr(nz) = iieptr(nz-1) + 1 for nz=2,nzone iieptr(nz) = iisptr(nz) - 1 - iipnts(nz) for nz=1,nzone The ibc array contains an index into the QBC array for each point in the ibpnts arrays. The interpolation and updating of chimera boundary points is done in the following manner: do i=1,iipnts(nz) s1 = q(ji(i) ,ki(i) ,li(i) ,n) s2 = q(ji(i)+1,ki(i) ,li(i) ,n) s3 = q(ji(i)+1,ki(i)+1,li(i) ,n) s4 = q(ji(i) ,ki(i)+1,li(i) ,n) s5 = q(ji(i) ,ki(i) ,li(i)+1,n) s6 = q(ji(i)+1,ki(i) ,li(i)+1,n) s7 = q(ji(i)+1,ki(i)+1,li(i)+1,n) s8 = q(ji(i) ,ki(i)+1,li(i)+1,n) c a1 = s1 a2 = -s1 +s2 a3 = -s1 +s4 a4 = -s1 +s5 a5 = s1 -s2 +s3 -s4 a6 = s1 -s2 -s5 +s6 a7 = s1 -s4 -s5 +s8 a8 = -s1 +s2 -s3 +s4 +s5 -s6 +s7 -s8 c qbc(iisptr-1+i) = a1 & + a2*dj(i) & + a3*dk(i) & + a4*dl(i) & + a5*dj(i)*dk(i) & + a6*dj(i)*dl(i) & + a7*dk(i)*dl(i) & + a8*dj(i)*dk(i)*dl(i) enddo do i=1,ibpnts(nz) q(jb(i),kb(i),lb(i)) = qbc(ibc(i)) enddo Multiblock Patched grids ------------------------ The multiblock, patched grid interface capability in INS3D requires that two neighboring grids overlap such that they share at least two coincident planes. The interzonal boundary condition information is passed into INS3D using the overset type of boundary condition file INTOUT. Additional Notes: ----------------- Steady-state convergence: The most important input parameters when trying to make the code converge are beta, dtau, and the line-relaxation sweeps parameters njswp, nkswp, and nlswp. Finding the best choice for these inputs is best done through numerical experiment. Typically, one should choose line-relaxation sweeps depending on which family of grid lines is aligned normal to the no-slip walls. For example, if the grid has a wall at an l=1 surface, it is ususally best to set njswp=nkswp=0, and set nlswp to a value between 2 and 10. Use a lower value for easier cases (ie laminar flow with larger spacings at the wall), and try a higher value for cases with finer grid spacings at the walls. For internal flows with multiple walls, it might be best to use sweeps on all directions, ie. njswp=nkswp=nlswp=2. 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. The metric subroutine which was added under version 1.80-37 can produce negative jacobians on grids which previously had only positive jacobians. This happens on meshes which have extreme stretching. The best thing to do is correct the mesh by adding more points or changing the spacing to reduce the stretching ratio. Alternatively, the older version of subroutine metric is included in the source file metric.f, and has been renamed metric2. This could be utilized in such cases.
Updated: Wednesday, 18-Oct-2000 06:53:32 PDT
WebWork: Stuart E. Rogers <rogers@nas.nasa.gov> NASA Responsible Official: Stuart Rogers |