Simulation device, simulation method, and program

Abstract

To provide a simulation device that keeps analysis accuracy high, and can reduce an amount of calculation in an analysis using a particle method.SOLUTION: A simulation condition is input to an input unit. A processing unit is configured to, on the basis of the simulation condition input to the input unit, subject a physical quantity and a position of each of a plurality of particles to time evolution, using a particle method. The processing unit is configured to divide an analysis space into a plurality of divided areas, and dispose a plurality of particles different in a size in each of the plurality of divided areas. Upon subjecting the positions of the plurality of particles to the time evolution as for each of the plurality of divided areas, the processing unit is configured to determine whether or not to add a new particle to the divided area so that a sum of volume of the particle within the divided area is even as for each of the plurality of divided areas; and add the new particle to the divided area having the addition of the new particle determined.SELECTED DRAWING: Figure 6

Description

[Technical field]

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

[0002] There is a known simulation method that analyzes the behavior of a fluid by approximating the movement of the fluid as the movement of a particle system. This simulation method is called a particle method. In the particle method, the fluid is represented by multiple particles, and the wall boundary placed in the fluid flow field is also generally represented by multiple particles. [Prior art documents] [Non-patent literature]

[0003] [Non-Patent Document 1] Takahiro Harada, Seiichi Koshizuka, "Improvement of the wall boundary calculation method in SPH", Transactions of Information Processing Society of Japan, Vol. 48, No. 4, pp. 1838-1846, 2007. Summary of the Invention [Problem to be solved by the invention]

[0004] In fluid analysis using the particle method, the spatial resolution is determined by the diameter of the particles. The particle diameter must be set to be sufficiently small compared to the length scale of the fluid motion to be analyzed. In general, the length scale of fluid motion varies depending on the location. For example, in areas far from the wall of an object placed in a fluid, a sufficiently accurate analysis can be performed even if the particle diameter is large. However, near a wall, the length scale of fluid motion is represented by the thickness of the boundary layer, so in order to accurately analyze the motion of the fluid near the wall, the particle diameter must be smaller than the thickness of the boundary layer.

[0005] In order to perform a highly accurate analysis in the entire analysis space, the diameter of the particles must be made small according to the thickness of the boundary layer, even in areas far from the wall, which increases the number of particles placed in the analysis space, resulting in an increase in the amount of calculations.

[0006] An object of the present invention is to provide a simulation device, a simulation method, and a program that are capable of maintaining high analytical accuracy and reducing the amount of calculation in an analysis using a particle method. [Means for solving the problem]

[0007] According to one aspect of the present invention, A simulation device for analyzing a fluid flow using a particle method, comprising: an input section for inputting simulation conditions; a processing unit that performs time evolution of each position of a plurality of particles using a particle method based on the simulation conditions inputted to the input unit; Equipped with The processing unit includes: Divide the analysis space into multiple sub-domains, A plurality of particles having different sizes are arranged in each of the plurality of divided regions; determining whether or not to add a new particle to each of the plurality of divided regions so that a sum of the volumes of the particles at the current time in the divided region is constant when the positions of the plurality of particles in each of the plurality of divided regions are evolved over time; There is provided a simulation device that adds new particles to a divided region where it has been determined that new particles should be added.

[0008] According to another aspect of the invention, A simulation method for analyzing a fluid flow using a particle method, comprising the steps of: Divide the analysis space into multiple sub-domains, A plurality of particles having different sizes are arranged in each of the plurality of divided regions; determining whether or not to add a new particle to each of the plurality of divided regions so that a sum of the volumes of the particles at the current time in the divided region is constant when the positions of the plurality of particles in each of the plurality of divided regions are evolved over time; A simulation method is provided for adding new particles to a divided region where it has been determined that new particles should be added.

