HK1224765A1 - Temperature coupling algorithm for hybrid thermal lattice boltzmann method - Google Patents

Description

Temperature coupling algorithm for hybrid thermal lattice boltzmann method

Priority protection

This application, which is incorporated herein by reference in its entirety, claims 35u.s.c. § 119(e) priority to U.S. provisional patent application No.61/860,392 filed 2013, 7, 31.

Background

Lattice Boltzmann Method (LBM) (or thermal Lattice Boltzmann Method (LBM)) is a class of Computational Fluid Dynamics (CFD) methods used for fluid simulation. Instead of solving the Navier-Stokes equations, discrete Boltzmann equations are solved to simulate the flow of Newtonian fluids using a collision model such as Bhatnagar-Gross-Krook (BGK). By simulating the fluidization and collision processes across a limited number of particles, the inherent particle interactions indicate a miniature of viscous flow behavior applicable across a greater mass.

Disclosure of Invention

In general, this document describes techniques for simulating the transport of particles in a fluid volume in a set of lattice velocities, where the transport causes collisions between the particles; and generating a distribution function for the transport of the particles, wherein the distribution function comprises a thermodynamic step and a particle collision step, and wherein the thermodynamic step is substantially independent of and separate from the particle collision step.

In some examples, the distribution function further comprises an advection step, and wherein the thermodynamic step is included in the distribution portion by enhancing the advection step with the thermodynamic step instead of enhancing the particle collision step with the thermodynamic step. In other examples, the thermodynamic step includes a temperature of the fluid volume during transport. In still other examples, generating includes: determining a post-collision distribution function f for a collision at a particular location x in a fluid volume at a particular time t_i' (x, t), wherein f_i′(x,t)＝f_i(x,t)+C_i(x, t) wherein C_iIs a collision operator, and f_iIs a distribution function for the particles prior to collision; from post-collision distribution function f_i' (x, t) deduction of fractional block g_i(x, t) to obtain a pre-advection particle density distribution function f_i″(x,t)＝f_i′(x,t)-g_i(x, t) wherein a part f of the particles_i′(x,t)-g_i(x, t) is advected to another location in the fluid volume, where g_i(x, t) represents the distribution of particles that are not advected; simulating advection of the particle portion at time t + Δ t to the other location in the fluid volume, the other location being denoted as (x + c)_iΔ t), wherein c_iIs the velocity vector of the particle before the collision, and Δ t is the interval between a certain time t and another point in time; obtaining the density distribution function of advected particles based on simulation of advectionWhereinAnd whereinIs at position x from position x-c_iOf advected particlesDistributing; block g to be deducted before_i(x, t) add back to density distribution functionTo form a density distribution function after advectionCalculating the mass, momentum and temperature of the particle at the location x at time t + Δ t; determining g using the calculated temperature, mass and momentum_i(x, t + Δ t); and the difference g_i(x,t+Δt)-g_i(x, t) addition to moving State f_i(x,t)+C_i(x,t)。

In some examples, g_iIs defined according to the following formula:

where ρ is the fluid density;

wherein T is₀Is the constant lattice temperature; wherein P is the pressure in the fluid volume; wherein T is the calculated temperature; and wherein w_iIs a constant weighting factor. The generated distribution function is based on:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)＋Ci(x，t)+[g_i(x+c_iΔt，t+Δt)-g_i(x，t)](ii) a Wherein x is the particle position within the volume; wherein t is a specific first point in time; where i is the index number of the lattice velocity in the set; wherein c is_iIs the velocity vector of the particle prior to collision; wherein C is_iIs a collision operator; wherein Δ t is the interval between the first point in time and the second point in time; wherein g is_iIs a thermodynamic step; and wherein f_iIs a distribution function for the particles at the location x at time t.

In some examples, the method conserves mass by modifying the stopping state to be according to:

in other examples, the particle collision step comprises an isothermal equilibrium distribution function. In other examples, the distribution function is for the location x + c_iThe distribution function of the particles at Δ t and at time t + Δ t and is expressed as:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t-g_i(x，t)](ii) a Wherein x is the particle position within the volume; wherein t is a specific first point in time; where i is the index number of the lattice velocity in the set; wherein C is_iIs a collision operator; wherein c is_iIs the velocity vector of the particle before collision; wherein Δ t is the interval between the first point in time and the second point in time; wherein g is_iIs a thermodynamic step; and wherein f_iIs a distribution function for the particles at the location x at time t. In still other examples, g_iIs defined according to the following formula:

where ρ is the fluid density; wherein T is₀Is the constant lattice temperature; wherein P is the pressure in the fluid volume; wherein T is the actual temperature of the fluid; and wherein w_iIs a constant weighting factor. In some examples, the set of lattice velocities is based on the lattice boltzmann method.

All or part of the foregoing may be implemented as a computer program product comprising instructions stored on one or more non-transitory machine-readable storage media and executable on one or more processing devices. All or a portion of the foregoing may be implemented as an apparatus, method, or electronic system that includes one or more processing devices and memory storing executable instructions that implement the described functionality.

The details of one or more implementations are set forth in the accompanying drawings and the description below. Other features, objects, and advantages will be apparent from the description and drawings, and from the claims.

These systems and methods and techniques may be implemented using various types of numerical simulation methods, such as the Shan-Chen method and the lattice boltzmann equation for multiphase flow. Further information regarding the lattice boltzmann equation will be described herein. However, the systems and techniques described herein are not limited to simulations using the lattice boltzmann equation, but may be applied to other numerical simulation methods.

These systems and techniques may be implemented using lattice gas simulations using the lattice boltzmann equation. Conventional lattice gas simulations assume a finite number of particles per lattice site, where the particles are represented by short vectors of bits. Each bit represents a particle moving in a particular direction. For example, a bit in the vector may indicate the presence (when set to 1) or absence (when set to 0) of a particle moving in a particular direction. Such a vector may have six bits, for example, a value of 110000 indicates that two particles are moving in opposite directions along the X axis, and that no particles are moving along the Y axis and the Z axis. A set of collision rules governs the behavior of collisions between particles at each site (e.g., a 110000 vector may become a 001100 vector, indicating that a collision between two particles moving along the X-axis results in two particles moving along the Y-axis). The rule is implemented by providing a state vector to a lookup table that performs a permutation on the bits (e.g., transform 110000 to 001100). The particles then move to contiguous sites (e.g., two particles moving along the Y-axis will move to adjacent sites to the left and right along the Y-axis).

In an enhanced system, the state vector at each lattice site includes many more bits (e.g., 54 bits for subsonic flow) to provide variation in particle energy and direction of movement, and a collision rule involving a subset of full state vectors is employed. In a further enhanced system, more than a single particle is allowed to exist in each momentum state at each lattice site or voxel (these two terms are used interchangeably throughout this document). For example, in an eight-bit implementation, 0-255 particles may move in a particular direction at a particular voxel. A state vector is a set of integers (e.g., a set of octets providing integers in the range of 0 to 255) rather than a set of bits, where each integer represents the number of particles in a given state.

In a further enhancement, the Lattice Boltzmann Method (LBM) uses a mesoscopic representation of the fluid to simulate a 3D unsteady compressible turbulent flow process in complex geometries at a deeper level than is possible with conventional computational fluid dynamics ("CFD") methods. A brief overview of the LBM method is provided below.

Boltzmann-hierarchical mesoscopic representation

It is well known in statistical physics that fluid systems can be represented by kinetic equations on a so-called "mesoscopic" level. At this level, detailed movements of individual particles need not be determined. In contrast, the properties of a fluid are represented by a particle distribution function defined using a single particle phase space, f ═ f (x, v, t), where x is the spatial coordinate and v is the particle velocity coordinate. Typical hydrodynamic quantities, such as mass, density, fluid velocity and temperature, are simple moments of the particle distribution function. The kinetics of the particle distribution function obey the boltzmann equation:

where F (x, t) represents the externally or self-consistent generated body force at (x, t). The collision term C represents the interaction of the particles at the respective velocities and positions. It is important to emphasize that without specifying a particular form for the collision term C, the above boltzmann equation applies to all fluid systems, and not just to the well-known case of lean gases (as originally constructed by boltzmann).

