General wall boundary condition treatment for k-omega turbulence model

Abstract

The present disclosure relates to a general wall boundary condition treatment for k-Omega turbulence models. Disclosed is a technique for simulating a physical process and determining a boundary condition for a specific energy dissipation rate of a k-Omega turbulence fluid model of a fluid flow by calculating a value of a specific energy dissipation rate of turbulence valid for a viscous sublayer, a buffer layer and a logarithmic region of a defined boundary in a simulation space from a cell center distance and a fluid flow variable. The value is determined by applying a buffer layer correction factor as a first boundary condition for the specific energy dissipation rate and by applying a viscous sublayer correction factor as a second boundary condition for the specific energy dissipation rate.

Description

Technical Field

[0001] This manual relates to computer simulations of physical processes such as physical fluid flow. Background Technology

[0002] Discrete solutions to the Navier-Stokes differential equations are generated by performing high-precision floating-point arithmetic operations on variables representing macroscopic physical quantities (e.g., density, temperature, flow velocity) at each of many discrete spatial locations, thereby simulating high Reynolds number flows. One approach is to simulate the problem of interest using the so-called SST k-Omega turbulence model, a widely used turbulence model that solves two partial differential equations (k) and (Omega). One of these equations (k) is for turbulent kinetic energy (TKE), and one of these equations (Omega) is for specific energy dissipation rate (SEDR). The combination of these equations defines the turbulent state of the fluid flow.

[0003] In principle, it is possible to solve the governing equations of fluid dynamics (e.g., the Navier-Stokes equations) to simulate the problem of interest. Unfortunately, most practical engineering problems involve flow conditions occurring under turbulent conditions at high Reynolds number flows. This necessitates the application of turbulence modeling techniques in computational fluid dynamics (CFD) simulation solvers to reduce the time and computational model size to realistic levels, enabling solutions to be obtained in such conditions at reasonable computational and time costs. Turbulence modeling is a method for approximating the behavior of turbulent fluid flows.

[0004] One such model is the shear stress transport formula of the k-Omega turbulence model, or simply the "SST k-Omega turbulence model". The SST k-Omega turbulence model involves the solution of two partial differential equations. Typically, for partial differential equations, the k and Omega partial differential equations require specification of boundary conditions in order to solve these equations.

[0005] Although the equation (k) for turbulent kinetic energy is governed by Dirichlet and Neumann boundary conditions that bind the equation to the specific energy dissipation rate (Omega partial differential equation), the Omega partial differential equation lacks Dirichlet or Neumann type boundary conditions at the wall boundary, and thus, it is not easy to combine the Omega partial differential equation with the k partial differential equation to provide a solution.

[0006] To provide efficient solutions to these equations, near-wall analysis is typically performed to simplify the Omega partial differential equations, thereby deriving ancillary equations that describe the behavior of the Omega partial differential equations as they asymptotically approach the wall boundary.

[0007] However, while the auxiliary equation is not a formal boundary condition that can be used to solve the Omega partial differential equation, forcing it near the wall boundary allows for a solution. The problem is how to enforce the auxiliary equation. Various solutions have been proposed. One proposed solution for enforcing the behavior of the auxiliary equation generates additional mesh points to recover the asymptotic behavior. Another proposed solution modifies the auxiliary equation to enforce it as a Dirichlet boundary condition. Summary of the Invention

[0008] As mentioned above, the drawback of the specific SST k-Omega turbulence model and the more general k-Omega model is that the Omega partial differential equation for the specific energy dissipation rate lacks Dirichlet or Neumann-type boundary conditions at the wall boundary, and therefore cannot be easily combined with the k equation to provide a solution.

[0009] This paper describes a general wall boundary condition treatment for Omega partial differential equations, used to specify boundary conditions in simulations using the SST k-Omega turbulence model, and more generally, the solution is generally applicable to k-Omega turbulence model formulas.

[0010] The general wall boundary conditions described in this paper automatically apply the correct auxiliary asymptotic relations required in the viscous-sublayer and logarithmic layer, thus eliminating the pathological mesh-dependent accuracy degradation exhibited when the viscous sublayer is not properly resolved by the mesh. The described general wall boundary conditions can eliminate the pathological mesh-dependent accuracy degradation exhibited when the first element leaving the wall boundary is located within a buffer layer. Furthermore, if the remainder of the boundary layer is properly resolved, the general wall boundary condition processing can accurately predict results, regardless of the position of the first element leaving the wall.

[0011] According to one aspect, a computer-implemented method for simulating fluid flow around a simulated physical object includes: receiving a model of a simulation space by one or more computing systems, the model comprising a grid defining a representation of physical objects in the simulation space, wherein the grid comprises a plurality of cells having a resolution describing the surfaces of the physical objects; determining boundary conditions for the specific energy dissipation rate of a k-Omega turbulent fluid flow model of the fluid flow by determining the cell centers of the cells in the grid and by the one or more computing systems calculating, based on the cell center distances and fluid flow variables, values ​​of the specific energy dissipation rate of turbulence effective for viscous layers, buffer layers, and logarithmic regions defined in the simulation space.

[0012] Other areas include computing systems and computer program products.

[0013] One or more of the above aspects may include any one or more of the following features or other features as disclosed below.

[0014] For cells located at the boundary position y+<3, where y is the cell at the boundary, some aspects apply a buffer layer correction factor as a boundary condition for the specific energy dissipation rate. Determining the boundary conditions involves one or more computational systems applying a buffer layer correction factor as a boundary condition for the specific energy dissipation rate of cells located at the boundary position y+<3, where y is the cell at the boundary, and the correction factor is given according to the following equation:

