![[Home]](/images/NASlogo.gif){kind=small}
![[NASA]](/images/NASAball.gif)
|
Updated: Wednesday, 18-Oct-2000 06:53:32 PDT
WebWork: Stuart E. Rogers <rogers@nas.nasa.gov> NASA Responsible Official: Stuart Rogers |
![[Help]](/images/help.gif)
|
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 |
|