JP5516949B2 - Numerical fluid calculation method, numerical fluid calculation device, program, and recording medium - Google Patents

Description

  The present invention relates to a numerical fluid calculation method and a numerical fluid calculation apparatus for solving a convection diffusion equation using a finite volume method. The present invention also relates to a program for operating a computer as such a numerical fluid computing device, and a recording medium on which such a program is recorded.

  With the advancement of computers, manufacturing technology has made tremendous progress. For example, a CAD (Computer Aided Design) system that supports design, a CAE (Computer Aided Engineering) system that supports analysis, a CAM (Computer Aided Manufacturing) system that supports processing, and a CAT (Computer Aided Testing) system that supports inspection , It has become an indispensable tool for making things.

  However, these systems are realized by separate software, and the data structure is not versatile. For example, a wire frame model widely used in a CAD system represents the surface shape of an object and does not represent the internal structure of the object. Therefore, the wire frame model cannot be used for the simulation as it is.

  A voxel model is known as a general-purpose modeling technique that can be used for both design and simulation. The voxel model is a modeling technique for expressing an object as a set of cubes called voxels. By making each voxel bear a physical quantity, the voxel model can be used for both design and simulation. However, when the object is expressed as a set of voxels, information on the surface of the object (for example, the direction of the normal line) is lost. Therefore, the boundary value problem cannot be solved with high accuracy by the numerical fluid calculation using the voxel model.

  As an improvement of the voxel model, a cut cube model is known in which a voxel that intersects the surface of an object is replaced with a cut cube (referred to as “Kitta Cube”). The cut cube is a (6 + N) plane obtained by cutting a cube having at most one cutting point on each side with N triangular cutting planes having three cutting points as vertices. The cutting point is provided at the intersection of each side of the base cube and the surface of the object. The cut cube model is suitable for the boundary value problem because it has a surface made of cut triangles that closely approximates the surface of the object.

However, there are 10 ⁶⁴ patterns of triangular cut surfaces in the cut cube, and 144 disposition patterns of cut points. Therefore, the numerical fluid calculation algorithm using the cut cube model includes many cases and it is not practical to code this. Even if coding can be performed, the amount of calculation becomes enormous, and the required processing speed cannot be achieved.

  An attempt to solve such a problem is described in Non-Patent Document 1. Specifically, it is expected that three-dimensional numerical fluid calculation will be possible by replacing the cut cube in the cut cube model with a cut cell (for example, replacing two triangular cut surfaces sharing one side with one cut surface). It is shown.

Thunder, 1 others, “Simulation of incompressible viscous fluids by VCAD”, Abstracts of 19th Symposium on Fluid Dynamics, December 2005, p. 51, CD-ROM NO. A1-3

  However, when the cut cube in the cut cube model is replaced with a cut cell as described in Non-Patent Document 1, the cut surfaces of all the cut cells cannot be aligned with the cut surfaces of the adjacent cut cells without contradiction. . That is, there is a hole in the surface of the cut cell model constituted by the cut surface of each cut cell, which causes a problem that the accuracy of the solution is greatly reduced.

  The present invention has been made in view of the above problems, and in the numerical fluid calculation method for solving the convection diffusion equation using the finite volume method, the number of necessary cases can be divided without replacing the cut cube with a cut cell or the like. It is to realize a numerical fluid calculation method in which the above is reduced.

  In order to solve the above problems, the numerical fluid method according to the present invention is an algebraic equation obtained by integrating a convection diffusion equation on each cell, and includes a physical quantity on each cell and an adjacent one adjacent to the cell. In a numerical fluid calculation method in which an equation system consisting of algebraic equations established between the above physical quantities on a cell is solved by a computer, a cell division step for dividing the calculation target region into a plurality of cells, the boundary of the calculation target region being defined As a boundary cell to be included, a hexahedron having at most one cutting point on each side is cut by N triangular cutting planes having three cutting points as vertices, and (6 + N) planes (N is equal to or less than Nmax) Cell division step using a natural number), a cutting area calculating step for calculating the area of the N triangular cutting planes for each of the boundary cells, and for each of the boundary cells, The integral value on each triangular cut surface of the convection flux and the integral value on each triangular cut surface of the diffusive flux included in the notation algebraic equation are calculated for each triangle calculated in the cutting area calculation step. It includes a boundary condition setting step that is calculated from the area of the cut surface and the boundary condition imposed on the boundary of the calculation target region.

  In order to solve the above problem, the numerical fluid calculation device according to the present invention is an algebraic equation obtained by integrating the convection diffusion equation on each cell, and includes a physical quantity on each cell, In a numerical fluid computing device that solves an equation system consisting of algebraic equations established between adjacent physical cells on adjacent cells, cell dividing means for dividing a calculation target region into a plurality of cells, the boundary of the calculation target region being defined As a boundary cell to be included, a hexahedron having at most one cutting point on each side is cut by N triangular cutting planes having three cutting points as vertices, and (6 + N) planes (N is equal to or less than Nmax) A natural number), a cutting area calculation means for calculating the area of the N triangular cutting planes for each of the boundary cells, and each of the boundary cells included in the algebraic equation. The integral value on each triangular cut surface of the convection flux, and the integral value on each triangular cut surface of the diffusion flux, the area of each triangular cut surface calculated by the cut area calculating means, and the above Boundary condition setting means for calculating from the boundary condition imposed on the boundary of the calculation target region.

  According to the above configuration, as a boundary cell including the boundary of the calculation target region, a hexahedron having at most one cutting point on each side is cut by N triangular cutting planes having three cutting points as vertices. A (6 + N) face body obtained as described above is used. Therefore, the boundary of the calculation target region can be approximated with high accuracy by these triangular cut surfaces. Since the convection term and diffusion term included in the algebraic equation to be solved are calculated using these triangular cut surfaces, the algebraic equation to be solved can be calculated with high accuracy. Therefore, the value of each physical quantity on each cell can be calculated with high accuracy.

  Moreover, the number of cut surfaces (number of triangular cut surfaces) N in each boundary cell is limited to Nmax or less. Therefore, it is possible to reduce the case division necessary for calculating the (algebraic coefficient) to be solved as compared with the case where the number of cut surfaces in each boundary cell is not limited.

  In the numerical fluid calculation method according to the present invention, for each of the boundary cells, a boundary area calculation step of calculating an area of a boundary surface between adjacent cells adjacent to the boundary cell; A coefficient calculation step that is a coefficient included in the algebraic equation and calculates a coefficient related to the physical quantity on an adjacent cell adjacent to the boundary cell by using the area of each boundary surface calculated in the boundary area calculation step It is preferable that these are further included.

  According to the above configuration, as a boundary cell including the boundary of the calculation target region, a hexahedron having at most one cutting point on each side is cut by N triangular cutting planes having three cutting points as vertices. A (6 + N) face body obtained as described above is used. Since the coefficients included in the algebraic equation to be solved are calculated using the area of the boundary surface of these boundary cells, the algebraic equation to be solved can be calculated with high accuracy. Therefore, the value of each physical quantity on each cell can be calculated with high accuracy.

  In the numerical fluid calculation method according to the present invention, it is preferable that the upper limit value Nmax of the number of triangular cut surfaces N included in each of the boundary cells is four.

  According to said structure, the division | segmentation required in order to calculate the area of the boundary surface of a boundary cell can be reduced significantly (refer FIGS. 5-8).

  The numerical fluid calculation apparatus according to the present invention may be realized by a computer. In this case, a program for realizing the numerical fluid computing device according to the present invention in the computer by operating the computer as each of the above means and a computer-readable recording medium recording the program also fall within the scope of the present invention. .

  As described above, according to the present invention, in the numerical fluid calculation method for solving the convection diffusion equation using the finite volume method, the numerical value obtained by reducing the number of necessary cases without replacing the cut cube with a cut cell or the like. A fluid calculation method can be realized.

