Improving performance and accuracy of stable explicit diffusion

Abstract

The present invention relates to improving the performance and accuracy of stable explicit diffusion. Methods, computer program products and systems can be used to simulate physical processes. A method includes determining an input flux to be applied to a first element. The method includes determining an applied flux, which is an amount of flux that can be applied to the first element without causing numerical instability. The method includes determining a balance flux, which is the difference between the input flux and the applied flux. The method also includes providing a balance flux to a second element.

Description

Background Art

[0001] High Reynolds number flows have been simulated by generating discretized solutions of the Navier-Stokes differential equations by performing high precision floating point arithmetic operations at each of many discrete spatial locations of variables representing macroscopic quantities (e.g., density, temperature, flow velocity). Another approach replaces the differential equations with an approach generally referred to as lattice gas (or cellular) automata, in which the macroscopic level simulation provided by solving the Navier-Stokes equations is replaced by a microscopic level model that performs operations on particles moving between sites on the lattice. Summary of the invention

[0002] In general, one innovative aspect of the subject matter described in this specification can be embodied in a method comprising an act of determining an input flux to be applied to a first element. The method comprises an act of determining an applied flux, the applied flux being an amount of flux that can be applied to the first element without causing numerical instability. The method comprises an act of determining a balance flux, the balance flux being a difference between the input flux and the applied flux. The method also comprises an act of providing the balance flux to a second element.

[0003] Other embodiments of this aspect include corresponding computer systems, devices, and computer programs recorded on one or more computer storage devices, each configured to perform the actions of the method. The system of one or more computers can be configured to perform specific operations or actions by installing software, firmware, hardware, or a combination thereof on the system, which software, firmware, hardware, or a combination thereof causes or causes the system to perform these actions in operation. One or more computer programs can be configured to perform specific operations or actions by including instructions that cause the device to perform these actions when executed by a data processing device.

[0004] The foregoing and other embodiments may each optionally include one or more of the following features, alone or in combination. A second element may be determined based on a direction of a scalar gradient. A method may include an act of providing at least a portion of a balance flux to a third element. Determining an input flux to be applied to a first element may include identifying a heat flux applied to each face of the first element. Determining an applied flux may include determining that applying a corresponding heat flux to at least one face will result in numerical instability. The physical process can be one of the following: heat flow in power generating equipment (such as engines); heat flow in turbomachinery; heat flow in electromagnetic machinery; waste heat management of electronic equipment; thermal management and protection of sensors and actuators; thermally driven stress and fatigue; thermally driven mechanical shock; thermally driven chemical changes in solids; thermally driven demagnetization; combined electrical heat generation and heat flow in conductors; heat generation and conduction in semiconductors; heat flow and current in thermoelectric devices; thermally driven dimensional changes; heat sinks; solid conduction in heat exchangers; thermal energy storage in single phase and phase change materials; detailed heat flow in complex structures such as PCBs, tires and reinforced concrete; electric heaters for engine blocks, sensors, reactors, steering wheels, car seats and batteries on automobiles; electric heaters for de-icing and defrosting car windshields and mirrors; heat conduction through vehicle structures during manufacturing and operation.

[0005] The details of one or more embodiments of the subject matter described in this specification are set forth in the accompanying drawings and the description below. Other features, aspects, and advantages of the subject matter will become apparent from the description, drawings, and claims. BRIEF DESCRIPTION OF THE DRAWINGS

[0006] Figure 1 and 2 The velocity components of the two LBM models are shown.

[0007] Figure 3 It is a flow chart of the process followed by a physical process simulation system.

[0008] Figure 4 It is a perspective view of a microscopic block.

[0009] Figure 5A and 5B yes Figure 3 Illustration of the lattice structure used in the system.

[0010] Figure 6 and 7 A variable resolution technique is shown.

[0011] Figure 8 An example of splitting the heat flux of an element into an applied flux and a balance flux is shown.

[0012] Fig. 9An example of sending balance flux sequentially across multiple elements is shown.

[0013] Fig.10 Element 1000 is shown during application of a heat flux.

[0014] Like reference numbers and names in the various drawings represent like elements. DETAILED DESCRIPTION

[0015] The present invention relates to computer simulation of physical processes, such as fluid flow.

[0016] 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 the following equation, where f i (0) is known as the equilibrium distribution function, which is defined as:

[0017]

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

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

[0020] According to this simulation, conventional fluid variables, such as mass ρ and fluid velocity u, are obtained as simple sums. Here, c i and w i The LBM model can be efficiently implemented on scalable computer platforms and runs with maximum robustness for time-invariant flows and complex boundary conditions.

[0021] A standard technique for obtaining the macroscopic equations of motion for a fluid system from the Boltzmann equation is the Chapman-Enskog method, in which successive approximations of the full Boltzmann equation are taken.

[0022] In a fluid system, small perturbations in density travel at the speed of sound. In a gas system, the speed of sound is generally determined by the temperature. The importance of the effects of compressibility in a fluid is measured by the ratio of the characteristic velocity to the speed of sound, which is called the Mach number.

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

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

[0025] More complex models may also be used, such as a 3D-2 model including 101 velocities and a 2D-2 model including 37 velocities. The velocities are more clearly described by the velocity components along each axis, as recorded in Tables 1 and 2, respectively.

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

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

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

[0029] refer to Figure 3, a physical process simulation system operates according to process 300 to simulate a physical process such as fluid flow. Prior to flow simulation, a simulation space is modeled (step 302) 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 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.

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

[0031] Re=uL / v

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

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

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

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

[0036] c i =(c ix , c iv , c iz )

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

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

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

[0040] 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}. Figure 4 Shown in.

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

[0042] S={F α}

[0043] Where α is an index that enumerates a particular facet. Facets are not restricted to the boundaries of voxels, but typically have a size that is on the order of the size of the voxels adjacent to the facet, or has a size that is slightly smaller than the size of the voxels adjacent to the facet, so that the facet affects a relatively small number of voxels. To implement surface dynamics, attributes are assigned to facets. Specifically, each facet F α With unit normal (n α ), surface area (A α )、Center position(x α ) and the surface element distribution function (fi (α)). The total energy distribution function q i (α) is treated in the same way as the flow distribution for the interaction of surface elements and volume elements.

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

[0045] Identify voxels affected by surfels

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

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

[0048] Execute simulation

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

[0050] Scalar diffusion

[0051] Numerical simulation of diffusion-dominated physical phenomena is very common due to applications in conduction heat transfer, mass diffusion, electrical conduction, etc. The governing equations for these phenomena are formulated as a set of partial differential equations (PDEs) that contain unstable diffusion and volume source terms. The numerical solution involves discretizing the spatial domain of interest and then advancing the solution in time using time integration techniques. Spatial discretization is usually done using highly automated mesh generation tools, while time discretization (time step size) needs to be carefully selected to ensure the stability and accuracy of the numerical solution at an acceptable numerical cost. Specifically, the stability characteristics of the time-marching scheme (Courant-Friedrichs-Lewy (CFL) constraint) determine the maximum time step size that can be used without making the solution unstable. Two types of time-marching schemes are commonly used - implicit and explicit. On the one hand, implicit methods satisfy the CFL constraint by construction, so large time steps can be used without making the solution unstable (however, too large time steps usually lead to inaccurate results). Implicit methods require solving large systems of matrix coefficients, making their implementation both non-trivial and computationally expensive. On the other hand, explicit methods are very simple to implement, computationally cheap (per iteration), and highly parallelizable, but require strict CFL constraints to be satisfied. This constraint on explicit diffusion schemes is given by The CFL number given is less than a certain limit (which is O(1)), where κ is the diffusion coefficient, Δ x is the size of the smallest spatial grid and Δ t is the time step size. In other words, if the spatial grid size Δ x Anywhere in the domain is reduced by a factor F, then the time step size Δ t Will have to follow F 2is reduced in order to maintain numerical stability. Therefore, for spatial grids with small-sized elements, explicit methods may require extremely small time step sizes, which seriously affects simulation performance. This is true even if the number of such small-sized elements in the simulation domain is very limited - the smallest element in the entire domain determines the CFL condition and therefore the time step size. For practical problems involving complex geometries, the use of irregular grids is inevitable for surface and volume discretization. On these grids, Δ x can vary significantly, and the use of explicit schemes can become very inefficient due to the extremely small time step sizes required by the CFL constraints. As a result, practitioners of explicit schemes spend a great deal of time and effort trying to improve the quality of the spatial mesh in an attempt to alleviate the problem. Even so, it is nearly impossible to remove all small-sized elements from any discretization of the true geometry, so a small time step size is the only way to make the solution (at least locally) stable.

