Stabilizing explicit numerical schemes

Abstract

The present disclosure relates to stabilizing explicit numerical schemes. A computer-implemented technique for simulating fluid flow around a solid surface is disclosed, comprising: receiving a model of a simulation space comprising a lattice structure represented as a collection of voxels and a representation of a physical object, wherein the voxels have an appropriate resolution for accounting for the surface of the physical object; simulating the movement of particles in the fluid volume, wherein the movement of the particles causes collisions between the particles; identifying, by a computing system, a face between two voxels, wherein at least one face violates a stability condition; computing, by the computing system, a modified heat flux using a spatially averaged temperature gradient near the two voxels, wherein at least one face violates the stability condition; and performing, by the computing system, an advection operation on the particles to a subsequent voxel.

Description

Background Art

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

[0002] Computational fluid dynamics is a branch of fluid mechanics that involves computer-implemented numerical analysis techniques to analyze and simulate fluid flows in physical environments.

[0003] The so-called "Lattice Boltzmann Method" (LBM) is a powerful 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 fundamental dynamical processes in a basic Boltzmann dynamics 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 based on their microscopic velocities.

[0004] A common aspect found in real computational fluid dynamics simulation problems involving the "lattice Boltzmann method" (LBM) are problems involving complex geometries, such as irregular grids for surface and volume discretization. Summary of the Invention

[0005] Numerical simulations of diffusion-dominated physical phenomena are commonly used for applications involving conductive heat transfer, mass diffusion, electrical conduction, and the like. The governing equations for these phenomena are formulated as a set of partial differential equations (PDEs) that include unsteady diffusion and volumetric source terms. The numerical solution involves discretizing the spatial domain of interest and applying time integration techniques to advance the solution in time. The spatial discretization process is typically implemented using highly automated grid generation tools, while the temporal discretization (for applying time integration techniques to advance the solution in time, i.e., the time-step) is chosen to ensure the stability and accuracy of the numerical solution at an acceptable computational cost.

[0006] In particular, the stability properties of the Courant-Friedrichs-Lewy (CFL) constraint, often referred to as a time marching scheme, determine the maximum time step that can be used without introducing significant instabilities to the solution. Two types of time marching schemes are commonly employed—implicit and explicit.

[0007] Implicit methods satisfy the CFL constraints by construction and can therefore use large time steps without making the solution unstable (however, too large a time step generally leads to inaccurate results). Implicit methods require the solution of a large system of matrix coefficients, making their practical implementation both non-trivial and computationally expensive.

[0008] On the other hand, explicit methods are very simple to implement, computationally cheap (per iteration), and highly parallelizable. However, explicit methods need to satisfy a CFL constraint. This CFL constraint for explicit diffusion schemes stipulates that the CFL number given by (κΔ_t) / (Δ_x^2) is less than a certain limit (which is O(1)), where κ is the diffusion rate, Δ_x is the size of the minimum spatial grid, and Δ_t is the time step. In other words, if the spatial grid size Δ_x decreases by a factor F anywhere in the domain, then the time step Δ_t will have to increase by a factor F. 2 Reduce to maintain numerical stability.

[0009] Therefore, for spatial grids with small-sized elements, explicit methods may require extremely small time steps, severely impacting simulation performance. This is true even if the number of such small-sized elements is very limited in the simulation domain, since the smallest element in the entire domain determines the CFL condition and therefore the time step.

[0010] For practical problems involving complex geometries, the use of irregular grids for surface and volume discretization is unavoidable. On these grids, Δ_x may vary significantly, and the use of explicit schemes may become very inefficient due to the extremely small time steps required by the CFL constraint.

[0011] Therefore, explicit solution implementers spend a lot of time and effort trying to improve the quality of the spatial grid in an attempt to alleviate this problem. Even so, it is almost impossible to remove all small-sized elements from any discretization of the actual spatial domain, and therefore, small time steps are (at least locally) the only way to make the explicit solution stable.

[0012] According to one aspect, a system, computer-implemented method, and computer program product for simulating fluid flow around a solid surface include receiving a model of a simulation space comprising a lattice structure and a representation of a physical object represented as a collection of voxels, wherein the voxels have an appropriate resolution to account for the surface of the physical object, simulating the motion of particles in the fluid volume, wherein the motion of the particles causes collisions between the particles, identifying faces (facets) between two voxels, wherein at least one of the faces violates a stability condition, calculating a modified heat flux using a spatially averaged temperature gradient near the two voxels, wherein the at least one of the faces violates the stability condition, and performing an advection operation on the particles to reach a subsequent voxel.

