Simulation device, simulation method, and program

Abstract

To provide a simulation device that eliminates the need for setting spatial resolution for each region within an analysis space as pretreatment analysis and can reduce a calculation load.SOLUTION: Fluid is expressed by multiple calculation particles, a physical quantity is imparted to each of the multiple calculation particles, the plurality of calculation particles are disposed in an analysis space obtained by modeling a flow path, and a physical quantity imparted to each of the multiple calculation particles and positions of the multiple calculation particles are developed temporarily by solving a governing equation of fluid. In processing for temporarily developing a physical quantity imparted to each of multiple calculation particles and positions of the multiple calculation particles, processing for comparing a characteristic dimension of a flow path in respective positions of the multiple calculation particles with a size of each of the multiple calculation particles, dividing each of the multiple calculation particles according to a comparison result, and generating smaller calculation particles, or processing for uniting some of the multiple calculation particles and generating larger calculation particles is performed.SELECTED DRAWING: Figure 2

Description

[Technical field]

[0001] The present invention relates to a simulation device, a simulation method, and a program. [Background technology]

[0002] A known fluid analysis technique is the particle method, which analyzes the time evolution of a flow field by discretizing a continuous field using a set of virtual particles (hereafter referred to as computational particles) and solving the governing equations. In this analysis, the spatial resolution depends on the diameter of the computational particles. In conventional techniques, the diameter of the computational particles is set in the pre-processing stage of the analysis, and multiple computational particles with the same diameter are placed in the analysis space to perform the analysis.

[0003] A method is known in which the diameter of a computational particle used in an analysis is set or a shape model is modified based on a given shape model and the minimum value of the diameter of the computational particle (Patent Document 1). In order to obtain sufficient spatial resolution, it is necessary to set the diameter of the computational particle according to the narrowest region of the flow path. When arranging computational particles having a single diameter in the analysis space, even if the narrow region of the flow path is small relative to the entire analysis space, the diameter of the computational particle must be set according to the narrow region of the flow path in the entire analysis space. In this case, in the region where the flow path is wide, the spatial resolution becomes excessive, and the analysis load increases.

[0004] A method is known in which an analysis space is divided into a plurality of regions and calculation particles having different diameters are arranged in each divided region (Patent Document 2). Also, a method is known in which an analysis space is divided into a plurality of regions according to the distance from a wall and calculation particles having different diameters are arranged in each divided region (Patent Document 3). [Prior art documents] [Patent documents]

[0005] [Patent Document 1] JP 2017-134479 A [Patent Document 2] JP 2007-140993 A [Patent Document 3] JP 2022-108658 A Summary of the Invention [Problem to be solved by the invention]

[0006] In the conventional method, the operator needs to know in advance which areas require high spatial resolution and which areas do not, which increases the man-hours required for pre-processing of the analysis. In addition, depending on the skill of the operator, the spatial resolution may be set inappropriately. The object of the present invention is to provide a simulation device, a simulation method, and a program that can reduce the calculation load and eliminate the need to set the spatial resolution for each area in the analysis space as pre-processing of the analysis. [Means for solving the problem]

[0007] According to one aspect of the present invention, A simulation device that analyzes a fluid flow inside a flow path using a particle method, comprising: an input section for inputting simulation conditions; a processing unit that performs calculations based on the simulation conditions inputted to the input unit; Equipped with The processing unit includes: a function of expressing a fluid by a plurality of computational particles, assigning a physical quantity to each of the plurality of computational particles, arranging the plurality of computational particles in an analysis space in which the flow path is modeled, and evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time by solving a governing equation for the fluid; a function of performing a process of comparing a characteristic dimension of the flow path at the position of each of the plurality of computational particles with a size of each of the plurality of computational particles in a process of evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time, and performing a process of dividing each of the plurality of computational particles to generate smaller computational particles or a process of combining some of the plurality of computational particles to generate larger computational particles according to a comparison result; A simulation device is provided having:

[0008] According to another aspect of the invention, A simulation method for analyzing a fluid flow inside a flow path using a particle method, comprising the steps of: A computer given the simulation conditions expressing a fluid using a plurality of computational particles based on given simulation conditions, assigning a physical quantity to each of the plurality of computational particles, arranging the plurality of computational particles in an analysis space in which the flow path is modeled, and evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time by solving a governing equation for the fluid; A simulation method is provided in which, in a process of evolving the physical quantities assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time, a characteristic dimension of the flow path at the positions of each of the plurality of computational particles is compared with the size of each of the plurality of computational particles, and, depending on the comparison result, each of the plurality of computational particles is divided to generate smaller computational particles, or some of the plurality of computational particles are combined to generate larger computational particles.