[0009] According to yet another aspect of the present invention, A step of dividing an analysis space into a plurality of divided domains; A step of arranging a plurality of particles having different sizes for each of the plurality of divided regions; a step of determining whether or not to add a new particle to each of the plurality of divided regions so that a sum of the volumes of the particles at the current time in the divided region becomes constant when the positions of the plurality of particles in each of the plurality of divided regions are evolved over time; a step of adding new particles to the divided area where it is determined that new particles should be added; A program for causing a computer to execute the above is provided. Effect of the Invention

[0010] By placing multiple particles of different sizes in each divided region, it is possible to make the spatial resolution different for each divided region. By making the particles smaller in divided regions where high spatial resolution is required, it is possible to perform analysis at the required high spatial resolution. By making the particles larger in divided regions where high spatial resolution is not required, the number of particles can be reduced. This makes it possible to reduce the amount of calculations in the simulation. By adding new particles so that the sum of the current volumes of the particles in the divided regions remains constant, it is possible to simulate a compressible flow field in which the number of particles in the divided regions fluctuates. [Brief description of the drawings]

[0011] [Figure 1] FIG. 1 is a perspective view showing an example of an analytical model to be analyzed by a simulation apparatus according to the present embodiment. [Diagram 2] FIG. 2 is a schematic diagram showing a part of the analytical space in two dimensions. [Diagram 3] FIG. 3A is a schematic diagram showing a portion of two adjacent divided regions, and FIG. 3B is a schematic diagram showing how physical quantities are exchanged between the two divided regions. [Figure 4] FIG. 4 is a schematic diagram showing a part of the analysis space in two dimensions. [Diagram 5] FIG. 5 is a block diagram of the simulation device according to this embodiment. [Figure 6] FIG. 6 is a flowchart showing the procedure of the simulation method of this embodiment. [Figure 7] 7A and 7B are diagrams showing the initial state of particles placed in the analysis space. [Figure 8] 8A and 8B show the instantaneous particle distribution and Mach number distribution when the flow field reaches a steady state. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

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

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

number

number

[0014] where x ij is the position vector x i and the position vector x j The absolute value of the difference between x and q = x ij / h, where 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 multiple particles with spheres of the same mass and physical density. In this specification, the kernel width h may be referred to as the diameter of a particle.

[0015] 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

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

number

number

[0017] 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

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

number

[0019] Next, we will explain the governing equations when analyzing compressible flows. An example of discretizing the governing equations for inviscid compressible flows using the SPH method is shown below.

number

number

number

[0020] Next, a simulation device and a simulation method according to an embodiment of the present invention will be described with reference to FIGS. 1 to 8B.

[0021] In fluid analysis using particle methods such as the SPH method, the spatial resolution of the analysis is determined by the diameter of the computational particles. The diameter of the computational particles must be set sufficiently small relative to the length scale of the fluid motion to be analyzed. In general, the length scale of the fluid motion varies depending on the location. For example, near the wall of an object, the length scale of the fluid motion is represented by the boundary layer thickness, and high spatial resolution is required. On the other hand, in areas away from the wall, the length scale of the fluid motion is often larger than the boundary layer thickness. Therefore, in areas away from the wall, analysis can be performed with a lower spatial resolution than in areas near the wall. Since numerical analysis must be performed with limited computational resources, it is necessary to set the computational particle diameter small in areas requiring high spatial resolution to ensure the necessary spatial resolution, while at the same time, to set the computational particle diameter large in areas not requiring high spatial resolution to save the number of computational particles required for analysis.

[0022] Conventionally, a method has been proposed that divides the analysis space into multiple rectangular parallelepiped regions and enables the spatial resolution to be set for each divided region. However, when analyzing a flow field with complex wall shapes, it is difficult to effectively divide the region into those that require high spatial resolution and those that do not, even if the region is divided into a simple shape such as a rectangular parallelepiped. In addition, the conventional method targets incompressible flows and cannot be used for highly compressible flows that involve sudden expansion or contraction of gas or the generation of shock waves.

