JP2020201961A - Computer simulation of physical fluid on irregular spatial grids for stabilizing explicit numerical diffusion problem - Google Patents

Abstract

To provide a computer simulation of physical fluid on irregular spatial grids for stabilizing explicit numerical diffusion problem.SOLUTION: Computer implemented techniques for simulating a fluid flow around a surface of a solid includes: receiving a model of a simulation space including a lattice structure represented as a collection of voxels and a representation of a physical object; simulating movement of particles in a volume of fluid; identifying faces between two voxels by a computing system; computing a modified flux by using a spatially averaged gradient in the vicinity of the two voxels; and performing by the computing system, advection operation on the particles to subsequent voxels.SELECTED DRAWING: Figure 1

Description

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

Computational fluid dynamics is a division of fluid dynamics that includes computerized numerical analysis techniques for analyzing and simulating fluid flow in a physical environment.

The so-called "Lattice Boltzmann method" (LBM) is an advantageous technique used in computational fluid dynamics. The dynamics underlying the Lattice Boltzmann system belong to the basic physics of kinetic theory, which includes the motion of many particles according to the Boltzmann equation. The basic Boltzmann motor system has two basic dynamic processes: collision and advection. The collision process involves interactions between particles that relax to equilibrium according to conservation law. The advection process involves modeling the movement of a particle from one location to another according to the microscopic velocity of the particle.

One common aspect of actual computational fluid dynamics simulation problems, including the "Lattice Boltzmann method" (LBM), is that of complex shapes, such as those involving irregular grids due to surface and volume discretization. ..

Numerical simulations of diffusion-dominant physical phenomena are commonly used in applications including conductive heat transfer, mass diffusion, electrical conduction, and the like. The governing equations for these phenomena are formulated as a set of PDEs that include the unsteady (unstable) diffusion term and the volume source term. Numerical solutions include discretizing the spatial domain of interest and applying time integration techniques to advance the solution in time. Spatial discretization procedures are usually accomplished using highly automated grid generation tools, while time discretization (time integration techniques are applied to advance the solution in time, i.e., in time step size). For) is chosen to guarantee the stability and accuracy of the numerical solution at an acceptable computational cost.

In particular, the stability property, commonly referred to as the Courant-Friedrichs-Lewy (CFL) constraint of the time forward scheme, determines the maximum time step size that can be used without introducing significant instability into the solution. Two types of time progression schemes, namely implicit and explicit schemes, are commonly used.

Implicit methods meet the CFL constraints by construction, so large time steps can be used without destabilizing the solution (however, too large time steps usually give inaccurate results). Implicit methods require solutions to large systems of matrix coefficients, thus making their actual implementation non-trivial and computationally expensive.

On the other hand, the explicit method is very simple to implement, computationally inexpensive (per iteration), and highly parallelizable. However, the explicit method needs to meet the CFL constraint. This CFL constraint of the explicit diffusion scheme requires that the number of CFLs given by (κΔ_t) / (Δ_x ^ 2) be less than a certain limit (O (1)), where κ is. Diffusivity, Δ_x is the size of the smallest space grid, and Δ_t is the time step. In other words, if the spatial grid size Δ_x is reduced by a factor F somewhere in the domain, then the time step Δ_t must be reduced by a factor F ² to maintain numerical stability.

Therefore, the explicit method may require extremely small time steps in a spatial grid with small size elements that significantly affect simulation performance. This is true even if the number of such small size elements is very limited within the simulation domain, because the smallest element within the entire domain determines the CFL condition, resulting in a time step. This is because it determines the size.

In real problems involving complex shapes, the use of irregular grids is unavoidable in surface and volume discretization. On these grids, Δ_x can vary significantly and the use of explicit methods can be very inefficient due to the extremely small time steps required by the CFL constraints.

Therefore, explicit practitioners spend a great deal of time and effort trying to improve the quality of the spatial grid and try to mitigate the problem. Even then, it is almost impossible to remove all small-sized elements from the discretization of realistic spatial domains, so that small time steps (at least locally) stabilize the explicit solution. The only way to get it done.

According to one aspect, systems, computer implementation methods, and computer program products for simulating fluid flow around the surface of a solid have a lattice structure represented as a collection of voxels and a representation of physical objects. Receiving a model of the simulation space that contains, the voxel has the appropriate resolution to describe the surface of a physical object, receiving and simulating the movement of particles in the volume of a fluid. Simulating that the movement of particles causes collisions between particles and identifying the faces (facets) between two voxels, identifying that at least one of the faces violates stability conditions. And to calculate the modified heat flux using the spatially averaged temperature gradient in the vicinity of the two voxels, and to calculate that at least one of the surfaces violates the stability condition. Includes performing a particle transfer operation to a subsequent voxel.

The disclosed technique introduces a modification to the heat flux calculation between two adjacent elements if at least one of the elements violates a constraint, eg, a Courant-Friedrichs-Lewy (CFL) constraint). .. These modifications to the heat flux calculation depend on the material and geometric properties of the two elements, as well as the existing state of the quantity of objects in the immediate vicinity of the elements, and stabilize the numerical solution regardless of the size of the two elements. Helps to ensure space-time accuracy. If the two adjacent elements are large (and therefore satisfy the CFL constraint), the newly proposed flux calculation becomes an explicit method embodiment, which is that this new method is consistent with the explicit method. It means that the above deficiencies are overcome by an explicit method.

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

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

It is a figure which shows the system for the simulation of a fluid flow. It is a flow chart which shows the operation simulation based on the lattice Boltzmann model using the heat flux calculation. FIG. 5 is a flow diagram comprising a mode of modification to heat flux calculation between two adjacent elements when at least one of the elements violates the constraint. FIG. 5 is a flow diagram comprising a mode of modification to heat flux calculation between two adjacent elements when at least one of the elements violates the constraint. FIG. 5 is a flow diagram comprising a mode of modification to heat flux calculation between two adjacent elements when at least one of the elements violates the constraint. FIG. 5 is a flow diagram comprising a mode of modification to heat flux calculation between two adjacent elements when at least one of the elements violates the constraint. FIG. 5 is a flow diagram comprising a mode of modification to heat flux calculation between two adjacent elements when at least one of the elements violates the constraint. It is a figure which shows the velocity component of two LBM models represented in Euclidean space (prior art). It is a figure which shows the velocity component of two LBM models represented in Euclidean space (prior art). It is a flow chart of the procedure followed by a physical process simulation system using a modified heat flux calculation between two adjacent elements. It is a perspective view of a microblock (prior art). It is a figure of the lattice structure used by the system of FIG. 1 (prior art). It is a figure of the lattice structure used by the system of FIG. 1 (prior art). It is a figure which shows the variable resolution technique (prior art). It is a figure which shows the variable resolution technique (prior art). It is a figure which shows the movement of a particle (prior art). It is a figure which shows the region affected by the facet of a surface (prior art). It is a figure which shows the surface dynamics (prior art). It is a flow chart of the procedure for carrying out surface dynamics.

Model simulation space In an LBM-based physical process simulation system, the fluid flow is represented by a distribution function value f _i evaluated by a set of discrete velocities c _i . The dynamics of the distribution function is described in Equation I. Dominated by 1.
f _i (x + c _i , t + 1) = f _i (x, t) + C _i (x, t) Equation (I.1)
This equation is a well-known lattice Boltzmann equation describing the time evolution of the distribution function f _i. The left side represents the change in distribution due to the so-called "flow process". The flow process is when a fluid pocket starts at one mesh position and then moves along one of a plurality of velocity vectors to the next mesh position. At that point, the "collision coefficient", i.e., the effect of the fluid's neighborhood pocket on the fluid's starting pocket, is calculated. The fluid can only move to different mesh positions, so proper selection of velocity vectors is required so that all components of all velocities are multiples of the common velocity.