[0013] The disclosed technique introduces modifications to the heat flux calculation between two adjacent elements when at least one element violates a constraint (e.g., a Counton-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 quantities of interest in the immediate vicinity of the elements, and help stabilize the numerical solution regardless of the size of the two elements and ensure spatio-temporal accuracy. When the two adjacent elements are large (and therefore satisfy the CFL constraint), the newly proposed flux calculation is implemented as an explicit solution, implicitly indicating that this new approach is consistent with the explicit method and also overcomes the aforementioned drawbacks of the explicit method.

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

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

[0016] Figure 1 Depicted is a system for simulating fluid flow.

[0017] Figure 2 A flow chart showing an operational simulation for heat flux calculations based on a lattice Boltzmann model is depicted.

[0018] Figure 3-Figure 7 A flow chart is shown involving aspects of modifying a heat flux calculation between two adjacent elements when at least one of the elements violates a constraint.

[0019] Figure 8 and 9Velocity components of two LBM models represented in Euclidean space are shown (prior art).

[0020] Figure 10 is a flow chart of the process followed by a physical process simulation system using a modified heat flux calculation between two adjacent elements.

[0021] Figure 11 is a perspective view of a microscopic block (prior art).

[0022] Figures 12A-12B yes Figure 1 Illustration of the lattice structure used in the system (prior art).

[0023] Figure 13 and 14 Variable resolution technology (prior art) is shown.

[0024] Figure 15 The motion of particles is shown (prior art).

[0025] Figure 16 The area affected by the bins of the surface is shown (prior art).

[0026] Figure 17 Surface dynamics (prior art) are shown.

[0027] Figure 18 is a flow chart of the process used to perform surface dynamics. DETAILED DESCRIPTION

[0028] Modeling the simulation space

[0029] 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),

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

[0031] This equation describes the distribution function f iThe 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 occurs when a fluid block starts at a mesh location and moves along one of the velocity vectors to the next mesh location. At that point, the "collision operator" is calculated, that is, the influence of the nearby fluid block on the starting fluid block. The fluid can only move to one mesh location, so the appropriate choice of velocity vector is necessary so that all components of all velocities are multiples of a common speed.

[0032] 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 a prescribed value given by the second equation, which is the "equilibrium" form,

[0033]

[0034] 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)} is constructed from the physical point of view of a well-defined local equilibrium given by:

[0035]

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

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

[0038] 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 equation of the following form:

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

[0040] 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 iThe (x,t)-derived hydrodynamic moments obey the Navier-Stokes equations in the macroscopic limit. These moments are defined as:

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

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

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

[0044] The standard technique for deriving the macroscopic equation 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 employed. In a fluid system, small perturbations in density travel at the speed of sound. In gaseous systems, the speed of sound is generally determined by the temperature. The significance of the effects of compressibility in a fluid is measured by the ratio of the characteristic velocity to the speed of sound, which is known as the Mach number.

[0045] 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.

[0046] refer to Figure 1 , shows a system 10 for simulating, for example, fluid flow about a representation of a physical object. In this implementation, the system 10 is based on a client-server architecture and includes a server system 12 implemented as a massively parallel computing system 12 and a client system 14, the server system 12 including a memory 18, a bus system 11, interfaces 20 (e.g., user interface / network interface / display or monitor interface, etc.), and a processing device 24, which together with a grid and a simulation engine 34 provide a simulation process 30.

[0047] In memory 18 are a grid preparation engine 32 and a simulation engine 34. Although Figure 1Grid preparation engine 32 is shown in memory 18, but may be a third-party application executing on a different system (e.g., system 14 or another system) than server 12. Regardless of whether grid preparation engine 32 executes in memory 18 or on a different system than server 12, grid preparation engine 32 receives user-supplied grid definition 28, and grid preparation engine 32 prepares the grid and sends the prepared grid to simulation engine 34.

[0048] The simulation engine 34 includes a particle collision interaction module 34a, a particle boundary model module 34b, and an advection module 34c that performs advection operations. The system 10 accesses a data repository 38 that stores 2D and / or 3D grids and libraries. The advection module 34c includes a submodule 36 that performs advection operations based on a correction to the flux calculation between two adjacent elements, as discussed below.

[0049] Before executing a simulation in a simulation engine, the simulation space is modeled as a collection of voxels. Typically, a computer-aided design (CAD) program is used to generate the simulation space. For example, a CAD program can be used to draw a microdevice positioned in a wind tunnel. The data generated by the CAD program is then processed to add a lattice structure with appropriate resolution and to account for objects and surfaces within the simulation space.

