CN111428423A - Lattice boltzmann solver for realizing total energy conservation - Google Patents

Abstract

In addition to a lattice boltzmann (L B) function for the fluid flow, the technique further includes modifying a set of state vectors of the particles by adding a particular total energy to a state of the particles to be advected at intervals of time and subtracting the particular total energy from a state of the particles that are not advected at intervals of time, and performing particle advection according to the modified set of states.

Description

Lattice boltzmann solver for realizing total energy conservation

Priority requirement

Priority of U.S. provisional patent application No. 62/790,528 entitled "L ATTICEBO L TZMANN SO L VER ENFORMCING TOTA L ENERGY CONSERVE SERVICE" (lattice boltzmann solver to achieve total ENERGY conservation), filed 2019, 1, 10.s.c. § 119(e), the entire contents of which are incorporated herein by reference.

Background

The present description relates to computer simulations of physical processes such as physical fluid flow.

The underlying dynamics of lattice Boltzmann systems lie in the fundamental physics of the kinetic theory, which involves the movement of many particles according to the Boltzmann equation there are two fundamental dynamic processes in the basic Boltzmann kinetic system-collisions and advection-collision processes involving interactions between particles that obey the law of conservation and relax to equilibrium-advection processes involving the modeling of the movement of particles from one location to another according to their microscopic velocities.

If a single distribution function is used to solve for all three quantities, the lattice used for L BM should satisfy moments (moment) up to the eighth order.

Disclosure of Invention

According to one aspect, a method for simulating fluid flow on a computer, wherein the simulation achieves total energy conservation, includes simulating activity of a fluid on a grid, the activity of the fluid being simulated to model movement of particles on the grid, storing, in a computer-accessible memory, a set of state vectors for each grid location in the grid, each state vector including a plurality of entries corresponding to a particular (particulate) one of the possible momentum states at the corresponding grid location, calculating, by the computer, a set of energy values for the grid locations in the grid, performing, by the computer, advection of the particles to subsequent grid locations for a time interval, and subtracting, by the computer, the particular total energy value from states of particles that are not advected for the time interval, the state vector of the particle is corrected.

Aspects also include a computer program product, one or more machine-readable hardware storage devices, an apparatus, and a computing system.

Among the features described herein, one or more of the above aspects may include one or more of the following features.

After the particles advection to subsequent grid locations according to the modified state set, this aspect includes adding local pressure terms back to the particle states to provide an appropriate pressure gradient prior to computing moments

The local pressure terms include theta-1 terms calculated by the computer. The act of simulating fluid flow includes simulating fluid flow based in part on the first set of discrete lattice rates. The method also includes simulating a time evolution of the scalar based in part on the second set of discrete lattice rates. The second set of discrete lattice rates is the same as the first set of discrete lattice rates. The second set of discrete lattice rates has a different number of lattice rates than the first set of discrete lattice rates. Simulating a temporal evolution of the total energy includes collecting incoming distributions from adjacent grid locations for the collision and energy operators, weighting the incoming distributions, determining an outgoing distribution that is a product of the collision and energy operators, and propagating the determined outgoing distribution.

Simulating the time evolution of the scalar comprises simulating the time evolution of the scalar based in part on a collision operator that effectively expresses the momentum flux resulting from the product of the collision and energy operators as a first order term. To recover the exact shear stress for any Prandtl (Prandtl) number, the energy collision term is given by:

wherein pi ═ Σ_ic_ic_i(f_i-f_i ^eqA filtered second moment that is an unbalanced component.

The aspect also includes applying a zero net surface flux boundary condition such that the incoming distribution is equal to the determined outgoing distribution. Determining the outgoing distribution includes determining the outgoing distribution to provide a zero surface scalar flux. The scalar comprises a scalar selected from the group consisting of temperature, concentration, and density. Computing a set of energy values for grid locations in a grid includes computing a set of energy values from E_t＝E+v²Total energy E given by/2_tWhere E is energy and v is velocity.

The simulated activity includes the computer's calculation of a specific total energy E_tApplying a balanced distribution and a second distribution function, wherein the second distribution is defined as a distribution f with the flow_iA specific scalar of advection. The second distribution function takes into account f_iAn unbalanced contribution to the energy equation to obtain the correct flow behavior near the boundary and across different grid resolutions, wherein the collision operator for the distribution function is given by:

wherein the term omega_f,Ω_qRepresenting the respective collision operator; the equilibrium distribution used in the above equation is:

total energy this term is added prior to advection so that pressure is convected and energy is removed from rest in order to compensate for the added energy. Removing the added energy conserves total energy and provides the correct pressure velocity term. By modifying the stop state

So that

Can be used forProviding for the removal of the added energy.

The techniques described herein for simulating fluid flow use the lattice Boltzmann (L B) method and scalar solver transfer equations.A L BM solver for compressible flows with total energy conservation and higher stability ranges is provided that can be used for real world applications.

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

Drawings

FIG. 1 depicts a system for fluid flow simulation.

Fig. 2-3 depict flow charts illustrating operations for complex fluid flow simulations and total energy calculations.

Fig. 4-6 are illustrative comparisons of predicted results based on calculated expectations.

Fig. 7 and 8 illustrate the velocity components of two L BM models (prior art).

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

Fig. 10 is a perspective view of a micro-block (prior art).

Fig. 11A and 11B are illustrations of lattice structures used by the system of fig. 1 (prior art).

Fig. 12 and 13 illustrate a variable resolution technique (prior art).

Fig. 14 illustrates a parallelepiped that helps to understand the interaction of particles with bins.

Figure 15 illustrates the movement of a particle from a voxel to a surface (prior art).

Fig. 16 illustrates the movement of particles from surface to surface (prior art).

FIG. 17 is a flow chart of a process for performing surface dynamics (prior art).

FIG. 18 is a process for generating a distribution function (prior art flow chart) with thermodynamic steps unrelated to the particle collision step.

Detailed Description

In the process discussed in fig. 9 below, a flow simulation process using a lattice boltzmann solver that achieves total energy conservation is described.

In FIGS. 7, 8, and 10-18, these figures are labeled as "Prior Art". The figures are labeled as Prior Art because they typically appear in the patents referenced above.

A. Method for solving scalar quantity

In the systems and methods described herein, an L BM-based physical process simulation system is based on modeling of scalars (as opposed to vectors) coupled with modeling of fluid flow.

The details of the technique of Coupling scalar modeling to simulated fluid flow are described in detail in co-pending U.S. patent publication US-2016 0188768-A1 entitled "Temperature Coupling Algorithm for Hybrid Thermal Method L" by Pradeep Gopalakrishnan et al, the entire contents of which are incorporated herein by reference.

For example, the system may be used to determine a convective temperature distribution within the system. For example, if the system (formed by the volume represented by the plurality of voxels) includes a heat source and there is an airflow within the system, some regions of the system will be hotter than others depending on the airflow and proximity to the heat source. To model this situation, the temperature distribution within the system may be represented as a scalar, where each voxel has an associated temperature.

In another example, the system may be used to determine chemical distribution within the system. For example, if the system (formed by the volume represented by the voxels) includes a source of contaminants, such as dirty bombs or chemicals or other particulates suspended in air or liquid, and there is a flow of air or liquid within the system, certain regions of the system will have a higher concentration than other regions based on the flow and proximity to the source. To model this situation, the chemical distribution within the system may be represented as a scalar, with each voxel having an associated concentration.

In some applications, multiple different scalars may be simulated simultaneously. For example, the system may simulate both temperature and concentration distributions in the system.

For example, the lattice boltzmann (L B) method for solving scalar transfer equations may be used to indirectly solve scalar transfers.

The method assigns a vector to each voxel in a volume to represent fluid flow and a scalar to each voxel in the volume to represent a desired scalar variable (e.g., temperature, density, concentration, etc.).

Unlike methods based on discretized macroscopical equations, L BM begins with "mesoscopic" Boltzmann dynamics equations to predict macroscopic fluid dynamics, the resulting compressible and unstable solution method can be used to predict various complex flow physics, such as aeroacoustics and purely acoustic problems.

Modeling simulation space

In an L BM-based physical process simulation system, fluid flow may be controlled by a set of dispersion velocities c_iDistribution function value f to be evaluated_iAnd (4) showing. The kinetics of the distribution function is governed by the equation (I.1), where f_i(0) Known as the equilibrium distribution function, is defined as:

wherein,

this equation describes the distribution function f_iThe well-known lattice-glass equation of time evolution. The left hand side shows the change in distribution due to the so-called "fluidization process". The fluidization process is a block of fluid that starts at a grid location and then moves to the next grid location along one of the velocity vectors. At that point, a "collision operator" is calculated, i.e., the effect of nearby fluid blocks on the originating fluid block. The fluid can only move to another grid location, so proper selection of velocity vectors is necessary so that all components of all velocities areMultiple of common rate (common speed).

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

The BGK operator is constructed from the physical argument that, regardless of the specifics of the collision, the distribution function is approximated by the collision

{f^eq(x, v, t) } given a well-defined local balance:

where the parameter τ represents the characteristic relaxation time to reach equilibrium via collision. When dealing with particles, such as atoms or molecules, the relaxation time is usually considered as a constant.

From this simulation, conventional fluid variables, such as mass ρ and fluid velocity u, are obtained as a simple sum in equation (I3):

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

Because of symmetry considerations, the set of velocity values is selected in such a way that they form a certain lattice structure when spanned in the configuration space the dynamics of this discrete system follow L BE having the form:

f_i(x+c_i,t+1)-f_i(x,t)＝C_i(x,t)

wherein the collision operator typically takes the form of BGK as described aboveForm (a). By appropriate choice of the form of the equilibrium distribution, it can be shown theoretically that the lattice boltzmann equation yields the correct fluid dynamics and thermo-fluid dynamics. I.e. from f_i(x, t) derived hydrodynamic moments obey the Navier-Stokes equation in macroscopic constraints. These moments are defined by equation (I3) above.

c_iAnd w_iThe L BM model can be efficiently implemented on scalable computer platforms and operates with maximum robustness to time-unstable flows and complex boundary conditions.

The standard technique for obtaining macroscopic equations for the motion of fluid systems from the boltzmann equation is the Chapman-Enskog method, in which successive approximations of the full boltzmann equation are taken in the fluid system small perturbations in density travel at the speed of sound in gas systems acoustic velocity is generally determined by temperature.

Referring to FIG. 1, a system 10 for simulating fluid flow (e.g., near a representation of a physical object) is shown. The system 10 in this implementation is based on a client-server architecture and includes a server system 12 and a client system 14 implemented as a massively parallel computing system 12. The server system 12 includes a memory 18, a bus system 11, an interface 20 (e.g., a user interface/network interface/display or monitor interface, etc.), and a processing device 24. Within memory 18 are a grid preparation engine 32 and a simulation engine 34. The system 10 accesses a data store 38 that stores 2D and/or 3D grids, coordinate systems, and libraries.

Although FIG. 1 shows grid preparation engine 32 in memory 18, the grid preparation engine may be a third party application executing on a different system than server 12. Whether grid preparation engine 32 executes in memory 18 or on a system other than server 12, grid preparation engine 32 receives user-provided grid definitions 30, and grid preparation engine 32 prepares the grid and sends the prepared grid to simulation engine 34. Simulation engine 34 includes a particle collision interaction module, a particle boundary modeling module, and an advection module that performs advection operations as described below. The simulation engine 34 also includes a total energy solver, as discussed below in FIG. 3.

Referring now to FIG. 2, a process 40 for simulating fluid flow proximate a representation of a physical object is shown, in examples that will be described herein, the physical object is an airfoil shaped body (airfoil), however, the use of airfoil shaped bodies is merely illustrative as the physical object may be any shape, and in particular may have a flat and/or curved surface(s). The process receives a grid for the physical object being simulated, e.g., from a client system 14 or by retrieving from a data repository 42. in other embodiments, an external system or server 12 generates a grid for the physical object being simulated based on user input.A process pre-computes 44 geometry from the retrieved grid and performs 46 dynamic lattice Boltzmann model simulation using the pre-computed geometry corresponding to the retrieved grid

The set of state vectors is modified by adding a specific total energy value to the state of the advected particle and subtracting the specific total energy value from the state of the particle that is not advected within the time interval.

Total energy conservation compressible flow solver

When performing complex fluid flow simulations, it may be beneficial to solve for the scalar (e.g., temperature distribution, concentration distribution, and/or density) simultaneously in conjunction with solving for the fluid flow. A particular scalar is distinguished from a scalar distribution, the scalar distribution being the difference between q (equation 10) and h (equation 2), and the particular scalar q is used as E _ t, and h is the density × E _ t. The method discussed is a scalar placement on top of the convected stream as described in the above-mentioned U.S. patent publication US-2016-0188768-a 1.

Processing of pressure terms in total energy conservation

Referring now to FIG. 3, a compressible L BM solver with a total energy conservation process 50 is now described, the compressible L BM solver having a better stability range than existing methods and being suitable for real world applications.

Wherein the total energy E_tFrom E_t＝E+v²Given by/2, ρ is density, v is velocity, p is pressure, k is thermal conductivity, T is temperature, and τ is shear stress. This is a common equation that all computational fluid dynamics tools for high velocity flows need to solve.

One of the complexities involved with the total energy solver is the addition of the pressure-convection term "p" in equation 1. like all computational fluid dynamics tools for high velocity flows, the lattice Boltzmann method (L BM) must solve this equation for total energy conservation.

The pressure convection problem of the L BM method can be explained by the fact that normal scalar transport has the following equation

The left hand side of the above equation includes a convection term, while the right hand side of the equation has a diffusion term. Unlike the total energy equation, all the terms in the above equation contain only one scalar variable S, which makes advection simpler. The process 50 simulates 52 the evolution of the particle distribution according to the selected form of the equilibrium distribution.

In addition to the equilibrium distribution (selected form of equation 1) 52, the use 54 for calculating the bits is also usedDetermining the total energy E_tProcess 50 includes advection 56 of the particles in the L BM to the next unit q. meeting energy conservation is achieved by the product of these two distributions_iA specific scalar of advection. Unlike prior methods that use a second distribution function independently to solve for energy, the second distribution function takes into account f_iUnbalanced contribution to the energy equations to obtain correct flow behavior near the boundary and across different grid resolutions.

The collision operator of the distribution function is given by:

f_i′(x，t)＝f_i(x，t)+Ω_f[f_i(x，t)，f_i ^eq(x，t)]equation 2

The term Ω in equations 2 and 3_f,Ω_qAre used to represent the respective collision operators. The collision operator Ω in equation 5 is discussed extensively in U.S. patent No. 9,576,087_fThe entire contents of this U.S. patent are incorporated herein by reference. The equilibrium distribution used in the above equation is

The observation that is apparent in the above-described distribution function equation 7 for a total energy solver is that the equilibrium distribution does not involve the pressure term (p) and the θ -1 term in equation 1, i.e., p ═ ρ RT, which is critical for heat flow. The second distribution is a specific total energy, which is the internal energy E ═ C_νT and kinetic energy v²The sum of/2.

Thus, the distribution function will always provide a positive value during a crash operation, and thus the distribution function provides a solver with a correspondingly high stability. In equation (2) the post-collision state is advected as is, the distribution function of equation 4 will only result in an isothermal pressure gradient in momentum, and no pressure convection term in the total energy equation. The pressure gradient and pressure convection terms are due to advection, while collisions do not affect these terms, in contrast to other methods that use a balanced distribution and include pressure terms during both advection and collisions, as shown below.

The term θ ═ RT/T₀The above balance is made negative corresponding to the pressure, and thus the prior art is highly unstable. It should also be noted that the second distribution function

Is used independently for the total energy, which results in noise across different grid resolutions.

The process 50 modifies 58 the advection state by adding a certain total energy to the advected state while removing the certain total energy from the stopped state (non-moving state).

More specifically, in this modification 58, pressure introduction is only performed during the advection step in order to recover the correct momentum and energy equations. Correcting motion states at grid locations prior to advection processes

The value, as shown in equation 8 below:

however, this term (the second equation in equation 5) is added to q_iAdditional energy is also added to the motion state. To compensate for the increased energy represented in this term, the same amount of energy is removed from the stopped state (the state where the particles do not convect). By doing so, the process saves energy and provides the correct pressure velocity term. In particular, by modifying the stop state

So that

To provide total energy conservation.

After advection, the local pressure terms are added back to the particle state before the moments are calculated, thereby creating the appropriate pressure gradient

The modifications before advection (equation 8) and after advection (equation 6) restore the correct equation of state (or correct pressure gradient) for momentum

) And also recover the correct pressure convection term for the total energy