The right side of the first equation is the "collision operator" described above, which represents a change in the distribution function due to a collision between fluid pockets. The particular form of the collision operator can be, but is not limited to, the form of the Bhatnagar, Gross, and Krook (BGK) operators. The collision operator brings the distribution function closer to the default value given by the second equation, which is the "equilibrium" form.

The BGK operator has a distribution function that collides, regardless of the details of the collision.

を介して、｛ｆ^eq（ｘ，ν，ｔ）｝で与えられる明確な局所平衡に近づくという物理的論拠によって構築され、ここで、パラメータτは、衝突を介した平衡までの固有緩和時間を表す。

Is constructed by the physical rationale of approaching the clear local equilibrium given by {f ^eq (x, ν, t)}, where the parameter τ sets the intrinsic relaxation time to equilibrium via collision. Represent.

From this simulation, conventional fluid variables such as mass ρ and flow velocity u are obtained as a simple sum in equation (I.3). See below.

Due to symmetry considerations, the set of velocity values is chosen to form a particular lattice structure across the arrangement space. Such discrete dynamics of the form _{f i (x + c i,} t + 1) = f i (x, t) + C i (x, t)
According to the LBM equation with, where the collision operator usually takes the BGK form as described above. With proper selection of equilibrium distribution forms, it can be theoretically shown that the lattice Boltzmann equation yields correct hydrodynamic and thermo-hydrodynamic results. That is, the hydrodynamic moment derived from f _i (x, t) is Navier a macroscopic extreme - according to Stokes equations. These moments
ρ (x, t) = Σ _i f _i (x, t); ρ (x, t) u (x, t) = Σ _i c _i f _i (x, t) Equation (I.3)
Defined as, where ρ and u are fluid densities and velocities, respectively.

The set values of c _i and ω _i define the LBM model. The LBM model can be efficiently implemented on an extensible computer platform and can be executed with excellent robustness against temporally unsteady flows and complex boundary conditions.

The standard technique for obtaining the macroscopic equation of motion of a fluid system from the Boltzmann equation is the Chapman-Enskog method, which incorporates a successive approximation of the complete Boltzmann equation. In a fluid system, low-density disturbances move at the speed of sound. In gas systems, the speed of sound is generally determined by temperature. The importance of the compressible effect on the flow is measured by the ratio of characteristic velocity to sound velocity, known as the Mach number.

For further description of conventional LBM-based physical process simulation systems, an interpretation is referred to in US Patent Application Publication No. 2016-0188768, the entire contents of which are incorporated herein by reference.

With reference to FIG. 1, a system 10 for simulating a fluid flow is shown, for example with respect to the representation of a physical object. The system 10 of this embodiment includes a server system 12 implemented as a large-scale parallel computing system 12 based on a client-server architecture, and a client system 14. The server system 12 includes a memory 18, a bus system 22, an interface 20 (eg, user interface / network interface / display or monitor interface, etc.), and a processing device 24 that performs a simulation process 30 together with a mesh and a simulation engine 34. including.

The memory 18 includes a mesh preparation engine 32 and a simulation engine 34. Although FIG. 1 shows the mesh preparation engine 32 in memory 18, the mesh preparation engine can be a third party application running on a system different from the server 12 (eg, system 14 or another system). .. Whether the mesh preparation engine 32 runs in memory 18 or on a different system than the server 12, the mesh preparation engine 32 receives a user-specified mesh definition 28, and the mesh preparation engine 32 receives the mesh. Prepare and send the prepared mesh to the simulation engine 34.

The simulation engine 34 includes a particle collision interaction module 34a, a particle boundary model module 34b, and an advection module 34c that performs an advection operation. System 10 accesses a data repository 38 that stores 2D and / or 3D meshes and libraries. The advection module 34c includes a submodule 36 that performs an advection operation according to a modification of the flux calculation between two adjacent elements, as discussed below.

Before running the simulation on the simulation engine, the simulation space is modeled as a set of voxels. Generally, the simulation space is generated using a computer-aided design (CAD) program. For example, a CAD program can be used to draw a microdevice located in a wind tunnel. The data generated by the CAD program is then processed to add a lattice structure with appropriate resolution to reveal objects and surfaces in the simulation space.

Next, with reference to FIG. 2, a process of simulating a fluid flow with respect to the representation of a physical object is shown. In the examples discussed herein, the physical object is a wing. However, the use of wings is merely exemplary, as physical objects can be of any shape and can have planes and / or curved surfaces in particular. The process receives a mesh for a physical object that is simulated, for example, from a client system 14 or by a search from a data repository 35a. In another embodiment, either the external system or the server 12 generates a mesh of simulated physical objects based on user input. The process precalculates the geometry from the retrieved mesh 35b and performs a dynamic lattice Boltzmann model simulation using the precomputed geometry corresponding to the retrieved mesh 35c. The Lattice Boltzmann model simulation is the simulation of the evolution 35d of the particle distribution 34a and the next cell q (bar above) in the LBM mesh by the boundary determination (not shown in FIG. 2) created from the engine 34b (FIG. 1). ) Includes the advection of particles 35e. The advection 35c process tests the CFL constraint violation 37a and, if found, applies the correction to the flux calculation 37b using the engine 36 (FIG. 1).

Stabilization of the explicit numerical scheme of the diffusion problem for irregular spatial grids For this explanation, an explicit Euler scheme and a finite volume formulation are assumed. In the following description, the quantity of interest is temperature and the governing equation is the heat conduction equation. Numerical schemes need to calculate the heat flux on all sides of the element. These fluxes are then summed and used to update the temperature of the element under consideration. Consider two face-sharing adjacent elements α and β. According to Fourier's law of heat conduction, the heat flux is

であり、ここで、

And here,

は、共通面の熱伝導率であり、

Is the thermal conductivity of the common surface,

は、共通面に垂直な温度勾配であり、「ｍ」は、量が時間ステップ「ｍ」で評価されることを明示するために使用される。熱流入α（熱流出αの代わりに）が考慮されるので、フーリエの法則の一般に使用される形式の負号は省かれる。使用される温度勾配は、特に、異なるサイズの要素が存在する状態で、粒子移流の滑らかさを確実するように計算される。２つの隣接する要素αおよびβがＣＦＬの制約を満たす場合、時間ステップｍからｍ＋１への間に共通面を横切るエネルギー移送は、熱流束に共通面の面積Ａ^αβおよび時間ステップサイズΔ_tを乗算することによって得られ、すなわち、

Is the temperature gradient perpendicular to the common plane and "m" is used to indicate that the quantity is evaluated in the time step "m". Since heat inflow α (instead of heat outflow α) is considered, the negative sign of the commonly used form of Fourier's law is omitted. The temperature gradients used are calculated to ensure the smoothness of particle advection, especially in the presence of elements of different sizes. If two adjacent elements α and β satisfy the CFL constraint, the energy transfer across the common plane between the time step m and m + 1 multiplies the heat flux by the area A ^{α β of the} common plane and the time step size Δ _t . Obtained by doing, i.e.

である。

Is.

In the conventional method, the final temperature of the element α at the end of the time step is calculated from the net energy transfer to α (the sum of the energy transfers from all planes).

The above equation (3), the temperature change is proportional to the net heat flux is inversely proportional to the size ∀ ^alpha elements, i.e., the small element, that the same net energy transfer is stated to bring a larger temperature change Please note.

With reference to FIG. 3, the flux correction calculation process 36 includes determining the time step. 36a where the time step is tested in the flux calculation process 36. If the time step Δ _t is not as large as 36b, which violates the CFL constraint on at least one of the two factors, then the above form (Equation 3) is less likely to result in numerical instability in the solution. , The above-mentioned form of flux calculation can be used. If the time step Δ _t is as large as 36c, which violates the CFL constraint on at least one of the two elements, the above-mentioned form (Equation 3) may result in numerical instability in the solution. In this case, the modified flux calculation method 36d is applied.

Without loss of generality, the modification method can be explained as follows.

