Computer simulation of physical fluids on a grid in an arbitrary coordinate system

Abstract

Disclosed is a computer simulation of a physical fluid on a grid in an arbitrary coordinate system. Specifically disclosed is a computer-implemented technique for simulating fluid flow near a solid surface, comprising: receiving a coordinate system for representing a curvilinear grid conforming to the solid surface; simulating the transport of particles in a volume of the fluid using a lattice velocity set, the transport causing collisions between particles; executing a distribution function for the transport of particles, the distribution function including particle collision determinations and changes in the particle distribution associated with the curvilinear grid; performing, by a computing system, an advection operation in the coordinate system under constraints applied to particle momentum values; and mapping, by a computing system, values ​​obtained from the simulation onto the curvilinear grid by converting particle momentum values ​​and spatial coordinates determined in the coordinate system to momentum and spatial values ​​in curvilinear space.

Description

Background Art

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

[0002] The so-called "Lattice Boltzmann Method" (LBM) is an advantageous technique used in computational fluid dynamics. The underlying dynamics of a lattice Boltzmann system lies in the basic physics of kinetic theory, which involves the motion of many particles according to the Boltzmann equation. There are two basic dynamical processes in the basic Boltzmann dynamical system - collision and advection. The collision process involves the interaction between particles that obey conservation laws and relax to equilibrium. The advection process involves modeling the movement of particles from one location to another according to their microscopic velocities.

[0003] In the standard LBM model, particle velocities have a set of discrete constant values ​​that form exact links from one lattice site to its neighboring sites on a simple Bravais lattice corresponding to a three-dimensional (3D) uniform cubic Cartesian grid.

[0004] Various attempts have been made to extend LBM to grids based on arbitrary coordinate systems (arbitrary (one or more) grids). One of the main representative methods is to relax the one-to-one advection mapping between a pair of lattice sites. On such arbitrary grids, particles after advection from their original grid sites usually do not land on a single neighboring site. In those solutions, the position of the particle is represented by those grid sites through interpolation. Summary of the invention

[0005] While important in many applications, uniform cubic Cartesian grids pose fundamental limitations for others. For example, in realistic fluid simulations, curved surfaces of solid geometry are often simulated. Cartesian grids do not exhibit smooth conformance to curved surfaces of solid geometry. Furthermore, realistic physical flows often have small structures concentrated in certain spatial regions and directions. For example, in a so-called turbulent boundary layer, the proportion of flow perpendicular to the wall is much smaller than in the tangential direction or in the bulk of the fluid region. Therefore, the spatial resolution requirement is significantly higher in the direction normal to the wall inside the boundary layer. Cubic Cartesian grids do not provide the flexibility to have different spatial resolutions in different directions.

[0006] According to one aspect, a computer-implemented method for simulating fluid flow near a solid surface comprises: receiving, by a computing system, a coordinate system for representing a curvilinear grid conforming to the solid surface; simulating, using a lattice velocity set, the transport of particles in a volume of a fluid, the transport causing collisions between particles; executing a distribution function for the transport of particles, the distribution function including particle collision determinations and changes in the particle distribution associated with the curvilinear grid; performing, by the computing system, an advection operation in the coordinate system under constraints applied to particle momentum values; and mapping, by the computer system, values ​​obtained from the simulation onto the curvilinear grid by converting particle momentum values ​​and spatial coordinates determined in the coordinate system into momentum and spatial values ​​in the curvilinear space.

[0007] Other aspects include a computer program product, one or more machine-readable hardware storage devices, an apparatus, and a computing system.

[0008] The method disclosed herein extends the current LBM-based simulation predicted on a Cartesian grid to a non-Cartesian grid framework in a curved space. The method disclosed herein is based on volume equations so that conservation of mass and momentum is satisfied. The resulting mass continuity equation will have the correct form in curvilinear coordinates, so the method does not need to introduce any artificial mass source terms to correct artifacts in the resulting fluid dynamics. In addition, as in the continuous dynamics theory on manifolds, in the disclosed extended LBM, the only external source term is the presence of inertial forces due to the curved space. This inertial force term does not contribute to mass and is constructed without relying on the analytical form in the continuous dynamics theory.

[0009] This inertial force term in discrete space and time asymptotically recovers the force in the continuous dynamic theory in the hydrodynamic limit. In addition, in the discrete spacetime LBM model, the inertial force enforces exact momentum conservation for the underlying Euclidean space. The force term is constructed so that it adds momentum to the system with the appropriate discrete time moments to produce the correct resulting Navier-Stokes hydrodynamics of the viscous order.

[0010] Other features and advantages will become apparent from the following description, including the drawings and claims. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 A system for simulating fluid flow is depicted.

[0012] Figure 2 A flow chart showing operations for formulating a lattice Boltzmann model in curvilinear coordinates is depicted.

[0013] Figure 3A flow chart illustrating operations using a stored curvilinear mesh that corresponds to a simulated physical object is depicted.

[0014] Figure 4 A flow chart showing simulation operations using a lattice Boltzmann model expressed in curvilinear coordinates is depicted.

[0015] Figure 5 A flow chart showing constraint operations applied to a distribution function is depicted.

[0016] Figure 6A-6B Representations of physical objects with shapes expressed in curvilinear and Cartesian spaces are depicted.

[0017] Figure 7A-7B A representation of the movement of particles in curved and Cartesian space is depicted.

[0018] Figures 8A-8B Depicted is an alternative representation of the movement of particles in curved and Cartesian space.

[0019] Fig. 9 and Fig.10 The velocity components of two LBM models represented in non-Euclidean space are illustrated.

[0020] Fig.11 It is a flow chart of the process followed by the physical process simulation system.

[0021] Fig.12 This is a perspective view of a microblock.

[0022] Figures 13A-13B yes Figure 1 Illustration of the lattice structure in a non-Euclidean space used by the system.

[0023] Fig.14 and 15 Variable resolution technology is illustrated.

[0024] Fig.16 The movement of particles is shown.

[0025] Fig.17 The area affected by the facets of the surface is illustrated.

[0026] Fig.18 A flow chart of surface dynamics is illustrated.

[0027] Fig.19 Illustration of the movement of particles from a voxel to a surface.

[0028] Fig. 20 is a flow chart of the process used to perform surface dynamics.

[0029] Fig.21 is a screenshot representing a fluid simulation. DETAILED DESCRIPTION

[0030] Modeling the Simulation Space

[0031] In the physical process simulation system based on LBM, the fluid flow can be represented by a set of discrete velocities c i The distribution function value f to be evaluated i The dynamics of the distribution function is governed by equation (I.1),

[0032] f i (x+c i , t+1) = f i (x,t)+C i (x, t) Equation (I.1)

[0033] This equation describes the distribution function f i The well-known lattice Boltzmann equation for the time evolution of . The left-hand side shows the changes in the distribution due to the so-called "fluidization process". A fluidization process is when a piece of fluid starts at a grid location and moves along one of the velocity vectors to the next grid location. At that point, the "collision operator" is calculated, i.e., the influence of the nearby fluid pieces on the starting fluid piece. The fluid can only move to one other grid location, so an appropriate choice of the velocity vectors is necessary so that all components of all velocities are multiples of a common speed.