1, showing an embodiment of the present invention, is a block diagram showing a configuration of a numerical fluid computing device. FIG. 1, showing an embodiment of the present invention, is a block diagram showing a hardware configuration of a numerical fluid computing device realized by using a computer. FIG. BRIEF DESCRIPTION OF THE DRAWINGS It is explanatory drawing which shows embodiment of this invention and shows the structure of a cutting | disconnection cube model. (A) is sectional drawing of a target object, (b) is a perspective view of an internal cell, (c) is a perspective view of a boundary cell. It is an explanatory view showing an embodiment of the present invention and showing classification of border cells. BRIEF DESCRIPTION OF THE DRAWINGS FIG. 1 illustrates an embodiment of the present invention and is a plan view illustrating a configuration of a boundary surface having no cutting point. FIG. 5 is a plan view showing an embodiment of the present invention and showing a configuration of a boundary surface including two cutting points. FIG. 4 is a plan view illustrating a configuration of a boundary surface including three cutting points according to the embodiment of the present invention. FIG. 5 is a plan view illustrating a configuration of a boundary surface including four cutting points according to the embodiment of the present invention. 1, showing an example of the present invention, is a schematic diagram showing a structure of a system to be calculated. FIG. 10 is a graph showing an evaluation result of average Nusselt in the system shown in FIG. 9. (A) shows the average Nusselt number on the small spherical surface, and (b) shows the average Nusselt number on the large spherical surface.

[Matters that form the basis of the present invention]
Prior to describing the embodiments of the present invention, matters serving as the basis of the present invention will be described. In the present specification, subscripts are represented by prefixing _, and superscripts are represented by prefixing ^. For example, a A _e represents a A_e, represents the L ² and L ^ 2.

