JP7496049B2 - A lattice Boltzmann solver that enforces total energy conservation - Google Patents

Description

[Claim of priority]
This application claims priority under 35 U.S.C. §119(e) to U.S. Provisional Patent Application No. 62/790,528, filed January 10, 2019, and entitled “LATTICE BOLTZMANN SOLVER ENFORCING TOTAL ENERGY CONSERVATION,” the contents of which are incorporated herein by reference in their entirety.

This specification relates to computer simulation of physical processes, such as physical fluid flows.

The so-called "Lattice Boltzmann Method" (LBM) is an advantageous technique used in computational fluid dynamics. The underlying dynamics of a lattice Boltzmann system resides in the fundamental physics of kinetic gases, with the motion of many particles following the Boltzmann equation. In a Boltzmann fundamental kinetic system, there are two fundamental dynamical processes: collision and advection. The collision process involves the interaction between particles following conservation laws, and relaxation to an equilibrium state. The advection process involves modeling the motion of particles from one location to another following microscopic particle velocities.

In compressible flows, a single distribution function is typically used to satisfy the conservation of mass, momentum and total energy. If a single distribution function is used to solve for all three quantities, the grid used for LBM should satisfy up to the eighth order of moments. In this situation, numerous grid velocities are required as the moment requirements increase and the resulting stencil size also increases. Other considerations with the use of a single particle distribution function include considerations regarding the stability of the distribution function, momentum and energy boundary conditions, body forces and heat sources, etc.

US Patent Application Publication No. 2016-0188768 U.S. Pat. No. 9,576,087

According to one aspect, a method of simulating fluid flow on a computer using a simulation that enforces total energy conservation includes simulating fluid activity across a mesh to model particle movement across the mesh; storing in a computer accessible memory a set of state vectors for each mesh location within the mesh, the set of state vectors including a plurality of entries, each entry corresponding to a particular one of the possible momentum states at the corresponding mesh location; the computer calculating a set of energy values for the mesh locations within the mesh; the computer performing advection of the particles to subsequent mesh locations over a time interval; and the computer modifying the state vectors of the particles by adding specific total energy values to the advected particle states and subtracting specific total energy values from the particle states that were not advected over the time interval.

Aspects also include computer program products, devices, and computer systems on one or more machine-readable hardware storage devices.

One or more of the above aspects may include one or more of the following features, among others, as described herein:

After the particles are advected to the subsequent mesh location according to the modified set of states, the aspect includes adding back a local pressure term to the particle state and then calculating the moment to result in a correct pressure gradient (∇p). The local pressure term includes a θ-1 term that is calculated by the computer. The step of simulating the fluid flow activity includes simulating the fluid flow based in part on the first set of discrete grid velocities, and the method further includes simulating the time evolution of the scalar quantity based in part on the second set of discrete grid velocities. The second set of discrete grid velocities is the same grid velocity as the first set of discrete grid velocities. The second set of discrete grid velocities has a number of grid velocities different from the number of discrete grid velocities in the first set. The step of simulating the time evolution of the total energy includes collecting inflow distributions from neighboring mesh locations for a collision operator and an energy operator, weighting the inflow distributions, determining an outflow distribution as a product of the collision operator and the energy operator, and propagating the determined outflow distribution.

ここでの

は、非平衡成分のフィルタ処理された２次モーメントである。 Simulating the time evolution of the scalar quantity includes simulating the time evolution of the scalar quantity based in part on a collision operator that effectively represents the momentum flux resulting from the product of a collision operator and an energy operator as a first order term. To recover the exact shear stress for any Prandtl number, the energy collision term is given by:

Here

is the filtered second moment of the unbalanced component.

The aspect further comprises applying a zero net surface flux boundary condition such that the inflow distribution is equal to the determined outflow distribution. Determining the outflow distribution comprises determining the outflow distribution to result in a zero surface scalar flux. The scalar quantity comprises a scalar quantity selected from the group consisting of temperature, concentration and density. Calculating a set of energy values for mesh locations within the mesh further comprises calculating a total energy _Et given by _Et = E + ^v2 /2, where E is energy and v is velocity.

ここでの項Ω_f、Ω_qは、それぞれの衝突演算子を表し、上記方程式において使用される平衡分布は以下の通りである、

The step of simulating the activity includes the step of the computer applying a second distribution function of the equilibrium distribution and the specific total energy _Et , where the second distribution is defined as a particular scalar that is advected with the flow distribution _fj . The second distribution function takes into account the non-equilibrium contribution of _fj to the energy equation to obtain accurate flow behavior over different grid resolutions near the boundary, where the collision operator of the distribution function is given by:

where the terms Ω _f , Ω _q represent the respective collision operators and the equilibrium distribution used in the above equation is:

を

となるように修正することによって行われる。 The total energy is calculated by adding a term before advection as pressure is convected, and then removing energy from the resting state to compensate for the additional energy. By removing the additional energy, the total energy is conserved, resulting in the correct pressure velocity term. Removing the additional energy is equivalent to removing the resting state

of

This is done by modifying it so that

In the method disclosed herein, solving the total energy uses a second distribution defined as a specific scalar that is advected with the flow distribution. The method for simulating fluid flow described herein uses the Lattice Boltzmann (LB) method and a scalar solver transport equation. An LBM solver for compressible flows with total energy conservation and high stability zone is provided, which can be used in practical applications. During collisions, the equilibrium distribution used in the solver is positive, and during advection, temperature coupling is performed.

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

FIG. 1 illustrates a fluid flow simulation system. 1 is a flowchart showing operations for complex fluid flow simulation and total energy calculation. 1 is a flowchart showing operations for complex fluid flow simulation and total energy calculation. FIG. 1 illustrates an exemplary computation-based comparison of expected and predicted outcomes. FIG. 1 illustrates an exemplary computation-based comparison of expected and predicted outcomes. FIG. 1 illustrates an exemplary computation-based comparison of expected and predicted outcomes. FIG. 2 illustrates velocity components of two LBM models (prior art). FIG. 2 illustrates velocity components of two LBM models (prior art). 1 is a flowchart of a procedure followed by a physical process simulation system (prior art). FIG. 2 is a perspective view of a microblock. FIG. 2 is a diagram of a lattice structure (prior art) used by the system of FIG. 1; FIG. 2 is a diagram of a lattice structure (prior art) used by the system of FIG. 1; FIG. 1 illustrates variable resolution technology (prior art). FIG. 1 illustrates variable resolution technology (prior art). FIG. 1 illustrates a parallelepiped that is useful in understanding particle-facet interactions. FIG. 1 illustrates the movement of particles from a voxel to a surface (prior art). FIG. 1 illustrates the movement of particles from surface to surface (prior art). 1 is a flow chart of a procedure for performing surface dynamics (prior art); 1 is a flow chart of a process for generating a distribution function using a thermodynamic step that is independent of the particle collision step.

The steps illustrated in Figure 9 below describe a flow simulation process using a lattice Boltzmann solver that enforces total energy conservation.

In Figures 7, 8, and 10-18, the figures are labeled "Prior Art." These figures are labeled as prior art because they are commonly found in the patents referenced above. However, these figures are as they appear in the patent applications referenced above. However, the patent applications referenced above do not describe a lattice Boltzmann solver that implements total energy conservation, as described below, and therefore these figures do not take into account any modifications that may be made to the flow simulation using this LB solver.

A. Methods for Solving for Scalar Quantities In achieving complex fluid flow simulations, it can be beneficial to simultaneously solve for scalar quantities, such as temperature distribution, concentration distribution and/or density, along with solving for the fluid flow. The systems and methods described herein couple the modeling of scalar quantities (as opposed to vector quantities) to the modeling of fluid flow based on an LBM-based physical process simulation system. Exemplary scalar quantities that can be modeled include temperature, concentration and density.

Specific details of the technique for coupling scalar modeling to simulated fluid flow are detailed in co-pending U.S. Patent Application Publication No. 2016-0188768, entitled "Temperature Coupling Algorithm for Hybrid Thermal Lattice Boltzmann Method," by Pradeep Gopalakrishman et al., the contents of which are incorporated herein by reference in their entirety. A collision operator based on a Galilean-invariant collision operator with a relatively high stability region is used. Specific details of this collision operator are described in U.S. Patent No. 9,576,087, the contents of which are incorporated herein by reference in their entirety.

For example, the system can be used to determine the convection temperature distribution within a system. For example, if a system (formed by a volume represented by multiple voxels) contains a heat source and there is an airflow within the system, some regions of the system will be hotter than other regions based on their proximity to the airflow and the heat source. To model such a situation, the temperature distribution within the system can be represented as a scalar quantity, with each voxel having an associated temperature.

In another example, the system can be used to determine the distribution of chemicals within a system. For example, if a system (formed of a volume represented by multiple voxels) contains a contamination source, such as a dirty bomb or chemicals or other particles suspended in the air or liquid, and there is an air or liquid flow within the system, some regions of the system will have higher concentrations than other regions based on the flow and proximity to the contamination source. To model such a situation, the distribution of chemicals within the system can be represented as a scalar quantity, with each voxel having an associated concentration.

In some applications, multiple different scalar quantities can be simulated simultaneously. For example, the system can simulate both the temperature distribution and the concentration distribution within the system.

Scalar quantities can be modeled in various ways. For example, the lattice Boltzmann (LB) method for solving the scalar transport equation can be used to indirectly solve the scalar transport. For example, the method described herein can provide an indirect solution of the second-order macroscopic scalar transport equation:

In addition to the lattice Boltzmann function for fluid flow, a second set of distribution functions is introduced for the transport scalars. The method assigns a vector to each voxel in the volume to represent the fluid flow, and a scalar to each voxel in the volume to represent the desired scalar variable (e.g., temperature, density, concentration, etc.). The method fully recovers the macroscopic scalar transport equations that satisfy the exact conservation laws. This method is believed to improve the accuracy of the calculated scalars compared to other non-LBM methods. The method is also believed to improve the ability to consider complex boundary shapes.

This method of modeling scalar quantities can be used with time-explicit CFD/CAA solvers based on the Lattice Boltzmann Method (LBM), such as the PowerFLOW system available from Exa Corporation, Burlington, Massachusetts. Unlike methods based on discretizing the macroscopic continuum equations, LBM predicts macroscopic fluid dynamics starting from the "mesoscopic" Boltzmann equations of motion. The resulting compressible and unsteady solution method can be used to predict the physics of a variety of complex flows, including aeroacoustic and pure acoustic problems. Below, we provide a general discussion of LBM-based simulation systems, followed by a description of scalar solution methods that can be used with fluid flow simulations to support such modeling methods.

