Turbulent flow analysis system and program

Abstract

To provide a turbulent flow analysis system and program.SOLUTION: The turbulent flow analysis system 100 of the present invention for analyzing the Navier-Stokes equations includes: definition means 101 for defining initial conditions for the pressure distribution and velocity distribution of an analysis model; first calculation means 102 for determining the velocity distribution using a forward method; second calculation means 103 for calculating an unsteady Reynolds stress term from the calculated velocity distribution; reduction means 104 for reducing the Reynolds number in a range including strong turbulence that causes divergence by using a control variable obtained by vectorizing the unsteady Reynolds stress term; and third calculation means 105 for converging the pressure distribution including an inertia term calculated from the velocity distribution, an additional term for converging the velocity outflow, and a gradient of the unsteady Reynolds stress term required to satisfy the law of conservation of mass. The analysis system iteratively calculates the pressure distribution until a predetermined criterion is reached.SELECTED DRAWING: Figure 1

Description

[Technical field]

[0001] The present invention relates to a turbulent flow analysis technique, and more particularly to a system and a program for numerically analyzing highly turbulent flows with large Reynolds numbers. [Background technology]

[0002] The Navier-Stokes equations, which describe the motion of fluids, are written as complex partial differential equations, and their analysis requires a great deal of effort. Therefore, CFD (Computational Fluid Dynamics) technology has been developed to approximately predict the behavior of fluids through numerical simulations. In particular, when analyzing the flow of incompressible fluids using the finite difference method, various methods have been developed to strictly establish the continuity equation.

[0003] For example, Japanese Patent Laid-Open Publication No. 3-229156 (Patent Document 1) discloses a technique for finding a numerical solution to the Navier-Stokes equations and the continuity equations using a simultaneous discretization method that maintains the same block length and a simultaneous linear equation with this as its coefficient, using a conjugate gradient system solution with incomplete triangular decomposition.

[0004] However, in the conventional techniques such as those disclosed in Patent Document 1, the solutions tend to diverge or convergence takes a long time, which is insufficient from the viewpoint of practical use. Therefore, there has been a demand for a further technique for stably converging the solutions of the Navier-Stokes equations.

[0005] In relation to this inconvenience, the present inventor proposed in Patent Specification No. 7299613 (Patent Document 2) a system, method and program in which even if a new velocity outflow occurs, the value is converged to 0 after a few steps, and the continuity equation is approximately satisfied, thereby easily obtaining a convergent solution using the simplest central difference method. [Prior art documents] [Patent documents]

[0006] [Patent Document 1] Japanese Patent Application Publication No. 3-229156 [Patent Document 2] Patent No. 7299613 specification Summary of the Invention [Problem to be solved by the invention]

[0007] Although Patent Document 2 can provide efficient convergence in regions where the Reynolds number is small and the turbulence in the flow is relatively weak, the solution is prone to divergence in flows with large Reynolds numbers and high turbulence, and further technology is needed to stably converge the solution to the Navier-Stokes equations. [Means for solving the problem]

[0008] According to this embodiment, 1. A system for analyzing the Navier-Stokes equations, comprising: A means for defining initial conditions for pressure and velocity distributions in the model of the analysis; a first calculation means for calculating a velocity distribution using a forward method; A second calculation means for calculating an unsteady Reynolds stress term from the calculated velocity distribution; a reduction means for reducing the Reynolds number in a range including strong turbulence that causes divergence, using a control variable obtained by vectoring the unsteady Reynolds stress term; A third calculation means for converging the pressure distribution including an inertia term calculated from the velocity distribution, an additional term for converging the velocity outflow, and a gradient of the unsteady Reynolds stress term; Including, A turbulent flow analysis system is provided which is characterized by iteratively calculating the pressure distribution until a predetermined criterion is reached.

[0009] Furthermore, according to this embodiment, a program for performing turbulent flow analysis is provided. Effect of the Invention

[0010] According to the present invention, it is possible to provide a turbulent flow analysis system and program capable of stably converging a strong flow with a high degree of turbulence and a large Reynolds number. [Brief description of the drawings]

