Reducing artificial compressibility in digital internal fluid flow simulation

Abstract

The invention relates to reducing artificial compressibility in digital internal fluid flow simulation. Systems and methods for digital simulation of fluid flow in a three-dimensional computer-aided design (CAD) model of a simulated space include receiving a digital representation of the simulated space based on the digital three-dimensional CAD model. The digital representation includes a plurality of voxels. The first fluid flow is digitally simulated by applying a driving force to the plurality of voxels in the digital representation of the simulation space to generate a pressure field. A volumetric force to be applied to the plurality of voxels is determined based on a pressure gradient of the pressure field. The second fluid flow is digitally simulated by applying a volumetric force to the plurality of voxels.

Description

[0001] Cross-reference to related applications

[0002] This application claims the benefit of U.S. Patent Application No. 18 / 749,798, filed June 21, 2024, which is incorporated herein by reference in its entirety. Background Technology

[0003] This manual relates to simulating physical processes, such as fluid flow.

[0004] Fluid flow can be simulated by performing high-precision floating-point arithmetic on variables representing macroscopic physical quantities (e.g., density, temperature, flow velocity) at each of many discrete spatial locations, generating discretized solutions to the Navier-Stokes differential equations. Another approach replaces the differential equations with a technique commonly known as lattice gas (or cellular) automata, where the macroscopic-level simulation provided by solving the Navier-Stokes equations is replaced by a microscopic-level model that operates on particles moving between points on the lattice. Summary of the Invention

[0005] Compressible fluid flow solvers (e.g., Lattice Boltzmann Model (LBM) solvers) can be used for digital simulation and research of fluid flows in the automotive, aerospace, and energy industries. LBM solvers can offer high robustness, accuracy, scalability, and versatility, enabling them to be used in many different scenarios to simulate fluid flows within and around physical objects with arbitrarily complex geometries. For example, LBMs can be applied to internal fluid flows in objects with tortuous structures, such as porous media (e.g., carbonate rock cores) or gas diffusion layers (GDLs) in fuel cells.

[0006] When generating a digital representation of a physical object or system for fluid flow simulations, configuration parameters (e.g., dimensionless numbers) are selected for the digital representation to match the physical object or system. Ideally, all configuration parameters will correspond perfectly to the physical parameters; however, in practice, assumptions and / or adjustments are made to reduce computational costs and time. For example, when simulating fluid flow through the interior of a pipe, the Mach number (the ratio of the fluid velocity to the speed of sound) may be very small, potentially requiring millions of time steps to perform the simulation, which can be daunting. Alternatively, the Mach number can be artificially increased for the simulation to reduce the number of time steps and the computational resources required. Artificially increasing the Mach number can have unintended consequences, such as introducing compressibility effects (e.g., artificial compressibility) into incompressible fluid flows. Artificial compressibility effects reduce the accuracy and quality of the simulation results.

[0007] The present disclosure describes a method for simulating internal fluid flow using a compressible flow solver. To counteract artificial compressibility effects, a data processing system can determine a conformal body force for driving the internal fluid flow to decouple pressure and density terms. The data processing system can determine the conformal body force by performing a first simulation (e.g., a preparatory simulation) using an artificially elevated Mach number. The conformal body force can be determined based on a pressure gradient from the preparatory simulation. The data processing system can perform a second simulation (e.g., a main simulation) of the fluid flow with the conformal body force applied to drive the fluid flow.

[0008] In an example implementation, a computer system for numerically simulating fluid flow in a three-dimensional computer-aided design (CAD) model of a simulation space includes one or more processors; and a memory including: a mesh preparation engine to generate and store a digital representation of the simulation space based on the digital three-dimensional CAD model, the digital representation including a mesh comprising a plurality of voxels; and a simulation engine to read the digital representation of the mesh in the simulation space from the mesh preparation engine. The simulation engine stores instructions for numerically simulating fluid flow that, when executed by the one or more processors, cause the one or more processors to perform operations including: reading the digital representation of the mesh in the simulation space from the mesh preparation engine; numerically simulating a first fluid flow by applying a driving force to the plurality of voxels in the digital representation of the mesh in the simulation space to generate a digital pressure field; determining a digital volumetric body force to be applied to the plurality of voxels based on a pressure gradient of the digital pressure field; and numerically simulating a second fluid flow by applying the digital volumetric body force to the plurality of voxels.

[0009] In another example implementation, a method implemented by a data processing system for numerically simulating fluid flow in a three-dimensional CAD model of a simulation space includes: receiving, by the data processing system, a digital representation of the simulation space based on a digital three-dimensional CAD model, the digital representation including a plurality of voxels; numerically simulating, by the data processing system, a first fluid flow by applying a driving force to the plurality of voxels in the digital representation of the simulation space to generate a pressure field; determining, by the data processing system, a volumetric body force to be applied to the plurality of voxels based on a pressure gradient of the pressure field; and numerically simulating, by the data processing system, a second fluid flow by applying the volumetric body force to the plurality of voxels.

[0010] In another example implementation, one or more non-transitory machine-readable storage devices store instructions for numerically simulating fluid flow in a three-dimensional CAD model of a simulation space, the instructions executable by one or more processors to cause performance of operations comprising: receiving a digital representation of the simulation space based on a digital three-dimensional CAD model, the digital representation comprising a plurality of voxels; numerically simulating a first fluid flow by applying a driving force to the plurality of voxels in the digital representation of the simulation space to generate a pressure field; determining a volumetric body force to be applied to the plurality of voxels based on a pressure gradient of the pressure field; and numerically simulating a second fluid flow by applying the volumetric body force to the plurality of voxels.

[0011] In aspects combinable with one, some, or all of the example implementations, the first fluid flow and the second fluid flow comprise incompressible, steady-state fluid flow.

[0012] In another aspect combinable with one, some, or all of the preceding aspects, the digital volumetric body force comprises a digital conformal volumetric body force aligned with a streamwise direction.

[0013] In another aspect combinable with one, some, or all of the preceding aspects, the second fluid flow has a smaller change in density than the first fluid flow.

[0014] Another aspect combinable with one, some, or all of the preceding aspects comprises determining an actual pressure force for the second fluid flow based on the driving force from the first fluid flow and the measured pressure from the second fluid flow.

[0015] In another aspect combinable with one, some, or all of the preceding aspects, a number of voxels used in simulating the first fluid flow is different than a number of voxels used in simulating the second fluid flow.

[0016] Another aspect combinable with one, some, or all of the preceding aspects comprises storing the pressure field and the volumetric body force in a hardware storage device for subsequent simulation of fluid flow.