方程式（Ｉ１）
式中、

である。 Modeling the Simulation Space In an LBM-based physical process simulation system, a fluid flow is represented by distribution function values f _i evaluated at a set of discrete velocities c _i . The dynamics of the distribution function follows equation I1, where f _i (0) is known as the equilibrium distribution function defined as:

Equation (I1)
In the formula,

It is.

This equation is the well-known lattice Boltzmann equation that describes the time evolution of the distribution function f _i . The left-hand side describes the change in distribution due to the so-called "flow process". The flow process is the point in time when a fluid pocket starts at one mesh location and then moves along one of several velocity vectors to the next mesh location. At this point, we calculate the "collision coefficient", i.e. the influence of neighboring fluid pockets on the starting fluid pocket. Fluids can only move to another mesh location, so the correct choice of velocity vectors is necessary so that all components of all velocities are multiples of a common velocity.

The right hand side of the first equation is the "collision operator" mentioned above, which describes the change in the distribution function due to collisions between fluid pockets. A particular form of the collision operator can be the Bhatnagar, Gross and Krock (BGK) operator. The collision operator drives the distribution function to a prescribed value given by the second equation, which is its "equilibrium" form.

方程式(Ｉ２)
ここでのパラメータτは、衝突を介した平衡までの特性緩和時間を表す。粒子（例えば、原子又は分子）の取り扱いでは、通常は緩和時間が定数とみなされる。 The BGK operator is constructed according to the physical description that, whatever the details of the collision, the distribution function approaches a well-defined local equilibrium given by {f _i ^eq (x,v,t)}{f ^eq (x,v,t)} through the collision,

Equation (I2)
The parameter τ here represents the characteristic relaxation time to equilibrium via collisions. In dealing with particles (e.g. atoms or molecules), the relaxation time is usually taken to be a constant.

方程式(Ｉ３)
ここでのρ、ｕ及びＴは、それぞれ流体の密度、速度及び温度であり、Ｄは、（必ずしも物理的空間次元とは等しくない）離散化速度空間の次元である。 From this simulation, traditional fluid variables such as mass ρ and fluid velocity u are obtained as simple sums in equation (I3) below:

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

Due to symmetry considerations, the set of velocity values is chosen to form a particular lattice structure when spanning the configuration space. The dynamics of such a discrete system obeys the LBE, which has the form
f _i (x+c _i ,t+1)-f _i (x,t)=c _i (x,t)
The collision operator here usually takes the form BGK as described above. With the correct choice of equilibrium distribution form, it can be theoretically shown that the lattice Boltzmann equations yield the correct hydrodynamic and thermohydrodynamic results. That is, the hydrodynamic moments derived from f _i (x,t) obey the Navier-Stokes equations in the macroscopic limit. These moments are defined by equation (I3) above.

The collective values of c _i and w _i define the LBM model, which can be efficiently implemented on scalable computing platforms and run with high robustness for time unsteady flows and complex boundary conditions.

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 a successive approximation of the full Boltzmann equation. In a fluid system, small disturbances in density travel at the speed of sound. In gas systems, the speed of sound generally depends on the temperature. The importance of the effect of compressibility of the flow is measured by the ratio of the characteristic velocity to the speed of sound, which is known as the Mach number. For a further description of conventional LBM-based physical process simulation systems, see U.S. Patent Application Publication No. 2016/0188768, incorporated by reference above.

1 shows a system 10 for simulating fluid flow, such as around a representation of a physical object. In this implementation, the system 10 is based on a client-server architecture and includes a server system 12 implemented as a massively parallel computer system 12 and a client system 14. 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 unit 24. Within the memory 18 resides a mesh preparation engine 32 and a simulation engine 34. The system 10 has access to a data repository 38 that stores 2D and/or 3D meshes, coordinate systems, and libraries.

1 shows the mesh preparation engine 32 in memory 18, the mesh preparation engine can be a third party application running on a different system than the server 12. Whether the mesh preparation engine 32 runs in memory 18 or on a different system than the server 12, it receives the user-specified mesh definition 30 and prepares and transmits the mesh to the simulation engine 34. As disclosed below, the simulation engine 34 includes a particle collision interaction module, a particle boundary model module, and an advection module that performs advection operations. The simulation engine 34 also includes a total energy solver as described below in FIG. 3.

への粒子の移流とのシミュレーションを含む。 2, a process 40 for simulating fluid flow around a representation of a physical object is shown. 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, specifically having planar and/or curved surface(s). The process receives (42), for example, a mesh of the physical object to be simulated, for example, by retrieving it from the client system 14 or from a data repository 42. In other embodiments, either an external system or the server 12 generates the mesh of the physical object to be simulated based on user input. The process pre-computes (44) geometric quantities from the retrieved mesh, and runs (46) a dynamic lattice Boltzmann model simulation using the pre-computed geometric quantities corresponding to the retrieved mesh. The lattice Boltzmann model simulation includes the evolution of the particle distribution and the next cell in the LBM mesh according to a set of state vectors that are modified by adding specific total energy values to the states of the advected particles and subtracting specific total energy values from the states of the non-advected particles over a period of time.

Simulations with particle advection into the void are included.

Total Energy Conserving Compressible Flow Solver In achieving complex fluid flow simulations, it can be beneficial to simultaneously solve for scalar quantities such as temperature distribution, concentration distribution and/or density along with solving the fluid flow. The distinction between a particular scalar and a scalar distribution is made by using a particular scalar q, which is the difference between q (Eq. 10) and h (Eq. 2), E_t, whereas h is density*E_t. The method described is such that the scalar exists on top of the advecting flow, as described in the above-mentioned U.S. Patent Application Publication No. 2016-0188768.

方程式１
ここでの全エネルギーＥ_tは、Ｅ_t＝Ｅ＋ν²／２によって与えられ、ρは密度であり、νは速度であり、ｐは圧力であり、ｋは熱伝導率であり、Ｔは温度であり、τは剪断応力である。この方程式は、高速流のための全計算流体力学ツールによって解く必要がある一般方程式である。 3, a compressible LBM solver using a total energy conservation process 50 that has a better stability zone than conventional methods and can be applied in practical applications is described. The total energy equation that needs to be solved is given by:

Equation 1
The total energy _Et here is given by _Et = E + ^v2 /2, where p is density, v is velocity, p is pressure, k is thermal conductivity, T is temperature, and τ is shear stress. This equation is the general equation that needs to be solved by all computational fluid dynamics tools for high speed flows.

One of the complications with the total energy solver is the addition of the pressure convection term "p" in Equation 1. The Lattice Boltzmann Method (LBM), as well as all computational fluid dynamics tools for high-speed flows, must solve this equation for total energy conservation. The LBM method has some shortcomings in solving this equation for all cases, especially for the pressure convection term p. The disclosed method can solve this equation for complex problems and has other advantages.

上記方程式の左辺は対流項を含み、右辺は拡散項を有する。全エネルギー方程式とは異なり、上記方程式の全ての項は１つのスカラー変数Ｓしか含んでおらず、このことが移流を単純化する。プロセス５０は、選択された平衡分布形態に従って粒子分布の発展をシミュレートする（５２）。 The pressure convection problem as it exists for the LBM method can be described as follows: The normal scalar transport takes the form of the following equation:

The left side of the above equation contains a convective term and the right side has a diffusive term. Unlike the total energy equation, all terms in the above equation contain only one scalar variable S, which simplifies advection. The process 50 simulates (52) the evolution of the particle distribution according to a selected equilibrium distribution regime.

In addition to the equilibrium distribution (a selected form of Eq. 1) 52, a second distribution function is used to calculate the specific total energy _Et (54). The process 50 includes advection of particles to the next cell q in the LBM (56). The realization of energy conservation is achieved by the product of two distributions. The product of these distribution functions is defined as a particular scalar that is advected with the flow distribution _fj . This second distribution function differs from existing methods that use a second distribution function alone to solve for the energy, and takes into account the non-equilibrium contribution of _fj to the energy equation to obtain accurate flow behavior over different grid resolutions near the boundaries.

方程式２

方程式３ The collision operator for the distribution function is given by:

Equation 2

Equation 3

方程式４ The terms _Ωf and _Ωq in Equation 2 and Equation 3 are used to represent the respective collision operators. The collision operator _Ωf in Equation 5 is explained extensively in U.S. Patent No. 9,576,087, the contents of which are incorporated herein by reference in their entirety. The equilibrium distributions used in the above equations are as follows:

Equation 4

Distribution Function for the Total Energy Solver An obvious observation in the above formulation of Equation 7 is that the equilibrium distribution does not include the pressure term (p) as in Equation 1, and the θ-1 term, p=ρRT, which is essential for heat flow. The second distribution is the specific total energy, which is the sum of the internal energy E=C _v T and the kinetic energy v ² /2.

Therefore, this distribution function always results in a positive value during the collision calculation, and therefore this distribution function provides the solver with a corresponding high stability. If the post-collision state in equation (2) is advected as is, the distribution function in equation 4 will only result in an isothermal pressure gradient of momentum and will not result in a pressure convection term in the total energy equation. In contrast to other methods that use an equilibrium distribution and include pressure terms during advection and collision such as the following, while the pressure gradient term and pressure convection term arise due to advection, collisions do not affect these terms.

The term θ=RT/T ₀ makes the above equilibrium negative in response to pressure, and therefore the existing technique is highly unstable. Also, the use of the second distribution function h _i ^eq alone for all energies introduces noise across different grid resolutions.

The process 50 modifies the advected states by adding specific total energy to the advected states while removing specific total energy from the rest (non-moving) states (58).

方程式５ Specifically, in this modification (58) that recovers the accurate momentum and energy equations, the introduction of pressure occurs only during the advection step. Prior to the advection process, the moving state (∀i≠0) values at the mesh sites are modified as in Equation 8 below:

Equation 5

を

となるように修正することによって行われる。 However, adding this term (the second equation in Equation 5) _to qi also adds additional energy to the moving state. To compensate for the additional energy represented by this term, we remove the same amount of energy from the resting state (where no particles are advecting). By doing this, the process conserves energy, resulting in the correct pressure velocity term. Specifically, the conservation of total energy is

