JP6576303B2 - Simulation apparatus, simulation apparatus control method, and control program - Google Patents

Description

  The present invention relates to a simulation apparatus for simulating a flow of an incompressible fluid or a pseudo-compressible Newtonian fluid.

  As interest in the protection of ecosystems in recent years has increased and awareness of ensuring the safety of mankind has increased, the importance of research on the aquatic environment in rivers, lakes, and oceans has increased. In the aquatic environment, it is difficult to directly identify the causal factors involved in various ecosystems. Observe spatio-temporal changes such as chemical concentration, temperature, and water conditions, which are factors that have a large impact on the aquatic environment. Or the change is estimated in a test.

  In addition, since the water condition, that is, the flow of transporting chemical substances, has a great influence on the aqueous environment, there are cases where measurements are made at multiple points over several days to several years. However, since it requires a large amount of money and labor, it is difficult to say that the actual measurement data is sufficiently collected.

  On the other hand, based on numerical simulation of fluids, attempts have been made to supplement and predict the aqueous environment by adding various models.

  As a fluid numerical simulation, for example, Non-Patent Document 1 discloses a typical ocean model approximating hydrostatic pressure. Non-Patent Document 2 discloses a fluid analysis model.

George L. Mellor, "USERS GUIDE for A THREE-DIMENSIONAL PRIMITIVE EQUATION NUMERICAL OCEAN MODEL" 2002, [online], October 2002, Princeton University [Searched on December 1, 2015] <URL: (http: // www.ccpo.odu.edu/POMWEB/UG.10-2002.pdf> Toshio Kobayashi [Others] "Computational Fluid Dynamics Handbook" Maruzen, March 2003, pp.19-65 Author Takahiko Tanahashi "Basic Theory of CFD" IPC, pp.201-220, 1990 Hirt, C.W., Nichols, B.D., Romero, N.C. (1975): SOLA-A Numerical Solution Algorithm for Transient Fluid Flow, Los Alamos Scientific Laboratory, LA-5852, pp. 1-50

  However, the above-described conventional technology cannot solve the problem of numerical simulation in an aqueous environment such as a river, a lake, and the ocean. In an aquatic environment such as a river, a lake, and the ocean, the vertical scale of the calculation grid used for the numerical simulation is smaller on the order of 1/10 to 1/1000 compared to the horizontal scale. In the method of Non-Patent Document 1, the momentum in the vertical direction is not taken into consideration, and there is a risk that the accuracy of the solution is lowered for a flow in which this becomes a driving force. On the other hand, the method of Non-Patent Document 2 generally deviates from the stability limit related to the analysis method, so that the simulation calculation does not converge and a numerical solution may not be obtained.

  The present invention has been made in view of the above-described problems, and an object thereof is to realize a simulation apparatus and the like that can be executed with high accuracy without diverging numerical simulation operations in an aqueous environment.

  In order to solve the above problems, a simulation apparatus according to the present invention is a simulation apparatus for an incompressible fluid or a pseudo-compressible Newtonian fluid that moves in a three-dimensional space, and includes a calculation grid set in the three-dimensional space. At least one coordinate scale is at least twice as large as at least one other coordinate scale, a partial differential equation derived from a basic equation is decomposed corresponding to each coordinate scale, and further an orthogonal decomposition equation is used. Using the equation corresponding to the at least one coordinate scale, the long wavelength component calculation unit for calculating the pressure and the velocity of the three-direction component, and using the equation corresponding to the other coordinate scale, the pressure and the small coordinate scale are calculated. A short wavelength component calculation unit that calculates the velocity of the corresponding direction component, the long wavelength component calculation unit, and the short wavelength component calculation A combining unit for adding the result of calculation, a structure comprising a.

  In order to solve the above problems, a control method for a simulation apparatus according to the present invention is a control method for a simulation apparatus for an incompressible fluid or a pseudo-compressible Newtonian fluid that moves in a three-dimensional space. At least one coordinate scale of the calculation grid set in the original space is at least twice that of the other coordinate scale, and the partial differential equation derived from the basic equation is decomposed corresponding to each coordinate scale, and further orthogonal decomposition Using the equation corresponding to the at least one coordinate scale, the long wavelength component calculating step for calculating the pressure and the velocity of the three-way component, and the equation corresponding to the other coordinate scale. A short wavelength component calculating step for calculating a velocity of a direction component corresponding to the pressure and the other coordinate scale, and the long The length component calculating step, and the method comprising, a synthetic step adding the calculation result of the short wavelength component calculating step.

  According to the above configuration, when the equation of motion has the ellipticity of a partial differential equation, the partial differential equation is decomposed into equations corresponding to each coordinate scale, and a plurality of equations obtained by orthogonal decomposition are used. Therefore, the fluid can be calculated for each scale that has been decomposed without being diverged. That is, it is possible to execute numerical simulations in an aqueous environment with high accuracy without diverging.

  An incompressible fluid is a fluid whose density does not change depending on pressure. The pseudo-compressible Newtonian fluid is a fluid whose density can be changed by a physical factor other than pressure. For example, the density can be changed depending on temperature and chemical species concentration.

  In the simulation apparatus according to the present invention, it is preferable that the position of the liquid surface at the current time is the pressure boundary condition at the next time.

  In the simulation apparatus according to the present invention, the at least one coordinate scale may be larger than 10 times the other coordinate scale.

  According to said structure, it can apply to the calculation of an ocean calculation and a hydrostatic pressure approximation.

  The simulation apparatus according to each aspect of the present invention may be realized by a computer. In this case, the simulation apparatus is realized by a computer by operating the computer as each unit (software element) included in the simulation apparatus. A control program for the simulation apparatus and a computer-readable recording medium on which the control program is recorded also fall within the scope of the present invention.

  According to the present invention, when the equation of motion has the ellipticity of a partial differential equation, the partial differential equation is decomposed into equations corresponding to each coordinate scale, and a plurality of equations obtained by orthogonal decomposition are used. The fluid can be calculated for each scale that has been decomposed without diverging. That is, according to the present invention, it is possible to execute a numerical simulation operation in an aqueous environment with high accuracy without diverging, and it is possible to perform a calculation that is difficult to obtain practicality in terms of calculation load.

It is a block diagram which shows the principal part structure of a simulation apparatus. It is a figure which shows the grid example of staggered arrangement | positioning. It is a figure for demonstrating the pressure boundary condition in a free surface. It is a flowchart which shows the flow of simulation in a simulation apparatus. It is a figure which shows the rectangular three-dimensional area | region made into calculation object. It is a figure which shows calculation conditions. It is a figure which shows a calculation result. It is a figure for demonstrating the example of a calculation at the time of using for the ocean, and is a figure which shows the depth information of the Seto Inland Sea. It is a figure for demonstrating the example of a calculation at the time of using for the ocean, and is a figure which shows the relationship between the time and speed of the tidal current of the Seto Inland Sea. It is a figure which shows the result of the example of calculation at the time of using for the ocean.

  The simulation apparatus 1 according to the present embodiment can perform a simulation by storing a predetermined formula in advance.

  First, a method for deriving equations used in the simulation method executed by the simulation apparatus 1 according to the present embodiment will be described.

  In the following, a method for deriving an expression when the vertical coordinate scale of the calculation grid set in the three-dimensional space is twice or more the two horizontal coordinate scales will be described.

  In the numerical calculation of the incompressible fluid in the three-dimensional orthogonal coordinate system, the following formula (1) (Courant condition) needs to be satisfied in order to ensure stability.

Here, x and y indicate horizontal coordinates, and z indicates vertical coordinates. Δx, Δy, and Δz indicate coordinate scales in the respective directions, and Δt indicates a time step size. U, v, and w indicate the speed in each direction. The subscript “max” represents the maximum value in the calculation area.

Here, in an aqueous environment such as a river, a lake, and the ocean, the vertical scale of the calculation grid used for the numerical simulation is smaller on the order of 1/10 to 1/1000 compared to the horizontal scale. Based on this point, when Δx≈Δy >> Δz, Δt substantially depends on w _max / Δz from equation (1). In the case of Δx≈Δy≈Δz, a numerical solution can be obtained if Expression (1) is satisfied, but an enormous number of calculation grids are required in the horizontal direction.

  Therefore, the following hydrostatic pressure approximation model has been proposed. In order to avoid the problem that the calculation of the simulation does not converge and diverges due to the limitations of the calculation grid scale as described above, first, the following expression that approximates the hydrostatic pressure is used.

Here, p is pressure, p ₀ is atmospheric pressure, h is water surface height, ρ is density, and g is gravitational acceleration. By using this approximate expression, the Courant condition is relaxed to the following expression.

  However, in the above equation (3), the dynamic behavior in the vertical direction cannot be calculated. In recent simulations, accuracy at a resolution of several hundred meters or less is required in the land near the sea and in the bay, and it is essential to consider the dependency of horizontal velocity on the water depth and the influence of vertical flow. Therefore, the above equation (3) is not necessarily valid.

  Therefore, in the present embodiment, equations are derived by the following method in order to avoid the problem that the calculation of the simulation does not converge and diverges due to the limitations of the calculation grid scale as described above. Specifically, first, a quasi-dynamic equation of motion with respect to the vertical direction is derived, and then a discretization method is applied for application to the Projection method. It is then expanded to incorporate a free surface.

(Derivation of formula)
As the fluid to be simulated, a pseudo-compressible Newtonian fluid having no phase change is assumed.

  First, a mass conservation equation (4) and momentum conservation equations (5), (6), and (7) are given as basic equations to the Cartesian three-dimensional orthogonal coordinate system as follows.

Here, t is time, D / Dt is substance differentiation, and μ is viscosity. Taking into account Δx≈Δy >> Δz, the equation (7) in the z direction can be approximated as the following equation (8).

  Next, considering the wavelength in the vertical direction, the variable w and the variable p are decomposed into two, a long wavelength component and a short wavelength component, as follows. In addition, in this specification, a letter with a hat (^) added to an alphabetic character shown in the formula is indicated with a hat next to the alphabetic character. The same applies to other characters.

Here, w ^ among u, v, and w in the above formulas (4) to (8) is a wavelength component that can be captured on a scale where the wavelength is Δx≈Δy >> Δz (hereinafter, a long wavelength component and ) Of each direction, and p ^ of p means the pressure of the long wavelength component. In addition, w ′ of w means the velocity in the z direction of a wavelength component (hereinafter also referred to as a short wavelength component) that can be captured on a Δz scale, and p ′ of p is a short wavelength. It means the pressure of the component. When the above equations (4) to (6) and (8) are rewritten using the equations (9) and (10), the following equations (11), (12), and (13) ) And equation (14) are obtained respectively.

  Next, we applied the Projection method (see Takahiko Tanahashi's "Basic Theory of CFD" IPC, pp. 201-220, 1990) as the numerical method, and the Euler explicit method (Toshio Kobayashi [others] edited "Computational Fluid Dynamics Handbook". Maruzen, March 2003, see pp. 19-65), the above equation (12), equation (13), and equation (14) are discretized in terms of time, and the following equation (15) Equation (16) and Equation (17) can be obtained.

Here, the superscript “n + 1” indicates a value at the next time (n + 1) Δt, and G _x , G _y , and G _z are variables summarizing known values obtained from the current time n. Indicates.

(1. Speed predictor)
In the above equations (15), (16), and (17), p ^ ^{n + 1} is unknown, so the speed predictor is excluded except for the pressure term.

Calculate Therefore, each of the speed predictors can be described as the following Expression (18), Expression (19), and Expression (20).

  First, derivation is advanced for u and v. By substituting the above equations (18) and (19) into the above equations (15) and (16), the following equations (21) and (22) are obtained.

  From this, the following equations (23) and (24) are obtained, respectively.

  On the other hand, for w, the following equation (25) is obtained by substituting the above equation (20) into the above equation (17).

  From this, the following equation (26) is obtained.

  Next, the following formula (27) is obtained by discretizing the formula (11) with respect to time.

(2. Expression of long wavelength component of pressure)
In consideration of the wavelength component that can be captured on the scale of Δx≈Δy in the above equation (27), w ′ = p ′ = 0 and orthogonal decomposition, the following Poisson equation (28) for p ^ ^{n + 1} (pressure length) Wavelength component formula) can be obtained.