[0050] Now refer to Figure 2 , a process for simulating fluid flow around a representation of a physical object is shown. In the examples to be discussed herein, the physical object is an airfoil. However, the use of an airfoil is merely illustrative, as the physical object may be of any shape and may in particular have (one or more) planar and / or curved surfaces. The process receives 35a a mesh for the physical object being simulated, for example from a client system 14 or by retrieval from a data repository. In other embodiments, an external system or server 12 generates a mesh for the physical object being simulated based on user input. The process precomputes 35b geometric quantities based on the retrieved mesh and performs 35c a dynamic lattice Boltzmann model simulation using the precomputed geometric quantities corresponding to the retrieved mesh. The lattice Boltzmann model simulation comprises performing a lattice Boltzmann model simulation based on the mesh retrieved from an engine 34b ( Figure 1 ) produces a boundary determination ( Figure 2 (not shown) the evolution 35d of the particle distribution 34a in the LBM grid and the particle to the next cell The advection 35c process tests 37a for CFL constraint violations and, if a violation is found, uses the engine 36 ( Figure 1 ) Applying the correction 37b to the flux calculation.

[0051] Stabilized explicit numerical scheme for diffusion problems on irregular spatial grids

[0052] For the purpose of this specification, 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. The numerical scheme requires the calculation of the heat flux on all faces of an element. These fluxes are then accumulated and used to update the temperature of the element under consideration. Consider two adjacent elements α and β that share a face. According to Fourier's law of heat conduction, the heat flux is:

[0053]

[0054] in is the thermal conductivity at the common surface, The temperature gradient is normal to this common surface, and "m" is used to specify that the quantity is evaluated at time step "m". The negative sign of the commonly used Fourier law is removed because it is the heat entering α that is considered (not the heat leaving α). The temperature gradient used is calculated to ensure smooth advection of the particles, especially in the presence of elements of different sizes. If two adjacent elements α and β satisfy the CFL constraint, then the heat flux is multiplied by the area of ​​the common surface A. αβ and the time step size Δ t , to obtain the energy transfer across the common surface during time steps m to m+1, that is

[0055]

[0056] In the traditional approach, the final temperature of element α at the end of the time step is calculated based on the net energy transfer to α (the sum of the energy transfers from all faces):

[0057]

[0058] Note that equation (3) above shows that the temperature change is proportional to the net heat flux, which is proportional to the size of the element Inversely proportional, that is, for small elements, the same net energy transfer results in a larger temperature change.

[0059] refer to Figure 3 , the flux correction calculation process 36 includes the determination of the time step. In the flux calculation process 36, the time step is tested 36a. When the time step is not large enough to violate the CFL constraint of at least one of the two elements 36b, the above form (Equation 3) will not lead to numerical instability of the solution, and thus the flux calculation of the above form can be used. However, when the time step is large enough to violate the CFL constraint of at least one of the two elements 36c, the above form (Equation 3) may lead to numerical instability of the solution. In this case, the modified flux calculation process 36d is applied.

[0060] Without loss of generality, the revised method can be explained as follows:

[0061] Assume that element α is smaller than adjacent element β and thus violates the CFL constraint at least for element α. This numerical instability arises because for element α (assuming its size is smaller), the assumption is used to calculate The temperature gradient over the entire time step Δ t As mentioned above, for the same net energy transfer, the temperature of small elements varies more, so the standard explicit time integration requires a smaller time step to ensure that the constant temperature gradient assumption is valid. For a given Δ t , this problem persists as long as the instability exists in the problem and disappears only in a stable state where all incoming and outgoing fluxes on each element are exactly balanced with each other.

[0062] refer to Figure 4 As part of the revised approach, the flux calculation for the flux correction will be as defined in equation (1) Subdivision 49a is divided into two parts: (1a) applied flux 49b, which will be used (in the above summation) for the temperature evolution of α, and (1b) the balance flux 49c, which will be transferred to the other side of the interface αβ without changing the temperature of the element, and they are expressed as follows:

[0063]

[0064] Among them, and ΔG is given by:

[0065]

[0066]

[0067] as well as

[0068]

[0069] In the above equation, these geometric features are represented by the distance d used in calculating the temperature gradient. αβ , element volume and and the area A of the common surface (shared by adjacent elements) αβ The material properties are represented by and (Here ρ represents density and Cp represents specific heat) and thermal conductivity To illustrate. The flux term in the summation and An estimate of the fluxes that may exist on the different faces of the element α and β, respectively, is provided.

[0070] Two flux terms and The physical explanation of is as follows. The applied flux represents the total flux (given by equation (1)) The term can be used for the temperature evolution of element α without introducing numerical instabilities. This form can be derived from first principles for an isolated system consisting of elements α and β, thus including an estimate of the effect of the continuously evolving temperature field near the system. To this end, the term shows the geometrical / thermal properties of the element and the interaction of these elements with their surroundings (ΔG).

[0071] It should be noted that in equation (4), still relies on the fluxes observed at other surfaces at the previous time step to provide an estimate of the ongoing interactions at those surfaces during the current time step. For strongly transient problems, this leads to and The mismatch between.

[0072] The second term is called the balance flux To account for this mismatch, the flux is transported across α to its neighboring element on the other side of the interface αβ. This balance flux is only used in the temperature evolution when it is deposited in an element large enough to satisfy the CFL constraint, after which it is transported sequentially along the flux direction.

[0073] The above scheme strictly ensures that the correct amount of total flux is bound at the interface between any two elements α and β At the same time, the amount of flux available for the temperature evolution of the small element α is precisely controlled. In summary, this scheme is able to maintain numerical stability as well as good spatial and temporal accuracy. Finally, it should be noted that in the steady state, the applied flux becomes equal to the full flux, And thus the balance flux becomes equal to zero,

[0074] Characteristics of the modified flux calculation process

[0075] As mentioned in the above description, and as Figure 5-7As shown in , the modified flux calculation process 36 requires several algorithmic processes. The modified flux calculation process includes identifying the faces to which the modified definition of heat flux is to be applied 52a. At all faces that do not violate the CFL condition, the standard definition of heat flux (Equation 3) is used 52b. At all faces that violate the CFL condition, for example, at any face between two elements (at least one of which violates the CFL condition criteria), a modified heat flux calculation 52c is used. The total heat flux is calculated using the spatially averaged temperature gradient near the element under consideration to ensure smoothness of the solution. (In contrast, the standard heat flux utilizes the temperature gradient calculated based on the traditional difference form). The modified heat flux calculation process 52 includes dividing 52d the modified heat flux into two parts - the applied flux term 49b and the balance flux term 49c. The flux calculation 52e is then performed based on this division.

[0076] Now refer to Figure 6 , the flux calculation (Equation 3) using the applied flux calculation 52f is always used in the temperature evolution equation for the element under consideration.

[0077] Now refer to Figure 7 , the balance flux may or may not be used for temperature evolution, depending on the size of the elements involved. If the element is small enough to violate the CFL constraint 52g, the balance flux is calculated 52h and transferred 52i to the adjacent element in the flux direction, and if the adjacent element is large enough not to violate the CFL constraint, the balance flux is used 52j for temperature evolution. The balance flux is transferred successively along the flux direction until the balance flux is finally transferred to a large enough element (large enough not to violate the CFL constraint), at which element the balance flux is applied to the temperature evolution.

[0078] Algorithm uniqueness and advantages

[0079] Several approaches are known to overcome the problem of numerical instability in diffusion problems on irregular grids with varying element sizes. The most common approach is to impose additional constraints on the grid generation tool to reduce such scenarios. Even then, since the problem cannot be completely avoided, it is common practice to use a global time step that is small enough to ensure stability, or to use local subloops when encountering small grid elements. The first approach (small global time step) substantially increases the computational cost (even if there is a single occurrence of a small-sized element anywhere on the spatial grid), while the second approach (local subloops) increases the complexity of the algorithm and its implementation. An alternative approach is to use an implicit solution instead of an explicit solution, or at least to restrict the use of the explicit solution to (one or more) local regions near the small elements. This implicit solution approach suffers from implementation complexity and the non-local nature of the solution, which provides the system with equations that are inconvenient for large-scale parallel computer implementation.

[0080] In contrast, the modified flux calculation method provides several significant 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 in the grid). For every conceivable situation, this is a huge benefit in terms of computational cost and ease of implementation. The modified flux calculation method has a dependency on the geometric properties of two adjacent elements, thus ensuring that the modified flux calculation method will work properly regardless of the size of the elements. Therefore, the usual constraints on the grid generation process (grid quality, size, etc.) can be relaxed to a large extent. The computational cost of calculating various terms is reasonable because the mathematical form of the term is simple and does not involve any iteration. This is in stark contrast to existing methods that may have higher computational costs (reducing the time step size or using a hybrid implicit-explicit scheme). In addition, due to the volumetric nature of the formula, the scheme maintains accurate conservation, which is required in many applications. The modified flux calculation is still explicit in nature and requires information from elements within a small distance from the element being considered. This explicitness means that minimal changes to the existing computing system from the original system implementation scheme are required. The parallelization characteristics of the original explicit method are retained, so the modified flux calculation method can be implemented using large-scale parallel computers.