of

This is done by modifying it so that

方程式６ After advection, we add the local pressure term back to the particle state before computing the moments, which gives us the correct pressure gradient (∇p).

Equation 6

The pre-advection (equation 8) and post-advection (equation 6) modifications restore the correct equation for the momentum state (or the correct pressure gradient (∇p)) and also the correct pressure convection term ∇(up) in the total energy. Because the pressure term is eliminated (equation 9), the state becomes positive during collisions, improving the stability region of the total energy solver.

As discussed below in FIG. 9 , advection is corrected by applying a total energy correction for advection in which the computer modifies the state vectors of the particles by adding a specific total energy value to the states of the advected particles and subtracting a specific total energy value from the states of the particles that were not advected over that time interval.

方程式７
ここでの方程式７は、全エネルギーに対するフィルタ処理された非平衡寄与である。エネルギー保存は、方程式（７）の右辺における第２項

である衝突項を以下の状態にする２つの状態の乗算によって満たされる。

方程式８ Regularized Collision Operator for Energy The collision operator _Ωf used for the flow state is discussed in the above mentioned U.S. Pat. No. 9,576,087. In the following, we consider the collision operator _Ωq used for the energy. The regularized collision with first moment is given by

Equation 7
Equation 7 here is the filtered non-equilibrium contribution to the total energy. Energy conservation is the second term on the right-hand side of equation (7)

The collision term is satisfied by the multiplication of two states, which results in the following state:

Equation 8

方程式９ The zeroth moment of the above term is zero, conserving energy. During the Chapman-Enskog expansion, the first moment of the collision term is involved in the derivation of thermal diffusion and is driven by shear stresses. The collision term can be expressed as follows, as a first order approximation:

Equation 9

を乗算した結果、運動量流束項ρυυが生じる。この項ρυυは、高速流では数値的不安定性をもたらすとともに、マッハ数依存の拡散ももたらす。これらの影響を軽減するために、衝突演算子を以下のように定義することによって項ρυυを削除する。

方程式１０ In the above equation, both states

The multiplication by results in a momentum flux term ρυυ. This term ρυυ can lead to numerical instabilities at high speeds as well as Mach number dependent diffusion. To reduce these effects, we eliminate the term ρυυ by defining the collision operator as

Equation 10

The above equation eliminates the momentum flux term ρυυ, improving numerical stability in high Mach flows.

Many methods are available for the collision operator. One of these methods is called "regularization". In the (pending) patent by Hudong and Raoyang on lattice Boltzmann based scalars, "regularization" is used in conjunction with the additional Galilean invariance feature. However, this regularized collision operator is not ideal for use with high speed flows. During the derivation, the operator ends up having higher order terms such as ρυυ (Eq. 12). In LBM, the velocity values are generally always less than 1, so for low speed flows these error terms become much smaller and can be ignored. However, for high speed flows the velocity values become high, almost 1, so the error term ρυυ becomes relatively large and is eliminated by Eq. 10.

方程式１１
ここでの

は、非平衡成分のフィルタ処理された２次モーメントである。 Viscous effects for arbitrary Prandtl numbers. The relaxation time used in the energy solver is the shear stress with _viscosity equal to the thermal diffusivity (τ _q −0.5ρRT), which leads to erroneous results. To recover the correct shear stress for arbitrary Prandtl numbers, we add an additional term to the energy impingement term:

Equation 11
Here

is the filtered second moment of the unbalanced component.

Total Energy Solver Capabilities The total energy solver can be included in existing flow solvers, such as PowerFLOW®, and can be used for a wide range of industrial applications, including those involving subsonic and transonic flows.

In this section, we present several benchmarks to demonstrate some of the potential advantages of the disclosed total energy solver, especially for compressible flows related to energy conservation aspects. Simulation results are compared to analytical solutions and to a typical finite difference PDE-based hybrid solver solution.

Figure 4 shows a graph of the total temperature for the flow across the plate at x = 1.5 m from the leading edge at a Mach speed (Ma) of Ma = 0.7. As shown in Figure 4, the total energy solver for the flow across the plate (solid line 62) accurately recovers the total temperature, whereas the hybrid solver (LBM flow solver and finite difference based entropy solver) (dotted line 64) does not preserve the total temperature.

5 shows a graph of static temperature at a Mach speed (Ma) of Ma=0.7 for a flat plate at x=1.5 m from the leading edge. The graph shows temperature T as T/ _T∞ (x-axis) versus u as u/ _u∞ (y-axis). The disclosed total energy solver shown as solid line 66 is compared to a hybrid solver (LBM flow solver and finite difference based entropy solver) (dotted line 68). This FIG. 1 shows that the static temperature at a Mach speed (Ma) of Ma=0.7 for a flat plate at x=1.5 m from the leading edge matches the Waltz equation (asterisk line 70).

FIG. 6 shows a graph of thermal simulation results for a Converging-Diverging Verification (CDV) nozzle with back pressure P _b =0:75P _t . The graph shows temperature T as T/T _∞ (x-axis) versus x (y-axis). The disclosed total energy solver, shown as a solid line 72, is compared to a hybrid solver (LBM flow solver and finite difference based entropy solver) (dotted line 74). The disclosed total energy solver is in good agreement with the dashed line 76 of the analytical decision plot. A comparison of static temperatures shows that the finite difference based compressible solver fails to recover temperatures after impact, while the disclosed total energy solver achieves the expected values.

A system using the disclosed total energy solver can be used to determine the convective temperature distribution within the system. For example, if a system (formed of a volume represented by multiple voxels) contains a heat source and there is an airflow within the system, some regions of the system will be hotter than other regions based on their proximity to the airflow and the heat source. To model such a system, the temperature distribution within the system can be represented as a scalar quantity, with each voxel having an associated temperature.

In another example, the system can be used to determine the distribution of chemicals within a system. For example, if a system (formed of a volume represented by multiple voxels) contains a contamination source, such as a dirty bomb or chemicals or other particles suspended in the air or liquid, and there is an air or liquid flow within the system, some regions of the system will have higher concentrations than other regions based on the flow and proximity to the contamination source. To model such a situation, the distribution of chemicals within the system can be represented as a scalar quantity, with each voxel having an associated concentration.

In some applications, multiple different scalar quantities can be simulated simultaneously. For example, the system can simulate both the temperature distribution and the concentration distribution within the system.

Scalar quantities can be modeled in various ways. For example, the lattice Boltzmann (LB) method for solving the scalar transport equation can be used to indirectly solve the scalar transport. For example, the methods described herein can provide an indirect solution of the following second-order macroscopic scalar transport equation mentioned above:

In such configuration simulations, in addition to the lattice Boltzmann functions for fluid flow, a second set of distribution functions is introduced for the transport scalars. The method assigns vectors to each voxel in the volume to represent the fluid flow, and assigns scalars to each voxel in the volume to represent the desired scalar variables (e.g., temperature, density, concentration, etc.). The method fully recovers the macroscopic scalar transport equations that satisfy the exact conservation laws. The method is believed to improve the accuracy of the calculated scalars compared to other non-LBM methods. The method is also believed to improve the ability to consider complex boundary shapes.

This method of modeling scalar quantities can be used with time-explicit CFD/CAA solvers based on the Lattice Boltzmann Method (LBM), such as the PowerFLOW system available from Exa Corporation, Burlington, Massachusetts. Unlike methods based on discretizing the macroscopic continuum equations, LBM predicts macroscopic fluid dynamics starting from the "mesoscopic" Boltzmann equations of motion. The resulting compressible transient solution can be used to predict a variety of complex flow physics, including aeroacoustic and pure acoustic problems. Below, we provide a general discussion of LBM-based simulation systems, followed by a description of scalar solvers that can be used with fluid flow simulations to support such modeling methods.

Referring to FIG. 7, the first model (2D-1) 100 is a two-dimensional model that includes 21 velocities. Of these 21 velocities, one (105) represents a particle that is not moving, three sets of four 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) (110-113), twice the normalized velocity (2r) (120-123), or three times the normalized velocity (3r) (130-133), and two sets of four velocities represent particles moving with normalized velocity (r) (140-143) or twice the normalized velocity (2r) (150-153) with respect to both the x and y lattice axes.

Also as shown in FIG. 8, the second model (3D-1) 200 is a three-dimensional model that includes 39 velocities, each velocity represented by one of the arrows in FIG. 8. Of these 39 velocities, one 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 grid, eight velocities represent particles moving at a normalized velocity (r) with respect to all three of the x, y, and z grid axes, and twelve velocities represent particles moving at twice the normalized velocity (2r) with respect to two of the x, y, and z grid axes.

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

In the three-dimensional model 3D-2, which consists of 101 velocities, one represents particles that are not moving (group 1), three sets of six velocities represent particles that are 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 the x, y, or z axis of the grid (groups 2, 4, and 7), and three sets of eight velocities represent particles that are moving at the normalized velocity (r), twice the normalized velocity (2r), or three times the normalized velocity (3r) with respect to all three of the x, y, and z grid axes (groups 1, 2, and 3). 3, 8, and 10), 12 velocities represent particles moving at twice the normalized velocity (2r) about two of the x, y, and z grid axes (Group 6), 24 velocities represent particles moving at normalized velocity (r) and twice the normalized velocity (2r) about two of the x, y, and z grid axes and not moving about the remaining axis (Group 5), and 24 velocities represent particles moving at normalized velocity (r) about two of the x, y, and z grid axes and three times the normalized velocity (3r) about the remaining axis (Group 9).

In the two-dimensional model 2D-2, which consists of 37 velocities, one represents particles that are not moving (group 1), three sets of four velocities represent particles that are 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 either the x or y axis of the lattice (groups 2, 4, and 7), two sets of four velocities represent particles that are moving at normalized velocity (r) or twice the normalized velocity (2r) with respect to both the x and y lattice axes, eight velocities represent particles that are moving at normalized velocity (r) with respect to one of the x and y lattice axes and twice the normalized velocity (2r) with respect to the other axis, and eight velocities represent particles that are moving at normalized velocity (r) with respect to one of the x and y lattice axes and three times the normalized velocity (3r) with respect to the other axis.