[0011] [Figure 1] FIG. 1 is a functional block diagram of an information processing device 100 according to the present embodiment. [Diagram 2] FIG. 2 is a flowchart of a computer-executable program implementing the above algorithm. [Diagram 3] FIG. 3 shows the result 300 of a calculation using the turbulent flow calculation program of Patent Document 2. [Figure 4] FIG. 4 shows the calculation results 400 of the average flow velocity at each cross section using a conventional method, with the vertical axis representing the average velocity and the horizontal axis representing the relative position in the main flow direction. [Diagram 5] FIG. 5 shows the result 500 of calculating the outlet mean velocity using the circular pipe turbulent flow model explained in FIG. 3 using the program of this embodiment, with Ω as a parameter. [Figure 6] FIG. 6 shows a result 600 obtained when the calculation algorithm of this embodiment is applied to the circular pipe turbulent flow model described with reference to FIG. 3 when Ω=2. [Figure 7] FIG. 7 shows a calculation result 700 corresponding to FIG. 4, obtained by taking the average velocity on the vertical axis and the relative position in the main flow direction on the horizontal axis, for the velocity calculated by this embodiment. [Figure 8] Figure 8 shows the calculation results 800, with the vertical axis representing the ratio of the calculated pipe friction coefficient to the pipe friction coefficient of a smooth circular pipe, and the horizontal axis representing the Reynolds number. [Figure 9] Figure 9 shows the results of a simulation in which random numbers between -0.4 and 0.4 were added to the velocity components in the three directions. [Figure 10] Figure 10(a) and (b) show the simulated velocity distribution around a sphere and a cylinder, respectively. [Figure 11] FIG. 11 shows the calculation results of the velocity obtained by the same embodiment as in FIG. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0012] The process described in Japanese Patent No. 7299613 will be described below in order to clarify the entire process of the present application. The Navier-Stokes equation of motion is written as the following equations (1) to (3).

[0013]

number

[0014]

number

[0015]

number

[0016] The above formula is made dimensionless with the characteristic length L and the characteristic velocity U as shown in the following formula (4). In formula (4), p represents the pressure distribution, ρ represents the fluid density, Re represents the Reynolds number, and μ represents the fluid viscosity coefficient. The program module that implements the above formulas (1) to (3) corresponds to the first calculation means of this embodiment.

[0017]

number

[0018] When the above formulas (1) to (3) are calculated using the formula (5), the following formula (6) is obtained.

[0019]

number

[0020]

number

[0021] In this embodiment, an initial condition is given, and the velocity field for the next step is obtained by an explicit method.

[0022] In the above formula (6), ∇u represents the velocity source. In order to keep the time change of the velocity source zero, that is, to maintain the equation of continuity in the next step, the following formula (7) must be satisfied.

[0023]

number

[0024] In order to satisfy the above formula (7), the pressure distribution needs to satisfy the following formula (8).

[0025]

number

[0026] Here, even if ∇u is not zero, in order for ∇u to converge quickly to zero, an additional term is added to the above formula (8) to give the following formula (9). That is, by adding an additional term obtained by multiplying ∇u, which is the velocity outflow, by a positive constant κ to the right-hand side as in formula (9), the value of ∇u at each location can be kept small, and the total amount of velocity outflow can be converged to zero. Note that the additional term is not limited to ∇u multiplied by a constant κ in formula (9), and various functions of ∇u can be used as the additional term, as described later.

[0027]

number

[0028] Here, by substituting the above formula (9) into formula (6), we obtain the following formula (10).

[0029]

number

[0030] As shown in the above formula (10), even if a velocity outflow exists, that is, when ∇u occurs, the additional term makes it disappear immediately, and a solution that approximately satisfies the continuity equation can be obtained. Note that hereafter, the left side of formula (10) is referred to as the additional term.

[0031] The above-mentioned solution is executed by an information processing device such as a computer or various computing systems. Here, the case where the solution is executed by an information processing device will be described as an example. Fig. 1 is a functional block diagram of an information processing device 100 according to this embodiment.

[0032] The information processing device 100 of this embodiment includes an initial condition definition unit 101, a velocity distribution calculation unit 102, and an unsteady Reynolds stress calculation unit 103.