(Cutting cube model)
A cut cube model used in the numerical fluid calculation method of the present invention will be described with reference to FIGS.

  The cut cube model is a modeling technique in the volume CAD system developed by the inventors. As shown in FIG. 3A, the hexahedral internal cell IC included in the object Ω and the object Ω This is a modeling technique for expressing the object Ω using a 6 + N-faced boundary cell BC including the boundary surface ∂Ω. Internal cells are sometimes referred to as “non-boundary cells”.

  In the cut cube model, a hexahedron (for example, a cube) as shown in FIG. 3B is used as the internal cell IC. In the internal cell IC, each of the six boundary surfaces forms a boundary with six adjacent cells (may be an internal cell or a boundary cell) adjacent to the internal cell IC. Further, the internal cell IC includes one lattice point. The position of the lattice point included in the internal cell IC is, for example, the position of the center of gravity of the internal cell IC.

  In this specification, a lattice point included in a cell of interest (hereinafter referred to as “cell of interest”) is denoted by P, and six boundary surfaces of the cell of interest are represented by e (x-axis positive direction) and w (x-axis negative). Direction), n (y-axis positive direction), s (y-axis negative direction), t (z-axis positive direction), b (z-axis negative direction). In addition, lattice points included in cells adjacent to the target cell (hereinafter referred to as “neighboring cells”) via the boundary surfaces e, w, s, n, t, and b are respectively E, W, S, and N. , T, B. In addition, the areas of the boundary surfaces e, w, s, n, t, and b are denoted as A_e, A_w, A_n, A_s, A_t, and A_b, respectively (see FIG. 3B).

  In the cut cube model, a (6 + N) face body as shown in FIG. 3C is used as the boundary cell BC. That is, a base hexahedron (for example, a cube congruent with the internal cell) having at most one cutting point on each side is cut by N triangular cutting planes having three cutting points as vertices. The (6 + N) face plane is used. The cutting point is provided at the intersection of each side of the base hexahedron and the boundary surface ∂Ω of the object Ω (if the size of each cell is sufficiently small, the number of cutting points provided on each side) Can be made 1 or less). The boundary cell BC illustrated in FIG. 3C is obtained by cutting along a triangular cut surface kt_1 having cut points K_1, K_2, and K_3 as vertices and a triangular cut surface kt_2 having cut points K_1, K_2, and K_4 as vertices. It is what was done.

  In the boundary cell BC, each of the six boundary surfaces excluding the triangular cut surface KT defines a boundary with six adjacent cells (which may be internal cells or boundary cells) adjacent to the boundary cell BC. Make it. On the other hand, the triangular cut surface KT constitutes a boundary between the object Ω and the outside world (solid). In addition, the boundary cell BC includes one lattice point, like the internal cell IC. The position of the grid point is, for example, the barycentric position of the boundary cell BC.

  It should be noted that the number N of triangular cut planes can be different for each boundary cell. However, in the present invention, an upper limit is set for the number N of cut surfaces of the boundary cell. In other words, the cell size is set to be sufficiently small so that no border cell is generated in which the number of cut surfaces N exceeds a predetermined upper limit value. In the present embodiment, the upper limit of the number N of cut surfaces of the boundary cell is four. In other words, the cell size is set to be sufficiently small so that no border cell having a cutting surface number N exceeding 4 is generated. Thereby, each boundary cell is classified into one of the four types shown in FIG.

(Convection diffusion equation)
The equation targeted by the present invention is the convection diffusion equation (1) followed by the physical quantity φ = φ (t, x, y, z) defined on the object Ω. The convection diffusion equation (1) may be referred to as “convection diffusion equation” or “transport equation”. As the physical quantity φ, each component of momentum per unit mass or energy per unit mass is assumed, but is not limited thereto.

  Here, U = (u, v, w) is a flow velocity (vector field), ρ is a density (scalar field), Γ is a diffusion coefficient (scalar field), and S is a generation term (scalar field). ρUφ is called convective flux, and Γ∇φ is called diffusion flux.

  In the following description, for the sake of simplicity, it is assumed that the density ρ and the diffusion coefficient Γ are constant regardless of time and position, but the present invention is not limited to this. That is, the density ρ and the diffusion coefficient Γ may depend on time or position. In the following description, for the sake of simplicity, it is assumed that the generated term S can be expressed by the sum of the term S_P × φ proportional to the physical quantity φ and the constant term S_C. However, the present invention is not limited to this. Absent. That is, the generation term S may include a second or higher order term of the physical quantity φ.

  In the numerical fluid calculation method of the present invention, the convection diffusion equation (1) is solved using a finite volume method. That is, an equation system (simultaneous equations) consisting of algebraic equations obtained by integrating the convection diffusion equation (1) on each cell is solved by a computer. Each algebraic equation is an algebraic equation that is established between a physical quantity φ_P at a lattice point P included in a target cell and a physical quantity φ_NB at a lattice point NBε {E, W, N, S, T, B} included in an adjacent cell. It is. These algebraic equations can be derived from the convection diffusion equation (1) as follows.

  That is, first, the convection diffusion equation (1) is integrated on each internal cell to obtain the equation (2). Further, the convection diffusion equation (1) is integrated on each boundary cell to obtain the equation (3).

  Here, Δt is a time step, and ΔV is the volume of the cell of interest. A_nb is the area of each boundary surface nbε {e, w, n, s, t, b} of the cell of interest, and A_kt is each triangular cut surface ktε {kt_1, kt_2,. kt_N} area. In the expressions (2) and (3), the area A_nb is treated as an area vector facing the outward normal direction of the boundary surface nb. In the expression (3), the area A_KT is treated as an area vector facing the outward normal direction of the triangular cut surface kt.

  Further, (ρφ) _P is a value of ρφ at the lattice point P included in the target cell. (ΡUφ) _nb and (Γ∇φ) _nb are respectively the convective flux ρUφ and the diffusion flux Γ at each boundary surface nbε {e, w, n, s, t, b} of the target cell. The value of ∇φ. The sum Σ_nb crosses all six boundary surfaces of the cell of interest.

  Further, (ρUφ) _kt and (Γ∇φ) _kt are the convective flux ρUφ and the diffusion flow at the triangular cut surfaces kt∈ {kt_1, kt_2,..., Kt_N} (representative points) of the cell of interest. The value of the bundle Γ∇φ. The sum Σ_kt crosses all N triangular cut surfaces of the cell of interest.

  Next, using a suitable difference method, the value φ_nb at the boundary surface nbε {e, w, n, s, t, b} of the physical quantity φ is converted into the lattice point NB∈ {E, W, N, The algebraic equation (4) is obtained by substituting the value φ_NB in S, T, B}.