(3. Formula of long wavelength component of velocity)
Similarly, when w ′ = p ′ = 0 in the above formula (21), formula (22), and formula (25) in consideration of the wavelength component that can be captured on the scale of Δx≈Δy, As shown in Expression (29) to Expression (31), a time update expression (expression of a long wavelength component of velocity) of u, v, and w ^ can be obtained.

(4. Formula of short wavelength component of pressure)
On the other hand, when the above equation (27) is solved by orthogonal decomposition with respect to the wavelength component that can be captured at the scale of Δz, the following Poisson equation (32) (expression of the short wavelength component of pressure) for p ′ ^{n + 1} can be obtained. .

(5. Formula of short wavelength component of velocity)
Similarly, when the above equation (26) is solved with respect to the wavelength component that can be captured on the scale of Δz, as shown in the following equation (33), the time update equation for w ′ (the equation for the short wavelength component of the velocity) is obtained. Can get.

Here, the calculation of w ′ requires the restriction of Δt due to the stability problem. Therefore, instead of Delta] t, for example, the following equation (34), and preferably used Delta] t _z calculated by the equation (35).
Δt _z = Δt / E _z (34)

Here, INT means rounding down to the integer.

[Example of numerical solution of Poisson equation]
(1. Long wavelength component of pressure)
A method for numerically solving the Poisson equation (28), which is an expression of the long wavelength component of pressure, will be described.

  Following the usual computational fluid dynamics method, the variable uses staggered arrangement. As shown in FIG. 2, the scalar is arranged at the lattice center and the vector is arranged at the lattice boundary surface. FIG. 2 is a diagram illustrating an example of a lattice with a staggered arrangement. In FIG. 2, it is drawn in two dimensions for easy viewing. This suppresses non-physical vibrations of variables. When the lattice positions of the x, y, and z coordinate axes are set to i, j, and k, respectively, the first and second terms (= Sijk) of the Poisson equation (28) are expressed by the following formula (36). It can be discretized as follows.

  Similarly, the third term (pressure term) of the Poisson equation (28) can be discretized as shown by the following equation (37). In the following, superscripts representing time are omitted for easy description.

  By substituting the above equations (36) and (37) into the above equation (28), an operation relating to p can be performed. The calculation method is roughly divided into the direct method and the iterative method. Here, the Gauss-Seidel method, which is an example of the iterative method (Toshio Kobayashi [etc.], “Computational Fluid Dynamics Handbook” Maruzen, March 2003 , P.19-65), the following equation (38) is derived.