[0015] ω′ Hyb =f blend ω Hyb

[0016] Where, ω′ Hyb It is the correction factor, f blend It is a blending function, and ω Hyb It is the viscous layer correction function.

[0017] Determining the boundary conditions includes applying a viscous subsurface correction factor as the boundary condition for the specific energy dissipation rate. Determining the boundary conditions also includes having one or more computational systems apply a viscous subsurface correction factor as the boundary condition for the specific energy dissipation rate, wherein the viscous subsurface correction factor is given by the following equation:

[0018]

[0019] in, It is a correction factor. It is the cell at position 1, and It is the cell at position 2.

[0020] Some aspects also include: accessing the k-Omega model; initializing the accessed k-Omega model with the determined boundary conditions; and performing the initialized k-Omega model to simulate fluid flow around the simulated physical object. Some aspects also include: accessing a k-Omega turbulent fluid flow model, wherein the k-Omega turbulent fluid flow model includes a first partial differential equation for determining the turbulent kinetic energy of the fluid flow and a second partial differential equation for determining the specific energy dissipation rate of the fluid flow in the simulation space.

[0021] Some aspects also include determining whether the location is at the buffer layer; and when it is at the buffer layer, applying a correction that increases the specific energy dissipation rate only at the buffer layer. Applying the correction also includes applying a mixing function that acts on the specific energy dissipation rate at the buffer layer and prevents the mixing function from affecting the value at the viscous layer, which is defined as a boundary in the simulation space.

[0022] One or more of the above aspects can provide one or more of the following advantages.

[0023] Compared to current boundary condition strategies, this method for specifying boundary conditions can significantly improve the accuracy of turbulent CFD simulations using any family of k-Omega models. Unlike current boundary condition strategies, this method removes grid dependencies found on grids without a properly analytical viscous underlayer from the predictions. Unlike current boundary condition strategies, this method removes or significantly reduces grid dependencies on grids where the first grid element leaving the wall resides in the buffer layer. Unlike current boundary condition strategies, this method provides boundary condition treatment throughout the inner layers of the turbulent boundary layer, independent of the location of the first grid point leaving the wall. Unlike current boundary condition strategies, this method can generalize the application of boundary conditions for Omega partial differential equations at the wall, thereby enforcing the asymptotic behavior of both the viscous and logarithmic underlayers. Attached Figure Description

[0024] Figure 1 A system for simulating fluid flow is described, which includes turbulent boundary layer models for compressible or incompressible flow.

[0025] Figure 2 A flowchart illustrating the operations used for the formulation of the k-Omega SST turbulence model is depicted.

[0026] Figure 3 A flowchart illustrating the operations used for correction to determine boundary conditions is depicted.

[0027] Figure 4 The graph illustrates the solution of the Omega PDE in turbulent channel flow.

[0028] Figure 5 It is a graph illustrating the combined auxiliary relationships of Omega at the viscous and logarithmic layers.

[0029] Figure 6 The graph illustrates the mesh correlation based on the resolution of the viscous sublayer.

[0030] Figure 7It is a graph illustrating the distribution of Omega and channel flow when the near-wall resolution is not fine enough in the gradation.

[0031] Figure 8 This is a diagram illustrating two cells near the boundary.

[0032] Figure 9 The graph (uncorrected) shows the average velocity distribution and friction factor for different near-wall resolutions, with multiple examples of resolution cases plotted, but for clarity, only selected examples are indicated by reference lines.

[0033] Figure 10 The graph (corrected) shows the average velocity distribution and friction factor for different near-wall resolutions, with several examples of resolution cases plotted, but for clarity, only selected examples are indicated by reference lines.

[0034] Figure 11 The graph illustrates the mismatch of Omega in the buffer layer.

[0035] Figure 12 This is a graph illustrating the performance of boundary layer conditions in the buffer layer, with several examples plotted for different resolutions, but for clarity, only selected examples are indicated by reference lines.

[0036] Figure 13 The graph illustrates the mixing function that improves accuracy in the buffer layer.

[0037] Figure 14 The graph illustrates the correction of Omega in the buffer layer.

[0038] Figure 15 The graph illustrates the results of applying boundary condition correction, showing several examples of resolution scenarios, but for clarity, only selected examples are indicated by reference lines. Detailed Implementation

[0039] The following describes a general auxiliary equation that is derived and correctly reproduces near-wall behavior and wall function constraints using a correction that eliminates mesh-dependent solutions in near-wall constraints, providing a practical implementation of auxiliary boundary condition treatment for the Omega partial differential equations in the k-Omega turbulence model used in computational fluid dynamics simulations.

[0040] In the following discussion, boundary condition processing is introduced for the Omega partial differential equation applicable to all k-Omega turbulence models. In the following discussion, the distance from the wall boundary is referred to as "y+n", where "y" is the boundary and "n" is the number of cells away from the boundary, and n is a number that takes an integer or a decimal value.

[0041] The boundary condition method discussed below refers to: the viscous sublayer nominally between y+0.1 and y+5; the log layer away from the wall as a fully turbulent region, for example, y+>25; and the buffer layer that is neither the viscous sublayer nor the fully turbulent log region and is generally defined by 5<y+<25. Other variations within these ranges are possible.

[0042] The boundary condition method is based on deriving a general auxiliary relation for ω, which describes the asymptotic behavior of the Omega partial differential equation in the viscous sublayer (e.g., the layer adjacent to the wall boundary, typically y+ ≤ 5) ω v and fully conforms to the turbulent equilibrium log layer (e.g., the layer away from the wall, typically y+ ≥ 25, but still within the boundary layer) ω L . This auxiliary equation method also includes a correction factor that remedies the accuracy degradation associated with relatively poor grid resolution in the viscous sublayer.