It is assumed that element α is smaller than adjacent element β and therefore the CFL constraint is violated at least at element α. This numerical instability is due to the element α (assumed to be small in size).

を計算するために使用される温度勾配が、時間ステップΔ_tの期間の全体を通して一定値のままであると仮定することが正しくないという理由で生じる。上記のように、同じ正味エネルギー移送では、小さい要素の温度変化は大きく、したがって、標準の陽的時間積分は、一定温度勾配の仮定が有効であることを保証するために時間ステップを減少させることを必要とする。所与のΔ_tでは、この問題は、問題に非定常性が存在する限り存続し、すべての要素に関して入って来るおよび出て行く流束がすべて互いに正確に釣り合うときの定常状態でのみ解決する。

It arises because it is incorrect to assume that the temperature gradient used to calculate the value remains constant throughout the period of time step _Δt . As mentioned above, for the same net energy transfer, the temperature change of the small elements is large, so the standard explicit time integral reduces the time step to ensure that the constant temperature gradient assumption is valid. Needs. At a given Δ _t , this problem persists as long as there is non-stationarity in the problem and is solved only in the steady state when all the incoming and outgoing fluxes for all elements are precisely balanced with each other. ..

Referring to FIG. 4, as part of the correction method, the flux-corrected flux calculation is a term defined as in equation (1).

を、以下の２つの部分に細分する４９ａ。２つの部分は、（１ａ）αの温度発展（上述の合計において）に向けて使用されることになる付加流束（ａｐｐｌｉｅｄ ｆｌｕｘ）４９ｂ

Is subdivided into the following two parts 49a. The two parts are (1a) an applied flux 49b that will be used for the temperature evolution of α (in the sum described above).

と、（１ｂ）要素αの温度変化なしに、インタフェースαβの反対側に送出されることになる平衡流束４９ｃ

And (1b) the equilibrium flux 49c that will be sent to the opposite side of the interface αβ without the temperature change of the element α.

とであり、それらは、以下のように、

And they are, as follows:

として表され、ここで、項

Expressed as, where the term

およびΔＧは、

And ΔG are

によって与えられ、そして、

Given by, and

として表される。

It is expressed as.

In the above formula, the geometric features, the distance d ^.alpha..beta used in the calculation of the temperature gradient, and element volume ∀ ^alpha and ∀ ^beta, by the area A ^.alpha..beta the common plane (shared by adjacent elements) expressed. Material properties are

および

and

（ここで、ρは密度を示し、Ｃ_pは比熱を示す）、ならびに熱伝導率

(Here, ρ indicates density and C _p indicates specific heat), as well as thermal conductivity.

によって説明される。加算における流束項

Explained by. Flux term in addition

および

and

は、それぞれ、要素αおよびβの異なる面に存在する可能性がある流束の推定量を提供する。

Provides an estimator of flux that may be present on different faces of elements α and β, respectively.

Two flux terms

および

and

の物理的解釈は、以下の通りである。付加流束

The physical interpretation of is as follows. Additional flux

は、数値不安定を導入することなく要素αの温度発展に向けて使用することができる全流束

Can be used for the temperature evolution of element α without introducing numerical instability

（式（１）により与えられる）の一部分を表す。この形態は、要素αおよびβからなる孤立系のための第１原理から導出され、この系の近くで連続的に発展する温度場の影響の推定量を含むことができる。この理由で、この項は、要素の幾何学的／熱的特性

Represents a part of (given by equation (1)). This form is derived from first principles for isolated systems consisting of elements α and β and can include an estimator of the effects of temperature fields that develop continuously near this system. For this reason, this term refers to the geometric / thermal properties of the element.

ならびに環境とのこれらの要素の相互作用（ΔＧ）の両方への依存を示す。

It also shows its dependence on both the interaction of these factors with the environment (ΔG).

In equation (4)

は、依然として、以前の時間ステップに他の面で観察された流束に依存しており、それによって、現在の時間ステップの間にそれらの面で継続している相互作用の推定量を提供することに留意すべきである。強硬な過渡的問題では、これは、

Still relies on the flux observed in other planes in the previous time step, thereby providing an estimator of the interactions that continue in those planes during the current time step. It should be noted that. In a tough transitional problem, this is

と

When

との間にずれをもたらす。

Brings a gap between and.

A second term called the equilibrium flux

は、このずれを説明しており、それはαを通ってインタフェースαβの反対側の隣接する要素に伝えられる。この平衡流束は、ＣＦＬ制約を満たすのに十分な大きさの要素に置かれるときのみ温度発展に付加され、それまで、流束方向に沿って連続的に送出される。

Explains this deviation, which is transmitted through α to the adjacent element on the opposite side of interface αβ. This equilibrium flux is added to the temperature evolution only when placed on an element large enough to satisfy the CFL constraint, and until then it is continuously delivered along the flux direction.

The scheme described above is the correct amount of total flux at the interface between the two elements α and β, while precisely controlling the amount of flux effective for temperature evolution of the small element α.

が組み入れられることを厳格に保証する。全体として、このスキームは、数値安定、ならびに良好な空間および時間精度を維持することができる。最後に、定常状態では、付加流束は全流束に等しくなり、すなわち、

Is strictly guaranteed to be incorporated. Overall, this scheme can maintain numerical stability as well as good spatial and temporal accuracy. Finally, in steady state, the additional flux is equal to the total flux, ie

であり、結果として、平衡流束は０に等しくなり、すなわち、

As a result, the equilibrium flux is equal to 0, i.e.

であることに留意すべきである。

It should be noted that

Characteristics of the Modified Flux Calculation Process As mentioned in the above description and as shown in FIGS. 5-7, the modified flux calculation method 36 requires several algorithmic processes. The modified flux calculation process involves identifying the plane to which the modified definition of heat flux should be applied 52a. 52b (Equation 3) uses the standard definition of heat flux for all aspects that do not violate Courant conditions. A modified heat flux calculation is used in all aspects that violate the CFL condition, eg, in any aspect between two elements that at least one violates the CFL condition criteria 52c. The total heat flux amount is calculated using a spatially averaged temperature gradient in the vicinity of the element under consideration to ensure the smoothness of the solution. (In contrast, the standard heat flux utilizes a temperature gradient calculated based on the conventional differential form.) The modified heat flux calculation process 52 applies the modified heat flux to two parts, namely the additional flux 49b. Includes 52d splitting into terms and equilibrium flux 49c terms. The flux calculation is then performed according to the division 52e.

Next, referring to FIG. 6, the flux calculation uses the additional flux calculation 52f (Equation 3), which is always used in the temperature evolution equation of the element under consideration.

Next, referring to FIG. 7, the equilibrium flux may or may not be used for temperature development, depending on the size of the elements involved. If the element is 52g small enough to violate the CFL constraint, the equilibrium flux is calculated 52h, sent to the element adjacent in the direction of the flux 52i, and if the adjacent element is large enough not to violate the CFL constraint, the equilibrium flux However, 52j is used for temperature development. The equilibrium flux is in the direction of the flux until the equilibrium flux is finally transferred to a sufficiently large element (not to violate the CFL constraint) that applies the equilibrium flux towards temperature evolution. It is sent continuously along.

Uniqueness and Advantages of Algorithms Several methods are known to overcome the problem of numerical instability in the diffusion problem of irregular grids with varying element sizes. The most commonly used technique is to impose additional constraints on the grid generation tool to reduce such scenarios. Even then, this problem cannot be completely avoided, so it is common to use global time steps small enough to guarantee stability, or to use local subcycling when encountering small grid elements. Method. The first method (small global time step) effectively increases the computational cost even if a single small size element occurs somewhere in the spatial grid, while the second method (local subcycling). Increases the complexity of the algorithm and its implementation. An alternative approach is to use the implicit method instead of the explicit method, or at least limit the use of the explicit method to local areas near small elements. This implicit method method suffers from the complexity of implementation and the nonlocal nature of the solution that results in simultaneous equations that are inconvenient for massively parallel processing computer implementation.