The LBM models described above provide a particular class of efficient and robust discrete velocity dynamics models for the numerical simulation of flows in both two and three dimensions. This type of model contains a particular set of discrete velocities and weights associated with these velocities. These velocities correspond to Cartesian grid points in velocity space (in a non-Euclidean space) that facilitate accurate and efficient implementation of discrete velocity models, particularly of the class known as Lattice Boltzmann models. Such models can be used to simulate flows with high fidelity.

Referring to FIG. 9, a physical process simulation system operates according to procedure 300 to simulate a physical process such as a fluid flow. Prior to simulation, the simulation space is modeled as a collection 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 depict a microdevice located within a wind tunnel. The data generated by the CAD program is then processed to add a grid structure of appropriate resolution to account for objects and surfaces within the simulation space.

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 (I4)

The characteristic length of an object represents a large-scale feature of the object. For example, when simulating flow around a microdevice, the height of the microdevice can be considered as the characteristic length. When one is interested in the flow around a small object region (e.g., a car side mirror), one can increase the resolution of the simulation or use a high-resolution region around the region of interest. As the grid resolution increases, the voxel dimension decreases.

The state space is represented as _f (x,t), where _f represents the number of elements or particles per unit volume of state i at a lattice site (i.e., the density of particles in state i) represented by the three-dimensional vector x at time t. At a known time increment, the number of particles is simply represented as _f (x). The combination of all states at a lattice site is represented as f(x).

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

Each state i represents a different velocity vector at a particular energy level (i.e., energy level 0, 1 or 2). The velocity c of each state is denoted with "velocity" in each of the three dimensions as follows:
c _i =(c _ix , c _{i v} , c _iz ) Equation (I5)

A state of energy level zero represents a stopped 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 zero velocity 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 zero velocity in the other two.

Generating all possible permutations of the three energy levels gives a total of 39 possible states (1 energy 0 state, 6 energy 1 states, 8 energy 3 states, 6 energy 4 states, 12 energy 8 states and 6 energy 9 states).

Each voxel (i.e., each lattice site) is represented by a state vector f(x). The state vector completely defines the state 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. Using this set of rates, the system can generate Maxwell-Boltzmann statistics for the achieved equilibrium state vector. For total energy, the second variance q _i (x,t) uses the same number of entries (39 states) to represent the specific energy. This is used together with the state vector f(x) to solve the total energy conservation.

Voxels are grouped into 2x2x2 volumes called microblocks for processing efficiency. These microblocks are organized to enable parallel processing of the voxels and minimize the overhead associated with the data structure. The abbreviation for a voxel within a microblock is defined as _Ni (n), where n represents the relative position of the lattice site within the microblock, nε{0, 1, 2,..., 7}. A microblock is shown in Figure 10.

11A and 11B, a surface S (FIG. 11A) is represented in a simulation space (FIG. 11B) as a collection of facets F _α .
S = {F _α } Equation (I6)
where α is an index that enumerates a particular facet. Facets are not restricted to voxel boundaries and typically have a size 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 purpose of implementing surface dynamics. Specifically, each facet _Fα has a unit normal ( _nα ), a surface area ( _Aα ), a center position ( _xα ), and a facet distribution function (f _i (α)) that describes the surface dynamics of the facet. The total energy dispersion function q _i (α) is treated in the same way as the flow dispersion for the facet-voxel interaction.

With reference to FIG. 12, different levels of resolution can be used in different regions of the simulation space to increase processing efficiency. Typically, the region 650 around the object 655 is of most interest, and therefore this region is simulated with the highest resolution. Since the effect of viscosity decreases with distance from the object, a decreasing level of resolution (i.e., an expanded voxel volume) is used to simulate regions 660, 665 located at increasing distances from the object 655. Similarly, as shown in FIG. 13, a region 770 around less important features of the object 775 can be simulated with a lower level of resolution, while a region 780 around the most important features of the object 775 (e.g., the leading and trailing faces) can be simulated with the highest level of resolution. Regions 785 far from the center are simulated with the lowest level of resolution and the most voxels.

C. Identifying Voxels Influenced by Facets Referring again to FIG. 9, once the simulation space has been modeled (step 302), voxels influenced by one or more facets are identified (step 304). A voxel can be influenced by a facet in several ways. First, a voxel intersected by one or more facets is affected in that it has a reduced volume compared to non-intersected voxels. This is because the facets and the material below the surface that the facets represent occupy a portion of the voxel. A fractional factor _Pf (x) indicates the portion of the voxel that is not influenced by the facets (i.e., the portion that can be occupied by the fluid or other material whose flow is to be simulated). For non-intersected voxels, _Pf (x) is equal to 1.

Voxels that interact with one or more facets by either transferring particles to the facet or receiving particles from the facet are also identified as voxels influenced by the facet. Every voxel that is 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, further voxels will also contain such states.

Referring to FIG. 14, for each state i with a nonzero velocity vector c, a facet _Fα receives particles from or transfers particles to a region defined by a parallelepiped _Gα having a height defined by the magnitude of the vector dot product (|c _i n _i |) of the velocity vector C _i and the unit normal n _α of the facet, and a base defined by the surface area _Aα of the facet, such that its volume V _iα is equal to:
V _iα = | c _i n _α | A _α Equation (I7)

A facet _Fα receives particles from the volume _Vα when the state's velocity vector points towards it (|c _i n _i |<0), and transfers particles into this region when the state's velocity vector points away from it (|c _i n _i |>0). As we will see later, this equation is modified when another facet occupies a portion of the parallelepiped _Gα , a situation that can occur near a non-convex feature such as an interior corner.

The parallelepiped G _iα of the facet F _α can overlap some or all of the voxels. The number of voxels or the fraction depends on the size of the facet relative to the size of the voxel, the energy of the states 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.

The fraction N(x) of a _voxel overlapped by a parallelepiped G is defined as _V (x). Using this term, the flux _Γ (x) of a particle in state i moving between voxel N(x) and facet _F is equal to the density of particles in state i at the voxel ( _N (x)) multiplied by the volume of the region of overlap of the voxels ( _V (x)).
Γ _iα (x) = N _i (x) V _iα (x) Equation (I8)

When one or more facets intersect the parallelepiped G _iα , the following conditions apply:
V _iα =ΣV _α (x) +ΣV _iα (β) Equation (I9)
where the first sum corresponds to all voxels where _G overlaps, and the second term corresponds to all facets that intersect _G. When the parallelepiped G does not intersect _with any other facets, this equation reduces to:
V _iα =ΣV _iα (x) Equation (I10)

D. Running the Simulation Once voxels influenced by one or more facets have been identified (step 304), a timer is initialized and the simulation begins (step 306). During each time increment of the simulation, the movement of particles from voxel to voxel is simulated by an advection phase (steps 308-316) that accounts for the interaction of the particles with the surface facets. Next, the interaction of particles within each voxel is simulated by a collision phase (step 318), based in part on the use of the entropy solver described above. The timer is then incremented (step 320). If the incremented timer does not indicate that the simulation is complete (step 322), 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. Surface Boundary Conditions In order to properly 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 the net mass flux at each energy level (i.e., energy levels 1 and 2) to be equal to zero.

The other two boundary conditions relate to the net momentum of the particles interacting with the facet. At surfaces without skin friction, referred to herein as slip surfaces, the net tangential momentum flux is equal to zero and the net normal momentum flux is equal to the local pressure on the facet. Thus, the components of the total momentum received perpendicular to the facet normal nα are equal to the components of the total momentum transferred perpendicular to the facet normal nα (i.e., the tangential components), and the difference between the components of the total momentum received parallel to the facet normal nα and the components of the total momentum transferred parallel to the facet normal nα (i.e., the normal components) is equal to the local pressure on the facet. At non-slip surfaces, surface friction reduces the total tangential momentum of the particles transferred by the facet relative to the total tangential momentum of the particles received by the facet by a factor related to the amount of friction.

2. Voxel to Facet Collection As a first step in simulating particle-surface interactions, particles are collected from voxels and provided to facets (step 308). As discussed above, the flux of particles in state i between voxel N(x) and facet F is:
Γ _iα (x) = N _i (x) V _iα (x) Equation (I11)

From this equation, for each state i pointing towards a facet (c _i n _α <0), the number of particles contributed by the voxel to the facet F _α is:
Γ _iαV → _F = Σ _x Γ _iα (x) = Σ _x N _i (x) V _iα (x) Equation (I12)

Only voxels for which V _iα (x) has a nonzero value need to be summed. As mentioned above, the facet size is chosen so that V _iα (x) has a nonzero value for only a small number of voxels. Since V _iα (x) and P _f (x) can also have non-integer values, Γ _α (x) is stored and processed as a real number.

3. Facet-to-Facet Transfer Next, the particles are transferred between facets (step 310). If a facet F _β intersects a parallelepiped G _iα where facet F _α is in the inflow state (c _i n _α <0), then a portion of the state i particles received by facet F _α will be inflowing from facet F _β . Specifically, facet F _α will receive a portion of the state i particles generated by facet F _β during the previous time increment. This relationship is illustrated in FIG. 17, where a portion 1000 of parallelepiped G _iα intersected by facet F _β is equal to a portion 1005 of parallelepiped G _iβ intersected by facet F _α . As mentioned above, the intersected portion is denoted as V _iα (β). Using this term, the flux of state i particles between facets F _β and F _α can be expressed as:
Γ _iα (β,t−1)=Γ _i (β) V _iα (β)/V _iα Equation (I.13)
where Γ _i (β,t−1) is a measure of the number of particles in state i produced by facet F _β during the previous time increment. From this equation, for each state i pointing towards facet F _α (c _i n _α <0), the number of particles donated to facet F _α by other facets is
Γ _iαF→F = Σ _β Γ _iα (β) = Σ _β Γ _i (β,t-1) V _iα (β)/V _iα Equation (I.14)
The total flux of particles of state i into the facet is
Γ _iIN (α) = Γ _iαF→F + Γ _iαF→F = Σ _x N _i (x) V _iα + Σ _β Γ _i (β,t−1) V _iα (β)/V _iα Equation (I.15)

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 N(α) of the facet distribution function are set equal to the flux of particles entering these states divided by the volume V _iα ,
When c _i n _α <0, we have the following:
N _i (α)=Γ _iIN (α)/V _iα Equation (I.16)