[0034] The right hand side of the first equation is the "collision operator" mentioned above, which represents the change in the distribution function due to collisions between fluid blocks. The specific form of the collision operator used here can be, but is not limited to, Bhatnagar, Gross and Krook (BGK). It forces the distribution function to go to a prescribed value given by the second equation, where the second equation is the "equilibrium" form,

[0035]

[0036] The BGK operator is based on the fact that (regardless of the details of the collision) the distribution function is approached by {f eq (x,v,t)} given by the physical point of view of a well-defined local equilibrium:

[0037]

[0038] The parameter τ represents the characteristic relaxation time for reaching equilibrium via collisions.

[0039] From this simulation, conventional fluid variables, such as mass p and fluid velocity u, are obtained as simple summations in equation (I.3) below.

[0040] Due to symmetry considerations, the set of velocity values ​​is chosen in such a way that they form a certain lattice structure when spanned in the configuration space. The dynamics of this discrete system follows the LBM equations of the following form:

[0041] f i (x+c i , t+1) = f i (x,t)+C i (x, t)

[0042] where the collision operator usually takes the BGK form as described above. By appropriate choice of the equilibrium distribution form, it can be shown theoretically that the lattice Boltzmann equation produces correct hydrodynamic and thermo-hydrodynamic results. That is, from f i The (x, t) derived hydrodynamic moments obey the Navier-Stokes equations in the macroscopic limit. These moments are defined as:

[0043] ρ(x, t) = ∑ i f i (x,t);ρ(x,t)u(x,t)=∑ i c i f i (x, t) Equation (I.3)

[0044] where ρ and u are the fluid density and fluid velocity, respectively.

[0045] c i and w i The LBM model can be efficiently implemented on scalable computer platforms and runs with maximum robustness for time-invariant flows and complex boundary conditions.

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

[0047] For further description of a conventional LBM-based physical process simulation system, reference is made to U.S. Patent Publication No. US-2016-0188768-A1, the entire contents of which are incorporated herein by reference.

[0048] refer to Figure 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 implemented as a massively parallel computing system 12 (standalone or cloud-based) and a client system 14. The server system 12 includes a memory 18, a bus system 11, an interface 20 (e.g., a user interface / network interface / display or monitor interface, etc.), and a processing device 24. In the memory 18 is a grid preparation engine 32 and a simulation engine 34.

[0049] Although Figure 1 A grid preparation engine 32 is shown in memory 18, but the grid preparation engine may be a third-party application executing on a different system than the server 12. Whether the grid preparation engine 32 is executed in memory 18 or on a different system than the server 12, the grid preparation engine 32 receives a user-provided grid definition 30, and the grid preparation engine 32 prepares the grid and (stores or) sends the prepared grid to a simulation engine 34. The simulation engine 34 includes a particle collision interaction module, a particle boundary modeling module, and advection operations. The system 10 accesses a data repository 38 that stores 2D and / or 3D grids (Cartesian and / or curvilinear), coordinate systems, and libraries.

[0050] Reference now Figure 2 , a process 40 for simulating fluid flow near a representation of a physical object is shown. In this example, the physical object is an airfoil. The use of an airfoil is merely illustrative, as the physical object can be of any shape, and in particular has (one or more) curved surfaces. The process receives 42 a mesh, such as a curvilinear mesh that conforms to the physical object being simulated, for example, from a client system 14 or by retrieval from a data repository 38. In other embodiments, an external system or server 12 based on user input generates a curvilinear mesh that conforms to the physical object being simulated. The process precomputes 44 geometric quantities based on the retrieved curvilinear mesh, and performs 46 a dynamic lattice Boltzmann model simulation using the precomputed geometric quantities corresponding to the retrieved curvilinear mesh.

[0051] Reference now Figure 3 , aspects of the construction 50 of the curvilinear lattice and the determination of the pre-calculated geometric quantities are explained below. A curvilinear coordinate system (corresponding to the physical object to be simulated) is selected so that x=x(q) is uniquely defined by q, where q≡(q 1 ,q 2 ,q 3 ) are the coordinates of three non-co-linear congruencies along the parameterized basis curve.

[0052] The system 10 starts 52 with a 3D curvilinear mesh, where the curvilinear mesh has a layout in a given space such that the positions of all its vertices in space are known and specified. Let x be any spatial point in 3-dimensional (3D) Euclidean space. A spatial point x is defined 54 only at a site (i.e., a vertex) of the curvilinear mesh. Once the curvilinear mesh is given (i.e., the "original" curvilinear mesh {x}), all spatial positions x on the mesh {x} are known. i .

[0053] The system 10 selects 56 a coordinate system {q} that is a one-to-one mapping between a point x in space and a site (q) on the curvilinear grid. For any site x on the curvilinear grid {x}, there is a unique value q associated with x, i.e., x=x(q), such that x is uniquely defined. The system 10 constructs the coordinate values ​​of q on the grid as follows:

[0054] For any site x(q), on the curve grid, the nearest neighboring site along the αth (α=1,2,3) coordinate curve in the positive or negative direction is the spatial point x ±i Since this unique mapping x ±α = x(q ±α ), where q ±α are unique coordinate values ​​of adjacent sites, so it is entirely possible to choose a spatial variation of q such that q ±α The αth coordinate component of q differs only by a constant distance (d0), and where (α, β = 1, 2, 3). In this example, the constant d0 is chosen to be the unit element (1) without loss of generality.

[0055] In this construction, the coordinate values ​​{q} provide a simple uniform 3D cubic Cartesian lattice structure with a lattice spacing of unit element (value d0). This variant "Cartesian" lattice {q} is produced by a deformation (bending, twisting and stretching / compressing) of the original curvilinear grid {x} in Euclidean space. Thus, the topology of the "Cartesian" lattice {q} is the same as the original curvilinear grid {x}, but the resulting Cartesian lattice is on a non-Euclidean space.

[0056] When a curved mesh is provided, the spatial locations of all vertices {x} on the mesh are specified, and the distance from any vertex to any other vertex on the curved mesh is also fully determined. The distance from a site x(q) to one of its neighbors x(q ±α )(α=1;2;3) distance vector D ±α (q) is defined as:

[0057] D ±α (q)≡x(q±α )-x(q);α=1,2,3 (Equation 1)

[0058] For example, D ±1 (q)≡x(q 1 ±1,q 2 ,q 3 )-x(q 1 ,q 2 ,q 3 ).

[0059] Due to the spatial inhomogeneity of general curvilinear grids, the spatial distance from a grid site to its nearest neighbor site usually varies with position. In other words, D ±α (q) is a function of q. In addition, the distance value in the positive direction along the αth coordinate curve is usually not equal to the distance value in the negative direction.

[0060] Specifically, in terms of the distance vector, D α (q)≠-D -α (q). For example, according to the definition given by (Equation 1),

[0061] D -1 (q) = x(q 1 -1,q 2 ,q 3 )-x(q 1 ,q 2 ,q 3 )=-D1(q 1 -1,q 2 ,q 3 )≠-D1(q)=-(x(q 1 +1,q 2 ,q 3 )-x(q 1 ,q 2 ,q 3 ))

[0062] Equation 2

[0063] This inequality holds only if the curvilinear grid is a uniform Cartesian lattice (such that |D α |=constant, independent of the spatial coordinate value q) becomes equal everywhere.

[0064] The system 10 may store 58 the constructed coordinate system for the selected curvilinear grid.

[0065] Based on the concepts of basic differential geometry, the system 10 constructs 60 tangent basis vectors along each coordinate direction