[0052] Improved Diffusion for Explicit Schemes

[0053] To overcome the deficiencies of the above-mentioned explicit schemes for diffusion problems on irregular grids, new modifications are introduced to the flux computation between two neighboring elements when at least one of the two elements would otherwise violate the CFL constraint. As described below, these modifications depend on the material and geometric properties of the two elements and the state of the quantities of interest in the vicinity of the elements, and they help to stabilize the numerical solution regardless of the size of the two elements and ensure spatiotemporal accuracy. When the two neighboring elements are large (thus satisfying the CFL constraint), the flux computation reduces to a textbook implementation, which means that the described method is a consistent extension of the standard method.

[0054] An explicit Euler scheme and a finite volume formulation are assumed. In the following examples, the quantity of interest is temperature and the governing equation is the heat conduction equation, however, it should be understood that any scalar quantity can be used with an appropriate governing equation. The numerical scheme requires the calculation of the heat fluxes at all faces of the element. These fluxes are then added and used to update the temperature of the element under consideration. Consider two faces sharing neighboring elements α and β. According to Fourier's law of heat conduction, the heat flux is:

[0055]

[0056] in is the thermal conductivity of the common surface, is the temperature gradient normal to the common face, and “m” is used to specify that these quantities are evaluated at time step “m”. Since we are considering the heat entering α (rather than the heat leaving α), the negative sign in the usual form of Fourier’s law is discarded. The temperature gradient used here is calculated to ensure smoothness, 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 face 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:

[0057]

[0058] 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):

[0059]

[0060] We note that equation (3) above states that the temperature change is proportional to the net heat flux and 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.

[0061] When the time step Δ t The above form may lead to numerical instabilities when is large enough to violate the CFL constraint for at least one of the two elements. In the case where element α is smaller than element β, the CFL constraint is violated at least for element α. This numerical instability arises because, for element α (which we assume is small), the assumption used to compute The temperature gradient in the time step Δ t is incorrect. As mentioned above, for the same net energy transfer, the temperature change of small elements is larger, so the standard explicit time integration requires a reduction in the time step to ensure that the constant temperature gradient assumption is valid. Obviously, for a given Δ t , this problem persists as long as there is an instability in the problem, and it disappears only in the steady state when all incoming and outgoing fluxes on each element are perfectly balanced with each other.

[0062] As part of our new approach, we propose to replace the term defined in equation (1) with It is divided into two parts: (1) Applying flux which will be used for the temperature evolution of α (in the above sum), and (2) the balance flux It will be transferred to the other side of the interface αβ without changing the temperature of element α.

[0063] Figure 8 An example of separating the heat flux of an element into an applied flux and a balance flux is shown. Element 800 receives a heat flux, as shown by arrow 802. Conventionally, heat flux 804 would be fully applied to element 800, causing a violation of the CFL constraint and resulting in instability. In this example, heat flux 804 is separated into applied flux 806 and balance flux 810. Conceptually, applied flux 806 is the amount of heat flux that can be applied to element 800 without violating the CFL constraint. Balance flux 808 is the difference between heat flux 804 and applied flux 806.

[0064] Applied flux 806 is applied to element 800. Upon receiving heat flux from element 800, the balance flux is transferred to downstream elements (as indicated by arrows 808a-c).

[0065] In some embodiments, the balance flux may be applied to the downstream element in a time step after the time step in which the apply flux was applied to element 800. That is, if the apply flux is applied to element 800 at time t=1, the balance flux may be applied to the downstream element at some time t>1.

[0066] Heat flux term can be expressed as follows:

[0067]

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

[0069]

[0070]

[0071] and

[0072]

[0073] In the above equation, the geometric feature is represented by the distance d used to calculate the temperature gradient. αβ , element volume and and the area of ​​the common surface A αβ The material properties are represented by and (Here ρ represents density and C p represents specific heat) and thermal conductivity Consider the flux term in the summation and Reasonable estimates of the fluxes that may exist on the different faces of elements α and β are given separately.

[0074] Two flux terms and The physical explanation of is as follows.