In general, C comprises a complex multidimensional integral of a two-point correlation function. To form a closed system with a single distribution function f and for the purpose of efficient computation, one of the most convenient and physically consistent forms is the well-known BGK operator. The BGK operator is approximated by { f via collisions according to a distribution function (regardless of the details of the collision) { F }^eq(x, v, t) } given the physical argument of well-defined local equilibrium:

where the parameter τ represents the characteristic relaxation time to reach equilibrium via collision. For particles (e.g., atoms or molecules), the relaxation time is typically taken to be constant. In the "mixed" (hydrodynamic) representation, this relaxation time is a function of hydrodynamic variables like strain rate, turbulent kinetic energy and others. Thus, turbulence can be expressed as a gas with turbulent particles ("vortices") of locally determined characteristic properties.

The numerical solution of the boltzmann-BGK equation has several computational advantages over the solution of the Navier-Stokes equation. First, it can be immediately recognized that there are no complex non-linear terms or higher order spatial derivatives in the equation, and thus there are few problems with advection instability. At this description level, the equation is local, as no processing pressure is required, which provides considerable advantages for algorithm parallelization. Another desirable feature of the linear advection operator, along with the absence of a diffusion operator having a second-order spatial derivative, is that it readily achieves physical boundary conditions, such as no-slip surfaces or slip surfaces, in a manner that mimics how particles actually interact with a real, solid surface, rather than achieving mathematical conditions for fluid partial differential equations ("PDEs"). One of the direct benefits is that there is no problem of handling the movement of interfaces on solid surfaces, which helps enable lattice-boltzmann-based simulation software to successfully simulate complex turbulent aerodynamics. Furthermore, certain physical properties from the boundary, such as a limited roughness surface, may also be incorporated in the force. Furthermore, the BGK collision operator is only local, and the calculation of self-consistent host forces can be done via neighbor information only. Therefore, the computation of the boltzmann-BGK equation can be efficiently adapted to parallel processing.

Lattice boltzmann equation

Solving successive boltzmann equations represents a significant challenge in that it involves numerical evaluation of integral-differential equations in position and velocity phase space when it is observed that not only position but also velocity phase space can be discretized, a great simplification occurs, which results in an efficient numerical algorithm for the solution of boltzmann equations_iI 1.., b }), which can also be systematically derived from the first principles as a discretization of the continuous boltzmann equation. As a result, LBE does not suffer from the well-known problems associated with the lattice gas approach. Thus, instead of processing a continuous distribution function f (x, v, t) in the phase space, only a limited set f of discrete distributions need to be tracked_i(x, t), where the subscript denotes the discrete speed index. A key advantage of dealing with this kinetic equation rather than a macroscopic description is that the increased phase space of the system is offset by the locality of the problem.

Due to symmetry considerations, the set of velocity values is selected in such a way that: in such a way that they form a certain lattice structure when spanned in the configuration space. The dynamics of such a discrete system follow a path having the form f_i(x+c_i,t+1)-f_i(x,t)＝C_i(x, t) LBE, where the collision operator typically takes the form of BGK as described above. By appropriate choice of the form of the equilibrium distribution, it can be shown theoretically that the lattice boltzmann equation yields the correct fluid dynamics and thermo-fluid dynamics. I.e. from f_i(x, t) derived hydrodynamic moments obey the Navier-Stokes equation in macroscopic constraints. These moments are defined as:

where ρ, u, and T are the fluid density, velocity, and temperature, respectively, and D is the dimension of the discretized velocity space (not equal to the physical space dimension at all).

Other features and advantages will be apparent from the following description, including the drawings and claims.

Drawings

Fig. 1 and 2 show the velocity components of two LBM models.

FIG. 3 is a flow chart of a procedure followed by the physical process simulation system.

Fig. 4 is a perspective view of a micro-tile.

Fig. 5A and 5B are illustrations of lattice structures used by the system of fig. 3.

Fig. 6 and 7 illustrate variable resolution techniques.

Fig. 8 shows the area affected by the facets of the surface.

Figure 9 shows the movement of a particle from a voxel to a surface.

Fig. 10 shows the movement of particles from surface to surface.

FIG. 11 is a flow chart of a procedure for performing surface dynamics.

FIG. 12 is a flow chart of a process for generating a distribution function having a thermodynamic step independent of the particle collision step.

FIG. 13 is a block diagram of components of a system for generating a distribution function having a thermodynamic step independent of a particle collision step.

Detailed Description

A. Overview

A system consistent with the present disclosure couples temperature to the LBM by removing the temperature coupling from the collision step and introducing the temperature coupling to the LBM as a separate thermodynamic step of binding temperature. By introducing a thermodynamic step into the distribution function (generated with LBM), this distribution function can be used to represent simulations for high mach and high temperature range applications (since temperature is not coupled into LBM by modifying the first order term that imposes the temperature range limitation).

B. Modeling simulation space

In an LBM-based physical process simulation system, fluid flow may be modeled by a set of discrete velocities c_iDistribution function value f to be evaluated_iAnd (4) showing. The dynamics of the distribution function are governed by equation 4, where f_i(0) Known as the equilibrium distribution function, which is defined as:

wherein

This equation describes the distribution function f_iWell-known lattice boltzmann equation of time evolution. The left hand side shows the change in distribution due to the so-called "fluidization process". The fluidization process is that a piece of fluid is at the grid positionFrom where it then moves along one of the velocity vectors to the next grid location. At that point, a "collision operator" is calculated, i.e., the effect of nearby fluid blocks on the originating fluid block. The fluid can only move to another grid location and therefore proper selection of the velocity vector is necessary so that all components of all velocities are multiples of a common rate (common).

The right hand side of the first equation is the "collision operator" mentioned above, which represents the change in the distribution function due to collisions between fluid blocks. The particular form of collision operator used here is attributed to Bhatnagar, Gross and Krook (BGK). It forces the distribution function to the prescribed value given by the second equation, which is a "balanced" form.

From this simulation, conventional fluid variables, such as mass ρ and fluid velocity u, are obtained as a simple sum in equation (3). Here, c_iAnd w_iThe set of values of (a) defines the LBM model. The LBM model can be efficiently implemented on a scalable computer platform and operates with maximum robustness to time-invariant flows and complex boundary conditions.

The standard technique for obtaining macroscopic equations for the motion of a fluid system from the boltzmann equation is the Chapman-Enskog method, in which a successive approximation of the full boltzmann equation is taken.

In a fluid system, small perturbations in density 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 in a fluid is measured by the ratio of the characteristic velocity to the speed of sound, which is referred to as mach number.