Since the pressure term (equation 9) is removed, the state becomes positive during the collision and the stability range of the total energy solver is improved.

The advection is corrected by applying a total energy correction to the advection in which the computer corrects the state vector of the particles by adding a specific total energy value to the state of the particles of the advection and subtracting the specific total energy value from the state of the particles that are not advected within the time interval, as discussed in fig. 9 below.

Energy regularized (regularized) collision operator

The collision operator Ω for flow states is discussed in the above-mentioned us patent 9,576,087_f. The collision operator Ω for energy is discussed below_q. A regular collision with a first moment is given by:

equation 7 is the filtered unbalanced contribution to the total energy. The conservation of energy is satisfied by the multiplication of these two states, which results in the collision term, the second term on the right hand side of equation (7)

Comprises the following steps:

the zero moment of the above term is zero and saves energy. In the Chapman-Enskog expansion process, the first moment of the collision term involves the derivation of thermal diffusion and work done by shear stress. As a first order approximation, the collision term can be expressed as:

in the above equation, the momentum flux term ρ ν is the multiplication of two states

The result of (1). This term ρ ν causes numerical instability for high velocity flows and also mach number dependent diffusion. To mitigate these effects, the term ρ ν is removed by defining the collision operator as:

the above expression omits the momentum flux term ρ ν and improves the numerical stability for high mach flow.

For L BM, the speed values will generally begin to be less than "1" (one) and therefore for lower speed flows, these error terms will be much lower and can be ignored, however, for higher speed flows, the speed values will be up to 1 order and therefore the error term ρ ν is relatively higher and removed by equation 10.

Viscous effect of arbitrary prandtl number

Relaxation time of τ used in energy solver_qThis will result in a viscosity equal to the thermal diffusivity (τ)_qThe shear stress term of-0.5 ρ RT) performs incorrect work. To recover the exact shear stress for any prandtl number, the following additional terms are added to the energy collision term.

Where pi ═ Σ_ic_ic_i(f_i-f_i ^eq) Is the second moment of the filtered unbalanced component.

Total energy solver capability

Total energy solvers can be included, for example

And the total energy solver can be used for a wide range of industrial applications, including those involvingAnd applications for high subsonic and transonic flows.

Several benchmarks are discussed in this section to show some of the potential advantages of the disclosed total energy solver for compressible streams, particularly in terms of energy conservation. The simulation results were compared to the analytical solution and the typical finite difference hybrid solver solution based on PDE.

FIG. 4 shows a graph of the total temperature of the flow above the plate at Mach velocity (Ma) for Ma 0.7 and x 1.5m from the leading edge, as shown in FIG. 4, the total energy solver (solid line 62) for the flow above the plate accurately recovers the total temperature, whereas the hybrid solver (L BM flow solver and finite difference based entropy solver) (dashed line 64) cannot maintain the total temperature.

Fig. 5 shows a graph of the static temperature of the plate at mach rate (Ma) for Ma 0.7 and x 1.5m from the leading edge. The graph is relative to as u/u_∞U (y-axis) is plotted as T/T_∞Temperature t (x-axis) the total energy solver of the present disclosure, represented by solid line 66, is compared to a hybrid solver (L BM flow solver and finite difference-based entropy solver) (dashed line 68) the figure illustrates that for Ma 0.7, the static temperature of the plate at mach rate (Ma) at x 1.5m from the leading edge matches the wolstz equation (asterisk line 70).

With reference to FIG. 6, the back pressure P is shown_b＝0:75P_tGraph of results of temperature simulation of a convergence-divergence verification (CDV) nozzle. The graph plots temperature T against temperature T on the x (y-axis)_∞(x-axis) the total energy solver of the present disclosure, represented by the solid line 72, is compared to a hybrid solver (L BM flow solver and a finite difference-based entropy solver) (dashed line 74.) the total energy solver of the present disclosure has good consistency with an analytically determined dotted line 76.

A system with the disclosed total energy solver may be used to determine the convective temperature distribution within the system. For example, if the system (formed by a volume represented by a plurality of voxels) includes a heat source and there is an airflow within the system, some regions of the system will be hotter than others depending on the airflow and proximity to the heat source. To model such a system, the temperature distribution within the system may be represented as a scalar, where each voxel has an associated temperature.

In another example, the system may be used to determine chemical distribution within the system. For example, if the system (formed by the volume represented by the voxels) includes a source of contaminants, such as dirty bombs or chemicals or other particulates suspended in air or liquid, and there is a flow of air or liquid within the system, certain regions of the system will have a higher concentration than other regions based on the flow and proximity to the source. To model this situation, the chemical distribution within the system may be represented as a scalar, with each voxel having an associated concentration.

In some applications, multiple different scalars may be simulated simultaneously. For example, the system may simulate both temperature and concentration distributions in the system.

For example, the lattice boltzmann (L B) method for solving scalar transfer equations may be used to indirectly solve scalar transfers.

In such an arrangement simulation, a second set of distribution functions is introduced for delivering scalars in addition to the lattice boltzmann function for fluid flow.

Unlike methods based on discretized macroscopical equations, L BM begins with "mesoscopic" Boltzmann dynamics equations to predict macroscopic fluid dynamics, the resulting compressible and unstable solution method can be used to predict various complex flow physics, such as aeroacoustics and purely acoustic problems.

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