[0075] Apply flux Represents the total flux The part of (given by equation (1)) that can be used for the temperature evolution of element α without introducing numerical instabilities. This form is derived from first principles for an isolated system consisting of elements α and β to include an estimate of the effect of the continuously evolving temperature field in the vicinity of the system. Therefore, this term shows an effect on the geometric / thermal properties of the element and the interaction of these elements with their surroundings (ΔG).

[0076] It should be noted that in equation (4), still relies on the fluxes observed on other faces at the previous time step to provide an estimate of the ongoing interactions at those faces during the current time step. For strongly transient problems, this leads to and The mismatch between This mismatch is taken into account by the second term of , which is transmitted through α to the neighboring element on the other side of the interface αβ. This balance flux is only applied to the temperature evolution when it is deposited in an element large enough to satisfy the CFL constraint (after it has been transported successively along the flux direction).

[0077] Fig. 9 An example of sequentially transmitting the balance flux across multiple elements is shown. In this example, input flux 900 is passed to a relatively small element 902 at a time step (in this example, time step t). The flux is such that applying the entire flux to element 902 would result in a violation of the CFL constraint. Therefore, a portion 912 of the flux is applied to element 902, and the balance flux 908 (i.e., the portion of the flux that cannot be applied to element 902 without violating the CFL constraint) is passed to a second element 904 (e.g., at time t+1).

[0078] At time step t+1, the simulation may attempt to apply the balance flux to the second element 904. Additional input flux 912 (e.g., originating from a different element) may also be simultaneously applied to the second element 904. This input flux 913 may be aggregated with the balance flux 908. The simulation may attempt to apply the resulting aggregated flux to the second element 904. In this example, applying the entire aggregated flux to the second element 904 would result in a violation of the CFL constraint. Therefore, a portion of the combination of input flux 915 and balance flux 908 is applied to the second element 904, and the new balance flux 910 is passed on to the larger third element 906 at time step t+2.

[0079] At time step t+2, balance flux 910 is aggregated with flux 914 (which is the flux passed to third element 906 as part of the simulation). In this example, the entire aggregate flux 916 (flux 914 and balance flux 910) can be applied to the third element without violating the CFL constraint, so no new balance flux is generated.

[0080] The above scheme strictly ensures that the correct amount of total flux is bound at the interface between any two elements α and β The flux available for the temperature evolution of the small element α is precisely controlled at the same time. Overall, the 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 As a result, the balance flux becomes equal to zero

[0081] As indicated in the above description, the process of stabilizing an explicit numerical scheme for diffusion problems on irregular spatial grids may include several algorithmic steps. Fig.10 Element 1000 is shown during application of a heat flux.

[0082] The process can identify faces to which a modified definition of heat flux is applied - any face between two elements where at least one of the elements violates the CFL condition. In some embodiments, on all other faces, the standard definition of heat flux can be used. For example, Fig.10 , in this example, heat fluxes are applied to three faces (face 1008, face 1010, and face 1012) of element 1000. In this example, the process may calculate that heat flux 1002 applied to face 1008 may use the standard definition, while heat flux 1004 applied to face 1010 and heat flux 1006 applied to face 1012 require modified heat flux definitions.

[0083] At the faces where modified heat flux is required, the process can calculate the total amount of heat flux using the spatially averaged temperature gradient near the element under consideration to ensure smoothness of the solution. Conversely, the standard heat flux can utilize the temperature gradient calculated based on the traditional difference form.

[0084] As mentioned above, the modified heat flux can be divided into two parts - the applied flux and the balance flux. The applied flux can be used in the temperature evolution equation of the element under consideration. The balance flux may or may not be used in the temperature evolution, depending on the size of the element. If the element is small enough to violate the CFL constraint, the balance flux can be transferred to the neighboring elements in the flux direction. On the other hand, if the element is large, the balance flux can be used for the temperature evolution.

[0085] As described above, the balance flux can be transported successively along the flux direction until it is finally delivered to a sufficiently large element (where it is applied to the temperature evolution).