[0023] FIG. 1 is a perspective view showing an example of an analytical model to be analyzed by the simulation device according to the present embodiment. The analytical model includes an analytical space 10 and a wall 11 arranged in the analytical space 10. The analytical space 10 is divided into a plurality of regions (hereinafter, referred to as divided regions Ω). In the example shown in FIG. 1, the analytical space 10 is divided into three divided regions Ω. 1 , Ω 2 , Ω 3 The innermost first divided region Ω 1 contains the wall 11, and the second divided region Ω 2 The second divided domain Ω of the analysis space 10 contains the first divided domain. 2 The area outside of is the third divided area Ω3 is equivalent to.

[0024] A plurality of particles 20 are arranged in each divided region Ω. Note that in FIG. 1, only one particle 20 is shown for each divided region Ω. The diameter of the particle 20 varies for each divided region Ω. For example, in the first divided region Ω, which is the innermost region, 1 The particle 20 in the innermost region is the smallest, and the third region Ω is the outermost region. 3 The particle 20 in the divided region Ω is the largest. The size of the particles 20 in one divided region Ω is the same. The size of the particles 20 is determined according to the spatial resolution required to analyze the movement of the fluid in each divided region Ω.

[0025] When applying the SPH method, the following three technical elements are required to divide the analysis space 10 into a plurality of divided regions Ω of arbitrary shapes and to change the spatial resolution for each divided region Ω. That is, first, a technology to divide the analysis space 10 into a plurality of divided regions Ω of arbitrary shapes, second, a technology to exchange physical quantities between adjacent regions of the divided region Ω, and third, a technology to control particles flowing into and out of each divided region Ω. These technologies are described below.

[0026] [Technology for dividing the analytical space into multiple divided domains] First, a technique for dividing an analysis space into a plurality of divided regions Ω will be described. For example, a signed distance function (SDF) is used as a method for dividing the analytical space 10 (Fig. 1) into multiple divided regions Ω of arbitrary shapes. The signed distance function φ is a continuous function defined in the entire analytical space 10, and has the shortest distance from a given arbitrary boundary S to an arbitrary position x as the function value. Moreover, the sign of the signed distance function φ is negative inside the boundary S and positive outside it.

[0027] To calculate the value of the signed distance function φ, we set φ=0 on the boundary S as the initial condition and solve the advection equation below.

number

[0028] The position and shape of a boundary S for dividing the analytical space 10 are determined, and a signed distance function φ is calculated using the boundary S. The analytical space 10 is divided into a plurality of divided regions Ω by at least one iso-value surface for which the value of the signed distance function φ is equal. The iso-value surface becomes the boundary surface between two mutually contacting divided regions Ω. The shape of the boundary S may be, for example, based on the shape of a wall 11 to be placed in the analytical space 10.

[0029] The shape of the boundary S can be arbitrarily determined by the user using CAD data, etc. Therefore, the analytical space 10 can be divided into a plurality of divided regions Ω of arbitrary shapes.

[0030] Next, with reference to FIG. 2, a specific example of finding the value of the signed distance function φ at a plurality of discretized lattice points in the analysis space 10 will be described.

[0031] FIG. 2 is a schematic diagram showing a partial region of the analytical space 10 in two dimensions. A boundary S is set within the analytical space 10. Furthermore, an orthogonal equal-distance grid is laid out within the analytical space 10. For each of multiple lattice points P in the vicinity of the boundary S, the exact value of the distance Lps from the lattice point P to the boundary S is calculated. This value corresponds to the value of the signed distance function φ at the position of the lattice point P. For lattice points P inside the boundary S, a negative value is given as the distance Lps. For lattice points P other than the lattice points P in the vicinity of the boundary S, an appropriate value larger than the distance Lps from the lattice points P in the vicinity of the boundary S to the boundary S is set as the initial value. At this time, a negative value is set for the lattice points P inside the boundary S, and a positive value is set for the lattice points P outside the boundary S.