[0017] In another aspect combinable with one, some, or all of the preceding aspects, the pressure field comprises a force balance based on driving forces in the simulation space, fluid inertia, and friction at solid boundaries.

[0018] Another aspect, which can be combined with one, some, or all of the preceding aspects, includes storing, in a memory, a digital pressure field generated by digitally simulating the first fluid flow by applying a driving force to a plurality of voxels in a digital representation of a mesh in an analog space; storing, in the memory, a digital volume body force to be applied to the plurality of voxels, the digital volume body force determined based on a pressure gradient of the digital pressure field; and storing, in the memory, a result of a digital simulation of the second fluid flow generated by applying the digital volume body force to the plurality of voxels.

[0019] One or more of the above aspects can provide one or more of the advantages disclosed herein. The method reduces the computational resources and computational complexity required to simulate incompressible fluid flow. The method corrects for artificially generated compressibility of the fluid in the internal flow simulation when using high simulation Mach numbers. High simulation Mach numbers enable the simulation to be performed using fewer time steps, thereby reducing the number of computations required to simulate the fluid flow. The compressible fluid solver uses fewer computations than an incompressible fluid solver because the compressible fluid solver neither iterates steps nor solves many sets of linear equations to determine a pressure field of the fluid flow. The method both reduces the computational steps and thus the computational resources required and enables the simulation to be parallelized on a parallel computing system.

[0020] This reduction in computational complexity saves computational resources because less processing power is required to perform the computations relative to the amount of processing power required for more complex computations. This reduction in computational complexity also increases the speed at which the processing device performs the computations. Generally, processing power includes the ability of a computer (or processing device) to process data. The method improves the accuracy of the incompressible fluid flow simulation by decoupling density and pressure to reduce artificially generated compressibility by applying conformal body forces to voxels of a simulation grid. By varying a global factor of the conformal body forces, the conformal body forces can be adapted to different sized boundary setups, such as different pressure driving forces. The reusability of the conformal body forces is advantageous because it improves computational efficiency for future digital simulations.

[0021] Other features and advantages of the method will be apparent from the following detailed description, and from the claims. BRIEF DESCRIPTION OF DRAWINGS

[0022] Figure 1 A system for digital simulation of fluid flow is depicted.

[0023] Figure 2A A flowchart showing operations for formulating a lattice-Boltzmann model simulation is depicted.

[0024] Figure 2BA flowchart of a method of depicting the effect of reducing artificial compressibility in internal fluid flow simulations.

[0025] Figure 3 and Figure 4 Illustrates velocity components of two LBM's represented in Euclidean space (prior art).

[0026] Figure 5 Is a flowchart of the process followed by a physical process simulation system using corrected computer aided design (CAD) drawings.

[0027] Figure 6 Is a perspective view of a microblock (prior art).

[0028] Figures 7A-7B Is Figure 1 A diagram of the lattice structure used by the system of

[0029] Figure 8 and Figure 9 Each illustrates a variable resolution technique (prior art).

[0030] Figure 10 Illustrates movement of particles (prior art).

[0031] Figure 11 Illustrates an area affected by a facet of a surface (prior art).

[0032] Figure 12 Illustrates movement of particles from surface to surface (prior art).

[0033] Figure 13 Is a flowchart of the process for performing surface dynamics operations.

[0034] Figures 14A-14F Shows simulation results for water flow driven by pressure through bipolar channels of a fuel cell.

[0035] The details of one or more embodiments of the application are set forth in the accompanying drawings and the description below. Other features, objects, and advantages of the application will be apparent from the description and drawings, and from the claims. DETAILED DESCRIPTION

[0036] One method for simulating fluid flow is the so-called Lattice Boltzmann method. LBM-based simulations have many advantages, including robustness, compatibility with multiphase and multi-component flows, localized computation, high scalability, and grid-independent solutions. In an LBM-based physical process simulation system, fluid flow is represented by distribution function values evaluated at a set of discrete velocities using well-known lattice Boltzmann equations that describe the time evolution of the distribution function. The distribution function involves two processes: a streaming process and a collision process.

[0037] Figure 1 A schematic diagram of an example data processing system 10 that performs a lattice Boltzmann (LB)-based simulation is illustrated. The system 10 in this implementation is based on a client-server or cloud-based architecture and includes a server system 12 implemented as a large-scale parallel computing system 16 (stand-alone or cloud-based) and a client system 14 coupled via a network 15. The server system 12 includes a memory 18, a bus system 22, 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 mesh preparation engine 32 and a simulation engine 34.

[0038] While Figure 1 The mesh preparation engine 32 in the memory 18 is shown, but the mesh preparation engine can be a third-party application executing on a system different from the server 12. Whether the mesh preparation engine 32 is executing in the memory 18 or on a system different from the server 12, the mesh preparation engine 32 receives a user-supplied mesh definition 30 based on a CAD-generated drawing 31 and then prepares the mesh and sends (and / or stores) the prepared mesh to the simulation engine 34.

[0039] The simulation engine 34 includes a collision interaction module 34a that includes a surface dynamics conversion 34b, a boundary handling module 34c, and an advection operation 34d. The system 10 accesses a data repository 38 that stores 2D and / or 3D meshes (Cartesian and / or curvilinear), coordinate systems, and libraries.

[0040] Referring Figure 2AFIG. 4 shows a process 40 for simulating fluid flow through a representation of a physical object. Examples of physical objects include porous media (such as carbonate rock), piping systems, GDLs of fuel cells, etc. The physical object can be of any shape, and in particular can have planar and / or curved surfaces. The process 40 receives 42 a grid (or mesh) for the physical object being simulated, e.g., from a client system 14 or retrieves the grid from a data repository 38. In other embodiments, an external system or server 12 generates the grid for the physical object being simulated based on user input. The process pre-computes 44 geometric quantities from the retrieved grid and performs a dynamic lattice Boltzmann model simulation 46 using the pre-computed geometric quantities corresponding to the retrieved grid. The lattice Boltzmann model simulation includes a simulation 46 of the evolution of a particle distribution, including surface dynamic conversion, boundary modeling, and particle advection transport into the next cell of the LBM grid.

[0041] Compressible fluid flow solvers, such as LBM solvers, can be used to numerically simulate fluid flow with high robustness, accuracy, scalability, and versatility, enabling LBM solvers to be used in many different scenarios to simulate fluid flow in and around physical objects with arbitrary complex geometrical configurations.