[0033] The initial condition definition unit 101 is a means for defining the initial conditions of the pressure distribution and velocity distribution of the model to be analyzed. The initial condition definition unit 101 can also define various conditions such as boundary conditions, and corresponds to the definition means in this embodiment.

[0034] The velocity distribution calculation unit 102 is means for calculating the velocity distribution at the next time based on the velocity distribution at a certain time and the pressure distribution calculated by the pressure distribution calculation unit 103 .

[0035] The unsteady Reynolds stress calculation unit 103 provides a function of calculating each stress for describing the pressure distribution in terms of unsteady Reynolds stress, and generating a control variable from the unsteady Reynolds stress term in a specified cell.

[0036] The information processing device 100 also includes a Reynolds number reducing unit 104 and a pressure distribution calculation unit 105. The Reynolds number reducing unit 104 reduces the Reynolds number of a cell in a range including strong turbulence that causes divergence based on the above control variables, and essentially suppresses divergence by increasing viscosity.

[0037] The pressure distribution calculation unit 105 is a means for calculating the pressure distribution using an inertia term calculated using the velocity distribution at a certain time, an additional term for converging the velocity outflow, and the gradient of the unsteady Reynolds stress term. The pressure distribution calculation unit 105 calculates the pressure distribution using an additional term for appropriately converging the solution and the gradient of the unsteady Reynolds number provided to satisfy the law of conservation of mass.

[0038] Detailed processing in the velocity distribution calculation unit 102, the unsteady Reynolds stress calculation unit 103, the Reynolds number reduction unit 104 and the pressure distribution calculation unit 105 will be described later.

[0039] The above-mentioned software blocks correspond to functional means realized by the CPU executing the program of the present embodiment and causing each piece of hardware to function. The functional means shown in each embodiment may be realized entirely in software, or part or all of them may be implemented as hardware that provides the same functions.

[0040] Furthermore, all of the above-mentioned functional means do not necessarily have to be included in the configuration shown in Fig. 1. For example, in another preferred embodiment, each functional means may be realized by cooperation of a plurality of information processing devices 100, for example, distributed processing, cloud computing, or the like, and a part or all of the software may be configured by a hardware circuit such as an FPGA.

[0041] The algorithm of the turbulent flow analysis process executed by the turbulent flow analysis system of this embodiment will be described below. In this embodiment, the time average value of the velocity average is defined as (11),

[0042]

number

[0043]

number

[0044] Then, the velocity component including the fluctuation is defined by the following equation (13).

[0045]

number

[0046] Here, the convection term included in the Navier-Stokes equation of motion is decomposed into a steady component and a fluctuating component, and the unsteady Reynolds stress term is described by equation (14).

[0047]

number

[0048] By explicitly including the unsteady Reynolds stress in equation (14), the convection term can be formulated as follows in equation (15):

[0049]

number

[0050] The program module that performs the above formulation corresponds to the second calculation means. It is as described in fluid mechanics textbooks, but here, the unsteady Reynolds stress term included in the inertia term of the Navier-Stokes equation of motion is extracted as a parameter, the unsteady Reynolds stress term is regarded as a vector, and its absolute value is used as the control variable of equation (16).

[0051]

number

[0052] In this embodiment, the values ​​of the control variables are constantly monitored for all cells, and the Reynolds number of a cell whose control variable exceeds a predetermined value and its surrounding cells is forcibly lowered, i.e., the viscosity of the cell is forcibly increased, thereby improving convergence. In other words, the technical idea of ​​preventing divergence by increasing the viscosity of a cell that is likely to diverge is introduced. Note that a cell is the smallest volume element used when calculating pressure distribution and velocity distribution.

[0053] An example of a control using a control variable is e.g. (1) Reduce the Reynolds number of cells with a control variable value above 6 to 7 (hard braking), (2) Reduce the Reynolds number of cells where the control variable value exceeds 4 to 1000 (deceleration brake) However, these values ​​are merely examples, and it will be understood that functions with other values ​​or control variables as independent variables can be appropriately set in order to optimize a specific embodiment. The program module that performs this processing corresponds to the reduction means for reducing the Reynolds number.