[0032] The above equation (11) is discretized to evolve the value of the signed distance function φ over time. The calculation ends when the values ​​of the signed distance function φ at all lattice points P converge. This converged value corresponds to the value of the signed distance function φ at each lattice point P. In this way, the value of the signed distance function φ at each position of multiple lattice points P defined in the analysis space 10 can be determined. A surface that includes multiple lattice points P with approximately the same value of the signed distance function φ is adopted as the boundary of the divided region Ω.

[0033] [Technology for exchanging physical quantities between two adjacent divided regions] Next, a technique for exchanging physical quantities between two adjacent divided regions Ω will be described with reference to FIGS. 3A and 3B.

[0034] FIG. 3A shows two adjacent divided regions Ω 1 and divided area Ω 2 1 is a schematic diagram showing a part of one divided region Ω 1 A plurality of particles 20A are arranged in the other divided region Ω 2 A plurality of particles 20B are arranged in one divided region Ω. 1 The particle 20A in the other divided region Ω 2 The particle in the region Ω is larger than 20B. 1 The diameter of particle 20A in h 1 and the division region Ω 2 The diameter of particle 20B in h 2 The division area Ω 1 The serial number of the particle 20A in the region is denoted as i, and the region Ω 2 The serial number given to the particle 20B in is denoted as j.

[0035] Figure 3B shows two divided regions Ω 1 , Ω 2 2 is a schematic diagram showing the exchange of physical quantities between two divided regions Ω 1 , Ω 2 At the boundary of b1 , Ω b2 Define the buffer region Ω b1 is one of the divided regions Ω 1The other division domain Ω 2 The inwardly expanding buffer region Ω b2 is one of the divided regions Ω 2 The other division domain Ω 1 The area expanded inwards. 1 The expanded buffer area Ω b1 is the divided area Ω 2 Similarly, the divided region Ω 2 The expanded buffer area Ω b2 is the divided area Ω 1 overlaps with a portion of.

[0036] Buffer area Ω b1 , Ω b2 The width Δx of the divided region Ω 1 The diameter h of particle 20A in 1 and divided area Ω 2 Diameter h of particle 20B in 2 The buffer area Ω should be about four times the diameter of the larger one. b1 Physical quantity A of particle 20A in i Instead of solving the governing equations, we use the buffer region Ω b1 The divided area Ω overlapping 2 The physical quantity A is calculated by kernel interpolation of the physical quantity of the particle 20B in the i is calculated using the following formula:

number

number

[0037] [Technology for controlling particles flowing into and out of divided regions] Next, a technique for controlling particles flowing into and out of the divided region Ω will be described with reference to FIGS.

[0038] In simulations using the SPH method, multiple particles move with the flow, so it is necessary to remove or add particles according to the particles flowing out of the divided domain Ω and the particles flowing into the divided domain Ω. If the particle removal and addition processes are not consistent with the physical quantities of the flow field, calculations will fail and analysis accuracy will decrease.

[0039] When the analysis space 10 (Figure 1) is divided into arbitrary shapes, the direction and speed of the flow at the boundary surface of the divided domain Ω cannot be determined in advance. Furthermore, when the flow is unsteady, the direction and speed of the flow at the boundary surface of the divided domain Ω changes from moment to moment. Therefore, a method is required to handle particles flowing into multiple divided domains Ω of arbitrary shapes and particles flowing out of the divided domains Ω in a manner that is consistent with the physical quantities of the flow field.

[0040] Furthermore, in flow fields where the compressibility of the fluid is strongly manifested, the local density of the flow field fluctuates from moment to moment, making it impossible to apply control techniques that are based on incompressible flow fields.

[0041] 4 is a schematic diagram showing a two-dimensional view of a part of the analysis space 10. Two divided regions Ω are defined in the analysis space 10. 1 , Ω 2 is defined. 1 From the adjacent split region Ω 2 Inward buffer region Ω b1 is expanding. Buffer region Ω b1 The outer surface of S 1 It is marked as such.