As also shown in FIG. 8, a second model (3D-1)200 is shown, comprising a three-dimensional model of 39 speeds, where each speed is represented by one of the arrows of FIG. 8. Of these 39 velocities, one velocity represents a particle that is not moving; three sets of six velocities represent particles moving at a normalized rate (r), twice the normalized rate (2r), or three times the normalized rate (3r) in either a positive or negative direction along the x, y, or z axis of the lattice; eight velocities represent particles moving at a normalized rate (r) with respect to all three of the x, y, z lattice axes; and twelve velocities represent particles moving at twice the normalized rate (2r) with respect to two of the x, y, z lattice axes.

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

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

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

The L BM model described above provides a specific class of efficient and robust discrete velocity dynamical models for numerical simulation of flows in two and three dimensions.

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

The Reynolds number relates to the viscosity (v) of the flow, the characteristic length (L) of the object in the flow, and the characteristic velocity (u) of the flow:

re u L/v equation (I4)

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

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

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

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

c_i＝(c_ix,c_iy,c_izequation (I5)

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

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

Each voxel (i.e., each lattice site) is represented by a state vector f (x). The state vector fully defines the state of a voxel and includes 39 entries. These 39 entries correspond to one energy zero state, 6 energy one states, 8 energy three states, 6 energy four states, 12 energy eight states, and 6 energy nine states. Using this velocity set, the system can generate maxwell-boltzmann statistics on the implemented equilibrium state vectors. For total energy, q for second allocation_i(x, t) the same number of terms (39 states) are used to represent a particular total energy. Along with the state vector f (x), this is used to solve for total energy conservation.

For processing efficiency, voxels are grouped in volumes of 2x2x2 called microscopic blocks. The micro-tiles are organized to allow parallel processing of voxels and to minimize overhead associated with the data structure. Shorthand notation for voxels in microscopic blocks is defined asN_i(n), where n represents the relative position of lattice sites in the micro-slabs and n {0,1,2,. 7 }. The microscopic blocks are shown in fig. 10.

Referring to fig. 11A and 11B, the surface S (fig. 11A) is represented as a bin F in the simulation space (fig. 11B)_αSet of (2):

S＝{F_αequation (I6)

A bin is not limited to a boundary of a voxel, but generally has a size on the order of the size of voxels adjacent to the bin, or slightly smaller than the size of voxels adjacent to the bin, such that the bin affects a relatively small number of voxels_αWith unit normal (n)_α) Surface area (A)_α) Center position (x)_α) And a bin distribution function (f) describing the surface dynamic properties of the bins_i(α)) total energy distribution function q_i(α) is processed in the same way as voxel interactions and flow distribution for face-elements.

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

C. Identifying voxels affected by bins (faces)

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

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

Referring to FIG. 14, for a vector c having a non-zero velocity_iEach state i, bin F of_αFrom parallelepiped G_iαA defined region for receiving or transmitting particles thereto, wherein the parallelepiped G_iαWith velocity vector c_iSum bin (| c)_in_i|) unit normal n_αHeight defined by the magnitude of the vector dot product of (a) and surface area of the bin (a)_αBase defined so that it is parallelepipedal G_iαVolume V of_iαEqual to:

V_ia＝|c_in_a|A_aequation (I7)

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

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

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

_ia(x)＝N_i(x)+V_ia(x) Equation (I8)

When parallelepipedon G_iαWhen intersected by one or more bins, the following condition is true:

V_ia＝∑V_a(x)+∑V_ia(β) equation (I9)

Wherein the first summation takes into account_iαAll voxels overlapping and the second term taking into account G_iαAll bins that intersect. When parallelepipedon G_iαWhen not intersected by another bin, this expression reduces to:

V_ia＝∑V_ia(x) Equation (I10)

D. Performing simulations

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

1. Boundary conditions for surfaces

In order to correctly model the interaction with the surface, four boundary conditions are met for each bin. First, the combined mass of particles received by a bin is equal to the combined mass of particles transmitted by the bin (i.e., the net mass flux to the bin is equal to zero). Second, the combined energy of the particles received by a bin is equal to the combined energy of the particles transmitted by the bin (i.e., the net energy flux to the bin is equal to zero). These two conditions may be satisfied by requiring the net mass flux at each energy level (i.e., energy level one and energy level two) to be equal to zero.

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

2. Clustering from voxels to bins

As a first step in simulating the interaction between the particles and the surface, the particles are aggregated from voxels and provided to bins (step 308). As indicated above, voxel N (x) and bin F_αThe flux of the particles in state i in between is:

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

In this view, for the directional bin F_αEach state i (c)_in_α<0) Is supplied from a body elementSurface element F_αThe number of particles of (a) is:

_iαV→F＝Σ_{X iα}(x)＝Σ_XN_i(x)V_iα(x) Equation (I12)

Only its V_iα(x) Voxels with non-zero values need to be summed. As indicated above, the size of the bins is chosen such that V_iα(x) Only for a few voxels have non-zero values. Because of V_iα(x) And P_f(x) May have non-integer values, therefore_α(x) Stored and processed as real numbers.

3. Bin to bin movement

Next, the particles are moved between bins (step 310). If used in bin F_α(ii) an incoming state (c)_in_α<0) Of the parallelepiped G_iαIs covered by another surface element F_βIntersect, then by F_αPart of the received state i particles will come from bin F_β. In particular, bin F_αReception by bin F during previous time increment_βPart of the resulting state i particles. This relationship is shown in FIG. 17, where the bin F is binned_βIntersecting parallelepipeds G_iαIs equal to quilt element F_αIntersecting parallelepipeds G_iβPortion 1005. As noted above, the intersecting portion is denoted V_iα(β) Using this term, bin F_βAnd surface element F_αThe flux of state i particles in between can be described as:

_iα(β，t-1)＝_i(β)V_iα(β)/V_iαequation (I.13)

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

_iαF→F＝Σ_{β iα}(β)＝Σ_{β i}(β，t-1)V_iα(β)/V_iαequation (I.14)

And the total flux to state i particles in the bin is:

_iIN(α）＝_iαF→F+_iαF→F＝Σ_xN_i(x)V_iα+Σ_{β i}(β，t-1)V_iα(β)/V_iαequation (I.15)

The state vector N (α) (also known as the bin distribution function) for a bin has M entries corresponding to the M entries of the voxel state vector M is the number of discrete lattice rates the input states of the bin distribution function N (α) is set equal to the flux of particles into those states divided by the volume V_iα：

For c_in_α<0，N_i(α)＝_iIN(α)/V_iαEquation (I.16)

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

for c_in_α≥0，N_i(α)＝_iOTHER(α)/V equation (I.17)

Wherein_iOTHER(α) is a method for generating using the above_iIN(α), but applying the technique to divide incoming states (c)_in_α<0) Out of state (c)_in_αNot less than 0). In an alternative method of the invention,_iOTHER(α) may utilize information from a previous time step_iOTHERA value of (α) is generated such that:

_iOTHER(α,t)＝_iOUT(α, t-1) equation (I.18)