[0066] g α (q)≡[D α (q)-D-α (q)] / 2; α=1,2,3 (Equation 3)

[0067] Through this construction, even-odd symmetry is achieved, so that g α (q) = -g -α (q). Unlike the Cartesian coordinate system in Euclidean space, the basis tangent vector g of the curvilinear coordinate system is α (q)(α=1,2,3) are usually not orthogonal. That is,

[0068] Therefore, the system 10 constructs 62 a corresponding metric tensor, which is defined based on the above-mentioned basis tangent vector as:

[0069] g αβ (q)≡g α (q)·g β (q),α,β=1,2,3 (Equation 4)

[0070] and the volume J of a cell centered at x(q) is defined as:

[0071] J(q)≡(g1(q)xg2(q))·g3(q) (Equation 5)

[0072] And choose an appropriate "handedness" so that J(q)>0, where J(q) is a constant in the space of a uniform Bravais lattice. ("Handedness" as used herein refers to the direction convention of the resulting vector as usually mentioned in vector analysis (usually the right-hand rule)). It can be verified by:

[0073] g(q)=det[g αβ (q)]=J 2 (q) (Equation 6)

[0074] where det[g αβ (q)] is the measure [g αβ (q)] tensor. Cotangent basis vector g α (q)(α=1,2,3) is constructed 64 as:

[0075] g 1 (q)≡g2(q)xg3(q) / J(q)

[0076] g 2 (q)≡g3(q)xg1(q) / J(q)

[0077] g 3 (q)≡g1(q)xg2(q) / J(q) (Equation 7)

[0078] The basis tangent and cotangent vectors are orthogonal to each other, where is the Kronecker delta function.

[0079] Similarly, the inverse metric tensor is defined as:

[0080] g αβ (q)≡g α (q)·g β (q),α,β=1,2,3 (Equation 8)

[0081] And the inverse metric tensor is the inverse of the metric tensor, [g αβ (q)]=[g αβ (q)] -1 or

[0082]

[0083] and det[g αβ (q)]=1 / det[g αβ (q)].

[0084] In the case of the basic geometric quantities defined above, the lattice Boltzmann velocity vector is introduced 66 on a general curvilinear grid, analogous to the velocity vector on a standard Cartesian lattice.

[0085]

[0086] constant is a positive or negative integer or zero, and it is an array of three-dimensional coordinates For example, in the so-called D3Q19, the Greek index ranges from 0 to 18,

[0087] c i ∈{(0, 0, 0), (±1, 0, 0), (0±1, 0), (0, 0, ±1), (±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1)}

[0088] A set of moment isotropy and normalization conditions are satisfied in order to recover the correct full Navier-Stokes fluid dynamics. These are that when there exists an appropriate set of constant weights {ω i When i=1,...b}, the set of lattice component vectors allows moment isotropy up to order 6, that is,

[0089]

[0090]

[0091]

[0092]

[0093] where T0 is a constant temperature and depends on the choice of a set of lattice vectors, and is the Kronecker delta function. Note that the three-dimensional array is c i , and e i (q) is the lattice vector of Equation 9.

[0094] A specific set of geometric quantities are defined for use with the LBM model as follows:

[0095]

[0096] When the grid is a uniform Cartesian grid, the term disappear.

[0097] Thus, once the curve mesh is specified, all of the above-mentioned geometric quantities are completely determined, and a curve mesh corresponding to the physical object to be simulated can be stored 68, and because the geometric quantities 0 are completely determined, they can be pre-computed before starting a dynamic LBM simulation (see Figure 2 ).

[0098] refer to Figure 4 , simulation process 46( Figure 2 ) simulates 80 the evolution of the particle distribution according to, for example, a modified lattice Boltzmann equation (LBE) adapted to a curvilinear space. The process performs advection 82 of the particle to the next cell q in the LBM non-Euclidean space, maps 84 the result values ​​obtained from the simulation onto the original curvilinear grid by converting the momentum and spatial coordinates of the coordinate system in the non-Euclidean space to a curvilinear system in the Euclidean space, and renders 86 the mapped result values ​​in the Euclidean space onto a display device, etc.

[0099] refer to Figure 5 , showing more details of the simulation process 46. The simulation process 46 ( Figure 2 ) is simulated according to, for example, a modified lattice Boltzmann equation (LBE) adapted to a curved space 80 ( Figure 4 ) evolution of the particle distribution, see Equation 12. This process performs the advection of particles to the next cell q in the LBM non-Euclidean space 82 ( Figure 4). The process imposes constraints 83 on the particle advection. The constraint process 83 constrains 83a the changes in the particle distribution due to the effective external forces associated with the curvature and inhomogeneities of the general curvilinear grid, imposes 83b constraints on the momentum flux, and uses a defined 83c equilibrium distribution function that recovers the correct Euler and Navier-Stokes equations in curvilinear coordinates in the hydrodynamic limit. Thus, the process maps 84 the resulting values ​​obtained from the simulation by converting the momentum and spatial coordinates in a coordinate system in a non-Euclidean space to a curvilinear system in Euclidean space. Figure 4 ) onto the original curve grid and render the resulting values ​​of the mapping in Euclidean space 86 ( Figure 4 ) to a display device, and so on.

[0100] Reference now Fig. 6A , shows an airfoil shaped body 92 in a general 3D curvilinear mesh 93 in Euclidean space (real world space), where only two dimensions are illustrated (for clarity). As shown, the mesh 93 generally conforms to the shape of the airfoil shaped body 92, and more specifically to the outer surface of the airfoil shaped body 92. The degree of conformity may be greater or less than that shown. Fig. 6A To a lesser extent than shown in .