The facet distribution function is a simulation tool for generating 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. Populating the outward states using the method for populating inward states described above results in:
If c _i n _α ≧0,
N _i (α) = Γ _iOTHER (α) / V _iα Equation (I.17)
where Γ _iOTHER (α) is determined using the method for generating Γ _iIN (α) described above, but applying this method to states (c _i n _α ≧0) other than the admission state (c _i n _α <0). Alternatively, the value of Γ _iOUT (α) from the previous time step can be used to generate Γ _iOTHER (α) as follows:
Γ _iOTHER (α,t) = Γ _iOUT (α,t-1) Equation (I.18)

In the equilibrium state (c _i n _α =0), both V _iα and V _iα (x) are zero. In the equation for _{N i} (α), V _iα (x) appears in the numerator (from the equation for Γ _iOTHER (α)) and V _iα appears in the denominator (from the equation for N _i (α)). Thus, the equilibrium state N _i (α) is determined as the limit of N _i (α) as V _iα and V _iα (x) approach zero. Values for 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 beginning of the simulation based on the initial conditions of temperature and pressure. These values are then adjusted over time.

方程式（Ｉ．１９）
この式から、以下のように法線運動量Ｐ_n（α）を求める。
Ｐ_n（α）＝ｎ_α・Ｐ（α） 方程式（Ｉ．２０） 4. Execution of facet surface dynamics Next, execute surface dynamics so that each facet satisfies the four boundary conditions mentioned above (step 312). Figure 18 shows the procedure for executing surface dynamics of a facet. First, determine the total momentum perpendicular to the facet _Fα (step 1105) by determining the total momentum P(α) of the particles at the facet as follows:
For all i,

Equation (I.19)
From this equation, the normal momentum P _n (α) is calculated as follows.
_Pn (α)= _nα ·P(α) Equation (I.20)

This normal momentum is then removed (step 1110) to produce N _n- (α) using a pushing/pulling method, whereby particles transition between states in a manner that only affects the normal momentum. The pushing/pulling method is described in U.S. Patent No. 5,594,671, which is incorporated by reference.

The particles of N _n− (α) are then collided to produce a Boltzmann distribution N _n−β (α), step 1115. As described below with respect to a fluid dynamics implementation, the Boltzmann distribution can be achieved by applying a set of collision rules to N _n− (α).

Next, the outgoing flux distribution of the facet F _α is calculated based on the incoming flux distribution and the Boltzmann distribution (step 1120). First, the difference between the incoming flux distribution Γ _i (α) and the Boltzmann distribution is calculated as follows.
ΔΓ _i (α) = Γ _iIN (α) - _{N n -βi} (α) V _iα Equation (I.21)

方程式（Ｉ．２４）
ここでのｊは、エネルギーレベル１状態では１に等しく、エネルギーレベル２状態では２に等しい。 Using this difference, the outflow flux distribution is
When n _α c _i >0, it becomes:
Γ _iOUT (α)=N _n−βi (α)V _iα −. Δ. Γ _i *(α) Equation (I.22)
where i* is the state with the opposite direction to state i. For example, if state i is (1,1,0,0), then state i* is (-1,-1,0,0). The outflow flux distribution can be further refined to account for skin friction and other factors:
When nαci>0, the following occurs:
Γ _iOUT (α) = N _n-Bi (α) V _{iα - ΔΓ i} * (α) + C _f (n _α · c _i ) - [N _n-βi * (α) - _{N n-βi} (α)] V _iα + (n _α · c _i ) (t _{1 α} · c _i ) ΔN _j,1 V _iα + (n _α · c _i ) (t _{2 α} · c _i ) ΔN _j,2 V _iα Equation (I.23)
where _C is a function of skin friction, _t is a first tangent vector perpendicular to _n , _t is a second tangent vector perpendicular to both _n and _t , and _ΔN and _ΔN are distribution functions corresponding to the energy of state i (j) and the indicated tangent vector. The distribution functions are determined according to the following equation:

Equation (I.24)
Here j is equal to 1 for the energy level 1 state and 2 for the energy level 2 state.

The functions of each term in the equation of Γ _iOUT (α) are as follows: The first and second terms enforce the normal momentum flux boundary condition as long as collisions are effective in producing a Boltzmann distribution, but include the exception of the tangential momentum flux. The fourth and fifth terms correct for exceptions that may occur due to discreteness effects or non-Boltzmann structures resulting from this insufficient collisions. Finally, the third term adds an amount of surface friction to enforce the desired change in tangential momentum flux on the surface. The generation of the friction coefficient C _f is explained below. Note that all terms involved in the vector manipulations are geometric factors that can be calculated before the start of the simulation.

方程式（Ｉ．２６） From this formula, the tangential velocity is calculated as follows:
u _i (α)=(P(α)-P _n (α)n _α )/ρ Equation (I.25)
where ρ is the density of the facet distribution such that

Equation (I.26)

Similarly to the above, the difference between the inflow flux distribution and the Boltzmann distribution is calculated as follows:
ΔΓ _i (α) = Γ _iIN (α) - _{N n -βi} (α) V _iα Equation (I.27)

Then, the outflow flux distribution is:
Γ _iOUT (α) = N _n-βi (α) V _iα - ΔΓ _i * (α) + C _f (nαci) [N _n-βi * (α) - _{N n-βi} (α)] V _iα Equation (I.28)
This equation corresponds to the first two rows of the exit flux distribution obtained by the previous method, but does not require correction for the anomalous tangential flux.

方程式（Ｉ．２９）
となり、ここでのｐ_αは、ファセットＦ_αにおける平衡圧であって、ファセットに粒子を提供するボクセルの平均密度及び温度値に基づき、ｕ_αは、ファセットにおける平均速度である。 Using either method, the resulting flux distribution satisfies all the momentum flux conditions, i.e.

Equation (I.29)
where p _α is the equilibrium pressure at the facet F _α , based on the average density and temperature values of the voxels providing particles to the facet, and u _α is the average velocity at the facet.

方程式（Ｉ．３０）
ここでの添字ｊは、状態ｉのエネルギーを示す。次に、このエネルギー差を用いて差分項（ｄｉｆｆｅｒｅｎｃｅ ｔｅｒｍ）を生成すると、
ｃｊｉｎα＞０の場合、以下のようになる。

方程式（Ｉ．３１）
この差分項を用いて流出流束を修正すると、流速は、
ｃ_jiｎ_α＞０の場合、以下のようになる。
Γ_αjiOUTf＝Γ_αjiOUT＋δΓ_αji 方程式（Ｉ．３２）
この演算は、接線運動量流束をそのままにした状態で質量及びエネルギー流束を訂正する。ファセットの近傍において流れがおおよそ均一であって平衡に近い場合には、この調整はわずかなものである。結果として得られる調整後の法線運動量流束は、近傍の平均特性に近傍の非均一性又は非平衡特性に起因する補正をプラスしたものに基づく平衡圧である値へとわずかに変化する。 To ensure that the mass and energy boundary conditions are satisfied, the difference between the input and output energies for each energy level j is measured as follows:

Equation (I.30)
Here, the subscript j denotes the energy of state i. Next, this energy difference is used to generate a difference term,
If cjinα>0, then:

Equation (I.31)
When the outflow flux is corrected using this difference term, the flow velocity becomes
When c _ji n _α >0, we have the following.
_ΓαjiOUTf = _ΓαjiOUT + _δΓαji Equation (I.32)
This operation corrects the mass and energy fluxes while leaving the tangential momentum flux unchanged. If the flow is roughly uniform and close to equilibrium in the vicinity of the facet, the adjustment is small. The resulting adjusted normal momentum flux changes slightly to a value that is an equilibrium pressure based on the average properties of the neighborhood plus a correction due to non-uniform or non-equilibrium properties of the neighborhood.

5. Voxel-to-Voxel Movement Referring again to Figure 9, particles are moved between voxels along a 3D rectilinear grid (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 near the surface). In a typical simulation, voxels that are not near enough to the surface to interact with the surface make up the majority of voxels.

Each of these discrete states 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 of movement of the particle along the corresponding axis.

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

This movement operation becomes a bit more complicated for voxels that interact with one or more surfaces. This can result in one or more fractional particles moving to the facet. Moving these fractional particles to the facet causes the voxel to remain filled with fractional particles. Move these fractional particles to the voxels occupied by the facet.

Also shown is total energy correction to advection 315, in which the computer corrects the state vectors of the particles by adding specific total energy values to the states of particles that were advected and subtracting specific total energy values from the states of particles that were not advected over the time interval (see FIG. 3).

方程式（Ｉ．３３）
ここでのＮ（ｘ）はソースボクセルである。 15, when a portion 900 of state i particles in voxel 905 have been moved to facet 910 (step 308), the remaining portion 915 is moved to voxel 920 in which facet 910 is located, from where state i particles are directed to facet 910. Thus, if the state population is equal to 25 and V _iα (x) is equal to 0.25 (i.e., one-quarter of the voxels intersect the parallelepiped G _iα ), then 6.25 particles are moved to facet F _α and 18.75 particles are moved to voxels occupied by facet F _α . Since a voxel may be intersected by multiple facets, the number N(f) of state i particles moving to voxels occupied by one or more facets is given by:

Equation (I.33)
where N(x) is the source voxel.

方程式（Ｉ．３４）
ここでのＰ_f（ｘ）は、部分的なボクセルの体積縮小に相当する。この式から、各状態ｉについて、ファセットからボクセルに向かう粒子の総数Ｎ_(x)は以下のようになる。

方程式（Ｉ．３５） 6. Scattering from Facets to Voxels The effluent particles from each facet are then scattered into voxels (step 316). Essentially, this step is the inverse of the collection step, where particles move from voxels to facets. The number of particles in state i that move from facet _Fα to voxel N(x) is given by:

Equation (I.34)
Here, P _f (x) corresponds to the partial voxel volume shrinkage. From this formula, for each state i, the total number of particles N _(x) going from the facet to the voxel is given by:

Equation (I.35)

After scattering particles from facets into voxels, combining these with particles advected from surrounding voxels, and integerizing the result, certain directions in certain voxels may either underflow (go negative) or overflow (go above 255 in an 8-bit implementation). This results in mass, momentum, and energy increasing or decreasing after truncating these quantities to fit within the range of acceptable values.