[0042] Buffer area Ω b1 Surface S1 Along the inside of k These lattice points P k For example, an orthogonal equidistant grid is laid out over the entire analysis space 10, and a buffer region Ω b1 on the surface S 1 Distance from 1 / 2 ×h 1 It can be placed by selecting a lattice point within h 1 is the buffer area Ω 1 is the diameter of the particle (kernel width) at the lattice point P k The interval is the kernel width h 1 In this case, the physical quantities at each grid point are the same as the surface S 1 The above physical quantities can be considered. k Any physical quantity A on k can be expressed by the following equation using the kernel function interpolation formula: where k is the lattice point P k It is a serial number assigned to the

number

[0043] Equation (14) is the divided area Ω 1 Buffer region Ω extending from b1 At a lattice point P along the surface of k The physical quantity of 1 The partitioned region Ω adjacent to 2 This means that the internal physical quantity is calculated by interpolating. Here, the physical quantity Ak to be interpolated is density ρ k and the velocity vector v k It is.

[0044] Surface S at the lattice point positions 1 Normal vector n k and the velocity vector v k Using and, the inner product n k ·v k By checking the sign of the surface S 1 Lattice point Pk In other words, it is possible to distinguish between the part where the fluid is flowing in and the part where the fluid is flowing out at the lattice point P k At the position, a buffer region Ω crosses the surface S1. b1 P is a lattice point into which particles flow and the dot product is positive. k At the position of b1 Particles are flowing out from the center of the tube. The volumetric flow rate is (n k ·v k )h 1 2 It can be estimated that:

[0045] Lattice point P where the inner product is negative k In this case, the buffer area Ω b1 Across the surface of the partition region Ω 1 The number of particles flowing into k / dt can be expressed by the following formula:

number

[0046] Divided area Ω 1 The position, number, and time of particles flowing into the flow path can be calculated by solving the differential equation given by equation (15), but due to errors associated with the discretization of the flow field, inconsistencies occur in the number of particles flowing in. Next, we will explain a method to correct this inconsistency.

[0047] Divided area Ω 1 The volume of V 1 The division domain Ω 1 is defined in the pre-processing stage of the simulation, so the volume V1 is constant throughout the simulation. Using the assumption that the analysis space 10 is always filled with fluid due to the continuity of the flow field, the partitioned domain Ω 1 The volume of particles present in m i / ρ i The sum of the partitioned region Ω 1 Volume V 1 That is, at any time t, the following equation must be satisfied:

number

[0048] From equation (16), the control volume V at time t is 0 The volume of the particle relative to V in (t) can be calculated using the following formula: 0 As the initial state, the partition domain Ω 1 It is advisable to adopt the volume of

number

[0049] Combining equations (15) and (17), for each lattice point P k For this reason, the differential equation that determines the time at which the fluid particles flow in is finally expressed as follows:

number

[0050] Equation (18) is the lattice point P k When particles are flowing into the division region Ω through k When particles are flowing out of the partition region Ω, the number of particles flowing in N k In this way, multiple lattice points P k For each of the divided regions, it is judged whether the particle is flowing in or out of the divided region, and the result of the judgment of inflow and outflow is taken into account and multiple lattice points P k Then, it is decided whether to add a new particle to each of the

[0051] By applying the above-mentioned technique to all divided domains Ω present in the analysis space 10, it is possible to connect flow fields with different spatial resolutions.

[0052] [Simulation device] 5 is a block diagram of the 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.

[0053] 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.

[0054] [Simulation method] Fig. 6 is a flowchart showing the procedure of the simulation method of this embodiment. Each step shown in Fig. 6 is realized by the processing unit 31 (Fig. 5) executing a program stored in the storage unit 33 (Fig. 5).

[0055] First, the processing unit 31 acquires the simulation conditions input to the input unit 30 (step S1). The simulation conditions include information defining the fluid to be simulated, initial conditions, boundary conditions, information for dividing the analysis space 10 into a plurality of divided regions Ω (FIG. 1), and the like. The information defining the fluid includes physical properties such as the density and viscosity of the fluid. The initial conditions include information for arranging a plurality of particles 20 (FIG. 1) in the analysis space 10, information specifying the initial velocity of the particles 20, and the like. The boundary conditions include information specifying the shape and size of the analysis space 10. The information for dividing the analysis space 10 into a plurality of divided regions Ω (FIG. 1) includes information specifying the shape and size of each of the divided regions Ω, and information specifying the kernel width of each of the divided regions Ω.