For parallel state (c)_in_α＝0)，V_iαAnd V_iα(x) Are all zero. In the use for N_i(α) wherein V is_iα(x) Present in the molecule (according to the application)_iOTHERExpression of (α)) And V is_iαAppear in the denominator (according to the use for N_iExpression of (α). thus, when V_iαAnd V_iα(x) Near zero, N for parallel state_i(α) is determined to be N_iThe values for the states with zero speed, i.e., the rest state and states (0,0,0,2) and (0,0,0, -2), are initialized at the beginning of the simulation based on initial conditions for temperature and pressure.

4. Executive surface element surface dynamics

Next, surface dynamics are performed on each bin to satisfy the four boundary conditions discussed above (step 312). A process for performing surface dynamics for the bins is shown in FIG. 18. Initially, the bin F is determined by determining the combined momentum P (α) of the particles at the bin_αCombined momentum of (step 1105):

for all of the i's, the average,

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

P_n(α)＝n_αequation P (α) (I.20)

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

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

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

Δ_i(α)＝_iIN(α)-N_n-βi(α)V_iαequation (I.21)

Using this difference, the outgoing flux distribution is:

for n_αc_i>0，_iOUT(α)＝N_n-βi(α)V_iα-.Δ._i*(α) equation (I.22)

And wherein i is a state having an opposite direction to state i. For example, if state i is (1, 1,0,0), then state i is (-1, -1,0, 0). To account for skin friction and other factors, the outgoing flux distribution can be further refined as:

for n_αc_i>0，

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

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

For_iOUTThe first and second terms implement normal momentum flux boundary conditions to the extent that collisions have effectively produced a Boltzmann distribution but include tangential momentum flux anomalies the fourth and fifth terms correct for such anomalies, which may be due to dispersion effects or non-Boltzmann structures due to insufficient collisions_fThe generation of (c) is described below. It should be noted that all items relating to vector manipulation are possibleGeometry factors calculated before starting the simulation.

Thus, the tangential velocity is determined as:

u_i(α)＝(P(α)-P_n(α)n_α) /[ rho ], equation (I.25)

Where ρ is the density of the bin distribution:

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

Δ_i(α)＝_iIN(α)-N_n-βi(α)V_iαequation (I.27)

Then, the outgoing flux distribution becomes:

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

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

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

wherein p is_αIs in the surface element F_αAnd based on the average density and temperature values of the volume elements providing particles to the surface element, and u_αIs the average velocity at that bin.

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

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

for c_jin_α>0，

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

for c_jin_α>0，_αjiOUTf＝_αjiOUT+_αjiEquation (I.32)

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

5. Moving from voxel to voxel

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

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

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

The moving operation becomes slightly more complex for voxels that interact with one or more surfaces. This may result in one or more fractional particles being transferred to the bin. This transfer of fractional particles to the voxel results in the fractional particles remaining in the voxel. These fractional particles are transferred into voxels that are binned.

Also shown is a total energy correction to advection 315, where the computer corrects the state vector of the particle by adding a specific total energy value to the state of the advected particle and subtracting the specific total energy value from the state of the particle that is not advected over a time interval.

Referring to FIG. 15, when a portion 900 of the state i-particle of a voxel 905 is moved to a bin 910 (step 308), the remaining portion 915 is moved to the voxel 920 in which the bin 910 is located and from which the state i-particle points to the bin 910. Thus, if the state population is equal to 25 and V_iα(x) Equal to 0.25 (i.e. one quarter of a voxel and parallelepiped G)_iαIntersect), then 6.25 particles will be moved to bin F_αAnd 18.75 particles will be moved to the quilt element F_αOccupied voxels. Since multiple bins will intersect a single voxel, the number of state i particles delivered to a voxel N (f) occupied by one or more bins is:

where N (x) is the source voxel.

6. Dispersion from bins to voxels

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

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

after dispersing the particles from bin to voxel, the particles are combined with the particles that have been advected from the surrounding voxels, and the result is integer-ized, with the possibility that some directions in some voxels may underflow (become negative) or overflow (in an eight-bit implementation, over 255). This will result in an increase or loss of mass, momentum and energy after they have been truncated to fit within the allowed value range.

To prevent this from happening, the out-of-bounds mass, momentum, and energy are accumulated prior to truncation of the violation state. For the energy to which the state belongs, an amount of quality equal to the value of gain (due to underflow) or loss (due to overflow) is added back to the randomly (or sequentially) selected state that has the same energy and is not itself subject to overflow or underflow. The additional momentum resulting from this addition of mass and energy is added up and added to the momentum from truncation. By adding mass only to the same energy state, both mass and energy are corrected when the mass counter reaches zero. Finally, the momentum is corrected using a push/pull technique until the momentum accumulator returns to zero.

7. Performing fluid dynamics

Fluid dynamics are performed (step 318) figure 9. This step may be referred to as micro-dynamics or intra-voxel manipulation. Similarly, the advection process may be referred to as inter-voxel operation. The micro-kinetic operations described below can also be used to collide particles in a bin to produce a boltzmann distribution.

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

E. Variable resolution

Referring to fig. 12, variable resolution (as discussed above) employs different sized voxels, hereinafter referred to as coarse voxels 12000 and fine voxels 1205. (the following discussion refers to voxels having two different sizes; it should be appreciated that the described techniques can be applied to voxels of three or more different sizes to provide higher levels of resolution). The interface between the coarse voxels and the regions of the fine voxels is referred to as a Variable Resolution (VR) interface 1210.

When using variable resolution at or near the surface, bins may interact with voxels on both sides of the VR interface. These bins are classified as VR interface bin 1215 (F)_αIC) Or VR Fine primitive 1220 (F)_αIF). VR interface bin 1215 is a bin located on the coarse side of the VR interface and has a coarse parallelepiped 1225 that extends into the fine voxel. (coarse parallelepiped is determined by the dimensions of the coarse voxels_iAnd the fine parallelepiped determines c from the dimensions of the fine voxel_iDimension (d). VR fine bin 1220 is a bin positioned on the fine side of the VR interface and has a fine parallel hexahedron 1230 that extends into the coarse voxel. Interface bin related processing may also involve coarse bin 1235 (F)_αC) And fine bin 1240 (F)_αF) The interaction of (a).

For both types of VR bins, surface dynamics were performed at a fine scale and operated as described above. However, VR bins differ from other bins in the way particles are advected to and from the VR bins.

The interaction with the VR bins is processed using the variable resolution process 1300 shown in fig. 13. Most of the steps of this process are performed using equivalent steps discussed above for interaction with non-VR bins. Process 1300 is performed during a coarse time-step (i.e., a time period corresponding to a coarse voxel) that includes two phases, each phase corresponding to a fine time-step. Bin surface dynamics are performed during each fine time step. Therefore, VR interface bin F_αICConsidered as two fine bins of equal size and orientation, called black bins F respectively_αICbAnd red bin F_αIC. Black surface element F_αICbAssociated with the first fine time step within the coarse time step, and the red bin F_αICrAssociated with a second fine time step within the coarse time step.

Initially, particles move (advect) between bins through a first surface-to-surface advection phase (step 1302). The particles are weighted by a weighting factor V_-αβFrom black surface element F_αICbMoving to rough surface element F_βCThe weight factor corresponding to the sub-bin F_αExtend out and are positioned in the surface element F_βVolume of the unobstructed part of the rear rough parallelepiped (fig. 12, 1225) minus the secondary bin F_αExtend out and are positioned in the surface element F_βThe rear fine parallelepiped unobstructed part (fig. 12, 1245). C of fine voxel_iIs the magnitude of the coarse voxel c_iHalf the magnitude of (c). As described above for bin F_αThe volume of the parallelepiped defined as:

V_iα＝|c_in_α|A_αequation (I.36)

Thus, due to the surface area A of the voxel_αDoes not vary between coarse and fine parallelepipeds and is due to the unit normal n_αIs always 1, and thus corresponds to a face elementThe volume of the fine hexahedron of (a) is half of the volume of the corresponding coarse parallelepiped of that bin.