[0043] Referring now to Figure 1 , system 10 includes a simulation engine 34 configured for CFD simulation using any model in the k-Omega turbulence model family. The simulation engine 34 includes a turbulence model module 34a that executes a specific k-Omega turbulence model and a boundary module 34b that automatically specifies boundary conditions for use with the turbulence model module 34a. The boundary conditions executed in the boundary module 34b are derived from a general solution for ω that satisfies turbulent behavior within the boundary layer. In other words, an auxiliary relation is provided that describes the asymptotic behavior of the Omega partial differential equation in the viscous sublayer (e.g., the layer adjacent to the wall boundary, for example, y+ ≤ 5) while fully conforming to the turbulent equilibrium log layer (e.g., the layer away from the wall, for example, y+ ≥ 25, while still within the boundary layer).

[0044] 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 massively parallel computing system 12 (stand-alone or cloud-based) and a client system 14. The server system 12 includes a memory 18, a bus system 22, an interface 20 (e.g., user interface / network interface / display or monitor interface, etc.) and a processing device 24. Stored in the memory 18 are a mesh preparation engine 32 and the simulation engine 34.

[0045] although Figure 1 A mesh preparation engine 32 in memory 18 is shown, but the mesh preparation engine could be a third-party application running on a different system than server 12. Regardless of whether the mesh preparation engine 32 runs in memory 18 or on a different system, it receives a user-supplied mesh definition 30, prepares a mesh based on the physical object being modeled, and sends (and / or stores) the prepared mesh to simulation engine 34 for simulation. System 10 accesses a data repository 38 storing 2D and / or 3D meshes (Cartesian and / or curves), coordinate systems, and libraries. Any number of physical objects or physical fluid flows can be represented in the simulation space used to represent physical objects or physical flows. The simulation can be used for a wide range of technical and engineering problems, including aerodynamic / aerospace analysis, weather, environmental engineering, industrial systems, biological systems, general fluid flow, and combustion systems.

[0046] Now refer to Figure 2 The diagram illustrates a process 40 for simulating fluid flow around a representation of a physical object. In the example discussed herein, the physical object is an airfoil. However, the use of an airfoil is merely illustrative, as physical objects can have any shape and, in particular, can have planar and / or curved (one or more) surfaces. Process 40 receives, for example, a mesh (or grid) for the physical object being simulated from either client system 14 or a data repository 38. In other embodiments, an external system or server 12 generates the mesh for the physical object being simulated based on user input.

[0047] The process pre-calculates 44 geometric quantities from the retrieved mesh. Process 40 determines 46 boundary conditions (e.g., in boundary module 34b) based on the retrieved mesh and flow field variables. Flow field variables such as velocity, omega, k, temperature, etc., are used to calculate ω of the object being simulated. HCorr This is used in boundary module 34b. Once the turbulence variables are known, the partial differential equations k and Omega are solved. The turbulence model provides eddy viscosity, which represents the effects of turbulence that the mesh cannot resolve and is used to advance the solution. The mesh includes the simulation space or divides the simulation space into multiple cells.

[0048] The process uses the selected k-Omega turbulence model to perform a computational fluid dynamics (CFD) simulation. The model is initialized with boundary conditions and executed using pre-computed geometric quantities corresponding to the retrieved mesh in the turbulence model module 34a.

[0049] Now refer to Figure 3The diagram illustrates a process 60 for simulating fluid flow around a simulated physical object. Process 60 includes receiving a model of a simulation space 62 from a computing system. This model includes a grid defining a representation of the physical object within the simulation space, wherein the grid comprises multiple cells with a resolution describing the surfaces of the physical object. Process 60 further includes determining boundary conditions for the specific energy dissipation rate of a k-Omega turbulent fluid flow model of the fluid flow by determining the cell centers of the cells in the grid and by one or more computing systems calculating, based on the cell center distances and fluid flow variables, 66 values ​​of the specific energy dissipation rate of turbulence effective for viscous layers, buffer layers, and logarithmic regions defined at the boundaries of the simulation space.

[0050] Process 60 further includes applying the 68a buffer layer correction factor as a boundary condition for the specific energy dissipation rate for cells located at the boundary position y+<3 (where y is the cell at the boundary). Process 60 also includes applying the 69b viscous sublayer correction factor as a boundary condition for the specific energy dissipation rate.

[0051] To derive the wall function assumptions discussed above, we analyze the fluid flow in a turbulent channel. The Reynolds-averaged governing equations for an incompressible turbulent channel are given by the following equation:

[0052]

[0053]

[0054] For a fully expanded flow and no spanwise velocity, the assumptions are as follows:

[0055]

[0056] The fact that the x and z components of turbulence are zero means that, on average, the flow does not change in the x and z directions. However, the flow can change instantaneously in parts of the geometric space in these directions.

[0057] Therefore, the continuity equation becomes

[0058]

[0059] This means that the advection is zero, therefore the steady-state governing equation is:

[0060]

[0061] The evaluation equation in the Z direction yields

[0062]

[0063] The equation is in equilibrium along the entire Z-axis, which means that the gradient of the Reynolds stress is constant and symmetrical along the Z-axis.

[0064]

[0065] This suggests that the Reynolds stress is antisymmetric.

[0066]

[0067] Therefore, the Reynolds stress at z = 0 is zero.

[0068]

