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

Description

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

In recent years, many efforts have been made to reproduce natural phenomena using computers. For example, numerical simulation such as fluid dynamics is an example, and the object is various such as liquid, gas, and electromagnetic wave.
When performing such a numerical simulation, a technique is often used in which a partial differential equation relating to two or more variables (hereinafter referred to as “multivariables”) is solved in terms of discretized space and time. . Examples of variables include speed, pressure, temperature, density, and the like.

In a numerical solution method (simulation) in a multivariable partial differential equation, high temporal stability (hereinafter, simply referred to as “stability”) is required for coupling between variables. This is because if the stability of coupling between variables is low, the calculation fails in a short time.
Therefore, conventionally, in the SMAC (Simplified Marker and Cell) method in fluid analysis and the FDTD (Finite Difference Time Domain) method in electromagnetic wave analysis, a staggered variable is set by shifting the definition points for each variable. Many lattices were used. This is because the staggered lattice has high coupling stability between variables.

  In addition, staggered grids have been often used in multi-moment schemes, which are calculation methods that use not only the values of each variable but also their spatial gradients (differential coefficients). However, the staggered grid has high stability in coupling between variables, but has a drawback in terms of calculation accuracy and ease of parallel calculation algorithm.

On the other hand, a collocated grid that defines all variables and spatial gradients on the same grid is superior to a staggered grid in terms of calculation accuracy and ease of parallel calculation algorithms.
For example, the calculation accuracy in the multi-moment scheme is second order of the staggered lattice and fourth order of the collocated lattice.
In addition, the collocated grid has the advantage that it can be easily combined with the AMR (Adaptive Mesh Refinement) method and the Cut-Cell method, which are expected to be developed in the future, compared to the staggered grid. Have.

  However, collocated grids can be used only under special conditions because the stability of coupling between variables is low and calculations often fail in a short time, and the application range is very limited. There was a problem of being.

For example, in an IDO (Interpolated Differential Operator) method that is a multi-moment scheme, a central interpolation function F (x) with respect to three points, j−1 points, j points, and j + 1 points having an interval Δx is expressed as j. Centering on the point, it can be expressed as a quintic function such as the following equation (1).
Note that f _j−1 , f _j , and f _{j + 1} represent the values of variables at j−1 point, j point, and j + 1 point, respectively, and fx _{, j−1} , fx _{, j} , fx _{, j + 1} represents the differential coefficient of the variable at j−1 point, j point, and j + 1 point, respectively. X represents coordinates.

F (x) = A ₅ x ⁵ + A ₄ x ⁴ + A ₃ x ³ + A ₂ x ² + f _{x, j} x + f _j (1)
Also, a differential expression of F (x) F _x (x) is expressed by equation (2).
F _x (x) = 5A ₅ x ⁴ + 4A ₄ x ³ + 3A ₃ x ² + 2A ₂ x + f _{x, j} (2)
The constraint conditions are F (Δx) = f _{j + 1} , F (−Δx) = f _j−1 , F _x (Δx) = fx _{, j + 1} , F _x (−Δx) = fx _{, j− 1} and these determine the values of A _{2 to} A ₅ .

Under these conditions, the spatial differential term in the governing equation (such as the partial differential equation used in the calculation) is expressed as follows by the equation (1), the equation (2), and the equation obtained by further differentiating the equation (2). f / ∂x = F _x (0) = fx _{, j} , ∂ ² f / ∂x ² = F _xx (0) = 2A ₂ ,... Note that F _xx (x) represents an expression obtained by differentiating F _x (x).

Here, the following equations (3) and (4), which are multivariable partial differential equations, are considered. F and u are variables, t is time, and x is coordinates.

On the collocated grid, spatial discretization such as the following Expression (5) and Expression (6) is performed using the interpolation function (Expression (1) and Expression (2)). Note that ∂f / ∂t | _j and ∂u / ∂t | _j represent differential coefficients at j points obtained from an interpolation function or the like.
As can be seen from the equations (5) and (6), only j-point variables appear in the discretization equation, and thus the coupling stability between the variables is poor (low).

  Accordingly, the present invention has been made in view of the above problems, and an object of the present invention is to provide a numerical calculation method with high stability of coupling between variables in a numerical solution of a partial differential equation using a collocated grid. And