[0086] Several approaches have been employed to overcome the problem of numerical instabilities for diffusion problems on irregular grids with varying element sizes. The most common approach is to enforce additional constraints on the grid generation tool to reduce such cases. Even so, 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) greatly increases the computational cost (even for small-sized elements that occur only once 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 abandon the explicit scheme altogether, or at least abandon the local region near the small elements, and use an implicit scheme instead. This approach suffers from the complexity of the implementation of the implicit scheme and the non-locality of the solution, which results in a system of equations that is not easily parallelizable.

[0087] In comparison, the new approach offers several distinct advantages:

[0088] (1) This approach allows the use of a single time step size that is determined based on temporal accuracy considerations (rather than the size of the smallest element in the grid). This is a huge benefit in terms of computational cost and ease of implementation for every conceivable scenario.

[0089] (2) The proposed method relies on the geometric properties of two adjacent elements, which ensures that the method works regardless of the size of the elements. Therefore, the usual constraints on the mesh generation process (mesh quality, size, etc.) can be relaxed to a large extent.

[0090] (3) The computational cost of computing the new terms is very reasonable. The mathematical form of the terms is simple and does not involve any iterations. This is in stark contrast to existing methods (reducing the time step size or using hybrid implicit-explicit schemes), which greatly increase the computational cost.

[0091] (4) Due to the volumetric nature of the formula, the scheme remains precisely preserved, which is an important requirement in many applications.

[0092] (5) The new variations are still explicit in nature and require information about elements within a small distance from the element under consideration, which means that only minor changes to the computational template in the original implementation are required. Therefore, the parallelization properties of the original explicit method are preserved and can be exploited in large-scale computations.

[0093] The new method can be used in a variety of different industrial applications. For example, the method can be used to simulate heat flow in power generation equipment (such as engines); heat flow in turbomachinery; heat flow in electromagnetic machinery (alternators, starters, traction motors, various actuators in the powertrain and in the cabin); waste heat management of electronic equipment; thermal management and protection of sensors and actuators; thermally driven stress and fatigue; thermally driven mechanical shock; thermally driven chemical changes in solids (such as plastic aging, glue curing, paint curing, heat curing); thermally driven demagnetization (such as the Curie temperature effect); combined electrical heat generation and heat flow in conductors (such as power cables); and thermal simulation of heat generation in semiconductors. Heat generation and conduction (e.g., across diodes, IGBTs, FETs); heat and current in thermoelectric devices (thermocouples, Seebeck effect); heat-driven dimensional changes; heat sinks; solid conduction in heat exchangers; thermal energy storage in single-phase and phase-change materials; detailed heat flow in complex structures such as PCBs, tires, and reinforced concrete; electric heaters for engine blocks, sensors, reactors, steering wheels, car seats, and batteries in automobiles; electric heaters for de-icing and defrosting automotive windshields and mirrors; heat conduction through vehicle structures during manufacturing and operation. Again, the above list is not exhaustive, but rather represents the kinds of applications in which the methods described herein may be used.

[0094] Embodiments of the subject matter and functional operations described in this specification may be implemented in digital electronic circuits, 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 computer program instruction modules encoded on a tangible, non-transitory program carrier for execution by a data processing device or for controlling the operation of a data processing device). Computer storage media 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.

[0095] 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 circuits (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 that constitutes processor firmware, a protocol stack, a database management system, an operating system, or a combination of one or more of them).

[0096] 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 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 stores 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 at multiple sites and interconnected by a data communication network.

[0097] The processes and logic flows described in this specification can be performed by one or more programmable computers executing one or more computer programs to perform functions by operating on input data and generating output. The processes and logic flows can also be performed by special purpose logic circuits (e.g., FPGAs (field programmable gate arrays) or ASICs (application specific integrated circuits)), and the apparatus can also be implemented as special purpose logic circuits.

[0098] A computer suitable for executing a computer program may be based on a general or special purpose microprocessor or both, or any other type of central processing unit. Typically, the central processing unit will receive instructions and data from a read-only memory or a random access memory or both. The basic elements of a computer are a central processing unit for executing or running instructions and one or more memory devices for storing instructions and data. Typically, a computer will also include or be operably coupled to receive data from or transfer data to one or more mass storage devices for storing data (e.g., magnetic disks, magneto-optical disks, or optical disks). Such devices are not required for a computer. In addition, a computer may be embedded in another device (e.g., a mobile phone, a personal digital assistant (PDA), a mobile audio or video player, a game console, a global positioning system (GPS) receiver, or a portable storage device (e.g., a universal serial bus (USB) flash drive), to name a few)).

[0099] 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 circuits.