To prevent this from happening, accumulate mass, momentum, and energy that are in excess of the limits before truncating the harmful state. For energy with states, mass equal to the value gained (due to underflow) or lost (due to overflow) is added back to a randomly (or sequentially) selected state with the same energy, which is itself not affected by overflow or underflow. The further momentum resulting from this mass and energy addition is accumulated and added to the momentum from the truncation. By adding mass only to states with the same energy, both mass and energy are corrected when the mass counter reaches zero. Finally, the momentum is corrected using a pushing/pulling method until the momentum accumulator returns to zero.

7. Perform Fluid Dynamics Perform fluid dynamics (FIG. 9, 318). This step can be referred to as a microdynamics operation or an intravoxel operation. Similarly, the advection procedure can be referred to as an intervoxel operation. The microdynamics operations described below can also be used to collide particles at the facets to generate a Boltzmann distribution.

The 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 dynamics of distribution in real fluid systems. The collision process can be fully 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 equilibrium distribution function, denoted by f ^eq in Equation (2), is fully specified by Equation (4). The selection of the velocity vector set c _i , weights, both listed in Table 1, together with Equation 2, ensures that the macroscopic behavior follows the correct fluid dynamics equations.

E. Variable Resolution Referring to Figure 12, variable resolution (as described above) employs different sized voxels, hereinafter referred to as coarse voxels 12000 and fine voxels 1205. (Although the following description refers to voxels having two different sizes, it should be understood that the described techniques may also be applied to three or more different sized voxels to provide additional levels of resolution.) The interface between the regions of coarse voxels and fine voxels is referred to as the variable resolution (VR) interface 1210.

Employing variable resolution at or near the surface allows facets to interact with voxels on either side of the VR interface. These facets are classified as VR interface facets 1215 (F _αIC ) or VR microfacets 1220 (F _αIF ). VR interface facets 1215 are facets located on the coarse side of the VR interface and have a coarse parallelepiped 1225 that extends into the microvoxel. (A coarse parallelepiped is one whose ci is dimensioned according to the dimensions of the coarse voxel, and a microparallelepiped is one whose ci is dimensioned according to the dimensions of the microvoxel.) VR microfacets 1220 are facets located on the micro side of the VR interface and have a microparallelepiped 1230 that extends into the coarse voxel. Processing related to interface facets can also involve interactions with coarse facets 1235 (F _αC ) and microfacets 1240 (F _αF ).

For both types of VR facets, surface dynamics are implemented at the small scale and behave as described above. However, VR facets differ from other facets in the way particles are advected to and from the VR facet.

Interaction with VR facets is performed using a variable resolution procedure 1300 shown in FIG. 13. Most steps of this procedure are performed using the equivalent steps described above for interaction with non-VR facets. Procedure 1300 is performed during a coarse time step (i.e., a period corresponding to a coarse voxel) that includes two phases, each corresponding to an atomic time step. Facet surface dynamics are performed during each atomic time step. Thus, the VR interface facet F _αIc is considered as two identically sized and identically oriented atomic facets, called the black facet F _αIcb and the red facet F _αICr , respectively. The black facet F _αICb is associated with the first atomic time step within the coarse time step, and the red facet F _αICr is associated with the second atomic time step within the coarse time step.

First, the particles are moved (advected) between the facets by a first surface-to-surface advection stage (step 1302). The particles are moved from the black facet F _αICb to the coarse facet F _βC with a weighting factor V _−αβ , which corresponds to the volume of the non-blocking part of the coarse parallelepiped (FIG. 12, 1225) extending from the facet F _α and lying behind _the facet F _β , minus the non-blocking part of the minute parallelepiped (FIG. 12, 1245) extending from the facet F _α and lying behind the facet F _β . The size of the c _i of the minute voxel is 1/2 the size of the c _i of the coarse voxel. As mentioned above, the volume of the parallelepiped of the facet F _α is defined as follows:
V _iα = | c _i n _α | A _α Equation (I.36)

Thus, since the surface area _Aα of a facet does not change between the rough and the microparallelopehedrons, and the unit normal _nα always has a magnitude of 1, the volume of the microparallelopehedron corresponding to a facet is 1/2 the volume of the rough parallelepiped corresponding to this facet.

The grains migrate from the rough facet FαC to the black facet _FβICb with a weighting factor of Vαβ, _which corresponds to the volume of the non-blocked portion of the microparallelope hexahedron extending from the facet _Fα and present behind the facet _Fβ .

The particles move from the red facet F _αICr to the coarse facet F _βC with a weighting factor of V _αβ , and from the coarse facet F _αC to the red facet F _βICb with a weighting factor of V _−αβ .

Particles move from the red facet F _αICr to the black facet F _βICb with a weighting factor of V _αβ . No advection from black to red occurs at this stage. Also, since the black and red facets represent successive time steps, advection from black to black (or from red to red) never occurs. For the same reason, particles move at this stage from the red facet F _αICr to the microfacet F _βIF or F _βF with a weighting factor of V _αβ , and from the microfacet F _αIF or F _αF to the black facet F _αICb with the same weighting factor.

Finally, the grains move from the microfacet F _αIF or F _αF to the other microfacet F _βIF or F _βF with the same weighting factor, and from the rough facet F _αC to the other rough facet F _C with a weighting factor of V _Cαβ , which corresponds to the volume of the non-block part of the rough _{parallelepiped} extending from the facet F _α and lying behind the facet F _β .

After the particles are advected between the surfaces, the particles are collected from the voxels in a first collection phase (steps 1304-1310). Particles of the microfacet F _αF are collected from the microvoxels using a microparallelopephron (step 1304), and particles of the coarse facet F _αC are collected from the coarse voxels using a coarse parallelephron (step 1306). Next, particles of the black facet F _αIRb and the VR microfacet F _αIF are collected from both the coarse and the microvoxels using a microparallelopephron (step 1308). Finally, particles of the red facet F _αIRr are collected from the coarse voxels using the difference between the coarse parallelephron and the microparallelopephron (step 1310).

Next, the coarse voxels that interact with the microvoxel or the VR facet are expanded into a set of microvoxels (step 1312). The coarse voxels are expanded in a state that moves particles to the microvoxel within a single coarse time step. For example, a coarse voxel in a suitable state where no facets are intersected is expanded into eight microvoxels oriented like the microblock in FIG. 4. A coarse voxel in a suitable state where one or more facets are intersected is expanded into a set of full and/or partial microvoxels that correspond to the portion of the coarse voxel that is not intersected by any facets. The particle densities N _i (x) are equal for the coarse voxel and the microvoxels resulting from the expansion of the coarse voxel, but the microvoxels may have a fractional coefficient P _f that is different from the fractional coefficient of the coarse voxel and the fractional coefficients of the other microvoxels.

Surface dynamics are then performed for the small facets F _αIF and F _αF (step 1314), and for the black facet F _αICb (step 1316). The dynamics are performed using the procedure shown in FIG.

Next, the particles are moved between the microvoxels, including the actual microvoxels and the microvoxels resulting from the unfolding of the coarse voxels (step 1318). After the particles are moved, the particles are scattered from the microfacets F _αIF and F _αF into the microvoxels (step 1320).

Particles are also scattered from the black facet F _αICb to tiny voxels (including tiny voxels resulting from the unfolding of the coarse voxels) (step 1322). Particles are scattered to tiny voxels if they would have received the particles if there was no surface present at that time. Specifically, particles are scattered to voxel N(x) when the voxel is an actual tiny voxel (as opposed to a tiny voxel resulting from the unfolding of the coarse voxel), when voxel N(x+c _i ), which is one velocity unit N(x) beyond voxel N(x), is an actual tiny voxel, or when voxel N(x+c _i ), which is one velocity unit N(x) beyond voxel N(x), is a tiny voxel resulting from the unfolding of the coarse voxel.

Finally, the first micro time step is completed by performing fluid dynamics on the micro voxels (step 1324). The voxels on which fluid dynamics is performed do not include the micro voxels that result from the unfolding of the coarse voxels (step 1312).

Procedure 1300 performs similar steps during the second microtime step. First, in a second surface-to-surface advection phase, particles move between surfaces (step 1326). Particles are advected from the black facets to the red facets, from the black facets to the microfacets, from the microfacets to the red facets, and from the microfacets to the microfacets.

After the particles are advected between the surfaces, the particles are collected from the voxels in a second collection phase (steps 1328-1330). Particles of the red facet F _αIRr are collected from the small voxels using a small parallelepiped (step 1328). Particles of the small facets F _αF and F _αIF are also collected from the small voxels using a small parallelepiped (step 1330).

Thus, as described above, surface dynamics are performed for the fine facets F _αIF and F _αF (step 1332), for the coarse facet F _αC (step 1134), and for the red facet F _αICR (step 1336).

Next, the particles are moved between voxels using the fine resolution such that the particles move between the fine voxels and the fine voxels that represent the coarse voxels (step 1338). Then, the particles are moved between voxels using the coarse resolution such that the particles move between the coarse voxels (step 1340).

Next, in a combination step, particles are scattered from the facets to the voxels while the microvoxels representing the coarse voxels (i.e., the microvoxels resulting from the expansion of the coarse voxels) are merged into the coarse voxels (step 1342). In this combination step, particles are scattered from the coarse facets to the coarse voxels using a coarse parallelepiped, particles are scattered from the microfacets to the microvoxels using a microparallelepiped, particles are scattered from the red facets to the micro or coarse voxels using a microparallelepiped, and particles are scattered from the black facets to the coarse voxels using the difference between the coarse parallelepiped and the microparallelepiped. Finally, fluid dynamics is performed on the micro and coarse voxels (step 1344).

F. Scalar Transport Solvers As mentioned above, to solve the fluid flow, various types of LBMs can be applied that act as background carriers for the scalar transport. During the simulation, both the fluid flow and the scalar transport are simulated. For example, the fluid flow is a flow that is simulated using the Lattice Boltzmann (LB) method, and the scalar transport is simulated using a set of distribution functions, referred to herein as the scalar transport equations.

方程式（Ｉ．３７） Although a detailed description of the LBM method for simulating fluid flow is provided herein, the following is an example of one fluid flow simulation method that can be used in conjunction with scalar simulation.

Equation (I.37)

方程式（Ｉ．３８）
ここではＴ₀＝１／３である。離散格子速度ｃ_iは以下の通りであり、

方程式（Ｉ．３９）
残りの粒子ではｗ₀＝１／３であり、デカルト方向の状態ではｗ_i＝１／１８であり、二重対角方向の状態ではｗ_i＝１／３６である。ｇ_i（ｘ，ｔ）は、外部体積力項である。流体力学的量￡及びｕは、粒子分布関数のモーメントである。