The particles are weighted by a weighting factor V_αβFrom coarse surface elements F_αCMove to black bin F_βICThe weight factor V_αβCorresponding to the slave bin F_αExtend out and are positioned in the surface element F_βVolume of the unobstructed portion of the latter fine parallelepiped.

The particles are weighted by a weighting factor V_αβFrom red bin F_αICrMoving to a coarse bin F_βCAnd with a weighting factor V_-αβFrom coarse surface elements F_αCShift to red bin F_βICr。

The particles are weighted by a weighting factor V_αβFrom red bin F_αICrMove to black bin F_βICb. At this stage, no black to red advection occurs. Furthermore, since the black and red bins represent consecutive time steps, black-to-black smoothing (or red-to-red smoothing) never occurs. For similar reasons, the particles in this phase are weighted by a factor V_αβFrom red bin F_αICrMove to fine flour element F_βIFOr F_βFAnd from the fine bin F with the same weighting factor_αIFOr F_αFMove to black bin F_αICb。

Finally, the particles are weighted by the same weight factor from the fine bin F_αIFOr F_αFMove to other bins F_βIFOr F_βFThen by a weighting factor V_CαβRough surface element F_αCMove to other rough surface elements F_CThe weight factor V of_CαβCorresponding to the slave bin F_αExtend out and are positioned in the surface element F_βVolume of the latter unblocked portion of the coarse parallelepiped.

After the particles advection between the surfaces, the particles are aggregated from the voxels in a first aggregation stage (step 1304-. For fine bin F_αFThe particles are collected from the fine voxels using fine parallelepipeds (step 1304), and for coarse bins F_αCCollecting particles from coarse voxels using coarse parallelepipeds (step)1306). Then for black bin F_αIRbAnd for VR fine bin F_αIFThe fine parallelepiped is used to cluster particles from both coarse and fine voxels (step 1308). Finally, the difference between coarse and fine parallelepipeds is used for the red bin F_αIRrThe particles are aggregated from the coarse voxels (step 1310).

Next, coarse voxels that interact with fine voxels or VR bins are decomposed into a set of fine voxels (step 1312). The state of the coarse voxels that will transport particles to the fine voxels in a single coarse time step is resolved. For example, the appropriate state of a coarse voxel that is not intersected by a voxel is decomposed into eight fine voxels that are oriented as the microblocks of FIG. 4. The appropriate state of coarse voxels intersected by one or more bins is decomposed into a set of full and/or partial fine voxels corresponding to portions of coarse voxels that are not intersected by any bin. Particle density N of coarse voxels and fine voxels resulting from decomposition thereof_i(x) Fractional factors P of equal, but fine voxels_fMay be different from the fractional factor for coarse voxels, or may be different from the fractional factor for other voxels.

Thereafter, the fine bin F is processed_αIFAnd F_αFSurface dynamics are performed (step 1314) and black bin F is binned_αICbSurface dynamics are performed (step 1316). The kinetics were performed using the process shown in fig. 11 and discussed above.

Next, the particle is moved between fine voxels including actual fine voxels and fine voxels decomposed from coarse voxels (step 1318). Once the particles have moved, the particles are removed from the fine bin F_αIFAnd F_αFDispersed to the fine voxels (step 1320).

The particles also extending from black surface element F_αICbDispersed to fine voxels (including fine voxels decomposed from coarse voxels) (step 1322). If a voxel receives a particle without a surface present at the time, the particle will disperse into a fine voxel. In particular, when the voxel is a real fine voxel (as opposed to a fine voxel decomposed from coarse voxels), when the hyper-voxel is a real oneGenerating voxel N (x + c) of 1 speed unit of voxel N (x)_i) For actual fine voxels, or when a voxel N (x) of 1 speed unit is exceeded (x + c)_i) For fine voxels derived from coarse voxel decomposition, the particles are dispersed into voxel N (x).

Finally, the first fine time step is completed by performing fluid dynamics on the fine voxels (step 1324). The voxel for which fluid dynamics are performed does not include fine voxels resulting from coarse voxel decomposition (step 1312).

Process 1300 performs similar steps during the second fine time step. Initially, in a second surface-to-surface advection stage, particles are moved between the surfaces (step 1326). The particles flow flat from the black bin to the red bin, from the black bin to the fine bin, from the fine bin to the red bin, and from the fine bin to the fine bin.

After the particles are advected between the surfaces, the particles are aggregated from the voxels in a second aggregation stage (step 1328-1330). For red bin F_αIRrThe particles are aggregated from the fine voxels using a fine parallelepiped (step 1328). Also for fine bins F_αFAnd F_αIFThe particles are aggregated from the fine voxels using the fine parallelepiped (step 1330).

Thereafter, as described above, for the fine bin F_αIFAnd F_αF(step 1332) for coarse bin F_αC(step 1134), and for the red bin F_αICr(step 1336) performing surface dynamics.

Next, the particle is moved between voxels using the fine resolution (step 1338), such that the particle moves back and forth between the fine voxel and the fine voxel representing the coarse voxel. The particles are then moved between voxels using the coarse resolution (step 1340) so that the particles move back and forth between the coarse voxels.

Next, in the combining step, particles are dispersed from bins to voxels, and fine voxels representing coarse voxels (i.e., fine voxels decomposed from coarse voxels) are agglomerated into coarse voxels (step 1342). In this combining step, the particles are dispersed from coarse bins to coarse voxels using coarse parallelepipeds, from fine bins to fine voxels using fine parallelepipeds, from red bins to fine or coarse bins using fine parallelepipeds, and from black bins to coarse voxels using the difference between coarse and fine parallelepipeds. Finally, fluid dynamics are performed on the fine voxels and the coarse voxels (step 1344).

F. Scalar transmission solver

As described above, various types of L BMs can be applied to solve for fluid flow that acts as a background carrier for scalar transfers.

Although a detailed description of the L BM method for simulating fluid flow is provided herein, the following is an example of one method of simulating fluid flow that may be used in conjunction with scalar simulation:

wherein f is_i(x, t) (i ═ 1.., 19.) is the particle distribution function, τ is the individual relaxation time,

is an equilibrium distribution function with a third order spread of fluid velocities,

wherein T is₀1/3. Discrete lattice velocity c_iComprises the following steps:

for stationary particles, w ₀1/3 for states in the cartesian direction, w _i1/18 for a state in the direction of the double diagonal, w_i＝1/36。g_i(x, t) is the external volume force term. The amount of hydrodynamic force and u are moments of the particle distribution function:

as described above, the fluid solver is used in conjunction with a scalar transmission solver that generates scalar transmission information. Thus, in addition to the flow solver, a separate set of distribution functions T is introduced_iFor scalar transmission. Thus, for each voxel in the system, the system simulates both fluid flow and scalar transmission to generate a state vector representing the fluid flow and a scalar quantity representing a scalar variable. These simulation results are stored as entries in a computer-accessible medium.

The set of scalar transfer functions provides an indirect way to solve the second order macroscopic scalar transfer equation. T is_iAn equation is provided that models the dynamic evolution of scalar quantities:

T_iis a scalar distribution function and T is the scalar to be solved for. Tau is_TIs the relaxation time corresponding to the scalar diffusivity. The relaxation time provides a measure of the time required for the system to relax to equilibrium. The term f is defined in equations (38), (39) and (41), respectively_i,ρ,T₀And u.