[0069] This means that the pressure does not change along the Z-axis and is sufficiently sufficient to “measure” the pressure at the center z = 0 in the spanwise direction. The evaluation equations in the Y-direction show that the wall normal Reynolds stress is only a function of the wall normal direction.

[0070]

[0071] This can be integrated to obtain...

[0072]

[0073] According to this equation, the pressure gradient in the streamwise direction is equal to the pressure gradient at the wall.

[0074]

[0075] Evaluate the X-momentum equation using equation (5.c) and apply the symmetry conditions to obtain...

[0076]

[0077]

[0078] Evaluating the equation at the center of the channel, y = H, yields...

[0079]

[0080] Therefore, the pressure drop in the channel can be calculated as

[0081]

[0082] The momentum equation can also be written based on the total stress as follows:

[0083]

[0084] Using equation (7.a), equation (7.c) can be written as follows:

[0085]

[0086] Therefore, the total stress in the channel can be written as a function of y+ (Equation 7d) and the friction Reynolds number (Equation 8.0).

[0087]

[0088] If this equation is expanded and the so-called "law-of-the-wall" is used to calculate the velocity gradient, Equation (8.0) can be written as:

[0089]

[0090] Wherein, the wall transfer law is a general function that correlates the value of the average velocity distribution with y+, and the velocity distribution consists of the following three regions: the viscous sublayer (wherein, for y+<5, the velocity follows a linear function), the logarithmic layer (wherein, for y+>25, the velocity follows a logarithmic function), and the buffer layer (for 5<y+<25, acting as a smooth transition function between the two layers), and at the wall it becomes:

[0091]

[0092] Therefore, the Reynolds stress can be calculated accurately

[0093]

[0094] The flow starts to become turbulent at Re τ >100, therefore y + <5, u + =y + can be written for the viscous sublayer as:

[0095]

[0096] This shows that in the viscous sublayer, the Reynolds stress is approximately zero. In the logarithmic region where y + >25, the Reynolds stress is approximated as

[0097]

[0098] Therefore, for a sufficiently large Re τ , the Reynolds stress is equal to the wall shear stress. In addition, even for a Reynolds stress Re τ =100, in the logarithmic region where y + =25, Equation (10.b) also provides (almost no turbulence), while for larger Reynolds numbers it is approximately equal to the turbulence at the wall or boundary.

[0099] By assuming the proximity is very close to the wall and eliminating the second term (for high Re), τ (It should be zero), and the method of generating the near wall is implemented by applying formula (9.c).

[0100]

[0101]

[0102] This could potentially correlate wall shear stress with turbulent variables. Turbulent kinetic energy budgeting predicts an equilibrium between generation and dissipation within the logarithmic region.

[0103]

[0104] Using the k-ε model, the eddy current viscosity is written as:

[0105]

[0106] Using equation (12.0) in (11.0), the wall shear can be written based on the turbulent kinetic energy in the logarithmic layer.

[0107]

[0108] Therefore, friction speed It can be written as

[0109] u τ =C μ 1 / 4 k 1 / 2 (14.a)

[0110] To calculate the value of ε in the viscous sublayer and the logarithmic region, the following relationship can be used. The diffusion of TKE balances the specific energy dissipation rate according to the TKE budget at the wall.

[0111]

[0112] This can be integrated as

[0113]

[0114] Given the boundary conditions k(y=0)=0 and C2=0, the turbulent kinetic energy can be written as:

[0115]

[0116] Since the turbulent kinetic energy must always be positive, this requires

[0117]

[0118] Further information can be obtained by finding the location of the critical point, since the point exists and is the minimum value according to equation (14.b).

[0119]

[0120] Since the critical point must lie within the domain y≥0, this imposes C1≤0. Applying these two conditions yields...

[0121]

[0122] Since this holds true for all y including y = 0, the condition in equation (19.0) becomes...

[0123] C1≥0, &C1≤0 (20.0)

[0124] The only way for equation 20.0 to be true is when C1 = 0, which implies that

[0125]

[0126] The relationship of dissipation rates in the logarithmic region can be obtained by substituting equation (10.b) into equation (11.0).

[0127]

[0128] Furthermore, by using equation (14.a), the logarithmic relation can be written as

[0129]

[0130] If we use Equation (23.0) to evaluate the eddy viscosity, we can find that the eddy viscosity is linear in the logarithmic region.

[0131] ν T =u τ κy (24.0)

[0132] Near-wall turbulence analysis of planar channel flow

[0133] The results predicted by the SST k-Omega model exhibit near-wall correlation. This correlation is due to the near-wall behavior of the Omega(ω) equation in the region near the wall.

[0134] Using the original model:

[0135]

[0136]

[0137] For a steady flow very close to the wall, equation (1.b) can be simplified to

[0138]

[0139] The equation has the following solution, which applies only to solutions very close to the wall.

[0140]

[0141] To develop the wall function formula, a turbulent channel flow analysis previously developed for the specific energy dissipation rate equation (22.0) was used. However, before using it, the process linked the ε and ω variables by using the definition of the corresponding eddy viscosity.

[0142]

[0143] Using this relationship, we can write equation (22.0) based on ω.

[0144]

[0145] Figure 4 The numerical solution for turbulent channel flow is shown to predict the two constraints presented in equations (27) and (29). The goal is to develop equations that incorporate these two constraints.

[0146]

[0147] Figure 5 It is shown that an analytical hybrid function (e.g., Equation 30.0) can predict numerical computations very closely. Figure 5 The viscous underlying layer y+<5(ω) shows the correct auxiliary relation that produces the correct relation as a function of y+. V ) and logarithmic layer y+>25(ω L The mixture of auxiliary relationships of ω at ) hyb The impact of ).

