JP2011186826A - Numerical calculation method, program and recording medium - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To achieve highly accurate calculation by using interpolation function in which a coefficient is determined by using proper information in downwind information in addition to upwind information, in a numerical calculation method used for computing partial differential equation. <P>SOLUTION: In the numerical calculation method, the interpolation function of the partial differential equation is a third-order polynomial, and a virtual point is set in the vicinity downstream of a reference point. A value of a variable at the virtual point is used for a restraint condition, to determine the coefficient of the interpolation function. Thus, the degree of the interpolation function is increased as well as non-physical numerical oscillation is reduced with respect to a high wave number, thereby enabling stable and highly accurate calculation. <P>COPYRIGHT: (C)2011,JPO&INPIT

Description

  The present invention relates to a numerical calculation method used for calculating a partial differential equation.

  In recent years, research and development relating to simulation of fluid physical phenomena using computers have become active. The object of simulation is various fluid motions caused by, for example, liquid, gas, or heat.

  When performing such numerical calculations, it is common to use a technique of solving partial differential equations relating to fluid variables (velocity, pressure, temperature, density, etc.) with respect to discretized space and time. This is because it is practically impossible to obtain a mathematical analytical solution (exact solution) of the partial differential equation. Therefore, it is important to develop a highly accurate numerical calculation method.

  In order to perform calculation with higher accuracy, not only the value of the variable but also the value of the differential coefficient (spatial gradient) of the variable or the interval integral average value (hereinafter also simply referred to as “integrated value”) is used as the dependent variable. The CIP method (Cubic-Interpolated Pseudo-Particle Scheme) and the IDO method (Interpolated Differential Operator Scheme), which are multi-moment methods, have been developed.

  The IDO method includes a conservative IDO method that strictly preserves mass, momentum, and energy that must be physically preserved, and a non-conservative IDO method that satisfies physical preservation within the range of high-precision numerical calculations. It can be divided roughly. In the storage type IDO method, for example, by setting the integral value of a variable as a dependent variable, a physical quantity related to the variable can be stored. The preservation type IDO method and the non-conservation type IDO method have advantages and disadvantages, respectively, and are selected and used according to the scene and situation.

  In addition, in the non-patent document 1 and the non-patent document 2 (substantially the same content as the non-patent document 1), the applicant of the present invention uses a calculation method using an interpolation function that locally approximates a partial differential equation for the conservative IDO method. Proposed. In the conservative IDO method, space is discretized by replacing a differential term in a partial differential equation with a differential value of an interpolation function.

  The technique disclosed in Non-Patent Document 1 and Non-Patent Document 2 is a high-precision numerical calculation method based on an interpolation function using the value of a variable on a lattice point (definition point) and an integral value in the conservative IDO method. It is. This method performs upwind interpolation used in compressible fluid analysis. Upwind interpolation is interpolation performed using information on the windward (original direction in which the fluid moves) rather than the grid point of interest (hereinafter also referred to as “reference point”). That is, when the reference point is j which is a grid point and the adjacent grid point in the windward direction is j-1, the value of the variable on the j point, the value on the j-1 point, and the integration between them The value was used to determine the coefficient of the interpolation function, which is a quadratic polynomial.

  This will be described with reference to FIG. Consider three points, a j−1 point, a j point, and a j + 1 point with a spatial interval Δx in the conservative IDO method. The reference point is j point, the j−1 point side is the windward side, and the j + 1 point side is the leeward side (opposite side on the windward side).

Assuming that the coordinates of the j point are x _j , the coordinates of the j−1 point and the j + 1 point are respectively x _j −Δx and x _j + Δx. In addition, the values of the variables at j−1 point, j point, and j + 1 point are assumed to be f _j−1 , f _j , and f _{j + 1} . In addition, an integral value from the j−1 point to the j point is ^x f _{j−1 / 2} , and an integral value from the j point to the j + 1 point is ^x f _{j + ½} .

Here, for example, an interpolation function F (x−x _j ) = a (xx _j ) ² + b (xx− _{j j} ) + c of a quadratic polynomial relating to a range from j−1 point to j point (a , B, c) is assumed. Since the interpolation function is a quadratic polynomial, the coefficient of the interpolation function can be determined if there are three constraint conditions. Therefore, three values f _j−1 , ^x f _{j−1 / 2} , and f _j are used. The constraint condition is as shown in the following formula (1).
_{^{Incidentally, ∫ xj xj-Δx F (}} x-x j) dx indicates that F a _{(x-x j)} is the integral of the section from the _x j _-Δx to _{x j} (hereinafter the same).
Therefore, (1 / Δx) ^∫xj _xj−Δx F (x−x _j ) dx represents an interval integral average value (an average value of the integrated values in the interval).

The value of the coefficient (a, b, c) of the interpolation function F (x−x _j ) can be determined by equation (1). Of these coefficients, b is particularly important. A value obtained by substituting x = _{x j} in the interpolation function F _{(x-x j)} the first derivative was _{∂F (x-x j) /} ∂x is b, the influence of this value gives the whole numerical calculation Because is big. Specifically, the value of b is expressed by the following equation (2).

  Here, in order to examine the accuracy (error) of the equation (2), when a Taylor variable is applied to the right side of the equation (2), the following equation (3) is obtained.