Referring to fig. 1, a first model (2D-1)100 is a two-dimensional model including 21 velocities. Of these 21 velocities, one velocity (105) represents a particle that is not moving; the three sets of four velocities represent particles moving at the normalized rate (r) (110-; and two sets of four velocities represent particles moving at the normalized rate (r) (140-.

As also shown in FIG. 2, the second model (3D-1)200 is a three-dimensional model that includes 39 velocities, where each velocity is represented by one of the arrows of FIG. 2. Of these 39 velocities, one velocity represents a particle that is not moving; three sets of six velocities represent particles moving at a normalized rate (r), twice the normalized rate (2r), or three times the normalized rate (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 rate (r) with respect to all three of the x, y, z lattice axes; and twelve velocities represent particles moving at twice the normalized rate (2r) with respect to two of the x, y, z lattice axes.

More complex models may also be used, such as a 3D-2 model comprising 101 velocities and a 2D-2 model comprising 37 velocities.

For the three-dimensional model 3D-2, of 101 velocities, one velocity represents a particle that does not move (group 1); three sets of six velocities represent particles moving at a normalized rate (r), twice the normalized rate (2r), or three times the normalized rate (3r) in either a positive or negative direction along the x, y, or z axis of the lattice (sets 2, 4, and 7); three sets of eight velocities represent particles moving at a normalized rate (r), twice the normalized rate (2r), or three times the normalized rate (3r) relative to all three of the x, y, z lattice axes (sets 3, 8, and 10); twelve velocities represent particles moving at twice the normalized rate (2r) relative to two of the x, y, z lattice axes (group 6); twenty-four velocities represent particles that move at the normalized rate (r) and twice the normalized rate (2r) with respect to two of the x, y, z lattice axes and do not move with respect to the remaining axes (group 5); and twenty-four velocities represent particles moving at a normalized rate (r) with respect to two of the x, y, z lattice axes and three times the normalized rate (3r) with respect to the remaining axes (group 9).

For the two-dimensional mode 2D-2, of the 37 velocities, one velocity represents a particle that is not moving (group 1); three sets of four velocities represent particles (sets 2, 4 and 7) moving at a normalized rate (r), twice the normalized rate (2r) or three times the normalized rate (3r) in the positive or negative direction along the x or y axis of the lattice; two sets of four velocities represent particles moving at a normalized rate (r) or twice the normalized rate (2r) with respect to both the x-lattice axis and the y-lattice axis; eight velocities represent particles moving at a normalized rate (r) with respect to one of the x-lattice axis and the y-lattice axis and at twice the normalized rate (2r) with respect to the other axis; and eight velocities represent particles moving at a normalized rate (r) with respect to one of the x-lattice axis and the y-lattice axis and three times the normalized rate (3r) with respect to the other axis.

The LBM model described above provides a specific class of efficient and robust discrete velocity dynamics models for numerical simulation of flow in two and three dimensions. This type of model includes a particular set of discrete velocities and weights associated with those velocities. These velocities coincide with grid points of cartesian coordinates in the velocity space, which facilitates accurate and efficient implementation of discrete velocity models, particularly of the kind known as lattice boltzmann models. With such a model, the flow can be simulated with high fidelity.

Referring to FIG. 3, a physical process simulation system operates in accordance with a protocol 300 to simulate a physical process, such as a fluid flow. Prior to simulation, the simulation space is modeled as a set of voxels (step 302). Typically, the simulation space is generated using a Computer Aided Design (CAD) program. For example, a CAD program may be used to render a microdevice located in a wind tunnel. Thereafter, the data generated by the CAD program is processed to add lattice structures with appropriate resolution and to illustrate objects and surfaces in simulation space.

The resolution of the lattice may be selected based on the reynolds number (reynolds number) of the system being simulated. Reynolds number relates to the viscosity (v) of the flow, the characteristic length (L) of the object in the flow and the characteristic velocity (u) of the flow:

Re-uL/v equation (6)

The characteristic length of the object represents a large scale characteristic of the object. For example, if the flow around a microdevice is being simulated, the height of the microdevice may be considered the feature length. When flow around a small region of the object (e.g., a side view mirror of an automobile) is of interest, the resolution of the simulation may be increased, or a region of increased resolution may be employed around the region of interest. The dimensions of the voxels decrease as the resolution of the lattice increases.

The state space is denoted f_i(x, t) wherein f_iThe number of elements or particles per unit volume in state i at a lattice site represented by a three-dimensional vector x at time t (i.e., the density of particles in state i) is represented. For a known time increment, the number of particles is referred to as f_i(x) In that respect The combination of all states of a lattice site is denoted as f (x).

The number of states is determined by the number of possible velocity vectors within each energy level. The velocity vector consists of integer linear velocities in space with three dimensions x, y, and z. For multi-species simulations, the number of states increases.

Each state i represents a different velocity vector at a particular energy level (i.e., energy level zero, one, or two). Speed of each state c_iWith its "rate" in each of the three dimensions indicated as follows:

c_i＝(c_i,x,c_i,y,c_i,z) Equation (7)

The energy level zero state represents a stopped particle that does not move in any dimension, i.e., c_stopped(0,0, 0). An energy level one state represents a particle having a velocity of ± 1 in one of the three dimensions and zero velocity in the other two dimensions. An energy level two state represents a particle with a velocity of ± 1 in all three dimensions, or a velocity of ± 2 in one of the three dimensions and zero velocity in the other two dimensions.

All possible permutations that generate three energy levels give a total of 39 possible states (one energy zero state, 6 energy one states, 8 energy three states, 6 energy four states, 12 energy eight states and 6 energy nine states).

Each voxel (i.e., each lattice site) is represented by a state vector f (x). The state vector fully defines the state of the voxel and includes 39 entries. These 39 entries correspond to one energy zero state, 6 energy one states, 8 energy three states, 6 energy four states, 12 energy eight states, and 6 energy nine states. Using this set of velocities, the system can generate maxwell-boltzmann statistics on the implemented equilibrium state vectors.

For processing efficiency, voxels are grouped in volumes of 2x2x2 called micro-tiles. The micro-tiles are organized to allow parallel processing of voxels and to minimize overhead associated with the data structure. A shorthand notation for a voxel in a micro-tile is defined as N_i(n), where n represents the relative position of lattice sites in the microblock and n ∈ {0,1, 2.., 7 }. the microblock is shown in FIG. 4.

Referring to fig. 5A and 5B, the surface S is represented in simulated space (fig. 5B) as facet F_αSet of (2):

S＝{F_αequation (8)

Where α is an index enumerating particular facets, facets are not limited to the boundaries of voxels, but typically have dimensions on the order of the dimensions of voxels adjacent to the facet, or have dimensions slightly smaller than the dimensions of voxels adjacent to the facet, so that the facet affects a relatively small number of voxels_αWith unit normal (n)_α) Surface area (A)_α) Center position (x)_α) And a facet distribution function (f) describing surface dynamic properties of the facets_i(α))。

Referring to FIG. 6, different resolution levels may be used in different regions of simulation space to improve processing efficiency. Typically, the region 650 around the object 655 is of most interest and is therefore modeled with the highest resolution. Because the effect of viscosity decreases with distance from the object, a decreasing resolution level (i.e., an expanding volume of the body element) is used to simulate regions 660, 665 that are spaced at increasing distances from the object 655. Similarly, as shown in fig. 7, a lower resolution level may be used to simulate an area 770 around less prominent features of the object 775, while a highest resolution level is used to simulate an area 780 around most prominent features (e.g., leading and trailing edge surfaces) of the object 775. The outlying region 785 is modeled with the lowest resolution level and the largest voxel.

C. Identifying voxels affected by facets

Referring again to FIG. 3, once the simulation space has been modeled (step 302), voxels affected by one or more facets are identified (step 304). Voxels can be affected by facets in a variety of ways. First, voxels intersected by one or more facets are affected by: the voxel has a reduced volume relative to non-intersecting voxels. This occurs because the facet and the material below the surface represented by the facet occupy a portion of the voxel. Fractional factor P_f(x) Indicating the portion of the voxel that is not affected by the facet (i.e., the portion that may be occupied by a fluid or other material for which flow is simulated). For non-intersecting voxels, P_f(x) Equal to 1.

Voxels that intersect one or more facets by transmitting particles to the facets or receiving particles from the facets are also identified as voxels affected by the facets. All voxels intersected by the facet will include at least one state that receives the particle from the facet and at least one state that transmits the particle to the facet. In most cases, additional voxels will also include this state.

Referring to FIG. 8, for a vector c having a non-zero velocity_iEach state i, facet F_αFrom parallelepiped G_iαThe defined region receiving or transmitting particles theretoWherein the parallelepiped G_iαWith velocity vector c_iAnd facet (| c)_in_i|) unit normal n_αAnd the height defined by the magnitude of the vector dot product of (A) and the surface area of the facet (A)_αBase defined so that it is parallelepipedal G_iαVolume V of_iαEqual to:

V_iα＝|c_in_α|A_αequation (9)

When the velocity vector of a state points to a facet (| c)_in_i|<0) Facet F_αFrom volume V_iαReceive particles, and when the velocity vector of the state points away from the facet (| c)_in_i|>0) And transporting the particles to the region. When the other facet occupies the parallelepiped G, as will be discussed below_iαThat is, a condition that may occur near a non-convex feature such as an interior corner, this expression must be modified.

Facet F_αOf the parallelepiped G_iαSome or all of the voxels may be overlapped. The number of voxels or portions thereof depends on the facet size relative to the voxel size, the energy of the state, and the orientation of the facet relative to the lattice structure. The number of voxels affected increases with the size of the facet. Thus, as noted above, the size of a facet is typically selected to be on the order of or smaller than the size of voxels located near the facet.

Quilt parallelepipeds G_iαThe overlapping portion of voxels N (x) is defined as V_iα(x) In that respect Using this term, at voxel V_iα(x) And facet F_αFlux of moving state i particles_iα(x) Equal to the density of state i particles in the voxel (N)_i(x) Multiplied by the volume (V) of the region overlapping the voxel_iα(x))：

