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