In other words, as can be seen from a comparison between Expression (2) and Expression (3), the error when Expression (2) is used for numerical calculation is equal to or less than the second term on the right side of Expression (3). That is, the calculation error is [Delta] x ² order, also, the numeric part of the coefficients in the top of the error term is 1/12.

Y. Imai, T. Aoki and K. Takizawa, “Conservative form of interpolated differential operator scheme for compressible and incompressible fluid dynamics”, Journal of Computational Physics, the Kingdom of the Netherlands, ELSEVIER, 2008, Vol. 227, Issue 4, p.2263-2285 Takayuki Aoki, Yosuke Imai, “Conservative IDO (Local Interpolation Differential Operator) Method”, Journal of Applied Mathematics, Iwanami “Applied Mathematics”, Iwanami Shoten, June 2008, Vol.18, No.2, p.32- 45

  Although the above-mentioned conventional method (method using equation (2)) can obtain a certain degree of calculation accuracy, it is the most important error in the calculation of numerical viscosity (fluid mechanics etc.) and eliminates the spatially fine profile. There is a problem that the numerical result becomes blurred (an error related to a spatially steep gradient becomes large) when the calculation is performed for a long time.

As a cause of this problem, for example, it is conceivable that information on the leeward side of the reference point is not used at all in determining the coefficient of the interpolation function. However, for example, in the method in which the degree of the interpolation function is increased by 1 to form a cubic polynomial, and the constraint value is set to 4 using the leeward integral value (“ ^x f _{j + 1/2} ” in the example of FIG. 12), In the region of high wavenumbers (space steep gradient), there are many calculation oscillations, which are not practical. In addition, it is not preferable to use information at a location far from the reference point because setting of the boundary condition becomes difficult or a calculation speed is delayed. Therefore, there has been a demand for a method for determining the coefficient of the interpolation function using appropriate information from the leeward information and performing high-accuracy calculation.

  Therefore, the present invention has been made in view of the above situation, and in the numerical calculation method used for the operation of the partial differential equation, in addition to the information on the windward, the coefficient is also calculated using appropriate information among the information on the windward. It is an object of the present invention to achieve highly accurate calculation by using an interpolation function that determines the above.

  In order to solve the above problems, in the present invention, the interpolation function of the partial differential equation is a cubic polynomial, a virtual point is assumed near the leeward side of the reference point, and the value of the variable at the virtual point is also used as a constraint condition. By doing so, the coefficient of the interpolation function is determined.

  According to the present invention, in the numerical calculation method used for the operation of the partial differential equation, in addition to the information on the windward, by using an interpolation function in which the coefficient is determined using appropriate information among the information on the windward, Highly accurate calculation can be realized.

It is a block diagram of the numerical calculation apparatus in this embodiment. It is a flowchart which shows the whole flow of the numerical calculation in this embodiment. It is explanatory drawing of the process of the determination of the 2nd interpolation function in 1st Embodiment. It is a flowchart showing the flow of the calculation for the determination of the 2nd interpolation function in 1st Embodiment. 3 is a graph showing an example of numerical calculation results in Example 1. It is a graph which shows the example of the numerical calculation result in Example 1, (a) is a graph which shows an amplitude ratio, (b) is a graph which shows a phase. It is a graph which shows the example of the numerical calculation result in Example 1, (a) is a graph which shows an amplitude ratio error, (b) is a graph which shows a phase error. It is a graph which shows the example of the numerical calculation result in Example 2, (a) is a graph which shows the example of the calculation result by a prior art (secondary upwind), (b) is the numerical calculation (3) of this embodiment. It is a graph which shows the example of the calculation result by the next upwind). It is a graph which shows the example of the numerical calculation result in Example 2, (a) is a graph which shows the example of the numerical calculation result by a prior art (secondary upwind), (b) is this embodiment (tertiary wind) It is a graph which shows the example of the numerical calculation result by the above. It is explanatory drawing of the process of the determination of the 2nd interpolation function in 2nd Embodiment. It is explanatory drawing of the process of the determination of the 2nd interpolation function in 3rd Embodiment. It is explanatory drawing of the process of the determination of the interpolation function in a prior art.

  Hereinafter, regarding a mode for carrying out the present invention (hereinafter referred to as “embodiment”), the first embodiment, the second embodiment, and the third embodiment (the generic name thereof is referred to as “present embodiment”). This will be described with reference to the drawings as appropriate. In the present embodiment, the case of the storage type IDO method will be described as an example. First, the configuration of a numerical calculation apparatus that executes the numerical calculation method according to the present embodiment will be described with reference to FIG.

  As shown in FIG. 1, a numerical calculation device 10 is a computer device that includes an input unit 1, a storage unit 2, a memory 3, a calculation unit 4, and an output unit 5. It can be realized by a computer or the like.