[0081] Therefore, the modified flux calculation method has advantages over existing methods. Although the above description of the modified flux calculation method is based on the finite volume formulation for heat conduction, the method is actually applicable to many diffusion-dominated problems.

[0082] Modeling the simulation space

[0083] In the physical process simulation system based on LBM, the fluid flow can be represented by a set of discrete velocities ci The distribution function value f to be evaluated i The dynamics of the distribution function is governed by equation (I.1), where f i (0) is known as the equilibrium distribution function, which is defined as:

[0084]

[0085] in,

[0086] 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 fluid block starts at a grid location and moves to the next grid location along one of the velocity vectors. At that point, the "collision operator" is calculated: that is, the influence of the nearby fluid block on the starting fluid block. Fluid can only move to one grid location, so the appropriate choice of velocity vectors is necessary so that all components of all velocities are multiples of a common speed.

[0087] 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 is Bhatnagar, Gross, and Krook (BGK). It forces the distribution function to a prescribed value given by the second equation, which is the "equilibrium" form.

[0088] The BGK operator is constructed based on the physical idea that, regardless of the details of the collision, the distribution function is close to the one given by {f eq (x,v,t)} gives a well-defined local equilibrium:

[0089]

[0090] where the parameter τ represents the characteristic relaxation time required to reach equilibrium via collisions. When dealing with particles (such as atoms or molecules), the relaxation time is usually considered a constant.

[0091] From this simulation, conventional fluid variables, such as mass ρ and fluid velocity u, are obtained as simple summations in equation (I3):

[0092]

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

[0094] 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 an LBE of the form:

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

[0096] 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 equations produce correct fluid dynamics and thermo-fluid dynamics. That is, from f i The (x, t)-derived hydrodynamic moments obey the Navier-Stokes equations in the macroscopic limit. These moments are defined by equation (I3) above.

[0097] c i and w i The LBM model can be implemented efficiently on scalable computer platforms and operates with maximum robustness to time-unstable flows and complex boundary conditions.

[0098] The standard technique for deriving the macroscopic equation 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 employed. In a fluid system, small perturbations in density travel at the speed of sound. In gaseous systems, the speed of sound is generally determined by the temperature. The significance of the effects of compressibility in a fluid is measured by the ratio of the characteristic velocity to the speed of sound, which is known as the Mach number.

[0099] The following provides a general discussion of LBM-based simulation systems that can be used in conjunction with CAD processes for fluid flow simulations.