[0101] Reference now Figure 6B , showing the ( Fig. 6A 92' (illustrated in two dimensions only). In this example, the representation 92' of the airfoil shaped body 92 is illustrated as a rectangular solid body (for illustration purposes only). As shown, the Cartesian grid 95 generally conforms to the shape of the representation 92' of the airfoil shaped body, and specifically to its outer surface. The degree of conformity may be greater or less than that shown.

[0102] Reference now Fig. 7A , Figure 7B , vectors 96 representing the momentum of particles of a fluid are shown in a general 3D curved grid 93 in Euclidean space (real world space, but only two dimensions are shown for clarity). As shown, vectors 96 represent particles of the fluid moving along straight lines through the curved space 96. However, in a variant Cartesian non-Euclidean space, vectors 96' representing the particles are seen to have curved motion.

[0103] Reference now Fig. 8A , Figure 8B ,and Fig. 7A , Figure 7BInstead, vectors 97 representing the momentum of particles of the fluid are shown in a general 3D curved grid 93 in Euclidean space (real world space, but only two dimensions are illustrated for clarity). As shown, vectors 97 represent particles of the fluid moving through curved space 96 in a curved motion. However, in a variant Cartesian non-Euclidean space, vectors 97' representing the particles are treated as having straight-line motion. The specific details of these processes will now be described.

[0104] Volume Lattice Boltzmann Modeling on Curvilinear Lattice

[0105] Described is a volume lattice Boltzmann modeling approach, which, while generally applicable to a variety of formulations, will be discussed as an illustrative example for a so-called "isothermal LBM" formulation.

[0106] The modified distribution is similar to the standard distribution of the evolution of the particle distribution. Volume lattice Boltzmann modeling on the curvilinear lattice is provided as:

[0107] N i (q+c i , t+1)=N i (q, t)+Ω i (q, t)+δN i (q, t) Equation 12

[0108] Where N i (q,t) is the unit element increment for t at time t in unit q belonging to discrete direction c i The number of particles without loss of generality. The term Ω in (Eq. 12) α (q,t) is the collision term that satisfies the local mass and momentum conservation equations:

[0109]

[0110]

[0111] The additional term δN in Equation 12 i (q,t) represents the change in particle distribution due to the effective external forces associated with the curvature and inhomogeneity of the general curvilinear grid. This additional term disappears in the standard LBM on a Cartesian lattice.

[0112] The particle density distribution function is related to N by the following equation i (q,t) related:

[0113] J(q)f i (q, t) = N i (q, t) Equation 14

[0114] As defined above, where J(q) is the volume of the cell centered at q. The basic fluid quantity is given by the standard hydrodynamic moment,

[0115]

[0116]

[0117] where ρ(q,t) and u(q,t) are the fluid density and fluid velocity at position q and time t.

[0118] Using the relationship in (Equation 9), the above velocity moment is rewritten as:

[0119]

[0120] And the velocity component values ​​in the curvilinear coordinate system are given by the following equations:

[0121]

[0122] To simplify the representation, the three-dimensional fluid velocity array is defined as U(q,t)≡(U 1 (q,t),U 2 (q,t)U 3 (q, t)), and Equation 17 is equivalently expressed as:

[0123]

[0124]

[0125] Equation 18 has the same form of fluid velocity as the LBM based on the standard Cartesian lattice. Similarly, the momentum conservation of the collision term in Equation 13 can also be written as:

[0126]

[0127] In LBM, the collision term is usually in a linearized form, where M ij and f i (q, t) denote the collision matrix and equilibrium distribution function, respectively. In particular, the so-called Bhatnagar-Gross-Krook (BGK) form corresponds to:

[0128]

[0129] where τ is the collision relaxation time. In order to recover the correct Navier-Stokes fluid dynamics, the collision matrix needs to satisfy additional conditions in addition to Eqs. 13 and 19.

[0130]

[0131] The BGK form satisfies this additional property. As mentioned above, the additional term δN in Eq. 12 i (q, t) represents the change in particle distribution due to the effective external forces associated with the curvature and inhomogeneity of the general curvilinear grid. This additional term vanishes in the standard LBM on a Cartesian grid. The advection process is an exact one-to-one hop from one site on the curvilinear grid to another site in the standard LBM, i.e.

[0132] N i (q+c i ,t+1)=N′ i (q,t) (Equation 22)

[0133] Where N′ i (q,t) is the post-collision distribution at (q,t). Since it is a curved grid, the number of particles advected from cell q is the same as that arriving at cell q+c. i The number of particles is exactly equal (see Eq. 22), but the momentum changes during advection.

[0134] Generally speaking,

[0135] (e i (q+c i )N i (q+c i ,t+1))≠e i (q)N′ i (q,t)

[0136] In the above equation, the left side of the inequality is the value of momentum at the end of the advection process, while the right side is the value of momentum at the beginning of the process. The inequality exists because the path of the particle is bent (and stretched or compressed) due to the curved grid, so that its velocity at the end of the advection is changed from its original value. This is fundamentally different from the situation on a uniform Cartesian lattice, because the particle has a constant velocity throughout the advection process.

[0137] Therefore, the total momentum value has the following inequality,

[0138]

[0139] where the right side of the inequality sign in (23) represents the total amount of momentum advected from all neighboring cells, while the left side is the total momentum reaching cell q after advection along the curved path. Therefore, from Equations 22 and 23, the net momentum change from all neighboring cells to cell q by advection is given by:

[0140]

[0141] Likewise, the net momentum change to all neighboring cells of cell q by advection out of cell q is given by:

[0142]

[0143] Then, if δN is i (q,t) imposes constraints:

[0144]

[0145]

[0146] Then for the underlying Euclidean space, conservation of mass is preserved and exact conservation of momentum is restored. Here

[0147] χ(q,t)=[χ I (q, t)+χ o (q, t)] / 2

[0148] x(q,t)=[x I (q,t)+x o (q,t)] / 2

[0149] More specifically, the first constraint in Eq. 26 does not pass through δN i (q, t) introduces a mass source. The second constraint in Eq. 26 introduces an “inertial force” that is exactly equal to the amount that satisfies conservation of momentum in the underlying Euclidean space at any lattice site and time t. This mechanism is similar to the description of continuum dynamics in curved space. Written in coordinate component form as F α (q,t)=χ(q,t)·g α (q), which, by directly placing the symbols from the preceding subsection, provides:

[0150]

[0151] The geometric function is defined in (Eq. 11). From (27) and (11), it can be seen that if the curvilinear grid is a regular uniform Cartesian grid, then F α (q, t) vanishes, as in conventional LBM. The second constraint in (26) can also be expressed in coordinate component form as:

[0152]

[0153] To recover the full viscous Navier-Stokes equations, an additional constraint on the momentum flux is imposed as follows,

[0154]

[0155] in

[0156]

[0157] The specific form chosen is:

[0158]

[0159] This equation satisfies the moment constraints of equations 26, 28, and 29. The form of the equilibrium distribution function is defined in order to recover the correct Euler equations and the Navier-Stokes equations in curvilinear coordinates in the hydrodynamic limit. In particular, the following basic conditions for the hydrodynamic moments are:

[0160]

[0161]

[0162]

[0163]

[0164]

[0165] Among them, above,

[0166] These basic conditions are met by the following equilibrium distribution:

[0167]

[0168] The above equilibrium distribution form is similar to the low Mach number expansion of the Maxwell-Boltzmann distribution, but expressed in curvilinear coordinates. In fact, if the curvilinear grid is a grid with g αβ =δ αβ The above equilibrium distribution form reduces to the standard LBM equilibrium distribution if the uniform Cartesian grid of q is used. With all the quantities and dynamics defined above, the lattice Boltzmann equation 12 simulated on the (non-Euclidean) uniform Cartesian lattice {q} obeys Navier-Stokes fluid dynamics in the case of curvilinear coordinates. Therefore, the resulting values ​​can be mapped onto the original curvilinear grid by the following simple transformation,

[0169] ρ(x,t)=ρ(q,t),

[0170] u(x,t)=U α (q,t)g α (q) Equation 34

[0171] Cartesian coordinates in non-Euclidean velocity space

[0172] The variant “Cartesian” lattice {q} resulting from the deformation (bending, twisting, and stretching / compressing) of the original curvilinear grid {x} in Euclidean space can be used to simulate fluid flow near a physical body in the same manner as the conventional Cartesian lattice (x) of LBM, where it is assumed that constraints are applied to particle momentum when advecting in the topology of the non-Euclidean space of the “Cartesian” lattice {q}, as discussed above, to return to a general curvilinear grid in Euclidean space.

[0173] refer to Fig. 9 , the first model (2D-1) 100 is a two-dimensional model including 21 velocities. Among the 21 velocities, one velocity (105) represents a particle that does not move; three sets of four velocities represent particles that move at a normalized velocity (r) (110-113), twice the normalized velocity (2r) (120-123), or three times the normalized velocity (3r) (130-133) in the positive or negative direction along the x or y axis of the lattice; and two sets of four velocities represent particles that move at a normalized velocity (r) (140-143) or twice the normalized velocity (2r) (150-153) relative to both the x and y lattice axes.

[0174] Also like Fig.10 As shown, a second model (3D-1) 200 is shown, which is a three-dimensional model including 39 speeds, each of which is represented by Fig. 9 Of the 39 velocities, one velocity represents a particle that is not moving; three sets of six velocities represent particles that move at a normalized velocity (r), twice the normalized velocity (2r), or three times the normalized velocity (3r) in the positive or negative direction along the x, y, or z axis of the lattice; eight velocities represent particles that move at a normalized velocity (r) relative to all three of the x, y, and z lattice axes; and twelve velocities represent particles that move at twice the normalized velocity (2r) relative to two of the x, y, and z lattice axes.

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

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

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

[0178] The LBM models described above provide a specific class of efficient and robust discrete velocity dynamic models for numerical simulation of flows in two and three dimensions. This type of model includes a specific set of discrete velocities and weights associated with those velocities. These velocities correspond to grid points of Cartesian coordinates in velocity space (in non-Euclidean space), which facilitates accurate and efficient implementation of discrete velocity models (especially the class known as lattice Boltzmann models). Using such models, flows can be simulated with high fidelity.

[0179] refer to Fig.11, a physical process simulation system operates according to process 300 to simulate a physical process such as fluid flow. Prior to simulation, the simulation space is modeled as a collection of voxels (step 302). Typically, the simulation space is generated using a computer-aided design (CAD) program. For example, a CAD program can be used to draw a microdevice located in a wind tunnel. Thereafter, the data generated by the CAD program is processed to add a lattice structure with appropriate resolution and to account for objects and surfaces within the simulation space.

[0180] The resolution of the lattice can be chosen based on the Reynolds number of the system being simulated. The Reynolds number is related to the viscosity of the flow (v), the characteristic length of objects in the flow (L), and the characteristic velocity of the flow (u):

[0181] Re=uL / v Equation (I.5)

[0182] The characteristic length of an object represents a large-scale feature of the object. For example, if the flow around a microdevice is being simulated, the height of the microdevice can be considered the characteristic length. When the flow around a small area of ​​an object (e.g., the side-view mirror of a car) is of interest, the resolution of the simulation can be increased, or a region of increased resolution can be employed around the area of ​​interest. The dimensionality of the voxel decreases as the resolution of the lattice increases.

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

[0184] The number of states is determined by the number of possible velocity vectors within each energy level. The velocity vectors consist of integer linear rates in a space with three dimensions x, y, and z. For multi-genus simulations, the number of states increases.

[0185] Each state i represents a different velocity vector at a particular energy level (i.e., energy level zero, one, or two). The velocity c of each state i The "velocity" in each of the three dimensions is indicated as follows:

[0186] c i =(c ix ,c iy ,c iz ,). Equation (I.6)

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

[0188] Generating all possible permutations of the three energy levels gives 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).