The input unit 1 is means for a user of the numerical calculation device 10 (hereinafter simply referred to as “user”) to input various types of information, such as a keyboard and a mouse.
The storage unit 2 is a means for storing various information, such as a ROM (Read Only Memory), a hard disk, and the like. Here, in order to perform numerical calculation (simulation) of the fluid or the like, the storage unit 2 is a definition point indicating a place where the numerical calculation is performed, a lattice point, a physical quantity (variable), initial conditions of its differential coefficient, A partial differential equation representing a physical phenomenon of fluid), an interpolation function (first interpolation function, second interpolation function), various programs, and the like are stored.

The memory 3 is temporary storage means used by the arithmetic unit 4 and is, for example, a RAM (Random Access Memory).
The computing unit 4 performs computation (processing) based on various information (programs, data, etc.) stored in the storage unit 2, and is, for example, a CPU (Central Processing Unit). In addition, although the calculation part 4 is set to one here, you may make it perform parallel calculation by providing multiple.
The output unit 5 outputs the result of the calculation by the calculation unit 4 and is, for example, a display.

Next, the overall flow of numerical calculation in the present embodiment will be described with reference to FIG. 2 (see FIG. 1 as appropriate).
As shown in FIG. 2, first, when the user performs an operation related to the start of numerical calculation by the input unit 1, the calculation unit 4 is a definition indicating a place where the calculation is performed based on information stored in the storage unit 2. A grid point as a point, an initial condition of each physical quantity, and the like are set (step S1). The grid points are set, for example, as “512” in the case of one dimension, “128 × 128” in the case of two dimensions, and “16 × 16 × 8” in the case of three dimensions. In addition, as initial conditions for each physical quantity, variable values at the respective lattice points, differential coefficients (spatial gradients), integral values, and the like of these values are set. Examples of variables include fluid speed, pressure, temperature, density, and the like.

  Next, based on the information stored in the storage unit 2, the calculation unit 4 sets a governing equation including a partial differential equation related to time and space, an interpolation function for correcting the value of each physical quantity, and the like (step) S2). Examples of the governing equation include a motion equation of compressible fluid and a mass conservation equation of fluid. The process of setting the interpolation function in step S2 is a feature of the present invention, and a detailed description thereof will be given later.

  Subsequently, the calculation unit 4 uses the various information set in step S1 and step S2, and based on the value of each physical quantity at each grid point at time “t”, each calculation at each grid point at time “t + 1”. The calculation of calculating the value of the physical quantity is performed for t = 0, 1, 2,... N (steps S3 to S5).

Next, the calculating part 4 outputs the value of each physical quantity calculated by step S3-step S5 to the output part 5 (step S6). Thus, the numerical calculation is completed.
Of the processes in steps S1 to S6, the interpolation function setting process in step S2 will be described below. Other processes are described in Japanese Patent No. 4304277 relating to the patent application previously filed by the present applicant. No. and JP 2008-197949 A, detailed description thereof is omitted here.

<First Embodiment>
Next, as the first embodiment, the interpolation function setting process in step S2 will be described. First, an image of setting an interpolation function will be described with reference to FIG. Here, a case will be described in which one value is calculated using the first interpolation function, and the second interpolation function is set (determining the coefficient) using that value and the other three values. The second interpolation function is a setting target interpolation function, and the first interpolation function is a means for creating the second interpolation function.

  As shown in FIG. 3, in the preserved IDO method, consider three points, j−1 point, j point, and j + 1 point, where the spatial interval is Δx. The reference point is j, the j-1 point is the windward, and the j + 1 point is the leeward.

Assuming that the coordinates of the j point are x _j , the coordinates of the j−1 point and the j + 1 point are respectively x _j −Δx and x _j + Δx. In addition, the values of the variables at j−1 point, j point, and j + 1 point are assumed to be f _j−1 , f _j , and f _{j + 1} . In addition, an integral value from the j−1 point to the j point is ^x f _{j−1 / 2} , and an integral value from the j point to the j + 1 point is ^x f _{j + ½} .

At this time, first, the first interpolation function H (x−x _j ) = a (x− _{j j} ) ² + b (x− _{j j} ) + c (... (4a) relating to the range from point j to point j + 1. )) Coefficients (a, b, c) are determined. Since the first interpolation function is a quadratic polynomial, the coefficient of the first interpolation function can be determined if there are three constraint conditions. Therefore, three values f _j , ^x f _{j + 1/2} and f _{j + 1} are used. The constraint condition is as the following equation (4b).

The value of the coefficient (a, b, c) of the first interpolation function H (x) is determined by the expression (4b) (the description of the value of the coefficient (a, b, c) is omitted). Using this first interpolation function H (x−x _j ), the value of f _{j + 1/2} can be obtained as in the following equation (5). That is, the first interpolation function H (x−x _j ) is a means for obtaining the value of f _{j + 1/2} .

  That is, j + 1/2 points (intermediate points between j and j + 1) that are not grid points are used as virtual points, and the value of f on the virtual points is calculated using the first interpolation function.

Next, the coefficients (c ₃ , c ₂ , c ₁ , c ₀ ) of the second interpolation function shown in the following equation (6) regarding the range from the j−1 point to the j + 1/2 point are determined.
^{_{F (x-x j) =}} c 3 (x-x j) 3 + c 2 (x-x j) 2 + c 1 (x-x j) + c 0
... Formula (6)