Here, C _C , C _R , C _L , C _T , C _B , C _F , and C _O are coefficients of p having the same subscript included in the formula (37). By using the above equation (37), a loop for i, j, k is assembled, and calculation is repeated until convergence, thereby obtaining a convergence value for the long wavelength component p ^ of pressure.

(2. Formula of short wavelength component of pressure)
Similarly, the Poisson equation (32), which is an equation for the short wavelength component of pressure, can also be discretized as shown by the following equation (39).

Similar to the long wavelength component p ^ of the pressure, the convergence value regarding the short wavelength component p 'of the pressure can be obtained by repeatedly calculating until convergence using the above formula (39).

(Pressure boundary condition)
Next, the pressure boundary condition considering the position of the liquid level will be described.

  When using a calculation grid of Δx≈Δy >> Δz as in the present embodiment, it is preferable to consider the position of the liquid surface because instability of the simulation and numerical error can be prevented. By calculating the free surface as a monovalent function in the horizontal direction, the position of the liquid surface can be taken into account. Examples of a method for calculating the free surface as a monovalent function in the horizontal direction include a SOLA-Surf (SOLution Algorithm-Surf) method.

  Therefore, in this embodiment, the SOLA-Surf method (Hirt, CW, Nichols, BD, Romero, NC (1975): SOLA-A Numerical Solution Algorithm for Transient Fluid Flow, Los Alamos Scientific Laboratory, LA-5852, pp. 1 -50), and calculate the free surface as a monovalent function h (the position of the liquid level from the seabed) in the horizontal direction. h can be obtained using the following equation (40). Any technique that can reproduce a free surface is acceptable. VOF (Volume of Fluid) method, ALE (Arbitrary Lagrangian-Eulerian) method, BFC (Boundary Fitted Coordinate) method (all of which are edited by Toshio Kobayashi [others] A method such as "Handbook" Maruzen, March 2003, pp. 19-65) is also applicable.

  In order to take into account the effect of the liquid level, it is necessary to impose pressure boundary conditions. FIG. 3 is a diagram for explaining the pressure boundary condition on the free surface.