_iα(x)＝N_i(x)V_iα(x) Equation (10)

When parallelepipedon G_iαIntersected by one or more facets, the following condition is true:

V_iα＝∑V_α(x)+∑V_iα(β) equation (11)

Wherein the first summation is denoted by G_iαAll voxels overlapping and the second term specification with G_iαAll facets intersecting. When parallelepipedon G_iαWithout being intersected by another facet, this expression reduces to:

V_iα＝∑V_iα(x) Equation (12)

D. Performing simulations

Once the voxels affected by the one or more facets are identified (step 304), a timer is initialized to begin the simulation (step 306). During each time increment of the simulation, the movement of the particle from voxel to voxel is simulated by an advection phase that accounts for the interaction of the particle with the surface facet (step 308-316). Next, the collision phase (step 318) simulates the interaction of the particles within each voxel. Thereafter, the timer is incremented (step 320). If the incremented timer does not indicate that the simulation is complete (step 322), then the advection and collision phases (steps 308-320) are repeated. If the incremented timer indicates that the simulation is complete (step 322), the results of the simulation are stored and/or displayed (step 324).

1. Boundary conditions for surfaces

In order to correctly simulate the interaction with the surface, each facet must satisfy four boundary conditions. First, the combined mass of the particles received by the facet must equal the combined mass of the particles transmitted by the facet (i.e., the net mass flux to the facet must equal zero). Second, the combined energy of the particles received by the facet must be equal to the combined energy of the particles transmitted by the facet (i.e., the net energy flux to the facet must be equal to zero). These two conditions may be satisfied by requiring the net mass flux at each energy level (i.e., energy level one and energy level two) to be equal to zero.

The other two boundary conditions are related to the net momentum of the particles interacting with the facets. For a surface without skin friction (referred to herein as a sliding surface), the net tangential momentum flux must equal zero and the net normal momentum flux must equal the local pressure at the facet. Hence, normal n to the facet_αThe components of the momentum received and transmitted by the perpendicular combination (i.e., the tangential components) must be equal to the normal n to the facet_αThe difference between the components of the momentum received and transmitted by the parallel combination (i.e., the normal components) must be equal to the local pressure at that facet. For a non-slip surface, the friction of the surface reduces the combined tangential momentum of the particles conveyed by the facet by a factor related to the amount of friction relative to the combined tangential momentum of the particles received by the facet.

2. Concentration from voxels to facets

As a first step in simulating the interaction between the particle and the surface, the particle is gathered from the voxel and provided to the facet (step 308). As indicated above, voxel n (x) and facet F_αThe flux of the particles in state i in between is:

_iα(x)＝N_i(x)V_iα(x) Equation (13)

In this view, for the directional facet F_αEach state i (c)_in_α<0) Supplied by voxels to facets F_αThe number of particles of (a) is:

only its V_iα(x) Voxels with non-zero values must be summed. As indicated above, the dimensions of the facets are selected such that V_iα(x) Only for a few voxels have non-zero values. Because of V_iα(x) And P_f(x) May have non-integer values, therefore_α(x) Stored and processed as real numbers.

3. Movement from facet to facet

Next, the particles are moved between the facets (step 310). If used for facet F_α(ii) an incoming state (c)_in_α<0) Of the parallelepiped G_iαIs faceted with another facet F_βIntersect, then by F_αA portion of the received state i particles will come from facet F_β. In particular, facet F_αWill receive during the previous time increment by facet F_βPart of the resulting state i particles. This relationship is shown in FIG. 10, where faceted F_βIntersecting parallelepipeds G_iαIs equal to the faceted surface F_αIntersecting parallelepipeds G_iβPortion 1005. As noted above, the intersecting portion is denoted V_iα(β) Using this term, facet F_βAnd facet F_αThe flux of state i particles in between can be described as:

_iα(β,t-1)＝_i(β)V_iα(β)/V_iαequation (15)

Wherein_i(β, t-1) is defined by facet F during a previous time increment_βMeasurement of the resulting state i particles. In this view, for the directional facet F_αEach state i (c)_in_α<0) Provided to facet F by other facets_αThe number of particles of (a) is:

and the total flux to the state i particles in the facet is:

the state vector N (α) (also known as the facet distribution function) for a facet has a correspondenceThe input states of the faceted distribution function N (α) are set equal to the flux of particles into those states divided by the volume V_iα：

For c_in_α<0，N_i(α)＝_iIN(α)/V_iαEquation (18)

The facet distribution function is a simulation tool for generating output flux from the facets and does not necessarily represent actual particles. To generate an accurate output flux, values are assigned to other states of the distribution function. The outward state is filled using the technique described above for filling the inward state:

for c_in_α≥0，N_i(α)＝_iOTHER(α)/V equation (19) wherein_iOTHER(α) is a method for generating using the above_iIN(α), but applying the technique to divide incoming states (c)_in_α<0) Out of state (c)_in_αNot less than 0). In an alternative method of the invention,_iOTHER(α) may utilize information from a previous time step_iOTHERThe value of (α) is generated such that:

_iOTHER(α,t)＝_iOUT(α, t-1) equation (20)