In order to solve the above problems, a numerical calculation method according to the present invention is a numerical calculation method by a numerical calculation device that performs numerical calculation related to a variable on each definition point, and the numerical calculation device includes a storage unit and a processing unit. The storage unit includes a value of each variable at three or more of the definition points, a first differential coefficient corresponding to the value of each variable of the three or more points, and a value of each variable. And an interpolation function determined by the first differential coefficient at a definition point that does not include a specific point among the three or more points, and the processing unit uses a second differential at the specific point by the interpolation function. A coefficient is calculated, and a third differential coefficient at the specific point is calculated by a weighted average of the second differential coefficient and the first differential coefficient at the specific point.
Moreover, a program for causing a computer to execute the numerical calculation method can be created, and further the program can be recorded on a recording medium.

  ADVANTAGE OF THE INVENTION According to this invention, the stability of the coupling between variables can be improved in the numerical solution of the partial differential equation using a collocated lattice.

Hereinafter, a numerical calculation method according to an embodiment of the present invention will be described with reference to the drawings as appropriate.
First, the configuration of a numerical calculation apparatus that performs a numerical calculation method will be described with reference to FIG. FIG. 1 is a configuration diagram of a numerical calculation device.

The numerical calculation device 10 is a computer device that includes an input unit 1, a storage unit (recording medium) 2, a memory 3, a processing unit 4, and an output unit 5, such as a personal computer, a workstation, or a supercomputer. Can be realized.
Note that although there is one processing unit 4 here, a plurality of processing units 4 may be provided to perform parallel computation.

The input unit 1 is means for a user of the numerical calculation device 10 (hereinafter simply referred to as “user”) to input various information, such as a keyboard and a mouse.
The storage unit 2 stores various types of information, such as a ROM (Read Only Memory) and a hard disk. Here, in order to perform a fluid dynamics numerical simulation, the storage unit 2 is a lattice point that is a definition point for performing a numerical calculation, initial conditions of physical quantities (variables, differential coefficients, etc.), a governing equation (a partial differential representing a natural phenomenon). Equations), interpolation functions, various programs, etc.

The memory 3 is temporary storage means used by the processing unit 4 and is, for example, a RAM (Random Access Memory).
The processing unit 4 performs arithmetic processing based on various types of information (programs, data, etc.) stored in the storage unit 2, and is, for example, a CPU (Central Processing Unit).
The output unit 5 outputs the result of the arithmetic processing by the processing unit 4 and is, for example, a display.

  Next, numerical simulation processing will be described with reference to FIG. 2 (see FIG. 1 as appropriate). FIG. 2 is a flowchart showing a rough flow of the numerical simulation process of fluid dynamics.

First, when the user performs an operation related to the start of a numerical simulation with the input unit 1, the processing unit 4 is based on information stored in the storage unit 2, and is a grid point that is a definition point for performing an operation, and an initial value of each physical quantity. Conditions and the like are set (step S1).
For example, the lattice points are set to 32 × 32 in the case of two dimensions and 16 × 16 × 8 in the case of three dimensions. As initial conditions for each physical quantity, variable values at the respective lattice points, spatial gradients (differential coefficients) of these values, and the like are set. Examples of variables include speed, pressure, temperature, density, and the like.

Next, the processing unit 4 sets, based on the information stored in the storage unit 2, 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).
As the governing equation, for example, an equation of motion of a compressible fluid can be cited, and details will be described later. Further, examples of the interpolation function include a cubic function, a rational function having a function in the denominator, and the details will be described later.

Subsequently, the processing 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, the physical quantity at each grid point at time (t + 1). A process of calculating a value is performed for t = 0, 1, 2,..., N-1 (steps S3 to S5).
And the process in step S4 in this is the point of this invention, The detail is demonstrated later with FIG.

  Next, the processing unit 4 outputs the value of each physical quantity calculated in steps S3 to S5 to the output unit 5 (step S6). Thus, the numerical simulation is completed.