[0148] If the grid resolution is not in y + <0.1, for example in Figure 7 In this study, when the near-wall resolution is not fine enough in the hierarchical classification, numerical calculations cannot accurately predict the stiffness variation of ω in the viscous sublayer. The numerical calculations illustrated with circles used a near-wall resolution of y+ = 0.05. The results show that the value of ω at the first cell does not fall within the expected theoretical value (ω...). v This leads to an incorrect distribution of ω at the remaining grid points located in the inner layer (see...). Figure 7 ), where the value of ω (solid circle) is higher than ω v Distribution. Therefore, the results show the grid dependence caused by the inability of diffusion discretization to properly resolve the steep gradient of ω as shown below.

[0149] For the viscous sublayer where y+ < 2, Equation 26 needs to be accurately discretized. Solving Equation 26 requires finding the layer at a distance of Y1 from the wall. Figure 8 Equation 30 is enforced at the first cell of the array. However, diffusion discretization does not adequately recover the analytical gradient of ω at the wall.

[0150]

[0151] It is assumed that specifying the value of ω at the first cell exiting the wall will be sufficient to capture the diffusion flux of the first element requiring diffusion flux. However, the first cell exiting the wall is not bound, and these values ​​are specified using Equation 30. The finite volume discretization used for diffusion cannot reproduce the steep change of ω, as shown below.

[0152]

[0153] Here, by replacing the viscous underlying behavior of ω and a certain algebraic manipulation in Equation 32.b, it is shown that the gradient of ω at the element surface is incorrect:

[0154]

[0155] That is, numerical discretization using a specified value of ω will always overpredict the diffusion flux, which is why Figure 9 The slope is steeper in the middle.

[0156] However, if the grid resolution is y + If <0.1, then the first grid point is close enough to the wall to allow the solution to eventually asymptotically approach the correct distribution of ω, even if incorrect flux is introduced by Equation 33.

[0157] Therefore, the global solution remains unchanged by any further mesh refinement. However, if the first grid point is located where y+>0.1, Equation 33 introduces a large numerical error that prevents the solution from asymptotically approaching the correct distribution of ω before it changes toward its logarithmic region value. Figure 6 As shown, this will generate mesh-dependent values ​​in the calculation and when the near-wall resolution is found to be 0.1 <y + The error is more pronounced when the value is less than 2.

[0158] It is possible to analytically derive the diffusion equation for ω, and thus a correction for the diffusion flux (Equation 34) can be introduced to account for the error in Equation 33.

[0159]

[0160] In the TKE equation, the corrected value of ω is used, thereby preventing the over production of associated TKE when ω is damped. It is very important to note that Equation 34 is only valid for y + <3.5 where Equation 33 holds. To show the accuracy improvement obtained by using Equation 34, numerical calculations of channel flow using Equation 30 are presented below, which show the predicted velocity distribution and friction factor.

[0161] Figure 9 shows the typical conventional results obtained when using near-wall resolution with y+<3, indicating that the results are grid-dependent. The lack of accuracy in friction factor and velocity distribution is obvious. That is, for grids with near-wall resolutions of Y+ 0.5 and 2.5, large errors are introduced (see Figure 9 ). Multiple resolutions are plotted, but for clarity, only selected resolutions are indicated by reference lines.

[0162] However, Figure 10 shows that when Equation 34 is used, grid dependency is eliminated and accuracy is greatly improved, as illustrated. The accuracy improvement of the friction factor and velocity distribution by Equation (34.0) is quite obvious. That is, the large errors introduced for grids with near-wall resolutions of Y+ 0.5 and 2.5, shown in Figure 9 , are eliminated for these resolutions, and accurate results are obtained, which confirms the importance of Equation 34 for improving calculation accuracy. Multiple resolutions are plotted, but for clarity, only selected resolutions are indicated by reference lines.

[0163] Although the error associated with the viscous discretization of Equation (33.0) can be corrected by Equation 34.0, unfortunately, when the near-wall resolution is located in the buffer layer, another error is encountered, where the value of ω satisfies neither the viscous sublayer nor the logarithmic layer.

[0164] Figure 11 shows that Equation 30 can predict the exact solution of the partial differential equation of ω (Equation 25.b) very closely in both the internal and logarithmic regions, while for Equation 30, the numerical value at the buffer layer is very close to the numerical solution of Equation 25.b, the values in the buffer layer are under-predicted.

[0165] Figure 12 shows the performance of the new boundary condition when the first grid point away from the wall is located on the buffer layer (ω hyb ). Because in the buffer layer (5<y+<25), ω v and ω LNeither of these are effective behaviors, so this behavior can be expected. The velocity distribution indicates underestimation of the logarithmic layer, which is associated with the increase in friction velocity relative to expectations. As a result, wall shear is greater, thus increasing the friction factor. Figure 10 As shown, the increase in the friction factor is due to the lower value of ω at the first cell, as predicted by Equation 34, when the first cell is located in the buffer layer. The lower value of ω results in a larger value of TKE, which in turn leads to greater wall shear. Multiple resolutions are plotted, but for clarity, only the selected resolution is indicated by the reference line.

[0166] In some cases, the viscous sublayer ω V With logarithmic layer ω L Both values ​​will underpredict the values ​​needed at the buffer layer, so neither of them, nor any of their averages, will be sufficient.

[0167] Therefore, it is necessary to develop a buffer layer correction that increases the value of ω only at the buffer layer. To provide such a correction function, equations 27, 29, and 30 will be used as follows:

[0168]

