JP5945183B2 - Simulation device, simulation method, and simulation program - Google Patents

Description

  The present invention relates to a simulation apparatus, and more particularly to a simulation apparatus, a simulation method, and a simulation program for simulating a free interface in a multiphase flow according to a volume tracking method.

  Conventionally, various proposals have been made on methods for analyzing free interfaces in multiphase flows. One of them, the volume tracking method, is a method of expressing an interface on a lattice using a volume occupancy.

  The volume occupation ratio is expressed by the following formula (1). In a two-phase flow system, the value of the function is 1 if a grid is completely occupied by one fluid, and the value of the function is 0 if it is completely occupied by the other fluid. . The value of this function takes a value from 0 to 1 depending on the proportion of one fluid.

  In the formula (1), t represents time, and α represents the occupation ratio of one fluid in each cell. Further, the v vector (represented by adding an arrow above v in the expression (1)) is the reconstructed flow rate u of the one fluid as shown in FIG. 28 of each cell. This is the vector quantity for the direction.

The occupation ratio of each cell will be described with reference to FIG. 29 and FIG.
FIG. 29 is a diagram schematically illustrating a state in which two fluids exist. In FIG. 29, the hatched portion indicates the presence area of one fluid, and the hatched portion indicates the presence area of the other fluid. FIG. 29 shows a state where 36 cells of 6 × 6 are defined. FIG. 30 shows the occupation ratio of the one fluid in each of the 36 cells in FIG.

  Various methods for reconstructing the interface as shown in FIG. 29 from the occupation ratio of a certain fluid in each cell as shown in FIG. 30 have been proposed.

  For example, Non-Patent Document 1 discloses a tangent of hyperbola for interface capturing (THINC) method. This method 1) accurately preserves the transferred flow, 2) efficiently eliminates numerical errors, and 3) avoids spurious oscillations in the vicinity of step changes in volume occupancy. For this purpose, a characteristic complementary function is adopted. In this method, a hyperbolic tangent function is used as a complementary function.

  Further, Non-Patent Document 2 discloses a weighted line interface calculation (THINC / WLIC) method in which unit normal vectors at both the x-axis and y-axis interfaces are weighted and used for reconstruction.

F. Xiao, Y. Honma, T. Kono, "A simple algebraic interface capturing scheme using hyperbolic tangent function", International Journal for Numerical Methods in Fluid 48 (2005) 1023-1040 K. Yokoi, "Efficient implementation of THINC scheme: A simple and practical smoothed VOF algorithm", Journal of Computational Physics 226 (2007) 1985-2002

  However, it cannot be said that the simulation of the flow of a complicated free interface can still accurately reproduce the actual flow, and there is still room for examination.

  The present invention has been conceived in view of such circumstances, and an object thereof is to provide a simulation apparatus, a simulation method, and a simulation program that can reconstruct a free interface more accurately.

  A simulation apparatus according to the present invention is a simulation apparatus for simulating a free interface in a multiphase flow composed of two or more kinds of fluids, wherein two or more kinds of fluids in each of a plurality of cells into which a simulation target is divided. The interface along the first axis of each cell according to the THINC method using the calculation means for calculating the occupancy ratio α of one of the fluids and the equations (A1) and (A2) including the variables β and γ reconstructing means for reconstructing χ (x),

The reconstruction means calculates the variable β in the direction along the first axis of each cell as the gradient of the normal vector of each cell with respect to the first axis.

  The simulation method according to the present invention is a simulation method executed by a computer for simulating a free interface in a multiphase flow composed of two or more kinds of fluids, and is a 2 in each of a plurality of cells into which a simulation target is divided. Along the first axis of each cell according to the THINC method using the steps (A1) and (A2) including the variables β and γ, and calculating the occupation ratio α of one of the fluids of the types or more Restructuring the interface,

Storing the constant set for each of two or more axes including the first axis defined in the simulation, which is a constant for obtaining the variable β in the formula (A1), and performing the reconstruction Calculates the variable β in the direction along the first axis of each cell as the gradient of the normal vector of each cell with respect to the first axis.

  The simulation program according to the present invention is a simulation program executed by a computer for simulating a free interface in a multiphase flow composed of two or more kinds of fluids. The simulation program is obtained by dividing a simulation target into the computer. Each cell according to the step of calculating the occupation ratio α of one of the two or more fluids in each of the plurality of cells, and the THINC method using the equations (A1) and (A2) including the variables β and γ Reconstructing the interface along the first axis of

A constant for obtaining the variable β in the formula (A1), the step of storing constants set for each of two or more axes including the first axis defined in the simulation, is executed, and reconstruction is performed. The step calculates the variable β in the direction along the first axis of each cell as the gradient of the normal vector of each cell with respect to the first axis.

  According to the present invention, it is possible to improve the accuracy of reconstruction of a free interface in simulation.