[0100] refer to Figure 8 The first model (2D-1) 200 is a two-dimensional model including 21 velocities. Of the 21 velocities, one velocity (205) represents a particle that is not moving; three sets of four velocities represent particles that move at a normalized velocity (r) (210-213), twice the normalized velocity (2r) (220-223), or three times the normalized velocity (3r) (230-233) 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) (240-243) or twice the normalized velocity (2r) (250-253) relative to both the x and y lattice axes.

[0101] refer to Figure 9, shows a second model (3D-1) 260 - a three-dimensional model including 39 velocities, each of which is represented by Figure 9 Of these 39 velocities, one velocity represents a particle that does not move; 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.

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

[0103] 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 of 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).

[0104] 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 and y lattice axes; eight velocities represent particles that move at a normalized velocity (r) relative to one of the x and y lattice axes 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 and y lattice axes and at three times the normalized velocity (3r) relative to the other axis.

[0105] The LBM models described above provide a specific class of efficient and robust discrete velocity dynamics models for numerical simulation of flows in two and three dimensions. This type of model consists of a special set of discrete velocities and weights associated with those velocities. These velocities correspond to Cartesian grid points in velocity space, which facilitates accurate and efficient implementation of discrete velocity models, particularly the class known as lattice Boltzmann models. Using such models, flows can be simulated with high fidelity.

[0106] refer to Figure 10 The physical process simulation system operates according to process 270 using the above-mentioned CAD process to simulate a physical process such as fluid flow. Prior to simulation, the simulation space is modeled as a collection of voxels (step 272). 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.

[0107] The physical process simulation system operates according to process 270 using the modified flux calculation process described above. The resolution of the lattice can be selected 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 the objects in the flow (L), and the characteristic velocity of the flow (u):

[0108] Re=uL / v Equation (I4)

[0109] 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 ​​the 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 used around the area of ​​interest. The dimensionality of the voxel decreases as the resolution of the lattice increases.

[0110] 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 as f(x).

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

[0112] Each state i represents a different velocity vector at a specific 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:

[0113] c ii =(c ix ,c iy ,c iz ,). Equation (I5)

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

[0115] 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).

[0116] 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, six energy one states, eight energy three states, six energy four states, 12 energy eight states, and six energy nine states. By using this set of velocities, the system can generate Maxwell-Boltzmann statistics for the achieved equilibrium state vector.

[0117] 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 N i (n), where n represents the relative position of the lattice site in the microscopic block and n{0,1,2,...,7}.

[0118] Microscopic blocks in Figure 11 Shown in.

[0119] refer to Figure 12A and 12B , surface S( Figure 12A ) in the simulation space ( Figure 12B ) is represented as a surface element F α A collection of:

[0120] S={F α} Equation (I6) where α is an index that enumerates a particular facet. Facets are not restricted to the boundaries of a voxel, but typically have dimensions on the order of the dimensions of the voxels adjacent to the facet, or slightly smaller than the dimensions 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 (f i (α)). Total energy distribution function q i (α) is treated in the same way as for the flow distribution of surfels and voxel interactions.

[0121] refer to Figure 13 Different resolution levels can be used in different regions of the simulation space to improve processing efficiency. Typically, the region 320 surrounding the object 322 is of greatest interest and is therefore simulated using the highest resolution. Because the effect of viscosity decreases with distance from the object, regions 324, 326 spaced at increasing distances from the object 322 are simulated using reduced resolution levels (i.e., enlarged voxel volumes).

[0122] Similarly, if Figure 14 As shown in , lower resolution levels can be used to model regions 340 around less prominent features of an object 342, while the highest resolution level is used to model regions 344 around the most prominent features (e.g., leading and trailing edge surfaces) of the object 342. Outlying regions 346 are modeled using the lowest resolution level and the largest voxels.

[0123] C. Identify voxels affected by facets

[0124] Reference again Figure 10 Once the simulation space has been modeled (step 272), the voxels that are affected by one or more bins are identified (step 274). Voxels can be affected by bins in a number of ways. First, a voxel that is intersected by one or more bins is affected in that the voxel has a reduced volume relative to non-intersecting voxels. This occurs because the bin and the material underlying the surface represented by the bin occupy a portion of the voxel. The fractional factor P f (x) indicates the portion of the voxel that is not affected by the surface element (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.

[0125] Voxels that intersect one or more bins by either sending particles to them or receiving particles from them are also identified as voxels affected by the bins. All voxels intersected by a bin will include at least one state for receiving particles from the bin and at least one state for sending particles to the bin. In most cases, additional voxels will also include such states.

[0126] refer to Figure 15 , for a non-zero velocity vector c i For each state i, the surface element F α From the parallelepiped G iα The area defined by the particle receiving or transmitting is 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:

[0127] V ia =|c i n a |A a Equation (I7)

[0128] When the velocity vector of the state points to the surface element (|c i n i |<0), panel 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), transmit particles to this area. As will be discussed below, when another surface element occupies the parallelepiped G iα This expression should be modified when θ is part of a convex feature, i.e., a condition that occurs near non-convex features such as interior corners.

[0129] Surf F α parallelepiped G iα Multiple voxels can be partially or completely overlapped. The number of voxels, or portions thereof, depends on the bin size relative to the voxel size, the energy of the state, and the orientation of the bin relative to the lattice structure. The number of voxels affected increases with the size of the bin. Thus, as noted above, the size of a bin is typically selected to be of the order of magnitude of, or smaller than, the size of voxels located near the bin.