[0169] Equation 35 is drawn in Figure 13 Above, it shows an increase of about 30% in the buffer layer, and it also exhibits asymptotic behavior extending beyond the buffer layer, which may be problematic. A new ω boundary condition equation ω is valid for viscous layers, buffer layers, and logarithmic regions. hyb It can be written as:

[0170] ω' Hyb =f blend ω Hyb (36.0)

[0171] Figure 14 The improvement at the buffer layer is shown when Equation 35 is used to correct Equation 30. The improvement, as illustrated in the diagrams for the buffer layer and at the buffer layer, is significant. Here, the behavior of ω inside the inner layer is general for wall function theory, and therefore Re is independent in the case of turbulent flow. These relations are developed under the assumption of equilibrium, and thus, can only be expected to be valid for simulated flows that satisfy the fundamental assumptions that give rise to the derivation of these equations.

[0172] Figure 14 The diagram illustrates the ω-relationship, which is related to ω. hyb Apply a mixing function to produce a value that is accurately called ω throughout the entire wall (all y+, the viscous sublayer, the buffer layer, and the logarithmic layer). HCorrThe function. While it can be said that the results of the mixing function will not apply when the flow is not in equilibrium, these are valid when applied to the entire wall function method. However, in general, these assumptions provide near-wall modeling for the idealized flow applied to derive the boundary condition equations, detaching their conditions. Therefore, Equation 36, which recovers the viscous sublayer, the logarithmic layer, and provides corrections for the mesh-independent results at the viscous sublayer and the buffer layer, can be regarded as the optimal boundary condition treatment for ω, and is valid as long as the wall function method is valid.

[0173] Equation 35 exhibits asymptotic behavior that penetrates well into the viscous sublayer and the logarithmic region. Since the correction to Equation 30 only needs to be applied in the buffer layer, the mixing function can be restricted to operate in the buffer layer to prevent it from affecting the values ​​of unnecessarily corrected values ​​at the inner and logarithmic layers, as in Equation 36.1.

[0174]

[0175] The results of using viscous underlayer correction (Equation 34) and buffer layer correction (Equation 36.1) significantly improved the accuracy at near-wall behavior, such as Figure 15 As shown in the image.

[0176] Figure 15 The diagram illustrates the application of the general boundary condition ω. HCorr The results show a significant improvement over the results obtained with the buffer layer. Multiple resolutions were plotted, but for clarity, only the selected resolution is indicated by the reference line. While the results with the buffer layer are greatly improved, the predicted friction factor also matches better and essentially eliminates the [previous issues]. Figure 11 The large error observed in the study.

[0177] Embodiments of the subject matter and functional operation described in this specification may be implemented as digital electronic circuit systems, tangibly implemented computer software or firmware, computer hardware (including the structures disclosed in this specification and their equivalents), or combinations thereof. Embodiments of the subject matter described in this specification may be implemented as one or more computer programs (i.e., one or more modules of computer program instructions encoded on a tangible, non-transitory program carrier for execution by a data processing apparatus or for controlling the operation of a data processing apparatus). Computer storage media may be machine-readable storage devices, machine-readable storage substrates, random or serial access memory devices, or combinations thereof.

[0178] The term "data processing apparatus" refers to data processing hardware and encompasses all kinds of devices, apparatuses, and machines for processing data, including (for example) programmable processors, computers, or multiple processors or computers. The apparatus may also be or further include special-purpose logic circuit systems (e.g., FPGAs (Field-Programmable Gate Arrays) or ASICs (Application-Specific Integrated Circuits)). In addition to hardware, the apparatus may optionally include code that creates the execution environment for computer programs (e.g., code constituting processor firmware, protocol stacks, database management systems, operating systems, or combinations thereof).