[0189] Each voxel (i.e., each lattice site) is represented by a state vector f(x). This state vector completely defines the state of the voxel and includes 39 entries. These 39 entries correspond to one energy zero state, 6 energy one states, 8 energy three states, 6 energy four states, 12 energy eight states, and 6 energy nine states. By using this set of velocities, the system can generate Maxwell-Boltzmann statistics for the achieved equilibrium state vector.

[0190] For processing efficiency, voxels are grouped in 2x2x2 volumes called microblocks. Microblocks are organized to allow parallel processing of voxels and minimize the overhead associated with data structures. The shorthand notation for voxels in a microblock is defined as Ni(n), where n represents the relative position of the lattice site in the microblock and n{0,1,2,...,7}. Microblocks are Fig.12 Shown in.

[0191] refer to Fig.13A and 13B , surface S( Fig.13A ) in the simulation space ( Fig. 13B ) is represented as a facet F α A collection of:

[0192] S={F α Equation (I.7)

[0193] Where α is an index that enumerates a particular facet. Facets are not restricted to the boundaries of voxels, but typically have a size that is on the order of the size of the voxels adjacent to the facet, or has a size that is slightly smaller than the size of the voxels adjacent to the facet, so that the facet affects a relatively small number of voxels. To implement surface dynamics, attributes are assigned to facets. Specifically, each facet F α With unit normal (n α ), surface area (A α )、Center position(x α ) and the surface element distribution function (fi (α)).

[0194] refer to Fig.14 , different resolution levels can be used in different regions of the simulation space to improve processing efficiency. Typically, the region 650 surrounding the object 655 is of greatest concern and is therefore simulated using the highest resolution. Because the effect of viscosity decreases with distance from the object, regions 660, 665 spaced at increasing distances from the object 655 are simulated using reduced resolution levels (i.e., enlarged voxel volumes). Similarly, as Fig.15 As shown in , lower resolution levels can be used to simulate regions 770 around less prominent features of object 775, while the highest resolution level is used to simulate regions 780 around the most prominent features (e.g., leading and trailing edge surfaces) of object 775. Outlying regions 785 are simulated using the lowest resolution level and the largest voxels.

[0195] C. Identify the voxels affected by the surface element

[0196] Reference again Fig.11 Once the simulation space has been modeled (step 302), the voxels affected by one or more facets are identified (step 304). Voxels can be affected by facets in a number of ways. First, a voxel intersected by one or more facets is affected in that the voxel has a reduced volume relative to non-intersecting voxels. This occurs because the facet and the material beneath the surface represented by the facet occupy a portion of the voxel. The fractional factor P f (x) indicates the portion of the voxel that is not affected by the facet (i.e., the portion that can be occupied by fluid or other material for which flow is simulated). For non-intersecting voxels, P f (x) is equal to 1.

[0197] Voxels that intersect one or more panels by transmitting particles to or receiving particles from the panels are also identified as voxels affected by the panels. All voxels intersected by a panel will include at least one state for receiving particles from the panel and at least one state for transmitting particles to the panel. In most cases, additional voxels will also include such states.

[0198] refer to Fig.16 , for a non-zero velocity vector c i For each state i, the face F α From the parallelepiped G iα The area defined by G receives particles or transmits particles to it, where the parallelepiped G iα With the velocity vector c i and facets (|c i n i |)'s unit normal n αThe height is defined by the magnitude of the vector dot product and the surface area A of the bin α Define the base so that the parallelepiped G iα The volume V iα equal:

[0199] V ia =|c i n a |A a Equation (I.8)

[0200] When the velocity vector of the state points to the surface element (|c i n i |<0), facet F α From the volume V iα Receive particles, and when the velocity vector of the state points away from the surface element (|c i n i |>0), and transmit particles to this region. As will be discussed below, when another facet occupies the parallelepiped G iα , i.e., a condition that occurs near non-convex features such as interior corners, this expression must be modified.

[0201] Panel F α The parallelepiped G iα Part or all of a plurality of voxels may be overlapped. The number of voxels or portions thereof depends on the size of the bin relative to the size of the voxel, 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 of the order of magnitude of or smaller than the size of a voxel located near the bin.

[0202] Parallelepiped G iα The portion of the voxel N(x) that overlaps is defined as V iα (x). Using this term, in voxel V iα (x) and panel F α The flux of particles moving between states i is Γ iα (x) is equal to the density of particles in state i in the voxel (N i (x)) multiplied by the volume of the region overlapping the voxel (V iα (x)):

[0203] Γ ia (x) = N i (x)+V ia (x) Equation (I.9)

[0204] When the parallelepiped G iα When intersected by one or more polygons, the following conditions are true:

[0205] V ia =∑V a (x)+∑V ia (β) Equation (I.10)

[0206] The first summation indicates that G iα All voxels overlap and the second term describes the iα All the face elements intersecting. When the parallelepiped G iα When not intersected by another face, this expression simplifies to:

[0207] V ia =∑V ia (x) Equation (1.11)

[0208] D. Execute simulation

[0209] Once the voxels affected by one or more facets are identified (step 304), a timer is initialized to begin the simulation (step 306). During each time increment of the simulation, the movement of particles from voxel to voxel is simulated by an advection phase that accounts for the interaction of particles with surface facets (steps 308-316). Next, a collision phase (step 318) simulates the interaction of 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), the advection phase and the collision phase (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).