The lattice positions in the x and z directions are numbered i and k, respectively. The pressure at the position (i, k) is expressed as p _ik . In accordance with the staggered arrangement, p _ik is defined by the lattice center position.

As an example, consider the case where there is a liquid level between the positions of the p _{ik + 1} of p _ik. The ratio of the distance between the liquid surface position and the position of p _{ik to} the cell width Δz is η. The following equation (41) can be obtained by giving the condition that the pressure on the atmosphere side becomes _patm at the liquid surface position and considering the linear interpolation relationship.

The pressure obtained by the above equation (41) is given as the boundary condition of the pressure equation (38) and the equation (39).

  In the Projection method, the pressure is directly and repeatedly calculated. Therefore, in addition to directly reflecting the conditions in consideration of the position of the liquid level, the solution method of the Poisson equation can be freely selected. Therefore, by adopting the Projection method as the numerical solution method as in this embodiment, it becomes easy to incorporate the position of the liquid level into the pressure boundary condition. Other numerical solutions include the SOLA method, SMAC (Simplified Maker and Cell) method, and the like. In particular, the SOLA method is a method of repeating speed and pressure at the same time, so that the position of the liquid level is easily incorporated as a pressure boundary condition.

  As will be described later, the simulation apparatus 1 according to the present embodiment can perform a simulation by storing the formula calculated above in advance.

  The apparatus configuration and processing flow of the simulation apparatus 1 according to this embodiment will be described in detail below.