1 is a schematic configuration diagram of a computer, which is a representative hardware configuration for realizing a simulation apparatus according to an embodiment of the present invention. It is a figure which shows the change of the one-dimensional occupation rate by a THINC method. It is a figure for demonstrating the calculation method of the variable in embodiment. It is a figure which shows the outline of a THINC / WLIC method. It is a figure for demonstrating the outline of a THAINC method. It is a figure for demonstrating the outline of the calculation method of this Embodiment. It is a flowchart of the process in the simulation of this Embodiment. It is a figure for demonstrating the processing content of FIG. It is a figure for demonstrating the processing content of FIG. It is a figure for demonstrating the processing content of FIG. It is a figure for demonstrating the processing content of FIG. It is a figure for demonstrating the processing content of FIG. It is a figure which shows the density of each cell corresponding to FIG. It is a figure which shows the density of each cell corresponding to FIG. (X, y, z) coordinate system centered coordinates in (i, j, k) and cell Ci, j, of k, it is a diagram for explaining the calculation of the x-axis direction of the normal vector n _x. (X, y, z) coordinate system centered coordinates in (i, j, k) and cell Ci, j, of k, it is a diagram for explaining the calculation of the x-axis direction of the normal vector n _x. (X, y, z) coordinate system centered coordinates in (i, j, k) and cell Ci, j, of k, it is a diagram for explaining the calculation of the x-axis direction of the normal vector n _x. FIG. 18 schematically shows a cross section taken along line L1 in FIG. FIG. 18 schematically shows a cross section taken along line L2 in FIG. It is a figure which shows typically an example of the lattice of the cell defined in a three-dimensional simulation. It is a figure which shows typically an example of the calculation result of the occupation rate of each cell in a three-dimensional simulation. It is a figure which shows typically an example of the calculation result of the density of each cell in a three-dimensional simulation. It is a figure which shows the result of the reconstruction of the interface at the time of performing the simulation which carries out the rotation movement of the bubble in a liquid phase 5 times. In a two-dimensional simulation, it is a figure which shows the result of the reconstruction of the interface at the time of moving the bubble in a liquid phase. It is a figure which shows the result of the reconstruction of the interface at the time of moving the bubble in a liquid phase in three-dimensional simulation. It is a figure for demonstrating the content of (beta) in the calculation method of this Embodiment. It is a figure for demonstrating the content of (beta) in the calculation method of this Embodiment. It is a figure for demonstrating the prior art. It is a figure for demonstrating the prior art. It is a figure for demonstrating the prior art.

  Embodiments of the present invention will be described in detail with reference to the drawings. Note that the same or corresponding parts in the drawings are denoted by the same reference numerals and description thereof will not be repeated.