In contrast, the modified flux calculation method offers some distinct advantages. The modified flux calculation method allows the use of a single time step size selected based on time accuracy considerations rather than the size of the smallest element of the grid. This is a great advantage in terms of computational cost and ease of implementation for all possible scenarios. The modified flux calculation method has a dependency on the geometric properties of the two adjacent elements, thereby ensuring that the modified flux calculation method works regardless of the size of the element. Therefore, the usual constraints on the grid generation process (grid quality, size, etc.) can be significantly relaxed. The computational cost of computing various terms is reasonable because the mathematical form of the terms is simple and does not include iterations. This is in sharp contrast to existing methods that can have high computational costs (reducing the time step size or using the hybrid implicit-explicit method). In addition, due to the volumetric properties of the formulation, this scheme maintains accurate storage, which is a requirement in many applications. Modified flux calculations are still explicit in nature and require information from elements within a small distance of the element under consideration. Yang means that minimal modifications to existing computing systems are required according to the original system embodiments. The parallelization characteristics of the original explicit method are preserved, and therefore the modified flux calculation method can utilize massively parallel processing computer embodiments.

Therefore, the modified flux calculation method has advantages over existing methods. Although the above description of the flux calculation correction method was based on the finite volume formulation of heat conduction, this method is practically applicable to many diffusion dominating problems.

Model simulation space In an LBM-based physical process simulation system, fluid flow is represented by a distribution function value f _i evaluated with a set of discrete velocities c _i . The dynamics of the distribution function is governed by equation I1, where _fi (0) is known as the equilibrium distribution function.

として定義され、ここで、

Defined as, here,

である。

Is.

This equation is a well-known lattice Boltzmann equation describing the time evolution of the distribution function f _i. The left side represents the change in distribution due to the so-called "flow process". In the flow process, the fluid pocket starts at one mesh position and then moves along one of the velocity vectors to the next mesh position. At that point, the "collision coefficient", i.e., the effect of the fluid's neighborhood pocket on the fluid's starting pocket, is calculated. The fluid can only move to different mesh positions, so proper selection of velocity vectors is required so that all components of all velocities are multiples of the common velocity.

The right side of the first equation is the "collision operator" described above, which represents a change in the distribution function due to a collision between fluid pockets. A particular form of collision operator is the form of the Bhatnagar, Gross, and Krook (BGK) operators. The collision operator brings the distribution function closer to the default value given by the second equation, which is the "equilibrium" form.

The BGK operator has a distribution function that collides, regardless of the details of the collision.

を介して、｛ｆ^eq（ｘ，ν，ｔ）｝によって与えられる明確な局所平衡に近づくという物理的論拠によって構築され、ここで、パラメータτは、衝突を介した平衡までの固有緩和時間を表す。粒子（例えば、原子または分子）を扱う場合、緩和時間は、一般に、定数と見なされる。

Constructed by the physical rationale of approaching the clear local equilibrium given by {f ^eq (x, ν, t)}, where the parameter τ sets the inherent relaxation time to equilibrium via collision. Represent. When dealing with particles (eg atoms or molecules), relaxation time is generally considered a constant.

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

ここで、ρ、ｕ、およびＴは、それぞれ、流体の密度、速度、および温度であり、Ｄは、離散化された速度空間次元である（物理的空間次元と必ずしも等しくない）。

Where ρ, u, and T are the density, velocity, and temperature of the fluid, respectively, and D is the discrete velocity space dimension (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 across the arrangement space. Such discrete dynamics of the form _{f i (x + c i,} t + 1) -f i (x, t) = C i (x, t)
According to the LBM equation with, where the collision operator usually takes the BGK form as described above. With proper selection of equilibrium distribution forms, it can be theoretically shown that the lattice Boltzmann equation results in correct fluid and thermo-fluid dynamics. That is, the hydrodynamic moment derived from f _i (x, t) is Navier a macroscopic extreme - according to Stokes equations. These moments are defined by the above equation (I3).

The set values of c _i and w _i define the LBM model. The LBM model can be efficiently implemented on an extensible computer platform and can be executed with excellent robustness against temporally unsteady flows and complex boundary conditions.

The standard technique for obtaining the macroscopic equation of motion of a fluid system from the Boltzmann equation is the Chapman-Enskog method, which incorporates a successive approximation of the complete Boltzmann equation. In a fluid system, low-density disturbances move at the speed of sound. In gas systems, the speed of sound is generally determined by temperature. The importance of the compressible effect on the flow is measured by the ratio of characteristic velocity to sound velocity, known as the Mach number.

A general discussion of LBM-based simulation systems that can be used with CAD processes to simulate fluid flow is provided below.

Referring to FIG. 8, the first model (2D-1) 200 is a two-dimensional model containing 21 velocities. Of these 21 velocities, one (205) represents a non-moving particle and three sets of four velocities are either positive or negative along either the x-axis or the y-axis of the lattice. Either the normalization speed (r) (210-213), twice the normalization speed (2r) (220-223), or three times the normalization speed (3r) (230-233). Representing a moving particle, two sets of four velocities are normalized velocities (r) (240-243) or twice the normalized velocity (2r) (2r) for both the x- and y-velocities. Represents a moving particle in 250-253).

With reference to FIG. 9, a second model (3D-1) 260, a three-dimensional model containing 39 velocities, is shown, each velocity being represented by one of the arrows in FIG. To. Of these 39 velocities, one represents a non-moving particle, and three sets of six velocities are either positive or negative along the x-axis, y-axis, or z-axis of the lattice. Crab represents a particle moving at any of the normalized velocity (r), twice the normalized velocity (2r), or three times the normalized velocity (3r), eight of which are x, y, z Represents particles moving at a normalization velocity (r) for all three of the lattice axes, with 12 being twice the normalization velocity for two of the x, y, and z lattice axes ( Represents a moving particle in 2r).

More complex models such as the 3D-2 model with 101 velocities and the 2D-2 model with 37 velocities can also be used. In the three-dimensional model 3D-2, one of the 101 velocities represents a non-moving particle (Group 1), and three sets of six velocities are the x-axis, y-axis, or z-axis of the lattice. Represents a particle moving along either in the positive or negative direction at either the normalization velocity (r), twice the normalization velocity (2r), or three times the normalization velocity (3r). (Groups 2, 4, and 7), the three sets of eight velocities are normalized velocities (r), twice the normalized velocities (2r), or twice the normalized velocities (2r) for all three of the x, y, z lattice axes. Representing particles moving at three times the normalization velocity (3r) (groups 3, 8, and 10), 12 are of the normalization velocity for two of the x, y, z lattice axes. Representing particles moving twice (2r) (group 6), 24 are the normalization velocity (r) and twice the normalization velocity for two of the x, y, z lattice axes. Represents particles that are moving in (2r) and not moving with respect to the remaining axes (Group 5), with 24 normalizing velocities for two of the x, y, and z lattice axes (group 5). In r), and with respect to the remaining axes, represent particles moving at three times the normalizing velocity (3r) (Group 9).

In the two-dimensional model 2D-2, one of the 37 velocities represents a non-moving particle (Group 1), and three sets of four velocities are on either the x-axis or the y-axis of the lattice. Represents a particle moving along either in the positive or negative direction at either the normalizing velocity (r), twice the normalizing velocity (2r), or three times the normalizing velocity (3r) (3r). Groups 2, 4, and 7), two sets of four velocities are moving at a normalized velocity (r) or twice the normalized velocity (2r) with respect to both the x and y velocities. Representing a particle, the eight velocities are moving at a normalized velocity (r) with respect to one of the x and y lattice axes and at twice the normalized velocity (2r) with respect to the other axis. Representing a particle, the eight velocities are moving at a normalized velocity (r) with respect to one of the x and y lattice axes and at three times the normalized velocity (3r) with respect to the other axis. Represents a particle.