[0009] According to another aspect of the invention, A program for causing a computer to realize a function of analyzing a fluid flow inside a flow path using a particle method, a function of expressing a fluid as a plurality of computational particles based on given simulation conditions, assigning a physical quantity to each of the plurality of computational particles, arranging the plurality of computational particles in an analysis space in which the flow path is modeled, and evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time by solving a governing equation for the fluid; a function of performing a process of comparing a characteristic dimension of the flow path at the position of each of the plurality of computational particles with a size of each of the plurality of computational particles in a process of evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time, and performing a process of dividing each of the plurality of computational particles to generate smaller computational particles or a process of combining some of the plurality of computational particles to generate larger computational particles according to a comparison result; A program for causing a computer to realize the above is provided. Effect of the Invention

[0010] Since the computational particles are combined or split according to the characteristic dimensions of the flow path, there is no need to set the spatial resolution for each region in the analysis space as a preprocessing step for the analysis. Furthermore, when the computational particles are combined, the computational load can be reduced. [Brief description of the drawings]

[0011] [Figure 1] FIG. 1A is a cross-sectional view showing an example of the shape of a flow path through which the fluid to be analyzed flows, FIG. 1B is a graph showing the absolute value of a signed distance function in shading, FIG. 1C is a schematic diagram for explaining a method for determining a characteristic dimension at an arbitrary position in the flow path, and FIG. 1D is a graph showing the distribution of the calculated characteristic dimensions in shading. [Diagram 2] FIG. 2 is a schematic diagram showing an example of the merging and division of computational particles during the analysis of a flow. [Diagram 3] FIG. 3A is a schematic diagram of a computational particle before division, and FIG. 3B is a schematic diagram of two computational particles generated by dividing one computational particle. [Figure 4]FIG. 4A is a schematic diagram of four computational particles generated by dividing one computational particle, and FIG. 4B is a schematic diagram of eight computational particles generated by dividing one computational particle. [Diagram 5] FIG. 5A is a schematic diagram showing an example of the distribution of a plurality of computational particles, and FIG. 5B is a schematic diagram showing an example of the distribution of computational particles after two computational particles are merged. [Figure 6] FIG. 6 is a block diagram of a simulation device according to the embodiment. [Figure 7] FIG. 7 is a flowchart showing the procedure of the simulation method according to the embodiment. [Figure 8] FIG. 8A is a diagram showing the distribution of computational particles at the start of calculation, and FIG. 8B is a diagram showing the distribution of computational particles in a steady state. [Figure 9] FIG. 9 is a graph showing the calculation results of the distribution of the time-averaged values ​​of the pressure and the mean radius of the droplets on the central axis of the flow channel in the steady state, compared with the experimental results presented in the paper. [Figure 10] 10A and 10B are diagrams showing a cross section perpendicular to the flow direction of a flow channel that is to be analyzed by a simulation method according to a modified example. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0012] [Overview of particle method] Before describing an embodiment of the present invention, a general SPH (Smoothed Particle Hydrodynamics) method will be briefly described. The SPH method is a Lagrangian type numerical calculation method that discretizes a continuum using a set of multiple calculation particles and calculates the time evolution of a flow field. The SPH method uses a continuous and differentiable kernel function as the mass distribution of a single particle, and also uses a discretization model in which the superposition of the differential coefficients of the kernel function is the spatial derivative.

[0013] In the SPH method, the mass m of a computational particle is distributed around it by a kernel function W, which is continuously distributed, and the density field is given as a superposition of the mass m and the density field of the i-th computational particle. i Density ρ i is expressed by the following formula:

number

[0014] For example, the following function can be used as the kernel function W. The following function is called the C2Wendland kernel.

number

[0015] h is the smoothing length, also called the kernel width. The kernel width h is equal to the diameter of an equivalent sphere obtained by replacing each of the multiple calculation particles with a sphere having the same mass and physical density. In this specification, the kernel width h may be referred to as the diameter of the calculation particle.

[0016] Any position x i The value of any scalar quantity A and vector quantity A in can be expressed by the following equation by superposition of kernel functions.

number

[0017] The gradient of any scalar quantity A is expressed by the following equation by superposing the derivatives of the kernel function:

number

number

[0018] Here, ∇W(x ij ,h) is the derivative of the kernel function, and when the C2Wendland kernel is used as the kernel function, it is expressed by the following equation.

number

[0019] The divergence for any vector quantity A can be expressed by the following equation, similar to equation (4).

number

[0020] The second derivative with respect to any scalar quantity A can be calculated, for example, using the following formula:

number