For parallel state (c)_in_α＝0)，V_iαAnd V_iα(x) Are all zero. In the use for N_i(α) wherein V is_iα(x) Present in the molecule (according to the application)_iOTHERExpression of (α) and V_iαAppear in the denominator (according to the use for N_iExpression of (α). thus, when V_iαAnd V_iα(x) Near zero, N for parallel state_i(α) is determined to be N_i(α).

The values of the states with zero velocity, i.e. the rest state and states (0,0,0,2) and (0,0,0, -2), are initialized at the beginning of the simulation based on initial conditions for temperature and pressure. These values are then adjusted over time.

4. Performing faceted surface dynamics

Next, surface dynamics are performed on each facet to satisfy the four boundary conditions discussed above (step 312.) the procedure for performing surface dynamics for the facet is shown in FIG. 11 Initially, the to facet F is determined by determining the combined momentum of the particle at the facet P (α)_αCombined momentum of (step 1105):

for all of the i's, the average,

from this equation, the normal momentum P_n(α) is determined as:

P_n(α)＝n_αequation (22) P (α)

This normal momentum is then eliminated (step 1110) using a push/pull technique to produce N_n-(α). according to this technique, particles are moved between states in a manner that affects only normal momentum.

Thereafter, N_n-(α) the particles are collided to generate a Boltzmann distribution N_n-β(α) (step 1115) — the boltzmann distribution may be transformed into N, as described below with respect to performing fluid dynamics_n-(α) applying a set of collision rules.

Then, based on the incoming flux distribution and Boltzmann distribution, a determination is made for facet F_αThe outgoing flux profile of (step 1120). First, the incoming flux distribution_iThe difference between (α) and the boltzmann distribution is determined as:

Δ_i(α)＝_iIN(α)-N_n-βi(α)V_iαequation (23)

Using this difference, the outgoing flux distribution is:

for n_αc_i>0，_iOUT(α)＝N_n-βi(α)V_iα-.Δ._i*(α) equation (24) and where i is the state with the opposite direction to state i, for example, if state i is (1,1,0,0), state i is (-1, -1,0, 0).

For n_αc_i>0，

Wherein C is_fIs a function of skin friction, t_iαIs and n_αPerpendicular first tangent vector, t_2αIs and n_αAnd t_1αAll perpendicular second tangent vectors, and Δ N_j,1And Δ N_j,2Is a distribution function of the energy (j) corresponding to the state i and the indicated tangent vector. The distribution function is determined according to the following equation:

where j is equal to 1 for level 1 state and 2 for level 2 state.

For_iOUTThe first and second terms implement the normal momentum flux boundary condition to the extent that a collision has effectively produced a Boltzmann distribution but includes tangential momentum flux anomalies_fThe generation of (c) is described below. It should be noted that all terms relating to vector manipulation are possibleThe geometry factor calculated before the start of the simulation.

In this view, the tangential velocity is determined as:

u_i(α)＝(P(α)-P_n(α)n_α) /ρ, equation (27)

Where ρ is the density of the facet distribution:

as before, the difference between the incoming flux distribution and the boltzmann distribution is determined as:

Δ_i(α)＝_iIN(α)-N_n-βi(α)V_iαequation (29)

Then, the outgoing flux distribution becomes:

_iOUT(α)＝N_n-βi(α)V_iα-Δ_i*(α)+C_f(n_αc_i)[N_n-βi*(α)-N_n-βi(α)]V_iαequation (30)

Which corresponds to the first two rows of outgoing flux distributions determined by previous techniques, but does not require correction for anomalous tangential flux.

Using either approach, the resulting flux distribution satisfies all momentum flux conditions, namely:

wherein p is_αIs at facet F_αAnd based on the average density and temperature values of voxels providing particles to the facet, and u_αIs the average velocity at the facet.

To ensure that the quality and energy boundary conditions are met, the difference between the input energy and the output energy is measured for each energy level j:

where the index j represents the energy of state i. This energy difference is then used to generate a difference term:

for c_jin_α>0，

This difference term is used to modify the outgoing flux so that it becomes:

for c_jin_α>0，_αjiOUTf＝_αjiOUT+_αjiEquation (34)

This operation corrects for mass and energy flux while leaving the tangential momentum flux unaltered. This adjustment is small if the flow is approximately uniform and near equilibrium in the vicinity of the facet. After the adjustment, the resulting normal momentum flux is slightly altered to a value that is the equilibrium pressure based on the neighborhood average property plus the correction due to the non-uniformity or non-equilibrium properties of the neighborhood.

5. Moving from voxel to voxel

Referring again to FIG. 3, the particles move between voxels along a three-dimensional rectilinear lattice (step 314). This voxel-to-voxel movement is the only movement operation performed on voxels that do not interact with facets (i.e., voxels that are not located near the surface). In a typical simulation, voxels that are not located close enough to the surface to interact with the surface constitute the majority of voxels.

Each separate state represents a particle moving along the lattice at an integer rate in each of the three dimensions x, y, and z. Integer rates include: 0, ± 1 and ± 2. The sign of the velocity indicates the direction in which the particle is moving along the corresponding axis.

For voxels that do not interact with the surface, the move operation is computationally very simple. The entire population of states (population) moves from its current voxel to the destination voxel during each time increment. At the same time, the particle of the destination voxel is moved from that voxel to its own destination voxel. For example, an energy level 1 particle moving in the +1x and +1y directions (1,0,0) moves from its current voxel to a voxel that is +1 in the x direction and 0 for the other directions. The particle terminates at its destination voxel in the state (1,0,0) it had before moving. Based on local interactions with other particles and the surface, the interaction within the voxel will likely change the particle count for that state. If not, the particles will continue to move along the lattice at the same rate and direction.

For voxels that interact with one or more surfaces, the move operation becomes slightly more complex. This results in one or more fractional particles being transported to the facet. This transport of the fractional particles to the facet causes the fractional particles to remain in the voxel. These fractional particles are transmitted into voxels occupied by facets. For example, referring to FIG. 9, when a portion 900 of the state i-particle of voxel 905 is moved to facet 910 (step 308), the remaining portion 915 is moved to the voxel 920 in which facet 910 is located and from which the state i-particle is directed to facet 910. Thus, if the state population is equal to 25 and V_iα(x) Equal to 0.25 (i.e. one quarter of a voxel and parallelepiped G)_iαIntersect), then 6.25 particles will be moved to facet F_αAnd 18.75 particles will be moved to be faceted F_αOccupied voxels. Since multiple facets will intersect a single voxel, the number of state i particles transmitted to a voxel N (f) occupied by one or more facets is:

where N (x) is the source voxel.

6. Dispersion from facets to voxels

Next, the outgoing particles from each facet are dispersed to the voxels (step 316). Essentially, this step is the inverse of the aggregation step that moves the particles from the voxel to the facet. From F_αThe number of state i particles that the facet moves to voxel N (x) is:

wherein P is_f(x) Illustrating the volume reduction of a portion of the voxel. In this view, for each state i, a voxel N is directed from the facet_(x)The total number of particles of (a) is:

after dispersing the particles from facet to voxel, the particles are combined with the particles that have been advected from the surrounding voxels, and the result is integer-ized, with the possibility that some directions in some voxels may underflow (become negative) or overflow (in an eight-bit implementation, over 255). This will result in gain or loss of mass, momentum and energy after they have been truncated to fit within the allowed value range. To prevent this from happening, the out-of-bounds mass, momentum, and energy are accumulated prior to truncation of the violation state. For the energy to which the state belongs, an amount of quality equal to the value of gain (due to underflow) or loss (due to overflow) is added back to the randomly (or sequentially) selected state that has the same energy and is not itself subject to overflow or underflow. The additional momentum resulting from this addition of mass and energy is added up and to the momentum from truncation. By adding mass only to the same energy state, both mass and energy are corrected when the mass counter reaches zero. Finally, the momentum is corrected using a push/pull technique until the momentum accumulator returns to zero.

7. Performing fluid dynamics

Finally, fluid dynamics are performed (step 318). This step may be referred to as micro-dynamics or intra-voxel manipulation. Similarly, the advection procedure may be referred to as inter-voxel operation. The micro-kinetic operations described below can also be used to impact particles at the facets to produce a boltzmann distribution.

Fluid dynamics are ensured in the lattice boltzmann equation model by a specific collision operator known as the BGK collision model. This collision model mimics the dynamics of the distribution in a real fluid system. The collision process can be well described by the right hand side of equations 1 and 2. After the advection step, the conservation quantities of the fluid system, specifically the density, momentum and energy, are obtained from the distribution function using equation 3. From these quantities, the equation (2) is given by f^eqThe indicated equilibrium distribution function is completely specified by equation (4). Velocity vector set c_iThe choice of weights are listed in table 1, along with equation 2 to ensure that the macroscopic behavior follows the correct fluid dynamics equations.

E. Mixed thermal lattice boltzmann method

Lattice Boltzmann Method (LBM) is used as a replacement for conventional Computational Fluid Dynamics (CFD) for a wide range of industrial applications within the limits of being almost incompressible. For thermal applications, LB solvers with energy conservation are highly unstable due to the lack of H-theorem and require higher order models to provide accurate results. Applying boundary conditions to energy in addition to mass and momentum further requires complex algorithms. To avoid these difficulties, most LBM methods use a separate finite difference solver or lattice boltzmann scalar solver for the energy. For convenience, the numbering of the following equations will start over at (1.0).

The standard isothermal LBM is given by the following equation:

f_i(x+c_iΔt,t+Δt)＝f_i(x,t)+C_i(x,t)(1.0)

wherein f is_i(x, t) is for speed value c at (x, t)_iAnd C is a particle distribution function of_i(x, t) is the particle collision operator. For example, the most common collision operator has the form BGK:

in equation 1.1, f_i(x, t) is the particle distribution function and τ is the relaxation time.

Standard equilibrium distribution f_i ^eqGiven by the following formula,

wherein the constant weighting factor w_iAnd lattice temperature T₀Is based on lattice assembly. The purpose of the weighting factors is to satisfy the order up to expectationThe moment of (a) is isotropic. For example, for D3Q19, T₀Is (1/3) and the weighting factor satisfies moment isotropy up to the second order. The above isothermal LBM results in the equation of state P ═ ρ T₀. In this example, P is pressure (e.g., in Nm)^-2In units), ρ is the density (e.g., in kgm)^-3In units) and T₀Is a temperature constant.

In the case of the hybrid thermal LBM solver, temperature coupling back to LBM is performed by modifying the first order term of the equilibrium distribution such that the second moment results in the actual pressure. The modified LBM equation is

Wherein

In equation 1.4 above, P satisfies the actual equation of state for an ideal gas, P ═ ρ RT and R is the gas constant factor β at rest equilibriumGiven, and q is the total number of mobile lattices.

The above equilibrium modifications do not alter the zeroth and first order moments, and result in the correct thermodynamic pressure for the second order moment,

this method requires only minimal modification to conventional LBM and the stress obtained as the second moment includes all lattice directions. A key drawback of this algorithm is the temperature range limitation it imposes due to changes in the dominant order term. For example, for an ideal gas, the distribution of the stop states for the D3Q19 lattice set becomes