The LBM model described above provides a specific class of efficient and robust discrete velocity dynamics models for numerical simulation of flow in both 2D and 3D. This type of model contains a specific set of discrete velocities and their associated velocities. The velocity coincides with the grid points of Cartesian coordinates in velocity space, which facilitates accurate and efficient implementation of the discrete velocity model, in particular the kind known as the Lattice Boltzmann model. Such a model can be used to simulate flow with high fidelity.

With reference to FIG. 10, the physical process simulation system operates according to procedure 270 for simulating a physical process such as fluid flow using a CAD process as discussed above. Prior to simulation, the simulation space is modeled as a set of voxels (step 272). Generally, the simulation space is generated using a computer-aided design (CAD) program. For example, a CAD program can be used to draw a microdevice located in a wind tunnel. The data generated by the CAD program is then processed to add a lattice structure with appropriate resolution to reveal objects and surfaces in the simulation space.

The physical process simulation system operates according to procedure 270 using the modified flux calculation process discussed above. The resolution of the grid can be selected based on the Reynolds number of the simulated system. The Reynolds number is related to the viscosity of the flow (ν), the characteristic length of the object in the flow (L), and the characteristic velocity of the flow (u).
Re = uL / ν equation (I4)
Is.

The characteristic length of an object represents a large-scale feature of the object. For example, when the flow around the microdevice is simulated, the height of the microdevice can be considered as the characteristic length. If the flow around a small area of the object (eg, a car side mirror) is of interest, the resolution of the simulation can be improved, or the improved resolution area can be used around the area of interest. .. Voxel dimensions decrease as the resolution of the grid increases.

State space, f _i (x, t) is represented as, where, f _i is the number of elements or particles per unit volume of the state i in the lattice sites represented by three-dimensional vector x at time t (i.e. , The density of particles in state i). For known time increment, the number of particles are referred to simply as f _i (x). The combination of all states of the lattice site is represented 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. The number of states increases in simulations of many species.

Each state i represents a different velocity vector at a particular energy level (ie, energy levels 0, 1, or 2). The velocity c _{i of} each state is indicated by each "velocity" in three dimensions as follows.
c _i = (c _ix , c _{i ν} , c _iz ) Equation (I5)

The state of energy level 0 represents a _stopped particle that has not moved to any dimension, i.e. c _stopped = (0,0,0). The state of energy level 1 represents a particle having a velocity of ± 1 in one of the three dimensions and a velocity of 0 in the other two dimensions. The state of energy level 2 represents a particle having a velocity of ± 1 in all three dimensions, or a velocity of ± 2 in one of the three dimensions and a velocity of 0 in the other two dimensions.

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

Each voxel (ie, each grid site) is represented by a state vector f (x). The state vector fully defines the voxel state and contains 39 entries. The 39 entries are in 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. Correspond. Using this speed set, the system can produce Maxwell-Boltzmann statistics for the equilibrium state vector achieved.

For processing efficiency, voxels are grouped into 2x2x2 volumes called microblocks. Microblocks are organized to allow parallel processing of voxels and minimize the overhead associated with data structures. Shorthand voxel microblocks is defined as N _i (n), where, n is represents the relative positions of the lattice sites in the macroblock, n∈ {0,1,2, ···, 7 } in is there.

The microblock is shown in FIG.

Referring to FIGS. 12A and 12B, the surface S (FIG. 12A) is represented as a set of facets F _alpha to the simulation space (FIG. 12B),
S = {F _α } formula (I6)
Where α is an index that lists a particular facet. Facets are generally sized to be equal to or slightly smaller than the size of voxels adjacent to the facets so that the facets affect a relatively small number of voxels, rather than being restricted to voxel boundaries. To. Properties are assigned to facets to perform surface dynamics. In particular, each facet F _alpha, unit normal (n _alpha), surface area (A _alpha), the center position (x _alpha), and facet facet distribution function describing the surface kinetics of the (f _i (α)) Have. The total energy distribution function q _i (α) is treated in the same way as the flow distribution of facet-voxel interactions.

With reference to FIG. 13, different levels of resolution can be used in different regions of the simulation space to improve processing efficiency. In general, the area 320 around the object 322 is the most important and therefore simulated with the highest resolution. Since the effect of viscosity decreases with distance from the object, a reduced level of resolution (ie, the expanded voxel volume) is used to simulate regions 324, 326 separated by an increased distance from the object 322. Will be done.

Similarly, as shown in FIG. 14, a low level of resolution may be used to simulate the area 340 around the less important features of object 342, while the highest level of resolution is object 342. It is used to simulate the area 344 around the most important features of (eg, top and end surfaces). The distant region 346 is simulated using the lowest level of resolution and the largest voxels.

C. Identifying Voxels Affected by Facets With reference to FIG. 10 again, after the simulation space is modeled (step 272), voxels affected by one or more facets are identified (step 274). Voxels can be affected in several ways by facets. First, voxels with one or more facets intersecting are affected in that they are reduced in volume compared to non-intersecting voxels. This occurs because the facets and the underlying material represented by the facets occupy a portion of the voxel. Fractional factor P _f (x) indicates the portion of the voxel that is unaffected by the facet (ie, the portion that can be occupied by the fluid or other substance whose flow is being simulated). For non-intersecting voxels, P _f (x) is equal to 1.

Voxels that interact with one or more facets by transferring the particles to or receiving the particles from the facets are also identified as voxels affected by the facets. All voxels at which the facets intersect will include at least one state of receiving the particles from the facets and at least one state of transferring the particles to the facets. In most cases, additional voxels will also include such a condition.

Referring to FIG. 15, for each state i having a non-zero velocity vector c _i , the facet F _α either receives a particle from the region defined by the parallelepiped G _{i α} or transfers the particle to it and its parallelepiped G _iα is the velocity vector c _i and the vector inner product of the unit normal n _alpha facet size chromatic and height defined by a bottom surface defined by the surface area a _alpha facet (| | c _i n _i) and, as a result, the volume V _i.alpha parallelepiped G _i.alpha is
V _{i α} = | c _i n _α | A _α formula (I7)
be equivalent to.

The facets F _alpha, the velocity vector of the state is directed towards the facet _{(| c i n i | <} 0) when receives the particles from the volume V _i.alpha, the velocity vector of the state is directed away from the facet and are _{(| c i n i |>} 0) when, for transporting the particles to the region. As discussed below, this formula, another facet occupies a portion of the parallelepiped G _i.alpha, i.e., when a condition for obtaining occurring in the vicinity of the non-convex surface features such as internal corner must be modified.

Parallelepiped G _i.alpha facet F _alpha may overlap with some or all of a number of voxels. The number of voxels or parts thereof depends on the size of the facets relative to the size of the voxels, the energy of the state, and the orientation of the facets with respect to the lattice structure. The number of voxels affected increases with facet size. As a result, the size of the facet is generally chosen to be equal to or smaller than the size of the voxels placed near the facet, as described above.

Portion of the voxel N (x) which overlap parallelepiped G _i.alpha is defined as V _{iα (x).} Using this term, the flux Γ _{i α} (x) of the particles in state i moving between the voxel N (x) and the facet F _α is the density of the particles in state i in the voxel (N _i (x). )) _Is equal to the product of the volume of the region of overlap with the voxel (V _iα (x)).
Γ _iα (x) = N _i (x) + V _iα (x) Equation (I8)

When parallelepiped G _i.alpha intersects one or more facets, following conditions are true.
_{V iα = ΣV α (x)} + ΣV iα (β) formula (I9)
Here, the first sum corresponds to all voxels in which G _iα overlaps, and the second term corresponds to all facets intersecting G _iα . If parallelepiped G _i.alpha does not intersect with another facet, the formula
_V iα = ΣV iα (x) formula (I10)
Is.