[0056] When the processing unit 31 acquires the simulation conditions, it divides the analysis space 10 into a plurality of divided regions Ω (step S2). For example, as described with reference to FIG. 2, the analysis space 10 is divided using a signed distance function φ, and the constant value surface of the signed distance function φ is set as the boundary surface of the plurality of divided regions Ω.

[0057] The processing unit 31 defines a buffer region for each of the multiple divided regions Ω (step S3). For example, as shown in FIG. 3B, 1 , Ω 2 One of the divided regions Ω 1 The other division domain Ω 2 Extend the buffer area Ω b1 Similarly, we define the partition domain Ω 2 Divide the area Ω 1 Extend the buffer area Ω b2 Next, the processing unit 31 defines the divided region Ω 1 , Ω 2 and the buffer region Ω b1 , Ω b2 In each of the divided regions Ω, a plurality of particles 20 (FIG. 1) having different diameters (kernel widths) are placed. Here, the size of the particles 20 is specified by the kernel width h in equation (2) (step S4).

[0058] After arranging the plurality of particles 20, the processing unit 31 b1 , Ω b2 For the particles 20, the buffer region Ω is calculated by the method described with reference to FIG. b1 , Ω b2 Then, the processing unit 31 determines physical quantities to be assigned to the particles 20 such as (step S5). After assigning physical quantities to all the particles 20, the processing unit 31 calculates the physical quantities to be assigned to the particles 20, i.e., the velocity and density, by solving the governing equations shown in formulas (8) to (10) for each of the divided regions Ω (step S6).

[0059] Next, the processing unit 31 controls the number of particles for each divided region Ω (step S7). The control of the number of particles is performed by applying the method described with reference to FIG. 4. Specifically, k (Figure 4) using differential equation (18) to calculate the number of particles flowing in, N k is evolved over time, and the number of particles flowing in, N k When the value of reaches the decision threshold, the grid point P k A particle is added to the position of the grid point P. kThe physical quantities of the particles in the divided domain Ω to which the particle belongs are determined by interpolation using a kernel function.

[0060] After controlling the number of particles, the physical quantity and position of each particle in the divided region Ω are updated (step S8). Steps S5 to S8 are repeated until the process of evolving the physical quantity and position over time is completed. When the process of evolving the physical quantity and position of the particle over time is completed, the analysis result is output to the output unit 32 (FIG. 5) (step S9).

[0061] [Excellent Effects of the Examples] Next, the excellent effects of the above embodiment will be described. In the above embodiment, a different spatial resolution, i.e., the diameter of the particles 20, can be set for each of the multiple divided regions Ω (FIG. 1) to perform a simulation. Therefore, the diameter of the particles in the region where high spatial resolution is required can be made small to satisfy the spatial resolution requirement. Also, in the region where high spatial resolution is not required, the diameter of the particles can be made large to reduce the amount of calculation. In this way, the required spatial resolution can be maintained while reducing the amount of calculation.

[0062] In the above embodiment, as described with reference to FIG. 2, a boundary S of any shape can be determined, and the boundary surfaces of a plurality of divided regions Ω can be determined using a signed distance function φ. For example, high spatial resolution is required near the wall surface of the fluid flow field and the surface of an object placed in the flow field. In this way, it is preferable to determine the region requiring high spatial resolution according to the shape of the wall surface placed in the flow field. In this embodiment, by setting the boundary S shown in FIG. 2 according to the shape of the wall surface, it is possible to set the divided region Ω having a shape that reflects the shape of the wall surface. For example, the simulation method according to this embodiment can be applied to the analysis of the flow field around an object with a complex shape or the flow field in a flow channel with a complex shape.

