Computer system for simulating physical processes using surface algorithms - Patents.com

Abstract

To provide computer implemented techniques for simulating elements of a fluid flow.SOLUTION: These techniques include the steps of: storing in a memory state vectors for a plurality of voxels, the state vectors comprising a plurality of entries that correspond to particular momentum states of a plurality of possible momentum states at a voxel; storing in a memory a representation of at least one surface that is sized and oriented independently of the size and orientation of the voxels; performing interaction operations on the state vectors, the interaction operations modelling interactions between elements of different momentum states; performing surface interaction operations on the representation of the surface, the surface interaction operations modelling interactions between the surface and substantially all elements of voxels; and performing move operations on the state vectors to reflect movement of elements to new voxels.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] The present description relates to simulating physical processes, for example fluid flow. [Background technology]

[0002] The Lattice Boltzmann Method (LBM) is used to simulate diverse and complex fluid flows around various geometric shapes such as cars and airplanes. One of the most important parts of the entire solver is its algorithms that handle the boundary conditions, i.e. the dynamics of the LBM particle distribution at the edge of the fluid computation domain, such as near the surface of a solid wall. How well such algorithms are formulated has a direct impact on the resulting accuracy of the simulation.

[0003] U.S. Patent No. 5,848,260 describes a novel technique for a volume-based formulation (the '260 slip algorithm) that ensures strict mass conservation and accurate application of momentum fluxes across interfaces. [Prior art documents] [Patent documents]

[0004] [Patent Document 1] U.S. Patent No. 5,848,260 Summary of the Invention

[0005] Although the patented approach has been very successful, it does have some limitations. For example, the resulting fluid quantities near the walls result in higher than desired artificial noise. Such artifacts manifest themselves as effective surface roughness, resulting in increased numerical dissipation and thicker than desired boundary layers. This can negatively impact the accuracy of the next level of the simulation, especially for very smooth streamlined enclosures. Secondly, the plateau of the patented process is not wide enough, thus limiting its capabilities for simulating faster fluid flows.

[0006] The process described in the '260 patent, i.e., the slip algorithm, is based on a volumetric representation of a lattice Boltzmann fluid domain along with a surface representation. In the '260 patent, the surface is described as a collection of surface elements, each having a specific surface area and a surface normal. In the '260 patent, boundary conditions are achieved by the interaction of adjacent surface elements with particles in the fluid domain.

[0007] However, in the novel slip algorithm described here, boundary conditions are provided by the interactions of all surface elements with particles in the fluid domain, as defined by three fundamental processes:

[0008] The incoming distribution functions (corresponding to particles moving towards the interface) from the fluid region near the interface are collected according to a so-called parallelogram / parallel hexahedron. Surface dynamics involves transforming the set of incoming distribution functions into a set of outgoing distribution functions (corresponding to particles reflected from the interface) that satisfy certain boundary conditions on the mass and momentum fluxes passing through the interface. The outgoing distribution functions are distributed in the fluid region near the interface according to a parallelogram / parallel hexahedron.

[0009] This approach largely addresses the limitations mentioned above. While the overall framework of the new approach is similar to that of the '260 patent, there is a significant change in the surface dynamics that transforms a set of inflow distribution functions into a set of outflow distributions. More specifically, the process used in the surface dynamics to achieve zero momentum flux through the surface is significantly modified. Achieving an exact zero tangential momentum flux is the basis of an accurate boundary condition algorithm so that specific amounts of momentum flux can be added later with precision. This modification helps overcome the long-standing limitations mentioned above.

[0010] According to one aspect, a computer-implemented method for simulating elements of a fluid flow includes storing in memory state vectors of a plurality of voxels, the state vector including a plurality of entries corresponding to particular momentum states of a plurality of possible momentum states at the voxel; storing in memory a representation of at least one surface sized and oriented independent of the size and orientation of the voxel; performing interaction operations on the state vectors, the interaction operations modeling interactions between elements of different momentum states; performing surface interaction operations on the representation of the surface, the surface interaction operations modeling interactions between the surface and substantially all elements of the voxel; and performing move operations on the state vectors to reflect movement of elements to new voxels.

[0011] According to a further aspect, a data processing system for simulating elements of a fluid flow includes instructions for causing the data processing system to store in a memory state vectors of a plurality of voxels, the state vector including a plurality of entries corresponding to particular momentum states of a plurality of possible momentum states for the voxel; store in the memory a representation of at least one surface sized and oriented independent of the size and orientation of the voxel; perform an interaction operation on the state vector, the interaction operation models interaction between elements of different momentum states; perform a surface interaction operation on the representation of the surface, the surface interaction operation models interaction between the surface and substantially all elements of the voxel; and perform a movement operation on the state vector to reflect movement of the elements to new voxels.

[0012] According to a further aspect, a non-transitory computer readable medium storing a computer program product for simulating elements of a fluid flow, the computer program product including instructions for causing a data processing system to: store in a memory state vectors of a plurality of voxels, the state vector including a plurality of entries corresponding to particular momentum states of a plurality of possible momentum states at the voxel; store in the memory a representation of at least one surface sized and oriented independent of the size and orientation of the voxel; perform an interaction operation on the state vector, the interaction operation models interaction between elements of different momentum states; perform a surface interaction operation on the representation of the surface, the surface interaction operation models interaction between the surface and substantially all elements of the voxel; and perform a movement operation on the state vector to reflect movement of the elements to new voxels.

[0013] One or more of the above aspects may include one or more of the following features, among others described herein.

[0014] The instructions for performing the surface interaction operation include instructions for collecting elements from a first set of at least one voxel that interacts with the facet, modeling interactions between the collected elements and all of the facets to generate a set of surface-interacting elements, and scattering the surface-interacting elements to a second set of at least one voxel that interacts with the facet. An entry of the state vector represents a density of elements per unit volume in a particular momentum state of the voxel. The state vector includes one or more of a plurality of integer and / or floating point values ​​that represent a density of elements per unit volume in a particular momentum state and have a greater number of possible values.

[0015] The surface interaction operation is performed using one or more of integer and / or floating point values, and the surface interaction operation is performed using values ​​representing real numbers. The surface interaction operation is performed using integer and / or floating point numbers. The surface interaction operation includes representing voxels that intersect with the surface as sub-voxels. The representation of the at least one surface includes a plurality of facets that are sized and oriented independent of the size and orientation of the voxels and represent the at least one surface, and the surface interaction operation further includes modeling an interaction between the facets and elements of at least one voxel near the facets.

[0016] An element represents a particle of the fluid, and a facet represents at least one surface through which the fluid flows.

[0017] Performing the surface interaction operation includes calculating the difference between the inflow distribution and the Boltzmann distribution, determining the total momentum from the differences of all state vectors, and generating the outflow distribution based on the determined differences.

[0018] One or more of the above aspects may provide one or more of the advantages disclosed herein.

[0019] In addition to handling boundary conditions on solid walls, the slip algorithm can also be applied to handle fluid boundaries between domains of two different lattice Boltzmann solvers and between domains in two different reference frames. Unlike solid walls, the mass flux across a boundary is generally non-zero. Also, the momentum flux includes an additional term that represents the convection of flow from one domain to another. All of this is a non-trivial but simple extension of the previous slip algorithm for solid surfaces ('260).

[0020] Other features and advantages of the present invention will become apparent from the following detailed description of the preferred embodiments, and from the claims. [Brief description of the drawings]

[0021] [Figure 1] FIG. 1 illustrates a system for simulating fluid flow including a novel surface dynamics transformation. [Diagram 2] 1 is a flow diagram illustrating operations for formulating a lattice Boltzmann model simulation with surface dynamics transformation. [Diagram 3] 1 is a flow diagram showing a simulation operation using a lattice Boltzmann model with surface dynamics transformation. [Figure 4] 1 is a flow diagram showing dynamics transformation. [Diagram 5] FIG. 1 illustrates velocity components of two LBM models represented in Euclidean space (prior art). [Figure 6] FIG. 1 illustrates velocity components of two LBM models represented in Euclidean space (prior art). [Figure 7] 4 is a flow diagram of the procedure followed by a physical process simulation system using a corrected CAD drawing. [Figure 8] FIG. 1 is a perspective view of a microblock (prior art). [Figure 9A]FIG. 2 is a diagram of a lattice structure used by the system of FIG. 1 (prior art). [Figure 9B] FIG. 2 is a diagram of a lattice structure used by the system of FIG. 1 (prior art). [Figure 10] FIG. 1 illustrates a variable resolution technique (prior art). [Figure 11] FIG. 1 illustrates a variable resolution technique (prior art). [Figure 12] FIG. 1 illustrates particle movement (prior art). [Figure 13] FIG. 2 shows the faceted area of ​​a surface (prior art). [Figure 14] 1 is a flow diagram of the novel surface dynamics. [Figure 15] 1 is a flow diagram of a procedure for performing surface dynamics (prior art). DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0022] The details of one or more embodiments of the invention are set forth in the accompanying drawings and the description below. Other features, objects, and advantages of the invention will become apparent from the description and drawings, and from the claims.

[0023] One of the methods for simulating fluid flows is the so-called Lattice Boltzmann Model (LBM). In an LBM-based physical process simulation system, the fluid flow is represented by a distribution function value and evaluated at a set of discrete velocities using the well-known lattice Boltzmann equation (see Equation 1 below), which describes the time evolution of the distribution function. The distribution function includes two processes: the flow process and the collision process.

[0024] In surface dynamics, a set of "inflow" particle distributions is mapped / transformed into a set of "outflow" particle distributions such that the exact zero tangential momentum flux boundary condition is achieved. As noted above, this specific part is where the method differs significantly from the '260 patent. Before describing the specific algorithmic differences in the following steps of surface dynamics, we first provide an overview of the concept of surface dynamics below.

[0025] For each surface element,

number

number

number

number

[0026] In surface dynamics, the value of the "outflow" distribution function of a surface element is determined based on the known value of the "inflow" distribution function.

number

number

[0027] Symbolically, surface dynamics can be expressed as the following logical transformation:

number

number

number

[0028] Fluid boundary conditions are typically defined such that mass M and momentum P obey specified values ​​on the interface of the fluid solver. Specifically, since particles cannot penetrate solid walls, a zero mass flux condition, M=0, applies. Zero tangential momentum flux through a surface is given by

number

[0029] For a solid wall, this is equivalent to having a frictionless boundary condition, where

number

number

[0030] The surface dynamics in the slip algorithm is a relationship that converts the values ​​of the inflow distribution function to values ​​of the outflow distribution such that the two flux conditions with specified values ​​are satisfied. After achieving the zero mass flux condition and the zero tangential momentum flux condition, the remaining process to satisfy the overall mass flux and momentum flux conditions with specified values ​​is the same between the '260 patent and the new approach, so the process will not be described here.

[0031] Since there are many distribution functions other than the two flux conditions, the value of the outflow distribution is not unique. In fact, the '260 patent and the current method give different values ​​because of different relationships linking the inflow distribution function with the outflow distribution function.

[0032] A surface dynamics algorithm for zero tangential momentum flux. The specific relationship in the '260 patent that satisfies zero mass flux and zero tangential momentum flux is the following formulation:

number

number

number

[0033] For example, for an isothermal LBM solver, the equilibrium distribution is given by Equation 5,

number

number

number

[0034] Substituting this equation (equation 4) along with (equation 5) into the definition of momentum flux (equation 3), we get

number

number

[0035] Here, the quantity p=ρT0 denotes the pressure value at the surface. More general boundary conditions can be realized by adding appropriate mass and momentum fluxes to the zero mass flux and zero tangential momentum flux conditions described here.

[0036] One defining feature of the '260 patent formulation is the runoff distribution, as shown in equation (4):

number

number

[0037] 1, a system 10 for performing Lattice Boltzmann (LB) based simulations includes a novel surface dynamics transformation 34b. The system 10 in this implementation is based on a client-server or cloud-based architecture and includes a server system 12 implemented as a massively parallel computing system 12 (standalone or cloud-based) and a client system 14 coupled via a network 15. The server system 12 includes a memory 18, a bus system 11, an interface 20 (e.g., a user interface / network interface / display or monitor interface, etc.), and a processing device 24. The memory 18 includes a mesh preparation engine 32 and a simulation engine 34.

[0038] 1 shows a mesh preparation engine 32 in memory 18, which may be a third-party application executing on a different system than server 12. Regardless of whether the mesh preparation engine 32 executes in memory 18 or on a different system than server 12, the mesh preparation engine 32 receives a user-provided mesh definition 30, the mesh preparation engine 32 prepares the mesh, and sends (and / or stores) the prepared mesh to the simulation engine 34.

[0039] The simulation engine 34 includes a collision interaction module 34a, which includes a novel surface dynamics transformation 34b, a boundary module 34c, and an advection particle collision interaction module 34d. The system 10 has access to a data repository 38 that stores 2D and / or 3D meshes (Cartesian and / or curvilinear), coordinate systems, and libraries.

[0040] 2, a process 40 is shown for simulating fluid flow for a representation of a physical object. In the example described herein, the physical object is an airfoil. However, the use of an airfoil is merely exemplary, as the physical object can be any shape, and in particular can have planar and / or curved surfaces. The process 40 receives 42 a mesh (or grid) for the physical object being simulated, for example, from the client system 14 or retrieves it from the data repository 38. In other embodiments, either an external system or the server 12 generates a mesh for the physical object being simulated based on user input. The process pre-computes 44 geometric quantities from the retrieved mesh and performs 46 a kinetic lattice Boltzmann model simulation using the pre-computed geometric quantities corresponding to the retrieved mesh. The lattice Boltzmann model simulation includes simulating 46 the evolution of particle distribution, including surface dynamics transformation, boundary modeling, and advection of particles to the next cell in the LBM mesh.

[0041] 3, a simulation process 46 simulates the evolution of particle distributions according to, for example, a Lattice Boltzmann Process (LBp) modified to fit the novel surface dynamics transformation 34b (FIG. 1). The process 46 (see FIG. 2) performs collision operations 46a and dynamics transformation 46b of the results of all collision operations for all particles, followed by boundary modeling 46c and advection 46d of the particles to the next cell in the LBM space.

[0042] 4, a simulation process 46 for simulating the time evolution of particles generally performs the conventional processes of collision calculation 46a, boundary modeling 46c, and advection 46d. However, the transformation of a set of inflow particle distributions to a set of outflow particle distributions described in the '260 patent is replaced by a kinetic transformation process that takes the results of all inflow particle distributions 46b to formulate each set of outflow particle distributions. The kinetic transformation 46b involves collecting 52 all inflow sets of distributions from a mesh location from the collision calculation, rather than collecting sets of inflow distributions from adjacent mesh locations. This new formulation avoids the one-to-one dependency in the '260 patent discussed above, and replaces that dependency with each "outflow" distribution function being dependent on all "inflow" distributions 54.

[0043] Instead of Equation 4, a new relationship is given in the novel approach, the following mathematical formulation:

number

number

number

[0044] In (Equation 9) the equilibrium distribution is defined in the same way as in (Equation 5).

number

number

number

number

number

[0045] Similar to the approach of the '260 patent, the novel approach defined by (Eqs. 7)-(Eqs. 9) can be shown to achieve boundary conditions of zero mass flux and zero tangential momentum flux.

[0046] Although seemingly a more complicated procedure, the main difference in the new formulation is that it avoids a one-to-one relationship between the outflow and inflow distributions. Indeed, as shown in (Eqs. 7)–(Eqs. 9), the outflow distribution

number

[0047] The new formulation largely overcomes the aforementioned shortcomings, providing less noisy surfaces and potentially supporting higher fluid velocities near the boundaries.

[0048] Detailed Example The procedure described below in Figure 7 describes a flow simulation process that uses a CAD drawing with identified free space to construct the simulation space. In Figures 5 and 6 above, as well as Figures 8, 9A-9B, 10 and 11, and Figures 12, 13, and 15, each of these figures is labeled as Prior Art because they appear in the above referenced patents.

[0049] However, the above referenced patent does not describe the dynamics conversion process 46b described herein, and therefore the diagrams that appear in the above referenced patent do not take into account modifications that may be made to the flow simulation using the dynamics conversion process 46b.

[0050] Modeling the Simulation Space In the LBM-based physical process simulation system, the fluid flow is modeled with discrete velocities c i The distribution function value f evaluated on the set i The dynamics of the distribution function is given by Equation I.1,

number

number

number

[0051] Equation (I.1) expresses the distribution function f iis the well-known lattice Boltzmann equation that describes the time evolution of . The left hand side describes the change in distribution due to the so-called "flow process". A flow process is when a pocket of fluid starts at one mesh location and then moves along one of several velocity vectors to the next mesh location. At that point we calculate the "collision factor", i.e. the influence of nearby pockets of fluid on the starting pocket of fluid. Fluids can only move to different mesh locations, so the velocity vectors need to be chosen appropriately so that all components of all velocities are multiples of a common velocity.

[0052] The right hand side of the first equation is the aforementioned "collision operator" which describes the change in the distribution function due to collisions between pockets of fluid. A particular form of the collision operator is that of the Bhatnagar, Gross, and Krook (BGK) operator. The collision operator forces the distribution function to a prescribed value given by the second equation, which is its "equilibrium" form:

[0053] The BGK operator ensures that the distribution function is consistent across the collisions, regardless of the collision details.

number

[0054] From this simulation, traditional fluid variables such as mass ρ and fluid velocity u are obtained as simple sums in the following equation (I.3):

number

[0055] By considering symmetry, the set of velocity values ​​is chosen to form a certain lattice structure when spanning the configuration space. The dynamics of such a discrete system can be expressed as a set of discrete lattices of the form f i (x+c i ,t+1)-f i (x,t)=C i (x,t) where the collision operator usually takes the BGK form as shown above. By appropriately choosing the equilibrium distribution form, it can be shown theoretically that the lattice Boltzmann equation produces the correct hydrodynamics and thermohydraulics. That is, i The hydrodynamic moments derived from (x,t) obey the Navier-Stokes equations in the macroscopic limit. These moments are defined by equation (I.3) above.

[0056] c i And w i The set of values ​​of defines the LBM model, which can be efficiently implemented on scalable computational platforms and run with good robustness for time-unsteady flows and complex boundary conditions.

[0057] The standard technique for obtaining the macroscopic equations of motion for a fluid system from the Boltzmann equation is the Chapman-Enskog method, which takes successive approximations of the full Boltzmann equation. In a fluid system, small density disturbances travel at the speed of sound. In gas systems, the speed of sound is generally determined by the temperature. The importance of the effect of compressibility of the flow is measured by the ratio of a characteristic velocity to the speed of sound, which is known as the Mach number.

[0058] A general discussion of an LBM-based simulation system, including a dynamics transform 46b for performing a fluid flow simulation, is provided below. For a further description of an LBM-based physics simulation system, see the above-referenced U.S. patents, which are incorporated by reference.

[0059] 5, the first model (2D-1) 200 is a two-dimensional model that includes 21 velocities. Of these 21 velocities, one velocity (205) represents a particle that is not moving, a set of three velocities represent particles moving in either the positive or negative direction along either the x or y axis of the lattice at either the normalized velocity (r) (210-213), twice the normalized velocity (2r) (220-223), or three times the normalized velocity (3r) (230-233), and a set of two velocities represent particles moving at either the normalized velocity (r) (240-243) or twice the normalized velocity (2r) (250-253) relative to both the x and y lattice axes.

[0060] Referring to Figure 6, a second model (3D-1) 260 is shown, which is a three-dimensional model that includes 39 velocities, each velocity represented by one of the arrows in Figure 10. Of these 39 velocities, one velocity represents a particle that is not moving, a trio of six velocities represent particles moving at either a normalized velocity (r), twice the normalized velocity (2r), or three times the normalized velocity (3r) in either a positive or negative direction along the x, y, or z axis of the lattice, eight velocities represent particles moving at a normalized velocity (r) for all three lattice axes, x, y, and z, and twelve velocities represent particles moving at twice the normalized velocity (2r) for two lattice axes, x, y, and z.

[0061] More complex models can also be used, such as the 3D-2 model, which contains 101 velocities, and the 2D-2 model, which contains 37 velocities. For the 3D model 3D-2, which consists of 101 velocities, one velocity represents particles that are not moving (group 1), a triad of six velocities represent particles moving at either normalized velocity (r), twice the normalized velocity (2r), or three times the normalized velocity (3r) in either the positive or negative direction along the x, y, or z axis of the lattice (groups 2, 4, and 7), and a triad of eight velocities represent particles moving at normalized velocity (r), twice the normalized velocity (2r), or three times the normalized velocity (3r) for all three lattice axes, x, y, and z (groups 1, 2, and 3). The velocities represent particles moving at twice the normalized velocity (2r) along two of the x, y, and z lattice axes (Group 6); the velocities represent particles moving at the normalized velocity (r) and twice the normalized velocity (2r) along two of the x, y, and z lattice axes and not moving along the remaining axes (Group 5); and the velocities represent particles moving at the normalized velocity (r) along two of the x, y, and z lattice axes and three times the normalized velocity (3r) along the remaining axis (Group 9).

[0062] For the two-dimensional model 2D-2, which consists of 37 velocities, one velocity represents particles that are not moving (group 1); three sets of four velocities represent particles moving at either the normalized velocity (r), twice the normalized velocity (2r), or three times the normalized velocity (3r) in either the positive or negative direction along either the x or y axis of the lattice (groups 2, 4, and 7); two sets of four velocities represent particles moving at the normalized velocity (r) or twice the normalized velocity (2r) for both the x and y lattice axes; eight velocities represent particles moving at the normalized velocity (r) for one of the x and y lattice axes and twice the normalized velocity (2r) for the other axis; and eight velocities represent particles moving at the normalized velocity (r) for one of the x and y lattice axes and three times the normalized velocity (3r) for the other axis.

[0063] The LBM model described above provides a particular class of efficient and robust discrete velocity kinetic models for the numerical simulation of flows in both two and three dimensions. This class of models includes a particular set of discrete velocities, and weights associated with those velocities. The velocities correspond to Cartesian grid points in velocity space, which facilitates accurate and efficient implementation, particularly of the class of discrete velocity models known as lattice Boltzmann models. Such models can be used to simulate flows to a high degree of completeness.

[0064] 7, the operation of the physics simulation system is described by procedure 270 for simulating a physics process such as a fluid flow. Prior to flow simulation, a simulation space is modeled as a collection of voxels (272) using the CAD drawing described above. The simulation space is generated using a computer-aided design (CAD) program and a gap correction process for the CAD-generated drawing. For example, a CAD program can be used to draw an airfoil placed in a wind tunnel.

[0065] The grid resolution can be selected based on the Reynolds number of the system being simulated, which is related to the viscosity of the flow (ν), the characteristic length of an object in the flow (L), and the characteristic velocity of the flow (u). Re=uL / ν Equation (I.4)

[0066] The characteristic length of an object represents a large-scale feature of the object. For example, if the flow around a microdevice is simulated, the height of the microdevice can be considered the characteristic length. If the flow around a small region of the object (e.g., a car side mirror) is of interest, the resolution of the simulation can be increased or regions of increased resolution can be used around the region of interest. As the grid resolution increases, the dimensions of the voxels decrease.

[0067] The state space is f i (x,t), where f irepresents the number of elements or particles of state i per unit volume (i.e., the density of particles in state i) at the lattice sites described by the three-dimensional vector x at time t. For a known time increment, the number of particles is simply f i The combination of all states of the lattice sites is denoted as f(x).

[0068] The number of states is determined by the number of possible velocity vectors in each energy level. A velocity vector is an integer linear velocity in a space with three dimensions: x, y, and z. In multi-species simulations, the number of states increases.

[0069] Each state i represents a different velocity vector at a particular energy level (i.e., energy level 0, 1, or 2). The velocity of each state c i is expressed using the "velocity" in each of the three dimensions as follows: c i =(c ix ,c iv ,c iz ) Equation (I.5)

[0070] The state of energy level 0 represents a stationary particle that is not moving in any dimension, i.e., c stopped =(0,0,0). A state of energy level 1 represents a particle with a velocity of ±1 in one of the three dimensions and a velocity of 0 in the other two. A state of energy level 2 represents a particle with a velocity of ±1 in all three dimensions, or a velocity of ±2 in one of the three dimensions and a velocity of 0 in the other two.

[0071] Generating all possible permutations of the three energy levels gives a total of 39 possible states (one energy 0 state, six energy 1 states, eight energy 3 states, six energy 4 states, 12 energy 8 states, and six energy 9 states).

[0072] Each voxel (i.e., each lattice site) is represented by a state vector f(x). The state vector completely defines the status of the voxel and contains 39 entries. The 39 entries correspond to one energy 0 state, six energy 1 states, eight energy 3 states, six energy 4 states, twelve energy 8 states, and six energy 9 states. By using this velocity set, the system can produce Maxwell-Boltzmann statistics for the achieved equilibrium state vector.

[0073] For processing efficiency, voxels are grouped into 2x2x2 volumes called microblocks. The microblocks are organized to enable parallel processing of the voxels and to minimize the overhead associated with data structures. The shorthand notation for the voxels in a microblock is N i (n), where n represents the relative position of the lattice site within the microblock, n ∈ {0, 1, 2, ..., 7}.

[0074] FIG. 8 shows a microblock.

[0075] 9A and 9B, a surface S (FIG. 9A) is defined in a simulation space (FIG. 9B) by a facet F α It is represented as a collection of. S={F α} Equation (I.6) where α is an index naming a particular facet. Facets are not restricted to voxel boundaries, but the size of a facet is typically comparable to or slightly smaller than the size of the voxels adjacent to the facet, so that the facet affects a relatively small number of voxels. The facets are assigned properties for the purposes of implementing surface dynamics. Specifically, each facet F α is the unit normal (n α ), surface area (A α ), center position (x α ), and the facet distribution function (f i(α) has a total energy distribution function q i (α) is treated in the same way as the flow distribution of facet-voxel interaction.

[0076] 10, to improve processing efficiency, different levels of resolution can be used in different regions of the simulation space. Typically, the region 320 around the object 322 is of greatest interest, and therefore this region is simulated with the highest resolution. Because the effect of viscosity decreases with distance from the object, decreasing levels of resolution (i.e., enlarged voxel volumes) are used to simulate regions 324, 326 located at intervals of increasing distance from the object 322.

[0077] 11, a lower level of resolution can be used to simulate regions 340 around less important features of the object 342, while the highest level of resolution is used to simulate regions 344 around the most important features (e.g., leading and trailing faces) of the object 342. An outer region 346 is simulated using the lowest level of resolution and the largest number of voxels.

[0078] C. Identification of voxels influenced by facets Referring again to FIG. 11, once the simulation space is modeled (272), voxels that are influenced by one or more facets are identified (274). Voxels may be influenced by facets in various ways. First, voxels that are intersected by one or more facets are affected in that the volume of the voxel is reduced compared to non-intersected voxels. This occurs because the facets and the subsurface material represented by the facets occupy a portion of the voxel. A fractional factor P f (x) indicates the portion of the voxel that is not influenced by the facet (i.e., the portion that may be occupied by the fluid or other material being simulated for flow). For non-intersecting voxels, P f(x) is equal to 1.

[0079] Voxels that interact with one or more facets by transferring particles to the facet or receiving particles from the facet are also identified as voxels influenced by the facet. Every voxel intersected by a facet will contain at least one state that receives particles from the facet and at least one state that transfers particles to the facet. In most cases, additional voxels will also contain such states.

[0080] Referring to FIG. 12, each state i has a non-zero velocity vector c i If F has α is a parallelepiped G iα and receiving particles from or transferring particles to a region defined by a parallelepiped G iα is the velocity vector c i and the unit normal of the facet n α Vector dot product with (|c i n α |) and the surface area of ​​the facet A α and thus the parallelepiped G iα Volume V iα teeth, V iα =|c i n α |A α Equation (I.7) is equal to.

[0081] Facet F α is the velocity vector of the state pointing towards the facet (|c i n α |<0), if the volume V iα receives a particle from the facet, and the velocity vector of the state is pointing away from the facet (|c i n α |>0), transfer the particle to that region. As explained below, if another facet is in the G iαThis formula is modified when

[0082] Facet F α of parallelepiped G iα may overlap some or all of several voxels. The number of voxels or parts thereof depends on the size of the facet relative to the size of the voxel, the energy of the state, and the orientation of the facet with respect to the lattice structure. The number of affected voxels increases with the size of the facet. Therefore, as mentioned above, the size of the facet is usually chosen to be comparable to or smaller than the size of the voxels located near the facet.

[0083] parallelepiped G iα The part of the voxel N(x) that overlaps with V iα (x). Using this term, we can distinguish between voxel N(x) and facet F α The flux Γ of particles in state i moving between iα (x) is the density of particles of state i in a voxel (N i (x)) and the volume of the area overlapping the voxel (V iα (x)). Gamma iα (x)=N i (x)V iα (x) Equation (I.8)

[0084] Parallelepiped G iα If intersects with one or more facets, then the following condition: V iα =ΣV iα (x) + ΣV iα (β) Equation (I.9) is true, where the first sum is G iα corresponds to all voxels that overlap with G, and the second term corresponds to iα corresponds to all facets that intersect with the parallelepiped G iα If does not intersect with another facet, then the expression V iα =ΣV iα (x) Equation (I.10) Scale down to.

[0085] D. Running the Simulation Once voxels influenced by one or more facets have been identified (274), a timer is initialized (276) to begin the simulation. During each time increment of the simulation, the movement of particles from voxel to voxel is simulated by an advection phase (278-286), which corresponds to the interaction of the particles with the surface facets. Next, a collision phase (288) simulates the interaction of particles within each voxel. The timer is then incremented (290). If the incremented timer does not indicate that the simulation is complete (292), the advection and collision phases (278-290) are repeated. If the incremented timer indicates that the simulation is complete (292), the results of the simulation are stored and / or displayed (294).

[0086] 1. Surface boundary conditions To correctly simulate the interaction with the surface, each facet satisfies four boundary conditions. First, the total mass of particles received by the facet is equal to the total mass of particles transferred by the facet (i.e., the net mass flux to the facet is equal to zero). Second, the total energy of particles received by the facet is equal to the total energy of particles transferred by the facet (i.e., the net energy flux to the facet is equal to zero). These two conditions can be satisfied by requiring that the net mass flux at each energy level (i.e., energy levels 1 and 2) is equal to zero.

[0087] The other two boundary conditions are related to the net momentum of the particle interacting with the facet. For a frictionless surface, referred to here as a slip surface, the net tangential momentum flux is equal to 0 and the net normal momentum flux is equal to the local pressure at the facet. Thus, the normal to the facet, n α and the component of the total momentum received perpendicular to the facet normal n αThe component of the total momentum transferred perpendicular to the facet (i.e., the tangential component) is equal to α The component of the total momentum received parallel to the normal of the facet, n α The difference between the component of the total momentum transferred that is parallel to (i.e., the normal component) is equal to the local pressure at the facet. For non-slip surfaces, surface friction reduces the total tangential momentum of the particles transferred by the facet compared to the total tangential momentum of the particles received by the facet by a factor related to the amount of friction.

[0088] 2. Voxel-to-facet collection In simulating the interaction between the particles and the surface, particles are collected from voxels and provided to facets (278). As described above, a voxel N(x) and a facet F α The flux of particles in state i between Gamma iα (x)=N i (x)V iα (x) Equation (I.11) It is.

[0089] From this equation, each state i is a facet F α (c i n α <0), the voxel is faceted by F α The number of particles provided to Gamma iαV→F =Σ X Gamma iα (x)=Σ X N i (x)V iα (x) Equation (I.12) It is.

[0090] V iα Only voxels with non-zero values ​​of (x) are summed. As mentioned above, the size of a facet is determined by V iα (x) is chosen so that it has non-zero values ​​for only a small number of voxels. iα (x) and P f Since (x) can have non-integer values, Γ α(x) is stored and manipulated as a real number.

[0091] 3. Moving from facet to facet The particle then moves between the facets (280). α is the inflow state (c i n α <0) iα is another facet F β If it intersects with facet F α Some of the particles in state i that are received by come from facet Fβ. Specifically, facet F α is the facet F during the previous time increment. β This relationship is shown in Figure 15, where the parallelepiped G iα Facet F β The intersection 380 is a parallelepiped G iβ Facet F α As mentioned above, the intersection is equal to the intersection 382 of V iα (β). Using this term, the facet F β and facet F α The flux of particles in state i between Gamma iα (β,t-1)=Γ i (β)V iα (β) / V iα Equation (I.13) where Γ i (β,t-1) is the facet F during the previous time increment. β From this equation, each state i is a measure of the particle in state i produced by facet F α (c i n α <0), facet F is determined by other facets. α The number of particles provided to Gamma iαF→F =Σ β Gamma iα (β)=Σ β Gamma i (β,t-1)V iα(β) / V iα Equation (I.14) and the total flux of particles of state i into the facet is Gamma iIN (α)=Γ iαF→F +Gamma iαF→F =Σ x N i (x)V iα +Σ β Gamma i (β,t-1)V iα (β) / V iα Equation (I.15) It is.

[0092] The facet state vector N(α), also called the facet distribution function, has M entries corresponding to the M entries of the voxel state vector, where M is the number of discrete lattice velocities. The input states of the facet distribution function N(α) are the flux of particles entering these states over a volume V iα is set equal to c divided by i n α If < 0, N i (α)=Γ iIN (α) / V iα Equation (I.16) It is.

[0093] The facet distribution function is a simulation tool for generating the output flux from a facet and does not necessarily represent an actual particle. To generate an accurate output flux, assign values ​​to the other states of the distribution function. Using the technique for assigning inward states described above, assigning the outward states results in c i n α If ≧0, N i (α)=Γ iOTHER (α) / V iα Equation (I.17) where Γ iOTHER (α) is the above-mentioned Γ iIN (α), but the inflow state (c i n α <0) other than the state (c i nα 0). In an alternative approach, we use the Γ iOUT Using the value of (α), iOTHER (α), Gamma iOTHER (α,t)=Γ iOUT (α,t-1) Equation (I.18) It can be generated as follows.

[0094] Parallel state (c i n α =0), V iα and V iα (x) is always 0. i In the formula for (α), V iα (x) is (Γ iOTHER (α) appears in the numerator of the equation, V iα is (N i (α) appears in the denominator of the equation for finding the parallel state. i (α) is V iα and V iα N when (x) approaches 0 i The values ​​for the states with zero velocity (i.e., the rest state, as well as states (0,0,0,2) and (0,0,0,-2)) are initialized at the start of the simulation based on the initial conditions of temperature and pressure. These values ​​are then adjusted over time.

[0095] 4. Running Faceted Surface Dynamics Next, surface dynamics is performed 282 for each facet such that it satisfies the boundary conditions above. Figure 14 shows the procedure 390 for performing surface dynamics for a facet.

[0096] 278, velocity and density from voxels to facets F α The velocity is then projected along the surface (i.e.

number

number

[0097] The momentum difference is projected along the tangent direction of the facet (398).

number

[0098] From the momentum difference and the Boltzmann distribution above, the zero tangent flux

number

number

[0099] To account for skin friction and other factors, the outflow flux distribution is further

number

[0100] 5. Voxel to Voxel Movement Referring again to Figure 7, particles are moved 284 from voxel to voxel along a 3D rectilinear grid. This voxel-to-voxel movement is the only movement operation performed for voxels that do not interact with facets (i.e., voxels that are not located near a surface). In a typical simulation, voxels that are not located near enough to a surface to interact with the surface make up the majority of voxels.

[0101] Each discrete state represents a particle moving with an integer velocity along the lattice in each of the three dimensions x, y, and z. Integer velocities include 0, ±1, and ±2. The sign of the velocity indicates the direction the particle is moving along the corresponding axis.

[0102] For voxels that do not interact with surfaces, the movement operation is computationally quite simple. During every time increment, the entire population of states is moved from its current voxel to its destination voxel. At the same time, the particles in the destination voxel are moved from that voxel to their own destination voxel. For example, a particle at energy level 1 moving in the +1x and +1y directions (1,0,0) will be moved from its current voxel to a voxel that is +1 in the x direction and 0 in the other directions. The particle will eventually arrive at its destination voxel in the same state (1,0,0) as it was before the movement. Interactions within the voxel may change the number of particles in this state based on local interactions with other particles and surfaces. If not, the particle will continue to move along the grid with the same speed and direction.

[0103] For voxels that interact with one or more surfaces, the translation operation becomes slightly more complex. This can result in one or more fractional particles being transferred to the facet. When such fractional particles are transferred to the facet, they will remain in the voxel. These fractional particles are transferred to the voxels that the facet occupies.

[0104] 13, when a portion 360 of state i particles in voxel 362 moves 278 to facet 364, a remaining portion 366 moves to voxel 368 in which facet 364 is located, and state i particles are directed from voxel 368 to facet 364. Thus, if the population of states is equal to 25, and V iα (x) is equal to 0.25 (i.e., one quarter of the voxels are parallelepiped G iα , 6.25 particles are intersected by facet F α 18.75 particles are moved to facet F α Since multiple facets may intersect a single voxel, the number of particles in state i that are transferred to a voxel N(f) occupied by one or more facets is given by

number

[0105] 6. Facet-to-voxel scattering The effluent particles from each facet are then scattered into voxels (286). Essentially, this scattering is the reverse of collection, where particles are transferred from voxels to facets. α The number of particles in state i that move from to voxel N(x) is

number

number

[0106] After scattering particles from facets into a voxel, combining them with particles advected from surrounding voxels, and integerizing the result, certain directions in certain voxels may underflow (go negative) or overflow (go beyond 255 in 8-bit implementations). This results in mass, momentum, and energy amounts being truncated to fit within the allowed range of values, resulting in an increase or decrease in mass, momentum, and energy. To protect against such events, the out-of-range mass, momentum, and energy are accumulated before truncating the problematic state. For the energy to which the state belongs, an amount of mass equal to the increased (by underflow) or decreased (by overflow) value is added to a randomly (or sequentially) selected state that has the same energy and is itself not subject to overflow or underflow. The additional momentum resulting from this addition of mass and energy is accumulated and added to the momentum from the truncation. By adding mass only to same-energy states, both mass and energy are corrected when the mass counter reaches 0. Finally, the momentum is corrected using a push / pull technique until the momentum accumulator returns to zero.

[0107] 7. Running Fluid Dynamics Perform fluid dynamics (288), Figure 7. This step is sometimes called microdynamics or intravoxel operations. Similarly, the advection step is sometimes called intervoxel operations. The microdynamics operations described below can also be used to collide particles at the facets to create a Boltzmann distribution.

[0108] Fluid dynamics is guaranteed in the Lattice Boltzmann Equation model by a specific collision operator known as the BGK collision model. This collision model mimics the distribution dynamics in real fluid systems. The collision process can be adequately described by the right-hand side of Equation 1 and Equation 2. After the advection step, the conserved quantities of the fluid system, specifically density, momentum, and energy, are obtained from the distribution function using Equation 3. From these quantities, the f in Equation (2) can be calculated aseq The equilibrium distribution function, denoted by is fully specified in equation (4). The set of velocity vectors, c i ,The ,choices ,of ,both ,weights ,are ,listed ,in ,Table ,1 ,and ,together ,with ,Equation ,2 ,ensure ,that ,the ,macroscopic ,behavior ,follows ,the ,correct ,hydraulic ,equations.

[0109] Variable Resolution Variable resolution (as described in US Patent Application Publication No. 2013 / 0151221) can also be used, resulting in the use of voxels of different sizes, e.g., coarse and fine voxels.

[0110] The subject matter and functional operations described herein may be implemented in digital electronic circuitry, tangibly embodied computer software or firmware, computer hardware (including the structures disclosed herein and their structural equivalents), or in a combination of one or more of these. The subject matter described herein may be implemented as one or more computer programs (i.e., one or more modules of computer program instructions encoded on a tangible, non-transitory program carrier for execution by or to control the operation of a data processing apparatus). The computer storage medium may be a machine-readable storage device, a machine-readable storage substrate, a random or serial access memory device, or a combination of one or more of these.

[0111] The term "data processing apparatus" refers to data processing hardware and encompasses any kind of apparatus, device, and machine for processing data, including, by way of example, a programmable processor, a computer, or multiple processors or computers. An apparatus can also be or further include special purpose logic circuitry (e.g., an FPGA (field programmable gate array) or an ASIC (application specific integrated circuit)). In addition to hardware, an apparatus can optionally include code that creates an execution environment for a computer program (e.g., code that constitutes a processor firmware, a protocol stack, a database management system, an operating system, or one or more combinations thereof).

[0112] A computer program, which may also be referred to or described as a program, software, software application, module, software module, script, or code, may be written in any form of programming language, including compiled or interpreted languages, or declarative or procedural languages, and may be deployed in any form, such as a stand-alone program or as a module, component, subroutine, or separate unit suitable for use in a computing environment. A computer program may correspond to a file in a file system, but need not. A program may be stored in part of a file that holds other programs or data (e.g., one or more scripts stored in a markup language document, in a single file dedicated to the program, or in multiple cooperating files (e.g., a file that stores one or more modules, subprograms, or portions of code)). A computer program may be deployed such that the program is executed on one computer or on multiple computers that are located at one site or distributed across multiple sites and interconnected by a data communications network.

[0113] Computers suitable for executing computer programs can be based on general purpose or dedicated microprocessors or both, or any other type of central processing unit. Computer-readable media suitable for storing computer program instructions and data include all forms of non-volatile memory on media and memory devices, including, by way of example, semiconductor memory devices (e.g., EPROM, EEPROM, and flash memory devices), magnetic disks (e.g., internal hard disks or removable disks), magneto-optical disks, and CD-ROM and DVD-ROM disks. The processor and memory may be supplemented by, or incorporated in, special purpose logic circuitry.

[0114] An embodiment of the subject matter described herein may be implemented in a computing system that includes a back-end component (e.g., as a data server), includes a middleware component (e.g., an application server), or includes a front-end component (e.g., a client computer having a graphical user interface or a web browser through which a user can interact with an implementation of the subject matter described herein), or any combination of one or more such back-end, middleware, or front-end components. The components of the system may be interconnected by any form or medium of digital data communication (e.g., a communication network). Examples of communication networks include a local area network (LAN) and a wide area network (WAN) (e.g., the Internet).

[0115] A computing system may include clients and servers. Clients and servers are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship. In some embodiments, a server sends data (e.g., HTML pages) to a user device acting as a client (e.g., for the purpose of displaying the data and receiving user input from a user interacting with the user device). Data generated at the user device (e.g., results of user interaction) can be received at the server from the user device.

[0116] Specific embodiments of the subject matter have been described. Other embodiments are within the scope of the following claims. For example, the actions recited in the claims can be performed in a different order and still achieve desirable results. As an example, the processes depicted in the accompanying figures do not necessarily require the particular order depicted or sequential order to achieve desirable results. In some cases, multitasking and parallel processing may be advantageous. [Explanation of symbols]

[0117] 10. System 11 Bus System 12 Large-scale parallel computing systems / server systems 15 Network 18 Memory 20 Interface 24 Processing Device 30 User-supplied mesh definition 32 Mesh Preparation Engines 34 Simulation Engine 34a Collision Interaction Module 34b Surface Dynamics Transformation 34c Boundary Module 34d Advection-Particle Collision Interaction Module 38 Data Repositories

Claims

1. 1. A computer-implemented method for simulating elements of a fluid flow about a representation of a surface of a physical object, comprising: storing in a memory a state vector of a plurality of voxels, the state vector including a plurality of entries corresponding to particular momentum states of a plurality of possible momentum states for the voxels; an entry of the state vector represents a density of elements per unit volume in a particular momentum state of said voxel; storing in a memory a representation of at least one surface of said physical object that is sized and oriented independent of the size and orientation of said voxels; the representation of the at least one surface includes a plurality of facets, each of the plurality of facets being assigned a surface characteristic; performing an interaction operation on the state vector, the interaction operation modeling interactions between elements of different momentum states; performing a surface interaction operation on the representation of the surface of the physics object, the surface interaction operation models interactions between the surface of the physics object and all elements of a voxel by sampling density from a fluid flow near the surface and rescaling the density such that a zero mass flux condition is satisfied on the surface; performing a translation operation on the state vector to reflect the translation of the element to a new voxel; 4. A computer-implemented method comprising:

2. The surface interaction calculation includes: collecting elements from a first set of at least one voxel that interacts with the facet; modeling interactions between the collected elements and all facets to generate a set of surface interaction elements; scattering the surface interaction element into a second set of at least one voxel that interacts with the facet; The method of claim 1 , comprising:

3. The method of claim 1 , further comprising sampling a velocity from the fluid stream near the surface.

4. 2. The method of claim 1, wherein the state vector comprises one or more of a plurality of integer and / or floating point values ​​that represent a density of elements per unit volume in a particular momentum state and have a larger number of possible values.

5. The method of claim 1 , wherein the surface interaction operations are performed using one or more of integer values ​​and / or floating point values, and the surface interaction operations are performed using values ​​representing real numbers.

6. The method of claim 1 , wherein the surface interaction operations are performed using integer and / or floating point numbers.

7. The method of claim 1 , further comprising representing voxels that intersect with a surface as subvoxels.

8. the representation of at least one surface includes a plurality of facets that are sized and oriented independent of the size and orientation of the voxels and represent the at least one surface; The method of claim 1 , wherein the surface interaction operation includes modeling an interaction between a facet and the element of at least one voxel near the facet.

9. The method of claim 1 , wherein the elements represent particles of a fluid and the facets represent at least one surface through which the fluid flows.

10. The step of performing a surface interaction operation comprises: calculating a difference between the inflow distribution and the Boltzmann distribution; determining the total momentum from the differences of all state vectors; generating a runoff distribution based on the determined difference; and The method of claim 1 further comprising:

11. 1. A data processing system for simulating elements of a fluid flow about a representation of a surface of a physical object, comprising: storing in a memory a state vector of a plurality of voxels, the state vector including a plurality of entries corresponding to particular momentum states of a plurality of possible momentum states at the voxels; an entry of the state vector represents a density of elements per unit volume in a particular momentum state of said voxel; storing in a memory a representation of at least one surface of said physical object that is sized and oriented independent of the size and orientation of said voxels; the representation of the at least one surface includes a plurality of facets, each of the plurality of facets being assigned a surface characteristic; performing an interaction operation on the state vector, the interaction operation modeling an interaction between elements of different momentum states; performing a surface interaction operation on the representation of the surface of the physics object, the surface interaction operation models an interaction between the surface of the physics object and all elements of a voxel by sampling density from a fluid flow near the surface and rescaling the density such that a zero mass flux condition is satisfied on the surface; performing a translation operation on the state vector to reflect the translation of the element to a new voxel; A data processing system comprising a memory having stored therein a program including instructions for causing the data processing system to perform the above.

12. The instructions for performing the surface interaction operation include: collecting elements from a first set of at least one voxel that interacts with the facet; modeling interactions between the collected elements and all facets to generate a set of surface interaction elements; scattering the surface interaction element into a second set of at least one voxel that interacts with the facet; 12. The data processing system of claim 11, further comprising instructions for:

13. 12. The data processing system of claim 11, wherein the state vector comprises one or more of a plurality of integer and / or floating point values ​​representing a density of elements per unit volume in a particular momentum state and having a greater number of possible values.

14. 12. The data processing system of claim 11, wherein the surface interaction operations are performed using one or more of integer values ​​and / or floating point values, and wherein the surface interaction operations are performed using values ​​representing real numbers.

15. 12. The data processing system of claim 11, wherein the surface interaction operations are performed using integer and / or floating point numbers.

16. The instructions for performing the surface interaction operation include:

12. The data processing system of claim 11, further comprising instructions for representing voxels that intersect with a surface as sub-voxels.

17. The representation of at least one surface includes a plurality of facets that are sized and oriented independent of voxel size and orientation and represent the at least one surface, and the instructions include: The data processing system of claim 11 , further comprising instructions for modeling an interaction between a facet and said element of at least one voxel near said facet.

18. 12. The data processing system of claim 11, wherein the elements represent particles of a fluid and the facets represent at least one surface through which the fluid flows.

19. The instructions for performing the surface interaction operation include: calculating a difference between the inflow distribution and the Boltzmann distribution; determining the total momentum from the differences of all state vectors; generating a runoff distribution based on the determined difference; and 12. The data processing system of claim 11, further comprising instructions for:

20. 1. A tangible, non-transitory computer readable medium storing a computer program product for simulating elements of a fluid flow about a representation of a surface of a physical object, the computer program product comprising: storing in a memory a state vector of a plurality of voxels, the state vector including a plurality of entries corresponding to particular momentum states of a plurality of possible momentum states at the voxels; an entry of the state vector represents a density of elements per unit volume in a particular momentum state of said voxel; storing in a memory a representation of at least one surface of said physical object that is sized and oriented independent of the size and orientation of said voxels; performing an interaction operation on the state vector, the interaction operation modeling an interaction between elements of different momentum states; the representation of the at least one surface includes a plurality of facets, each of the plurality of facets being assigned a surface characteristic; performing a surface interaction operation on the representation of the surface of the physics object, the surface interaction operation models an interaction between the surface of the physics object and all elements of a voxel by sampling density from a fluid flow near the surface and rescaling the density such that a zero mass flux condition is satisfied on the surface; performing a translation operation on the state vector to reflect the translation of the element to a new voxel; 23. A computer readable medium comprising instructions for causing a data processing system to:

21. The instructions for performing the surface interaction operation include: collecting elements from a first set of at least one voxel that interacts with the facet; modeling interactions between the collected elements and all facets to generate a set of surface interaction elements; scattering the surface interaction element into a second set of at least one voxel that interacts with the facet; 21. The computer readable medium of claim 20, further comprising instructions for:

22. 21. The computer-readable medium of claim 20, wherein the state vector comprises one or more of a plurality of integer and / or floating point values ​​that represent a density of elements per unit volume in a particular momentum state and have a greater number of possible values.

23. 21. The computer-readable medium of claim 20, wherein the surface interaction operations are performed using one or more of integer values ​​and / or floating point values, and wherein the surface interaction operations are performed using values ​​representing real numbers.

24. The instructions for performing a surface interaction operation include: calculating a difference between the inflow distribution and the Boltzmann distribution; determining the total momentum from the differences of all state vectors; generating a runoff distribution based on the determined difference; and 21. The computer readable medium of claim 20, further comprising instructions for:

Images