[0042] In generating a digital representation of a physical object and fluid flow configurations associated with the physical object, dimensionless numbers such as Reynolds number (ratio of inertial forces to viscous forces), Bond number (ratio of gravitational forces to surface tension forces), and Mach number (ratio of flow velocity to speed of sound in the fluid) can be used to characterize fluid forces that control fluid dynamics. Dimensionless numbers are particularly useful for scaling fluid dynamics to match between various length scales, fluid properties, and flow velocities. In an ideal situation, dimensionless numbers can be matched without assuming control dynamics. However, in most situations, this can be very difficult to achieve due to limitations of computer resources. For example, to simulate water flow at 0.8 meter / second through a circular pipe with length 1 m and radius 0.01 m, the Mach number is 5.33e-4. At this Mach number, it would take 1 million to 2 million time steps to make the fluid go through the pipe, with radius and speed of sound being 10 grid (e.g., each voxel in the digital representation has a side length of 0.001 m), and 1.0 grid per time step (e.g., 6.67e-7 seconds per time step). Such a small time step can result in unacceptably huge computational cost. More complex geometrical configurations of the physical object further exacerbate this problem. A common workaround to solve this problem is to increase the simulation Mach number within a reasonable range. However, increasing the simulation Mach number can result in increased numerical compressibility effects.

[0043] An increase in the simulated Mach number results in higher fluid velocities, which in turn results in higher frictional forces and larger pressure changes throughout the conduit. In compressible flow solvers, such as LBM, the pressure of the fluid is related to the density. Higher pressure changes result in density changes that do not reflect real-world fluid flow. The excessive compressibility effects resulting from an increased simulated Mach number can cause problems because the compressibility effects reduce the accuracy of fluid flow field predictions.

[0044] To overcome the artificial compressibility effects introduced by increasing the simulated Mach number, a conformal body force can be used to drive the fluid flow and decouple the density and pressure changes. A conformal body force can be generated to replicate the fluid flow dynamics of a pressure driven flow while significantly reducing the artificial compressibility effects. Two simulation approaches (preparation simulation and main simulation) can be used to develop an appropriate conformal body force (e.g., a force acting on a fluid element along the flow direction of the fluid flow). In the preparation simulation, only a single phase pressure driven flow can be simulated to determine the pressure profile along the flow path. For example, the pressure profile can be determined by a force balance between the pressure, inertial forces, and friction from solid boundaries (e.g., pipe walls, pore walls). A volume conformal body force field is generated based on the pressure gradient from the pressure profile of the preparation simulation. In the main simulation, the previously generated conformal body force field can be applied to the fluid flow in place of the pressure driving force to serve as the main driving force for the fluid. The resulting pressure is represented by the volume conformal body force plus a density dependent pressure term. This decoupling can reduce the size of the density dependent pressure term such that the density changes along the flow path are reduced, thereby reducing the artificial compressibility effects.

[0045] Using this conformal volume body force approach, simulating incompressible fluid flow using a compressible fluid flow solver improves the accuracy of fluid flow predictions compared to simulating using only an artificially increased Mach number. At the same time, compared to simulating incompressible fluid flow at the actual fluid flow Mach number, this approach reduces the computational cost by reducing the number of calculations needed to accurately simulate the fluid flow. This approach solves the compressible Navier-Stokes equations or an asymptotically equivalent equation.

[0046] Figure 2B A flowchart depicting a process 100 for simulating internal fluid flow is depicted. A data processing system receives (102) a digital representation of a simulation space, the digital representation containing a plurality of voxels. For example, the data processing system receives the digital representation from the mesh preparation engine 32 or from the data repository 38. In some implementations, the data processing system generates the digital representation of the simulation space based on a CAD model or drawing representing a physical object.

[0047] The data processing system simulates (104) a first fluid flow by applying driving forces to a plurality of voxels in a digital representation of a simulation space to generate a pressure field. Driving forces include forces that cause fluid to flow from one location to another, such as pressure, gravity, and inertial forces. The data processing system can apply boundary conditions that reflect real-world boundary conditions and run the simulation until the simulation results have converged. For example, the data processing system can simulate a pressure-driven flow by applying pressure driving forces to the voxels in the simulation space.

[0048] The data processing system determines (106) a volumetric body force to be applied to the plurality of voxels based on pressure gradients of the pressure field. For example, the data processing system determines a conformal body force acting in a flow direction of the fluid flow by determining pressure gradients from the pressure field of the simulation of the first fluid flow. The data processing system can store the pressure field and the volumetric body force in a data repository or other hardware storage device for access at a later time.

[0049] The data processing system simulates (108) a second fluid flow by applying the volumetric body force to the plurality of voxels. The simulation of the first fluid flow and the simulation of the second fluid flow occur in the same geometry. In some embodiments, the simulation of the first fluid flow and the simulation of the second fluid flow can have different resolutions (e.g., different numbers of voxels). For example, the first fluid flow can be simulated at a lower resolution than the second fluid flow. The volumetric body force can be scaled up or down to the resolution of the second fluid flow (e.g., by bi-cubic or linear interpolation). In some implementations, the data processing system scales the magnitude of the volumetric body force to achieve a desired fluid flow.

[0050] The data processing system can determine an actual pressure for the second fluid flow based on the driving forces from the first fluid flow and the measured pressure from the second fluid flow. The measured pressure can be a density-dependent pressure. Both the first fluid flow and the second fluid flow can be incompressible steady-state fluid flows.

[0051] In some implementations, the data processing system stores the volumetric body force and the pressure field in a hardware storage device for use in additional (e.g., subsequent) simulations of fluid flows with the same geometry.

[0052] Detailed Examples

[0053] In the processes discussed below in Figure 5 , CAD drawings are used to describe the flow simulation process, with identified void spaces configuring the simulation space. In the previously discussed Figure 3 and Figure 4 and Figure 6 , Figure 7A , Figure 7B, Figure 8 , Figure 9 , Figure 10 , Figure 11 and Figure 12 In these figures, each of these figures is labeled as prior art because these figures appear in U.S. Patent 5,848,260 (the '260 patent) or U.S. Patent 11,847,391 (the '391 patent), both of which are hereby incorporated in their entirety herein.

[0054] However, the figures appearing in the above patents do not take into account any modification made to the effects of artificial compressibility for fluid simulation because the processes described herein are not described in the above cited patents.

[0055] Modeling the simulation space