[0063] Furthermore, in the above embodiment, the differential equation (18) is solved to obtain the lattice points P kThe number of particles flowing into the partition region Ω every k is calculated and new particles are added to the divided domain Ω. In this way, the number of particles is controlled so that the sum of the current volume of particles in the divided domain Ω remains constant. This makes it possible to analyze highly compressible flow fields in which the number of particles in the divided domain Ω fluctuates.

[0064] [Verification of effectiveness] In order to verify the excellent effects of the above embodiment, a simulation of a flow field simulating a turbine blade row was performed using the simulation method according to the above embodiment. This simulation will be described below with reference to Figures 7A to 8B.

[0065] 7A and 7B are diagrams showing the initial state of particles arranged in the analysis space 10. The blade chord length Lv was set to 135.7 mm, the pitch Ly was set to 87.59 mm, and the blade axial length Lx was set to 102.5 mm. The cross-sectional shape and dimensions of the blade were the same as those used in the experiment disclosed in "White, AJ, Young, JB and Walters, PT, Experimental validation of condensing flow theory for a stationary cascade of steam turbine blades, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 354, 1704 (1996), pp. 59-88."

[0066] The height of the analysis space 10 was set equal to the pitch Ly of the blade row. The distance from the inflow boundary Sin to the leading edge of the blade, and the distance from the trailing edge of the blade to the outflow boundary Sout were each set equal to the axial length Lx of the blade.

[0067] FIG. 7A shows the kernel width h 0FIG. 7B shows a case where the simulation method according to the above embodiment is applied to divide the analysis space 10 into four divided regions Ω 0 , Ω 1 , Ω 2 , Ω 3 The divided region Ω is shown. 0 , Ω 1 , Ω 2 are the regions within 4 mm, 16 mm, and 40 mm from the blade surface, respectively. 3 is the region from 40 mm upstream from the leading edge of the blade to 40 mm downstream from the trailing edge of the blade. 0 , Ω 1 , Ω 2 , Ω 3 The kernel width h 0 , h 1 , h 2 , h 3 were set to 1 mm, 2 mm, 4 mm, and 8 mm, respectively.

[0068] A uniform inflow condition was applied to the inflow boundary Sin, and a Neumann condition was imposed on the outflow boundary Sout. Periodic boundary conditions were imposed in the vertical direction of the analysis space 10. Adiabatic and slip conditions were imposed on the wall surfaces. The stagnation pressure and stagnation temperature of the flow were set to 41.9 kPa and 378.5 K, respectively. As the initial conditions for the calculation, particles with the kernel width set for each divided domain Ω were placed at equal intervals throughout the entire flow field, the velocity was set to zero, and the stagnation pressure and stagnation temperature were given. The time step size for the analysis was 2.67 × 10 -7 It was seconds.

[0069] The number of particles in the initial state of the calculation was 199,592 under the conditions shown in Fig. 7A, and 47,216 under the conditions shown in Fig. 7B. That is, in this embodiment, the number of particles could be reduced to about 1 / 4 compared to the method in which the analysis space 10 is not divided. The time required for calculation was also reduced to about 1 / 4.

[0070] Figures 8A and 8B are diagrams showing instantaneous values ​​of particle distribution and Mach number distribution when the flow field has reached a steady state. Particles with relatively small Mach numbers are shown in a relatively light color. Figure 8A shows the analysis results when the analysis space 10 is not divided (Figure 7A), and Figure 8B shows the analysis results when it is divided into divided regions Ω (Figure 7B). The flow expands as it passes through the blade row, and becomes supersonic. As can be seen from Figure 8A, a region where the Mach number changes suddenly extends linearly from near the trailing edge of the blade, and it can be seen that a shock wave is generated at the boundary. Figure 8B shows the analysis results when the divided region Ω, which has a relatively small kernel width, is used. 0 , Ω 1 There is a region where the Mach number changes suddenly.