By c_iThe set of lattice velocities represented may be a set of discrete lattice rates for scalar simulation. In general, lattice velocities for scalar distributionsFor example, when a 19 rate L B model is used for fluid simulation, a 6 rate L B model may be used for scalar simulation, since the 19 rate L B model has a higher order of lattice symmetry than the 6 rate L B model, the same 19 rate model for scalars is used in the examples provided below.

The standard well-known BGK (e.g., as described above) includes non-equilibrium moments of all orders. It is believed that including all of the moments of non-equilibrium is not necessarily isotropic, hydrodynamic, or physically meaningful. Therefore, a BGK regularization/filtering collision operator form is used. Collision operator phi_iRepresenting a future collision factor. The collision operator extracts the unbalanced scalar attribute only at the relevant supported order (e.g., only the first order). The operator also retains and releases the mode of interest while the unbalanced properties associated with non-supported/undesired higher order modes are removed. This projection is sufficient to recover scalar transport physics (e.g., convection and diffusion). Using such a future collision operator is relevant to significantly reduce noise, provide better advection performance (e.g., may be more galileo invariant), and are considered more stable than other solutions of the well-known BGK operator. Such a form (e.g., as shown in equation 43) ensures that only the first order imbalance moment contributes to scalar spread over the hydrodynamic range. By this collision process, all higher order imbalance moments can be filtered out. It is believed that the collision operator Φ is used as described above_iBenefits are provided including elimination of digital noise and improved robustness exhibited in BGK. The scalar T acts as its own balance and does not require a complex expression of scalar balancing distribution functions. Collision operator phi_iThe overall calculation of (a) is quite efficient. It is believed that filtering higher order imbalance moments may additionally provide for reducing imbalance moments at higher ordersBalancing the advantages of aliasing that may exist in the solution.

It can be shown that the collision in (42) obeys the law of scalar conservation. Multiplication by f on both sides of equation (42)_i′(x,t)＝f_i(x+c_iT +1) and note:

the result is

Wherein T is_i' (x, t) denotes the right side of equation (42). Thus, the scalar collision operator preserves the local ρ T, which means that local energy conservation is achieved if the scalar is considered to be temperature. Due to T_iAnd f_iTravel together, so the energy distribution E-f during advection_iT_iIs completely retained. Thus, global conservation of ρ T is achieved. Furthermore, most notably, this scheme preserves the exact invariance of the scalar T uniformity. Thus, it is intuitively obvious that if there is any place

Then at all times thereafter, regardless of the background flow field, phi (x, t) is 0 and

none of the previous lattice boltzmann scalar models exhibit this fundamental property.

Using the Chapmann-Enskog expansion, it can be shown that equation (42) recovers the following second order macroscopic scalar transfer equation:

where κ ═ (τ)_T-1/2)T₀. The condition of constant uniformity ensures that rho is within

And (c) out.

Boundary condition

L BM one significant advantage is the ability to handle complex geometries in the generalized volume L B surface algorithm for achieving frictionless/frictionless Boundary Conditions (BC) on any geometry, mass conservation, and tangential and normal momentum fluxes on the boundary are precisely achieved.

During particle advection, as shown in FIG. 14, each surface element collects incoming particles from its neighboring fluidic cells (step 1410). In contrast to other point-by-point L B, the boundary condition is enforced on a discretized set of surface elements

Weighting is performed (step 1412). After the input quantity is received, the following surface scalar algorithm is applied (step 1414). The outgoing distribution of the surface is not determined,

wherein,

P_i(≡|n·c_i|A)

c_i·nn<0,c_in＝-c_i·P_i＝P_i·n

where n isTo the surface normal of the fluid domain, and c_i·nn<0,c_in＝-c_i。 P_i(≡|n·c_i| A) is the particle direction c associated with the surface normal n and the area A of a given surface element_iVolume of parallelepiped, obviously P_i＝P_i*. The outgoing profile propagates from the surface elements back to the fluid cells according to the same surface advection process (step 1416). It is readily demonstrated that the above described surface scalar collisions achieve an accurate zero surface scalar flux. Summing in the outgoing direction, the outgoing scalar flux is:

note P_i＝P_i*And T in equation (49)ⁱⁿThe second summation term on the right hand side is zero. In addition, due to conservation of mass flux,

the total outgoing scalar flux is the same as the total incoming scalar flux:

therefore, zero net surface flux (adiabatic) BC is fully satisfied on arbitrary geometries. If an external scalar source Q (t) is specified on the surface, the source term can be added directly to the formula (48)

If the boundary condition has a specified scalar quantity T_w(e.g., surface temperature), then the surface heat flux can be calculated accordingly:

Q(t)＝ρC_pκ(T_ω-Tⁱⁿ) Equation (I.52)

Numerical verification

15-18 show four sets of simulation results that demonstrate the capability of the L B scalar solver in terms of its numerical accuracy, stability, Galileo invariance, grid orientation independence, etc. by comparison, the results of using two different second order FD schemes, a flux limiter scheme of type van L eer and a direct hybrid scheme (a hybrid of the central and first order windward schemes) are also presented.

A. Shear wave attenuation

The first test case was a constant, uniform fluid flow carrying a temperature shear wave attenuation, the initial temperature distribution was uniform, plus a lattice wavelength L of 16, amplitude

Is sinusoidally varied. T is_AThe velocity of the background mean flow is 0.2, the thermal diffusivity κ is 0.002, so low resolution and κ can prove numerical stability and accuracy well, L B scalar solver and finite difference methods show excellent agreement with theory for temperature decay without background flow in the case of non-zero background mean flow L B scalar solver can still compare exactly with theory, however, FD results show significant numerical error, temperature curves at lattice time step 81 are plotted in fig. 15.

B. Inclined channel with heat source

A second test case was to simulate the temperature distribution in channel flows with different lattice orientations. The channel walls are free to slide, so the fluid flow remains uniform, as a result of which U₀0.2. The channel width was 50 (lattice spacing) and the flow rate Re was 2000. The thermal diffusivity, κ, was 0.005. Temperature on the wall at T _W1/3 fixed, constant volume heat source q 5 × 10 is used in the bulk fluid domain^-6. The flow has periodic boundary conditions in the flow direction, which in the case of lattice alignmentThis is easily achieved when the channel (light) tilt is tilted as shown in FIG. 16, the input and output channel boundaries are perfectly matched in the coordinate direction, thus again achieving periodic boundary conditions in the flow direction, to demonstrate lattice independence, we chose a tilt angle of 26.56 degrees, like the first test case, when the channel is lattice aligned, the temperature distribution using the L B scalar solver and the two FD schemes is very well matched with the analytical solution.

C. Temperature propagation in inclined channels

In addition, the computation of scalar advection may be further compromised for FD methods due to strong dependence on windward information.instead, as described above, the boundary processing of the L B scalar solver achieves precise local scalar flux conservation.scalar advection in such near-wall regions may be precisely computed.As an example, high temperature convection is performed in a channel tilted 30 degrees₀0.0909. thermal diffusivity, k, is 0.002. initially, the temperature was 1/3 everywhere except at the inlet, T, 4/9. then, the temperature front should convect to downstream locations at a later time without distortion due to uniform background fluid flow. the calculated temperature front distribution across the channel at lattice time step 2000 is shown in fig. 17. the temperature front of the L B scale solver remains an approximately straight curveIt is also worth mentioning that the L B scalar results show the thinnest temperature front, which means that the L B scalar solver has less numerical spread.