[0056] In an LBM-based physical process simulation system, fluid flow is represented by distribution function values evaluated at a set of discrete velocities. The dynamics of the distribution function is governed by the lattice Boltzmann equation that relates the change of distribution due to a so-called "streaming process" to the change of the distribution function due to a "collision process". The streaming process is when a fluid pocket starts at a grid location and then moves to the next grid location along one of a number of velocity vectors. At this time, a "collision factor" is calculated, i.e., the effect of nearby fluid pockets on the starting fluid pocket. The fluid can only move to another grid location, so it is necessary to select the velocity vectors correctly so that all components of all velocities are multiples of a common velocity. The collision process uses a "collision operator" to represent the change of the distribution function due to collisions between fluid pockets. A particular form of the collision operator is the Bhatnagar, Gross, and Krook (BGK) operator. The collision operator forces the distribution function to tend toward a prescribed value.

[0057] The BGK operator is constructed on the basis of the following physical argument: regardless of the details of the collision, the distribution function approaches a well-defined local equilibrium via the collision, reaching equilibrium via the collision according to a characteristic relaxation time. When dealing with particles (e.g., atoms or molecules), the relaxation time is often taken as a constant.

[0058] From this simulation, the regular fluid variables, such as mass and fluid velocity, are obtained by simply summing the products of the distributions. Due to symmetry considerations, a set of velocity values is chosen such that they form certain lattice structures when striding across the configuration space. The dynamics of this discrete system obeys the LBE, where the collision operator is usually in the BGK form as described above. By choosing the equilibrium distribution form properly, it can be shown theoretically that the lattice Boltzmann equation produces the correct fluid dynamics and thermal hydrodynamics. That is, the fluid momentum derived from the distribution function obeys the Navier-Stokes equations in the macroscopic limit.

[0059] The collective values of the lattice velocities and associated weights define the LBM. The LBM can be implemented efficiently on scalable computer platforms and runs with high robustness for time- dependent flow and complex boundary conditions.

[0060] The standard technique to obtain the macroscopic equations of motion for fluid systems from the Boltzmann equation is the Chapman-Enskog method, where a successive approximation of the full Boltzmann equation is taken. In a fluid system, a small perturbation in density travels at the speed of sound. In a gas system, the speed of sound is generally determined by the temperature. The importance of compressibility effects in the flow is measured by the ratio of a characteristic velocity, called the Mach number, to the speed of sound.

[0061] A general discussion of the LBM-based simulation system, including the dynamics transform 34b to perform fluid flow simulation, is provided below. For further explanation of the LBM-based physical process simulation system, the reader is referred to the '260 patent.

[0062] Reference Figure 3 The first model (2D-1) 200 is a two-dimensional model that includes 21 velocities. Of these 21 velocities, one velocity (205) represents a particle that is not moving; three sets of four velocities represent particles that are moving in the positive or negative direction along the x-axis or y-axis of the lattice at a normalized rate (r) (210-213), twice the normalized rate (2r) (220-223), or three times the normalized rate (3r) (230-233); and two sets of four velocities represent particles that are moving relative to both the x and y lattice axes at the normalized rate (r) (240-243) or twice the normalized rate (2r) (250-253).

[0063] Reference Figure 4 A second model (3D-1) 260 is illustrated - a three-dimensional model that includes 39 velocities, where each velocity is represented by Figure 4One of the arrows represents. Out of the 39 velocities, one velocity represents a particle that is not moving; three sets of six velocities represent particles that are moving in the positive or negative direction along the x, y, or z axis of the lattice at a normalized rate (r), twice the normalized rate (2r), or three times the normalized rate (3r); eight velocities represent particles that are moving at the normalized rate (r) relative to all three of the x, y, z lattice axes; and twelve velocities represent particles that are moving at twice the normalized rate (2r) relative to two of the x, y, z lattice axes.

[0064] More complex models can also be used, such as a 3D-2 model including 101 velocities and a 2D-2 model including 37 velocities. For the three-dimensional model 3D-2, out of the 101 velocities, one velocity represents a particle that is not moving (set 1); three sets of six velocities represent particles that are moving in the positive or negative direction along the x, y, or z axis of the lattice at a normalized rate (r), twice the normalized rate (2r), or three times the normalized rate (3r) (sets 2, 4, and 7); three sets of eight velocities represent particles that are moving at the normalized rate (r), twice the normalized rate (2r), or three times the normalized rate (3r) relative to all three of the x, y, z lattice axes (sets 3, 8, and 10); twelve velocities represent particles that are moving at twice the normalized rate (2r) relative to two of the x, y, z lattice axes (set 6); twenty-four velocities represent particles that are moving at the normalized rate (r) and twice the normalized rate (2r) relative to two of the x, y, z lattice axes and not moving relative to the remaining axis (set 5); and twenty-four velocities represent particles that are moving at the normalized rate (r) relative to two of the x, y, z lattice axes and at three times the normalized rate (3r) relative to the remaining axis (set 9).

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

[0066] The LB models described above provide a particular class of efficient and robust discrete velocity kinetic models for numerical simulation of flow in both two and three dimensions. This type of model includes a particular set of discrete velocities and weights associated with the velocities. The velocities correspond to grid points of a Cartesian coordinate in a velocity space, which facilitates accurate and efficient implementation of the discrete velocity model, particularly that known as the lattice Boltzmann model. Using such models, flow can be simulated with high fidelity.

[0067] Reference is made to Figure 5 A physical process simulation system is described that operates according to the process 270 to simulate a physical process such as fluid flow. Prior to the flow simulation, the simulation space is modeled (272) as a collection of voxels using CAD drawings discussed above. The simulation space is generated using a computer aided design (CAD) program and gap correction processing of CAD generated drawings. For example, the CAD program can be used to draw an airfoil located in a wind tunnel.

[0068] 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, the characteristic length of the objects in the flow, and the characteristic velocity of the flow.

[0069] The characteristic length of the objects represents the large scale features of the objects. For example, if a flow around a micro device is being simulated, the height of the micro device can be considered the characteristic length. When the flow around a small region of an object is of interest (e.g., the side mirror of a car), the resolution of the simulation can be increased or a region of increased resolution can be employed around the region of interest. The size of the voxels decreases as the resolution of the lattice is increased.

[0070] The state space is represented as a distribution function of particles or particles per unit volume at a given state at a lattice site represented by a spatial vector at a given time. The number of states is determined by the number of possible velocity vectors within each energy level. The velocity vectors are integer linear velocities in a space having three dimensions, x, y, and z. For multi-species simulations, the number of states is increased. Each state represents a different velocity vector at a particular energy level (i.e., energy level zero, one, or two). The velocity of each state is indicated by its "velocity" in each of the three dimensions.