Hybrid Method When the hybrid method is used as the difference method, the coefficient a_NB in the algebraic equation (4) for the internal cell is given by the equations (5a) to (5h). The constant b in the algebraic equation (4) for the internal cell is given by the equation (5i).

  Here, the parentheses || · || indicate the largest element among the elements in the parentheses. F_nb is an amount defined by the equation (6), and D_nb is an amount defined by the equation (7). In Expression (7), r_ (P → NB) represents a vector having the lattice point P of the target cell as the start point and the lattice point NB of the adjacent cell as the end point. e_ (P → NB) is a unit vector having the same direction as r_ (P → NB).

  On the other hand, when the hybrid method is used as the difference method, the coefficient a_NB in the algebraic equation (4) for the boundary cell is given by the equations (5a) to (5h) as in the algebraic equation (4) for the inner cell. The constant b in the algebraic equation (4) for the boundary cell is given by the equation (8) including the integral of the convective flux ρUφ and the diffusion flow velocity Γ∇φ on the triangular cut surface.

  The integral of the convection flux ρUφ and the diffusion flow velocity Γ∇φ on the triangular cut surface can be calculated from the values of φ and ∇φ given as boundary conditions on the boundary surface ∂Ω.

Upwind Difference Method When the upwind difference method is used as the difference method, the coefficient a_NB in the algebraic equation (4) for the internal cell is given by the equations (9a) to (9b). The constant b in the algebraic equation (4) for the internal cell is given by the equation (9c). F_nb is an amount defined by equation (6) described above, and D_nb is an amount defined by equation (7).

  Here, the upwind difference method means that the value φ_nb at the boundary surface nbε {e, w, n, s, t, b} of the physical quantity φ is represented by the lattice point NB∈ {E, It means replacing with the value φ_NB in W, N, S, T, B}. In Expression (10), Δr_ (P → NB) represents a vector starting from the lattice point P of the target cell and ending at the lattice point NB of the adjacent cell, and Δr_ (NB → P) is a lattice of the adjacent cell. This represents a vector starting from the point NB and ending at the lattice point P of the cell of interest. M_nb ^ + and M_nb ^-are defined by the equation (11).

  When the integration of the time term in the target cell in equation (1) is discretized as in equation (12), the coefficient a_t in equation (9b) and the constant b_temporal in equation (9c) ).

  Further, the constant b_CT and the constant b_DT in the equation (9c) are given as in the equations (14a) to (14b).

  On the other hand, when the upwind difference method is used as the difference method, the coefficient a_NB in the algebraic equation (4) for the boundary cell is given by the equations (9a) to (9b) similarly to the algebraic equation (4) for the inner cell. . Also, the constant b in the algebraic equation (4) for the boundary cell is given by the equation (15) including the integral of the convective flux ρUφ and the diffusion flow velocity Γ∇φ on the triangular cut surface.

  Here, the hybrid method and the finite volume method using the algebraic equation (4) derived from the upwind difference method have been described, but the scope of application of the present invention is not limited to this. That is, the finite volume method using the algebraic equation (4) derived from other discretization methods such as the central difference method, the power method, and the exponent method (exact method) is also included in the scope of application of the present invention.

  Further, here, the flow velocity U is known and the transport of the physical quantity φ in the flow field defined by the flow velocity U has been described, but the present invention is not limited to this. That is, the numerical fluid calculation method for calculating the flow velocity U * of the unsteady incompressible viscous fluid by combining the algebraic equation (4) and the pressure correction equation is also included in the scope of the present invention.

Embodiment
Hereinafter, an embodiment of the present invention will be described with reference to the drawings.