[0021] [Governing equations for wet steam flow] The governing equations for compressible flow with phase change, assuming inviscid, adiabatic, homogeneous flow, are expressed as follows:

number

number

number

number

number

[0022] When the above-mentioned equations (9) to (13) are discretized using the SPH method, they can be expressed as follows:

number

number

number

number

number

number

number

number

[0023] [Calculation of characteristic dimensions of the flow path] In this embodiment, calculation particles are divided or combined during calculation based on the characteristic dimensions of the flow path through which the fluid flows. Next, the characteristic dimensions of the flow path will be described with reference to Figs. 1A to 1D.

[0024] FIG. 1A is a cross-sectional view showing an example of the shape of a flow channel 10 through which a fluid to be analyzed flows. The shape and size of the flow channel 10 are determined by an object Ω. The shape and size of the object Ω are given by, for example, CAD data. The surface of the flow channel 10 defined by the object Ω is denoted as ∂Ω. In FIG. 1A, the object Ω is placed in a two-dimensional analysis space, and a two-dimensional flow channel 10 is defined.

[0025] In FIG. 1A, a fluid flows in from the left end of the flow channel 10 and flows out from the right end. The flow channel width at the inflow end of the flow channel 10 is labeled Wi, the flow channel width at the outflow end is labeled Wo, and the length of the flow channel 10 is labeled L. The flow channel has a constant flow channel width portion near the inflow end, and the flow channel width gradually narrows toward the downstream to form a throat portion, and then gradually widens. The length of the portion where the flow channel width gradually narrows is shorter than the length of the portion where the flow channel width gradually widens. The flow channel width is constant near the outflow end. The flow channel width Wo at the outflow end is smaller than the flow channel width Wi at the inflow end. The flow channel width at the throat portion is labeled Wt.

[0026] In fluid analysis using the particle method, when analyzing the flow around an object, it is necessary to calculate the distance from each of the positions of multiple calculation particles to the object surface in order to provide boundary conditions on the object surface.

[0027] To calculate the distance from an object surface to an arbitrary position in the analytical space, for example, a signed distance function can be used. If an arbitrary position in the analytical space is denoted as x and the value of the signed distance function at position x is denoted as f(x), f(x) is defined as follows:

number

[0028] Figure 1B shows the absolute value of the signed distance function in shading. The value of the signed distance function is zero at the surface ∂Ω of the object Ω, and the absolute value of the signed distance function increases (becomes darker in Figure 1B) as one moves away from the surface ∂Ω. Inside the flow path 10 (outside the object Ω), the value of the signed distance function is positive, and outside the flow path 10 (inside the object Ω), the value of the signed distance function is negative. The dashed lines in Figure 1B show contour lines of the values ​​of the signed distance function.

[0029] The thin solid lines in Figure 1B represent an orthogonal uniform grid stretched in the analysis space. The value of the signed distance function is calculated at each of the multiple grid points of the orthogonal uniform grid. Note that the grid spacing of the orthogonal uniform grid used in the actual analysis is narrower than that of the orthogonal uniform grid shown in Figure 1B.

[0030] FIG. 1C is a schematic diagram for explaining a method for determining a characteristic dimension at an arbitrary position of the flow channel 10. The characteristic dimension is calculated at each position of the lattice points of an orthogonal equally spaced lattice laid out in the analytical space. First, a gradient vector of the value of the signed distance function is calculated at each position of the lattice points. In FIG. 1C, the gradient vector at each position of some lattice points is represented by an arrow. The gradient vector at an arbitrary position in the analytical space is perpendicular to the isopleth of the value of the signed distance function and points in the direction of the larger value of the signed distance function.

[0031] Next, a local maximum value of the signed distance function is searched for along the direction of the gradient vector. For example, the step size of the multiple calculation positions on the line along the gradient vector when searching for the local maximum value is set to be equal to the lattice interval of the orthogonal equal-spaced lattice. The starting point of the arrow indicating the gradient vector shown in FIG. 1C corresponds to the position where the characteristic dimension is searched for, and the end point corresponds to the position where the value of the signed distance function is the local maximum. A value twice the searched local maximum value is adopted as the characteristic dimension at each position of the lattice point.

[0032] For grid points with negative values ​​of the signed distance function (grid points inside the object Ω), the above calculation is not performed and zero is assigned as the characteristic dimension.