方程式（Ｉ．４０） where f _i (x,t) (i=1, . . . , 19) are particle distribution functions, τ is the unity relaxation time, and f _i ^eq (x,t) is the equilibrium distribution function containing the third order evolution of the fluid velocity.

Equation (I.38)
Here, T ₀ = 1/3. The discrete grid velocities c _i are:

Equation (I.39)
For the remaining particles w ₀ = 1/3, for the Cartesian states w _i = 1/18 and for the bidiagonal states w _i = 1/36. g _i (x,t) are the external body force terms. The hydrodynamic quantities £ and u are the moments of the particle distribution function.

Equation (I.40)

As described above, the fluid solver is used in conjunction with the scalar transport solver that generated the scalar transport information. Thus, for scalar transport, a separate set of distribution functions T _i is introduced in addition to the fluid solver. Thus, the system simulates both fluid flow and scalar transport for each voxel of the system to generate state vectors representing fluid flow and scalar quantities representing scalar variables. These simulation results are stored as entries in a computer accessible medium.

方程式（Ｉ．４１）

方程式（Ｉ．４２）

方程式（Ｉ．４３）
Ｔ_iはスカラー分布関数であり、Ｔは解かれるスカラーである。τ_Tは、スカラー拡散率に対応する緩和時間である。緩和時間は、システムが緩和して平衡になるまでにどれほど掛かるかについての尺度を示す。項ｆ_i、ρ、Ｔ₀及びｕは、それぞれ方程式（３８）、（３９）及び（４１）で定義される。 The set of scalar transport functions provides an indirect solution of the second-order macroscopic scalar transport equation: _{T i} provides the equations that model the dynamic evolution of the scalar quantities.

Equation (I.41)

Equation (I.42)

Equation (I.43)
T _i is a scalar distribution function and T is the scalar to be solved. τ _T is the relaxation time corresponding to the scalar diffusivity. The relaxation time provides a measure of how long it takes for a system to relax to equilibrium. The terms f _i , ρ, T ₀ and u are defined in equations (38), (39) and (41), respectively.

The grid velocity set represented by c _i can be a discrete set of grid velocities used in the scalar simulation. In general, the grid velocity set for the scalar distribution does not need to be the same as the grid velocity set for the fluid distribution, since the scalar solver is an additional system attached to the basic fluid solver. For example, a simulation of scalar evolution can use a smaller number of grid velocities than a simulation of fluid flow. A different grid velocity set for the scalar can be applied as long as the scalar grid velocity set is a subset of the fluid grid velocity set. For example, if a 19-velocity LB model is used for the fluid simulation, a 6-velocity LB model can be used for the scalar simulation. The 19-velocity LB model has a higher order of lattice symmetry than the 6-velocity LB model, so the examples below use the same 19-velocity lattice model for scalars.

Standard well-known BGK (e.g., as described above) includes all orders of non-equilibrium moments. It is believed that including all non-equilibrium moments is not necessarily isotropic, hydrodynamically or physically meaningful. Therefore, a form of the BGK regularized/filtered collision operator is used. The collision operator Φ _i represents the future collision coefficient. This collision operator extracts non-equilibrium scalar features only in the relevant supported orders (e.g., only the first order). This operator preserves and relaxes the modes of interest, while the non-equilibrium features associated with unsupported/undesired higher order modes are eliminated. This scheme is sufficient for recovering the physical properties of scalar transport (e.g., advection and diffusion). The use of this future collision operator is believed to significantly reduce noise, provide better advection behavior (e.g., can enhance Galilean invariance), and be stable compared to other solutions of the well-known BGK operator. Such a form (e.g., as shown in Eq. 43) ensures that only first order non-equilibrium moments contribute to scalar diffusion in the hydrodynamic range. Any higher-order non-equilibrium moments are filtered out by this collision process. The use of the collision operator Φ _i described above is believed to provide advantages including the elimination of numerical noise and improved robustness as shown in BGK. The scalar T acts as a unique equilibrium, and no complex expressions for the scalar equilibrium distribution functions are required. The overall computation of the collision operator Φ _i is fairly efficient. It is believed that filtering higher-order non-equilibrium moments can provide the additional advantage of reducing aliasing that may be present in higher-order equilibrium solutions.

方程式（Ｉ．４４）
となり、この結果、

方程式（Ｉ．４５）
となることが分かり、ここでのＴ’_i（ｘ，ｔ）は方程式（４２）の右辺を示す。従って、スカラーを温度とみなした場合、スカラー衝突演算子は局所的なρＴを保存し、このことは局所的エネルギー保存の実現を意味する。Ｔ_iはｆ_iと共に伝播するので、移流中にエネルギー分布Ｅ＝ｆ_iＴ_iが完全に維持される。従って、ρＴの大域的保存が達成される。さらに、最も注目すべき点として、このスキームは、スカラーＴの均一性について厳密な不変性を維持する。これについては、あらゆる場所で

である場合、背景の流れ場にかかわらずその後の全ての時点においてあらゆる場所でφ（ｘ，ｔ）＝０かつ

であることが容易に分かる。この基本特性は、過去のどの格子ボルツマンスカラーモデルでも実証されていない。 It can be seen that the collisions in equation (42) obey the scalar conservation law. Multiplying both sides of equation (42) by _f'i (x,t)= _fi (x+ _cj ,t+1), we obtain

Equation (I.44)
As a result,

Equation (I.45)
where T' _i (x,t) denotes the right hand side of equation (42). Thus, if we consider the scalar as temperature, the scalar collision operator conserves the local ρT, which means realizing local energy conservation. Since T _i propagates with f _i , the energy distribution E = f _i T _i is perfectly preserved during advection. Thus, global conservation of ρT is achieved. Moreover, and most notably, this scheme maintains strict invariance for the uniformity of the scalar T, which can be seen everywhere.

If φ(x,t)=0 everywhere at all subsequent times, regardless of the background flow field,

It is easy to see that this fundamental property has not been demonstrated in any previous lattice Boltzmann scalar model.

方程式（Ｉ．４６）
κ＝（τ_T－１／２）Ｔ₀であることが分かる。この均一性不変条件は、ρが∇Ｔの範囲外にあることを確実にする。 Using the Chapman-Enskog expansion, equation (42) recovers the second-order macroscopic scalar transport equation:

Equation (I.46)
It can be seen that κ=(τ _T −1/2)T _0. This uniformity invariant ensures that ρ is outside the range of ∇T.

Boundary Conditions One substantial advantage of LBM is its ability to accommodate complex geometries. In a generalized volumetric LB surface algorithm for achieving frictionless/frictionless boundary conditions (BC) of arbitrary geometries, mass is conserved and the tangential and normal momentum fluxes on the boundary are realized exactly. Local detailed balance is completely satisfied. As a direct extension of this method, adiabatic (zero scalar flux) BC of any geometry of scalars can be derived. Once the adiabatic BC is realized, a prescribed finite flux BC can be achieved.

方程式（Ｉ.４７）
であり、ここでは、

方程式（Ｉ.４８）
Ｐ_i（≡｜ｎ・ｃ_i｜Ａ）
ｃ_i・nｎ＜０，ｃ_iｎ＝－ｃ_i・Ｐ_i＝Ｐ_i・nである。 Unlike other point-wise LB, boundary conditions are enforced on a discretized set of surface elements. These piece-wise flat surface elements together represent a curved shape. As shown in FIG. 14, during particle advection, each surface element collects inflow particles from adjacent fluid cells (step 1410). The inflow distributions f _i ⁱⁿ , T _i ⁱⁿ are weighted by the volumetric overlap of the parallelepiped from the underlying surface element and the cell in the particle movement direction (step 1412). After receiving the inflow, the following surface scalar algorithm is applied (step 1414): To obtain the outflow distribution from the surface,

Equation (I.47)
where:

Equation (I.48)
P _i (≡ |n·c _i |A)
c _{i ·n} n < 0, c _i n = -c _{i ·} P _i =P _{i ·n} .

方程式（Ｉ．４９） where n is the surface normal pointing towards the fluid domain, c _{i ·n} n<0, c _{i ·} P _i (≡|n ·c _i |A) is the volume of the parallelepiped in particle direction ci associated with surface normal n and area A of a given surface element, obviously P _i =P _{i *} . The outflow distribution is back-propagated from the surface elements to the fluid cells following the same surface advection process (step 1416). It is not difficult to show that the above surface scalar collisions achieve a strictly zero surface scalar flux. Summing over the outflow direction, the outflow scalar flux becomes:

Equation (I.49)

であるため、総流出スカラー流束は総流入スカラー流束と同じである。

方程式（Ｉ．５０） Note that P _i =P _{i *} and the definition of T ⁱⁿ equation (49), the second summation term on the right hand side is zero. Also, due to the conservation of mass flux,

Therefore, the total outflow scalar flux is the same as the total inflow scalar flux.

Equation (I.50)

方程式（Ｉ．５１） Therefore, the zero net surface flux (adiabatic) BC is perfectly satisfied for any geometry. If an external scalar source Q(t) is specified on the surface, one can add the source term directly to equation (48).

Equation (I.51)

方程式（Ｉ．５２） If the boundary conditions have a prescribed scalar quantity T _w (eg, surface temperature), then the surface heat flux can be calculated according to:

Equation (I.52)

Numerical Validation Figures 15-18 show four sets of simulation results that prove the capabilities of the LB scalar solver in terms of its numerical accuracy, stability, Galilean invariance, and grid orientation independence, etc. For comparison, results with two different second-order FD schemes are also shown: a van Leer type flux limiter scheme, and a direct mixed scheme (mixture of central and first order upwind schemes).

である。ＴＡは定数である。背景平均流の速度は０．２であり、熱拡散率κは０．００２である。このような低分解能及びκでは、数値安定性及び精度が十分に有効となり得る。背景流を伴わない温度減衰では、ＬＢスカラーソルバ及び有限差分法の両方が理論との優れた一致を示す。非ゼロの背景平均流では、ＬＢスカラーソルバが依然として正確に理論に匹敵することができる。しかしながら、ＦＤ結果は顕著な数値誤差を示す。図１５には、格子時間ステップ８１における温度プロファイルを示す。流束制限ＦＤスキームでは数値拡散がはっきりと見られるが、混合ＦＤスキームでは正確な温度プロファイルもその位置も維持することができない。 A. Shear Wave Damping The first test case is thermal shear wave damping carried by a constant uniform fluid flow. The initial temperature distribution is a uniform distribution plus a spatial sinusoidal variation with a grating wavelength L = 16 and a magnitude of