[0071] From the simulation results, it was confirmed that the simulation method according to the above embodiment is capable of efficiently simulating a flow field in which the compressibility of the fluid is strongly manifested, accompanied by a sudden expansion of the fluid and the generation of shock waves.

[0072] In the simulation results shown in FIG. 8B, the divided region Ω 2 , Ω 3 It appears that the change in Mach number is gentler in . This is thought to be the result of the numerical viscosity increasing due to the large kernel width, which in turn blunts the shock wave. If a more detailed simulation is required, it is sufficient to change the setting of the division domain Ω and the kernel width and perform the simulation again.

[0073] The above-mentioned embodiment is merely an example, and the present invention is not limited to the above-mentioned embodiment. For example, it will be obvious to those skilled in the art that various modifications, improvements, combinations, etc. are possible. In the above-mentioned embodiment, the SPH method is used as the particle method, but the method according to the above-mentioned embodiment can also be applied to fluid analysis using other particle methods, such as the dissipative particle method and the MPS method. [Explanation of symbols]

[0074] 10 Analysis space 11. Wall 20, 20A, 20B particles 30 Input section 31 Processing section 32 Output section 33 Storage section

Claims

1. A simulation device that analyzes a fluid flow using a particle method, comprising: an input section for inputting simulation conditions; a processing unit that performs time evolution of each position of a plurality of particles using a particle method based on the simulation conditions inputted to the input unit; Equipped with The processing unit includes: Divide the analysis space into multiple sub-domains, A plurality of particles having different sizes are arranged in each of the plurality of divided regions; determining whether or not to add a new particle to each of the plurality of divided regions so that a sum of the volumes of the particles at the current time in the divided region is constant when the positions of the plurality of particles in each of the plurality of divided regions are evolved over time; A simulation device that adds new particles to a divided region where it has been determined that new particles should be added.

2. The processing unit includes: In the determination of whether or not to add new particles to the plurality of divided regions, 2. The simulation device according to claim 1, wherein, for each of a plurality of lattice points defined along the boundaries of each of the plurality of divided regions, the number of particles per unit time flowing through the lattice point into the divided region is evolved over time, and when the number of flowing-in particles exceeds a judgment threshold, particles are added to the lattice point where the number of flowing-in particles exceeds the judgment threshold.

3. The processing unit includes:

3. The simulation device according to claim 2, wherein in determining whether or not to add new particles to a divided region, it is determined for each of the plurality of lattice points whether particles are flowing in or out of the divided region, and the determination result of inflow or outflow is taken into account to determine whether or not to add new particles to each of the plurality of lattice points.

4. The processing unit includes:

4. The simulation device according to claim 1, wherein, when analyzing behavior of particles using a particle method for each of the plurality of divided regions, when focusing on one of the plurality of divided regions, a buffer region is defined by extending the divided region of interest toward an adjacent divided region, particles having the same size as the particles placed in the divided region of interest are placed in the buffer region, and physical quantities assigned to particles in the buffer region are calculated based on physical quantities assigned to particles in an adjacent divided region that includes the buffer region.

5. A simulation method for analyzing a fluid flow using a particle method, comprising the steps of: Divide the analysis space into multiple sub-domains, A plurality of particles having different sizes are arranged in each of the plurality of divided regions; determining whether or not to add a new particle to each of the plurality of divided regions so that a sum of the volumes of the particles at the current time in the divided region is constant when the positions of the plurality of particles in each of the plurality of divided regions are evolved over time; A simulation method for adding new particles to a divided region where it has been determined that new particles should be added.

6. A step of dividing an analysis space into a plurality of divided domains; a step of arranging a plurality of particles having different sizes for each of the plurality of divided regions; a step of determining whether or not to add a new particle to each of the plurality of divided regions so that a sum of the volumes of the particles at the current time in the divided region becomes constant when the positions of the plurality of particles in each of the plurality of divided regions are evolved over time; a step of adding new particles to the divided area where it is determined that new particles should be added; A program that causes a computer to execute the following.

Images