[0033] FIG. 1D is a diagram showing the distribution of the calculated characteristic dimensions in shading. In FIG. 1D, the characteristic dimensions at the lattice points in region A1 are relatively large, and the characteristic dimensions at the lattice points in region A3 are relatively small. The characteristic dimensions at the lattice points in regions A2 and A4 are intermediate between the sizes in region A1 and region A4. That is, in the two-dimensional flow channel 10, the characteristic dimensions at the lattice points in the region where the flow channel width is relatively wide are relatively large, and the characteristic dimensions at the lattice points in the region where the flow channel width is relatively narrow are relatively small. That is, it can be said that the characteristic dimensions reflect the value of the flow channel width.

[0034] [Merging and splitting of computational particles] Next, the conditions for merging and splitting computational particles will be explained. A characteristic dimension of the flow path at each position of the multiple computational particles is calculated. For example, the characteristic dimension of the flow path at each position of the multiple computational particles can be calculated by performing a linear interpolation calculation using the characteristic dimension values ​​of the flow path assigned to each lattice point of an orthogonal equally spaced lattice stretched in the analysis space. In this case, if the value of the signed distance function assigned to a lattice point is negative, the lattice point is excluded from the target of the interpolation calculation.

[0035] The serial number given to the calculation particle of interest is denoted as i, and the characteristic dimension of the flow channel at the position of the i-th calculation particle is denoted as w i and the kernel width (diameter) of the i-th calculation particle is denoted as h i It is marked as such.

[0036] If the following equation (23) is satisfied, it is determined that the spatial resolution is insufficient, and the calculation particles are divided.

number

[0037] If the following equation (24) is satisfied, it is determined that the spatial resolution is excessive, and the calculation particles are merged.

number

[0038] Fig. 2 is a schematic diagram showing an example of the merging and division of computational particles during the analysis of a flow field. As described with reference to Fig. 1D, the characteristic dimension is relatively large in the region A1, relatively small in the region A3, and intermediate between the characteristic dimensions in the regions A2 and A4.

[0039] A computational particle 20 moves from the inlet end (the left end in FIG. 2) of the flow channel 10 toward the outlet end (the right end in FIG. 2). When the computational particle 20 enters the region A2 from the region A1, the characteristic dimension value w i When equation (23) is satisfied, the computational particle 20 is divided into a plurality of computational particles 20, and the diameter of each of the computational particles 20 is reduced. Similarly, when the computational particle 20 enters the region A3 from the region A2, the characteristic dimension value w iWhen equation (23) is satisfied, the computational particle 20 is further divided, and the diameter of the computational particle 20 becomes smaller.

[0040] When the computational particle 20 enters the region A4 from the region A3, the characteristic dimension value w i As a result, when equation (24) is satisfied, multiple computational particles 20 are combined into one computational particle 20, and the diameter of the computational particle 20 becomes larger.

[0041] Next, a method for dividing computational particles will be described with reference to FIGS. 3A to 4B. Fig. 3A is a schematic diagram of a computational particle 20a before division. Fig. 3B is a schematic diagram of two computational particles 20a1 and 20a2 generated by dividing one computational particle 20a (Fig. 3A). Fig. 4A is a schematic diagram of four computational particles 20a1-20a4 generated by dividing one computational particle 20a (Fig. 3A). Fig. 4B is a schematic diagram of eight computational particles 20a1-20a8 generated by dividing one computational particle 20a (Fig. 3A).

[0042] First, as shown in FIG. 3B, a case where the particle is divided into two computational particles 20a1 and 20a2 will be described. The position of the computational particle 20a before division is denoted as x1, the kernel width as h1, the mass as m1, and the volume as V1. The kernel width of each of the two computational particles 20a1 and 20a2 after division is denoted as h2, the mass as m2, and the volume as V2. The positions of the two computational particles 20a1 and 20a2 after division are denoted as x 21 , x 22 Here, V1=h1 3 , V2=h2 3 It is.

[0043] The mass m2, volume V2, and kernel width h2 of each of the computational particles 20a1 and 20a2 after division can be calculated by the following equations.

number

[0044] The positions x of the two computational particles 20a1 and 20a2 after division21 , x 22 can be given, for example, by the following formula:

number

[0045] Next, as shown in FIG. 4A, a case where the particle is divided into four computational particles 20a1 to 20a4 will be described. The kernel width of each of the four computational particles 20a1 to 20a4 after division is denoted as h2, the mass as m2, and the volume as V2. The positions of the four computational particles 20a1 to 20a4 after division are respectively denoted by x 21 ~x 24 It is marked as such.

[0046] The mass m2, volume V2, and kernel width h2 of each of the computational particles 20a1 to 20a4 after division can be calculated by the following formulas.

number

[0047] The positions x of the four computational particles 20a1 to 20a4 after division 21 ~x 24 can be given, for example, by the following formula:

number