[0071] The energy level zero state represents a stopped particle that is not moving in any dimension, i.e., the particle has a velocity of zero in each dimension. The energy level one state represents a particle that has a ±1 velocity in one of the three dimensions and zero velocity in the other two dimensions. The energy level two state represents a particle that has a ±1 velocity in all three dimensions or a ±2 velocity in one of the three dimensions and zero velocity in the other two dimensions.

[0072] 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). Each voxel (i.e., each lattice site) is represented by a state vector. The state vector completely defines the state of the voxel and includes 39 entries. The 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 produce Maxwell-Boltzmann statistics for the realized equilibrium state vectors.

[0073] For processing efficiency, the voxels are grouped in 2x2x2 volumes called microblocks. The microblocks are organized to allow parallel processing of the voxels and to minimize the overhead associated with the data structure.

[0074] Figure 6 A microblock is illustrated in FIG. 3. The voxels are represented at the corners of the microblock.

[0075] Reference is made to Figure 7A and Figure 7B , the surface S( Figure 7A ) is represented in the simulation space ( Figure 7B ) as a collection of facets F α , where a is an index that enumerates a particular facet. The facets are not limited to voxel boundaries, but are generally sized on the order of or slightly smaller than the size of the voxels adjacent to the facet, such that the facet influences a relatively small number of voxels. For the purpose of implementing surface dynamics, the facets are assigned properties. In particular, each facet F α has a unit normal (n α ), a surface area (A α ), a center position (x α ), and a facet distribution function (f i (a)) that describes the surface dynamic properties of the facet. The total energy distribution function is treated in the same manner as the flow distribution for the facet and voxel interactions.

[0076] Reference is made to Figure 8 , different levels of resolution can be used in different regions of the simulation space to improve processing efficiency. Typically, the region 320 surrounding the object 322 is of the most interest and, therefore, is simulated at the highest resolution. Since the effects of viscosity diminish with distance from the object, a reduced level of resolution (i.e., an expanded voxel volume) is employed to simulate regions 324, 326 spaced at increasing distances from the object 322.

[0077] Similarly, as Figure 9As shown, a lower level of resolution can be used to simulate the region 340 around less important features of object 342, while the highest level of resolution can be used to simulate the region 344 around the most important features of object 342 (e.g., leading edge and trailing edge faces). The lowest level of resolution and the largest voxels are used to simulate the peripheral region 346.

[0078] C. Identify voxels affected by the surface patch

[0079] Refer again Figure 5 Once the simulation space has been modeled (272), voxels affected by one or more patches are identified (274). Voxels can be affected by patches in several ways. First, voxels intersected by one or more patches are affected in that they have a reduced volume relative to non-intersecting voxels. This occurs because the patches, and the material beneath the surface they represent, occupy a portion of the voxel. The fractional factor indicates the portion of the voxel not affected by the patches (i.e., the portion that can be occupied by the fluid or other material being simulated for its flow). For non-intersecting voxels, the fractional factor is equal to one.

[0080] Voxels that interact with one or more patches by transferring particles to or receiving particles from patches are also identified as voxels affected by patches. All voxels intersected by patches will include at least one state that receives particles from a patch and at least one state that transfers particles to a patch. In most cases, additional voxels will also include such states.

[0081] refer to Figure 10 For a non-zero velocity vector c i For each state i, the patch F α From the parallelepiped G iα The defined region receives particles or is directed towards the parallelepiped G. iα Defined region-transfer particle, parallelepiped G iα Having a unit normal n of the surface α and velocity vector c i The magnitude of the vector dot product is defined by the height and the surface area A of the patch. α The defined base makes the parallelepiped G iα Volume V iα It equals the base multiplied by the height.

[0082] noodle F α When the velocity vector of the state points towards the surface, from volume V iα Receive particles and transfer them to the region when the state's velocity vector points away from the facet. As will be discussed below, when another facet occupies parallelepiped G... iαWhen a surface is part of a face (as can be the case for non-convex features of a surface such as internal corners), this relationship is modified.

[0083] Patch F α Parallelepiped G iα Can overlap some or all of the voxels. The number of voxels or portions thereof depends on the size of the patch relative to the size of the voxels, the energy of the state, and the orientation of the patch relative to the lattice structure. The number of voxels affected increases with the size of the patch. Thus, as noted above, the size of the patch is typically chosen to be on the order of or smaller than the size of the voxels located near the patch.

[0084] The flux of particles of a given state moving between a voxel and a patch F α is equal to the density of particles of that state in that voxel multiplied by the volume of the region overlapping that voxel.

[0085] When a parallelepiped G iα is intersected by one or more patches, the volume V iα of the parallelepiped is equal to the sum of the volumes associated with each voxel overlapping G iα and the sum of the volumes associated with all patches intersecting G iα .

[0086] D. Performing the simulation

[0087] Referring to Figure 5 , once the voxels affected by one or more patches are identified (274), a timer is initialized to start the simulation (276). During each time increment of the simulation, the movement of particles from voxel to voxel is simulated by considering an advection phase (278-286) of the interaction of the particles with the surface patches. Next, a collision phase (288) simulates the interaction of the particles within each voxel. Thereafter, the timer is incremented (290). If the incremented timer does not indicate that the simulation is complete (292), the advection and collision phases (278-288) are repeated. If the incremented timer indicates that the simulation is complete (292), the results of the simulation are stored and / or displayed (294).

[0088] 1. Boundary conditions for surfaces

[0089] To properly model the interaction with the surface, each patch satisfies four boundary conditions. First, the combined mass of the particles received by the patch is equal to the combined mass of the particles transmitted by the patch (i.e., the net mass flux to the patch is equal to zero). Second, the combined energy of the particles received by the patch is equal to the combined energy of the particles transmitted by the patch (i.e., the net energy flux to the patch is equal to zero). These two conditions can be satisfied by requiring that the net mass flux at each energy level (i.e., energy level one and energy level two) be equal to zero.

[0090] The other two boundary conditions relate to the net momentum of the particles interacting with the patch. For a surface with no surface friction (referred to herein as a slip surface), the net tangential momentum flux is equal to zero and the net normal momentum flux is equal to the local pressure at the patch. Thus, the difference between the combined, received, and transmitted momentum of the particles with the normal n α to the patch that is perpendicular (i.e., the tangential component) is equal, while the difference between the combined, received, and transmitted momentum of the particles with the normal n α to the patch that is parallel (i.e., the normal component) is equal to the local pressure at the patch. For a non-slip surface, the friction of the surface causes the combined tangential momentum of the particles transmitted by the patch to be reduced relative to the combined tangential momentum of the particles received by the patch by a factor related to the amount of friction.