D. Running Simulation After the voxels affected by one or more facets have been identified (step 274), a timer is initialized to start the simulation (step 276). During each time increment of the simulation, the movement of the particles from voxel to voxel is simulated by an advection step (steps 278-286) that causes the interaction of the particles with the surface facets. The collision step (step 288) then simulates the interaction of particles within each voxel. The timer is then incremented (step 290). If the incremented timer does not indicate the completion of the simulation (step 292), the advection and collision steps (steps 278-290) are repeated. If the incremented timer indicates the completion of the simulation (step 292), the simulation results are stored and / or displayed (step 294).

1. 1. Surface Boundary Conditions Each facet must meet four boundary conditions in order to correctly simulate surface interaction. First, the total mass of particles received by the facet must be equal to the total mass of particles transferred by the facet (ie, the net mass flux to the facet must be equal to 0). Second, the total energy of the particles received by the facet must be equal to the total energy of the particles transferred by the facet (ie, the net energy flux to the facet must be equal to zero). These two conditions can be met by requiring that the net mass flux at each energy level (ie, energy levels 1 and 2) be equal to zero.

The other two boundary conditions relate to the net momentum of the particles interacting with the facets. On a surface without surface friction, referred to herein as a slip surface, the net tangential momentum flux must be equal to zero and the net normal momentum flux must be equal to the local pressure of the facet. Therefore, the component of the total momentum received and the total momentum transferred (ie, the tangent component), which is perpendicular to the facet normal n _α , must be equal, while parallel to the facet normal n _α. The difference between the total momentum received and the component of the total momentum transferred (ie, the normal component) must be equal to the local pressure at the facet. On 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 friction.

2. 2. Voxel-to-facet collection As a first step in simulating the interaction between particles and surfaces, particles are collected from the voxels and fed to the facets (step 278). As described 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)
Is.

From this equation, for each state directed to the facet _{F α i (c i n α} <0), the number of voxels facet F _alpha to supply particles,
Γ _{iα V → F} = Σ _X Γ _{i α} (x) = Σ _X N _i (x) V _{i α} (x) Equation (I12)
Is.

Only voxels with a non-zero value for V _iα (x) must be summed. As mentioned above, the facet size is chosen such that V _iα (x) has a non-zero value only in a small number of voxels. Since V _{i α} (x) and P _f (x) can have non-integer values, Γ _α (x) is stored and processed as a real number.

3. 3. Moving from facet to facet Next, particles are moved between facets (step 280). When the parallelepiped G _{i α in} the inflow state of facet F _α (c _i n _α <0) intersects another facet F _β , some of the particles in state i received by facet F _α are from facet F _β. Will come. In particular, the facet F _α will receive some of the particles in state i created by the facet F _β during the previous time increment.

Next, with reference to FIG. 17, the relationship between the particles in state i created by facet F _β during the previous time increment is shown. 17, a portion 380 of the parallelepiped G _i.alpha intersecting the facet F _beta is, it is shown equal to the part 382 of the parallelepiped G _i.beta intersecting the facet F _alpha. As mentioned above, the intersection is represented as _Viα (β). Using this term, the flux of particles in state i between facet F _β and facet F _α is
_{Γ iα (β, t-1} ) = Γ i (β) V iα (β) / V iα formula (I.13)
Where Γ _i (β, t-1) is the amount of particles in state i produced by facet F _β during the previous time increment. The number from this equation, facet _{F α (c i n α <} 0) for each state i, which is directed to, particles fed to the facet F _alpha by another facet,
Γ _{iα F → F} = Σ _β Γ _{i α} (β) = Σ _β Γ _i (β, t-1) V _{i α} (β) / V _{i α} equation (I.14)
And the total flux of the particles of state i to 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)
Is.

The facet state vector N (α), also called the facet distribution function, has M entries corresponding to the M entries of the voxel state vector. M is the number of discrete lattice velocities. Input state of the facet distribution function N (alpha) is set equal to those obtained by dividing the flux of particles to their state in the volume V _i.alpha, against c _i n _{α <0,}
N _i (α) = Γ _iIN (α) / V _{i α} equation (I.16)
Is.

The facet distribution function is a simulation tool for generating the output flux from the facet and does not necessarily represent the actual particle. Assign values to other states of the distribution function to generate an accurate output flux. Using the above techniques for introducing inward state, outwards state is turned on, the _{c i} n α ≧ 0,
_{N i (α) = Γ iOTHER} (α) / V iα formula (I.17)
Here, Γ _iOTHER (α) uses the above technique for generating Γ _iIN (α), but states other than the inflow state (c _i n _α <0) (c _i n _α ≧). It is determined by applying this technique to 0). In an alternative approach, Γ _iOTHER (α) can be generated using the value of Γ _iOUT (α) from the previous time step, and as a result,
Γ _iOTHER (α, t) = Γ _iOUT (α, t-1) Equation (I.18)
Will be.

In the parallel state (c _i n _α = 0), both V _{i α} and V _{i α} (x) are 0. In the equation of N _i (α), V _{i α} (x) appears in the numerator (from the equation of Γ _iOTHER (α)) and V _i α appears in the denominator (from the equation of N _i (α)). As a result, the parallel state N _i (alpha) is determined as a limit of N _i (alpha) when V _i.alpha and V _i.alpha (x) approaches zero. The values of the states with zero velocity (ie, the stationary state, and the states (0,0,0,2) and (0,0,0,-2)) are the beginning of the simulation based on the initial conditions of temperature and pressure. Is initialized to. These values are then adjusted over time.

4. Execution of facet surface dynamics Next, surface dynamics are executed for each facet so as to satisfy the four boundary conditions discussed above (step 282). The procedure for performing surface dynamics on the facet is shown in FIG. First, for all i, the total momentum perpendicular to the facet F _α determines the total momentum P (α) of the particles in the facet.

として決定される（ステップ３９２）。この式から、法線運動量Ｐ_n（α）は、
Ｐ_n（α）＝ｎ_α・Ｐ（α） 式（Ｉ．２０）
として決定される。

Is determined as (step 392). From this equation, the normal momentum P _n (α) is
P _n (α) = n _α · P (α) equation (I.20)
Is determined as.

This normal momentum is then removed using a push / pull technique (step 394) to produce N _n- (α). By this technique, particles are moved between states so that they only affect normal momentum. The push / pull technique is described in US Pat. No. 5,594,671 which is incorporated by reference.

Then, the particles of N _n- (α) are collided to create a Boltzmann distribution N _n-β (α) (step 396). The Boltzmann distribution can be achieved by applying a set of collision rules to N _n- (α), as described below for the practice of fluid dynamics.

Effluent flux distribution of facet F _alpha is flowing flux distribution, corrected flux calculation to take into account the CFL constraint violations, and on the basis of the Boltzmann distribution is determined (step 398).

First, the difference between the inflow flux distribution Γ _i (α) and the Boltzmann distribution is
ΔΓ _i (α) = Γ _iIN (α) −N _{n-β i} (α) V _{i α} equation (I.21)
Is determined as.