For high temperatures, RT>1.5T₀The stall state equilibrium becomes negative and thus causes instability. Furthermore, for low temperature RT<T₀The positive (probability) range of the mobile state decreases, which also leads to instability at high speeds. Due to this drawback, most hybrid thermal LBM models are limited to low temperature range and low mach number applications.

F. Methods utilizing force protocols

In another method, the difference between the pressures resulting from the isothermal LBM and the actual equation of state is calculated as follows:

the above item F is applied as an external main force. To avoid the discrete lattice effect of the host force, higher order force terms with modifications to the velocity definition must be employed. Since the force term added is a gradient in pressure, the positive range of equilibrium is not adversely affected. However, the calculation of such additional force terms requires a finite difference calculation that affects the isotropic properties of the LBM. The calculation of the gradient near the wall is also challenging, which can result in non-zero normal velocities for high heat flux simulations.

G. Temperature coupling algorithm

As discussed in the previous section, the conventional coupling method has its own advantages and disadvantages. In the case of a change of the dominant order of the equilibrium (section E), the pressure gradient obtained includes all the lattice directions and is easy to achieve. However, altering the balance reduces the positive range and therefore reduces the availability range of the LB solver. The introduction of the subject force term (part F) does not affect balance and therefore provides stability for a wide range of rates and temperatures. However, the computation of gradients with finite templates affects the isotropic nature of the LB scheme and introduces errors like grid dependency, lattice orientation dependency and complexity when applying the correct boundary conditions.

In this section, a temperature coupling algorithm is described to overcome the shortcomings of the prior methods. Linear stability analysis shows that whenever the equilibrium is negative, collisions lead to numerical instability. Additionally, Chapmann multiscale expansion shows that the pressure gradient is the result of an advection procedure. It is therefore possible to remove the temperature coupling from the collision step and to introduce it during advection as follows

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t)-g_i(x，t)](3.1)

Wherein the equilibrium distribution is the same as equation (1.2) and the pressure coupling term is given by

In this example, where ρ is the fluid density; t is₀Is the constant lattice temperature; p is the pressure in the fluid volume; t is the temperature of the fluid during transport; and w_iIs a constant weighting factor. In equation 3.1 above, the first part is the mathematical expression for the crash process (i.e., f)_i(x,t)+C_i(x, t)), wherein C_iIs a collision operator. The second term of this equation is the thermodynamic step (e.g., g)_i(x,t+Δt)-g_i(x, t)). By performing a multiscale analysis, it can be shown that the second moment of the coupling term results in the correct pressure gradient without affecting the other quantities in the mass and momentum conservation equations.

Another method of temperature coupling is to add a time variation to the coupling as well

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t+Δt)-g_i(x，t)](3.3)

The difference between equation [3.1] and equation [3.3] is the new force term g calculated at t + Δ t. This introduces an additional time derivative term to the conservation of mass. This can be eliminated simply by changing the stop state as follows

These techniques for introducing temperature during advection require minimal changes to the conventional LBM method. The algorithm for temperature coupling versus time, equation [3.3], is described below with respect to FIG. 12.

Referring to fig. 12, the system implements a process 1200 in generating a distribution function for particle transport, where the thermodynamic step is separated from the collision step. In this example, the distribution function includes an advection step, and the thermodynamic step is included in the distribution portion by enhancing the advection step with the thermodynamic step rather than enhancing the particle collision step with the thermodynamic step. Generally, advection involves the transport of particles (e.g., from one region to another in a horizontal direction). In this example, the thermodynamic step includes the temperature of the fluid volume during transport.

In operation, the system simulates (1202) the transport of particles in a fluid volume in a set of lattice velocities, where the transport causes collisions between the particles. The system also generates (1203) a distribution function for the transport of the particles, wherein the distribution function comprises a thermodynamic step and a particle collision step, and wherein the thermodynamic step is substantially independent of and separate from the particle collision step. The generated process includes acts 1204, 1206, 1208, 1210, 1212, 1214, 1216, 1218.

In the example of FIG. 12, the system determines (1204) a post-collision distribution function f for determining a collision at a particular location x in the fluid volume at a particular time t_i' (x, t), wherein f_i′(x,t)＝f_i(x,t)+C_i(x, t) wherein C_iIs a collision operator, and f_iIs a distribution function for the particles before collision.

In this example, the system distributes the function f from after-collision_i' (x, t) subtracting (1206) fractional block g_i(x, t) to obtain a pre-advection particle density distribution function f_i″(x,t)＝f_i′(x,t)-g_i(x, t) in which g_i(x, t) represents the portion of the particle that must be retained at the same location during advection (e.g., the distribution function for the portion of the particle that is not advected during advection). In this example, the system separates the particles into a first part f that is advected as described below_i(x,t)-g_i(x, t) and a second portion g not advected during the advection step_i(x,t)。

The system simulates 1208 advection of a second portion of the particles at time t + Δ t to another location (e.g., an adjacent cell) in the fluid volume, denoted as (x + c)_iΔ t), wherein c_iIs the velocity vector of the particle before the collision and at is the interval between a certain time t and another point in time. In the example, the particles are advected to adjacent cells in the direction of the discrete particle velocity. The advected portion of the particle is a second portion of the particle (e.g., f)_i′(x,t)-g_i(x, t)). The system obtains (1210) a density distribution function of advected particles based on an advection simulationWhereinIs at position x from position x-c as described below_iIs advected. In this example, f_i″(x-c_iΔ t, t) denotes the position x-c at time t_iParticle density distribution function before advection at Δ t. Item x-c_iAt represents another location (e.g., in an adjacent cell) from which the particle is advected to location x. The other position x-c_iAt is based on the original position x and the distance the particle can move within the time interval at (based on the velocity value c)_i). Equation ofRepresents the distribution of particles at position x at time t + Δ t and at position x-c at time t_iThe advection-front particle density distribution function at Δ t is the same because at position x-c_iThe particle at Δ t is advected to position x.

Block g that the system will deduct before_i(x, t) adding (1212) back to the density distribution function (advected from neighboring cells)To form a density distribution function after advection:

f_i(x，t+Δt)＝f_i(x，t+Δt)+g_i(x，t)＝f_i(x-c_iΔt，t)+g_i(x, t) which is also denoted as:

in this example, the previously subtracted block g_i(x, t) addition back to the density distribution function results in conservation of mass. Moreover, the addition restores equation 3.1, as it can be accomplished by adding c_iAdded to the position vector x of all the terms in the above equation.

The system calculates (1214) the mass, momentum, and temperature of the particle at position x at time t + Δ t. The system determines (1216) g using the calculated temperature, mass, and momentum_i(x, t + Δ t). System will be poor g_i(x,t+Δt)-g_i(x, t) adding (1218) to each mobility state f_i(x, t + Δ t) to illustrate the time variation of the temperature coupling. This action recovers equation [3.3] as follows]

f(x，t＋Δt)＝f(X-c_iΔt，t)＋C_i(X-c_iΔt，t)+[g(x，t)-g_i(x-c_iΔt，t)]+[g_i(x，t＋Δt)-g_i(x，t)]

＝f(x-c_iΔt，t)+C_i(x-c_iΔt，t)+[g(x，t+Δt)-g_i(x-c_iΔt，t)]

The number of particles added to the move described in the above state is subtracted from the stopped state to conserve mass according to equation (3.4). The sum of the g terms in equation (3.4) represents the difference in population between the g groups of particles starting from the non-stopped state at time t and at time t + Δ t. The system introduces a summation term in the equation for the stopped state to supplement the difference, thereby ensuring overall mass conservation. The system repeats acts 1204, 1206, 1208, 1210, 1212, 1214, 1216, 1218. The whole cycle repeats from one time step t to a subsequent time step t + Δ t and so on.

In the example, equations 3.1 and 3.3 are at x + c_iHas a position of, and a speed of c_iIs advected from position x, rather than having a velocity c at position x_iFrom position x-c_iAnd (4) smoothing flow. To obtain the formula for 3.1, except that g is at position x + c_iAt is outside of the time t calculation (at act 1216), the system repeats the acts included in process 1200.