[0091] 2. Aggregation from voxel to patch

[0092] The interaction between the simulation particles and the surface is modeled by aggregating particles from voxels and providing them to patches (278). As noted above, the flux of particles between a given state voxel and a patch is related to the degree of overlap of the parallelepiped of the patch with the voxel for the given state. Only voxels whose parallelepipeds have a non-zero value are summed. As noted above, the size of the patches is chosen so that the parallelepiped volume of a patch has a non-zero value for only a few voxels. The volume of the parallelepiped and the flux of particles can have non-integer values, and thus the data processing system stores and processes these quantities as real numbers.

[0093] 3. Movement from patch to patch

[0094] Next, the particles move between patches (280). If the parallelepiped G α for the incident state of patch F iα is intersected by another patch F β , then a portion of the particles of the incident state received by patch F α will come from patch F β . In particular, patch F α will receive a portion of the particles of the given state produced by patch F β during the previous time increment. Figure 12 This relationship is illustrated in FIG. 3, where the parallelepiped Giα of the face F β of the intersecting portion 380 equals the parallelepiped G iβ of the face F α of the intersecting portion 382. For each state of the face F α of the face F α of the face F α of the face F α of the face F

[0095] For a given state, the total flux of particles into the face F α equals the sum of particles from all voxels that overlap the parallelepiped of the face and the sum of fluxes from all faces that have a volume overlap with the parallelepiped volume of the face F α . The state vector of a face (also referred to as a face distribution function) has a number of entries corresponding to the number of discrete lattice velocities in the LBM, where the number of entries is the number of discrete lattice velocities in the LBM. The input states of the face distribution function are set equal to the flux of particles into those states divided by the volume V iα .

[0096] The face distribution function is a simulation tool for generating output flux from a face and does not necessarily represent actual particles. To generate accurate output flux, values are assigned to the other states of the distribution function. The outwards states are populated using the techniques described above for populating the inwards states. In an alternative approach, the values of the outwards flux from the previous time step can be used to generate the flux from states other than the incident state.

[0097] For a parallel state (e.g., a state whose velocity is parallel to the face), the volume of the associated parallelepiped is zero. The face distribution function for a parallel state is determined as the limit of the face distribution function as the volume of the parallelepiped and any overlapping parallelepipeds approaches zero. The values for states whose velocity is zero (i.e., the rest state and states (0, 0, 0, 2) and (0, 0, 0, -2)) are initialized at the start of the simulation based on the initial conditions for temperature and pressure. These values are then adjusted over time.

[0098] 4. Perform face surface dynamics operations

[0099] Next, for each face, surface dynamics operations are performed 282 to satisfy the boundary conditions. Figure 13 A process 390 for performing surface dynamics operations for a face is illustrated in FIG. 3B.

[0100] During the gather step 278 Figure 5 ) from the voxels to the facets F α The velocities are then projected along the surface. A Boltzmann equilibrium distribution is computed (392) based on the density of the samples and the projected velocities, and the density and velocities are scaled to satisfy the constraints. A new Boltzmann equilibrium distribution is then computed (392) using the resulting density and velocities. The difference between the incident distribution and the new Boltzmann distribution is determined (394), and the combined momentum difference between all incident states and the corresponding Boltzmann distribution for the facet is computed (396).

[0101] The momentum difference is projected (398) along the tangential direction of the facet. From the momentum difference and the Boltzmann distribution, the outgoing flux of the facet is computed to satisfy the perfect slip boundary condition by satisfying zero tangential flux (399).

[0102] To account for surface friction and other factors, a surface friction component can be added to the outgoing flux distribution based on the product of the surface friction coefficient of the surface and the difference between the outgoing and incident particle fluxes. More detailed descriptions of applying surface friction and corrections to different lattice levels required for perfect mass and energy conservation are presented in the '260 patent and the '391 patent.

[0103] 5. From voxel to voxel

[0104] Referring again to Figure 5 , particles are moved between voxels (284) along the three-dimensional straight-line lattice. This voxel-to-voxel movement is the only movement operation performed for voxels that do not interact with a facet (i.e., voxels that are not located near a surface). In a typical simulation, voxels that are located too far from a surface to interact with the surface constitute the majority of the voxels.

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

[0106] For voxels that do not interact with surfaces, the move operation is computationally simple. During each time increment, the entire fill of a state is moved from its current voxel to its destination voxel. At the same time, the particles of that destination voxel are moved from that voxel to their own destination voxels. For example, a level 1 particle that is moving in the +1x and +1y directions (1,0,0) is moved from its current voxel to a voxel that is +1 in the x direction and 0 for the other directions. The particle terminates in its destination voxel with the same state that it had before it moved (1,0,0). Interactions within the voxel will have the potential to change the particle count of that state based on local interactions with other particles and surfaces. If not, the particle will continue to move along the lattice at the same rate and direction.

[0107] For voxels that interact with one or more surfaces, the move operation becomes slightly more complex. This can result in one or more fractional particles being transferred to a face. This transfer of fractional particles to a face results in fractional particles remaining in the voxel. These fractional particles are transferred to the voxel occupied by the face.

[0108] Reference Figure 11 When a portion 360 of the particles of a given state for a voxel 362 moves to a face 364 (278), the remaining portion 366 moves to the voxel 368 in which the face 364 is located, and the particles of that state are directed from the voxel 368 to the face 364. Thus, if the state fill is equal to 25, and the volume of the parallelepiped that overlaps the voxel is equal to 0.25 (i.e., one quarter of the voxel intersects the parallelepiped), then 6.25 particles will move to the face 364, and 18.75 particles will move to the voxel 368 occupied by the face 364. Because multiple faces can intersect a single voxel, the number of particles of a given state that are transferred to the voxel occupied by one or more faces is related to the fraction of the volume of the voxel that does not overlap the parallelepiped of the intersecting faces (e.g., the number of particles transferred to the voxel is equal to the total number of incident particles from the source voxel multiplied by the fraction of the voxel that does not overlap the parallelepiped of the intersecting faces).

[0109] 6. Scattering from faces to voxels

[0110] Next, the outgoing particles from each face are scattered to voxels (286). Essentially, this scattering is the inverse of the aggregation of particles moving from voxels to faces. The number of particles of a given state that move from a face to a voxel is proportional to the outgoing flux of particles from that face multiplied by the ratio of the volume of the parallelepiped that overlaps the voxel to the total volume of the parallelepiped. A scaling factor can also be included to account for the volume reduction of partial voxels. For each state, the total number of particles directed from a face to a voxel is determined based on the sum of the particles moving to that voxel from each face.