Using this difference, the outflow flux distribution
For n _α c _i > 0
Γ _iOUT (α) = N _{n-β i} (α) V _{i α} −. Δ. Γ _{i *} (α) Equation (I.22)
Here, i ^* is a state having a direction opposite to that of the state i. For example, when the state i is (1,1,0,0), the state i ^* is (-1, -1,0,0). To take into account surface friction and other factors, the outflow flux distribution can be further refined as follows. For n _α c _i > 0
_{Γ 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α formula (I.23)
, And the Here, C _f, is a function of surface friction, t _i.alpha is, n is a first tangent vector perpendicular to the _alpha, t _2.alpha is, n _alpha and t second both to a vertical _1α The tangent vector of, ΔN _{j, 1} and ΔN _{j, 2} are distribution functions corresponding to the state i and the energy (j) of the displayed tangent vector. The distribution function is

に従って決定され、ここで、ｊは、エネルギーレベル１の状態では１に等しく、エネルギーレベル２の状態では２に等しい。

Here, j is equal to 1 in the state of energy level 1 and equal to 2 in the state of energy level 2.

The role of each term in the equation of Γ _iOUT (α) is as follows. The first and second terms enforce boundary conditions for normal momentum flux to the extent that collisions are effective in creating a Boltzmann distribution, but include tangential momentum flux anomalies. The fourth and fifth terms correct for this anomaly that may occur due to discrete effects or non-Boltzmann structures due to inadequate collisions. Finally, the third term adds a specific amount of surface fraction to force the desired change in the tangential momentum flux of the surface. The generation of the friction coefficient C _f will be described below. Note that all terms, including vector manipulation, are geometric factors that can be calculated before starting the simulation.

From this equation, the tangential velocity is
u _i (α) = (P (α) -P _n (α) n _α ) / ρ equation (I.25)
Determined as, where ρ is the density of the facet distribution

である。

Is.

As before, the difference between the inflow flux distribution and the Boltzmann distribution is
ΔΓ _i (α) = Γ _iIN (α) −N _{n-β i} (α) V _{i α} equation (I.27)
Is determined as.

Next, the outflow flux distribution is
_{Γ iOUT (α) = N n} -βi (α) V iα -ΔΓ i * (α) + C f (n α c i) [N n-βi * (α) -N n-βi (α)] V iα Equation (I.28)
This corresponds to the first two lines of the outflow flux distribution determined by previous techniques, but does not require correction to anomalous tangential flux.

Regardless of which method is used, the resulting flux distribution meets all of the momentum flux conditions, i.e.

であり、ここで、ｐ_αは、ファセットＦ_αの平衡圧力であり、ファセットに粒子を供給するボクセルの平均密度および温度値に基づき、ｕ_αは、ファセットにおける平均速度である。

Here, p _α is the equilibrium pressure of facet F _α , and u _α is the average velocity at facet, based on the average density and temperature values of the voxels that feed the facets.

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

として測定され、ここで、インデクスｊは、状態ｉのエネルギーを表す。次いで、このエネルギー差を使用して、差の項を生成し、ｃ_jiｎ_α＞０に対して、

Where the index j represents the energy of the state i. This energy difference is then used to generate a term for the difference, for c _j n _α > 0.

である。この差の項を使用して、ｃ_jiｎ_α＞０に対して、
Γ_αjiOUTf＝Γ_αjiOUT＋δΓ_αji 式（Ｉ．３２）
になるように流出流束を修正する。この操作は、接線運動量流束を不変のままにしながら、質量およびエネルギー流束を補正する。この調節は、流れがファセットの近傍でほぼ一様で平衡に近い場合、わずかである。調節の後に結果として生じる法線運動量流束は、近傍の平均特性に、近傍の非均一または非平衡特性による補正を加えたものに基づいて、平衡圧力である値までわずかに変更される。ＣＦＬ制約が違反されている場合、プロセスは、修正流束計算手法を、図１０のプロセスに含まれるいずれかの流束計算に適用する２８５。

Is. Using this difference term, for c _ji n _α > 0,
Γ _αjiOUTf = Γ _αjiOUT + δΓ _αji equation (I.32)
Correct the outflow flux so that This operation corrects the mass and energy flux while leaving the tangential momentum flux unchanged. This adjustment is slight when the flow is nearly uniform and close to equilibrium near the facets. The resulting normal momentum flux after adjustment is slightly modified to the equilibrium pressure value based on the mean characteristic of the neighborhood plus correction by the non-uniform or non-equilibrium characteristic of the neighborhood. If the CFL constraint is violated, the process applies the modified flux calculation method to any of the flux calculations included in the process of FIG. 10 285.

5. Moving from voxel to voxel With reference to FIG. 10 again, the particles are moved between voxels along a three-dimensional linear lattice (step 284). This voxel-to-voxel movement is the only movement operation performed on voxels that do not interact with the facets (ie, voxels that are not located near the surface). In a typical simulation, voxels that are not placed close enough to interact with the surface make up the majority of voxels.

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

For voxels that do not interact with the surface, the movement operation is computationally very simple. The entire population of states is moved from the current voxel to the destination voxel in time increments. At the same time, particles of the destination voxel are moved from that voxel to its own destination voxel. For example, energy level 1 particles moving in the + 1x and + 1y directions (1,0,0) are moved from the current voxel to a voxel at +1 in the x direction and 0 in the other direction. The particles eventually arrive at the destination voxel in the same state they had before the move (1,0,0). Interactions within voxels will change the number of particles for that state based on local interactions with other particles and surfaces. Otherwise, the particles will continue to move along the grid at the same velocity and direction.

The move operation is slightly more complicated for voxels that interact with one or more surfaces. This will transfer one or more fractionated particles to the facet. By transferring such fractionated particles to the facet, the fractionated particles will remain in the voxels. These fractionated particles are transferred to voxels occupied by facets.

Referring to FIG. 16, when a part 360 of the particles in the state i of the voxel 362 is moved to the facet 364 (step 278), the remaining part 366 is moved to the voxel 368 where the facet 364 is located. From there, the particles in state i are directed to the facet 364. Therefore, the state population is equal to 25, V _i.alpha (x) is equal to 0.25 (i.e., a quarter of voxel intersects the parallelepiped G _i.alpha) case, 6.25 pieces of particles facet F _alpha It would be moved to, so that 18.75 pieces of particles are moved to the voxels occupied by facet F _alpha. Since many facets intersect a single voxel, the number of particles in state i transferred to voxel N (f) occupied by one or more facets is

であり、ここで、Ｎ（ｘ）はソースボクセルである。

Where N (x) is the source voxel.

6. Scattering from facets to voxels Next, outflow particles from each facet are scattered to voxels (step 286). In essence, this step is the reverse of the collection step in which particles are moved from voxels to facets. The number of particles in state i moving from facet F _α to voxel N (x) is

であり、ここで、Ｐ_f（ｘ）は、部分ボクセルの体積減少に相当する。この式から、状態ｉごとに、ファセットからボクセルＮ（ｘ）に向けられる粒子の総数は、

Where P _f (x) corresponds to the volume reduction of the partial voxels. From this equation, the total number of particles directed from facets to voxels N (x) for each state i is

である。

Is.

After scattering particles from facets to voxels, combining them with particles advected from surrounding voxels, and integerizing the results, certain directions of a particular voxel either underflow (negative) or become negative. It can overflow (more than 255 for 8-bit implementations). This will result in either an increase or decrease in mass, momentum, and energy after these amounts have been truncated to fit acceptable values. To prevent the occurrence of such phenomena, off-boundary mass, momentum, and energy are stored before truncating the violating state. An amount of mass equal to an increased value (due to underflow) or a decreased value (due to overflow) with respect to the energy to which the state belongs has the same energy and is itself affected by overflow or underflow. Not received Added to randomly (or sequentially) selected states. The additional momentum resulting from this mass and energy addition is accumulated and added to the momentum due to 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 the push / pull technique until the momentum accumulator returns to zero.

7. Execution of Fluid Dynamics Execution of fluid dynamics (step 288), FIG. This step can also be called microdynamics or intra-voxel manipulation. Similarly, the advection procedure can also be referred to as an intervoxel operation. The microdynamic operations described below can also be used to collide particles with facets to create a Boltzmann distribution.

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 the distribution of an actual fluid system. The collision process can be fully described by the right side of Equations 1 and 2. After the advection step, the conserved amount of the fluid system, specifically the density, momentum, and energy, is obtained from the distribution function using Equation 3. From these quantities, the equilibrium distribution function represented by f ^eq in equation (2) is completely defined by equation (4). The velocity vector set c _i , weights, and selection of both are listed in Table 1, and together with Equation 2, ensure that the macroscopic behavior follows the correct hydrodynamic equation.

Variable resolution Variable resolution (as discussed in US Patent Application Publication No. 2013/0151221) may also be used, resulting in the use of voxels of different sizes, such as coarse and fine voxels. ..

By leveraging the intrinsic transient lattice Boltzmann-based physics, the system can perform simulations that accurately predict real-life conditions. For example, engineers evaluate product performance early in the design process before the prototype is built, if the impact of the change is of paramount importance to design and cost. The system uses CAD geometry to perform aerodynamic, air-acoustic, and thermal management simulations accurately and efficiently. The system can run simulations to address the following uses: Aerodynamics (aerodynamic efficiency; vehicle operation; dirt and water management; panel deformation; driving dynamics), aerodynamics (greenhouse wind noise; bottom wind noise; gap / seal noise; mirror, whistle, and tonal noise; sunroof and windows Buffeting; Passing / Urban noise; Cooling fan noise), Thermal management (Cooling airflow; Thermal protection; Brake cooling; Operation cycle simulation; Key-off and soak; Electronic equipment and battery cooling; ROA / Intake port), Air conditioning (Indoor comfort) Gender; HVAC unit & distribution system performance; HVAC system and fan noise; defrosting and thawing), power train (drive train cooling; exhaust system; cooling jacket; engine block), dirt and water management (pillar overflow, dust and dust accumulation) , Tire spray) etc.