[0048] Next, as shown in FIG. 4B, a case where the particle is divided into eight computational particles 20a1 to 20a8 will be described. The kernel width of each of the eight computational particles 20a1 to 20a8 after division is denoted as h2, the mass as m2, and the volume as V2. The positions of the eight computational particles 20a1 to 20a8 after division are respectively denoted by x 21 ~x 28 It is marked as such.

[0049] The mass m2, volume V2, and kernel width h2 of each of the computational particles 20a1 to 20a8 after division can be calculated by the following formulas.

number

[0050] The positions x of the eight computational particles 20a1 to 20a8 after division 21 ~x 28 can be given, for example, by the following formula:

number

[0051] Any physical quantity A of the sth computational particle 20s after division 2s can be given by the following equation using the gradient of the physical quantity A1 of the computational particle 20 before division, for example.

number

[0052] Even when the analysis space is two-dimensional, the positions of the computational particles 20 after division and the physical quantities of the computational particles 20 can be obtained in a similar manner.

[0053] Next, a method for merging computational particles will be described with reference to Figures 5A and 5B. When merging computational particles, it is determined whether or not there is a partner particle that constitutes a nearest pair. A "nearest pair" refers to a pair of computational particles that have a nearest relationship with each other.

[0054] The nearest pair will be specifically described with reference to Fig. 5A. Fig. 5A is a schematic diagram showing an example of the distribution of a plurality of computational particles 20. The arrows shown in Fig. 5A point from each of the computational particles 20 to the nearest computational particle 20. For example, the nearest particle of the computational particle 20a is the computational particle 20b, and the nearest particle of the computational particle 20b is the computational particle 20a. Therefore, the two computational particles 20a and 20b form a nearest pair.

[0055] The nearest particle of the computational particle 20c is the computational particle 20d. However, the nearest particle of the computational particle 20d is the computational particle 20e, not the computational particle 20c. Therefore, the computational particle 20c does not form a nearest pair with another computational particle 20. The nearest particle of the computational particle 20e is the computational particle 20d. Therefore, the computational particles 20d and 20e form a nearest pair.

[0056] When there is a partner computational particle 20 forming a nearest pair, it is determined whether or not the spatial resolution at the position of the partner particle forming the nearest pair is excessive. If it is determined that the spatial resolution of the partner particle is excessive, the two computational particles 20 forming the nearest pair are merged to generate one large computational particle having a kernel width larger than the kernel width of either of the two computational particles 20. In other words, when the spatial resolutions of the two computational particles 20 forming the nearest pair are both excessive, the two computational particles 20 are merged.

[0057] 5A, the computational particles 20 determined to have excessive spatial resolution are hatched. That is, the computational particles 20a, 20b, 20c, and 20d are determined to have excessive spatial resolution, while the other computational particles 20 are not determined to have excessive spatial resolution.

[0058] The computational particles 20a and 20b constituting the nearest pair are both determined to have excessive spatial resolution, so the computational particles 20a and 20b are merged. The computational particles 20d and 20e constituting the nearest pair are not merged because the spatial resolution of one of the computational particles 20e is not excessive.

[0059] 5B is a schematic diagram showing an example of the distribution of the computational particles 20 after the computational particles 20a and 20b are merged. The computational particles 20a and 20b are merged to generate a large computational particle 20ab. In FIG. 5B, the computational particle 20ab after the merger is indicated by relatively light hatching.

[0060] Next, a method for combining two computational particles 20 will be described. The position of one of the computational particles 20a before merging is x 11 , the kernel width is h 11 , mass m 11 , volume is V 11 and the position of the other computational particle 20b is denoted as x 12 , the kernel width is h 12 , mass m 12 , volume is V 12 The mass m2 and kernel width h2 of the computational particle 20ab (FIG. 5B) after merging can be calculated using the following equations.

number

[0061] The position x2 of the computational particle 20ab after merging is the position x of the two computational particles 20a and 20b before merging. 11 , x 12 can be determined using a mass-weighted average as follows:

number

[0062] Similarly, the arbitrary physical quantity A2 of the computational particle 20ab after the merging is the physical quantity A 11 , A 12 It can be calculated using the following formula:

number

number

[0063] 6 is a block diagram of a simulation device according to the embodiment. The simulation device according to the embodiment includes an input unit 30, a processing unit 31, an output unit 32, and a storage unit 33. Simulation conditions and the like are input from the input unit 30 to the processing unit 31. Furthermore, various commands and the like are input from a user to the input unit 30. The input unit 30 is composed of, for example, a communication device, a removable media reader, a keyboard, a pointing device, and the like.