[0111] After dispersing the particles from the patches to the voxels, combining them with the particles that have been advected in from the surrounding voxels, and integerizing the result, it is possible that some of the directions in some of the voxels can either underflow (become negative) or overflow (exceed 255 in an eight-bit implementation). This will cause a gain or loss in mass, momentum, and energy after these quantities are 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. 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) chosen state of the same energy and which itself did not experience an overflow or underflow. The additional momentum resulting from this addition of mass and energy is accumulated and added to the momentum from the truncated state. By adding mass only to the same energy state, when the mass counter reaches zero, both the mass and energy are corrected. Finally, the momentum is corrected using a push / pull technique until the momentum accumulator returns to zero.

[0112] 7. Performing fluid dynamics operations

[0113] Performing fluid dynamics operations (288) Figure 5 This can be referred to as micro-dynamics or intra-voxel operations. Similarly, the advection process can be referred to as inter-voxel operations. The micro-dynamics operations described below can also be used to collide particles at the patches to produce a Boltzmann distribution.

[0114] Fluid dynamics operations are ensured in the lattice Boltzmann equation model by a particular collision operator known as the BGK collision model. This collision model mimics the dynamics of the distribution in a real fluid system. After the advection step, the conserved quantities of the fluid system are obtained from the distribution function, specifically the density, momentum, and energy. From these quantities, the equilibrium distribution function can be completely specified. The choice of the set of velocity vectors and weights and the lattice Boltzmann equation ensure that the macroscopic behavior obeys the correct fluid dynamics equations.

[0115] Variable resolution

[0116] Variable resolution (as discussed in U.S. Patent 10,360,324, which is hereby incorporated in its entirety herein) can also be employed and will use different sizes of voxels, e.g., coarse voxels and fine voxels, for different regions of the grid depending on the region of interest. For example, the region of interest can contain fine voxels in order to resolve the fluid dynamics at a smaller scale and higher resolution. Outside the region of interest, coarse voxels can be used in the grid to reduce the amount of computation without reducing the accuracy of the fluid flow simulation in the region of interest. During the simulation of the physical process, particles can be advected between regions of different resolution.

[0117] Example

[0118] Figures 14A-14F An example of pressure-driven water flow through a bipolar channel of a fuel cell is shown. Figure 14A An oblique view of the channel geometry is shown, and Figure 14B A top view of the geometry is shown. The channel height is 0.01 m. Figure 14C and Figure 14D Density and velocity contours in steady state are shown, respectively, with an inlet velocity of 0.8 m / s and flow from left to right. The simulation Mach number was set to 0.01, which is artificially increased with respect to the actual Mach number of 5.33e-4. As shown in Figure 14C , the density almost doubles through the channel. This is a clear example of artificial compressibility effects due to the artificially increased simulation Mach number. Under the constraint of mass flux conservation, the density variation leads to a velocity variation, as shown in Figure 14D . The velocity field along the flow direction is expected to be almost constant along the long channel; however, the over-compressibility effects due to the increased simulation Mach number lead to an increase in the velocity along the flow direction over the length of the channel. These artificial compressibility effects lead to less accurate flow field predictions.

[0119] Figure 14E and Figure 14F Improved simulation results for flow through a serpentine channel are shown. The improved results are achieved by applying the process 100 Figure 2B . The density variation Figure 14E and the velocity variation Figure 14F are much smaller than the density and velocity variations in Figure 14C and Figure 14D . The density Figure 14E is almost constant, as expected for incompressible fluid flow. Likewise, the velocity Figure 14F is almost constant along the flow direction, with the expected transverse velocity variation due to the no-slip boundary condition imposed on the channel walls.

[0120] Referring to Figure 2B , Figure 5 and Figure 13 , the results of the processes 100, 270 and 390 can be used to manufacture physical structures, for example, using one or more computer-controlled manufacturing systems. For example, the simulations can be used to optimize the shape of a physical object to improve fluid flow around or through the object. The optimized shape can be used as path planning data for a computer-aided manufacturing process to produce the optimized shape. Examples of shapes that can be optimized include flow paths for internal duct or channel flow (e.g., for liquid cooling systems), nozzles and diffusers for wind tunnels, and the like.

[0121] Alternatively or additionally, the results of the fluid flow simulation of process 270 can be used to predict the head loss and / or velocity profile of a fluid flowing through a pipe or channel. The results of the fluid flow simulation of process 270 can also be used to predict fluid flow through complex geometrical configurations, such as cooling systems, to predict thermal performance of a car’s cooling system; to predict channels in a GDL of a fuel cell to predict fuel cell performance; or to predict a porous medium such as carbonate rock to predict production of hydrocarbons from a subterranean reservoir. The techniques described herein improve the accuracy of such fluid flow simulations in predicting the head loss and / or velocity profile of a fluid flowing through a pipe or channel and / or predicting the thermal performance of a car’s cooling system.

[0122] Embodiments of the subject matter and the functional operations described in this specification can be implemented in digital electronic circuitry, in tangibly-embodied

[0123] The term“data processing apparatus” refers to data processing hardware and encompasses all kinds of apparatus, devices, and machines for processing data, including by way of example a programmable processor, a computer, or multiple processors or computers. The apparatus can also be, or further include, special purpose logic circuitry, e.g., an FPGA (field programmable gate array) or an ASIC (application-specific integrated circuit). The apparatus can optionally include, in addition to hardware, code that creates an execution environment for computer programs, e.g., code that constitutes processor firmware, a protocol stack, a database management system, an operating system, or a combination of one or more of them.

[0124] A computer program (which can also be referred to or described as a program, software, a software application, a module, a software module, a script, or code) can be written in any form of programming language, including compiled or interpreted languages, or declarative or procedural languages, and can be deployed in any form, including as a stand-alone program or as a module, component, subroutine, or other unit suitable for use in a computing environment. A computer program may, but need not, correspond to a file in a file system. A program can 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, in a single file dedicated to the program in question, or in multiple coordinated files, e.g., files that store one or more modules, sub programs, or portions of code). A computer program can be deployed for use by one computer or by multiple computers that are situated at one site or distributed across multiple sites and interconnected by a data communication network.