[0054] By the way, when a value α that cannot be approximated is added to the right side of equation (10), the following equation (17) is obtained.

[0055]

number

[0056] When this constant α is added, the conventional method cannot converge to the following equation (18).

number

[0057] Therefore, in order to converge the above formula (18), the most direct way is to subtract the value corresponding to α added to the additional term in formula (10) to obtain formula (19).

number

[0058] The next problem is how to set the value of α. After investigation, the inventors found that the unsteady Reynolds stress term steadily generates a negative value of the following equation (20).

[0059]

number

[0060] In the process of deriving equation (6) by the calculation of equation (5), the Reynolds stress term is included in the form of its gradient. Therefore, in equations (8), (9), and (10), the Reynolds stress term is included in the form of its gradient. In other words, it has been found that the error in equation (10) arising from the unsteady Reynolds stress term can be offset by subtracting a constant multiple of the gradient of the unsteady Reynolds stress term from the additional term on the right side of equation (10). With this constant as Ω, in this embodiment, the above equation (10) is modified as shown in the following equation (21) using the gradient of the unsteady Reynolds stress term.

[0061]

number

[0062] In order to realize the above-mentioned modified formula (21), it is necessary to include the gradient of the Reynolds stress term in formula (9) to obtain formula (22). That is, formula (21) can be realized by substituting formula (22) into formula (6). Note that a program module that performs iterative calculations using the calculation of formula (22) below corresponds to the third calculation means of this embodiment. Formula (22) is essentially formula (9) with a term associated with Ω added, but formula (9) is simplified using the continuity equation, whereas formula (22) is the complete form before simplification.

[0063]

number

[0064] That is, in this embodiment, the velocity distribution is first obtained to calculate the unsteady Reynolds stress term, the unsteady Reynolds stress term is determined from the velocity distribution, the control variable is obtained, the Reynolds number of the cell where this value exceeds a certain value is reduced, and then the pressure distribution is calculated. By using the above algorithm, even if the flow has a high Reynolds number and a high turbulence degree, the Navier-Stokes equation and the law of conservation of mass are satisfied, and a scalable simulation is possible that does not diverge due to turbulence. It has been confirmed that the value of Ω is appropriate for the circular pipe turbulence model and the external flow around the bluff body, but if problems occur in other flow fields, the value of Ω can be appropriately determined so that the law of conservation of mass is satisfied.

[0065] FIG. 2 shows a flowchart of a computer executable program for implementing the above algorithm. The following process is performed for all cells set in the calculation. The process shown in FIG. 2 starts from step S200, and in step S201, the initial velocity distribution and pressure distribution are set. In step S202, the velocity distribution is obtained using the forward method using equations (1) to (3).

[0066] Next, in step S203, the unsteady Reynolds stress term in each cell is found, and in step S204, the absolute value of the unsteady Reynolds stress term, i.e., the control variable, is calculated for all cells. It is determined whether the control variable of a cell is 4 or more or 6 or more, and if the control variable is 4 or more, the Reynolds number of that cell and neighboring cells is reduced to 1000, and if the control variable is 6 or more, the Reynolds number of that cell and neighboring cells is reduced to 7, and in step S206, the pressure distribution is found from equation (12). Note that in this embodiment, the neighborhood refers to the length of one side of a cell being δ,

[0067]

number

[0068] The above example is just one example of a condition set to prevent the solution from diverging even if there is a large disturbance. In general, it is possible to give the solution as a function with the cell where the control variable is maximized as the center and the distance from that cell as variables.

[0069] After step S205, in step S206, equation (22) is repeatedly used to iteratively calculate the pressure distribution until a predetermined criterion is reached, and when the predetermined criterion is reached, the process proceeds to step S202 to find the velocity distribution for the next time step.

[0070]

number

[0071]

number

[0072] After the pressure distribution of all cells is calculated in S206, it is judged in S207 whether the calculation of the speed and cloth pressure for the required set time has been completed. If not, the process returns to S202 and moves to the calculation of the next time step. If the calculation has been completed, the process ends in S208.