[Configuration of Simulation Device 1]
First, with reference to FIG. 1, the principal part structure of the simulation apparatus 1 is demonstrated. FIG. 1 is a block diagram showing a main configuration of the simulation apparatus 1. As shown in FIG. 1, an input device 2 and a display device 3 are connected to the simulation device 1.

  The input device 2 is a device for operating the simulation device 1. The input device 2 only needs to accept an input operation for operating the simulation device 1 and transmit operation data indicating the content of the accepted input operation to the simulation device 1. For example, a mouse, a keyboard, or the like is applied. You can also.

  The display device 3 is a device that displays an image using display data output from the simulation device 1. The display device 3 may be any device that can display an image, and a liquid crystal display device, an organic EL (Electro Luminescence) display device, or the like can also be applied. On the display device 3, an operation screen for operating the simulation device 1, a simulation result, and the like are displayed.

  The simulation apparatus 1 includes a control unit 10 and a storage unit 20.

  The control unit 10 controls the operation of the simulation apparatus 1 in an integrated manner. The control unit 10 performs the above-described control by reading a program stored in the storage unit 20 into, for example, a RAM (Random Access Memory) and executing the program.

  The storage unit 20 stores data, programs, and the like necessary for the simulation apparatus 1 to operate. In addition, intermediate results calculated in the course of simulation, final results of simulation, and the like are also stored in the storage unit 20. Here, in the storage unit 20, parameters (density, viscosity, heat capacity, thermal conductivity, etc.), initial values, boundary conditions (pressure on the wall, tide flow, Etc.), variables, calculation instruction order, etc. are registered in advance and are read at the start of simulation or at the time of reference.

  The control unit 10 includes a simulation calculation unit 30, a simulation control unit 50, and a display control unit 60. The simulation calculation unit 30 includes a speed predictor calculation unit 31, a pressure long wavelength component calculation unit (long wavelength component calculation). Part) 32, pressure short wavelength component calculation part (short wavelength component calculation part) 33, pressure synthesis part (synthesis part) 34, speed long wavelength component calculation part (long wavelength component calculation part) 35, speed short wavelength component calculation part ( A short wavelength component calculation unit) 36, a velocity synthesis unit (synthesis unit) 37, and a pressure boundary condition calculation unit 38 are included.

  The velocity predictor calculation unit 31 calculates velocity predictors in the x, y, and z directions for each calculation grid position (i, j, k) on the x, y, and z coordinate axes. Specifically, for each calculation grid, calculation is performed using the above-described Expression (18), Expression (19), and Expression (20).

  The pressure long wavelength component calculation unit 32 calculates a long wavelength component of pressure for each calculation grid position (i, j, k) of the x, y, z coordinate axes. Specifically, for each calculation grid, the above-described equation (38) is used to repeatedly calculate until convergence, and the convergence value is notified to the pressure synthesis unit 34 as the calculation result of the long wavelength component of pressure.

  The pressure short wavelength component calculator 33 calculates a short wavelength component of pressure for each calculation grid position (i, j, k) of the x, y, z coordinate axes. Specifically, for each calculation grid, the above-described equation (39) is used to repeatedly calculate until convergence, and the convergence value is notified to the pressure synthesis unit 34 as the calculation result of the short wavelength component of pressure.

  The pressure synthesizer 34 calculates the long wavelength component of the pressure calculated by the pressure long wavelength component calculator 32 and the pressure short wavelength component calculator 33 for each calculation grid position (i, j, k) of the x, y, z coordinate axes. Is added to the short wavelength component of the calculated pressure, and the total value is notified to the simulation control unit 50 as the pressure of the calculation grid.

  The velocity long wavelength component calculator 35 calculates a velocity long wavelength component for each calculation grid position (i, j, k) of the x, y, z coordinate axes. Specifically, the long wavelength component of the velocity in the x, y, and z directions is calculated for each calculation grid using the above-described Expression (29), Expression (30), and Expression (31). Of the calculation results, the long wavelength component of the velocity in the x and y directions is notified to the simulation control unit 50 as the velocity in the x and y directions of the calculation grid. Further, the long wavelength component of the velocity in the z direction is notified to the velocity synthesis unit 37.

  The velocity short wavelength component calculator 36 calculates a short wavelength component of velocity in the z direction for each calculation grid position (i, j, k) of the x, y, z coordinate axes. Specifically, for each calculation grid, the short wavelength component of the z-direction velocity is calculated using the above-described equation (33), and the calculation result is notified to the velocity synthesis unit 37.

  The velocity synthesizer 37 has a long wavelength component of velocity in the z direction and a short wavelength component of velocity calculated by the velocity long wavelength component calculator 35 for each calculation grid position (i, j, k) of the x, y, z coordinate axes. The short wavelength component of the velocity in the z direction calculated by the calculation unit 36 is added, and the total value is notified to the simulation control unit 50 as the velocity in the z direction of the calculation grid.

  For each current time t, the pressure boundary condition calculation unit 38 calculates the liquid level position h at the current time as a boundary condition used for calculation of the simulation at the next time t + Δt using the above equation (40). Further, the pressure boundary condition is calculated using the above equation (41).

  The simulation control unit 50 uses the simulation calculation unit 30 to calculate each time and each calculation grid position, thereby executing a fluid simulation and transmitting the execution result to the display control unit 60. The simulation control unit 50 may store the simulation execution result in the storage unit 20.

  More specifically, when the simulation calculation unit 30 finishes calculating the pressure and velocity at a certain calculation grid position, the simulation control unit 50 updates the target calculation grid position, and the pressure and velocity at the updated calculation grid position. The simulation calculation unit 30 is caused to execute the speed calculation. This is repeatedly executed by the simulation calculation unit 30 until the calculation is completed for all calculation grid positions (all combinations of i, j, and k). Next, the simulation control unit 50 updates the time, and causes the simulation calculation unit 30 to execute the above-described calculation of pressures and velocities at all calculation grid positions. The simulation control unit 50 repeats the above processing while updating the time until the simulation end time.

  The display control unit 60 displays the simulation result executed by the simulation control unit 50 on the display device 3. Further, the display control unit 60 reads display data to be displayed from the storage unit 20 and causes the display device 3 to display it. That is, the display control unit 60 causes the display device 3 to display the simulation result of the above-described simulation device 1, reads operation screen data for operating the simulation device 1 from the storage unit 20, and causes the display device 3 to display the operation screen data. .

[Processing flow in the simulation apparatus 1]
Next, the flow of simulation in the simulation apparatus 1 will be described with reference to FIG. FIG. 4 is a flowchart showing the flow of simulation in the simulation apparatus 1.

  As shown in FIG. 4, first, the speed predictor calculation unit 31 calculates a speed predictor using the above equations (18), (19), and (20) (S101).