[0210] 1. Boundary conditions for surfaces

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

[0212] The other two boundary conditions are related to the net momentum of the particle interacting with the surface element. For a surface without skin friction (referred to herein as a sliding surface), the net tangential momentum flux must be equal to zero and the net normal momentum flux must be equal to the local pressure at the surface element. Therefore, the normal n to the surface element is α The components of the combined received and transmitted momentum that are perpendicular (i.e., the tangential components) must be equal and have a direction perpendicular to the surface element normal n.α The difference between the parallel components of the combined received and transmitted momentum (i.e., the normal components) must equal the local pressure at that cell. For a non-slip surface, the friction of the surface reduces the combined tangential momentum of the particles transmitted by the cell relative to the combined tangential momentum of the particles received by the cell by a factor related to the amount of friction.

[0213] 2. Aggregation from voxels to surfaces

[0214] As a first step in simulating the interaction between particles and surfaces, particles are aggregated from volume elements and provided to surface elements (step 308). As noted above, the volume element N(x) and the surface element F α The flux of particles in state i is:

[0215] Γ iα (x) = N i (x)V iα (x) Equation (1.12)

[0216] From this it can be seen that for the pointing surface element F α Each state i(c i n α <0), provided by the voxel to the surface element F α The number of particles is:

[0217] Γ iαV→F =∑ X Γ iα (x) = ∑ X N i (x)V iα (x) Equation (I.13)

[0218] Only its V iα (x) The voxels with non-zero values ​​must be summed. As noted above, the size of the bins is chosen so that V iα (x) has non-zero values ​​only for a few voxels. iα (x) and P f (x) can have non-integer values, so Γ α (x) is stored and processed as a real number.

[0219] 3. Movement from panel to panel

[0220] Next, the particles are moved between the bins (step 310). α The incoming state (c i n α <0) parallelepiped G iα By the other side element F β If they intersect, then F αThe part of the particles in state i received will come from the surface element F β In particular, the face element F α will receive the bin F during the previous time increment β The fraction of particles in state i that are generated. This relationship is Fig.17 As shown in the figure, the surface element F β Intersecting parallelepipeds G iα The part 1000 is equal to the surface element F α Intersecting parallelepipeds G iβ As noted above, the intersecting portion is denoted as V iα (β). Using this term, the panel F β With face element F α The flux of particles between states i can be described as:

[0221] Γ iα (β, t-1) = Γ i (β)V iα (β) / V iα Equation (I.14)

[0222] where Γ i (β,t-1) is the value of the bin F during the previous time increment. β The measurement of the state i particle produced. From this, it can be seen that for the pointing surface element F α Each state i(c i n α <0), provided by other bins to bin F α The number of particles is:

[0223] Γ iαF→F =∑ β G iα (β)=∑ β Γ i (β, t-1)V iα (β) / V iα Equation (I.15)

[0224] And the total flux of particles of state i into the bin is:

[0225] Γ iIN (α) = Γ iαF→F +Γ iαF→F =∑ x N i (x)V iα +∑ β Γ i (β, t-1)V iα (β) / V iα

[0226] Equation (I.16)

[0227] The state vector N(α) for the bin (also called the bin distribution function) has M entries corresponding to the M entries of the voxel state vector. M is the number of discrete lattice rates. The input states to the bin distribution function N(α) are set equal to the flux of particles into those states divided by the volume V iα :

[0228] For c i n α <0,N i (α) = Γ iIN (α) / V iα Equation (I.17)

[0229] The bin distribution function is a simulation tool used to generate output flux from bins, and does not necessarily represent actual particles. In order to generate an accurate output flux, values ​​are assigned to the other states of the distribution function. The outward states are filled using the technique described above for filling the inward states:

[0230] For c i n α ≥0,N i (α) = Γ iOTHER (α) / V Equation (I.18)

[0231] where Γ iOTHER (α) is the method used to generate Γ iIN (α), but applies the technique to all states except the incoming state (c i n α <0) (c i n α ≥0). In an alternative approach, Γ iOTHER (α) can use Γ from the previous time step iOTHER The value of (α) is generated so that:

[0232] Γ iOTHER (α,t)=Γ iOUT (α,t-1) Equation (I.19)

