 |
Description  |
|
|
BACKGROUND OF THE INVENTION
This invention relates to computer systems for simulating physical
processes, e.g., fluid flow.
The conventional approach to simulating high Reynolds number flow has been
to generate discretized solutions of the Navier-Stokes differential
equations, in which high-precision floating point arithmetic operations
are performed at each of many discrete spatial locations on variables
representing the macroscopic physical quantities (e.g., density,
temperature, flow velocity). The fastest and most powerful computers
available are used, and yet very limited and inexact results have been
achieved. To keep run times manageable, very coarse grid resolutions are
used, and even at those coarse resolutions there are unacceptable errors
in the solutions due to accumulated round off errors inherent in
performing successive floating point arithmetic operations.
There has long been an effort to replace the differential equation approach
with what is generally known as lattice gas (or cellular) automata, in
which the macroscopic-level simulation provided by solving the
Navier-Stokes equations is replaced by a microscopic-level model that
performs operations on particles moving between sites on a lattice. The
goal has long been to find a microscopic-level model of particle
interactions and movement that would produce the correct macroscopic
results (i.e., variations in density, temperature, etc. as prescribed by
the Navier Stokes equations).
The traditional lattice gas simulation assumes a limited number of
particles at each lattice site, with the particles being represented by a
short vector of bits. Each bit represents a particle moving in a
particular direction. For example, one bit in the vector might represent
the presence (when set to 1) or absence (when set to 0) of a particle
moving along a particular direction. Such a vector might have six bits,
with, for example, the values 110000 indicating two particles moving in
opposite directions along the X axis, and no particles moving along the Y
and Z axes. A set of collision rules governs the behavior of collisions
between particles at each site (e.g., a 110000 vector might become a
001100 vector, indicating that a collision between the two particles
moving along the X axis produced two particles moving away along the Y
axis). The rules are implemented by supplying the state vector to a lookup
table, which performs a permutation on the bits (e.g., transforming the
110000 to 001100). Particles are then moved to adjoining sites (e.g., the
two particles moving along the Y axis would be moved to neighboring sites
to the left and right along the Y axis).
Molvig et al. taught an improved lattice gas technique in which, among
other things, many more bits were added to the state vector at each
lattice site (e.g., 54 bits for subsonic flow) to provide variation in
particle energy and movement direction, and collision rules involving
subsets of the full state vector were employed. Molvig et al
PCT/US91/04930; Molvig et al., "Removing the Discreteness Artifacts in 3D
Lattice-Gas Fluids", Proceedings of the Workshop on Discrete Kinetic
Theory, Lattice Gas Dynamics, and Foundations of Hydrodynamics, World
Scientific Publishing Co., Pte., Ltd., Singapore (1989); Molvig et al.,
"Multi-species Lattice-Gas Automata for Realistic Fluid Dynamics",
Springer Proceedings in Physics, Vol. 46, Cellular Automata and Modeling
of Complex Physical Systems, Springer-Verlag Berlin, Heidelberg (1990)
(all hereby incorporated by reference). These improvements and others
taught by Molvig et al. produced the first practical lattice-gas computer
system. Discreteness artifacts that had made earlier lattice gas models
inaccurate at modeling fluid flow were eliminated.
In another approach to avoiding discreteness artifacts, referred to as the
lattice-Boltzmann model, the boolean variables of lattice gas techniques
are replaced by real variables. Chen et al., "Lattice Boltzmann Model for
Simulation of Magnetohydrodynamics," Physical Review letters, Vol. 67, n.
27, 30 Dec. 1991; Qian et al., "Lattice BGK Models for Navier-Stokes
Equation," Europhysics Letters, Vol. 17, pp. 479-484, 1 Feb. 1992. Rather
than monitoring the presence of individual particles for each state of
each site in the lattice, the lattice-Boltzmann model monitors the
particle distribution function for each said state.
As described by Qian et al., the real numbers used in the lattice-Boltzmann
approach permit application of the method of relaxation, which can be
described simply as:
N(t+1)=(1-.omega.)N(t)+.omega.N.sub.e,
where N(t) is a quantity at time t, N.sub.e is the quantity's equilibrium
(Boltzmann) value, and .omega. is a relaxation parameter having a value
between 0 and 2. The method is referred to as subrelaxation when
0<.omega.<1, and over-relaxation when 1<.omega.<2. (When .omega.=1, N(t+1)
simply equals N.sub.e, and no relaxation occurs.) Qian et al. also noted
that the shear viscosity of a simulated fluid can be reduced by increasing
.omega..
SUMMARY OF THE INVENTION
The invention features modifying the viscosity of a lattice gas system in
which each lattice site (voxel) is represented by a state vector that
includes an integer value for each state (e.g., in each voxel, 0-255
elements can be moving in a particular direction with a particular
energy). Decreasing the viscosity of the system increases the effective
density of the lattice and thereby dramatically increases the resolution
at which the system is able to simulate a given physical process.
Increasing the viscosity of the system permits simulation of high
viscosity fluids.
Viscosity is a measure of a fluid's resistance to a shear force (i.e., a
force which acts parallel to the direction of fluid flow). In an actual
fluid, viscosity results from interactions between neighboring particles
in the fluid that cause the velocities of the particles to gravitate
toward an average value. In a lattice system, viscosity results from
interactions between particles positioned in specific voxels that cause
the net velocity of the particles positioned in a voxel to gravitate
toward the net velocity of the particles positioned in neighboring voxels.
Because each voxel in a lattice system represents a region of simulated
space that is substantially larger than the physical space that would be
occupied by an actual particle, the viscosity resulting from interactions
between voxels is substantially greater than that resulting from molecular
particle interactions in real fluids (i.e., the "averaging" resulting from
each voxel interaction affects a substantially larger region of space than
that resulting from each molecular particle interaction).
Typically, the viscosity of the system is modified using relaxation
techniques. Because the relaxation technique described by Qian et al. is
directed to a system that uses real numbers, its applicability to an
integer-based system is not readily apparent. Indeed, direct application
of the relaxation technique described by Qian et al. to a lattice gas
system using state vectors comprised of integer values would not work
because multiplication of an integer value by the real-valued relaxation
parameter would result in a real value rather than an integer value. While
the real value produced by the multiplication could be rounded to an
integer value, such rounding would result in the system no longer
conserving mass, momentum and energy. The invention permits the
integer-valued state vectors to be modified by real-valued relaxation
parameters, and does so while conserving mass, momentum and energy, and
without introducing roundoff error.
Viscosity in a lattice system can be reduced by increasing the density of
the lattice (i.e., by decreasing the quantity of simulated space that is
represented by each voxel), and can also be reduced through use of
over-relaxation. Viscosity, .nu., can be expressed in terms of .omega.,
the relaxation parameter:
##EQU1##
where T is the temperature of the fluid. Thus, for example, relative to a
relaxation parameter of one (.nu.=T/2), a relaxation parameter of 1.8
(.nu.=T/18) will reduce the viscosity in the lattice by a factor of nine.
The invention promises to substantially improve the ability of lattice
systems to simulate physical processes. Because over-relaxation has the
same effect as increasing the density of the lattice (i.e., reducing the
viscosity of the lattice system), use of over-relaxation effectively
increases the density of the lattice. Use of over-relaxation therefore has
a dramatic effect on the processing necessary to simulate a physical
system with a particular resolution (or the resolution with which a
particular processor can simulate a physical system). For example, a
tenfold increase in the effective density of a three dimensional lattice
reduces the processing required to simulate a physical system with the
lattice to a particular level of resolution by a factor of almost ten
thousand (i.e., ten cubed less the additional processing required to
implement over-relaxation and multiplied by a tenfold decrease in the time
required to simulate a fluid of a given velocity).
The invention is implemented in a new computer system for simulating
physical processes. Instead of the lattice gas model in which at each
lattice site (voxel) there is at most a single element in any momentum
state (e.g., at most a single element moving in a particular direction
with a particular energy), the invention uses a multi-element technique in
which, at each voxel, multiple elements can exist at each of a number of
states (e.g., 0-255 elements can be moving in a particular direction). The
state vector, instead of being a set of bits, is now a set of integers
(e.g., a set of eight bit bytes providing integers in the range of 0 to
255), each of which represents the number of elements in a given state.
To model interactions between elements of different momentum states within
a voxel, the computer system performs interaction operations on the state
vectors. Typically, these interaction operations are performed for each
voxel by applying a set of rules to the state vector of the voxel, where
each rule modifies a particular set of integers from the state vector.
To change the viscosity of the simulated physical process, the computer
system performs viscosity modification operations on the state vectors.
These operations are typically performed after the interaction operations
and apply a set of rules that are similar to, or the same as, the rules
applied during the interaction operations. Where the same rules are
applied, the rules modify the state vectors by a first amount during the
interaction operations and a second amount during the viscosity
modification operations, where the first amount is related to the second
amount by a relaxation parameter. Because the rules used in the
interaction operations conserve mass, momentum and energy, this approach
ensures that these properties will be conserved during the viscosity
modification operations.
When the second amount is determined by multiplying the first amount by a
real number derived from the relaxation parameter, the result of the
multiplication is truncated to ensure that the second amount is an integer
value. To prevent the truncation operation from introducing statistical
bias into the system, a random value between zero and one is added to the
result of the multiplication prior to truncation.
The viscosity of the lattice system is reduced by using a relaxation
parameter having a value greater than one and less than two. As the
relaxation parameter approaches two, the viscosity of the simulated system
approaches zero and the system becomes unstable. Viscosity, which is
essentially a form of friction, tends to damp out fluctuations in the
system. Thus, instability occurs when there is no viscosity because these
fluctuations are allowed to spread unchecked through the system. It has
been found that instability can generally be avoided by using a relaxation
parameter that is less than or equal to 1.8.
After performing the interaction and viscosity modification operations, the
computer system performs move operations on the state vectors that reflect
movement of elements to new voxels.
BRIEF DESCRIPTION OF THE DRAWINGS
FIG. 1 is a flow chart of a procedure followed by a physical process
simulation system.
FIG. 2 is a perspective view of a microblock.
FIG. 3 is a flow chart of a procedure for performing collisions with
over-relaxation.
FIG. 4 is a block diagram of a system for performing slip surface dynamics.
FIG. 5 is an illustration of specular reflection.
FIG. 6 is a block diagram of a functional unit of a physical process
simulation system.
FIG. 7 is a block diagram of a microdynamics unit of the system of FIG. 6.
FIG. 8 is a block diagram of a single-voxel data path of the microdynamics
unit of FIG. 7.
DESCRIPTION OF THE PREFERRED EMBODIMENTS
Referring to FIG. 1, a physical process simulation system operates
according to a procedure 100. At startup, a timer is initialized (step
102). Next, a voxel counter, which points to a particular voxel (or
location) within the lattice, is initialized to point to the first voxel
in the lattice (step 104).
After initialization, the system loads the state vector corresponding to
the voxel designated by the voxel count (step 106). The state vector
completely defines the status of the voxel, and includes 49 or more
multi-bit entries, each of which corresponds to an integer value. These 49
entries correspond to a rest state, 24 directional vectors at a first
energy level, and 24 directional vectors at a second energy level. Though
only 49 entries are required, preferred embodiments provide for six rest
states and therefore use 54 entries. Six rest states are employed to
ensure that there are a sufficient number of rest "slots". Of course, this
same effect could be achieved by increasing the number of bits in the
single entry corresponding to the rest state in a 49 entry embodiment. By
using multi-bit entries, the system offers substantially improved
performance over systems that used single-bit entries to define voxel
status. In particular, unlike single-bit systems that could only produce
Fermi-Dirac statistics, which are unsuitable for many applications, the
system can produce Maxwell-Boltzmann statistics.
After loading the state vector, the system performs all intravoxel
operations on the state vector (step 108). Intravoxel operations are those
that do not require information about other voxels. For example, in a
fluid simulation system, intravoxel operations would account for
collisions between particles within a voxel.
Upon completing the intravoxel operations, the system increments the voxel
counter (step 110). If the new voxel count does not exceed the number of
voxels in the lattice (step 112), the system loads the state vector of the
next voxel (step 106) and continues processing.
If the new voxel count exceeds the number of voxels in the lattice (step
112), the system performs intervoxel operations (step 114). Intervoxel
operations are those that require information from more than one voxel.
For example, in a fluid simulation system, intervoxel operations would
account for movement of particles between voxels. After performing
intervoxel operations, the system increments the time (step 116),
reinitializes the voxel counter to point to the first voxel (step 104),
and continues processing.
Operation of a preferred embodiment of the system is described in detail
below. For clarity, the system described above has been described as
operating serially. However, as noted below, the system, like other
lattice systems, is ideally suited for parallel operations. For example,
intravoxel operations could be performed on multiple voxels
simultaneously. Similarly, as long as intravoxel operations on all of the
voxels involved in an intervoxel operation are complete, the intervoxel
operation could be performed simultaneously with other intravoxel
operations.
The disclosures of U.S. application Ser. No. 08/030,573, filed Mar. 12,
1993; PCT application Ser. No. PCT/US91/04930, filed Jul. 12, 1991; U.S.
application Ser. No. 07/812,881, filed Dec. 20, 1991; U.S. application
Ser. No. 07/555,754, filed Jul. 12, 1990; and U.S. application Ser. No.
08/165,293, filed Dec. 10, 1993 are all hereby incorporated by reference.
Before any of the computational operations are described, it is necessary
to briefly describe the elementary data structure that comprises the basic
state vector for each voxel. This is the basic element upon which the
majority of required computations operate. Each lattice site, or voxel
(these two terms are used interchangeably throughout this document),
contains 54 states for subsonic mono-species simulations. The number of
states will be lengthened for transonic flows or multiple-species
simulations.
In this document the state space is represented with the following
notation:
N.sub.i (x, t)
N.sub.i represents the number of particles in state i, at the lattice site
denoted by the 3-dimensional vector x, at time-step t.
The number of states is determined by the number of possible velocity
vectors within each energy level. The velocity vectors consist of integer
linear speeds in 4-dimensional space: x, y, z, and w. The 4'th dimension,
w, is projected back onto 3-dimensional space and thus does not indicate
an actual velocity in the 3-dimensional lattice. For subsonic flows, i
ranges from 0 to 53.
Each state represents a different velocity vector at a specific energy
level. The velocity of each state is indicated with its "speed" in each of
the 4 dimensions as follows:
c.sub.i =(c.sub.x, c.sub.y, c.sub.z, c.sub.w)
The energy level 0 state is known as a stopped particle, they are not
moving in any dimension, i.e. c.sub.stopped =(0, 0, 0, 0). Energy level 1
states have a .+-.1 in two of the four dimensions and a 0 velocity in the
other two. And energy level 2 states have either a .+-.1 in all four
dimensions, or a .+-.2 in one of the 4 dimensions and a 0 velocity in the
other three. Generating all of the possible permutations of these 3 energy
levels gives a total of 49 possible states (1 energy 0 state, 24 energy 1
states, 24 energy 2 states). In addition, the subsonic flow state space
maintains a total of 6 stopped states as opposed to 1, giving a total
state count of 54 instead of 49.
The voxels are grouped in to small 2.times.2.times.2 volumes that are
called microblocks. The microblocks are organized to optimize parallel
processing of the lattice sites as well as to minimize the overhead
associated with the data structure. A short-hand notation for the 8
lattice sites in the microblock is defined below and is used throughout
this document.
N.sub.i (x)
Where x.epsilon.{0, 1, 2, . . . , 7}
Where x represents the relative position of the lattice site within the
microblock. A microblock is illustrated in FIG. 2.
Microdynamics (Intravoxel Operations)
The microdynamics operations are those set of physical interactions that
occur purely within a voxel. This class of operations allows for the
permutation of the fluid state space to account for the physical
interactions between fluid particles and various types of object surfaces.
Normal Collisions
Normal collisions are operations that allow for particles to collide with
each other, thus changing their velocity and direction. A change in a
particle's velocity and direction is accomplished by moving that particle
into a different state, since it is the state that a particle is in that
determines its velocity vector.
The typical collision operation consists of two pairs of input state
vectors (4 in total) and likewise two pairs of output state vectors. The
basic collision operation collides two "incoming" particles and changes
their state into two "outgoing" particles. The incoming and outgoing pairs
must always conserve mass, momentum, and energy. Therefore, not all
possible quartets within the 54 states are "legal" collision sets.
The basic collision operator is bi-directional in nature, thus the
"incoming" and "outgoing" states are determined at the time the collision
takes place depending on the local collision states' populations. Two
pairs are selected and depending on local density, one pair (incoming)
will be the source of particles into the other (outgoing) pair.
The basic collision operation is described below.
C=SignOf[(N.sub.i .multidot.N.sub.j)-(N.sub.k .multidot.N.sub.l)]
NScatt=C.multidot..delta.
Where .delta..epsilon.{1,2,4,8}
N.sub.i =N.sub.i -NScatt
N.sub.j =N.sub.j -NScatt
N.sub.k =N.sub.k +NScatt
N.sub.l =N.sub.l +NScatt
Where SignOf is a function that returns only the sign (.+-.1) of the
operations enclosed in brackets. The SignOf operator returns a +1 if the
value is a 0. The number of particles scattered from states i and j to
states k and l, NScatt, is determined by multiplying the sign of the
collision operation, C, by a small positive constant .delta. (delta). The
.delta. is specified along with the state indices in the collision rule
list. The pair of states, i and j or k and l, with the larger product
becomes the source of the particles to the pair with the lower product. If
NScatt is negative, then particles are transferred from states k and l
into state i and j.
All 4 of the states represent particles at the same voxel. All of the
collisions only depend on the state information local to that particular
site. The state indices i, j, k, and l are determined such that a particle
from each of states i and j have the same total momentum and energy as a
particle from each of states k and l. All 4 of the indices must represent
4 different states and all 4 of the states must be at the same energy
level.
As an example of a normal collision, the following initial state is
proposed.
##EQU2##
As can be seen from the state selection for i, j, k, and l, the pair of i
and j have a net momentum of +2 in the x dimension and 0 in y, z, and w.
In addition, both particles are at energy level 1. The same is true of the
pair k and l. After the first step in the collision process described
above, C is equal to -1 (SignOf[(25.multidot.40)-(30.multidot.50)]=-1). A
negative sign indicates that the k and l pair are sourcing particles to
the i and j pair.
NScatt is calculated to be -4, based on a delta of 4 that was specified in
the collision rule. NScatt is then substracted from states i and j, and
added to k and l. The collision operation is now complete and four new
output state populations have been created. Where now:
N.sub.i =29
N.sub.j =44
N.sub.k =26
N.sub.l =46
By moving the same number of particles out of states k and l and into
states i and j, mass is also conserved.
There does exist a potential for overflow and underflow of a state's
particle count in the collision operation described above. It should be
noted that conservation of mass, momentum, and energy is paramount in this
simulation environment and that an overflow of a state would result in a
loss of mass, as well as momentum and energy, if it went unchecked.
Likewise, it is also possible to have underflow occur, in which case mass
would be created, not destroyed. Thus, it is required that the collision
operations preserve the mass of the quartet of states that it operates on.
This is accomplished by preventing any exchange of particles if the
operation would cause either an overflow or underflow in any of the states
involved in the collision.
Energy Exchanging Collisions
Energy exchanging collisions are performed just like the normal collisions
described in the previous section, except that the two outgoing particles
are at different energy levels than the two incoming particles. For
subsonic flows there are only 3 energy levels: 0 (stopped), 1, and 2. In
order to conserve energy, the only possible energy exchanging collisions
occur with one pair consisting of two energy 1 particles while the other
pair is comprised of an energy 2 particle and a stopped particle. What
makes these collisions special is that these collisions do not happen at
the full rate of normal collisions. Therefore, the class of energy
exchanging collisions must only be allowed to occur at a specified rate.
The rate is specified in two components, a forward and inverse rate. These
rates can be specified as integers and implemented in the collision
process as follows.
##EQU3##
Where are R.sub.1.fwdarw.2 and R.sub.2.fwdarw.1 are rates that control the
exchange of particles between states of different energy levels. The
energy exchange collisions are structured such that two energy 1 particles
(state i and j) are collided to produce an energy 2 particle and an energy
0 particle (states k and l), or vice versa. As with normal collisions,
mass, momentum, and energy are always conserved. The determination of
these rates is based on temperature.
The two rates, R.sub.1.fwdarw.2 and R.sub.2.fwdarw.1, that are used in the
energy exchange collisions are calculated based on the temperature of the
fluid. The temperature of the fluid, however, is not necessarily constant
over the length of a simulation, especially for simulations involving heat
transfer. These rates will have to be updated dynamically during the
simulation to reflect changes in the local temperature.
The two rates above are not independently calculated. The only relevant
information provided by the two rates is in their ratio. The overall Rate
is the ratio:
##EQU4##
Where the Rate is calculated based from the temperature, T, as described
below.
##EQU5##
Expanding the product terms out we get the following equation:
##EQU6##
The temperature range supported for subsonic flows is between 1/3 and 2/3.
If the temperature, T, is less than 0.5 the Rate is less than 1 and is
greater than 1 for temperatures above 0.5. From this Rate, the two energy
rates can then be determined. However, the energy exchange rates have to
be scaled within the allowable precision.
The same concerns regarding overflow and underflow conditions also apply to
these collision operations. These conditions are possible with energy
exchanging collisions and must also be prevented here as well.
Over-relaxation (Viscosity Reduction)
Referring to FIG. 3, the physical process simulation system performs
collisions with viscosity reduction for each fluid voxel according to a
procedure 120. Procedure 120 includes a collision stage (steps 122-134),
in which the system performs normal collision operations for the voxel to
place the voxel in the equilibrium (Boltzmann) state, and a viscosity
reduction stage (steps 136-148), in which the system performs the
viscosity reduction.
At the beginning of the collision stage, the system initializes a rule
counter to a value of one (step 122). Next, the system determines which
states of the state vector are affected by the rule corresponding to the
rule counter, and, assuming that the rule affects four states, loads the
identity of those states into variables i, j, k, and l (step 124). (Though
the following discussion assumes that each rule affects exactly four
states; the system accommodates rules affecting different number of
states.) After loading i, j, k, and l, the system determines N.sub.Scatt,
the number of particles scattered, as a function of N.sub.i, N.sub.j,
N.sub.k, and N.sub.l using either the normal collision function or the
energy exchanging collision function discussed above (step 126). The
system then subtracts N.sub.Scatt from N.sub.i and N.sub.j and adds
N.sub.Scatt (or 0 if the original change was omitted to avoid underflow or
overflow) to N.sub.k and N.sub.l (step 128). As discussed above, the
system does not perform the exchange of particles if the exchange would
result in an underflow or overflow condition in any of the affected
states. Next, the system stores N.sub.Scatt in an array entry (N.sub.Scatt
[r]) corresponding to the value of the rule counter (r) (step 130), and
increments the rule counter (step 132). If the incremented rule counter
does not exceed the maximum number of rules (step 134), which is 276 in
the preferred embodiment, the system loads values into i, j, k, and l
based on the rule corresponding to the incremented rule counter (step 124)
and continues with the collision stage as described above.
If the incremented rule counter exceeds the maximum number of rules (step
134), the system enters the viscosity reduction stage. First, the system
reinitializes the rule counter to a value of one (step 136). Next, the
system determines N.sub.VR, the number of particles scattered due to
viscosity reduction, as:
N.sub.VR =.left brkt-bot.(.omega.-1)N.sub.Scatt [r]+noise.right brkt-bot.
where .omega. is the relaxation parameter, N.sub.Scatt [r] is the array
entry corresponding to the rule counter, and "noise" is a random value
between zero and one (step 138). Use of "noise" ensures that the
truncation operation, which forces N.sub.VR to be an integer value, will
not statistically bias N.sub.VR in a particular direction. After
determining N.sub.VR, the system loads the identity of the states affected
by the rule corresponding to the rule counter into variables i, j, k, and
l (step 140). The system then subtracts N.sub.VR from N.sub.i and N.sub.j
and adds N.sub.Scatt to N.sub.k and N.sub.l (step 142). A | | |