Since the second interpolation function is a cubic polynomial, if there are four constraint conditions, the coefficient of the second interpolation function can be determined. Therefore, four values of f _j−1 , ^x f _{j−1 / 2} , f _j , and f _{j + 1/2} are used. The constraint condition is as the following formula (7).

These determine the values of the coefficients (c ₃ , c ₂ , c ₁ , c ₀ ) of the second interpolation function F (x). Of these factors, c ₁ is particularly important. Value when substituting x = _{x j} in the second interpolation function F _{(x-x j)} the first derivative was _{∂F (x-x j) /} ∂x is _{c 1,} the value is numerical This is because the influence on the whole is large. The value of c ₁ is specifically by the following equation (8).

Here, in order to investigate the accuracy (error) of Expression (8), when a Taylor variable is substituted for the right side of Expression (8), the following Expression (9) is obtained.

In other words, as can be seen from a comparison between Expression (8) and Expression (9), the error when Expression (8) is used for numerical calculation is equal to or less than the second term on the right side of Expression (9). That is, the calculation error is [Delta] x ³ order, also, the numeric part of the coefficients in the top of the error term is 1/90. As described above, in the interpolation function quadratic polynomial in the prior art (equation (2), (3)), and calculation error is [Delta] x ² order, also, the numeric part of the coefficients in the top of the error term Compared with 1/12, it can be seen that the calculation accuracy when using the equation (8) is significantly higher than in the case of the prior art. In other words, when comparing the second term or less on the right side of Equation (3) with the second term or less on the right side of Equation (9), the differential rank in each of the most significant terms is different. Nonetheless, as described above, it can be said that the calculation error in the case of using Equation (2) is obviously larger due to the difference between the order and the coefficient.

Similarly, the values of other coefficients are as shown in the following equation (10).