, where TA is a constant. The background mean flow velocity is 0.2 and the thermal diffusivity κ is 0.002. At such low resolution and κ, numerical stability and accuracy can be sufficiently effective. For temperature decay without background flow, both the LB scalar solver and the finite difference method show excellent agreement with theory. For non-zero background mean flow, the LB scalar solver can still accurately match the theory. However, the FD results show significant numerical errors. Figure 15 shows the temperature profile at grid time step 81. Numerical diffusion is clearly visible in the flux-limited FD scheme, but the mixed FD scheme cannot maintain the exact temperature profile or its location.

B. Inclined Channel with Volumetric Heat Source The second test case is a simulation of the temperature distribution of the channel flow with different lattice orientations. The channel walls are free-slip and the fluid flow remains uniform, resulting in U ₀ =0.2. The channel width is 50 (lattice spacing) and the flow rate Re is 2000. The thermal diffusivity κ is 0.005. The wall temperature is fixed at T _w =1/3. A constant volumetric heat source q=5×10 ^-6 is applied to the bulk fluid domain. The flow has periodic boundary conditions in the other flow directions, which is easy to achieve in the lattice-aligned regime. When the channel (light color) is inclined as shown in Figure 16, the periodic boundary conditions in the other flow directions are again achieved by perfectly matching the inner and outer channel boundaries in the coordinate directions. To demonstrate lattice independence, the inclination angle is chosen to be 26.56 degrees. Similar to the first test case, when the channel is lattice aligned, the temperature distribution using the LB scalar solver and the two FD schemes agree very well with the analytical solution. However, when the channel is tilted, the results from the FD schemes deviate significantly from theory. Figure 16 shows the simulation results of the temperature distribution across the tilted channel. It is clearly shown that the LB results are independent of the lattice orientation. Errors from the FD method also arise from the fundamental difficulty in handling the gradient calculations for tilted boundary orientations, thus introducing additional numerical artifacts. The LB scalar particle advection is accurate in the BC presented here, so that a scalar evolution independent of the lattice orientation can be achieved.

C. Temperature Propagation in an Inclined Channel Due to the absence of neighborhood information in non-grid aligned near-wall regions on a Cartesian grid, it is extremely difficult to accurately estimate the local gradients, which are essential for FD-based methods. Furthermore, the strong dependence on upwind information in FD methods can make the computation of scalar advection even more difficult. In contrast, the boundary treatment of the LB scalar solver achieves accurate local scalar flux conservation as discussed above. Scalar advection in such near-wall regions can be computed accurately. As proof, we implement high temperature convection in a 30 degree inclined channel. With free-slip and adiabatic BC implementation at the solid walls, the fluid velocity is constant U ₀ =0.0909 along the channel. The thermal diffusivity κ is 0.002. Initially, the temperature is 1/3 everywhere except at T=4/9 at the inlet. This temperature front should then be convected undistorted to downstream locations by the uniform background fluid flow. Figure 17 shows the computed temperature front distribution across the channel at grid time step 2000. The temperature front of the LB scalar solver maintains an almost straight shape, while the temperature fronts from the two FD schemes show significant distortion in the near-wall region. It can be said that the LB scalar result shows the thinnest temperature front, which means that the numerical diffusion of the LB scalar solver is reduced.

D. Rayleigh-Benard Natural Convection Rayleigh-Benard (RB) natural convection is a classical benchmark for the accuracy verification of numerical solvers. RB natural convection is a complex physical phenomenon, although the case setup is simple. When the Rayleigh number Ra exceeds a certain critical value, the system transitions from non-flow to convection. The current study, as shown in Figure 18, is carried out under the Boussinesq approximation, where the buoyancy force acting on the fluid is defined as αρg(T-T _m ), where α is the thermal expansion coefficient, g is gravity, and T _m is the average temperature value of the upper and lower boundaries. The most unstable wave number is κ _c = 3.117 when _Ra exceeds the critical value Ra = 1707.762, and this study uses a resolution of 101 × 50. The P _r used here is 0.71. Once the R _B RB convection is established, the heat transfer between the two plates is greatly enhanced. This heat transfer enhancement is expressed by Nu=1+<u _y T>H/κΔT.

Figure 19 shows a screenshot of a fluid simulation. Fluid simulation using the total energy conservation described above (rather than the conventional method) results in a similar representation of the fluid simulation and any conventional corresponding computational data output. However, such a fluid simulation using the method described above may perform the fluid simulation faster and/or using fewer computational resources than other methods when modeling an object, such as an actual physical object.

The equilibrium distribution described here does not contain a pressure term, and therefore its value is always positive during collision calculations, resulting in relatively strong stability. After the collision state is advected as is, an isothermal pressure gradient appears in the momentum, and no pressure convection term appears in the total energy equation. The conservation of the total P-energy is achieved by modifying the resting state, and after advection, the local pressure term is added back to the particle state before computing the moment to result in the correct pressure gradient.

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

Claims (18)

1. A method of simulating a fluid flow on a computer using a simulation that implements total energy conservation, comprising:
simulating fluid activity across the mesh to model particle movement across the mesh;
storing in a computer accessible memory a set of state vectors for each mesh location within said mesh, the state vectors including a plurality of entries, each entry corresponding to a particular one of the possible momentum states at the corresponding mesh location;
the computer modifies the set of state vectors of the particles by adding to the states of the advected particles a value of specific total energy obtained by calculating the specific total energy _Et = E + v2 ^/ 2, where E is the specific internal energy and v is the velocity, and subtracting the value of specific total energy from the states of the non-advected particles over a time interval;
The method according to claim 1, further comprising: After the particles have been advected to subsequent mesh locations according to the modified set of states, the method further comprises adding back a local pressure term to the particle states and then computing moments to yield a correct pressure gradient (∇p).
The method of claim 1. The local pressure terms include a θ-1 term calculated by the computer.
The method of claim 2. simulating the fluid flow activity includes simulating the fluid flow based in part on a first set of discrete grid velocities, the method further including simulating the time evolution of a scalar quantity based in part on a second set of discrete grid velocities that is the same grid speed as the first set of discrete grid velocities or that has a different number of grid speeds than the first set of discrete grid velocities;
The method of claim 2. the computer performing advection of the particles to subsequent mesh locations over a period of time;
The step of simulating the time evolution of the total energy includes:
collecting inflow distributions from neighboring mesh locations for a collision operator and an energy operator;
weighting the inflow distribution;
determining an outflow distribution as a product of the collision operator and an energy operator;
propagating the determined outflow distribution;
The method of claim 1 , comprising: To recover the exact shear stress for any Prandtl number , the energy impingement term is given by,

Here

is the filtered second moment of the non-equilibrium component,
The method of claim 1. applying a zero net surface flux boundary condition such that the inflow distribution is equal to the determined outflow distribution;
The method according to claim 5. The step of simulating the activity includes the step of the computer applying a second distribution function of the equilibrium distribution and the specific total energy _Et , the second distribution function being defined as a particular scalar that is advected with the flow distribution _f ;
The method of claim 1. The second distribution function takes into account the non-equilibrium contribution of f _i to the energy equation to obtain accurate flow behavior over different grid resolutions near the boundary, and the collision operator of the distribution function is given by:

where the terms Ω _f , Ω _q represent the respective collision operators and the equilibrium distribution used in the above equation is:

The method according to claim 8 . Simulating a physical process fluid flow at run -time and providing a computer system
simulating fluid activity across the mesh to model the movement of particles across the mesh;
storing in a computer memory a series of state vectors for each mesh location within the mesh, the state vectors including a plurality of entries, each entry corresponding to a particular one of the possible momentum states at the corresponding mesh location;
modifying said set of state vectors of the advected particles by adding to the states of the advected particles a value of specific total energy obtained by calculating the specific total energy _Et = E + v2 ^/ 2 , where E is the specific internal energy and v is the velocity, and subtracting said value of specific total energy from the states of the non-advected particles over a time interval;
A computer program characterized by causing a computer to perform the following: the value of the specific total energy is added before advection so that pressure is convected, and the value of the added specific total energy is removed from the rest state to correct the value of the added specific total energy;
11. A computer program product as claimed in claim 10 . Total energy is conserved by removing the added specific total energy value , resulting in the correct pressure velocity terms.
12. A computer program product as claimed in claim 11 . The removal of the added specific total energy value is

of

This is done by modifying
13. A computer program product as claimed in claim 12 . 1. A computer system for simulating a physical process fluid flow, comprising:
one or more processor units;
a computer memory coupled to the one or more processor units;
A computer readable medium having instructions stored thereon;
the instructions, when executed, simulate a physical process fluid flow and provide the computer system with:
simulating fluid activity across the mesh to model the movement of particles across the mesh;
storing in a computer memory a series of state vectors for each mesh location within said mesh, said state vectors including a plurality of entries, each entry corresponding to a unique momentum state among the possible momentum states at the corresponding mesh location;
modifying said state vectors of advected particles by adding to the states of the particles a value of specific total energy obtained by calculating the specific total energy _Et = E + ^v2 /2, where E is the specific internal energy and v is the velocity, and subtracting said value of specific total energy from the states of the particles not advected over a time interval;
To make something happen,
1. A computer system comprising: the system further includes instructions for adding back a local pressure term to the particle states after the particles have been advected to subsequent mesh locations according to the modified set of states, and then computing moments to result in a correct pressure gradient (∇p).
The system of claim 14 . the local pressure terms include a θ-1 term calculated by the system;
The system of claim 15 . The instructions for simulating the activity of the fluid flow include:
performing advection of said particles to subsequent mesh locations;
simulating the fluid flow based in part on a first set of discrete grid velocities; and simulating the time evolution of a scalar quantity based in part on a second set of discrete grid velocities that is the same grid speed as the first set of discrete grid velocities or that has a different number of grid speeds than the first set of discrete grid velocities.
The system of claim 15 . and instructions for causing the system to apply a second distribution function of the equilibrium distribution and specific total energy _Et , the second distribution function being defined as a unique scalar that is advected with the flow distribution _f ;
The system of claim 14.

Images