Next, the process in step S4 of FIG. 2 will be described with reference to FIGS. 3 and 4 (see FIG. 1 as appropriate). FIG. 3 is a diagram showing the relationship between each grid point and a variable, first-order differential coefficient (hereinafter simply referred to as “differential coefficient”).
As shown in FIG. 3, three of the lattice points, j−1 points, j points (specific points), and j + 1 points are arranged at an interval Δx.

At the j−1 point, j point, and j + 1 point, the coordinates are −Δx, 0, Δx, and the variable (velocity) u values are u _j−1 , u _j , u _{j + 1} , and the variable (velocity) u. The differential coefficients are ux _{, j-1} , ux _{, j} , ux _{, j + 1} , the variable (pressure) p values are _pj-1 , _pj , _{pj + 1} , respectively, and the variable (pressure) ) The differential coefficients of p are p _{x, j−1} , p _{x, j} , p _{x, j + 1} , respectively.

  FIG. 4 is a flowchart showing details of the contents of the process in step S4 of FIG. The physical quantities handled in steps S41 to S46 are all at time t, and the values of the physical quantities calculated in step S47 are those at time t + 1.

  Here, the interpolation function U (x) of the variable (velocity) u is set as a cubic function that is one of the Hermite interpolation functions as shown in the following equation (7). The Hermite interpolation function is a generic name for interpolation functions determined from the values of variables on the definition points and their differential coefficients, and is suitable as the interpolation function here.

U (x) = C ₃ x ³ + C ₂ x ² + C ₁ x + C ₀ (7)
U _x (x), which is a differential expression of U (x), is expressed as formula (8).
U _x (x) = 3C ₃ x ² + 2C ₂ x + C ₁ (8)

The constraint conditions are U (Δx) = u _{j + 1} , U (−Δx) = u _j−1 , U _x (Δx) = ux _{, j + 1} , U _x (−Δx) = ux _{, j− 1,} these, the value of C ₀ -C ₃ is determined. Here, the feature is that the differential coefficient at point j is not used as a constraint, that is, U _x (0) = ux _{, j} is not used. If U _x (0) = u _{x, j} is used, C ₁ acquired in step S42 becomes equal to u _{x, j,} and the meaning of taking a weighted average in step S43 is lost.

Returning to FIG. 4, under the above conditions, the processing unit 4 first acquires the original differential coefficient ux _{, j} (first differential coefficient) (step S41).
Next, the processing unit 4 acquires a differential coefficient C ₁ (second differential coefficient) obtained from the interpolation function (Equation (7) and Equation (8)) (Step S42). That is, the differential coefficient at the point j is obtained as C ₁ by the following equation (9).
∂f / ∂x = U _x (0) = C ₁ (9)

Subsequently, the processing unit 4 calculates an approximate differential coefficient ∂u / ∂x | _j (third differential coefficient) by the following equation (10) (step S43). Note that ∂u / ∂x | _j represents a differential coefficient at j point obtained from an interpolation function or the like.

  The stability of the numerical simulation can be enhanced by substituting an appropriate value for α in the equation (10) according to the numerical simulation conditions.

  For example, it has been demonstrated that, under many conditions, the stability of numerical simulations can be increased by substituting α = 2/3. A demonstration example (application example) will be described later, but here, a mathematical basis considered as a reason why α = 2/3 is often suitable will be described.

First, the equation of Taylor expansion around j point with respect to the physical quantity U (x) is shown in the following equation (11).

Then, using the equation (11) under the constraint that U (x + Δx) = u _{j + 1} , U (x) = u _j , U (x−Δx) = u _j−1, the secondary center When discrete approximation using the difference is performed, Expression (12) is obtained, and Expression (13) is obtained by rearranging Expression (12).

Further, an equation for Taylor expansion around the j point with respect to the spatial gradient U _x (x) is shown in the following equation (14).

Next, under the constraint that U _x (x + Δx) = ux _{, j + 1} , U _x (x) = ux _{, j} , U _x (x−Δx) = ux _{, j−1} , By performing discrete approximation by forward difference and backward difference using Expression (14) and rearranging the expressions, a third-order differential term can be obtained as shown in the following Expression (15).

Then, by substituting equation (15) into equation (13), the following equation (16) can be obtained.