Next, the pressure long wavelength component calculator 32 follows the Gauss-Seidel method to calculate the long wavelength component p ^ ^{n + 1} of the pressure using the above equation (38) (S102, long wavelength component calculation step). In parallel with this, the pressure short wavelength component calculation unit 33 calculates the short wavelength component p ′ ^{n + 1} of the pressure using the above equation (39) (S103, short wavelength component calculation step). Here, the pressure obtained by the above equation (41) is given as a boundary condition of the above equations (38) and (39). Next, the speed long wavelength component calculation unit 35 substitutes the long wavelength component p ^ ^{n + 1} of the pressure calculated in step S102 into the above formulas (29), (30), and (31), so that the speed length Wavelength components u ^{n + 1} , v ^{n + 1} , and w ^ ^{n + 1} are calculated (S105, long wavelength component calculation step). In parallel with this, the velocity short wavelength component calculation unit 36 substitutes the short wavelength component p ′ ^{n + 1} of the pressure calculated in step S103 into the above equation (33) to calculate the short wavelength component w ′ ^{n + 1} of the velocity. (S106, short wavelength component calculation step).

Next, the pressure synthesizer 34 calculates p from the long wavelength component p ^ ^{n + 1} calculated by the pressure long wavelength component calculator 32 and the short wavelength component p ′ ^{n + 1} calculated by the pressure short wavelength component calculator 33. ^ ^{N + 1} + p′n ^{+ 1} is calculated (S104, composition step).

The speed synthesizing unit 37, among the speed-length long wavelength components of the wavelength velocity component calculating unit 35 has calculated ^{u n + 1, v n +} 1, w ^ n + 1, the long wavelength component ^{w ^ n + 1} in the z-direction, velocity short wavelength component calculating From the short wavelength component w ′ ^{n + 1} of the velocity in the z direction calculated by the unit 36, w ^ ^{n + 1} + w′n ^{+ 1} is calculated (S107, synthesis step).

  Then, the simulation control unit 50 determines whether or not the calculation in steps S101 to S107 has been completed for all the calculation grid positions (S108). If not completed (NO in S108), the simulation grid position is updated. The process returns to step S101.

  On the other hand, when it is determined in step S108 that the calculation in steps S101 to S107 has been completed for all the calculation grid positions (YES in S108), the simulation control unit 50 initializes the positions (S109) and ends the simulation. It is determined whether or not the calculation has been completed up to the time (S110). When the calculation is not completed until the end time (NO in S110), the simulation control unit 50 causes the pressure boundary condition calculation unit 38 to calculate the liquid level position h at the current time (S112) and updates the time (S113). ). Then, the process returns to step S101.

  On the other hand, if it is determined in step S110 that the calculation has been completed up to the simulation end time (YES in S110), the simulation calculation ends.

  The current time is the time when steps S101 to S108 and S111 are executed. Then, calculation has been completed for all positions in step S108 (YES in S108), and when the positions are initialized in step S109, the process proceeds to the next time (next time when the simulation is executed).

  Here, in FIG. 4, it has been described that the operations of steps S102 and S103 are executed in parallel. That is, steps S102 and S103 are steps that can be independently calculated without dependency on each other after the speed predictor is calculated in step S101. Of course, steps S102 and S103 can be calculated one step at a time, and in this case, any step may be executed.

  As described above, the simulation method in the simulation apparatus 1 is an approach in which a long wavelength equation and a short wavelength equation are separately prepared and calculated. In addition, since the above equation (32) (that is, the above equation (39)) is limited to a component having a short scale, as shown in an example (calculation example) described later, calculation necessary for convergence is performed. The number of times does not increase extremely.

  Here, the size of the calculation grid used in the simulation apparatus 1 is such that the vertical scale (Δz) is smaller than ½ compared to the horizontal scale (Δx = Δy). However, the simulation can be executed even if the vertical scale (Δz) is approximately the same as the horizontal scale (Δx = Δy). The relationship between the lattice size ratio and the effect will be described later.

  Since the simulation apparatus 1 is based on the assumption that Δx≈Δy >> Δz, a third-order accuracy CIP (Constrained Interpolation Profile Scheme) method is applied to the calculation of the convection terms of u and v, and the spatial resolution is reduced. Improvements can be made. Further, the first-order upwind difference method can be applied to the calculation of the convection term of w to further improve the stability. Furthermore, the Poisson equation having the highest calculation load among the calculation steps can be accelerated by parallel calculation using the Red-Black SOR method.

〔Example〕
Next, with reference to FIGS. 5 to 10, the results of the simulation performed by the above-described method will be described.

[Influence of aspect ratio of calculation grid]
First, with reference to FIGS. 5 to 7, it will be demonstrated that, even if the aspect ratio of the calculation grid is changed, the calculation result can be obtained without divergence and the stability can be maintained by the method according to the present embodiment. FIG. 5 is a diagram showing a rectangular three-dimensional area as a calculation target. FIG. 6 is a diagram showing calculation conditions. Moreover, FIG. 7 is a figure which shows a calculation result.