Embodiments of the subject matter and functional operation described herein include digital electronic circuits, tangibly embodied computer software or firmware, computer hardware, including the structures and structural equivalents disclosed herein. It can be realized with or in combination of one or more of these. Embodiments of the subject matter described herein are tangible non-temporary program carriers for execution by one or more computer programs (ie, for execution by a data processor or for controlling the operation of a data processor). It can be realized as one or more modules of computer program instructions encoded in. The computer storage medium can be a machine-readable storage device, a machine-readable storage substrate, a random or serial access memory device, or a combination of one or more thereof.

The term "data processor" refers to data processing hardware, eg, programmable processors, computers, or all types of devices, devices, and machines for processing data, including multiple processors or computers. Include. The device can be a dedicated logic circuit (eg, FPGA (Field Programmable Gate Array) or ASIC (Application Specific Integrated Circuit)), or can further include it. In addition to the hardware, the device contains the code that creates the execution environment for the computer program (eg, the code that makes up the processor firmware, protocol stack, database management system, operating system, or one or more combinations thereof). Can be included as an option.

Computer programs, which can also be referred to or written as programs, software, software applications, modules, software modules, scripts, or code, are any form of programming language, including compiled or interpreted languages, or declarative or procedural languages. It can be written in and deployed in any form, including as a stand-alone program or as a module, component, subroutine, or separate unit suitable for use in a computing environment. Computer programs can, but do not have to, handle files in the file system. A program can be another program or data (eg, in a markup language document, in a single file dedicated to that program, or in a number of collaborative files (eg, one or more modules, subprograms, or code). It can be stored in a part of the file that holds one or more scripts) stored in the file that stores the part. A computer program is such that the program runs on one computer, is located at one site, or is distributed across many sites and runs on many computers interconnected by a data communication network. Can be deployed.

A computer suitable for executing a computer program can be based on a general purpose or dedicated microprocessor or both, or other types of central processing units. Computer-readable media suitable for storing computer program instructions and data include, for example, semiconductor memory devices (eg, EPROM, EEPROM, and flash memory devices), magnetic disks (eg, internal hard disks or removable disks), magneto-optical. Includes discs and all forms of non-volatile memory on media and memory devices, including CD-ROMs and DVD-ROM disks. The processor and memory may be supplemented by or incorporated into a dedicated logic circuit.

Embodiments of the subject matter described herein include back-end components (eg, as data servers), middleware components (eg, application servers), or front-end components (eg, application servers). For example, in a computing system that includes a graphical user interface or web browser that allows the user to interact with embodiments of the subject matter described herein, or in one or more such backs. It can be implemented with any combination of end, middleware, or front end components. The components of the system can be interconnected by digital data communication (eg, a communication network) of any form or medium. Examples of communication networks include local area networks (LANs) and wide area networks (WANs) (eg, the Internet).

The computing system can include clients and servers. Clients and servers are usually separated from each other and generally interact through a communication network. The client-server relationship arises from computer programs that run on their respective computers and have a client-server relationship with each other. In some embodiments, the server sends data (eg, an HTML page) to a user device that acts as a client (eg, to display data to a user interacting with the user device and to receive user input from the user). To do. Data generated on the user device (eg, the result of user interaction) can be received from the user device on the server.

Specific embodiments of the subject were described. Other embodiments are within the scope of the following claims. For example, the actions detailed in the claims can be performed in a different order and still achieve the desired result. As an example, the process shown in the attached figure does not necessarily require the specific order or sequential order shown to achieve the desired result. In some cases, multitasking and parallel processing may be advantageous.

10 System 12 Server System 14 Client System 18 Memory 20 Interface 22 Transfer Operation 24 Processing Device 28 Mesh Definition 30 Simulation Process 32 Mesh Preparation Engine 34 Simulation Engine 34a Particle Collision Interaction Module 34b Particle Boundary Model Module 34c Transfer Module 36 Submodule, Engine 38 Data repository

Claims (14)

A computer-implemented method for simulating fluid flow around the surface of a solid, said method.
It is a step of receiving a model of a simulation space containing a lattice structure represented as a set of voxels and a representation of a physical object by one or more computing systems, wherein the voxels cover the surface of the physical object. The receiving step and the receiving step, which have the appropriate resolution to explain,
A step of simulating the movement of particles in a volume of fluid by the one or more computer systems, wherein the movement of the particles causes collisions between the particles.
A step of identifying a face between two voxels by the computing system, and a step of identifying at least one of the faces violating stability conditions.
A step of calculating a modified heat flux by the computing system using a spatially averaged temperature gradient in the vicinity of the two voxels, wherein at least one of the surfaces violates the stability condition. There are steps to calculate and
A computer implementation method comprising the step of performing an advection operation of the particles into a subsequent voxel by the computing system. The computer method of claim 1, further comprising calculating a temperature gradient calculated based on a temperature difference form on adjacent surfaces where neither of the surfaces violates the stability conditions. The temperature difference form determines the net energy transfer to each of the elements corresponding to the sum of the energy transfers from all aspects of the element in order to calculate the final temperature of the element at the end of the time step. Item 2. The computer method according to item 2. The step of calculating the modified heat flux is
The step of calculating the additional flux by the computing system,
The computer method of claim 1, further comprising the step of calculating the equilibrium flux by the computing system. The computer method of claim 4, wherein the calculated additional flux is used for the temperature evolution of the element under consideration. The computer method of claim 4, wherein the equilibrium flux is used for temperature development depending on the size of the element. The computer method of claim 6, wherein the equilibrium flux is used for the temperature development if the size of the element is small enough to violate the constraint. The computer method of claim 4, further comprising sending the equilibrium flux to one or more adjacent elements in the direction of the flux by the computer system. The equilibrium flux is continuously delivered along the direction of the flux until the equilibrium flux is finally transferred to a sufficiently large element, and the equilibrium flux is at the sufficiently large element. The computer method according to claim 4, which is applied toward the temperature development of the above. The computer method of claim 9, wherein if the element is large enough, the equilibrium flux is used for temperature evolution. The computer method according to claim 1, wherein the conditions include stability characteristics. The computer method of claim 1, wherein the condition is a Courant-Friedrichs-Lewy (CFL) constraint of a time progression scheme that determines the maximum time step size that can be used to maintain a stable distribution. A device for simulating a fluid flow around the surface of a solid, said device.
With memory
A device comprising one or more processor devices configured to perform the function according to any one of claims 1-10. On the computer
On one or more machine-readable hardware storage devices, including executable instructions that can be executed by one or more processing devices to perform the function according to any one of claims 1-10. A tangible computer program product stored in.

Images