[0064] The processing unit 31 executes a simulation using the particle method based on the input simulation conditions and commands. Furthermore, it outputs the simulation results to the output unit 32. The simulation results include information representing the state of particles in the particle system that is the simulation target, and the change over time in the physical quantity of the particle system. The processing unit 31 includes, for example, a central processing unit (CPU) of a computer. A program for causing the computer to execute the simulation using the particle method is stored in the storage unit 33. The output unit 32 includes a communication device, a removable media writing device, a display, and the like.

[0065] 7 is a flowchart showing the procedure of the simulation method according to this embodiment. Each step of the flowchart shown in FIG.

[0066] First, the processing unit 31 acquires the simulation conditions inputted by the operator to the input unit 30 (step S1). The simulation conditions include information defining the fluid to be simulated, initial conditions, boundary conditions, etc. The information defining the fluid includes physical property values ​​such as the density and viscosity of the fluid. The initial conditions include information for arranging a plurality of computational particles in the analysis space, a method for specifying the initial physical property values ​​of the computational particles, etc. The boundary conditions include information specifying the shape and size of the analysis space, information specifying the shape of the object that defines the flow path, etc.

[0067] Next, the processing unit 31 calculates the value of the signed distance function at each position of the lattice points of the orthogonal equally spaced lattice for the given object shape (step S2). Furthermore, the processing unit 31 calculates the characteristic dimension of the flow channel 10 at each position of the lattice points using the value of the signed distance function (step S3). For this calculation, the method described with reference to FIG. 1C can be used. By this calculation, the distribution of the characteristic dimension shown in FIG. 1D is obtained.

[0068] Based on the input simulation conditions and the value of the signed distance function, a plurality of computational particles are placed in the analysis space (step S4). Specifically, a plurality of computational particles are placed in the region (inside the flow channel 10) where the value of the signed distance function is positive.

[0069] Next, the governing equations (14) to (19) are numerically solved for each of the calculation particles to calculate the physical quantities assigned to each of the calculation particles, and the physical quantities and positions of each of the calculation particles are updated (step S5). After that, for each of the calculation particles, the spatial resolution (diameter of the calculation particle) relative to the characteristic dimension of the flow path at the position of the calculation particle is evaluated (step S6). Specifically, it is determined whether or not the formulas (23) and (24) hold. If the formula (23) holds, it is determined that the spatial resolution is insufficient, and if the formula (24) holds, it is determined that the spatial resolution is excessive.

[0070] If the spatial resolution is insufficient for the characteristic dimensions of the flow channel, the calculation particles are divided (step S7). The calculation particles are divided, for example, by the method described with reference to FIGS. 3A to 4B. The number of divisions is determined by the kernel width h i and flow path characteristics and dimensions w i It is preferable to determine the spatial resolution based on the ratio of the spatial resolution to the characteristic dimensions of the flow channel. Furthermore, if the spatial resolution is excessive with respect to the characteristic dimensions of the flow channel, the calculation particles are combined (step S8). The calculation particles are combined, for example, by the method described with reference to Figs. 5A to 5B. Note that step S8 may be performed before step S7.

[0071] The procedure from step S5 to step S8 is repeated until the calculation end condition is satisfied (step S9). If the calculation end condition is satisfied, the calculation result is output to the output unit 32 (step S10).

[0072] Next, the results of simulating the flow of fluid in a flow channel using the simulation method according to this embodiment will be described with reference to Figures 8A to 9. The flow channel to be analyzed was the same two-dimensional flow channel as that shown in Figure 1A. The flow channel width Wi at the inflow end was 40 mm, the flow channel width Wo at the outflow end was 18 mm, the flow channel width Wt at the throat was 10 mm, and the flow channel length L was 200 mm.

[0073] The cross-sectional shape of the channel is the same as that used in the experiment described in CA Moses, GD Stein, “On the Growth of Steam Droplets Formed in a Laval Nozzle Using Both Static Pressure and Light Scattering Measurement”, Journal of Fluids Engineering, 100(3), 311-322, September 1987, ISSN 0098-2202, doi: 10.1115 / 1.3448672, https: / / doi.org / 10.1115 / 1,3448672. The flow conditions were the stagnation pressure p0 of 67.66 kPa and the stagnation temperature T0 of 376.7 K, following the experiment described in the above paper.