Here, as described above, the rectangular three-dimensional region shown in FIG. The rectangular three-dimensional area shown in FIG. 5 has L _x = L _y = 20 m and L _z = 5 m, and the aspect ratio (L _x / L _z ) of the calculation area is 4. The calculation grid size in the horizontal direction is set to Δx = Δy = 1 m. The inlet boundary and the outlet boundary (both surfaces of 2 m × 2.5 m) are arranged so as to be diagonally opposite to each other so that the flow causes a three-dimensional motion, and the rest is a slip wall. The inlet and outlet boundaries had the same normal as the yz plane and the velocity u was given by u = 1.0 m / s.

  With this setting, the entrance boundary has momentum only in the x direction, but the momentum is also transmitted to the y direction and the z direction via pressure.

  As an initial condition, the omnidirectional speed was set to 0 m / s. In addition, the time interval Δt was set to 0.1 second, and transient calculation was performed up to 400 seconds corresponding to the residence time. Further, as shown in FIG. 6, the vertical calculation lattice size Δz was changed to three types of 1 m, 0.1 m, and 0.01 m, and the influence on the velocity field was compared.

  The calculation results are shown in FIG. (A) of FIG. 7 shows the result of execution number 1, (b) shows the result of execution number 2, and (c) shows the result of execution number 3. Further, FIG. 7 shows a velocity field on the xy section at a position of z = 2.5 m after 400 seconds from the start of calculation. Referring to FIG. 7, it can be seen that in all the execution results, the tendency that the direction of the speed changes via the pressure near the wall surface is consistent. The time required for calculation using a personal computer (CPU: Intel Xeon 3.4Ghz, Memory: 32 GB) was about 70 hours in the execution number 3 with the highest load. In addition, when using the conventional calculation method (the method described in Toshio Kobayashi [others] “Numerical Fluid Mechanics Handbook” Maruzen, March 2003, pp. 19-65), time step 100 except for execution number 1 It diverged below and the calculation result could not be obtained.

[Calculation example when used in the ocean]
Next, a calculation example when the present embodiment is used in the ocean will be described with reference to FIGS. FIG. 8 is a diagram for explaining a calculation example when used in the ocean, and is a diagram illustrating depth information of the Seto Inland Sea. FIG. 9 is a diagram for explaining a calculation example when used in the ocean, and is a diagram showing a relationship between the tidal current time and speed in the Seto Inland Sea. FIG. 10 is a diagram illustrating a result of a calculation example when used in the ocean.

  In this calculation example, the tidal current of the Seto Inland Sea was calculated as an example of the ocean. The Seto Inland Sea that touches Kagawa Prefecture, Ehime Prefecture, and Hiroshima Prefecture shown in FIG. The seafloor shape was input from the water depth data of the Japan Marine Data Center of the Japan Coast Guard.

  The calculation grid has a size of Δx = Δy = 500 m and Δz = 5 m, and the total number of grids is 255000. The eastern end of the Seto Inland Sea was set as the entrance boundary and the west end as the exit boundary. Time average value 4.97m / s calculated based on the observation period (June 2, 1992 4:00 to June 17, 3:00 3:00) data shown in Fig. 9 as the inlet flow velocity which is a boundary condition. Gave.

  The calculation results are shown in FIG. As shown in FIG. 8, the calculation target area includes many islands including the archipelago and has a water depth of 50 m or less. However, in this embodiment, as shown in FIG. It can be seen that can be calculated appropriately.