[0073] In this iterative process, the processing of this embodiment prevents a decrease in calculation efficiency due to divergence of the solution, and makes it possible to simulate a fluid system using the Navier-Stokes equations even in the case of a highly turbulent flow with a large Reynolds number.

[0074] The present invention will be described below with reference to specific embodiments. (1) Conventional algorithm Ω is a steadily occurring component of equation (10).

[0075]

number

[0076] Figure 4 shows the calculation result 400 of the average flow velocity at each cross section by the conventional method, where the vertical axis represents the average velocity and the horizontal axis represents the relative position in the main flow direction. As shown in Figure 4, in the conventional method, the average flow velocity decreases between the inlet and the outlet, which indicates that the law of conservation of mass is not satisfied. This means that a negative error occurs in the velocity outflow term, and the velocity outflow cannot be converged to 0.

[0077] (2) Algorithm of this embodiment FIG. 5 shows the result 500 of calculating the outlet mean velocity with Ω as a parameter, which is the result of calculating the turbulent pipe model described in FIG. 3 using the program of this embodiment. The result of FIG. 5 shows that when this embodiment is applied to the turbulent pipe model, the law of conservation of mass can be satisfactorily satisfied by setting Ω=2 and subtracting the gradient of the unsteady Reynolds stress term from the velocity source equation (10). We have confirmed that Ω=2 is appropriate for pipe flows with various Reynolds numbers and external flows around bluff bodies, but if problems arise, the value of Ω can be determined in a similar manner.

[0078] Fig. 6 shows a result 600 obtained by applying the calculation algorithm of this embodiment to the circular pipe turbulent flow model described in Fig. 3 when Ω = 2. In Fig. 6, line 610 represents the calculated value of velocity distribution, curve 620 represents the theoretical value of laminar velocity distribution, curve 630 represents the turbulent velocity distribution of a smooth circular pipe, and curve 640 represents the turbulent velocity distribution of a rough circular pipe. It can be seen that the flow velocity distribution in the radial direction shown by line 610 of the calculation result shown in Fig. 6 has a more suppressed flow velocity decrease than the result shown in Fig. 3.

[0079] 7 shows a calculation result 700 corresponding to FIG. 4, obtained by using the velocity calculated by this embodiment with the average velocity on the vertical axis and the relative position in the main flow direction on the horizontal axis. As shown in FIG. 7, according to the calculation method of this embodiment, the average velocity is well preserved along the length of the tube, indicating that the law of conservation of mass is well satisfied.

[0080] (3) Pipe friction coefficient To study the effect of the pipe friction factor, the circular pipe turbulent flow model explained in Figure 3 was used, with the unit of non-dimensional time T being the time it takes to pass through the pipe diameter at the average flow velocity, and while gradually increasing the velocity from T = 0 to 0.5, a ±5% disturbance was applied to each velocity component at all points from T = 0.5 to 4. After that, periods when a similar disturbance was applied and periods when it was not applied were alternated every 0.25 time, and the pipe friction factor during the period when no disturbance was applied was calculated.

[0081] The length of the calculation domain was set to 10 times the diameter, and the calculation was continued while occasionally replacing the inlet velocity distribution with the average value of the velocity distribution in the center. Figure 8 shows the calculation result 800, in which the vertical axis indicates the ratio of the calculated value of the pipe friction coefficient to the pipe friction coefficient of a smooth circular pipe, and the horizontal axis indicates the Reynolds number. Here, δ is the length of one side of the cubic cell, and d is the diameter. As shown in Figure 8, the algorithm of this embodiment shows that the smooth pipe friction coefficient is calculated when δ / d is 0.03 or less.

[0082] (4) Lower limit of the critical Reynolds number Next, a simulation was performed by continuously adding random numbers of -0.4 to 0.4 to the velocity components in three directions for the fluid phase in the center part of the radius 0 to 0.2 excluding the inlet part. FIG. 9 shows the results. The vertical axis of FIG. 9 indicates the value of unsteady Reynolds stress, and the horizontal axis indicates the radial position (dimensionless). At a Reynolds number of 2200, the turbulence in the center part is transmitted to the vicinity of the wall, but there is no sign of turbulence being stimulated by it and induced from the vicinity of the wall. On the other hand, at a Reynolds number of 2500, it is maximum near the wall, and it can be seen that turbulence is induced from the vicinity of the wall. It is shown that the lower limit of the critical Reynolds number when the initial turbulence is large is near 2300. As described above, it has been shown that the algorithm of this embodiment can simulate turbulent flows with a turbulence degree of up to about 2300 at a critical Reynolds number with a sufficiently high accuracy.