[0100] To provide interaction with a user, embodiments of the subject matter described in this specification may be implemented on a computer having a display device (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor) for displaying information to the user and a keyboard and a pointing device (e.g., a mouse or trackball) through which the user can provide input to the computer. Other types of devices may also be used to provide interaction with a user; for example, the feedback provided to the user may be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback), and input from the user may be received in any form, including acoustic, voice, or tactile input. In addition, a computer may interact with a user by sending documents to and receiving documents from a device used by the user, for example, by sending web pages to a web browser on the user's device and responding to requests received from the web browser.

[0101] Embodiments of the subject matter described in this specification may 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 may 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).

[0102] A computing system may include a client and a server. The client and the server are usually remote from each other and usually interact through a communication network. The relationship between the client and the server is generated by means of computer programs that run on each computer and have 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 user interaction) can be received from the user device at the server.

[0103] Although this specification contains many specific implementation details, these should not be interpreted as limitations on the scope of any invention or the scope that can be claimed, but rather as descriptions of features specific to particular embodiments of particular inventions. Certain features described in this specification in the context of separate embodiments may also be implemented in combination in a single embodiment. Conversely, various features described in the context of a single embodiment may also be implemented in multiple embodiments individually or in any suitable sub-combination. In addition, although the above features may be described as working in certain combinations and even initially stated as such, one or more features from the claimed combination may be cut out of the combination in some cases, and the claimed combination may be directed to a sub-combination or a variation of the sub-combination.

[0104] Similarly, although operations are depicted in a particular order in the accompanying drawings, this should not be understood as requiring that the operations be performed in the particular order shown or in sequence, or that all of the operations shown be performed, to achieve the desired results. In some cases, multitasking and parallel processing may be advantageous. In addition, the separation of various system modules and components in the above-described embodiments should not be understood as requiring such separation in all embodiments, and it should be understood that the described program components and systems can generally be integrated into a single software product or packaged into multiple software products.

[0105] 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 some cases, multitasking and parallel processing may be advantageous.

Claims

1. A computer-implemented method for simulating a physical process, comprising the steps of: determining an input heat flux to be applied to a first element of a spatial grid, including identifying a corresponding heat flux to be applied to each of corresponding faces of the first element; determining whether applying the corresponding heat flux to at least one of the faces of the first element will result in numerical instability under a Courant-Friedrichs-Lewy constraint; In response to determining that applying the corresponding heat flux to at least one of the faces of the first element will result in numerical instability under the Courant-Friedrichs-Lewy constraint, splitting the corresponding input heat flux into an applied heat flux and a balance heat flux, the applied heat flux being an amount of heat flux that can be applied to the temperature evolution of the first element without causing the numerical instability under the Courant-Friedrichs-Lewy constraint; determining the balance heat flux, the balance heat flux being the difference between the corresponding input heat flux and the applied heat flux; as well as providing the balance heat flux to a second element of the spatial grid, The applied heat flux is determined based on the geometric characteristics, material properties and interaction of the first element with its surroundings.

2. The method according to claim 1, wherein: The second element is determined based on the direction of the scalar gradient.

3. The method of claim 1, further comprising providing at least a portion of the balance heat flux to a third element of the spatial grid.

4. The method according to claim 1, wherein: The physical process is one of: heat flow in power generating equipment such as engines; heat flow in turbomachinery; heat flow in electromagnetic machinery; waste heat management from electronic equipment; thermal management and protection of sensors and actuators; thermally driven stress and fatigue; Thermally driven mechanical shock; Thermally driven chemical changes in solids; Thermally driven demagnetization; Combined electrical heat generation and heat flow in conductors; Heat generation and conduction in semiconductors; heat flow and current in thermoelectric devices; Thermally driven dimensional changes; heat sinks; solid conduction in heat exchangers; thermal energy storage in single-phase and phase change materials; detailed heat flow in complex structures such as PCBs, tires, and reinforced concrete; electric heaters for engine blocks, sensors, reactors, steering wheels, car seats, and batteries in cars; electric heaters for de-icing and defrosting car windshields and mirrors; and heat conduction through vehicle structures during manufacturing and operation.