[0130] Parallelepiped G iα The portion of the voxel N(x) that overlaps is defined as V iα (x). Using this term, in the 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 element (N i (x)) multiplied by the volume of the region overlapping with the voxel (V iα (x)):

[0131] Γ ia (x) = N i (x)+V ia (x) Equation (I8)

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

[0133] V ia =∑V a (x)+∑V ia (β) Equation (I9)

[0134] The first summation takes into account the iα All voxels that overlap and the second term takes into account the iα All the face elements intersecting. When the parallelepiped Giα When not intersected by another polygon, this expression simplifies to:

[0135] V ia =∑V ia (x) Equation (I10)

[0136] D. Execute simulation

[0137] Once the voxels affected by one or more facets are identified (step 274), a timer is initialized to begin the simulation (step 276). 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 278-286). Next, a collision phase (step 288) simulates the interaction of particles within each voxel. Thereafter, the timer is incremented (step 200). If the incremented timer does not indicate that the simulation is complete (step 202), the advection phase and collision phase (steps 278-200) are repeated. If the incremented timer indicates that the simulation is complete (step 202), the results of the simulation are stored and / or displayed (step 204).

[0138] 1. Boundary conditions for surfaces

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

[0140] 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 should be equal to zero and the net normal momentum flux should be equal to the local pressure at the surface element. Therefore, the normal to the surface element, n α The components of the combined received and transmitted momentum that are perpendicular (i.e., the tangential components) should be equal and have a direction perpendicular to the normal n of the surface element. α The difference between the parallel components of the combined received and transmitted momentum (i.e., the normal components) should be equal to the local pressure at that cell. For a non-sliding 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.

[0141] 2. Aggregation from voxels to surfaces

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

[0143] Γ iα (x) = N i (x)V iα (x) Equation (I11)

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

[0145] Γ iαV→F =∑ X Γ iα (x)=∑ X N i (x)V iα (x) Equation (I12)

[0146] Only its V iα (x) The voxels with non-zero values ​​should 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.

[0147] 3. Moving from panel to panel

[0148] Next, the particles are moved between the bins (step 280). α The incoming state (c i n α <0) parallelepiped G iα By the other side element F β Intersect, then F α A portion of the particles in the received state i will come from the surface element F β In particular, the bin F α will receive the bin F during the previous time increment β The fraction of particles in state i that are generated.

[0149] Now refer to Figure 17 , which shows the bin F during the previous time increment β The relationship between the state i particle generated. Figure 17In the figure, the element F is shown. β Intersecting parallelepipeds G iα The part 380 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 the face element F α The flux of particles between states i can be described as:

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

[0151] 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:

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

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

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

[0155] Equation (I.15)

[0156] 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α :

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

[0158] The bin distribution function is a simulation tool used to generate output flux from bins and does not necessarily represent actual particles. To generate 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:

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

[0160] 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 utilize Γ from the previous time step iOTHER The value of (α) is generated so that:

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