[0233] For the parallel state (c i n α =0), V iα and V iα (x) are all zero. i In the expression of (α), V iα (x) appears in the numerator (according to the iOTHER The expression of (α) and V iα appears in the denominator (according to thei (α) expression). Therefore, when V iα and V iα When (x) is close to zero, N for parallel state i (α) is determined as N i 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 initialized at the beginning of the simulation based on the initial conditions for temperature and pressure. These values ​​are then adjusted over time.

[0234] 4. Performing panel surface dynamics

[0235] Next, surface dynamics are performed on each surface element to satisfy the four boundary conditions discussed above (step 312). The procedure for performing surface dynamics for a surface element is described in Fig.18 Initially, the combined momentum P(α) of the particle at the surface element is determined to determine the surface element F α Combined momentum (step 1105):

[0236] For all i,

[0237] According to this equation, the normal momentum P n (α) is determined as:

[0238] P n (α) = n α ·P(α) Equation (I.21)

[0239] This normal momentum is then removed using a push / pull technique (step 1110) to produce N n- (α). According to this technique, particles are moved between states in a manner that affects only the normal momentum. The push / pull technique is described in U.S. Patent No. 5,594,671, which is incorporated herein by reference.

[0240] Afterwards, N n- (α) particles are collided to produce the Boltzmann distribution N n-β (α) (step 1115). As described below with respect to implementing fluid dynamics, the Boltzmann distribution can be obtained by adding N n- (α) This is achieved by applying a set of collision rules.

[0241] Then, based on the incoming flux distribution and the Boltzmann distribution, the α The outgoing flux distribution of (step 1120). First, the incoming flux distribution Γ i The difference between (α) and the Boltzmann distribution is determined as:

[0242] ΔΓi (α) = Γ iIN (α)-N n-βi (α)V iα Equation (I.22)

[0243] Using this difference, the outgoing flux distribution is:

[0244] For n α c i >0,Γ iOUT (α) = N n-βi (α)V iα -.Δ.Γ i* (α) Equation (I.23)

[0245] and where i* is the state with the opposite direction to state i. For example, if state i is (1,1,0,0), then state i* is (-1,-1,0,0). To account for skin friction and other factors, the outgoing flux distribution can be further refined as:

[0246] For n α c i >0,

[0247] Among them C f is a function of skin friction, t iα is with n α The vertical first tangent vector, t 2α is with n α and t 1α are perpendicular to the second tangent vector, and ΔN j,1 and ΔN j,2 is the distribution function of the energy (j) corresponding to state i and the indicated tangent vector. The distribution function is determined according to the following formula:

[0248]

[0249] where j is equal to 1 for the level 1 state and equal to 2 for the level 2 state.

[0250] For Γ iOUT The function of each term in equation (α) is as follows. The first and second terms enforce the normal momentum flux boundary conditions to the extent that the collisions have effectively produced a Boltzmann distribution but include anomalies in the tangential momentum flux. The fourth and fifth terms correct for such anomalies, which may arise from discrete effects or non-Boltzmann structures due to insufficient collisions. Finally, the third term adds a specified amount of skin friction to implement the desired change in the tangential momentum flux at the surface. The friction coefficient C f The generation of is described below. Note that all terms involved in vector manipulation are geometric factors that can be calculated before starting the simulation.

[0251] From this, the tangential velocity is determined as:

[0252] u i (α)=(P(α)-P n (α)n α ) / ρ, Equation (I.26)

[0253] where ρ is the density of the bin distribution:

[0254]

[0255] As before, the difference between the incoming flux distribution and the Boltzmann distribution is determined as:

[0256] ΔΓ i (α) = Γ iIN (α)-N n-βi (α)V iα Equation (I.28)

[0257] Then, the outgoing flux distribution becomes:

[0258] Γ iOUT (α) = N n-βi (α)V iα -ΔГ i* (α)+C f (n α c i )[N n-βi *(α)-N n-βi (α)]V iα Equation (I.29)

[0259] It corresponds to the first two rows of the outgoing flux distribution determined by the previous technique, but without the need for correction for the anomalous tangential flux.

[0260] Using either method, the resulting flux distribution satisfies all the momentum flux conditions, namely:

[0261]

[0262] Among them, p α It is in the face element F α The equilibrium pressure at the bin and based on the average density and temperature values ​​of the volume element providing the particle to the bin, and u α is the average velocity at that bin.

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

[0264]

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

[0266] For c ji n α >0,

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

[0268] For c ji n α >0,Γ αjiOUTf =Γ αjiOUT +δΓ αji Equation (I.33)

[0269] This operation corrects the mass and energy fluxes while leaving the tangential momentum flux unchanged. If the flow in the neighborhood of the panel is roughly uniform and close to equilibrium, this adjustment is small. After the adjustment, the resulting normal momentum flux is slightly altered to a value that is the equilibrium pressure based on the neighborhood average properties plus a correction due to non-uniform or non-equilibrium properties of the neighborhood.

[0270] 5. Movement from voxel to voxel

[0271] Reference again Figure 3 , particles are moved between voxels along the three-dimensional rectilinear lattice (step 314). This voxel-to-voxel movement is the only movement operation performed on voxels that do not interact with surface elements (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 make up the majority of the voxels.

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

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

[0274] For voxels that interact with one or more surfaces, the move operation becomes slightly more complex. This results in one or more fractional particles being teleported to the surface. This fractional particle to surface teleportation results in fractional particles remaining in the voxel. These fractional particles are teleported to the voxels occupied by the surface.

[0275] refer to Fig.19 , when a portion 900 of the state i-particle of voxel 905 is moved to bin 910 (step 308), the remaining portion 915 is moved to voxel 920 in which bin 910 is located and from which the state i-particle points to bin 910. Thus, if the state population equals 25 and V iα (x) is equal to 0.25 (i.e., one quarter of the volume element and the parallelepiped G iα intersect), then 6.25 particles will be moved to the surface element F α And 18.75 particles will be moved to the surface element F α Since multiple facets can intersect a single facet, the number of particles in state i that are teleported to a facet N(f) occupied by one or more facets is:

[0276]

[0277] Where N(x) is the source voxel.

[0278] 6. Dispersion from surface elements to volume elements

[0279] Next, the outgoing particles from each bin are scattered to the voxels (step 316). Essentially, this step is the inverse of the aggregation step of moving particles from voxels to bins. α The number of particles in state i that move from the surface element to the volume element N(x) is:

[0280]

[0281] Where Pf (x) indicates that the volume of some voxels is reduced. Therefore, for each state i, the direction from the face element to the voxel N (x) The total number of particles is:

[0282]

[0283] After dispersing particles from bins to voxels, combining the particles with those that have been advected from surrounding voxels, and integerizing the results, it is possible that some directions in some voxels may underflow (become negative) or overflow (exceed 255 in an eight-bit implementation). This will result in a gain or loss of mass, momentum, and energy after these quantities have been truncated to fit within the allowed range of values. To prevent this from happening, the out-of-bounds mass, momentum, and energy are accumulated before the truncation of the offending state. For the energy to which the state belongs, an amount of mass equal to the value of the gain (due to underflow) or loss (due to overflow) is added back to a randomly (or sequentially) selected state with the same energy and which is not itself subject to overflow or underflow. The additional momentum resulting from this addition of mass and energy is accumulated and added to the momentum from the truncation. By only adding mass 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.

[0284] 7. Perform fluid dynamics

[0285] Perform fluid dynamics (step 318) Fig.11 This step can be called microdynamics or intra-voxel operations. Similarly, the advection procedure can be called inter-voxel operations. The microdynamics operation described below can also be used to collide particles in bins to produce Boltzmann distributions.

[0286] The fluid dynamics are ensured in the lattice Boltzmann equation model by a special collision operator known as the BGK collision model. This collision model mimics the dynamics of the distribution in a real fluid system. The collision process can be well described by equations I.1 and the right-hand side of equation I.2. After the advection step, the conserved quantities of the fluid system, specifically density, momentum and energy, are obtained from the distribution function using equation I.3. From these quantities, the f in equation (I.2) eq The equilibrium distribution function indicated is completely governed by equation (I.4). The velocity vector set c i The choices of , weights are listed in Table 1, together with Equation I.2 to ensure that the macroscopic behavior follows the correct fluid dynamics equations.

[0287] E. Variable resolution

[0288] refer to Fig. 20, variable resolution (as discussed above) employs voxels of different sizes, hereinafter referred to as coarse voxels 12000 and fine voxels 1205. (The following discussion relates to voxels having two different sizes; it will be appreciated that the described techniques may be applied to voxels of three or more different sizes to provide higher levels of resolution). The interface between the regions of coarse and fine voxels is referred to as a variable resolution (VR) interface 1210.

[0289] When variable resolution is used at or near a surface, the facets can interact with the voxels on both sides of the VR interface. These facets are classified as VR interface facets 1215 (F αIC ) or VR fine element 1220 (F αIF ). VR interface facet 1215 is a facet located on the rough side of the VR interface and having a rough parallelepiped 1225 extending into the fine voxel. (The rough parallelepiped is determined based on the dimensions of the rough voxel c i The dimension of the fine parallelepiped is determined by the dimension of the fine volume element c i The VR fine facet 1220 is a facet located on the fine side of the VR interface and has a fine parallelepiped 1230 extending into the coarse facet. Processing related to the interface facet may also involve processing related to the coarse facet 1235 (F αC ) and fine bin 1240 (F αF ) interactions.

[0290] For both types of VR surfaces, surface dynamics are performed at a fine scale and operate as described above. However, VR surfaces differ from other surfaces in the way particles are advected to and from the VR surface.

[0291] Fig.21 A screen shot of a fluid simulation is shown. A fluid simulation using the curve method described above (rather than the conventional Cartesian method or possible other methods for a curved mesh) will provide a similar description of the fluid simulation and any customized corresponding computational data output. However, when an object having a curved surface (e.g., an actual physical object) is modeled, such a fluid simulation using the curve method described above can perform the fluid simulation faster and / or by using fewer computational resources than other methods.

[0292] Many implementations have been described. However, it will be appreciated that various modifications may be made without departing from the spirit and scope of the claims. Therefore, other implementations are within the scope of the appended claims.

Claims

1. A computer-implemented method for simulating fluid flow near a surface of a solid, the method comprising: receiving, by a computing system, a Cartesian coordinate system in a non-Euclidean space representing a curvilinear mesh that conforms to a surface of the solid; simulating the transport of particles in a volume of a fluid using a set of lattice velocities according to a modified lattice Boltzmann equation adapted to represent a Cartesian coordinate system in a non-Euclidean space of a curvilinear grid, wherein the transport causes collisions between the particles; executing a distribution function for transport of the particles, wherein the distribution function includes particle collision determinations and changes in particle distribution associated with the curvilinear grid; performing, by a computing system, an advection operation in the coordinate system under constraints applied to the particle momentum values; and Values ​​obtained from the simulation are mapped onto the curvilinear grid by a computing system by converting particle momentum values ​​and spatial coordinates determined in the Cartesian coordinate system in non-Euclidean space to momentum and spatial values ​​in the curvilinear coordinate system of the curvilinear grid.

2. The computer-implemented method of claim 1, wherein: Solid surfaces are curved solid surfaces.

3. The computer-implemented method of claim 1, wherein: The grid is a curvilinear grid in Euclidean space.

4. The computer-implemented method of claim 1, wherein: The coordinate system is a non-Euclidean space defined in Riemannian geometry.

5. The computer-implemented method of claim 1, wherein: A general curvilinear coordinate system has a one-to-one mapping to a coordinate system of choice, which is a Cartesian coordinate system in a non-Euclidean space.

6. The computer-implemented method of claim 3, wherein: Particle advection is performed such that at each time step, the particle distribution function is uniquely advected from one voxel to another.

7. The computer-implemented method of claim 1, wherein: The change in particle distribution causes the computing system to simulate the application of forces applied in each voxel at each time increment in the curvilinear grid to enforce exact momentum conservation in the coordinate system.

8. The computer-implemented method of claim 3, wherein: A volume representation is used so that each voxel in a general curvilinear grid has a defined volume, and mass conservation is enforced, and mass and momentum densities are defined.

9. The computer-implemented method of claim 1, wherein: The mapping also includes: Determine momentum states that are constrained due to the representation of meshes and solid surfaces in Cartesian coordinate systems in non-Euclidean spaces.

10. An apparatus for simulating fluid flow near a surface of a solid, the apparatus comprising: Memory; One or more processor devices configured to: receiving a Cartesian coordinate system in a non-Euclidean space for representing a curvilinear mesh that conforms to a surface of the solid; simulating the transport of particles in a volume of a fluid using a set of lattice velocities according to a modified lattice Boltzmann equation adapted to represent a Cartesian coordinate system in a non-Euclidean space of a curvilinear grid, wherein the transport causes collisions between the particles; executing a distribution function for transport of the particles, wherein the distribution function includes particle collision determinations and changes in particle distribution associated with the curvilinear grid; performing an advection operation in said coordinate system, subject to constraints applied to the momentum values ​​of the particles; and The values ​​obtained from the simulation are mapped onto the curvilinear grid by converting the particle momentum values ​​and spatial coordinates determined in the Cartesian coordinate system in non-Euclidean space to momentum and spatial values ​​in the curvilinear coordinate system of the curvilinear grid.

11. The device of claim 10, wherein: Solid surfaces are curved solid surfaces.

12. The device of claim 10, wherein: The grid is a curvilinear grid in Euclidean space.

13. The device of claim 10, wherein: The coordinate system is a non-Euclidean space defined in Riemannian geometry.

14. The device of claim 10, wherein: A general curvilinear coordinate system has a one-to-one mapping to a coordinate system of choice, which is a Cartesian coordinate system in a non-Euclidean space.

15. The device of claim 12, wherein: Particle advection is performed such that at each time step, the particle distribution function is uniquely advected from one voxel to another.

16. The device of claim 10, wherein: The change in particle distribution causes the one or more processor devices to simulate the application of forces on each voxel at each time increment in the curvilinear grid to enforce exact momentum conservation in the coordinate system.

17. The device of claim 12, wherein: A volume representation is used so that each voxel in a general curvilinear grid has a defined volume, and mass conservation is enforced, and mass and momentum densities are defined.

18. The device of claim 10, wherein: The mapping also includes: Determine momentum states that are constrained due to the representation of meshes and solid surfaces in Cartesian coordinate systems in non-Euclidean spaces.

19. One or more machine-readable hardware storage devices storing instructions executable by one or more processing devices to cause a computer to: receiving a Cartesian coordinate system in a non-Euclidean space for representing a curvilinear mesh, the curvilinear mesh conforming to a surface of a solid; simulating the transport of particles in a volume of a fluid using a set of lattice velocities according to a modified lattice Boltzmann equation adapted to represent a Cartesian coordinate system in a non-Euclidean space of a curvilinear grid, wherein the transport causes collisions between the particles; executing a distribution function for transport of the particles, wherein the distribution function includes particle collision determinations and changes in particle distribution associated with the curvilinear grid; performing an advection operation in said coordinate system, subject to constraints applied to the momentum values ​​of the particles; and The values ​​obtained from the simulation are mapped onto the curvilinear grid by converting the particle momentum values ​​and spatial coordinates determined in the Cartesian coordinate system in non-Euclidean space to momentum and spatial values ​​in the curvilinear coordinate system of the curvilinear grid.

20. The one or more machine-readable hardware storage devices of claim 19, wherein: The grid is a curvilinear grid in Euclidean space.

21. The one or more machine-readable hardware storage devices of claim 19, wherein: Solid surfaces are curved solid surfaces.

22. The one or more machine-readable hardware storage devices of claim 19, wherein: The coordinate system is a non-Euclidean space defined in Riemannian geometry.

23. The one or more machine-readable hardware storage devices of claim 19, wherein: A general curvilinear coordinate system has a one-to-one mapping to a coordinate system of choice, which is a Cartesian coordinate system in a non-Euclidean space.

24. The one or more machine-readable hardware storage devices of claim 20, wherein: Particle advection is performed such that at each time step, the particle distribution function is uniquely advected from one voxel to another.

25. The one or more machine-readable hardware storage devices of claim 19, wherein: The change in particle distribution causes the computer to simulate the application of forces on each voxel at each time increment in the curvilinear grid to enforce exact momentum conservation in the coordinate system.

26. The one or more machine-readable hardware storage devices of claim 20, wherein: A volume representation is used so that each voxel in a general curvilinear grid has a defined volume, and mass conservation is enforced, and mass and momentum densities are defined.

27. The one or more machine-readable hardware storage devices of claim 19, wherein: The mapping also includes: Determine momentum states that are constrained due to the representation of meshes and solid surfaces in Cartesian coordinate systems in non-Euclidean spaces.

Images