D. Rayleigh-Bernard natural convection

The Rayleigh Bernard (RB) natural convection is the classical reference for numerical solver accuracy verification it has a simple case set-up, but has complex physical phenomena when the Rayleigh number Ra exceeds a certain threshold, the system undergoes a transition from no flow to convection, as shown in FIG. 18, the current study is conducted under the Businessq approximation where the buoyancy acting on the fluid is defined as α ρ g (T-T)_m) Wherein α is the thermal expansion coefficient, g is the gravity, T_mThe average temperature values for the top and bottom boundaries. Since when R is_aExceeds a critical value R_a1707.762, the most unstable wavenumber is κ_c3.117, the resolution used in the study was therefore 101 × 50. Pr used here is 0.71. after RB convection was established, the heat transfer between the two plates was greatly enhanced, Nu 1+<u_yT>H/κ Δ T describes the enhancement of heat transfer.

Fluid simulation using the above total energy conservation (rather than the conventional approach) would provide a similar description of the fluid simulation and any customized corresponding computational data output. However, when an object (e.g., an actual physical object) is modeled, such fluid simulation using the curve method described above may be performed faster than other methods and/or by using less computational resources than other methods.

The equilibrium distribution discussed herein does not involve a pressure term and therefore its value is always positive during crash operations, thereby providing relatively greater stability. After the collision state convects as it is, an isothermal pressure gradient in momentum is caused, and no pressure convection term exists in the total energy equation. Conservation of total pDuty is accomplished by modifying the stall state, such as adding the local pressure term back to the particle state after advection, before calculating the moment that provides the appropriate pressure gradient.

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

Claims (20)

1. A method for simulating fluid flow on a computer, wherein the simulation achieves total energy conservation, the method comprising:

simulating the movement of the fluid on the grid, the movement of the fluid being simulated to model the movement of the particles on the grid;

storing, in a computer-accessible memory, a set of state vectors for each grid location in the grid, each state vector comprising a plurality of entries corresponding to a particular momentum state of the possible momentum states at the corresponding grid location;

calculating, by a computer, a set of energy values for grid locations in a grid;

performing, by the computer, advection of the particles to subsequent grid locations for the time interval; and

modifying, by the computer, the set of state vectors for the particles by adding the particular total energy value to the states of the advected particles and subtracting the particular total energy value from the states of the particles that were not advected for the time interval.

2. The method of claim 1, wherein after advecting particles to subsequent grid locations according to the modified state set, the method further comprises:

adding local pressure terms back to the particle state prior to computing the moment to provide an appropriate pressure gradient

3. The method of claim 2, wherein the local pressure terms comprise θ -1 terms calculated by the computer.

4. The method of claim 2, wherein the activity of simulating fluid flow comprises simulating fluid flow based in part on a first set of discrete lattice rates, and the method further comprises:

the time evolution of the scalar is modeled based in part on a second set of discrete lattice rates, the second set of discrete lattice rates being the same as the first set of discrete lattice rates, or the second set of discrete lattice rates having lattice rates differing in number from the first set of discrete lattice rates.

5. The method of claim 1, wherein simulating the time evolution of total energy comprises:

collecting incoming distributions from adjacent grid locations for collision and energy operators;

weighting the incoming distribution;

determining an outgoing distribution as a product of collision and energy operators; and

propagating the determined outgoing distribution.

6. The method of claim 1, wherein to recover accurate shear stress for an arbitrary Prandtl (Prandtl) number, the energy collision term is given by:

wherein, pi is ∑_ic_ic_i(f_i-f_i ^eqA filtered second moment that is an unbalanced component.

7. The method of claim 5, further comprising:

a zero net surface flux boundary condition is applied such that the incoming distribution is equal to the determined outgoing distribution.

8. The method of claim 1, wherein computing a set of energy values for grid locations in a grid further comprises computing a mean by E_t＝E+v²Total energy E given by/2_tWhere E is energy and v is velocity.

9. The method of claim 1, wherein simulating an activity comprises simulating an activity of the patient

By computer for a particular total energy E_tApplying a balanced distribution and a second distribution function, wherein the second distribution is defined as a distribution f with the flow_iA specific scalar of advection.

10. The method of claim 9, wherein the second distribution function takes into account f_iAn unbalanced contribution to the energy equation to obtain the correct flow behavior near the boundary and across different grid resolutions, wherein the collision operator for the distribution function is given by:

wherein the term omega_f,Ω_qRepresenting the respective collision operator; the equilibrium distribution used in the above equation is:

11. a computer program product tangibly embodied in a computer-readable medium, the computer program product including instructions that, when executed, simulate a physical process fluid flow and cause a computing system to:

simulating the movement of the fluid on the grid, the movement of the fluid being simulated to model the movement of the particles on the grid;

storing, in computer memory, a set of state vectors for each grid location in the grid, each state vector comprising a plurality of entries corresponding to a particular momentum state of the possible momentum states at the corresponding grid location;

computing a set of energy values for grid locations in the grid;

performing advection of the particles to subsequent grid locations for the time interval; and

the set of state vectors for the particles is modified by adding a specific total energy value to the states of the advected particles and subtracting the specific total energy value from the states of the particles not advected in the time interval.

12. The computer program product of claim 11, wherein total energy adds the term before advection such that pressure is convected and energy is removed from rest to compensate for the added energy.

13. The computer program product of claim 12, wherein removing the added energy conserves total energy and provides the correct pressure velocity term.

14. The computer program product of claim 13, wherein the stop state is modified by

So that

Providing for the removal of the added energy.

15. A computer system for simulating physical process fluid flow, the system comprising:

one or more processor devices;

computer memory coupled to the one or more processor devices; and

a computer-readable medium storing instructions that, when executed, simulate a physical process fluid flow and cause a computing system to:

simulating the movement of the fluid on the grid, the movement of the fluid being simulated to model the movement of the particles on the grid;

storing, in computer memory, a set of state vectors for each grid location in the grid, each state vector comprising a plurality of entries corresponding to a particular momentum state of the possible momentum states at the corresponding grid location;

computing a set of energy values for grid locations in the grid;

performing advection of the particles to subsequent grid locations for the time interval; and

the state vector of the particle is corrected by adding a specific total energy value to the state of the advected particle and subtracting the specific total energy value from the state of the particle not advected in the time interval.

16. The system of claim 15, wherein the system further comprises instructions to:

after the particles advect to subsequent grid locations according to the modified state set, local pressure terms are added back to the particle states before moments are calculated to provide the appropriate pressure gradients

17. The system of claim 16, wherein the local pressure terms comprise θ -1 terms calculated by the system.

18. The system of claim 16, wherein the instructions to simulate the activity of the fluid flow comprise instructions to:

simulating fluid flow based in part on the first set of discrete lattice rates, an

The time evolution of the scalar is modeled based in part on a second set of discrete lattice rates, the second set of discrete lattice rates being the same as the first set of discrete lattice rates, or the second set of discrete lattice rates having lattice rates differing in number from the first set of discrete lattice rates.

19. The system of claim 15, wherein the instructions to compute the set of energy values for the grid locations in the grid comprise computing a mean from E_t＝E+v²Total energy E given by/2_tWhere E is energy and v is velocity.

20. The system of claim 15, further comprising instructions that cause the system to:

for a specific total energy E_tApplying a balanced distribution and a second distribution function, wherein the second distribution is defined as a distribution f with the flow_iA specific scalar of advection.

Images