[0162] 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 (α) expression) and V iα appears in the denominator (according to the i (α) expression). Therefore, when V iαand V iα When (x) is close to zero, N for parallel state i (α) is determined to be N i The values ​​of the states with zero velocity (i.e., the static 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.

[0163] 4. Performing Surfactant Surface Dynamics

[0164] Next, surface dynamics is performed on each of the bins to satisfy the four boundary conditions discussed above (step 282). The process for performing surface dynamics for a bin is described in Figure 18 Initially, the combined momentum P(α) of the particles at the bin is determined to determine the bin F α Combined momentum (step 392):

[0165] For all i, According to this equation, the normal momentum P n (α) is determined as:

[0166] P n (α) = n α P(α) Equation (I.20)

[0167] This normal momentum is then removed using a push / pull technique (step 394) 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.

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

[0169] Then, the flux distribution for the panel F is determined based on the incoming flux distribution, the modified flux calculation taking into account the CFL constraint violation, and the Boltzmann distribution. α The outgoing flux distribution of (step 398). First, the incoming flux distribution Γ i The difference between (α) and the Boltzmann distribution is determined as:

[0170] ΔΓ i (α) = Γ iIN (α)-N n-βi (α)Viα Equation (I.21)

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

[0172] For n α c i >0,Γ iOUT (α) = N n-βi (α)V iα -.Δ.Γ i* (α) Equation (I.22) 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:

[0173] For n α c i >0,

[0174] Among them C f is a function of skin friction, t iα is with n α The first perpendicular 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:

[0175]

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

[0177] 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 effectively produce 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 enforce the desired change in the tangential momentum flux at the surface. The friction coefficient C is 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.

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

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

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

[0181]

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

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

[0184] The outgoing flux distribution then becomes:

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

[0186] 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.

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

[0188]

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

[0190] 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:

[0191]

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

[0193] For cji n α >0,

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

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

[0196] This operation corrects the mass and energy fluxes while leaving the tangential momentum flux unchanged. If the flow in the neighborhood of the bin is roughly uniform and close to equilibrium, this adjustment is small. After the adjustment, the resulting normal momentum flux is slightly modified to the equilibrium pressure value based on the neighborhood average properties plus corrections due to non-uniform or non-equilibrium properties of the neighborhood. If the CFL constraint is violated, the process applies a modified flux calculation method to Figure 10 Any flux calculations involved in the process.

[0197] 5. Moving from voxel to voxel

[0198] Reference again Figure 10 , particles are moved between voxels along the three-dimensional rectilinear lattice (step 284). 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 constitute the majority of the voxels.

[0199] 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.

[0200] For voxels that do not interact with surfaces, the move operation is computationally very simple. The entire population of a state moves from its current voxel to the destination voxel during each time increment. At the same time, the particles of the destination voxel move from that voxel to their 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 in 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.

[0201] 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 transfer results in fractional particles remaining in the voxel. These fractional particles are then teleported to the voxel occupied by the surface.

[0202] refer to Figure 16 , when a portion 360 of the state i-particle of voxel 362 is moved to bin 364 (step 278), the remaining portion 366 is moved to voxel 368 in which bin 364 is located and from which the state i-particle is directed to bin 364. Thus, if the state population is equal to 25 and V iα (x) is equal to 0.25 (i.e., one quarter of the volume element and the parallelepiped G iα intersect), 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 transferred to a facet N(f) occupied by one or more facets is:

[0203]

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

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

[0206] Next, the outgoing particles from each bin are scattered to the voxels (step 286). Essentially, this step is the inverse of the aggregation step that moves 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:

[0207]

[0208] Among them, P f (x) shows that the volume of some voxels is reduced. Therefore, for each state i, the direction from the surface element to the voxel N is (x) The total number of particles is:

[0209]

[0210] After dispersing particles from bin to voxel, combining them with particles that have been advected from surrounding voxels, and integerizing the result, 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 gains or losses in mass, momentum, and energy after these quantities have been truncated to fit within the allowed value range.

[0211] To prevent this from happening, out-of-bounds mass, momentum, and energy are accumulated before the offending state is truncated. For the energy to which the state belongs, an amount of mass equal to the value gained (due to underflow) or lost (due to overflow) is added back to a randomly (or sequentially) selected state of the same energy that 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 states of the same energy, both mass and energy are corrected when the mass counter reaches zero. Finally, momentum is corrected using a push / pull technique until the momentum accumulator returns to zero.

[0212] 7. Perform fluid dynamics

[0213] Perform fluid dynamics (step 288) Figure 10 This step can be called microdynamics or intra-voxel operations. Similarly, the advection process can be called inter-voxel operations. The microdynamics operation described below can also be used to collide particles in the bin to produce a Boltzmann distribution.

[0214] 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 simulates the dynamics of distributions in real fluid systems. The collision process can be well described by the right-hand side of Equations 1 and 2. After the advection step, the conserved quantities of the fluid system, specifically density, momentum and energy, are obtained from the distribution function using Equation 3. From these quantities, the f in Equation (2) is obtained. eq The equilibrium distribution function is fully defined by equation (4). The velocity vector set c i The choices of , weights are listed in Table 1, together with Equation 2 to ensure that the macroscopic behavior follows the correct fluid dynamics equations.

[0215] Variable resolution

[0216] Variable resolution (as discussed in US 2013 / 0151221 A1) may also be employed and voxels of different sizes would be used, eg, coarse and fine voxels.

[0217] By leveraging unique transient physics principles based on lattice Boltzmann theory, the system can perform simulations that accurately predict real-world conditions. For example, engineers can evaluate product performance early in the design process, before building any prototypes, when the impact of changes is significant to the design and budget. The system can accurately and efficiently perform aerodynamic, aeroacoustic, and thermal management simulations using CAD geometry. The system can perform simulations to address the following applications: Aerodynamics (aerodynamic efficiency; vehicle handling; pollution and water management; panel deformation; driving dynamics), Aeroacoustics (greenhouse wind noise; underbody wind noise; gap / seal noise; mirror, whistle and tonal noise; sunroof and window buffeting; bypass / neighborhood noise; cooling fan noise), Thermal Management (cooling airflow; thermal protection; brake cooling; drive cycle simulation; shutdown and soak; electronics and battery cooling; ROA / air intake), Climate Control (cabin comfort; HVAC unit and power distribution system performance; HVAC system and fan noise; defrost and defogging), Powertrain (powertrain cooling; exhaust system; cooling jacket; engine block), Wastewater and Water Management (column overflow, dirt and dust accumulation, tire spray).

[0218] Embodiments of the subject matter and functional operations described in this specification may be implemented in digital electronic circuitry, tangibly embodied computer software or firmware, computer hardware (including the structures disclosed in this specification and their structural equivalents), or a combination of one or more of them. Embodiments of the subject matter described in this specification may be implemented as one or more computer programs (i.e., one or more modules of computer program instructions encoded on a tangible, non-transitory program carrier for execution by a data processing apparatus or for controlling the operation of a data processing apparatus). A computer storage medium may 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 of them.

[0219] The term "data processing apparatus" refers to data processing hardware and includes all types of apparatus, devices, and machines for processing data, including, for example, a programmable processor, a computer, or multiple processors or computers. The apparatus may also be or further include special-purpose logic circuitry (e.g., an FPGA (field programmable gate array) or an ASIC (application-specific integrated circuit)). In addition to hardware, the apparatus may also include code that creates an execution environment for a computer program (e.g., code constituting processor firmware, a protocol stack, a database management system, an operating system, or a combination of one or more of these).

[0220] A computer program may also be referred to or described as a program, software, software application, module, software module, script or code and may be written in any form of programming language, including compiled or interpreted languages ​​or declarative or procedural languages, and it may be deployed in any form, including as a standalone program or as a module, component, subroutine or other unit suitable for a computing environment. A computer program may, but need not, correspond to a file in a file system. A program may be stored in a portion of a file that holds other programs or data (e.g., one or more scripts stored in a markup language document, a single file dedicated to the program in question, or in multiple coordinated files (e.g., files storing one or more modules, subroutines or code portions)). A computer program may be deployed so that the program is executed on one computer or on multiple computers located at one site or distributed across multiple sites and interconnected by a data communications network.

[0221] Computers suitable for executing computer programs may be based on general or special purpose microprocessors, or both, or any other type of central processing unit. Computer-readable media suitable for storing computer program instructions and data include all forms of non-volatile memory on media and memory devices, including, for example, semiconductor memory devices (e.g., EPROM, EEPROM, and flash memory devices), magnetic disks (e.g., internal hard disks or removable disks), magneto-optical disks, CD ROMs, and DVD-ROM disks. The processor and memory may be supplemented by, or incorporated in, special purpose logic circuitry.

[0222] Embodiments of the subject matter described in this specification can be implemented in a computing system that includes a back-end component (e.g., as a data server), or includes a middleware component (e.g., an application server), or includes a front-end component (e.g., a client computer with a graphical user interface or a web browser through which a user can interact with an implementation of the subject matter described in this specification), or one or more such back-end, middleware, or front-end components. The components of the system can be interconnected by any form or medium of digital data communication (e.g., a communication network). Examples of communication networks include local area networks (LANs) and wide area networks (WANs) (e.g., the Internet).

[0223] A computing system may include a client and a server. The client and the server are typically remote from each other and typically interact via a communication network. The relationship between the client and the server arises by means of computer programs running on respective computers and having a client-server relationship with each other. In some embodiments, the server sends data (e.g., an HTML page) to a user device (e.g., for the purpose of displaying data to a user interacting with the user device and receiving user input), and the user device acts as a client. Data generated at the user device (e.g., the result of a user interaction) may be received at the server from the user device.

[0224] Specific embodiments of the subject matter have been described. Other embodiments are within the scope of the following claims. For example, the actions recited in the claims can be performed in a different order and still achieve the desired results. As an example, the processes depicted in the accompanying drawings do not necessarily require the particular order shown, or sequential order, to achieve the desired results. In certain circumstances, multitasking and parallel processing may be advantageous.

Claims

1. A computer-implemented method for simulating fluid flow around a solid surface, the method comprising: Receiving, by one or more computing systems, a model of a simulation space, the simulation space comprising a lattice structure and a representation of a physical object, the lattice structure represented as a collection of voxels, wherein the voxels have a resolution selected based on a surface of the physical object; simulating, by the one or more computing systems, motion of particles in a fluid volume, wherein the motion of the particles causes collisions between the particles; identifying, by the one or more computing systems, faces between two voxels, wherein at least one of the faces violates a Courant-Friedrichs-Lewy constraint, wherein a standard heat flux definition applies when the Courant-Friedrichs-Lewy constraint is not violated and a modified heat flux definition applies when the Courant-Friedrichs-Lewy constraint is violated; According to a modified heat flux definition, the one or more computing systems calculate a modified heat flux using a spatially averaged temperature gradient near the two volume elements, wherein the at least one of the surfaces violates the Courant-Friedrichs-Lewy constraint, wherein calculating the modified heat flux comprises: calculating an applied flux by the one or more computing systems, wherein the calculated applied flux is used for the temperature evolution of the considered voxel, calculating, by the one or more computing systems, a balance flux as the difference between the total flux and the applied flux, and The one or more computing systems transfer the balance flux to one or more adjacent voxels along the direction of flux flow until the balance flux is finally transferred to a specific voxel that is large enough to satisfy the Courant-Friedrichs-Lewy constraint, and the balance flux is applied to the temperature evolution at the specific voxel; performing, by the one or more computing systems, an advection operation on the particles to reach a subsequent voxel; and Results of the simulation are stored and / or displayed by the one or more computing systems.

2. The computer-implemented method of claim 1 , further comprising: At adjacent faces where none of the faces violates the Courant-Friedrichs-Lewy constraint, a temperature gradient based on the temperature difference is calculated.

3. The computer-implemented method of claim 2, wherein: The temperature difference represents the net energy transfer to each of the voxels to calculate the final temperature of the voxel at the end of a time period having a time step, the net energy transfer corresponding to the sum of the energy transfers from all faces of the voxel.

4. An apparatus for simulating fluid flow around a solid surface, the apparatus comprising: Memory; One or more processor devices configured to: Performing a computer-implemented method according to any one of claims 1-3.

5. A computer program product tangibly stored on one or more machine-readable hardware storage devices, comprising executable instructions executable by one or more processing devices to cause a computer to perform the computer-implemented method according to any one of claims 1 to 3.

Images