[0083] (5) Versatility Furthermore, the inventors have examined the applicability of the algorithm of this embodiment not only to turbulent pipe flows but also to flows with large turbulence around spheres and cylinders with large Reynolds numbers, with Ω=2. The inventors have examined whether the results can be applied to Figs. 10 and 11. In the simulation, the Reynolds number is 10,000, the calculation domain is 3 x 3 x 6 (length x width x depth = diameter), and the ratio of the length of one side of the cell to the diameter is 0.034. Figs. 10(a) and 10(b) show the velocity distributions by simulation around a sphere and a cylinder, respectively. Fig. 11 shows the same calculation results as Fig. 7, and in both cases, there is no change in the flow rate in the flow direction, and the law of conservation of mass is satisfied. Thus, it has been confirmed that this analysis method is effective not only for turbulent pipe flows but also for general flows with large turbulence.

[0084] As described above, according to this embodiment, it is possible to provide a turbulence analysis system and program that can obtain a solution that 1. satisfies the Navier-Stokes equations of motion, 2. largely satisfies the law of conservation of mass, and 3. does not diverge even in the presence of turbulence.

[0085] Although the present invention has been described above with reference to an embodiment, the present invention is not limited to the embodiment shown in the drawings, and may be modified within the scope of what a person skilled in the art can conceive, including other embodiments, additions, modifications, deletions, etc., and any aspect is within the scope of the present invention as long as it achieves the functions and effects of the present invention. [Explanation of symbols]

[0086] 100: Information processing device 101: Initial condition definition section 102: Speed ​​distribution calculation section 103: Speed ​​distribution calculation section

Claims

1. 1. A system for analyzing the Navier-Stokes equations, comprising: A means for defining initial conditions for pressure and velocity distributions in the model of the analysis; a first calculation means for determining a velocity distribution using a forward motion method; A second calculation means for calculating an unsteady Reynolds stress term from the calculated velocity distribution; a reduction means for reducing the Reynolds number in a range including strong turbulence that causes divergence, using a control variable obtained by vectoring the unsteady Reynolds stress term; a third calculation means for converging the pressure distribution including an inertia term calculated from the velocity distribution, an additional term for converging the velocity outflow, and a gradient of the unsteady Reynolds stress term required to satisfy the law of conservation of mass; Including, A turbulent flow analysis system, characterized in that the pressure distribution is calculated iteratively until a predetermined criterion is reached.

2. 2. The turbulent flow analysis system according to claim 1, wherein the means for reducing the Reynolds number varies a value of the reduced Reynolds number depending on the magnitude of the absolute value of the unsteady Reynolds stress term.

3. The turbulence analysis system according to claim 1 , wherein the predetermined criterion is the number of iterations or the convergence of the pressure distribution.

4. 1. A computer-executable program for solving the Navier-Stokes equations, comprising: Defining initial conditions for the pressure and velocity distributions in the model of the analysis; determining a velocity distribution using a forward method; calculating unsteady Reynolds stress terms from the calculated velocity distribution; reducing the Reynolds number in a range that includes strong turbulence causing divergence using the control variables that vectorize the unsteady Reynolds stress terms; Converging the pressure distribution including an inertia term calculated from the velocity distribution, an additional term for converging the velocity outflow, and a gradient of the unsteady Reynolds stress term required to satisfy the law of conservation of mass; Including, A program for repeatedly calculating the pressure distribution until a predetermined criterion is reached.

5. 5. The program according to claim 4, wherein the means for reducing the Reynolds number varies a value of the reduced Reynolds number depending on the magnitude of the absolute value of the unsteady Reynolds stress term.

6. The program according to claim 4 , wherein the predetermined criterion is the number of iterations or the convergence of the pressure distribution.

Images