[0125] Computers suitable for the execution of a computer program can be based on general or special purpose microprocessors or any other kind of central processing unit. The computer-readable medium suitable for storing the computer program instructions and data includes all forms of non volatile memory, including by way of example semiconductor memory devices, e.g., EPROM, EEPROM, and flash memory devices; magnetic disks such as internal hard disks or removable disks; magneto optical disks; and CD-ROM and DVD-ROM disks. The processor and the memory can be supplemented by, or incorporated in, special purpose logic circuitry.

[0126] 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 that includes a middleware component, e.g., an application server, or that includes a front end component, e.g., a client computer having 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 any combination of 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 a local area network (LAN) and a wide area network (WAN), e.g., the Internet.

[0127] The computing system can include clients and servers. A client and server are generally remote from each other and typically interact through a communication network. The relationship of client and server arises by virtue of computer programs running on the respective computers and having a client-server relationship to each other. In some embodiments, a server transmits data (e.g., an HTML page) to a user device (e.g., for purposes of displaying data to and receiving user input from a user interacting with the user device), which acts as a client. Data generated at the user device (e.g., a result of the user's interaction with the user device) can be received from the user device at the server.

[0128] Particular embodiments of the subject matter have been described. Other embodiments are within the scope of the following claims. For example, the acts recited in the claims can be performed in a different order and still achieve desirable results. As one example, the processes depicted in the accompanying figures do not necessarily require the particular order shown, or sequential order, to achieve desirable results. In some instances, multitasking and parallel processing can be advantageous.

Claims

1. A computer system for numerically simulating fluid flow in a three-dimensional computer-aided design (CAD) model of a simulation space, the computer system comprising: one or more processors; and a memory, the memory comprising: a mesh preparation engine to generate and store a digital representation of the simulation space based on a digital three-dimensional CAD model, the digital representation comprising a mesh comprising a plurality of voxels; and a simulation engine to read the digital representation of the mesh in the simulation space from the mesh preparation engine, wherein the simulation engine stores instructions for numerically simulating fluid flow that, when executed by the one or more processors, cause the one or more processors to perform operations comprising: reading the digital representation of the mesh in the simulation space from the mesh preparation engine; numerically simulating a first fluid flow by applying a driving force to the plurality of voxels in the digital representation of the mesh in the simulation space to generate a digital pressure field; determining a digital volume body force to be applied to the plurality of voxels based on a pressure gradient of the digital pressure field; and numerically simulating a second fluid flow by applying the digital volume body force to the plurality of voxels.

2. The computer system of claim 1, wherein the first fluid flow and the second fluid flow comprise incompressible steady-state fluid flow.

3. The computer system of claim 1, wherein the digital volume body force comprises a digital conformal volume body force that is aligned with a flow direction.

4. The computer system of claim 1, wherein a density variation of the second fluid flow is less than a density variation of the first fluid flow.

5. The computer system of claim 1, wherein the instructions further comprise determining an actual pressure for the second fluid flow based on the driving force from the first fluid flow and a measured pressure from the second fluid flow.

6. The computer system of claim 1, wherein a number of voxels used in simulating the first fluid flow is different than a number of voxels used in simulating the second fluid flow.

7. The computer system of claim 1, wherein the instructions further comprise storing the pressure field and the volume body force in a hardware storage device for subsequent simulation of fluid flow.

8. The computer system of claim 1, wherein the pressure field comprises a force balance based on a driving force in the simulation space, fluid inertia, and friction at solid boundaries.

9. The computer system of claim 1, further comprising: storing, in the memory, a digital pressure field generated by numerically simulating the first fluid flow by applying a driving force to the plurality of voxels in the digital representation of the mesh in the simulation space; storing, in the memory, a digital volume body force to be applied to the plurality of voxels, the digital volume body force determined based on a pressure gradient of the digital pressure field; and storing, in the memory, results of a digital simulation of the second fluid flow generated by numerically simulating the second fluid flow by applying the digital volume body force to the plurality of voxels. ​ ​ 10. A method implemented by a data processing system for numerically simulating fluid flow in a three-dimensional computer-aided design (CAD) model of a simulated space, the method comprising: receiving, by the data processing system, a digital representation of the simulated space based on the digital three-dimensional CAD model, the digital representation comprising a plurality of voxels; numerically simulating, by the data processing system, a first fluid flow by applying a driving force to the plurality of voxels in the digital representation of the simulated space to generate a pressure field; determining, by the data processing system, a volumetric body force to be applied to the plurality of voxels based on a pressure gradient of the pressure field; and numerically simulating, by the data processing system, a second fluid flow by applying the volumetric body force to the plurality of voxels.

11. The method of claim 10, wherein the first fluid flow and the second fluid flow comprise incompressible, steady-state fluid flow.

12. The method of claim 10, wherein the volumetric body force comprises a conformal volumetric body force aligned with a flow direction.

13. The method of claim 10, wherein a density variation of the second fluid flow is less than a density variation of the first fluid flow. determining an actual pressure for the second fluid flow based on the driving force from the first fluid flow and a measured pressure from the second fluid flow.

14. The method of claim 10, further comprising:

15. The method of claim 10, wherein a number of voxels used in simulating the first fluid flow is different than a number of voxels used in simulating the second fluid flow. storing the pressure field and the volumetric body force in a hardware storage device for subsequent simulation of fluid flow.

16. The method of claim 10, further comprising:

17. One or more non-transitory machine-readable storage devices storing instructions for numerically simulating fluid flow in a three-dimensional computer-aided design (CAD) model of a simulated space, the instructions executable by one or more processors to cause performance of operations comprising: receiving a digital representation of the simulated space based on the digital three-dimensional CAD model, the digital representation comprising a plurality of voxels; numerically simulating a first fluid flow by applying a driving force to the plurality of voxels in the digital representation of the simulated space to generate a pressure field; determining a volumetric body force to be applied to the plurality of voxels based on a pressure gradient of the pressure field; and numerically simulating a second fluid flow by applying the volumetric body force to the plurality of voxels.

18. The one or more non-transitory machine-readable storage devices of claim 17, wherein the first fluid flow and the second fluid flow comprise incompressible, steady-state fluid flow, and wherein a density variation of the second fluid flow is less than a density variation of the first fluid flow.

19. The one or more non-transitory machine-readable storage devices of claim 17, wherein the volumetric body force comprises a conformal volumetric body force aligned with a flow direction. ​ ​ 20. The one or more non-transitory machine-readable storage devices of claim 17, wherein the instructions further comprise determining an actual pressure for the second fluid flow based on a driving force from the first fluid flow and a measured pressure from the second fluid flow.

Images