Utilizing the techniques described herein, the distribution is balancedIn a distribution function (e.g., f)_i(x+c_iΔ t, t + Δ t)) and thus the stability range remains the same as that of isothermal LBM. This enables simulation of high Mach number and high temperature range applicationsCan be used. Furthermore, the pressure gradient calculated as a result of advection involves all lattice directions, which preserves the isotropic nature of the LBM solver. Conservation of local momentum is ensured because no force is added that is approximately calculated with a finite difference. The velocity definition is exact, unlike force methods that require alternating velocity definitions to remove the discrete lattice effect. The technique described herein is computationally inexpensive because it does not involve extensive pressure gradient calculations and complex boundary processing. These techniques also have simpler boundary conditions that ensure accurate pressure gradients in the near-wall region, resulting in a correct velocity profile. This is crucial for the prediction of heat flux for high temperature applications. These techniques can be readily extended to any LBM application involving a conservative force, such as magnetic forces, gravity, inter-component forces in multiphase applications, and the like.

Fig. 13 is a block diagram of components of a network environment 1300. The network environment 1300 also includes a system 1302 (for implementing the techniques described herein), the system 1302 including a memory 1304, a bus system 1306, and a processor 1308. The memory 1304 may include a hard disk drive and a random access memory storage device, such as a dynamic random access memory, a machine readable hardware storage device, a machine readable medium, or other type of non-transitory machine readable storage device. A bus system 1306, including, for example, a data bus and a motherboard, may be used to establish and control data communication between components of system 1302. Processor 1308 may include one or more microprocessors and/or processing devices. In general, processor 1308 can include any suitable processor and/or logic capable of receiving and storing data and capable of communicating via a network (not shown).

The system 1302 may be any of a variety of computing devices capable of receiving data, such as a server, a distributed computing system, a desktop computer, a laptop computer, a cell phone, a rack-mounted server, and so forth. The system 1302 may be a single server or a group of servers at the same location or at different locations. The illustrated system 1302 can receive data via an input/output ("I/O") interface 1310. I/O interface 1310 may be any type of interface capable of receiving data via a network, such as an ethernet interface, a wireless networking interface, a fiber optic networking interface, a modem, and so forth. The system 1302 is configured for communication with a data store 1312, which data store 1312 may be configured to store speed models, simulation data, and the like.

Embodiments may be implemented in digital electronic circuitry, or in computer hardware, firmware, software, or in combinations of them. Apparatus of the techniques described herein may be implemented in a computer program product tangibly embodied or stored in a machine-readable medium (e.g., a hardware storage device) for execution by a programmable processor; and method acts may be performed by a programmable processor executing a program of instructions to perform operations of the techniques described herein by performing operations on input data and generating output. The techniques described herein may be implemented in one or more computer programs that are executable on a programmable system including at least one programmable processor coupled to receive data and instructions from, and to transmit data and instructions to, a data storage system, at least one input device, and at least one output device. Each computer program may be implemented in a high level procedural or object oriented programming language, or in assembly or machine language if desired; and in any case, the language may be a compiled or interpreted language.

Suitable processors include, by way of example, both general and special purpose microprocessors. Generally, a processor will receive instructions and data from a read-only memory and/or a random access memory. Generally, a computer will include one or more mass storage devices for storing data files; such devices include magnetic disks, such as internal hard disks and removable disks; magneto-optical disks; and an optical disc. Storage devices suitable for tangibly embodying computer program instructions and data include all forms of non-volatile memory, including by way of example semiconductor memory devices, such as EPROM, EEPROM, and flash memory devices; magnetic disks, such as internal hard disks and removable disks; magneto-optical disks; and a CD-ROM disk. Any of the foregoing may be supplemented by, or incorporated in, ASICs (application-specific integrated circuits).

Many implementations have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the claims. Accordingly, other implementations are within the scope of the following claims.

Claims (33)

1. A method, comprising:

simulating transport of particles in a fluid volume in a set of lattice velocities, wherein the transport causes collisions between the particles; and

generating a distribution function for the transport of the particles, wherein the distribution function comprises a thermodynamic step and a particle collision step, and wherein the thermodynamic step is substantially independent of and separate from the particle collision step.

2. The method of claim 1, wherein the distribution function further comprises an advection step, and wherein the thermodynamic step is included in the distribution portion by enhancing the advection step with the thermodynamic step instead of enhancing the particle collision step with the thermodynamic step.

3. The method of claim 1, wherein the thermodynamic step comprises a temperature of the fluid volume during the transporting.

4. The method of claim 1, wherein generating comprises:

determining a post-collision distribution function f 'for a collision at a particular location x in the fluid volume at a particular time t'_i(x, t), wherein f'_i(x,t)＝f_i(x,t)+C_i(x, t) wherein C_iIs a collision operator, and f_iIs a distribution function for the particles prior to collision;

from post-collision distribution function f'_i(x, t) deducting a fraction block g_i(x, t) to obtain a pre-advection particle density distribution function f_i″(x,t)＝f_i′(x,t)-g_i(x, t), wherein a portion f of the particles_i′(x,t)-g_i(x, t) is advected to another location in the fluid volume, where g_i(x, t) represents the distribution of particles that are not advected;

simulating advection of the portion of the particle at time t + Δ t to the another location in the fluid volume, the another location being represented as (x + c)_iΔ t), wherein c_iIs the velocity vector of the particle before the collision, and Δ t is the interval between a certain time t and another point in time;

obtaining a density distribution function of the advected particles based on a simulation of the advectionWhereinAnd whereinIs at position x from position x-c_iDistribution of advected particles;

block g to be deducted before_i(x, t) add back to density distribution functionTo form a density distribution function after advection

Calculating the mass, momentum and temperature of the particle at the location x at time t + Δ t;

determining g using the calculated temperature, mass and momentum_i(x, t + Δ t); and

the difference g_i(x,t+Δt)-g_i(x, t) addition to moving State f_i(x,t)+C_i(x,t)。

5. The method of claim 4, wherein g_iIs defined according to the following formula:

where ρ is the fluid density;

wherein T is₀Is the constant lattice temperature;

wherein P is the pressure in the fluid volume;

wherein T is the calculated temperature; and

wherein w_iIs a constant weighting factor.

6. The method of claim 4, wherein the generated distribution function is according to the following equation:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t+Δt)-g_i(x，t)]；

wherein x is a particular location within the volume;

wherein t is a specific first point in time;

wherein i is the index number of the lattice velocities in the set;

wherein c is_iIs the velocity vector of the particle before the collision;

wherein C is_iIs a collision operator;

wherein Δ t is the interval between the first and second points in time;

wherein g is_iIs a thermodynamic step; and

wherein f is_iIs a distribution function for the particles at the location x at time t.

7. The method of claim 4, further comprising:

the mass is conserved by modifying the stopped state according to:

8. the method of claim 1, wherein the particle collision step comprises an isothermal equilibrium distribution function.

9. The method of claim 1, wherein the distribution function is for a location x + c_iThe distribution function of the particles at Δ t and at time t + Δ t and is expressed as:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t)-g_i(x，t)]；

wherein x is a particular location within the volume;

wherein t is a specific first point in time;

wherein i is the index number of the lattice velocities in the set;

wherein C is_iIs a collision operator;

wherein c is_iIs the velocity vector of the particle before the collision;

wherein Δ t is the interval between the first point in time and the second point in time;

wherein g is_iIs a thermodynamic step; and

wherein f is_iIs a distribution function for the particles at the location x at time t.

10. The method of claim 9, wherein g_iIs defined according to the following formula:

where ρ is the fluid density;

wherein T is₀Is the constant lattice temperature;

wherein P is the pressure in the fluid volume;

wherein T is the actual temperature of the fluid; and wherein w_iIs a constant weighting factor.

11. The method of claim 1, wherein the set of lattice velocities is based on the lattice boltzmann method.

12. One or more machine-readable hardware storage devices storing instructions, wherein the instructions are executable by one or more processing devices, the instructions comprising:

simulating transport of particles in a fluid volume in a set of lattice velocities, wherein the transport causes collisions between the particles; and

generating a distribution function for the transport of the particles, wherein the distribution function comprises a thermodynamic step and a particle collision step, and wherein the thermodynamic step is substantially independent of and separate from the particle collision step.

13. The one or more machine readable hardware storage devices of claim 12, wherein the distribution function further comprises an advection step, and wherein the thermodynamic step is included in the distribution portion by enhancing the advection step with the thermodynamic step instead of enhancing the particle collision step with the thermodynamic step.