[0179] Computer programs (also referred to or described as programs, software, software applications, modules, software modules, scripts, 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 standalone programs or as modules, components, subroutines, or other units suitable for use in a computing environment. A computer program may, but does not need to, correspond to a file in a file system. A program may be stored as a portion of a file that holds other programs or data (e.g., in a markup language document, in a single file dedicated to the program in question, or in one or more scripts within multiple coordinating files (e.g., files storing portions of one or more modules, subroutines, or code). Computer programs can be deployed such that they execute on one computer or on multiple computers located at one site or distributed across multiple sites and interconnected by a data communication network.

[0180] The processes and logic flows described in this specification can be executed by one or more programmable computers, which execute one or more computer programs to perform functions by manipulating input data and generating outputs. The processes and logic flows can also be executed by a special-purpose logic circuit system (e.g., a FPGA (Field-Programmable Gate Array) or an ASIC (Application-Specific Integrated Circuit)), and the apparatus can also be implemented as a special-purpose logic circuit system (e.g., a FPGA or an ASIC).

[0181] A computer suitable for executing computer programs can be based on a general-purpose or special-purpose microprocessor, or both, or any other type of central processing unit. Typically, the central processing unit receives instructions and data from read-only memory or random access memory, or both. The basic components of a computer are the central processing unit for executing or fulfilling instructions and one or more memory devices for storing instructions and data. Typically, a computer will also include one or more mass storage devices (e.g., magnetic, magneto-optical, or optical discs) for storing data, or will be operatively coupled to receive data from or transfer data to such mass storage devices, or both; however, a computer does not need to have such devices. Furthermore, a computer can be embedded in another device (e.g., to name just a few, mobile phones, personal digital assistants (PDAs), mobile audio or video players, game consoles, GPS receivers, or portable storage devices (e.g., Universal Serial Bus (USB) flash drives)).

[0182] 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, and CD-ROM and DVD-ROM disks. Processors and memory may be supplemented by or contained within dedicated logic circuitry systems.

[0183] To provide interaction with the user, embodiments of the subject matter described in this specification can 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 pointing device (e.g., a mouse or trackball), through which the user provides input to the computer. Other types of devices can also be used to provide interaction with the user; for example, feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback), and input from the user can be received in any form including acoustic, speech, or tactile input. Additionally, the computer can interact with the user by sending and receiving documents to and from the device used by the user (e.g., by sending a webpage to a webpage in response to a request received from a webpage on the user's device).

[0184] Embodiments of the subject matter described in this specification may be implemented in a computing system that includes back-end components (e.g., as a data server), middleware components (e.g., an application server), or front-end components (e.g., a client computer having a graphical user interface or 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 components, middleware components, or front-end components. The components of the system may be interconnected via 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).

[0185] A computing system may include clients and servers. Clients and servers are often geographically separated and typically interact via a communication network. The client-server relationship arises from computer programs running on their respective computers and having a client-server relationship with each other. In some embodiments, the server sends data (e.g., HTML pages) to a user device acting as a client (e.g., for the purpose of displaying data to a user interacting with the user device and receiving user input from the user). Data generated at the user device (e.g., the result of user interaction) may be received at the server from the user device.

[0186] While this specification contains numerous details of specific implementations, these should not be construed as limiting the scope of any invention or the scope that can be claimed, but rather as descriptions of features that may be specific to particular embodiments of a particular invention. 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 separately or in any suitable sub-combination in multiple embodiments. Furthermore, while features may be described above as functioning in certain combinations and even initially claimed in this way, in some cases, one or more features from the claimed combination may be removed from the combination, and the claimed combination may involve sub-combinations or variations thereof.

[0187] Similarly, although the operations are depicted in a specific order in the accompanying drawings, this should not be construed as requiring these operations to be performed in the specific order shown or in sequential order, or to perform all the illustrated operations to achieve the desired result. In some cases, multitasking and parallel processing can be advantageous. Furthermore, the separation of the various system modules and components in the above embodiments should not be construed as requiring such separation in all embodiments, and it should be understood that the described program components and systems can generally be integrated together in a single software product or packaged into multiple software products.

[0188] 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 result. As an example, the processes depicted in the drawings do not necessarily require the specific order or sequence shown to achieve the desired result. In some cases, multitasking and parallel processing can be advantageous.

Claims

1. A computer-implemented method for simulating fluid flow around a simulated physical object, the method comprising: A model of a simulation space is received by one or more computing systems. The model of the simulation space includes a mesh that defines a representation of a physical object being simulated in the simulation space. The mesh includes a plurality of cells having a resolution that describes the surface of the physical object, and the mesh defines a boundary having a viscous layer, a buffer layer, and a logarithmic region of the boundary. For the first element leaving the boundary, the boundary conditions for the specific energy dissipation rate of the k-Omega turbulent fluid flow model are determined by the following operation: Determine the cell center of each cell in the grid; The one or more computing systems calculate the specific energy dissipation rate of the effective turbulence for the viscous layer, the buffer layer, and the logarithmic region of the boundary defined in the simulation space, based on the determined cell center and fluid flow variables, to substantially eliminate the grid dependency on the grid where the first grid element leaving the boundary resides in the buffer layer.

2. The method according to claim 1, wherein, For a cell located at position y+<3 on the boundary, where y is the cell at the boundary and y+ is the distance from the boundary, the method further includes: The buffer layer correction factor is applied by the one or more computing systems as a boundary condition for the specific energy dissipation rate.

3. The method according to claim 1, wherein, Determining the boundary conditions also includes: The one or more computing systems apply a buffer layer correction factor as a boundary condition for the specific energy dissipation rate of the cell located at the boundary position y+<3, where y is the cell at the boundary and y+ is the distance from the boundary, and the correction factor is given by the following formula: in, It is the buffer layer correction factor. It is a mixture function, and It is the viscous layer correction function.

4. The method according to claim 1, further comprising: The one or more computing systems apply a viscous bottom-level correction factor as a boundary condition for the specific energy dissipation rate.

5. The method according to claim 1, wherein, Determining the boundary conditions also includes: The viscous sub-base correction factor is applied by the one or more computing systems as a boundary condition for the specific energy dissipation rate, wherein the viscous sub-base correction factor is given by the following formula: in, It is a correction factor for viscous substrate. y 1 It is the cell at position 1, and y 2 It is the cell at position 2. ω ν It is the specific energy dissipation rate in the viscous sublayer. ω v d yes ω ν The discretized form, and ω v f express ω v d The correction form.

6. The method according to claim 3, wherein, Determining the boundary conditions also includes: The viscous sub-base correction factor is applied by the one or more computing systems as a boundary condition for the specific energy dissipation rate, wherein the viscous sub-base correction factor is given by the following formula: in, It is a correction factor for viscous substrate. y 1 It is the cell at position 1, and y 2 It is the cell at position 2. ω ν It is the specific energy dissipation rate in the viscous sublayer. ω v d yes ω ν The discretized form, and ω v f express ω v d The correction form.

7. The method according to claim 1, further comprising: Access the k-Omega model; Initialize the visited k-Omega model with the determined boundary conditions; as well as Execute the initial k-Omega model to simulate fluid flow around the simulated physical object.

8. The method according to claim 1, further comprising: Access the k-Omega turbulent fluid flow model, wherein the k-Omega turbulent fluid flow model includes The first partial differential equation used to determine the turbulent kinetic energy of the fluid flow; as well as The second partial differential equation is used to determine the specific energy dissipation rate of the fluid flow in the simulated space.

9. The method according to claim 1, further comprising: Determine if the location is at the buffer layer; and if it is at the buffer layer... The application is a correction that increases the specific energy dissipation rate only at the buffer layer.

10. The method according to claim 9, wherein, The application of the correction also includes: A mixing function is applied that acts on the specific energy dissipation rate at the buffer layer and prevents the mixing function from affecting the value at the viscous layer of the boundary defined in the simulation space.

11. A system for simulating physical process flows around a simulated physical object, the system comprising: One or more processor devices; A memory operatively coupled to the one or more processor devices; A storage medium storing a computer program, the computer program including instructions for causing the system to perform the following operations: A model of a simulation space is received, the model of the simulation space comprising a mesh defining a representation of a physical object being simulated in the simulation space, wherein the mesh comprises a plurality of cells having a resolution describing the surface of the physical object, and the mesh defines a boundary having a viscous layer, a buffer layer and a logarithmic region of the boundary. For the first element leaving the boundary, the boundary condition for the specific energy dissipation rate of the k-Omega turbulent fluid flow model is determined by instructions for causing the system to perform the following operations: Determine the cell center of each cell in the grid; as well as The specific energy dissipation rate of the effective turbulence for the viscous layer, the buffer layer, and the logarithmic region defined by the boundary in the simulation space is calculated based on the determined cell center and fluid flow variables, so as to substantially eliminate the grid dependence of the first grid element residing on the grid at the buffer layer.

12. The system of claim 11, further configured to: The cell is located at the boundary at a position y + < 3, where, y is the cell at the boundary, and y+ is the distance from the boundary; and A buffer layer correction factor is applied as a boundary condition for the specific energy dissipation rate, wherein the correction factor is given by the following formula: in, It is the buffer layer correction factor. It is a mixture function, and It is the viscous layer correction function.

13. The system according to claim 11, further configured to: The viscous subsurface correction factor is applied as the boundary condition for the specific energy dissipation rate, wherein, The viscous sublayer correction factor is given by the following formula: in, It is a correction factor for viscous substrate. y 1 It is the cell at position 1, and y 2 It is the cell at position 2. ω ν It is the specific energy dissipation rate in the viscous sublayer. ω v d yes ω ν The discretized form, and ω v f express ω v d The correction form.

14. The system of claim 12, further configured to: The viscous subsurface correction factor is applied as the boundary condition for the specific energy dissipation rate, wherein, The viscous sublayer correction factor is given by the following formula: in, It is a correction factor for viscous substrate. y 1 It is the cell at position 1, and y 2 It is the cell at position 2. ω ν It is the specific energy dissipation rate in the viscous sublayer. ω v d yes ω ν The discretized form, and ω v f express ω v d The correction form.

15. The system according to claim 12, further configured to: Access the k-Omega model; Initialize the visited k-Omega model with the determined boundary conditions; and Execute the initial k-Omega model to simulate fluid flow around the simulated physical object.

16. A computer program product for simulating a physical process, the computer program product being tangibly stored on a non-transitory computer-readable storage medium, the computer program product including instructions for causing a system to perform the following operations: A model of a simulation space is received, the model of the simulation space comprising a mesh defining a representation of a physical object being simulated in the simulation space, wherein the mesh comprises a plurality of cells having a resolution describing the surface of the physical object, and the mesh defines a boundary having a viscous layer, a buffer layer and a logarithmic region of the boundary. For the first element leaving the boundary, the boundary condition for the specific energy dissipation rate of the k-Omega turbulent fluid flow model is determined by instructions for causing the system to perform the following operations: Determine the cell center of the cells in the grid; and The specific energy dissipation rate of the effective turbulence for the viscous layer, the buffer layer, and the logarithmic region defined by the boundary in the simulation space is calculated based on the determined cell center and fluid flow variables, so as to substantially eliminate the grid dependence of the first grid element residing on the grid at the buffer layer.

17. The computer program product of claim 16, further comprising instructions for causing the system to perform the following operations: The cell is located at the boundary at a position y + < 3, where, y is the cell at the boundary, and y+ is the distance from the boundary; and A buffer layer correction factor is applied as a boundary condition for the specific energy dissipation rate, wherein the correction factor is given by the following formula: in, It is the buffer layer correction factor. It is a mixture function, and It is the viscous layer correction function.

18. The computer program product of claim 16, further comprising instructions for causing the system to perform the following operations: The viscous subsurface correction factor is applied as the boundary condition for the specific energy dissipation rate, wherein, The viscous sublayer correction factor is given by the following formula: in, It is a correction factor for viscous substrate. y 1 It is the cell at position 1, and y 2 It is the cell at position 2. ω ν It is the specific energy dissipation rate in the viscous sublayer. ω v d yes ω ν The discretized form, and ω v f express ω v d The correction form.

19. The computer program product of claim 17, further comprising instructions for causing the system to perform the following operations: The viscous subsurface correction factor is applied as the boundary condition for the specific energy dissipation rate, wherein, The viscous sublayer correction factor is given by the following formula: in, It is a correction factor for viscous substrate. y 1 It is the cell at position 1, and y 2 It is the cell at position 2. ω ν It is the specific energy dissipation rate in the viscous sublayer. ω v d yes ω ν The discretized form, and ω v f express ω v d The correction form.

20. The computer program product of claim 16, further comprising instructions for causing the system to perform the following operations: Access the k-Omega model; Initialize the visited k-Omega model with the determined boundary conditions; and Execute the initial k-Omega model to simulate fluid flow around the simulated physical object.

Images