Further, when the windward and leeward are interchanged, the virtual point can be calculated in the same manner as described above, with the virtual point being j-1 / 2 point (the intermediate point between j-1 point and j). Then, in the second interpolation function, c _{1 'is} a coefficient corresponding to c ₁ is as shown in the following equation (11).

Incidentally, _{similarly, c 3, c 2, c} 3 is a coefficient corresponding to _{c 0 ', c 2',} c 0 ' , the value of the other coefficients, so that the following equation (12) .

  Next, the calculation (processing) for determining the coefficient of the second interpolation function in the first embodiment (processing for setting the interpolation function in step S2 in FIG. 2) will be described with reference to FIG. (See FIGS. 1 and 3 as appropriate). The calculation unit 4 of the numerical calculation device 10 performs each calculation while using the storage unit 2 and the memory 3.

First, the calculation unit 4 acquires three values f _j , ^x f _{j + 1/2} , and f _{j + 1} (step S21).
Next, the calculating part 4 determines the coefficient of a 1st interpolation function using the acquired three values (step S22).

Subsequently, the calculation unit 4 calculates a value of f _{j + 1/2} using the first interpolation function determined by the coefficient (step S23).
Next, the calculation unit 4 acquires four values of f _j , ^x f _{j + 1/2} , f _{j + 1} , and f _{j + 1/2} (step S24).

  Finally, the calculation unit 4 determines the coefficient of the second interpolation function using the acquired four values (step S25).

As described above, according to the first embodiment, in the numerical calculation method used for the calculation of the partial differential equation by the numerical calculation device 10, f _j and ^x f _{j + 1/2} conventionally used for upwind interpolation. , F _{j + 1} , and the coefficient of the second interpolation function is determined using the value f _{j + 1/2} at the virtual point j + 1/2 calculated using the first interpolation function. Therefore, the numerical calculation device 10 can realize high-accuracy calculation by using the second interpolation function in which the coefficient is determined using appropriate information on the downwind in addition to the information on the upwind.

(Fourier analysis)
Next, with reference to FIG. 5 to FIG. 7, verification of accuracy and stability using Fourier analysis will be described regarding an example of the numerical calculation method of the first embodiment. In the first embodiment, Fourier analysis is performed by substituting discrete variables that are Fourier-transformed into a linear advection equation, and the relationship between the wave number, the Couran number, the phase, and the amplitude ratio is graphed. For the time integration, the fourth-order Runge-Kutta method was used.

  Regarding Example 1, “3.Fourier analysis (p.2272-2276)” in Non-Patent Document 1 and “5 Fourier analysis (p.38-41)” in Non-Patent Document 2 described above are used. Since there is a detailed description, detailed description is omitted here. Compared with those described in Non-Patent Documents 1 and 2, the difference in the first embodiment is that the interpolation function is a cubic polynomial and the method for determining the coefficient of the interpolation function that is the cubic polynomial. It is.

  In the graph shown in FIG. 5, the eigenvalues of the matrix describing the time evolution of the partial differential equation are displayed on the coordinates where the vertical axis is the Courant number and the horizontal axis is the wave number. In such a graph, the numerical calculation is more stable as the line indicating the eigenvalue “1.0” is positioned upward. Here, it can be seen that the numerical value calculation method according to the first embodiment is stable because the line L51 indicating the eigenvalue “1.0” is located considerably above.

  In the graph shown in FIG. 6A, the ordinate is the amplitude ratio, the abscissa is the wave number, and the exact solution and the conservative IDO method (secondary upwind: a conventional interpolation function of a second order polynomial is used. Upwind interpolation using the value of the upwind interpolation method and the conservative IDO method (third-order upwind: third-order polynomial interpolation function in the first embodiment (here, downwind information is also used) The value obtained by the above method) is displayed.

  In such a graph, naturally, the closer to the exact solution, the higher the accuracy of numerical calculation. 6A, the storage type IDO method (third-order upwind) that is the method of the first embodiment is more strict than the storage-type IDO method (secondary upwind) that is the conventional method. It can be seen that the numerical calculation for the amplitude ratio is highly accurate and is close to the solution.

  In the graph shown in FIG. 6B, the vertical axis is the phase, the horizontal axis is the wave number, the exact solution, and the conservative IDO method (center: center interpolation using a conventional quadratic function interpolation function (3 Interpolation based on the center of the two grid points), a value based on the conservative IDO method (secondary upwind), and a value based on the conservative IDO method (third order upwind). And are displayed.

  Here, there is no significant difference in the difference between the conservative IDO method (center), the conservative IDO method (secondary upwind), and the conservative IDO method (third-order upwind). That is, it can be seen that the conservative IDO method (third-order upwind), which is the method of the first embodiment, is not particularly different from other conventional methods with respect to the accuracy of numerical calculation for the phase.

  Further, in the graph shown in FIG. 7A, the vertical axis is the amplitude ratio error, the horizontal axis is the wave number, and the value by the conservative IDO method (secondary upwind) and the conservative IDO method (third-order upwind). The value by and is displayed. In view of this, the conservative IDO method (third-order upwind) that is the method of the first embodiment has a significantly smaller amplitude ratio error than the conservative IDO method (second-order upwind) that is the conventional method. It can be seen that the numerical calculation of the amplitude ratio is highly accurate.

  In the graph shown in FIG. 7B, the vertical axis is the phase error, the horizontal axis is the wave number, the value by the conservative IDO method (center), the value by the conservative IDO method (secondary upwind), and the conservative form. The value by the IDO method (third-order upwind) is displayed.

  Here, there is no significant difference between the preserved IDO method (center), the preserved IDO method (secondary upwind), and the preserved IDO method (third-order upwind). That is, it can be seen that the conservative IDO method (third-order upwind), which is the method of the first embodiment, is not particularly different from other conventional methods with respect to the accuracy of numerical calculation for the phase.

  Thus, as can be seen from Example 1, although the numerical calculation method of the first embodiment has almost the same calculation accuracy for the phase as the conventional technology, the calculation accuracy for the eigenvalue and the amplitude ratio is higher than that of the conventional technology. Significantly high and comprehensive calculation accuracy is significantly high.

(Growth process simulation of compressible Rayleigh-Taylor instability)
Next, with reference to FIGS. 8 and 9, an example in which the numerical calculation method of the first embodiment is used for growth process simulation of compressible Rayleigh-Taylor instability will be described. In this example 2, as a calculation example of a compressible fluid, when a heavy (high density) fluid is placed on a light (low density) fluid, a disturbance applied to the interface between the two fluids is shown. It is a graph of two fluids when simulating the process of growth due to Rayleigh-Taylor instability. In addition, about the two fluids at that time, FIG. 8 shows a state seen from the side (a state seen in the horizontal direction), and FIG. 9 shows a state seen from the top (a state seen vertically downward). ing. Here, the numerical calculation method of the first embodiment is applied to the advection term calculation part in the compressible fluid calculation.

  Since this Example 2 has a detailed description in “6 Calculation Example 6.1 Compressible Rayleigh-Taylor Instability (p.41-43)” in Non-Patent Document 2 described above, a detailed explanation is given here. Omitted. Compared with the description in Non-Patent Document 2, the only difference in the second embodiment is that the interpolation function is a cubic polynomial and the method for determining the coefficient of the interpolation function that is the cubic polynomial.

  The graph shown in FIG. 8A shows the state of the two fluids when the calculation is performed by the conservative IDO method (secondary upwind). The graph shown in FIG. 8B shows the state of the two fluids when the calculation is performed by the conservative IDO method (third-order upwind). Both graphs have a resolution of 256 × 256.

  From these, it can be seen that, in the graph shown in FIG. 8A, the details are blurred, whereas in the graph shown in FIG. 8B, the details can be clearly displayed. That is, it can be seen that the storage type IDO method (third-order upwind), which is the method of the first embodiment, has higher numerical calculation accuracy than the storage-type IDO method (secondary upwind) that is the conventional method.

  Similarly, the graph shown in FIG. 9A shows the state of two fluids when the calculation is performed by the conservative IDO method (secondary upwind). The graph shown in FIG. 9B shows the state of the two fluids when the calculation is performed by the conservative IDO method (third-order upwind). Both graphs have a resolution of 128 × 128.

  From these, it can be seen that, in the graph shown in FIG. 9A, the details are blurry, whereas in the graph shown in FIG. 9B, the details can be clearly displayed. That is, it can be seen that the storage type IDO method (third-order upwind), which is the method of the first embodiment, has higher numerical calculation accuracy than the storage-type IDO method (secondary upwind) that is the conventional method.

Thus, as can be seen from Example 2, the numerical calculation method of the first embodiment has higher calculation accuracy than the prior art. The reason will be described as follows. That is, in such a simulation, initial linear growth γ = √gka (g: gravitational acceleration, k: wave number of disturbance, a = (ρ _H −ρ _L ) / (ρ _H + ρ _L ), ρ _H : heavy fluid Past the density of ρ _L : the density of light fluid), it moves to non-linear growth with Kelvin-Helmholtz instability, and the growth rate of long-wave waves increases, while forming a very complex interface, Free Fall (H to 0.06 agt ² , h: falling speed, t: time).

  And when calculating with the conventional method (conservative IDO method (secondary upwind)) with a large numerical viscosity, the wave with a large wave number is numerically attenuated, and the transition to Free Fall and the fall acceleration are slowed down. May end up. On the other hand, in the numerical calculation method according to the first embodiment (preservation type IDO method (third-order upwind)), such a delay is small and calculation can be performed with high accuracy.

In the numerical calculation method of the first embodiment, the value of f _{j + 1/2} is calculated by the _first interpolation function, and the coefficient of the second interpolation function is determined using the value of f _{j + 1/2} as the constraint condition. . As a result, not only the order of the interpolation function is increased by one, but also non-physical numerical vibration is reduced with respect to the high wave number, and stable and highly accurate calculation can be performed.

  A person skilled in the art creates a program for causing a computer to execute the numerical calculation method according to the first embodiment, and further stores the program in a recording medium such as a CD (Compact Disc) or a DVD (Digital Versatile Disk). Can be recorded.

Second Embodiment
Next, as a second embodiment, an interpolation function setting process in step S2 (see FIG. 2) will be described. In the second embodiment, overlapping description of the same points as in the first embodiment will be omitted as appropriate, and differences from the first embodiment will be mainly described.

An image of setting an interpolation function will be described with reference to FIG. The difference from the first embodiment is that the value calculated using the first interpolation function is f ^tmp and the coordinates of the f ^tmp are x _j + α · Δx (0 <α <1).

Similarly to the first embodiment, when the first interpolation function (formula (4a)) and the second interpolation function (formula (6)) are used, c _0, c ₁ , c ₂ , c ₃ The value is obtained as follows.

Here, when considered in the same manner as in equations (8) and (9), the exact solution that is the basis of c ₁ in equation (13) can be expressed as in the following equation (14).

Further, the equation (14) is modified into the following equation (15).

Here, in Equation (15), the third-order differential term affects the wave phase velocity, and the fourth-order differential term affects the wave viscosity dissipation. And in the conservative IDO method, the phase velocity of the wave is slower than the exact solution, so it is better to correct the phase velocity to increase, that is, it is preferable to satisfy E ₃ <0. Further, it is preferable to introduce a stable numerical viscosity, that is, it is preferable to satisfy E ₄ > 0. Therefore, if it calculates so that it may satisfy | fill, it will be understood that (alpha) should just satisfy | fill following Formula (16).
(√17-3) / 4 <α ≦ 1/2 Formula (16)

Of α satisfying equation (16), as can be seen from equations (14) and (15), when α = 1/2, E ₃ = 0, and the term of Δx ² in equation (15) is since disappear, c ₁ is the tertiary precision (magnitude of the error is proportional to [Delta] x ^3). Then, when α ≠ 1/2, c ₁ has a secondary accuracy (the magnitude of the error is proportional to Δx ² ). Therefore, in the first embodiment, the first interpolation is performed by using j + 1/2 point (middle point between j point and j + 1. That is, the point when α = 1/2 in the second embodiment) as a virtual point. It also shows the effectiveness of calculating (interpolating) the value of f on the virtual point using a function. The same applies to the accuracy of c ₂ and c ₃ .

As described above, according to the second embodiment, in the numerical calculation method used for the calculation of the partial differential equation by the numerical calculation device 10, f _j and ^x f _{j + 1/2} that have been conventionally used for upwind interpolation. , F _{j + 1} , and the coefficient of the second interpolation function is determined by using the value f _{j + α · Δx} at the virtual point j + α · Δx calculated using the first interpolation function. Therefore, the coefficient of the second interpolation function is determined using the information about the leeward point determined by setting α within the range satisfying the expression (16) in addition to the intermediate point between the j point and the j + 1 point, and the second By using the interpolation function of 2, high-accuracy calculation can be realized.

<Third Embodiment>
Next, as a third embodiment, an interpolation function setting process in step S2 (see FIG. 2) will be described. In the third embodiment, overlapping description of the same points as at least one of the first embodiment and the second embodiment will be omitted as appropriate, and differences from the first embodiment and the second embodiment will be mainly described. .

  An image of setting an interpolation function will be described with reference to FIG. The difference from the second embodiment is that the range of the first interpolation function is from j−1 point to j + 1 point, and the first interpolation function is not the second-order polynomial of equation (4a), but is expanded. 4th order polynomial (detailed explanation is omitted), and accordingly, the range of α is −1 <α <0, 0 <α <1 (that is, the virtual point is j−1 points) To j + 1 point, except for j point).

As in the first embodiment, the values of c _0, c ₁ , c ₂ , and c ₃ can be obtained by using the first interpolation function of the fourth-order polynomial and the second interpolation function (equation (6)). Is obtained by the following equation (17).

Here, when considered in the same manner as in equations (8) and (9), the exact solution that is the source of c ₁ in equation (17) can be expressed as in the following equation (18).

Further, the equation (18) is transformed into the following equation (19).

Here, considering E ₄ in Equation (19), since the denominator includes the term (2a + 1), when the denominator approaches 0 (α = −1 / 2), the value of E ₄ becomes infinite ( Asymptotically close to (infinity or small), the calculation is not stable. In consideration from the viewpoint of numerical viscosity, the value of E _4, it is desirable to satisfy the E _4> 0. In the range of −3/5 <α <−1/2, E ₄ > 0 is satisfied. However, since the value of E ₄ may gradually approach infinity, it is not preferable from the viewpoint of the stability of calculation. Therefore, the range of 0 <α <1 satisfies E ₄ > 0, and the value of E ₄ does not asymptotically approach infinity, which is preferable for calculation stability. This also applies to c ₂ and c ₃ .

  As described above, according to the third embodiment, the range of the first interpolation function is expanded from j−1 point to j + 1 point, and the first interpolation function is a quartic polynomial, so that the value of α is obtained. Regardless of the calculation, the third order accuracy can be obtained. In this case, the range of α is preferably 0 <α <1, and thus, the effectiveness of using appropriate downwind information, which is one of the features of the present invention, has been demonstrated.

When generality is derived from the above-described three embodiments, (1) a first interpolation function is created using five or more physical quantities, or (2) f _j as in the first embodiment. calculates the _{f j + 1/2} in a virtual point j + 1/2-point using a ^{_{x f j + 1/2,}} f j + 1, by any of the condition that, it can be said that the _{c 1} may be a tertiary accuracy.

Although description of this embodiment is finished above, the aspect of the present invention is not limited to these.
For example, the present invention is not limited to the above-described embodiments, and can be widely applied to general simulation of fluid physical phenomena using partial differential equations.

Further, the number of dimensions of the space in the partial differential equation to be used may be any number of one dimension or more.
Moreover, the physical variable regarding the fluid to which the numerical calculation method of the present embodiment is applied may be any variable such as density, speed, pressure, temperature, and the like.

Further, the case where the j-1 point side is taken up and the j + 1 point side is taken down has been mainly described, but conversely, the j-1 point side may be taken down and the j + 1 point side may be taken up. In other words, when upwind and downwind are interchanged in successive calculations, the interpolation function can be used properly in each scene.
Further, the present invention can be applied not only to the storage type IDO method but also to other calculation methods such as the CIP method.

Further, when determining the coefficient of the first interpolation function, three values f _j , ^x f _{j + 1/2} and f _{j + 1} are used, but other values may be used.
Further, the first interpolation function in the first embodiment and the second embodiment described above and the first interpolation function of the fourth-order polynomial in the third embodiment are examples, and other interpolation functions are used. Also good.

In addition, as the coefficient of the second interpolation function, the values according to the equations (8) and (10) to (12) may be used instead of the values themselves. The value according to them may be, for example, a value obtained by slightly changing the numerical part on the right side in the equations (8) and (10) to (12). For example, it is clear that even if “5f _j−1 ” is replaced with “4.99999f _j−1 ” on the right side of Expression (8), the calculation accuracy is not significantly affected. Also, for example, when a virtual point is considered between j point and j + 1, a point exactly in the middle of j point and j + 1 is set as a virtual point, but other than that, a point slightly deviated from the intermediate point may be set as a virtual point. Good. In that case, the coefficient of the second interpolation function is slightly different from the values shown in the equations (8) and (10) to (12), but if the range does not significantly affect the accuracy of the calculation, Those slightly different values also fall under “similar values”.

  In addition, the specific configuration of the numerical calculation device and each flowchart can be appropriately changed without departing from the spirit of the present invention.

DESCRIPTION OF SYMBOLS 1 Input part 2 Memory | storage part 3 Memory 4 Operation part 5 Output part 10 Numerical calculation apparatus

Claims (8)

A governing equation including a partial differential equation representing a physical phenomenon of a fluid and a cubic that locally approximates the partial differential equation with respect to the value of one or more types of fluid at each of a plurality of definition points in space. A numerical calculation method by a numerical calculation device that updates in a predetermined minute time unit using an interpolation function that is a polynomial and a coefficient is determined for each definition point,
The numerical calculation device includes a storage unit and a processing unit,
The storage unit stores the governing equation, the interpolation function, the plurality of definition points, and the value of the variable for each definition point,
Of the defined points, three points are j−1 points, j points, and j + 1 points, a reference point is a j point, and a j−1 point side is an upwind that is the original direction in which the fluid moves. , J + 1 point side is the leeward direction opposite to the windward side,
For each of the defined points j−1, j, and j + 1, the coordinates are x _j −Δx, x _j , x _j + Δx, and the value of one type of variable among the variables for each defined point Is f _j−1 , f _j , f _{j + 1} ,
The integral value of the variable from the defined point j-1 to the j point is ^x f _{j−1 / 2} , and the integral value of the variable from each defined point j to the j + 1 point is ^x f _{j + 1/2} ,
When the interpolation function is the following equation (6):
^{_{F (x-x j) =}} c 3 (x-x j) 3 + c 2 (x-x j) 2 + c 1 (x-x j) + c 0
... Formula (6)
The processor is
When the interpolation function represented by the equation (6) is used for the j point that is the reference point with reference to the storage unit, a value represented by the following equation (8) or a value equivalent thereto is used as the coefficient c _1. A numerical calculation method characterized by the use.
The processor is
With reference to the storage unit, when using the interpolation function represented by the equation (6) for the reference point j, it is represented by the following equation (10) as the coefficients c ₃ , c ₂ , c _0. 2. The numerical calculation method according to claim 1, wherein a value or a value equivalent thereto is used.
When the windward and leeward are swapped, the j-1 point side is leeward and the j + 1 point side is leeward,
The processor is
With reference to the storage unit, when using the interpolation function represented by the equation (6) for the reference point j, the value of c ₁ ′ represented by the following equation (11) or the coefficient c ₁ The numerical calculation method according to claim 1, wherein a value equivalent thereto is used.
The processor is
When the interpolation function represented by the above equation (6) is used for the j point that is the reference point with reference to the storage unit, it is represented by the following equation (12) as the coefficients c ₃ , c ₂ , and c _0. The numerical calculation method according to claim 3, wherein values of c ₃ ′, c ₂ ′, and c ₀ ′ or values equivalent thereto are used.
A governing equation including a partial differential equation representing a physical phenomenon of a fluid and a cubic that locally approximates the partial differential equation with respect to the value of one or more types of fluid at each of a plurality of definition points in space. A numerical calculation method by a numerical calculation device that updates in a predetermined minute time unit using an interpolation function that is a polynomial and a coefficient is determined for each definition point,
The numerical calculation device includes a storage unit and a processing unit,
The storage unit stores the governing equation, the interpolation function, the plurality of definition points, and the value of the variable for each definition point,
Of the defined points, three points are j−1 points, j points, and j + 1 points, a reference point is a j point, and a j−1 point side is an upwind that is the original direction in which the fluid moves. , J + 1 point side is the leeward direction opposite to the windward side,
For each of the defined points j−1, j, and j + 1, the coordinates are x _j −Δx, x _j , x _j + Δx, and the value of one type of variable among the variables for each defined point Is f _j−1 , f _j , f _{j + 1} ,
The integral value of the variable from the defined point j-1 to the j point is ^x f _{j−1 / 2} , and the integral value of the variable from each defined point j to the j + 1 point is ^x f _{j + 1/2} ,
When the interpolation function is the following equation (6):
^{_{F (x-x j) =}} c 3 (x-x j) 3 + c 2 (x-x j) 2 + c 1 (x-x j) + c 0
... Formula (6)
The processor is
When the interpolation function represented by the equation (6) is used for the j point that is the reference point with reference to the storage unit, a value represented by the following equation (13) or a value equivalent thereto is used as the coefficient c _1. A numerical calculation method characterized by the use.
However, (√17-3) / 4 <α ≦ 1/2 A governing equation including a partial differential equation representing a physical phenomenon of a fluid and a cubic that locally approximates the partial differential equation with respect to the value of one or more types of fluid at each of a plurality of definition points in space. A numerical calculation method by a numerical calculation device that updates in a predetermined minute time unit using an interpolation function that is a polynomial and a coefficient is determined for each definition point,
The numerical calculation device includes a storage unit and a processing unit,
The storage unit stores the governing equation, the interpolation function, the plurality of definition points, and the value of the variable for each definition point,
Of the defined points, three points are j−1 points, j points, and j + 1 points, a reference point is a j point, and a j−1 point side is an upwind that is the original direction in which the fluid moves. , J + 1 point side is the leeward direction opposite to the windward side,
For each of the defined points j−1, j, and j + 1, the coordinates are x _j −Δx, x _j , x _j + Δx, and the value of one type of variable among the variables for each defined point Is f _j−1 , f _j , f _{j + 1} ,
The integral value of the variable from the defined point j-1 to the j point is ^x f _{j−1 / 2} , and the integral value of the variable from each defined point j to the j + 1 point is ^x f _{j + 1/2} ,
When the interpolation function is the following equation (6):
^{_{F (x-x j) =}} c 3 (x-x j) 3 + c 2 (x-x j) 2 + c 1 (x-x j) + c 0
... Formula (6)
The processor is
When the interpolation function represented by the equation (6) is used for the reference point j, with reference to the storage unit, a value represented by the following equation (17) or a value equivalent thereto is used as the coefficient c _1. A numerical calculation method characterized by the use.
However, 0 <α <1   The program for making a computer perform the numerical calculation method as described in any one of Claims 1-6.   A recording medium on which the program according to claim 7 is recorded.