<Hardware configuration>
FIG. 1 is a schematic configuration diagram of a computer 100 which is a representative hardware configuration for realizing a simulation apparatus according to an embodiment of the present invention.

  Referring to FIG. 1, a computer 100 includes a computer main body 101 on which an FD (Flexible Disk) driving device 111 and a CD-ROM (Compact Disk-Read Only Memory) driving device 113 are mounted, a monitor 102, a keyboard 103, Mouse 104.

  In addition to the FD driving device 111 and the CD-ROM driving device 113, the computer main body 101 is connected to each other via a bus, which is a CPU (Central Processing Unit) 105, a memory 106, and a fixed storage device. A disk 107 and a communication interface 109 are included.

  The simulation according to the present embodiment is realized by the CPU 105 executing a program using computer hardware such as the memory 106. In general, such a program is stored in a recording medium such as the FD 112 or the CD-ROM 114, or distributed via a network or the like. Such a program is read from the recording medium by the FD driving device 111 or the CD-ROM driving device 113 or received by the communication interface 109 and stored in the fixed disk 107. Further, such a program is read from the fixed disk 107 to the memory 106 and executed by the CPU 105.

  CPU 105 is an arithmetic processing unit that performs various numerical logic operations, and provides a control function according to the present embodiment by sequentially executing programmed instructions. The memory 106 stores various information according to the program execution of the CPU 105.

  The monitor 102 is a display unit for displaying information output by the CPU 105, and includes, for example, an LCD (Liquid Crystal Display) or a CRT (Cathode Ray Tube). That is, the simulation result or the like is displayed on the monitor 102.

  The mouse 104 receives a command from a user corresponding to an operation such as click or slide. The keyboard 103 receives a command from the user corresponding to the input key.

  The communication interface 109 is a device for establishing communication between the computer 100 and another device, and can accept various data from the outside.

  Instead of the computer as described above, part or all of the computer may be realized by dedicated hardware.

  Further, a recording medium on which a program that causes the computer 100 to function in accordance with the present embodiment is not limited to the above-described FD or CD-ROM. Other examples of the recording medium include DVD-ROM (Digital Versatile Disk-Read Only Memory), USB (Universal Serial Bus) memory, memory card, hard disk, magnetic tape, cassette tape, MO (Magnetic Optical Disc), MD ( Examples thereof include a medium for storing a program in a nonvolatile manner, such as a Mini Disc, an IC (Integrated Circuit) card (excluding a memory card), an optical card, a mask ROM, an EPROM, and an EEPROM (Electronically Erasable Programmable Read-Only Memory).

<Summary of simulation method of this embodiment>
As described above, the THINC method and the THINC / WLIC method have been proposed for the reconstruction of the free interface. On the other hand, in the simulation method of the present embodiment, the constant (β) whose value is fixed in the THINC method is set as the gradient for each transport direction from the normal vector.

  Therefore, in this specification, the simulation method of the present embodiment will be described first after describing the THINC method and the THINC / WLIC method which are the premise of the simulation method of the present embodiment. Further, in this specification, a THAINC (tangent of hyperbola with adaptive function for interface capturing) method for setting the constant (β) in a mode different from the present embodiment will be described as a reference.

  In the present specification, the simulation method of the present embodiment is appropriately referred to as a “THAINC_M (tangent of hyperbola with adaptive function for interface capturing_modified) method”.

  Further, in this specification, the description will be given mainly assuming two axes of the x axis and the y axis as axes defined in the simulation. That is, two-dimensional simulation is mainly described. However, the simulation method disclosed in the present embodiment can be extended to three dimensions according to a known technique. In addition, about the method of extending to three dimensions, since a well-known technique is employable, description is not repeated here.

<THINC method>
FIG. 2 is a diagram showing a change in one-dimensional occupancy by the THINC method. The THINC method uses a smoothed hyperbolic tangent function as shown in FIG. In addition, the broken line of FIG. 2 has shown the change by the volume tracking method.

  In the THINC method, the following expression (2) is used as a characteristic function.

In equation (2), γ takes a value of “1” or “−1” depending on the value of _{nx, i} , as represented by the following equation (3).

N _{x, i in} Equation (3) is the x component of the normal vector, and the occupancy ratio α of the remaining nine cells when the calculation target cell is placed at the center of a 3 × 3 grid Is calculated based on Specifically, as shown in FIG. 3, 9 cells of a 3 × 3 grid are assumed. In the upper stage, the occupation ratio of each cell is α _{i−1, j + 1} , α _{i, j + 1} , α _{i + 1, j + 1} in this order from the left. In the middle row, α _i−1,1 , α _{i, j} , and α _{i + 1, j} are sequentially from the left. In the lower row, from the left, α _{i−1, j−1} , α _{i, j−1} , α _{i + 1, j−1} are in order. Based on these occupancy rates, nx _{, i} is calculated from the following equations (4-1) to (4-4) and (5). In these formulas, “NXRT”, “NXLT”, “NXRB”, and “NXLB” are the occupations when the four vertices of the center cell shown in FIG. 3 are moved by a minute distance in the x-axis direction. It shows the rate change. In Expression (4-1), “NXRT” represents the change amount of the occupation ratio in the upper right portion of the center cell. In the formula (4-2), “NXLT” represents the change amount of the occupation ratio of the upper left portion of the center cell. In the formula (4-3), “NXRB” represents the change amount of the occupation ratio in the lower right portion of the center cell. In Expression (4-4), “NXLB” represents the amount of change in the occupation ratio of the lower left portion of the center cell.

  Returning to Equation (2), β is a parameter that determines the slope of the hyperbolic tangent function shown in FIG. In general, a fixed value of 3.5 is employed as β.

“X _i ˜” in equation (2) is calculated according to the following equation (6) using the average of χ (x) in equation (2). In the present specification, the tilde symbol is represented by “˜”.

That is, χ (x) in Expression (2) is constructed using γ and β set as described above, and “x _i ˜” is calculated based on Expression (6).

<THINC / WLIC method>
In the THINC / WLIC method, the flux in the x direction is calculated according to equation (7). In addition, "RF" in Formula (7) means Right Face. That is, it means the rightmost surface (the surface on which the vertical line where the black circle shown by u in FIG. 5 described later exists), and more specifically, the cell number “i + 1/2”.

  On the right side of Equation (7), the first term is the flux according to the THINC method, and the second term is the flux according to the UPWIND method (primary upwind method).

In equation (7), ω _{x, i} is calculated according to the following equation (8).

In Expression (8), _{nx, i} is calculated from Expression (4-1) to Expression (4-4) and Expression (5) as described above. Also, _{ny, i} is calculated from the following formulas (9-1) to (9-4) and (10) based on the occupancy rate of nine cells as shown in FIG.

That is, in the THINC / WLIC method, as shown in FIG. 4, when ny _{, i} is 0, ω _{x, i} becomes 1 according to the equation (8), and the flux in the x direction is It is calculated based only on the first term on the right side of Equation (7). When n _{x, i} is 0, ω _{x, i} is 0 according to the equation (8), so that the flux in the x direction is based only on the second term on the right side of the equation (7). Calculated.

In the state between these cases, ω _{x, i} is calculated according to the ratio of the x component and the y component of the velocity vector in equation (8). Thereby, the flux in the x direction is calculated by combining the first term and the second term on the right side of Equation (7) according to Equation (7) based on the calculated values of ω _{x, i} . The first term on the right side of Equation (7) corresponds to the calculation of the flux by the THINC method, and the second term corresponds to the calculation of the flux by the UPWIND method.

<THAINC method>
The THAINC method, which is a simulation method of a comparative example of the present embodiment, will be described. The THAINC method uses constants (β) whose values have been fixed in the THINC method, constants set for each axis (βu, βl described later), and further, for each axis in the velocity vector of the fluid taking account of the occupation ratio. A method of setting based on the contribution degree (ω described later) is given.

FIG. 5 is a diagram for explaining the outline of the THAINC method.
Referring to FIG. 5, the THAINC method is characterized by a method for setting a variable β in the special function of the THINC method represented by Expression (2).

Specifically, the variable β is calculated and set using β _u , β _l, and ω _{x, i} , as shown in Expression (11) in FIG. The values of β _u and β _l are stored in the fixed disk 107, for example.

β _u is a value that changes the vicinity of the interface more rapidly than the value 3.5 conventionally used as β, and is set to 10, for example. Note that β _u in FIG. 5 shows the result when γ ≧ 0 in the equation (2). β _l is a value that makes the change in the vicinity of the interface more gradual than the value 3.5 conventionally used as β, and is, for example, 0.001. FIG. 5 shows a change in the occupation ratio in the vicinity of the interface when each of β _u and β _l is used as β in the equation (2).

According to the THAINC method, β is set according to the magnitude relationship of normal vectors (in the case of two-dimensional simulation, the magnitude relationship between nx _{, i} and ny _{, i} ).

Specifically, when _{nx, i} is larger than _{ny, i} , that is, when the x component of the velocity vector of the cell is larger than the y component, in the reconstruction of the cell, β is set so that the change in the occupation ratio of the fluid at the interface in the x direction becomes more rapid. This is because, in the setting of β according to the equation (11), β _u is reflected at a rate larger than β ₁ .

On the other hand, when _{nx, i} is smaller than _{ny, i} , that is, when the x component of the velocity vector of the cell is smaller than the y component, the cell is reconstructed in the x direction more than before. Β is set so that the change in the occupation ratio of the fluid at the interface becomes more gradual. This is because in the setting of β according to the equation (11), β ₁ is reflected at a rate larger than β _u .

That is, in the reconstruction of the interface according to the present embodiment, the contribution degree of the component in the reconstruction axial direction in the velocity vector in the field for a plurality of axes (in this embodiment, _{nx, i} And the degree of contribution of _{nx, i to} the sum of _{ny, i} ) will be reflected more greatly.

  The calculation flow in the THAINC method can be the same as the calculation flow in the simulation of the present embodiment described with reference to FIG.

  That is, first, in step S10, the initial values of the occupation ratio α and the velocity v of each cell are read. In step S20, the velocity v of each cell and the normal vectors in the x and y directions of the velocity are calculated. In step S30, “∇α” for the latest occupancy α is calculated. In step S40, “∇αv” for the product of the occupancy α and the speed v of each cell is calculated. In step S50, “の αv” is calculated. Occupancy rate α after Δt is calculated with respect to occupancy rate α, and in step S60, γ is calculated according to equation (3) based on the latest occupancy rate α, and β is calculated according to equation (11). By calculating “xi˜” according to the equations (2) and (6), the interface of each cell is reconstructed.

In step S70, when the simulation period (sum of Δt) _Tsum exceeds the specified period T, the process ends.

Since the THAINC method can be extended to a three-dimensional simulation by adopting a well-known technique, detailed description will not be repeated here. For example, in the three-dimensional simulation, ω _{x, i} is a normal vector n _{x in} the x-axis direction, a normal vector n _{y in the} y-axis direction, and a z-axis direction as shown in the following expression (12). Of the normal vector _nz .

<Calculation method of the present embodiment>
In the calculation method according to the present embodiment, each cell is basically reconstructed by a calculation method according to the THAINC method.

FIG. 6 is a diagram for explaining the outline of the calculation method of the present embodiment.
Referring to FIG. 6, the calculation method of the present embodiment is characterized by a variable β setting method in the special function of the THINC method expressed by the equation (2). Specifically, the variable β is calculated and set as shown in Expression (13) in FIG.

  Equation (13) shows β used for reconstruction in the k-axis (one of the axes defined in the simulation) direction. j means an axis other than the k-axis of the orthogonal coordinate axis defined in the simulation only in Expression (13). “K” is used as a character representing an axis only in equation (13).

  In the following formulas (14) to (16), β is calculated for each of the simulation models of the one-dimensional model (1D), the two-dimensional model (2D), and the three-dimensional model (3D). An expression is shown.

  First, when the simulation model is a one-dimensional model, since the j-axis does not exist, β is a constant (3.5) as shown in Expression (14).

When the simulation model is a two-dimensional model, as shown in equation (15), beta is a normal vector n _x in the x-axis _direction, and, using the normal vector n _y in the y-axis direction Table Is done. Specifically, when the x axis and the y axis are defined in the two-dimensional simulation, if the k axis is the x axis, the j axis is the y axis, and if the k axis is the y axis, the j axis is the x axis. It is. Then, beta is used to reconstruct the x-axis direction (i.e., beta _x) is calculated n _x as a value obtained by dividing n _y. Furthermore, beta is used to reconstruct the y-axis direction (i.e., beta _y) is calculated n _y as divided by the n _x.

When the simulation model is a three-dimensional model, as shown in Expression (16), β is a normal vector n _{x in} the x-axis direction, a normal vector n _{y in the} y-axis direction, and a z-axis direction. Of the normal vector _nz . Specifically, when the x-axis, y-axis, and z-axis are defined in the three-dimensional simulation, if the k-axis is the x-axis, the j-axis is the y-axis and the z-axis, and the k-axis is the y-axis For example, the j axis is the z axis and the x axis, and if the k axis is the z axis, the j axis is the x axis and the y axis. Then, beta is used to reconstruct the x-axis direction (i.e., beta _x) is calculated n _x as a value obtained by dividing by the square root of the sum of the squares of n _y and n _z. is used to reconstruct the y-axis direction beta (i.e., beta _y) is calculated n _y as the square root divided by the sum of the squares of n _z and n _x. Then, beta is used to reconstruct the z-axis direction (i.e., beta _z) is calculated n _z as the square root divided by the sum of the squares of n _x and n _y.

<Calculation flow in the present embodiment>
FIG. 7 is a flowchart of processing in the simulation of the present embodiment.

  In the following description, the movement of a gas in a liquid will be exemplified as a simulation target. Thereby, the occupation ratio of the gas in each cell is illustrated as occupation ratio (alpha). Note that the combination of phases in the multiphase flow to be simulated in the present embodiment is not limited to this. As long as two or more different fluids are present, a system having two or more different fluids may be used.

  Referring to FIG. 7, when an instruction to start calculation is input to computer 100, CPU 105 first reads initial values of occupancy α and speed v of each cell in step S10, and proceeds to step S20. To proceed. These values are stored in the fixed disk 107, for example.

  In step S20, the CPU 105 calculates the velocity v of each cell and the normal vectors in the x and y directions of the velocity, and proceeds to step S30.

  In step S30, the CPU 105 calculates “∇α” for the latest occupancy α and advances the process to step S40.

  In step S40, “∇αv” for the product of the occupation ratio α and the speed v of each cell is calculated, and the process proceeds to step S50.

  In step S50, the occupation rate α after Δt is calculated with respect to the occupation rate α of each cell, and the process proceeds to step S60.

In step S60, γ is calculated according to equation (3) based on the latest occupation ratio α, β is calculated according to equation (13), and “x _i ˜” is calculated according to equation (2) and equation (6). As a result, the interface of each cell is reconstructed, and the process proceeds to step S70.

  In step S70, the CPU 105 determines whether or not the simulation period (sum of Δt) Tsum exceeds the specified period T. If not, the process returns to step S20, and if it exceeds, the process is performed. Terminate.

  As described above, the occupancy rate of each cell is updated every minute time Δt by the process described with reference to FIG. 7, and the interface of each cell is described with reference to FIG. It is reconstructed according to the calculated calculation method of the present embodiment.

  In the simulation of the present embodiment, a cell lattice as shown in FIG. 8 (100 × 100 cell lattice in FIG. 8) is defined as a simulation target region. Each time the processing from step S20 to step S60 is performed once, the calculation result of the occupation ratio α of each cell is updated, and the interface in each cell is reconstructed.

  9 and 10 show the calculation results of the occupancy α of each cell for a part of the lattice-like cell shown in FIG. 8 (10 × 10 cell lattice). In addition, FIG. 10 shall show the calculation result after time (DELTA) T with respect to the calculation result of FIG.

  FIGS. 11 and 12 are diagrams showing the state of interface reconstruction corresponding to the calculation results of FIGS. 9 and 10, respectively. That is, FIG. 12 shows the result after time ΔT with respect to the result shown in FIG.

  In FIG. 9, in the cells other than the 16 cells in the 3rd to 6th columns and the 3rd to 6th rows, the occupation ratio α is set to “0.0”. On the other hand, in these 16 cells, a value of 0.4, 0.8, or 1.0 is calculated as the occupation ratio α. And the interface is reconstructed as shown in FIG. 11 corresponding to each occupation ratio α.

  In FIG. 10, with respect to FIG. 9, 16 cells each having a value other than 0.0 in the occupation ratio α are moved to the right by one column. In FIG. 12, correspondingly, the reconstructed interfaces are moved to the right with respect to those shown in FIG.

Further, the density of each cell is calculated according to the following equation (14) based on the occupation ratio α of each cell as shown in FIG. 11 and FIG. FIG. 13 shows the density calculation result corresponding to FIG. 11, and FIG. 14 shows the density calculation result corresponding to FIG. In equation (14), ρ _G is the gas density, and ρ _l is the liquid density. Further, in FIGS. 13 and 14, is 1 kg / m ³ as an example of the value of [rho _G, is 1000 kg / m ³ as an example of the value of [rho _l, are employed respectively.

<Extension to 3D>
In the simulation method (THAINC) of the present embodiment, as described with reference to the equation (13) in FIG. 6, when the simulation is performed in two or more dimensions, each axis of the cell to be reconstructed The normal vector of was used. The calculation of the normal vector of each cell in the three-dimensional simulation will be described.

15, (x, y, z) the center coordinates in a coordinate system (i, j, k) and cell C _{i which, j,} of _k, for explaining the calculation of the normal vector n _x in the x-axis direction FIG.

In FIG. 15, C _{i, j, k} is shown as a cube, and the amount of change in the occupancy when moving a minute distance in the x-axis direction in the vicinity of the eight vertices of the cube is “NXRTF”, “ NXLTF ”,“ NXRBF ”,“ NXLBF ”,“ NXRTB ”,“ NXLTB ”,“ NXRBB ”, and“ NXLBB ”.

  These change amounts are calculated according to the following equations (16-1) to (16-8).

  As understood from the equations (16-1) to (16-8), the occupancy rate α of the cell to be calculated and the cell adjacent to the cell is used for the calculation of the above eight change amounts. The

FIG. 16 shows cells used for the calculation of the above eight variations.
In FIG. 16, 27 cells are shown that are arranged three by three in the three-axis direction around the calculation target cell C _{i, j, k} . For the 27 cells, in FIG. 16, the x coordinates of the three cells arranged in the x-axis direction are indicated by “i−1”, “i”, and “i + 1”. In addition, y coordinates of three cells arranged in the y-axis direction are indicated by “j−1”, “j”, and “j + 1”. In addition, z coordinates of three cells arranged in the z-axis direction are indicated by “k−1”, “k”, and “k + 1”. That is, the three-dimensional coordinates of the 27 cells are (i-1, j-1, k-1) (i-1, j, k-1) (i-1, j + 1, k-1) ( i-1, j-1, k) (i-1, j, k) (i-1, j + 1, k) (i-1, j-1, k + 1) (i-1, j, k + 1) (i-1, j + 1, k + 1) (i, j-1, k-1) (i, j, k-1) (i, j + 1, k-1) (i , j-1, k) (i, j, k) (i, j + 1, k) (i, j-1, k + 1) (i, j, k + 1) (i, j + 1, k + 1) (i + 1, j-1, k-1) (i + 1, j, k-1) (i + 1, j + 1, k-1) (i + 1, j-1, k) (i + 1, j, k) (i + 1, j + 1, k) (i + 1, j-1, k + 1) (i + 1, j, k + 1) (i + 1 , j + 1, k + 1).

Then, the normal vector _nx is calculated according to the following equation (19) using these eight variations.

That is, the normal vector _nx is calculated as an average value of the above eight variations.
Returning to the calculation of the eight change amounts, each change amount is calculated using the occupancy ratio of the eight cells located in the vicinity of the vertex. FIG. 16 shows “NXLTF” which is one of the eight vertices of the cell C _{i, j, k} .

For example, as shown in FIG. 17, the amount of change in the occupancy ratio in the minute region in the x-axis direction near the vertex indicated by “NXLTF” is calculated in the cell C _{i, j, k} to be calculated, 7 cells (C _{i, j + 1, k} , C _{i, j + 1, k−1} , C _{i, j, k−1} , C _{i−1, j + 1, k} , C _{i−1) , j, k-1} ) is used. The subscript for C represents the coordinates of each cell. The occupancy of each of these cells is α _{i, j, k} , α _{i, j + 1, k} , α _{i, j + 1, k−1} , α _{i, j, k−1} , α _{i−1, j + 1, k} , α _{i−1, j + 1, k−1} , α _{i−1, j, k} , α _{i−1, j, k−1} .

  FIG. 18 schematically shows a cross section taken along line L1 in FIG. The cross section corresponds to the cross section of the eight cells in FIG. 17 at z = k. FIG. 19 schematically shows a cross section taken along line L2 in FIG. 18 and 19 show the occupancy rates of the respective cells. The cross section corresponds to a cross section at z = k−1 of the eight cells in FIG.

Then, “NXLTF” is an occupancy ratio of two cells adjacent in the direction of the axis (x-axis) for obtaining the normal vector for the eight cells of FIG. 17 as shown in the equation (16-2). It is calculated as an average of values obtained by dividing the difference by a minute distance (Δx) in the x-axis direction. That is, in these eight cells, there are four combinations of cells adjacent in the x-axis direction (C _{i, j + 1, k−1} and C _{i−1, j + 1, k−1} , C _{i, j, k-1} and C _{i-1, j, k-1} , C _{i, j + 1, k} and C _{i-1, j + 1, k} , C _{i, j, k} and C _{i-1, j , k} ). Therefore, in the numerator of the equation (16-2), the difference in occupancy in each of these four groups is expressed, so that in the denominator of the equation (16-2) “4” is described together with the minute distance (Δx).

Similarly, in the other formulas (16-1) to (16-8), the occupancy ratio of two cells adjacent in the x-axis direction in the vicinity of the eight vertices of the cells C _{i, j, k} An average of values obtained by dividing the difference by a minute distance (Δx) in the x-axis direction is calculated. The results calculated by these equations is utilized in accordance with the equation (19), the normal vector n _x in the x-direction is calculated.

Further, the normal vector n _y and z direction of the normal vector n _z in the y-direction, respectively, is calculated the following equation (20) in accordance with equation (21).

“NYRTF”, “NYLTF”, “NYRBF”, “NYLBF”, “NYRTB”, “NYLTB”, “NYRBB”, “NYLBB” in Formula (20), and “NZRTF”, “NZLTF” in Formula (21) The eight numerical values of “NZRBF”, “NZLBF”, “NZRTB”, “NZLTB”, “NZRBB”, and “NZLBB” are respectively the eight cells C _{i, j, k to} be reconstructed as shown in FIG. This is the rate of change of the occupancy rate in the minute area near the apex. Incidentally, 8 numbers in the formula (20) are the numerical values for determining the normal vector n _y in the y-direction, y-axis the difference in occupancy of four sets of two cells adjacent to the y-axis direction It is calculated as the average of the values divided by the minute distance (Δy) in the direction. In addition, since the eight numerical values in the equation (21) are numerical values for obtaining the normal vector _nz in the z direction, the difference between the occupancy ratios of four sets of two cells adjacent in the z axis direction is expressed as the z axis. It is calculated as the average of the values divided by the minute distance (Δz) in the direction.

As described above, in the three-dimensional simulation, and the normal vector n _z of the normal vector n _y and z direction of the normal vector n _x and y direction of the x-direction is calculated. These are used in the equation (16) to calculate β for each of the x direction, the y direction, and the z direction, thereby realizing reconstruction of each cell based on the three-dimensional simulation.

  In the three-dimensional simulation, a lattice of cells arranged in a three-dimensional direction (100 × 100 × 100 cell lattice in FIG. 20) is defined as a simulation target region. Then, each time the process from step S20 to step S60 in FIG. 7 is performed once, the calculation result of the occupancy rate α of each cell is updated as shown in FIG. Used to reconstruct the interface in each cell.

Further, the density of each cell is calculated according to the above equation (17) based on the occupation ratio α of each cell as shown in FIG. FIG. 22 is a diagram illustrating a density calculation result corresponding to FIG. In Figure 22, the 1 kg / m ³ as an example of the value of [rho _G, is 1000 kg / m ³ as an example of the value of [rho _l, are employed respectively.

<Comparison of calculation results>
Next, the result of the reconstruction of the interface according to the present embodiment is compared with the result of the reconstruction by various methods so far.

  FIG. 23 to FIG. 25 are diagrams showing simulation results by the THINC / WLIC method, the THAINC method, and the THAINC_M method for different simulation contents. In each figure, a graph in which the result of the simulation is plotted, a table indicating the numerical value of the result, and an explanatory diagram of the contents of the simulation are shown in order from the top.

In the graph, for each bubble, the error (ERR) relative to the initial shape for each bubble is plotted against the Courant number (CO). The error is shown as a percentage. The results of the THINC / WLIC method are plotted as diamonds, the results of the THAINC method as squares, and the results of the THAINC_M method as triangles. This error is calculated as an error (E _SP ) according to the following equation (22).

  In the table, the Courant number (CO) and the percentage value of the error corresponding to CO are shown.

(Rotate 5 times)
FIG. 23 is a diagram illustrating a result of the reconstruction of the interface when a simulation is performed in which bubbles in the liquid phase are rotated and moved five times. In this simulation, as shown in the explanatory diagram, the bubble in the initial state indicated by the bubble BU is rotated five times in the liquid LQ in the rotation direction indicated by the arrow R1 and returned to the original position.

  As understood from the graph and table of FIG. 23, the THAINC_M method has an error of about 1% lower than that of the THINC / WLIC method over the almost entire range of the indicated Couran number. Even the error is low.

(Movement in 2D simulation)
FIG. 24 is a diagram showing a result of reconstruction of the interface when bubbles in the liquid phase are moved in the two-dimensional simulation. In this simulation, as shown in the explanatory diagram, the bubble in the initial state indicated by the bubble BU is moved in the liquid LQ according to the arrow R2.

  As understood from the graph and table of FIG. 24, the THAINC_M method has a lower error than the THINC / WLIC method and the THAINC method over almost the entire range of the indicated Couran number. In particular, the THAINC_M method is about 1.5% when CO is 0.5, close to 3% when CO is 0.25, and when CO is 0.1. The error is low by 3% or more.

(Movement in 3D simulation)
FIG. 25 is a diagram showing a result of reconstruction of the interface when bubbles in the liquid phase are moved in the three-dimensional simulation. In this simulation, as shown in the explanatory diagram, the bubble in the initial state indicated by the bubble BU is moved in the three-dimensional space in the liquid LQ according to the arrow R3. In FIG. 25, the bubble after the movement is indicated by a bubble BU01.

  As can be understood from the graph and table of FIG. 25, the magnitude relationship of the error is switched between the THINC / WLIC method and the THAINC method depending on the value of CO. However, the THAINC_M method Over almost the entire area, the error is about 3% lower than that of the THINC / WLIC method and THAINC method.

<Summary of Embodiment>
The THAINC_M method of the present embodiment described above is characterized by the method of calculating β of the special function of the THINC method shown in Equation (2) (see FIG. 6) when two or more axes are defined. . Specifically, β is calculated using the normal vector for each axis of each cell. More specifically, as shown in Equation (15), Equation (16), and the like, β used for reconstruction in the first axial direction is the normal vector of the first axis, Calculated as the ratio of the axis normal vector to the combined vector.

In FIG. 26, the x component (vector n _x ) and y component (vector n _y ) of a normal vector for a certain cell in the two-dimensional simulation are schematically shown along with the interface SB and the normal vector (vector n). It is shown. In the example shown in FIG. 26, the value of the fluid occupancy α that is focused on the left and right of the interface SB changes. Α = 0 on the left side of the interface SB, and α = 1 on the right side. In this case, as described with reference to the equation (15), β (that is, β _x ) used for reconstruction in the x-axis direction is a gradient of the normal vector n with respect to the x-axis direction. Calculated. Further, β (that is, β _y ) used for reconstruction in the y-axis direction is calculated as a gradient of the normal vector n with respect to the y-axis direction.

FIG. 27 shows an x component (vector n _x ) of a normal vector for a certain cell in a three-dimensional simulation, a normal vector (vector n), and a combined vector (vector n _yz ) of its y component and z component. Are schematically shown. In FIG. 27, the y component of the normal vector n is indicated by “n _y ” on the y axis, and the z component of the normal vector n is indicated by “n _z ” on the z axis. Further, the scalar quantity of the combined vector (vector n _yz ) is expressed by the following equation (23).

β (that is, β _x ) used for reconstruction in the x-axis direction is calculated as a gradient of the normal vector n with respect to the x-axis direction, that is, as shown in Expression (16).

  According to the THAINC_M method of the present embodiment described above, in a simulation of a free interface in a multiphase flow composed of two or more fluids, even if the simulation dimension changes, the error is suppressed to be lower than that of the conventional calculation method. be able to.

<Other variations, etc.>
In the present embodiment described above, the CPU 105 that executes step S50 constitutes a calculation means for calculating the occupancy of each cell. The CPU 105 that executes step S60 constitutes a reconstructing unit that reconstructs the interface of each cell.

  100 computer, 101 computer main body, 102 monitor, 103 keyboard, 104 mouse, 106 memory, 107 fixed disk, 109 communication interface, 111 FD drive device, 113 CD-ROM drive device.

Claims (6)

A simulation device for simulating a free interface in a multiphase flow composed of two or more kinds of fluids,
A calculating means for calculating an occupancy ratio α of one of the two or more types of fluid in each of the plurality of cells into which the simulation target is divided;
Reconstructing the interface χ (x) along the first axis of each cell according to the tangent of hyperbola for interface capturing (THINC) method using equations (A1) and (A2) including variables β and γ Restructuring means for

Said reconstruction means, said variable β in the direction along the first axis of each said cell is calculated using the formula (A3), the simulation device.

(In the formula (A3), the letter k represents the first axis, the letter _nk represents the normal vector in the first axial direction, and the letter j represents the first coordinate axis defined in the simulation. An axis other than the axis is represented, and the letter n _j represents a normal vector of an axis other than the first axis.) A simulation device for simulating a free interface in a multiphase flow composed of two or more kinds of fluids,
A calculating means for calculating an occupancy ratio α of one of the two or more types of fluid in each of the plurality of cells into which the simulation target is divided;
Reconstructing the interface χ (x) along the first axis of each cell according to the tangent of hyperbola for interface capturing (THINC) method using equations (A1) and (A2) including variables β and γ Restructuring means for

The reconstructing means calculates the variable β with respect to the direction along the first axis of each cell by using the formula (A4) or the formula (A5).

(In Formula (A4) and Formula (A5), the letter n _x Represents the normal vector in the x-axis direction and the letter n _y Represents the normal vector in the y-axis direction and the letter n _z Represents a normal vector in the z-axis direction. ) A computer-implemented simulation method for simulating a free interface in a multiphase flow comprising two or more fluids,
Calculating an occupancy ratio α of one of the two or more fluids in each of a plurality of cells into which a simulation target is divided;
Reconstructing the interface along the first axis of each said cell according to a tangent of hyperbola for interface capturing (THINC) method using equations (A1) and (A2) comprising variables β and γ;

A constant for obtaining the variable β in the formula (A1), the constant being set for each of two or more axes including the first axis defined in the simulation,
The step of performing the rebuilding is a simulation method in which the variable β for the direction along the first axis of each of the cells is calculated using Expression (A3) .

(In the formula (A3), the letter k represents the first axis, the letter _nk represents the normal vector in the first axial direction, and the letter j represents the first coordinate axis defined in the simulation. An axis other than the axis is represented, and the letter n _j represents a normal vector of an axis other than the first axis.) A computer-implemented simulation method for simulating a free interface in a multiphase flow comprising two or more fluids,
Calculating an occupancy ratio α of one of the two or more fluids in each of a plurality of cells into which a simulation target is divided;
Reconstructing the interface χ (x) along the first axis of each cell according to the tangent of hyperbola for interface capturing (THINC) method using equations (A1) and (A2) including variables β and γ Steps,

A constant for obtaining the variable β in the formula (A1), the constant being set for each of two or more axes including the first axis defined in the simulation,
The step of performing the reconstruction is a simulation method in which the variable β with respect to the direction along the first axis of each cell is calculated using the formula (A4) or the formula (A5).

(In Formula (A4) and Formula (A5), the letter n _x Represents the normal vector in the x-axis direction and the letter n _y Represents the normal vector in the y-axis direction and the letter n _z Represents a normal vector in the z-axis direction. ) A simulation program executed by a computer for simulating a free interface in a multiphase flow composed of two or more fluids,
The simulation program is stored in the computer.
Calculating an occupancy ratio α of one of the two or more fluids in each of a plurality of cells into which a simulation target is divided;
Reconstructing the interface along the first axis of each said cell according to a tangent of hyperbola for interface capturing (THINC) method using equations (A1) and (A2) comprising variables β and γ;

A constant for obtaining the variable β in the formula (A1), the constant being set for each of two or more axes including the first axis defined in the simulation,
The step of performing the reconstruction is a simulation program for calculating the variable β with respect to the direction along the first axis of each of the cells using the formula (A3) .

(In the formula (A3), the letter k represents the first axis, the letter _nk represents the normal vector in the first axial direction, and the letter j represents the first coordinate axis defined in the simulation. An axis other than the axis is represented, and the letter n _j represents a normal vector of an axis other than the first axis.) A simulation program executed by a computer for simulating a free interface in a multiphase flow composed of two or more fluids,
Before Symbol simulation program, in the computer,
Calculating an occupancy ratio α of one of the two or more fluids in each of a plurality of cells into which a simulation target is divided;
Reconstructing the interface χ (x) along the first axis of each cell according to the tangent of hyperbola for interface capturing (THINC) method using equations (A1) and (A2) including variables β and γ Steps,

A constant for obtaining the variable β in the formula (A1), the constant being set for each of two or more axes including the first axis defined in the simulation,
The step of performing the reconstruction is a simulation program for calculating the variable β with respect to a direction along the first axis of each of the cells using the formula (A4) or the formula (A5).

In (Equation (A4) and formula (A5), the character _{n x} represents the normal vector of the x-axis direction, the character _{n y} represents a normal vector of the y-axis direction, the character _{n z} in the z-axis direction normal Represents a vector.)

Images