5. A non-transitory computer-readable medium storing instructions that, when executed, cause a computer processor to perform operations for simulating a physical process, the operations comprising the steps of: determining an input heat flux to be applied to a first element of a spatial grid, including identifying a corresponding heat flux to be applied to each of corresponding faces of the first element; determining whether applying the corresponding heat flux to at least one of the faces of the first element will result in numerical instability under a Courant-Friedrichs-Lewy constraint; In response to determining that applying the corresponding heat flux to at least one of the faces of the first element will result in numerical instability under the Courant-Friedrichs-Lewy constraint, splitting the corresponding input heat flux into an applied heat flux and a balance heat flux, the applied heat flux being an amount of heat flux that can be applied to the temperature evolution of the first element without causing the numerical instability under the Courant-Friedrichs-Lewy constraint; determining the balance heat flux, the balance heat flux being the difference between the corresponding input heat flux and the applied heat flux; as well as providing the balance heat flux to a second element of the spatial grid, The applied heat flux is determined based on the geometric characteristics, material properties and interaction of the first element with its surroundings.

6. The non-transitory computer readable medium of claim 5, wherein: The second element is determined based on the direction of the scalar gradient. 7 . The non-transitory computer readable medium of claim 5 , further comprising providing at least a portion of the balance heat flux to a third element of the spatial grid.

8. The non-transitory computer readable medium of claim 5, wherein: The physical process is one of: heat flow in power generating equipment such as engines; heat flow in turbomachinery; heat flow in electromagnetic machinery; waste heat management from electronic equipment; thermal management and protection of sensors and actuators; thermally driven stress and fatigue; Thermally driven mechanical shock; Thermally driven chemical changes in solids; Thermally driven demagnetization; Combined electrical heat generation and heat flow in conductors; Heat generation and conduction in semiconductors; heat flow and current in thermoelectric devices; Thermally driven dimensional changes; heat sinks; solid conduction in heat exchangers; thermal energy storage in single-phase and phase change materials; detailed heat flow in complex structures such as PCBs, tires, and reinforced concrete; electric heaters for engine blocks, sensors, reactors, steering wheels, car seats, and batteries in cars; electric heaters for de-icing and defrosting car windshields and mirrors; and heat conduction through vehicle structures during manufacturing and operation.

9. A system for simulating a physical process, comprising one or more processing devices and one or more hardware storage devices storing instructions, the instructions being operable when executed by the one or more processing devices to cause the one or more processing devices to perform operations comprising: determining an input heat flux to be applied to a first element of a spatial grid, including identifying a corresponding heat flux to be applied to each of corresponding faces of the first element; determining whether applying the corresponding heat flux to at least one of the faces of the first element will result in numerical instability under a Courant-Friedrichs-Lewy constraint; In response to determining that applying the corresponding heat flux to at least one of the faces of the first element will result in numerical instability under the Courant-Friedrichs-Lewy constraint, splitting the corresponding input heat flux into an applied heat flux and a balance heat flux, the applied heat flux being an amount of heat flux that can be applied to the temperature evolution of the first element without causing the numerical instability under the Courant-Friedrichs-Lewy constraint; determining the balance heat flux, the balance heat flux being the difference between the corresponding input heat flux and the applied heat flux; as well as providing the balance heat flux to a second element of the spatial grid, The applied heat flux is determined based on the geometric characteristics, material properties and interaction of the first element with its surroundings.

10. The system according to claim 9, wherein: The second element is determined based on the direction of the scalar gradient.

11. The system of claim 9, further comprising providing at least a portion of the balance heat flux to a third element of the spatial grid.

12. The system according to claim 9, wherein: The physical process is one of: heat flow in power generating equipment such as engines; heat flow in turbomachinery; heat flow in electromagnetic machinery; waste heat management from electronic equipment; thermal management and protection of sensors and actuators; thermally driven stress and fatigue; Thermally driven mechanical shock; Thermally driven chemical changes in solids; Thermally driven demagnetization; Combined electrical heat generation and heat flow in conductors; Heat generation and conduction in semiconductors; heat flow and current in thermoelectric devices; Thermally driven dimensional changes; heat sinks; solid conduction in heat exchangers; thermal energy storage in single-phase and phase change materials; detailed heat flow in complex structures such as PCBs, tires, and reinforced concrete; electric heaters for engine blocks, sensors, reactors, steering wheels, car seats, and batteries in cars; electric heaters for de-icing and defrosting car windshields and mirrors; and heat conduction through vehicle structures during manufacturing and operation.

Images