[Summary]
As described above, in this embodiment, a quasi-dynamic equation of motion with respect to the vertical direction is derived and discretized for application to the Projection method. Furthermore, the calculation method is expanded by incorporating the free surface into the boundary condition. As a result of the calculation using the calculation grid of Δx≈Δy >> Δz, the following conclusion is obtained.
(1) Even when the aspect ratio of the lattice size is away from 1, calculation is possible in this embodiment.
(2) For example, a complex tidal current can be calculated in the present embodiment even in a sea area including many islands including several islands including the archipelago such as the Seto Inland Sea.
(3) The relationship between the lattice size ratio and the effect of the present embodiment is specifically shown as follows. When the lattice size ratio (Δx / Δz) is “100 or more”, there is a possibility that the calculation does not converge with the conventional method, but in the present embodiment, the calculation is appropriately executed and the simulation is performed as described above. be able to. Note that the method according to the present embodiment is advantageous because the lattice size ratio (Δx / Δz) is likely to be 100 or more in many ocean calculations and calculations with hydrostatic approximation.

  In addition, this embodiment does not limit the scope of the present invention, and various modifications are possible within the scope of the present invention. For example, the present embodiment can be configured as follows.

  A calculation method according to the present invention is a calculation method for an incompressible fluid or a pseudo-compressible Newtonian fluid that moves in a three-dimensional space, and when the spatial region is divided so as to have a required resolution, the cubic is obtained. In the case where at least one coordinate scale of the original region is intended for an element that is more than twice the other coordinate scale, preferably larger than 10 squared, and the equation of motion has ellipticity of a partial differential equation, Even if the partial differential equation is decomposed into equations corresponding to each coordinate scale, and by using multiple equations obtained for pressure by further orthogonal decomposition, it is possible to calculate for each decomposed scale without causing divergence Good. The pseudo-compressible Newtonian fluid refers to a fluid whose fluid density can be changed by physical factors other than pressure. For example, the density can be changed depending on temperature and chemical species concentration.

  Furthermore, the calculation method according to the present invention may be a method of calculating a stable hydrology in a state having boundary conditions relating to a pressure such as a liquid level.

  In the present embodiment, the case where Δx≈Δy >> Δz has been described. However, the present invention only requires that at least one coordinate scale is at least twice as large as at least one other coordinate scale. For example, Δx >> Δy≈Δz, Δx << Δy≈Δz, Δx >> Δy >> It can also be applied to the case of Δz. In other words, if the scale of the grid to be decomposed in the coordinate direction is different, if the equation is decomposed according to each scale, the calculation of the numerical simulation in the water environment can be executed with high accuracy without diverging, and the calculation with improved stability Is possible.

  The present invention is also applicable to orthogonal systems such as polar coordinates. Furthermore, by applying to coordinate transformation methods such as the BFC (Boundary Fitted Coordinate) method, the non-Euclidean coordinate system (Toshio Kobayashi [etc.] “Computational Fluid Dynamics Handbook” Maruzen, March 2003, pp.19 (See -65 etc.).

[Example of software implementation]
The control blocks (in particular, the simulation calculation unit 30, the simulation control unit 50, and the display control unit 60) of the simulation apparatus 1 may be realized by a logic circuit (hardware) formed in an integrated circuit (IC chip) or the like. However, it may be realized by software using a CPU (Central Processing Unit).

  In the latter case, the simulation apparatus 1 includes a CPU that executes instructions of a program that is software that implements each function, a ROM (Read Only Memory) in which the program and various data are recorded so as to be readable by a computer (or CPU), or A storage device (these are referred to as “recording media”), a RAM (Random Access Memory) that expands the program, and the like are provided. And the objective of this invention is achieved when a computer (or CPU) reads the said program from the said recording medium and runs it. As the recording medium, a “non-temporary tangible medium” such as a tape, a disk, a card, a semiconductor memory, a programmable logic circuit, or the like can be used. The program may be supplied to the computer via an arbitrary transmission medium (such as a communication network or a broadcast wave) that can transmit the program. The present invention can also be realized in the form of a data signal embedded in a carrier wave in which the program is embodied by electronic transmission.

  The present invention is not limited to the above-described embodiments, and various modifications are possible within the scope shown in the claims, and embodiments obtained by appropriately combining technical means disclosed in different embodiments. Is also included in the technical scope of the present invention. Furthermore, a new technical feature can be formed by combining the technical means disclosed in each embodiment.

DESCRIPTION OF SYMBOLS 1 Simulation apparatus 2 Input apparatus 3 Display apparatus 10 Control part 20 Memory | storage part 30 Simulation calculation part 31 Speed predictor calculation part 32 Pressure long wavelength component calculation part (long wavelength component calculation part)
33 Pressure short wavelength component calculation unit (short wavelength component calculation unit)
34 Pressure synthesis unit (synthesis unit)
35 Velocity Long Wavelength Component Calculation Unit (Long Wavelength Component Calculation Unit)
36 Speed short wavelength component calculation unit (short wavelength component calculation unit)
37 Speed synthesis unit (synthesis unit)
38 Pressure boundary condition calculation unit 50 Simulation control unit 60 Display control unit

Claims (5)

A simulation device for incompressible fluid or pseudo-compressible Newtonian fluid moving in three-dimensional space,
At least one coordinate scale of the calculation grid set in the three-dimensional space is at least twice that of the other coordinate scale, and the partial differential equation derived from the basic equation is decomposed corresponding to each coordinate scale, and further orthogonal Using the decomposed equations,
A long wavelength component calculator that calculates the pressure and velocity of the three-way component using an equation corresponding to the at least one coordinate scale;
Using an equation corresponding to the other coordinate scale, a short wavelength component calculation unit that calculates the pressure and the velocity of the direction component corresponding to the other coordinate scale,
A simulation apparatus comprising: the long wavelength component calculation unit; and a synthesis unit that adds the calculation results of the short wavelength component calculation unit.   The simulation apparatus according to claim 1, wherein the position of the liquid level at the current time is set as a pressure boundary condition at the next time.   The simulation apparatus according to claim 1, wherein the at least one coordinate scale is larger than 10 square times the other coordinate scale. A control method of a simulation device for an incompressible fluid or a pseudo-compressible Newtonian fluid moving in a three-dimensional space,
At least one coordinate scale of the calculation grid set in the three-dimensional space is at least twice that of the other coordinate scale, and the partial differential equation derived from the basic equation is decomposed corresponding to each coordinate scale, and further orthogonal Using the decomposed equations,
A long wavelength component calculating step of calculating a pressure and a velocity of a three-way component using an equation corresponding to the at least one coordinate scale;
A short wavelength component calculating step for calculating a pressure and a velocity of a direction component corresponding to the other coordinate scale using an equation corresponding to the other coordinate scale;
A simulation apparatus control method comprising: the long wavelength component calculation step; and a synthesis step of adding the calculation results in the short wavelength component calculation step.   A control program for causing a computer to function as the simulation apparatus according to claim 1, wherein the control program causes the computer to function as the long wavelength component calculating unit, the short wavelength component calculating unit, and the combining unit.

Images