14. The one or more machine-readable hardware storage devices of claim 12, wherein the thermodynamic step comprises a temperature of the fluid volume during the transporting.

15. The one or more machine-readable hardware storage devices of claim 12, wherein generating comprises:

determining a post-collision distribution function f 'for a collision at a particular location x in the fluid volume at a particular time t'_i(x, t), wherein f'_i(x,t)＝f_i(x,t)+C_i(x, t) wherein C_iIs a collision operator, and f_iIs a distribution function for the particles prior to collision;

from post-collision distribution function f'_i(x, t) deducting a fraction block g_i(x, t) to obtain a pre-advection particle density distribution function f_i″(x,t)＝f_i′(x,t)-g_i(x, t), wherein a portion f of the particles_i′(x,t)-g_i(x, t) is advected to another location in the fluid volume, where g_i(x, t) represents the distribution of particles that are not advected;

simulating advection of the portion of the particle at time t + Δ t to the another location in the fluid volume, the another location being represented as (x + c)_iΔ t), wherein c_iIs the velocity vector of the particle before the collision, and Δ t is the specific time t andan interval between another point in time;

obtaining a density distribution function of the advected particles based on a simulation of the advectionWhereinAnd whereinIs at position x from position x-c_iDistribution of advected particles;

block g to be deducted before_i(x, t) add back to density distribution functionTo form a density distribution function after advection

Calculating the mass, momentum and temperature of the particle at the location x at time t + Δ t;

determining g using the calculated temperature, mass and momentum_i(x, t + Δ t); and

the difference g_i(x,t+Δt)-g_i(x, t) addition to moving State f_i(x,t)+C_i(x,t)。

16. The one or more machine-readable hardware storage devices of claim 15, wherein g is_iIs defined according to the following formula:

where ρ is the fluid density;

wherein T is₀Is the constant lattice temperature;

wherein P is the pressure in the fluid volume;

wherein T is the calculated temperature; and

wherein w_iIs a constant weighting factor.

17. The one or more machine readable hardware storage devices of claim 15, wherein the generated distribution function is according to:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t+Δt)-g_i(x，t)]；

wherein x is a particular location within the volume;

wherein t is a specific first point in time;

wherein i is the index number of the lattice velocities in the set;

wherein c is_iIs the velocity vector of the particle before the collision;

wherein C is_iIs a collision operator;

wherein Δ t is the interval between the first and second points in time;

wherein g is_iIs a thermodynamic step; and

wherein f is_iIs a distribution function for the particles at the location x at time t.

18. The one or more machine-readable hardware storage devices of claim 15, wherein the operations further comprise:

the mass is conserved by modifying the stopped state according to:

19. the one or more machine readable hardware storage devices of claim 12, wherein the particle collision step comprises an isothermal equilibrium distribution function.

20. The one or more machine-readable hardware storage devices of claim 12, wherein the distribution function is for at location x + c_iThe distribution function of the particles at Δ t and at time t + Δ t and is expressed as:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t)-g_i(x，t)]；

wherein x is a particular location within the volume;

wherein t is a specific first point in time;

wherein i is the index number of the lattice velocities in the set;

wherein C is_iIs a collision operator;

wherein c is_iIs the velocity vector of the particle before the collision;

wherein Δ t is the interval between the first point in time and the second point in time;

wherein g is_iIs a thermodynamic step; and

wherein f is_iIs a distribution function for the particles at the location x at time t.

21. The one or more machine-readable hardware storage devices of claim 20, wherein g is_iIs defined according to the following formula:

where ρ is the fluid density;

wherein T is₀Is the constant lattice temperature;

wherein P is the pressure in the fluid volume;

wherein T is the actual temperature of the fluid; and wherein w_iIs a constant weighting factor.

22. The one or more machine readable hardware storage devices of claim 12, wherein the set of lattice velocities is based on a lattice boltzmann method.

23. A system, comprising:

one or more processing devices; and

one or more machine-readable hardware storage devices storing instructions, wherein the instructions are executable by one or more processing devices, the instructions comprising:

simulating transport of particles in a fluid volume in a set of lattice velocities, wherein the transport causes collisions between the particles; and

generating a distribution function for the transport of the particles, wherein the distribution function comprises a thermodynamic step and a particle collision step, and wherein the thermodynamic step is substantially independent of and separate from the particle collision step.

24. The system of claim 23, the distribution function further comprising an advection step, and wherein the thermodynamic step is included in the distribution portion by enhancing the advection step with the thermodynamic step instead of enhancing the particle collision step with the thermodynamic step.

25. The system of claim 23, wherein the thermodynamic step comprises a temperature of the fluid volume during the transporting.

26. The system of claim 23, wherein generating comprises:

determining a post-collision distribution function f 'for a collision at a particular location x in the fluid volume at a particular time t'_i(x, t), wherein f'_i(x,t)＝f_i(x,t)+C_i(x, t) wherein C_iIs a collision operator, and f_iIs a distribution function for the particles prior to collision;

from post-collision distribution function f'_i(x, t) deducting a fraction block g_i(x, t) to obtain a pre-advection particle density distribution function f_i″(x,t)＝f_i′(x,t)-g_i(x, t), wherein a portion f of the particles_i′(x,t)-g_i(x, t) is advected to another location in the fluid volume, where g_i(x, t) represents the distribution of particles that are not advected;

simulating advection of the portion of the particle at time t + Δ t to the another location in the fluid volume, the another location being represented as (x + c)_iΔ t), wherein c_iIs the velocity vector of the particle before the collision, and Δ t is the interval between a certain time t and another point in time;

obtaining a density distribution function of the advected particles based on a simulation of the advectionWhereinAnd whereinIs at position x from position x-c_iDistribution of advected particles;

block g to be deducted before_i(x, t) add back to density distribution functionTo form a density distribution function after advection

Calculating the mass, momentum and temperature of the particle at the location x at time t + Δ t;

determining g using the calculated temperature, mass and momentum_i(x, t + Δ t); and

the difference g_i(x,t+Δt)-g_i(x, t) addition to moving State f_i(x,t)+C_i(x,t)。

27. The system of claim 26, wherein g_iIs defined according to the following formula:

where ρ is the fluid density;

wherein T is₀Is the constant lattice temperature;

wherein P is the pressure in the fluid volume;

wherein T is the calculated temperature; and

wherein w_iIs a constant weighting factor.

28. The system of claim 26, wherein the generated distribution function is according to the following equation:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t+Δt)-g_i(x，t)]；

wherein x is a particular location within the volume;

wherein t is a specific first point in time;

wherein i is the index number of the lattice velocities in the set;

wherein c is_iIs the velocity vector of the particle before the collision;

wherein C is_iIs collision operator；

Wherein Δ t is the interval between the first and second points in time;

wherein g is_iIs a thermodynamic step; and

wherein f is_iIs a distribution function for the particles at the location x at time t.

29. The system of claim 26, wherein the operations further comprise:

the mass is conserved by modifying the stopped state according to:

30. the system of claim 23, wherein the particle collision step comprises an isothermal equilibrium distribution function.

31. The system of claim 23, wherein the distribution function is for a location x + c_iThe distribution function of the particles at Δ t and at time t + Δ t and is expressed as:

f_i(x+c_iΔt，t+Δt)＝f_i(x，t)+C_i(x，t)+[g_i(x+c_iΔt，t)-g_i(x，t)]；

wherein x is a particular location within the volume;

wherein t is a specific first point in time;

wherein i is the index number of the lattice velocities in the set;

wherein C is_iIs a collision operator;

wherein c is_iIs the velocity vector of the particle before the collision;

wherein Δ t is the interval between the first point in time and the second point in time;

wherein g is_iIs a thermodynamic step; and

wherein f is_iIs a distribution function for the particles at the location x at time t.

32. The system of claim 23, wherein g_iIs defined according to the following formula:

where ρ is the fluid density;

wherein T is₀Is the constant lattice temperature;

wherein P is the pressure in the fluid volume;

wherein T is the actual temperature of the fluid; and wherein w_iIs a constant weighting factor.

33. The system of claim 32, wherein the set of lattice velocities is based on a lattice boltzmann method.