[0074] At the initial condition of the simulation, the kernel width h=4 / (2 1 / 2 The calculation particles were placed with a diameter of 1 mm and spaced 4 mm apart. The pressure and temperature of each calculation particle were set to the stagnation pressure and stagnation temperature, respectively, and the initial velocity was set to zero.

[0075] FIG. 8A is a diagram showing the distribution of the calculation particles 20 at the start of the calculation. -6 The time period was set to 100 s, and the positions and velocities of the computational particles were allowed to evolve over time from the initial state until the flow field reached a steady state.

[0076] Fig. 8B is a diagram showing the distribution of the computational particles 20 in the steady state. In Fig. 8B, the shading of the computational particles 20 represents the kernel width. In the region A1 where the characteristic dimension of the flow path is relatively large, the kernel width is relatively large, and in the region A3 where the characteristic dimension of the flow path is relatively small, the kernel width is relatively small. The kernel width in the regions A2 and A4 is intermediate between the kernel width in the region A1 and the kernel width in the region A3.

[0077] Figure 9 is a graph showing the calculation results of the distribution of the time-averaged values ​​of pressure and the mean radius of droplets on the central axis of the channel in the steady state, compared with the experimental results shown in the above paper. The horizontal axis represents the distance from the throat of the channel in units of [cm], the left vertical axis represents the normalized pressure normalized by the stagnation pressure p0, and the right vertical axis represents the mean radius of droplets in units of [nm]. The mean diameter r of the droplets was calculated using the following formula.

number

[0078] In Fig. 9, the solid line indicates the calculation results by the simulation method according to this embodiment, and the circle symbols indicate the experimental results described in the above paper. It can be seen that the results obtained by the calculation are almost consistent with the experimental results in terms of both the pressure distribution and the distribution of the average radius of the droplets. From the calculation results shown in Fig. 9, it was confirmed that it is possible to analyze the wet steam flow by the simulation method according to this embodiment.

[0079] Next, the excellent effects of this embodiment will be described. In this embodiment, in steps S7 and S8 (FIG. 7), a necessary and sufficient spatial resolution is automatically determined according to the characteristic dimensions of the flow path. Since the operator does not need to define an area requiring high spatial resolution and an area that may be analyzed with low spatial resolution, the operator is not required to have high skill in determining an appropriate spatial resolution according to the location in the flow path.

[0080] Furthermore, there is no need to perform preprocessing such as dividing the analysis space into multiple regions according to the required spatial resolution before performing the analysis, which reduces the number of steps required for the analysis.

[0081] Next, a modified example of the above embodiment will be described with reference to Figs. 10A and 10B. In the above embodiment, the flow field in a two-dimensional flow channel is analyzed, but in the modified example described with reference to Figs. 10A and 10B, the flow field in a three-dimensional flow channel is analyzed. Figs. 10A and 10B are diagrams showing a cross section (hereinafter simply referred to as "flow channel cross section") perpendicular to the flow direction of the flow channel to be analyzed by the simulation method according to the modified example. In the case of a three-dimensional flow channel, the value of the signed distance function can be calculated as described with reference to Fig. 1B, and the characteristic dimensions of the flow channel can be calculated as described with reference to Fig. 1C.

[0082] In the modified example shown in FIG. 10A, the flow channel 10 is rotationally symmetric about the central axis C, and the cross section of the flow channel at any position in the flow direction is circular. The diameter of the cross section of the flow channel at a certain position in the flow direction is denoted as D. The contour lines of the signed distance function defined for a flow channel with a circular cross section are concentric circles. The value of the signed distance function has a maximum value at the central axis C in the radial direction. Therefore, when the maximum value of the signed distance function is searched for from an arbitrary point P in the direction of the gradient vector, the value of the signed distance function at the position of the central axis C of the flow channel cross section is detected as the maximum value. Twice this maximum value is adopted as the characteristic dimension of the flow channel at the position of point P.

[0083] In the modified example shown in Fig. 10B, the cross section of the flow channel 10 at any position in the flow direction is an equilateral triangle, and when the cross section of the flow channel is moved in the flow direction, the center of gravity of the cross section of the flow channel moves along a straight line. The straight line followed by the center of gravity of the cross section of the flow channel is referred to as the central axis C of the flow channel. When a maximum value of the signed distance function is searched for from any point P in the direction of the gradient vector, a maximum value appears at a point Pm where the function intersects with a plane formed by the bisectors of each vertex of the cross section of the flow channel at any position in the flow direction (hereinafter referred to as "bisect plane B"). Twice this maximum value is adopted as the characteristic dimension of the flow channel at point P.

[0084] In this way, it is possible to define a characteristic dimension of a three-dimensional channel at any point in the channel, for example, the characteristic dimension of the channel can be said to be a quantity that depends on the degree of expansion of the channel at the position of the calculation particle.

[0085] For example, in the flow channel having an equilateral triangular cross section shown in Fig. 10B, near the midpoint of one side, the degree of expansion of the flow channel is small in the direction toward the nearest flow channel wall, but is large in the direction away from the flow channel wall and in two directions parallel to the flow channel wall. Near the apex of the flow channel cross section, the degree of expansion of the flow channel is small in three directions, the direction toward the nearest flow channel wall, the direction parallel to the flow channel wall toward the apex, and the direction away from the flow channel wall, and is large only in one direction, the direction parallel to the flow channel wall toward the opposite side from the apex. It can be said that the degree of expansion of the flow channel near the apex is smaller than the degree of expansion of the flow channel near the midpoint of each side.

[0086] The characteristic dimension of the flow channel near the apex, where the degree of spatial spread of the flow channel is relatively small, is smaller than the characteristic dimension of the flow channel near the midpoints of each side, where the degree of spatial spread of the flow channel is relatively large. In this way, the characteristic dimension of the flow channel at any position can be said to be an index representing the degree of spatial spread of the flow channel at that position.

[0087] Each embodiment is merely an example, and it goes without saying that partial replacement or combination of the configurations shown in different embodiments is possible. Similar effects due to similar configurations of multiple embodiments are not mentioned in each embodiment. Furthermore, the present invention is not limited to the above-described embodiments. For example, it will be obvious to those skilled in the art that various modifications, improvements, combinations, etc. are possible. [Explanation of symbols]

[0088] 10 Flow Path 20, 20a, 20a1, 20a2, 20ab, 20b, 20c, 20d, 20e, 20s calculation particles 30 Input section 31 Processing section 32 Output section 33 Storage section

Claims

1. A simulation device that analyzes a fluid flow inside a flow path using a particle method, comprising: an input section for inputting simulation conditions; a processing unit that performs calculations based on the simulation conditions inputted to the input unit; Equipped with The processing unit includes: a function of expressing a fluid by a plurality of computational particles, assigning a physical quantity to each of the plurality of computational particles, arranging the plurality of computational particles in an analysis space in which the flow path is modeled, and evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time by solving a governing equation for the fluid; a function of performing a process of comparing a characteristic dimension of the flow path at the position of each of the plurality of computational particles with a size of each of the plurality of computational particles in a process of evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time, and performing a process of dividing each of the plurality of computational particles to generate smaller computational particles or a process of combining some of the plurality of computational particles to generate larger computational particles according to a comparison result; A simulation device having the above configuration.

2. The simulation device according to claim 1 , wherein a characteristic dimension of the flow path at each position of the plurality of computational particles is determined based on a degree of spatial spread of the flow path at each position of the plurality of computational particles.

3. a value of a shortest distance to a wall surface of the flow path is assigned to each of a plurality of lattice points of a lattice laid out in the analysis space; 3. The simulation device according to claim 1, wherein the processing unit further has a function of determining the characteristic dimension at each position of the plurality of computational particles based on the value of the shortest distance assigned to a lattice point located in the vicinity of each position of the plurality of computational particles.

4. The processing unit includes:

4. The simulation device according to claim 3, further comprising: a search for a maximum value of the shortest distance value in a direction representing a gradient of the shortest distance value assigned to the plurality of lattice points from each position of the plurality of computational particles; and determining the characteristic dimension based on the maximum value obtained by the search.

5. A simulation method for analyzing a fluid flow inside a flow path using a particle method, comprising the steps of: A computer given the simulation conditions expressing a fluid using a plurality of computational particles based on given simulation conditions, assigning a physical quantity to each of the plurality of computational particles, arranging the plurality of computational particles in an analysis space in which the flow path is modeled, and evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time by solving a governing equation for the fluid; A simulation method in which, in a process of evolving the physical quantities assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time, a characteristic dimension of the flow path at the positions of each of the plurality of computational particles is compared with the size of each of the plurality of computational particles, and, depending on the comparison result, each of the plurality of computational particles is divided to generate smaller computational particles, or some of the plurality of computational particles are combined to generate larger computational particles.

6. A program for causing a computer to realize a function of analyzing a fluid flow inside a flow path using a particle method, a function of expressing a fluid as a plurality of computational particles based on given simulation conditions, assigning a physical quantity to each of the plurality of computational particles, arranging the plurality of computational particles in an analysis space in which the flow path is modeled, and evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time by solving a governing equation for the fluid; a function of performing a process of comparing a characteristic dimension of the flow path at the position of each of the plurality of computational particles with a size of each of the plurality of computational particles in a process of evolving the physical quantity assigned to each of the plurality of computational particles and the positions of the plurality of computational particles over time, and performing a process of dividing each of the plurality of computational particles to generate smaller computational particles or a process of combining some of the plurality of computational particles to generate larger computational particles according to a comparison result; A program that enables a computer to achieve the above.

Images