Further, when C ₁ in the equation (10) is calculated using the equations (7) and (8), the following equation (17) is obtained.
Further, when Expression (17) and α = 2/3 are substituted into Expression (10), the following Expression (18) is obtained.

As can be seen from the above description, the expression (18) matches the left two terms on the right side of the expression (16), and the quaternary accuracy (Δx in terms of the third and subsequent items from the left of the right side of the expression (16) Since the lowest order is 4), higher accuracy than the conventional calculation methods such as secondary accuracy and other calculation methods can be realized.
This is considered to be one of the mathematical grounds showing the advantage of substituting α = 2/3 in equation (10).

However, since the numerical simulation conditions vary widely, α = 2/3 is not always the best, and α may be set to a value other than 2/3 as necessary (2 / Substituting a value close to 3 or a value smaller than 2/3).
For example, if the wave wavelength of the object of numerical simulation (gas, liquid, electromagnetic wave, etc.) is relatively longer than the relative distance of the lattice, calculate by substituting a numerical value smaller than 2/3 for α. May be stabilized. This is because when the wavelength of the target wave is relatively long compared to the lattice interval (several times to several hundred times the lattice interval, etc.), when calculating the physical quantity of the j point, the adjacent lattice point (j−1 This is because the necessity of correction may be reduced depending on the physical quantity of point, j + 1 point, and the like.
That is, the feature of the present invention is to calculate an approximate differential coefficient from two values of an original differential coefficient and a differential coefficient obtained from an interpolation function, as shown in Expression (10), and α = 2 / 3 is one example.

Returning to FIG. 4, the processing unit 4 calculates the approximate differential coefficient ∂p / ∂x | _j with respect to the variable (pressure) p through steps S44 to S46, but steps S44 to S46 are steps S41 to S43. In this case, since u is replaced with p in FIG.

Next, the processing unit 4 calculates a governing equation using ∂u / ∂x | _j calculated in step S43 and ∂p / ∂x | _j calculated in step S46, and each physical quantity (at time t + 1) ( The values of speed, pressure and their spatial gradient are calculated (step S47).
For example, when the object of the numerical simulation is a compressible fluid, the governing equation is expressed by the following equation (19). Note that ρ is density.

If this equation (19) is integrated over time, the value of each physical quantity at time t + 1 can be calculated. For time integration, the Runge-Kutta method or the like may be used. Further, ∂u / ∂x | _j and ∂p / ∂x | if the _j [rho were calculated using alpha = 2/3 is, [rho are similarly calculated from the interpolation function for calculating the [rho You may make it use the value calculated | required by the weighted average of 2/3: 1/3 (2: 1) based on the value and the original value.
Further, instead of the interpolation function of Expression (7), a rational function having a function in the denominator, such as the following Expression (20), may be used as the interpolation function.

Although the numerator is a cubic function and the denominator is a linear function here, other functions may be used. Needless to say, this equation (20) can be applied not only to pressure but also to speed and other variables.
As described above, even when a rational function such as Expression (20) is used as the interpolation function, the method of calculating the approximate differential coefficient ∂u / ∂x | _j is changed from Expression (7) to Expression (20). Except for the replacement, the procedure is the same as that in steps S41 to S43 in the flowchart of FIG.

  In the early stage of numerical simulation, in the case of calculations related to blasts or shock waves with density ratios or pressure ratios of several hundreds or thousands of times, it is possible to improve the stability of the calculation by using a rational function such as equation (20). Can be increased.

Next, four application examples regarding the numerical calculation method by the numerical calculation apparatus 10 of the present embodiment will be described.
<Application example 1: Shallow water wave equation>
The dynamic process in the hydrostatic atmospheric model is described by the shallow water wave equation. When the shallow water wave equation is used for a numerical simulation, a highly accurate and stable coupling is required for the gravitational potential and the flow velocity.

  In the past, Temperton et al. Performed high-resolution calculations for the shallow water wave equation (the literature “An efficient two-time-level semi-Lagrangian semi-implicit integration scheme” (Temperton, C., and Stainforth, A , 1987, "Quart. J. Roy. Meteor. Soc. (113 edition) pages 1025-1039"). The calculation by Temperton et al. Is performed with a lattice number of 800 × 800.

In this case, the governing equation is a shallow water wave equation on the rotating sphere. This shallow water wave equation is expressed by the following equations (21) to (23) using Cartesian coordinates (x, y) on the polar stereo projection plane.
Where t is the time, u and v are the x-direction component and y-direction component of the flow velocity, φ is the geopotential height, m is the map factor, Ω is the angular velocity of the earth's rotation, and lat is the latitude. = U / m, V = v / m, S = m ² , f = 2Ωsin (lat) (Coriolis parameter).

  The calculation target area is a square area centered on the North Pole, and the reference latitude of the map factor is 60 degrees north latitude. As a boundary condition, a fixed wall condition is given near the equator. As an initial condition, a flow velocity distribution and a geopotential height distribution are given based on an analysis with respect to a pressure of 500 mb at 12:00 on February 28, 1984, and integration for two days is performed at 5-minute intervals.

  Under such conditions, FIG. 5A is a diagram showing an error of the geopotential height when calculated by a conventional calculation method that does not perform coupling between variables, and FIG. It is the figure showing the error of the geopotential height at the time of calculating by the numerical calculation method of the form, and the number of grid points is 100 × 100.

  In each figure, the center of the figure represents the North Pole, the curves 501 and 511 represent the equator, the curves 502 and 512, etc. represent the points where the geopotential heights obtained by the calculation of Temperton et al. Are equal. . Reference numerals 503 and 513 represent the relationship between the magnitude of the geopotential height error and the color.

  As can be seen by comparing these two figures, the calculation result obtained by the numerical calculation method of the present embodiment is much more accurate than the conventional calculation method, and is calculated by Temperton et al. And almost the same calculation results are obtained.

<Application example 2: Turbulent flow calculation>
For the calculation of turbulent flow, a spectral method in which each physical quantity is Fourier-expanded and calculation is performed in a wave number space is mainly used, and is known as the most accurate calculation method.
On the other hand, the numerical calculation method of this embodiment uses the governing equations shown in the following equations (24) to (26). Equations (24) and (25) are two-dimensional incompressible Navier-Stokes equations, and equation (26) is a continuous equation.
Note that u and v are the x direction component and y direction component of the flow velocity, x and y are coordinates, t is time, P is pressure, and Re is Reynolds number.

  Then, a predetermined turbulent flow profile (here, vorticity ω = ∂v / ∂x−∂u / ∂y) is given to a square region having the number of lattice points of 512 × 512. At this time, if the equations (24) to (26) are solved using the periodic boundary condition, a result as shown in FIG. 6 is obtained.

FIG. 6 is a diagram showing the relationship between dissipation and ETT (Eddy Turnover Time). Dissipation is a fluid statistic that represents the amount of kinetic energy converted to thermal energy per unit time, where Φ is the dissipation, ε is the entropy, and Re is the Reynolds number. = 2ε / Re. It should be noted that, if the vorticity and ω, is ε = ω ^2/2.
ETT is a dimensionless time, and is expressed as 1ETT = L / U _rms where L is the integral length in the initial stage and U _rms is the standard deviation of speed.

As shown in FIG. 6, the calculation result by the spectral method (curve 601: Spectral) and the calculation result by the numerical calculation method of the present embodiment (dots in the figure: IDO-SC) almost coincide.
Although the turbulent flow phenomenon is irregular in terms of time and space, according to the numerical calculation method of this embodiment, not only dissipation but also standard deviation of velocity, skewness, flatness, enstrophy, etc. As for other statistics, high-stability calculations similar to the spectral method can be realized.

<Application example 3: Two-dimensional cavity flow problem>
The problem of cavity flow is widely known as a benchmark for incompressible fluid calculations, and conventionally, the results of the multigrid method by Gear et al. Are considered to be the most accurate (reference: High-Resolution for incompressible flow using the Navier-Stokes equations and a multitigrid method "(Ghia, U., Ghia, KN, and Shin, CT, 1982," J. Comput. Phys. (48th edition) pages 387-411 ").

As for the problem of the cavity flow, since the coupling between the pressure and the flow velocity is important, the conventional IDO method cannot correctly evaluate the flow velocity. FIG. 7A is a diagram showing a result when calculation is performed by the conventional IDO method with the number of lattice points of 40 × 40.
Note that the governing equations used for the calculation by the conventional IDO method and the calculation by the numerical calculation method of the present embodiment are the above-described equations (24) to (26).

  In FIG. 7A, x and y represent coordinates. Further, v represents the flow velocity in the x direction corresponding to the x coordinate at a point where y = 0.5, and the dot (Ghia et al) in the vicinity of the curve 701 (IDO) is the calculation result of Gear et al. 256), a curve 701 is a calculation result (number of grid points 40 × 40) by a conventional calculation method. Further, u represents the flow velocity in the y direction corresponding to the y coordinate at a point where x = 0.5, and the dot (Ghia et al) in the vicinity of the curve 702 (IDO) is the calculation result of Gear et al. 256), a curve 702 is a calculation result (number of grid points 40 × 40) by a conventional calculation method. As can be seen from FIG. 7A, the accuracy of the conventional calculation results is not good.

On the other hand, FIG. 7B is a diagram showing a calculation result by the numerical calculation device of the present embodiment. In FIG. 7B, x and y represent coordinates.
Further, v represents the flow velocity in the x direction corresponding to the x coordinate at the point where y = 0.5, and the dot (Ghia et al) in the vicinity of the curve 711 (IDO-SC) is the calculation result (number of grid points). 256 × 256) and the curve 711 are the calculation results (number of grid points 40 × 40) by the numerical calculation method of the present embodiment. Furthermore, u represents the flow velocity in the y direction corresponding to the y coordinate at the point where x = 0.5, and the dot (Ghia et al) in the vicinity of the curve 712 (IDO-SC) indicates the calculation result (number of grid points). 256 × 256) and the curve 712 are the calculation results (number of grid points 40 × 40) by the numerical calculation method of the present embodiment.

  If FIG.7 (b) is seen, the calculation result by the numerical calculation method of this embodiment will be substantially in agreement with the calculation result of the gears of the number of lattice points 256x256, although it is very small number of lattice points of 40x40. It is clear that

<Application example 4: Shock wave problem>
In a shock wave problem in which a large pressure ratio or density ratio is given as an initial condition, in the calculation by the conventional IDO method or the like, the coupling of pressure, density and flow velocity is poor, and the calculation becomes unstable from the beginning and breaks down.

  Here, two fluids having different densities and pressures exist in a sealed tube with a diaphragm separated, and when the diaphragm is instantaneously removed, the high pressure fluid flows into the low pressure fluid side, and the low pressure fluid A case will be described in which a numerical simulation is performed on a one-dimensional shock tube problem in which a shock wave is generated due to rapid compression.

In this case, the governing equation is a one-dimensional compressible Navier-Stokes equation and a continuous equation. When expressed in a non-conservative form, the following equations (27) to (29) are obtained.
Here, t is time, x is coordinates, ρ is density, p is pressure, q is artificial viscosity, u is flow velocity, and e is internal energy.

FIG. 8A is a relationship diagram of pressure and coordinates when the pressure ratio is 10: 1, and FIG. 8B is a relationship diagram of pressure and coordinates when the pressure ratio is 10000: 1.
In FIG. 8A, the horizontal axis represents coordinates, and the vertical axis represents fluid pressure. Here, initially, a fluid having a pressure of 1 exists on the left side of x = 1.0, and a fluid of a pressure of 0.1 exists on the right side of x = 1.0. Under these conditions, x = 1 The pressure at each coordinate of the fluid at a certain time after a short time after the diaphragm provided at the point of 0.0 is instantaneously removed is shown.

  A curve 801 (Exact) is an analytical solution (exact solution) obtained logically, and a calculation result by the numerical calculation method of the present embodiment indicated by dots (IDO-SC) is almost in agreement with the curve 801. I understand that.

In FIG. 8B, the horizontal axis represents coordinates, and the vertical axis (logarithm display) represents fluid pressure (Pressure). Here, initially, a fluid having a pressure of 1 (10 ⁰ ) is present on the left side of x = 1.0, and a fluid having a pressure of 10 ⁻⁴ is present on the right side of x = 1.0. , X = 1.0, the pressure at each coordinate of the fluid at a certain time after a minute time has elapsed after the diaphragm provided at the point of 1.0 is instantaneously removed.

A curve 802 (Exact) is an analytical solution (exact solution) obtained logically, and a calculation result by the numerical calculation method of the present embodiment indicated by dots (IDO-SC) is almost in agreement with the curve 802. I understand that.
According to the numerical calculation method of the present embodiment, not only a general shock tube problem as in FIG. 8A, but also a shock tube problem with a larger pressure ratio as in FIG. 8B. Even for problems involving interactions with walls and walls, stable and highly accurate calculations can be performed. When the density ratio or pressure ratio is extremely large, the interpolation function is not a cubic function but a rational function having a function in the denominator, for example, a rational function having a linear function in the denominator and a cubic function in the numerator. It is desirable to use etc.

  As described above, according to the numerical calculation method of the present embodiment, stable coupling between variables can be realized in a multivariable partial differential equation using a collocated grid. In addition, a program for causing a computer to execute the numerical calculation method can be created, and the program can be stored (recorded) in the storage unit (recording medium) 2 (see FIG. 1).

This is the end of the description of the embodiments, but the aspects of the present invention are not limited to these.
For example, the present invention can be widely applied to general reproduction of natural phenomena using partial differential equations, such as numerical simulation related to electromagnetics and analysis of Rayleigh Taylor instability.
In addition, specific configurations such as the numerical calculation device and each flowchart can be appropriately changed without departing from the spirit of the present invention.

It is a block diagram of a numerical calculation apparatus. It is a flowchart showing the rough flow of the process of numerical simulation. It is a figure showing the relationship between each lattice point, a variable, and a differential coefficient. It is a flowchart showing the content of the process in step S4 of FIG. It is a figure showing the calculation result of the application example 1, (a) is a figure showing the error at the time of calculating with the conventional calculation method, (b) is the error at the time of calculating with the numerical calculation apparatus of this embodiment. FIG. It is a figure showing the calculation result of the application example 2, and is the figure which showed the relationship between dissipation and ETT. It is a figure showing the calculation result of the application example 3, (a) is the figure which showed the calculation result by the conventional IDO method, (b) is the figure which showed the calculation result by the numerical calculation apparatus of this embodiment. . It is a figure showing the calculation result of the application example 4, (a) is a relationship figure of the pressure in the case of 10: 1 pressure ratio, and a coordinate diagram, (b) is the relationship between the pressure in the case of a pressure ratio of 10000: 1. FIG.

Explanation of symbols

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

Claims (7)

To correct the value of one or more physical variables at each of a plurality of definition points in space, the governing equation including a partial differential equation representing a predetermined physical phenomenon, and the value of the variable used in the governing equation And an interpolation function in which a coefficient is determined based on the values of two or more variables other than its own defined point for each of the defined points, and is updated in a predetermined minute time unit. A numerical calculation method by a numerical calculation device,
The numerical calculation device includes a storage unit and a processing unit,
The storage unit stores the governing equation and the interpolation function,
When the processing unit updates the value of the variable for a specific point among the plurality of definition points,
Calculating a first derivative for the variable using the governing equation;
Calculating a second derivative for the variable using the interpolation function;
Calculating a third derivative for the variable by a weighted average of the first derivative and the second derivative ;
A numerical calculation method , wherein the value of the variable is updated to a value after the predetermined minute time by performing an operation based on the third differential coefficient and the governing equation . The numerical calculation method according to claim 1, wherein the plurality of definition points are collocated grid points in which all variables and their differential coefficients are defined on the same grid .   The numerical calculation method according to claim 1, wherein the interpolation function is a Hermite interpolation function.   The numerical calculation method according to claim 1, wherein the interpolation function is a rational function having a cubic function in a numerator. The weighted average according to claim 1 or claim 2, wherein the ratio of said first derivative and the previous SL second derivative is 1/3 to 2/3 or a value close to Numerical calculation method.   The program for making a computer perform the numerical calculation method as described in any one of Claims 1-5. A recording medium on which the program according to claim 6 is recorded.