(Configuration of numerical fluid computing device)
The configuration of the numerical fluid calculation apparatus 1 according to the present embodiment will be described with reference to FIG. FIG. 1 is a block diagram showing the configuration of the numerical fluid calculation apparatus 1. The numerical fluid computing device 1 includes a modeling unit 11, a cut surface count unit 12 (abbreviated as “count unit” in the figure), a boundary cell information generation unit 13, an internal cell information generation unit 14, a boundary cell coefficient A setting unit 15, an internal cell coefficient setting unit 16, a matrix calculation unit 17, a convergence determination unit 18, a UI unit 19, a parameter storage unit 21, and a mesh data storage unit 22 are provided.

  The modeling unit 11 is means for dividing the object Ω into a plurality of cells. Specifically, the object Ω is divided into a plurality of cells according to the cut cube model, and mesh data representing the vertex coordinates of each cell is generated. As described above, the modeling unit 11 uses a hexahedron as an internal cell and uses a (6 + N) face as a boundary cell. The mesh data generated by the modeling unit 11 is provided to the cut surface counting unit 12.

  The cut surface number counting unit 12 refers to the mesh data provided from the modeling unit 11, counts the cut surface number N in each boundary cell, and determines whether the maximum cut surface number is 4 or less.

  When the maximum number of cut surfaces is 4 or less, the cut surface number counting unit 12 returns “OK” to the modeling unit 11. When “OK” is returned from the cut surface count unit 12, the modeling unit 11 stores mesh data representing the vertex coordinates of each cell in the mesh data storage unit 22.

  On the other hand, when the maximum number of cut surfaces is larger than 4, the cut surface number counting unit 12 returns “NG” to the modeling unit 11. When “NG” is returned from the cut surface count unit 12, the modeling unit 11 subdivides the object Ω into smaller-sized internal cells and boundary cells. The modeling unit 11 repeats subdivision until the maximum number of cut surfaces becomes 4 or less (until “OK” is returned from the cut surface number counting unit 12), and finally the mesh whose maximum number of cut surfaces becomes 4 or less. Data is stored in the mesh data storage unit 22.

  The boundary cell information generation unit 13 determines (1) volume ΔV, (2) coordinates of lattice point P, and (3) lattice point NBε {E for each boundary cell from the mesh data stored in the mesh data storage unit 22. , W, N, S, T, B} and the grid point P | NB-P |, (4) coordinates of representative points of each boundary surface nbε {e, w, n, s, t, b} (5) area A_nb of each boundary surface nbε {e, w, n, s, t, b}, (6) coordinates of representative points of each triangular cut surface ktε {kt_1, kt_2, ..., kt_N}, (7) The area A_kt of each triangular cut surface ktε {kt_1, kt_2,..., Kt_N} and (8) the normal vector n_kt of each triangular cut surface ktε {kt_1, kt_2,. Various quantities calculated by the boundary cell information generation unit 13 are provided to the boundary cell coefficient setting unit 15 as boundary cell information.

  Note that the coordinates of the lattice point P of the boundary cell can be easily calculated by, for example, dividing the sum of the vertex coordinates of the boundary cell by the number of vertices. The coordinates of the representative point of each boundary surface nb can be easily calculated by dividing the sum of the vertex coordinates of the boundary surface nb by the number of vertices, for example. The area A_nb of each boundary surface nb can be calculated using a quadrature formula (described later) derived for each configuration of the boundary surface. The same applies to the volume ΔV of the boundary cell.

  Further, the coordinates of the representative point of each triangular cut surface kt can be easily calculated by dividing the sum of the vertex coordinates of the triangular cut surface kt by the number of vertices 3, for example. The area A_kt of each triangular cut surface kt is, for example, A_kt = | (β−α) × (γ−, where α, β, and γ are the position vectors of three vertices (cut points) of the triangular cut surface kt. α) | / 2 (“×” is an outer product). The normal vector n_kt of the triangular cut surface kt can be calculated by, for example, n_kt = (β−α) × (γ−α) / | (β−α) × (γ−α) | X is an outer product).

  From the mesh data stored in the mesh data storage unit 22, the internal cell information generation unit 14 determines (1) the volume ΔV, (2) the coordinates of the lattice point P, and (3) each lattice point NBε { E, W, N, S, T, B} and grid point P distance | NB-P |, (4) of representative points of each boundary surface nbε {e, w, n, s, t, b} Coordinates (5) The area A_nb of each boundary surface nbε {e, w, n, s, t, b} is calculated. The various amounts calculated by the internal cell information generation unit 14 are provided to the internal cell coefficient setting unit 16 as internal cell information.

  Note that the coordinates of the lattice point P of the internal cell can be easily calculated by dividing the sum of the vertex coordinates of the internal cell by the number of vertices 8, for example. The coordinates of the representative point of each boundary surface nb can be easily calculated by dividing the sum of the vertex coordinates of the boundary surface nb by the number of vertices 4, for example. Further, the area A_nb of each boundary surface nb can be easily calculated by multiplying the lengths of two adjacent sides, for example. For example, if the internal cell is a cube with the length of each side being L, Amb = L ^ 2 (constant). Further, the volume ΔV of the internal cell can be easily calculated by multiplying the lengths of the three adjacent sides, for example. For example, if the internal cell is a cube with side length L, ΔV = L ^ 3 (constant).

  The values of parameters {Δt, ρ, Γ, S, U = (u, v, w), φ_kt, (∇φ) _kt} necessary to construct the equation system consisting of the algebraic equation (4) It is provided by the user via 19 and stored in the parameter storage unit 21. The boundary cell coefficient setting unit 15 and the internal cell coefficient setting unit 16 can freely refer to these parameters stored in the parameter storage unit 21.

  The boundary cell coefficient setting unit 15 calculates the integration sum Σ (ρUφ) _kt · A_kt of the convection flux ρUφ on each triangular cut surface kt and the integral on each triangular cut surface kt of the diffusion flow velocity Γ∇φ. The sum is calculated with reference to the boundary values φ_kt and (Δφ) _kt stored in the parameter storage unit 21.

  Specifically, (1) the density ρ, the flow velocity U_kt, and the boundary condition φ_kt are read from the parameter storage unit 21, (2) the normal vector n_kt and the area A_kt are acquired from the boundary cell information generation unit 13, and (3 ) Calculate inner product U_kt · n_kt between flow velocity U_kt and normal vector n_kt, and (4) calculate (ρUφ) _kt · A_kt by multiplying inner product U_kt · n_kt by density ρ, boundary condition φ_kt, and area A_kt. To do. Then, Σ (ρUφ) _kt · A_kt is obtained by adding the obtained (ρUφ) _kt · A_kt.

  Similarly, (1) the diffusion coefficient Γ and the boundary condition (∇φ) _kt are read from the parameter storage unit 21, (2) the normal vector n_kt and the area A_kt are acquired from the boundary cell information generation unit 13, (3 ) Calculate the inner product (∇φ) _kt · n_kt of the boundary condition (∇φ) _kt and the normal vector n_kt, and (4) multiply the inner product (∇φ) _kt · n_kt by the diffusion coefficient Γ and area A_kt, (Γ∇φ) _kt · A_kt is calculated. Then, Σ (Γ∇φ) _kt · A_kt is obtained by adding the obtained (Γ∇φ) _kt · A_kt.

  The boundary cell coefficient setting unit 15 further includes {ρ, U, Γ, S_P} read from the parameter storage unit 21 and {A_nb, | NB−P |, ΔV} acquired from the boundary cell information generation unit 13. The coefficients a_P, a_E, a_W, a_N, a_S, a_T, and a_B are calculated according to equations (5a) to (5h). Further, {Σ (Γ∇φ) _kt · A_kt, Σ (ρUφ) _kt · A_kt} calculated in advance, {Δt, S_C, ρ, φ_P ^ T} read from the parameter storage unit 21, and boundary cell information With reference to ΔV acquired from the generation unit 13, the constant b is calculated according to the equation (8).

  The internal cell coefficient setting unit 16 refers to {ρ, U, Γ, S_P} read from the parameter storage unit 21 and {A_nb, | NB−P |, ΔV} acquired from the boundary cell information generation unit 13. The coefficients a_P, a_E, a_W, a_N, a_S, a_T, and a_B are calculated according to the equations (5a) to (5h). Further, by referring to {Δt, S_C, ρ, φ_P ^ T} read from the parameter storage unit 21 and ΔV acquired from the border cell information generation unit 13, the constant b is calculated according to the equation (5i).

  Note that φ_P ^ t referred to by the boundary cell coefficient setting unit 15 and the internal cell coefficient setting unit 16 to calculate the constant b is not set by the user via the UI unit 19, but is set by the matrix calculation unit 17. It is φ_P calculated last time. However, φ_P ^ t referred to in the first calculation is an initial value φ_P ^ 0 set by the user via the UI unit 19.

  The matrix calculation unit 17 is means for solving the algebraic equation (4) defined by the coefficient a_NB and the constant b set by the boundary cell coefficient setting unit 15 and the internal cell coefficient setting unit 16. As algorithms for solving the algebraic equation (4), the Jacobian iteration method, the Gauss-Seidel iteration method, the successive over relaxation method (SOR), and the like are known. Use sequential overrelaxation method.

  The matrix calculation unit 17 provides the updated solution to the convergence determination unit 18 every time the solution of the algebraic equation (4) (the value of the physical quantity φ at each lattice point) is updated. The convergence determination unit 18 determines whether or not the difference (norm) between the pre-update solution and the post-update solution is equal to or less than a predetermined threshold.

  When the difference between the solution before the update and the solution after the update is equal to or less than a predetermined threshold, the convergence determination unit 18 returns “OK” to the matrix calculation unit 17. When “OK” is returned from the convergence determination unit 18, the matrix calculation unit 17 outputs the updated solution to the outside and stores the updated solution in the parameter storage unit 21. The solution stored in the parameter storage unit 21 is referred to as the solution φ_P ^ t at time t in order to calculate the physical quantity φ_P ^ t + Δt at time t + Δt.

  If the difference between the pre-update solution and the post-update solution is greater than a predetermined threshold, the convergence determination unit 18 returns “NG” to the matrix calculation unit 17. When “NG” is returned from the convergence determination unit 18, the matrix calculation unit 17 updates the solution again. The matrix calculation unit 17 repeatedly updates the solution until the difference between the solution before the update and the solution after the update is equal to or less than a predetermined threshold (until “OK” is returned from the convergence determination unit 18), and finally A convergence solution is obtained in which the difference from the solution before update is equal to or less than a predetermined threshold.

(Configuration example using a computer)
The numerical fluid computing device 1 can be realized using, for example, a computer (electronic computer). FIG. 2 is a block diagram illustrating a hardware configuration of the numerical fluid calculation apparatus 1 realized using a computer.

  As shown in FIG. 2, the numerical fluid calculation apparatus 1 includes an arithmetic device 120, a main storage device 130, an auxiliary storage device 140, and an input / output interface 150 connected to each other via a bus 110. . An example of a device that can be used as the arithmetic device 120 is a CPU (Central Processing Unit). Further, as a device that can be used as the main storage device 130, for example, a semiconductor RAM (random access memory) can be cited. An example of a device that can be used as the auxiliary storage device 140 is a hard disk drive.

  As shown in FIG. 2, the input device 200 and the output device 300 are connected to the input / output interface 150. The input device 200 is a means for the user to input the values of the parameters {Δt, ρ, Γ, S, U = (u, v, w), φ_kt, (∇φ) _kt} to the numerical fluid calculation device 1. Yes, for example, a keyboard. The computing device 120 can refer to the parameters {Δt, ρ, Γ, S, U = (u, v, w), φ_kt, (を φ) _kt} input via the input device 200. Stored in the main storage device 130. That is, the main storage device 130 is used as the parameter storage unit 21.

  On the other hand, the output device 300 is a means for outputting the solution of the algebraic equation (4) (the value of the physical quantity φ at each lattice point), for example, a monitor. Instead of outputting the solution of the algebraic equation (4) via the output device 300, the solution of the algebraic equation (4) may be stored in the auxiliary storage device 140.

  The auxiliary storage device 140 stores a numerical fluid calculation program for operating the computer as the numerical fluid calculation device 1. The numerical fluid calculation program includes a modeling program, a cutting surface count program, a boundary cell information setting program, a boundary cell coefficient setting program, an internal cell coefficient setting program, a matrix calculation program, a convergence determination program, and a UI program. Including. The auxiliary storage device 140 is also used as the mesh data storage unit 22 that stores mesh data.

  The arithmetic device 120 causes the computer to function as the modeling unit 11 by executing instructions included in the modeling program developed on the main storage device 130. Similarly, included in the cutting plane number count program, the boundary cell information setting program, the boundary cell coefficient setting program, the internal cell coefficient setting program, the matrix operation program, the convergence determination program, and the UI program developed on the main storage device 130 By executing the instructions, the computer 1 can be connected to the cutting plane number counting unit 12, the boundary cell information generation unit 13, the internal cell information generation unit 14, the boundary cell coefficient setting unit 15, the internal cell coefficient setting unit 16, and the matrix calculation unit 17. , And function as the convergence determination unit 18.

  An object of the present invention is to supply a numerical fluid calculation apparatus 1 with a recording medium in which program codes (execution format program, intermediate code program, source program) of each of the above programs are recorded so as to be readable by a computer. This can also be achieved by reading and executing the program code recorded on the recording medium.

  Examples of the recording medium include a tape system such as a magnetic tape and a cassette tape, a magnetic disk such as a floppy (registered trademark) disk / hard disk, and an optical disk such as a CD-ROM / MO / MD / DVD / CD-R. Card system such as IC card, IC card (including memory card) / optical card, or semiconductor memory system such as mask ROM / EPROM / EEPROM / flash ROM.

  Further, the numerical fluid calculation apparatus 1 may be configured to be connectable to a communication network, and the program code may be supplied to the numerical fluid calculation apparatus 1 via the communication network. The communication network is not particularly limited. For example, the Internet, intranet, extranet, LAN, ISDN, VAN, CATV communication network, virtual private network, telephone line network, mobile communication network, satellite communication. A net or the like is available. Further, the transmission medium constituting the communication network is not particularly limited. For example, even in the case of wired such as IEEE 1394, USB, power line carrier, cable TV line, telephone line, ADSL line, etc., infrared rays such as IrDA or remote control, Bluetooth ( (Registered trademark), 802.11 wireless, HDR, mobile phone network, satellite line, terrestrial digital network, and the like can also be used. The present invention can also be realized in the form of a computer data signal embedded in a carrier wave in which the program code is embodied by electronic transmission.

(Quadrature formula)
As described above, the area A_nb of the boundary surface nb in the boundary cell can be calculated using the quadrature formula derived for each configuration of the boundary surface. Hereinafter, this quadrature formula will be described with reference to FIGS. In the following description, the ratio of the boundary surface nb occupying the surface of the basic cube is referred to as an aperture ratio A_fluid.

  From the condition that “at most one cutting point is provided on each side of the base cube”, the combinations to be considered are (1) when there is no cutting point, (2) when there are two cutting points, (3) When there are three cut points, it is classified as either (4) there are four cut points. Hereinafter, the cases (1) to (4) will be described in order.

(1) When there is no cutting point When there is no cutting point on the boundary surface nb, the entire boundary surface nb is included in the inside (fluid) of the object Ω (see FIG. 5A) or the boundary surface nb Either the whole is included in the outside (solid) of the object Ω (see FIG. 5B). In the former case, the aperture ratio A_fluid = 1, and in the latter case, the aperture ratio A_fluid = 0.

(2) When there are two cutting points When the boundary surface nb includes two cutting surfaces, the configuration of the boundary surface nb is one of the three types shown in FIGS. 6 (a) to 6 (c). Become. In this case, the aperture ratio A_fluid is calculated by Expression (16).

(3) When there are three cutting points When there are three cutting points on the boundary surface nb, the configuration of the boundary surface nb is one of the five types shown in FIGS. 7 (a) to 7 (e). become. In the case of the configuration shown in FIGS. 7A to 7D, the aperture ratio A_fluid is calculated by Expression (17).

  On the other hand, the configuration of FIG. 7E is not desirable because the problem is degraded in two dimensions. For this reason, when the arrangement of the cutting points shown in FIG. 7E occurs, it is desirable to recut the mesh (subdivide the object Ω).

(4) When there are four cutting points When there are four cutting points on the boundary surface nb, the configuration of the boundary surface nb is one of the four types shown in FIGS. 8 (a) to 8 (d). become. In the case of the configuration shown in FIGS. 8A to 8B, the aperture ratio A_fluid is calculated by Expression (18).

  Note that the configuration shown in FIGS. 8C to 8D is not desirable because the problem is degraded two-dimensionally. For this reason, when arrangement | positioning of the cutting | disconnection point shown in FIG.8 (c)-FIG.8 (d) arises, it is desirable to recut a mesh (subdividing object Ω).

〔Example〕
Finally, an embodiment in which the present invention is applied to the convection diffusion equation according to the temperature T will be described.

  A system to be calculated in the present embodiment is shown in FIG. The system to be calculated in this embodiment is a space between a small spherical surface (radius d_i) and a large spherical surface (radius d_o) having a common center (d_i <d_o), as shown in FIG. As boundary conditions, the temperature of the small spherical surface (inner boundary) = T_i and the temperature of the large spherical surface (outer boundary) = T_o are imposed (T_i> T_o).

  After setting the Rayleigh number Ra to each value in the table below, the value of the temperature T in each cell was calculated according to the numerical fluid calculation method of the present invention. The Rayleigh number Ra is an amount defined by the equation (19).

  In order to check the calculation accuracy, the average Nusselt number (hereinafter abbreviated as “Nu number”) on each spherical surface was evaluated. Note that the average Nusselt number on the small sphere is an amount defined as in Expression (20a), and the average Nusselt number on the large sphere is an amount defined as in Expression (20b).

  The evaluation results are shown in FIG. FIG. 10A shows the average Nusselt number on the small spherical surface, and FIG. 10B shows the average Nusselt number on the large spherical surface. In the figure, the average Nusselt number described in the following reference is shown as a comparative example. According to the numerical fluid calculation method of the present invention, it can be seen that highly accurate numerical fluid calculation is realized for any Rayleigh number shown in Table 1.

  [References] Greg V K, "Natural convection between concentric spheres", International Journal of Heat and Mass Transfer, 35 (8), 1991, pp1935-1945

  The present invention can be widely used for numerical fluid calculations. In particular, (V) it can be suitably used for numerical fluid calculation in cooperation with a CAD system.

DESCRIPTION OF SYMBOLS 1 Computational fluid calculation apparatus 11 Modeling part 12 Section number counting part 13 Boundary cell information generation part 14 Internal cell information generation part 15 Boundary cell coefficient setting part 16 Internal cell coefficient setting part 17 Matrix calculation part 18 Convergence determination part 19 UI part 21 Parameter storage unit 22 Mesh data storage unit IC internal cell BC boundary cell e, w, n, s, t, b boundary surface kt triangular cut surface P lattice point included in the target cell E, W, N, S, T, B Lattice points included in neighboring cells

Claims (7)

An algebraic equation obtained by integrating the convection diffusion equation on each cell, and an equation system comprising an algebraic equation that holds between a physical quantity on each cell and the physical quantity on an adjacent cell adjacent to the cell. in computational fluid calculation method by a computer solved,
A cell dividing step for dividing a calculation target area into a plurality of cells, and a hexahedron having at most one cutting point on each side as a boundary cell including the boundary of the calculation target area is apex with three cutting points A cell division step using a (6 + N) face (N is a natural number equal to or less than Nmax) obtained by cutting with N triangular cut planes;
For each of the boundary cells, a cutting area calculating step for calculating an area of the N triangular cutting planes;
For each of the boundary cells, the integrated value on each triangular cut surface of the convection flux and the integrated value on each triangular cut surface of the diffusion flux included in the algebraic equation are calculated in the cut area calculation step. And a boundary condition setting step for calculating from the area of each triangular cut surface calculated in (1) and the boundary condition imposed on the boundary of the calculation target region,
A numerical fluid calculation method characterized by the above. For each of the boundary cells, a boundary area calculation step for calculating an area of a boundary surface between adjacent cells adjacent to the boundary cell;
For each of the boundary cells, a coefficient included in the algebraic equation, the coefficient relating to the physical quantity on an adjacent cell adjacent to the boundary cell, the area of each boundary surface calculated in the boundary area calculation step A coefficient calculation step of calculating using
The numerical fluid calculation method according to claim 1, wherein: The upper limit Nmax of the triangular cut surface number N included in each of the upper Kisakai field cell is 4,
The numerical fluid calculation method according to claim 1 or 2, characterized in that: The convective diffusion equation, the physical quantity phi, the density [rho, flow velocity U = (u, v, w), the diffusion coefficient gamma, a product term as _{S = S P φ + S C} , is represented by the following formula (1) ,
The algebraic equation expresses the physical quantity at the lattice point P included in the cell as φ _P and the adjacent cell adjacent to the cell via the boundary surface nb (nb = e, w, n, s, t, b). The physical quantity at the included grid point NB (NB = E, W, N, S, T, B) is represented by the following equation (2), where φ _NB is
The area of the boundary surface nb A _nb, the flow velocity at the boundary surface nb (u _nb, v _nb, w _nb) and F _nb defined by the following equation (3) as, the diffusion coefficient in the boundary surface nb Γ _nb , the area vector of the boundary surface nb facing the outward normal direction is A _nb , the vector starting from the lattice point P and ending at the lattice point NB is r _{P → NB} , and the vector r _{P → NB} And D _nb defined by the following equation (4) where e _{P → NB} is the unit vector in the same direction as :
The coefficients a _NB and a _P included in the algebraic equation are expressed by the following equation (5), where the volume of the cell is ΔV,
When the cell is an internal cell, the coefficient b included in the algebraic equation is set as ΔV for the volume of the cell, Δt for the time step, or an initial value, or the previously calculated lattice point P The physical quantity is represented by the following equation (6), where φ _P ^t :
When the cell is a boundary cell, the coefficient b included in the algebraic equation is set as ΔV for the volume of the cell, Δt for a time step, or an initial value, or the lattice point P calculated at the previous time. The physical quantity is φ _P ^t , the area vector of the triangular cut surface kt (kt = kt ₁ , kt ₂ ,..., Kt _N ) facing the outward normal direction is A _kt , and the convection flow velocity at the triangular cut surface kt is (ρUφ ) _Kt , where the diffusion flow velocity at the triangular cut surface kt is (Γ∇φ) _{kt and} is expressed by the following equation (7)
The coefficient calculated in the coefficient calculation step is a coefficient a _NB (NB = E, W, N, S, T, B) included in the algebraic equation when the cell is a boundary cell .
The numerical fluid calculation method according to claim 2 , wherein:

An algebraic equation obtained by integrating the convection diffusion equation on each cell, and an equation system comprising an algebraic equation that holds between a physical quantity on each cell and the physical quantity on an adjacent cell adjacent to the cell. In the numerical fluid computing device to solve,
A cell dividing means for dividing a calculation target area into a plurality of cells, and as a boundary cell including a boundary of the calculation target area, a hexahedron having at most one cutting point on each side is apexed at three cutting points. Cell dividing means using a (6 + N) plane (N is a natural number equal to or less than Nmax) obtained by cutting with N triangular cut planes,
Cutting area calculating means for calculating the area of the N triangular cutting planes for each of the boundary cells;
For each of the boundary cells, the integrated value of each of the convective fluxes on each triangular cut surface and the integrated value of each of the diffusion fluxes on each triangular cut surface included in the algebraic equation are calculated as the cut area calculation means. Boundary condition setting means for calculating from the area of each triangular cut surface calculated in the above and the boundary condition imposed on the boundary of the calculation target region,
A numerical fluid computing device characterized by the above.   A program for causing a computer to operate as the numerical fluid calculation apparatus according to claim 5, wherein the computer functions as each unit included in the numerical fluid calculation apparatus.   A computer-readable recording medium on which the program according to claim 6 is recorded.

Images