Simulation method, simulation device, and program

Abstract

To provide a simulation method that can keep analysis accuracy high, and reduce an amount of calculation in an analysis using a particle method.SOLUTION: A simulation method is configured to: represent a fluid by a plurality of calculation particles; give each of the plurality of calculation particles a physical quantity; arrange the plurality of calculation particles in an analysis space; develop the physical quantity given to each of the plurality of calculation particles and a position of the plurality of calculation particles over time by resolving a dominant equation; evaluate the excess or deficiency of spatial resolution relying upon a size of the plurality of calculation particles as for each of the plurality of calculation particles in line with a state of a flow field represented by the plurality of calculation particles; and adjust the spatial resolution in line with an evaluation result of the excess or deficiency of the spatial resolution.SELECTED DRAWING: Figure 2

Description

[Technical Field]

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

[0002] A simulation method is known 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 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, spatial resolution is determined by the diameter of the particles. The particle diameter must be set small enough relative 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 the 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 fluid near the wall, the particle diameter must be sufficiently smaller than the thickness of the boundary layer.

[0005] To perform a highly accurate analysis across the entire analysis space, the particle diameter must be reduced according to the thickness of the boundary layer, even in regions far from the wall. This requires a large number of particles to be placed in the analysis space, which results in an increase in the amount of calculations.

[0006] An object of the present invention is to provide a simulation method, a simulation device, and a program that can maintain high analytical accuracy and reduce 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 method comprising: representing a fluid using 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; and evolving the physical quantity assigned to each of the plurality of computational particles and positions of the plurality of computational particles over time by solving a governing equation, evaluating, for each of the plurality of computational particles, whether or not a spatial resolution is sufficient or insufficient, depending on the size of the plurality of computational particles, according to the state of the flow field represented by the plurality of computational particles; A simulation method is provided for adjusting the spatial resolution depending on the result of the evaluation of whether the spatial resolution is excessive or insufficient.

[0008] According to another aspect of the present invention, A simulation device that analyzes fluid flow using a particle method, an input section for inputting simulation conditions; a processing unit that uses a particle method to evolve the positions of each of a plurality of particles over time based on the simulation conditions input to the input unit; Equipped with The processing unit a step of representing a fluid using 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, 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; a step of evaluating whether a spatial resolution, which depends on the size of the plurality of calculation particles, is excessive or insufficient for each of the plurality of calculation particles according to a state of a flow field represented by the plurality of calculation particles; The procedure for adjusting the spatial resolution according to the results of the evaluation of whether the spatial resolution is excessive or insufficient. A simulation device is provided that executes the above.

[0009] According to yet another aspect of the present invention, A program that causes a computer to execute a simulation that analyzes a fluid flow using a particle method, a step of representing a fluid using 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, 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; a step of evaluating whether a spatial resolution, which depends on the size of the plurality of calculation particles, is excessive or insufficient for each of the plurality of calculation particles according to a state of a flow field represented by the plurality of calculation particles; The procedure for adjusting the spatial resolution according to the results of the evaluation of whether the spatial resolution is excessive or insufficient. A program for executing the above is provided. [Effects of the Invention]

[0010] By adjusting the spatial resolution depending on the results of the evaluation of whether the spatial resolution is excessive or insufficient, analysis can be performed with an appropriate spatial resolution throughout the entire analysis space. [Brief explanation of the drawings]

[0011] [Figure 1] FIG. 1 is a block diagram of a simulation device according to an embodiment. [Figure 2] FIG. 2 is a flowchart showing the procedure of the simulation method according to this embodiment. [Figure 3] FIG. 3 is a schematic diagram showing an analytical space and a plurality of computational particles arranged in the analytical space. [Figure 4]FIG. 4 is a flowchart showing the steps of a method for evaluating spatial resolution. [Figure 5] FIG. 5 is a schematic diagram showing the positional relationship between the wall surface and the calculation particles to be evaluated. [Figure 6] FIG. 6 is a flowchart showing the procedure of the spatial resolution adjustment method executed in step SA6 (FIG. 2). [Figure 7] FIG. 7A is a schematic diagram of a computational particle before division, and FIG. 7B is a schematic diagram of two computational particles generated by dividing one computational particle. [Figure 8] FIG. 8A is a schematic diagram of four computational particles generated by dividing one computational particle, and FIG. 8B is a schematic diagram of eight computational particles generated by dividing one computational particle. [Figure 9] FIG. 9A is a schematic diagram showing an example of the distribution of a plurality of computational particles, and FIG. 9B is a schematic diagram showing an example of the distribution of computational particles after two computational particles are merged. [Figure 10] Figure 10A is a schematic diagram showing the distribution of multiple computational particles before the computational particles are divided, Figure 10B is a schematic diagram showing the distribution of multiple computational particles after one computational particle is divided, and Figure 10C is a schematic diagram showing the distribution of multiple computational particles after a certain time step has passed from the state shown in Figure 10B. [Figure 11] Figure 11A is a schematic diagram showing the distribution of multiple computational particles before two computational particles are merged, Figure 11B is a schematic diagram showing the distribution of multiple computational particles after two computational particles are merged, and Figure 11C is a schematic diagram showing the distribution of multiple computational particles at a point in time after a certain time step has passed from the state shown in Figure 11B. [Figure 12] FIG. 12 is a flowchart showing the procedure of a simulation method according to another embodiment. [Figure 13] FIG. 13 is a schematic diagram of the analysis model used in the actual simulation. [Figure 14] 14A and 14B are schematic diagrams showing the distribution of calculation particles at the start and end of the calculation, respectively. [Figure 15]FIG. 15 is a graph showing the velocity distribution at the end of the calculation and the velocity distribution in the initial state. DETAILED DESCRIPTION OF THE INVENTION

[0012] Before describing an embodiment of the present invention, we will briefly explain the general SPH (Smoothed Particle Hydrodynamics) method. The SPH method is a Lagrangian-type numerical calculation method that discretizes a continuum using a set of multiple computational 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 its periphery by a kernel function W, and the density field is given as a superposition of these. That is, the position x of the sth computational particle is given as s Density ρ at s 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 known as 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 also be referred to as the particle diameter.

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

number

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

number

number

[0018] Here, ∇W(x sm , 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 (6).

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 of flow] Next, we will explain the governing equations when analyzing compressible flows. As a means of analyzing high Reynolds number compressible flows with limited computational resources, for example, Reynolds-Averaged Navier-Stokes (RANS) analysis can be used. This is a technology that analyzes the mean field of high Reynolds number flows with limited computational resources by decomposing the governing equations of the flow into a mean field and a fluctuating field using Reynolds decomposition, modeling the contribution of the fluctuating field to the mean field, and numerically solving only the governing equations of the mean field.

[0022] The RANS analysis can be formulated using the Shear Stress Transport (SST) model, as explained in FR MENTER, "Two-equation eddy-viscosity turbulence models for engineering applications," AIAA Journal Vol. 32 No. 8, (1994), pp. 1598-1605, _eprint: https: / / doi.org / 10.2514 / 3.12149.

[0023] The governing equations for compressible flow using the SST model can be written, for example, as follows:

number

number

number

number

number

[0024] σ ij represents the viscous stress tensor and is given by the following equation:

number

number

[0025] where μ t is the eddy viscosity coefficient. The pressure p can be calculated using, for example, the following equation of state for an ideal gas:

number

[0026] The other variables in Equations (9) to (13) are constants used in the SST model. Since the SST model is widely used, a description of these constants will be omitted here.

[0027] By discretizing and numerically solving equations (9) to (13) using equation (4), compressible flows at high Reynolds numbers can be analyzed using the SPH method.

[0028] The governing equations discretized by the SPH method are written as follows:

number

number

number

number

number

number

[0029] Artificial advection term δv i For example, the following formula can be used:

number

number

[0030] [Spatial resolution] Next, we will explain the difference between excessive and insufficient spatial resolution when analyzing flow fields using particle methods. In fluid analysis using particle methods such as the SPH method, the spatial resolution of the analysis is determined by the diameter of the calculation particles (e.g., kernel width). In general, the length scale of fluid motion varies depending on the location. For example, near the wall of an object, the length scale of fluid motion is represented by the boundary layer thickness, and a relatively high spatial resolution is required.

[0031] On the other hand, in regions far from the wall, the length scale of fluid motion is often larger than the boundary layer thickness, so analyses can be performed with lower spatial resolution in regions far from the wall than in regions near the wall.

[0032] If calculation particles with a large diameter are placed in an area where analysis must be performed with high spatial resolution, the spatial resolution of the analysis will be insufficient. Conversely, if calculation particles with a small diameter are placed in an area where analysis can be performed with low spatial resolution, the spatial resolution of the analysis will be excessive. In the embodiment described below, the excess or deficiency of the spatial resolution of the analysis can be reduced, and analysis can be performed with the necessary and sufficient spatial resolution.

[0033] [Simulation procedure for example] The simulation device and simulation method according to this embodiment will be described with reference to FIGS.

[0034] 1 is a block diagram of a simulation device according to this embodiment. The simulation device according to this 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 instructions (commands) and the like are input from the 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.

[0035] 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, changes in the physical quantities of the particle system over time, etc. 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, etc.

[0036] FIG. 2 is a flowchart showing the procedure of the simulation method according to this embodiment. First, the processing unit 31 acquires the simulation conditions input to the input unit 30 (step SA1). The simulation conditions include, for example, information defining the fluid to be simulated, boundary conditions, initial conditions, etc. The information defining the fluid includes physical property values ​​such as the density and viscosity of the fluid. The boundary conditions include information specifying the shape and size of the analysis space. The initial conditions include a method for specifying the initial physical property values ​​of the fluid, an initial value of the kernel width of the calculation particles, and information for arranging multiple calculation particles in the analysis space.

[0037] The simulation conditions further include information for specifying an evaluation index for whether the spatial resolution is excessive or insufficient, a method for dividing the calculation particles, and a maximum kernel width. These pieces of information are referenced in the simulation of this embodiment, which will be described below.

[0038] Next, the processing unit 31 places a plurality of computational particles in the analysis space based on the acquired simulation conditions (step SA2).

[0039] 3 is a schematic diagram showing an analysis space 10 and a plurality of computational particles 20 arranged in the analysis space 10. As an example, the analysis space 10 has a rectangular parallelepiped shape, and conditions such as a no-slip condition and a periodic boundary condition are imposed on each face of the analysis space 10.

[0040] Next, the governing equations shown in Equations (17) to (22) are numerically solved for each of the multiple computational particles 20 to calculate the physical quantities assigned to each of the computational particles 20, and the physical quantities and positions of each of the computational particles 20 are updated (step SA3).

[0041] Next, it is determined whether or not to evaluate the spatial resolution (step SA4). For example, the processing unit 31 determines to evaluate the spatial resolution when a certain number of time steps have passed since the most recent time step in which the spatial resolution was evaluated.

[0042] When it is determined that the spatial resolution should be evaluated, the processing unit 31 evaluates the spatial resolution (step SA5). For example, if the diameter of the calculation particle 20 is too large compared to the spatial resolution required for the analysis, it is determined that the spatial resolution is insufficient. Conversely, if the diameter of the calculation particle 20 is too small compared to the spatial resolution required for the analysis, it is determined that the spatial resolution is excessive.

[0043] After evaluating whether the spatial resolution is excessive or insufficient, the processing unit 31 adjusts the spatial resolution based on the evaluation result (step SA6). For example, in a region where the spatial resolution is determined to be insufficient, the spatial resolution is increased by making the diameter of the calculation particles 20 smaller than the current diameter. At this time, the number of calculation particles 20 is increased in response to the reduction in the diameter of the calculation particles 20. In a region where the spatial resolution is determined to be excessive, the spatial resolution is decreased by making the diameter of the calculation particles 20 larger than the current diameter. At this time, the number of calculation particles 20 is decreased in response to the increase in the diameter of the calculation particles 20.

[0044] If it is determined in step SA4 that the spatial resolution is not to be evaluated, steps SA5 and SA6 are not executed.

[0045] The processes from step SA3 onwards are repeatedly executed until an analysis termination condition is met (step SA7). The analysis termination condition is specified, for example, by the number of time steps to be executed. When the analysis is completed, the analysis results are output to the output unit 32 (step SA8).

[0046] [Evaluation of spatial resolution] Next, an example of a method for evaluating the spatial resolution executed in step SA5 (FIG. 2) will be described with reference to FIGS.

[0047] 4 is a flowchart showing the steps of a method for evaluating spatial resolution, which is performed for each of a plurality of computational particles 20.

[0048] First, it is determined whether the computational particle 20 to be evaluated is close to a wall surface (step SB1). Here, among the wall surfaces of the analytical space 10 shown in Fig. 3, wall surfaces to which a periodic boundary condition is applied are not included in the wall surfaces to be determined. For example, the determination is made for wall surfaces to which a no-slip condition is applied.

[0049] 5 is a schematic diagram showing the positional relationship between the wall surface 11 and the calculation particle 20s to be evaluated. An xyz Cartesian coordinate system is defined in which the surface of the wall surface 11 is the xz plane and the normal direction of the wall surface 11 is the y direction. The diameter of the calculation particle 20s is defined as h s and the distance from the wall 11 to the center of the calculation particle 20s is denoted as y s As an example, y s <h s The computational particle 20s that satisfies the above condition is determined to be a particle located near the wall surface 11. The component of the fluid velocity vector parallel to the wall surface 11 is defined as v h In the vicinity of the wall surface 11, the velocity v h becomes larger.

[0050] For the calculation particle 20s that is not close to the wall, the Reynolds number Re defined by the shear strength γ at the center position of the calculation particle 20s is γs is calculated using the following formula (step SB2).

number

number

[0051] The processing unit 31 calculates the Reynolds number Re at the position of the calculation particle 20s. γs After calculating the Reynolds number Re γsIt is determined whether the Reynolds number Re is greater than 10 (step SB4). γs If the Reynolds number Re is greater than 10, it is determined that the spatial resolution at the position of the calculation particle 20s is insufficient (step SB6). γs When the Reynolds number Re is 10 or less, the processing unit 31 γs It is determined whether the Reynolds number Re is smaller than 2.5 (Step SB8). γs If is smaller than 2.5, it is determined that the spatial resolution at the position of the calculation particle 20s is excessive (step SB10). In this specification, the spatial resolution at the position of the calculation particle may be simply referred to as the "spatial resolution of the calculation particle."

[0052] When the computational particle 20s is close to the wall surface 11, the processing unit 31 calculates the dimensionless distance y s + is calculated using the following formula (step SB3).

number

number

number

number

[0053] Calculation particle 20s for the dimensionless distance y s + Once calculated, the processing unit 31 calculates the dimensionless distance y s +It is determined whether or not the dimensionless distance y is greater than 3 (step SB5). s + If y is greater than 3, it is determined that the spatial resolution of the position of the calculation particle 20s is insufficient (step SB7). s + If y is 3 or less, the processing unit 31 calculates the dimensionless distance y s + It is determined whether or not the dimensionless distance y is smaller than 1 (step SB9). s + If is smaller than 1, the processing unit 31 determines that the spatial resolution of the computational particles 20s is excessive (step SB11).

[0054] The above step of determining whether the spatial resolution is excessive or insufficient is repeated for all the calculation particles 20 (step SB12).

[0055] [Spatial resolution adjustment method] Next, an example of a method for adjusting the spatial resolution executed in step SA6 (FIG. 2) will be described with reference to FIGS. 6 to 9B.

[0056] 6 is a flowchart showing the procedure of the method for adjusting the spatial resolution executed in step SA6 (FIG. 2). This adjustment is executed for each of the plurality of computational particles 20.

[0057] Based on the evaluation result performed in step SA5, the processing unit 31 determines whether the spatial resolution of the computational particle 20 of interest is excessive (step SC1). If it is determined that the spatial resolution is not excessive, the processing unit 31 determines whether the spatial resolution of the computational particle 20 is insufficient (step SC5). If the spatial resolution of the computational particle 20 is insufficient, the processing unit 31 divides the computational particle 20 (step SC6) to generate multiple computational particles with smaller kernel widths.

[0058] When dividing one computational particle 20 into multiple computational particles 20, it is necessary to preserve physical quantities such as the mass and volume of the computational particle 20. Furthermore, the number and arrangement of the multiple computational particles 20 after division are arbitrary.

[0059] Fig. 7A is a schematic diagram of computational particle 20a before division. Fig. 7B is a schematic diagram of two computational particles 20a1 and 20a2 generated by dividing one computational particle 20a. Fig. 8A is a schematic diagram of four computational particles 20a1 to 20a4 generated by dividing one computational particle 20a. Fig. 8B is a schematic diagram of eight computational particles 20a1 to 20a8 generated by dividing one computational particle 20a.

[0060] First, as shown in FIG. 7B, 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 is.

[0061] The mass m2, volume V2, and kernel width h2 of each of the divided computational particles 20a1 and 20a2 can be calculated using the following equations.

number

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

number

[0063] Next, as shown in FIG. 8A, 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 by h2, the mass by m2, and the volume by 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 follows.

[0064] The mass m2, volume V2, and kernel width h2 of each of the divided computational particles 20a1 to 20a4 can be calculated using the following formulas.

number

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

number

[0066] Next, as shown in FIG. 8B, 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 by h2, the mass by m2, and the volume by 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 follows.

[0067] The mass m2, volume V2, and kernel width h2 of each of the divided computational particles 20a1 to 20a8 can be calculated using the following formulas.

number

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

number

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

[0070] 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 using a similar method.

[0071] If it is determined in step SC1 (FIG. 6) that the spatial resolution is excessive, it is determined whether or not there are any partner particles that form a nearest neighbor pair (step SC2). A "nearest neighbor pair" refers to a pair of computational particles 20 that have a relationship of nearest neighbor to each other.

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

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

[0074] If it is determined in step SC2 (FIG. 6) that there is a partner computational particle 20 that forms a nearest pair, it is determined whether the spatial resolution of the partner particle that forms the nearest pair is excessive (step SC3). If it is determined that the spatial resolution of the partner particle is excessive, the two computational particles 20 that form the nearest pair are merged (step SC4) 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, if the spatial resolutions of the two computational particles 20 that form the nearest pair are both excessive, the two computational particles 20 are merged.

[0075] Of the computational particles 20 shown in Fig. 9A, 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. In Fig. 9A, the computational particles 20a, 20b, 20c, and 20d are determined to have excessive spatial resolution, and are hatched.

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

[0077] 9B 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 larger computational particle 20ab. In FIG. 9B, the merged computational particle 20ab is indicated by relatively light hatching.

[0078] 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 , kernel width is h 11 , mass m 11 , the volume is V 11 and the position of the other computational particle 20b is denoted by x 12 , kernel width is h 12, mass m 12 , the volume is V 12 The mass m2 and kernel width h2 of the combined computational particle 20ab (FIG. 9B) can be calculated using the following equations:

number

[0079] 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 weighted average based on mass as follows:

number

[0080] Similarly, the arbitrary physical quantity A2 of the computational particle 20ab after merging is calculated by the physical quantity A of the computational particles 20a and 20b before merging. 11 , A 12 It can be calculated using the following formula:

number

number

[0081] The procedure from step SC1 is repeated until the process of adjusting the spatial resolution is completed for all the computational particles 20 (step SC7).

[0082] Next, the excellent effects of this embodiment will be described. In this embodiment, the spatial resolution is evaluated in step SA5 shown in Figure 2, and the spatial resolution is adjusted in step SA6, so that analysis can be performed at an appropriate spatial resolution throughout the entire analysis space. For example, analysis is performed at a high spatial resolution near the wall, allowing sufficient information about the flow field to be obtained. Furthermore, analysis is performed at a relatively low spatial resolution in areas far from the wall, reducing the computational load.

[0083] One possible method for performing an analysis with different spatial resolutions depending on the location within the analysis space is to divide the analysis space into multiple regions with different spatial resolutions and place computational particles of different sizes in each divided region. However, this method requires the exchange of physical quantities between the divided regions, which places a heavy computational load on the calculations. Increasing the number of divisions further increases the computational load. Furthermore, in complex flow fields, when computational particles move across divided regions, inconsistencies in the law of conservation of mass can occur, making the analysis numerically unstable. Furthermore, because the division of the analysis space is performed manually by the operator, there is a risk that an inappropriate spatial resolution will be set depending on the operator's skill.

[0084] In contrast to this, in this embodiment, the analysis space is not divided into a plurality of regions, and therefore the above-mentioned problem does not arise.

[0085] In this embodiment, the thresholds used as the criteria for judgment in steps SB4 and SB8 shown in Fig. 4 are set to 10 and 2.5, and the thresholds used as the criteria for judgment in steps SB5 and SB9 are set to 3 and 1, but other values ​​may be used as these thresholds. However, it is preferable to set the threshold used as the criteria for judging that the spatial resolution is excessive lower than the threshold used as the criteria for judging that the spatial resolution is insufficient, thereby creating a difference between the two. Setting the thresholds in this manner prevents the computational particles 20 from frequently repeating division and merging.

[0086] Furthermore, if the computational particles 20 repeatedly split and merge during a small number of time steps, inconsistencies may occur in the distribution of the multiple computational particles 20 that fill the analysis space 10 (FIG. 3), making the calculation numerically unstable. To avoid this, it is preferable to execute the process of adjusting the spatial resolution (step SA6) at a fairly long time interval. The number of time steps corresponding to this time interval is set to N. r and the maximum kernel width is h max When the representative sound speed of the fluid is denoted as c0 and the time step width of the calculation is denoted as Δt, the number of time steps N r It is recommended to set it as follows:

number

[0087] [Reducing errors caused by particle division and merging] Next, a method for reducing errors that may occur due to division and merging of the computational particles 20 will be described with reference to FIGS. 10A to 11C.

[0088] 10A is a schematic diagram showing the distribution of a plurality of computational particles 20 before the computational particles 20a are divided. In FIG. 10A, the computational particles 20a to be divided are hatched. Before division, the plurality of computational particles 20, including the computational particles 20a, are distributed almost uniformly.

[0089] 10B is a schematic diagram showing the distribution of a plurality of computational particles 20 after dividing the computational particle 20a (FIG. 10A). One computational particle 20a (FIG. 10A) is divided into two computational particles 20a1 and 20a2. In FIG. 10B, the two divided computational particles 20a1 and 20a2 are hatched.

[0090] By dividing one computational particle 20a (FIG. 10A) into two computational particles 20a1 and 20a2, a bias occurs in the distribution of the plurality of computational particles 20 including the two computational particles 20a1 and 20a2.

[0091] 11A is a schematic diagram showing the distribution of the plurality of computational particles 20 before the two computational particles 20a and 20b are merged. In FIG. 11A, the two computational particles 20a and 20b before merging are hatched. Before merging, the plurality of computational particles 20, including the computational particles 20a and 20b, are distributed almost uniformly.

[0092] 11B is a schematic diagram showing the distribution of multiple computational particles 20 after the computational particles 20a and 20b (FIG. 11A) are merged. Two computational particles 20a and 20b (FIG. 11A) are merged to generate one computational particle 20ab. In FIG. 11B, the computational particle 20ab generated by the merger is hatched.

[0093] When two computational particles 20a and 20b (FIG. 11A) are combined to generate one computational particle 20ab, a bias occurs in the distribution of the multiple computational particles 20 including the combined computational particle 20ab.

[0094] As shown in FIGS. 10B and 11B, when the distribution of the computational particles 20 becomes biased, the artificial advection term δv i By the action of (a), an artificial velocity is added to each of the computational particles 20 in a direction that eliminates the bias in distribution. The arrows shown in Fig. 10B and Fig. 11B indicate the direction of the artificially added velocity.

[0095] FIG. 10C is a schematic diagram showing the distribution of a plurality of computational particles 20 at a time when a certain time step has passed from the state shown in FIG. 10B. FIG. 11C is a schematic diagram showing the distribution of a plurality of computational particles 20 at a time when a certain time step has passed from the state shown in FIG. 11B. The artificial advection term δv in Equation (17) i By this action, the bias in the distribution of the plurality of computational particles 20 is eliminated.

[0096] The motion of the computational particles 20 brought about by this artificial advection term is unrelated to the fluid motion. Therefore, when the bias in the distribution of the computational particles 20 is eliminated, a distortion occurs in the distribution of the physical quantities that the flow field originally had. Below, we will explain how to suppress this distortion.

[0097] The governing equations shown in equations (17) to (22) include a Lagrangian derivative term. In the method described below, this Lagrangian derivative term is modified to include an artificial advection term, thereby reducing the distortion of physical quantities.

[0098] The artificial advection term δv in Eq. (17) i In order to suppress distortion of physical quantities due to unphysical motion caused by the wind, the Lagrangian derivative dA / dt, which takes into account the artificial advection term, is defined as follows: Note that DA / Dt is the normal Lagrangian derivative without taking into account the artificial advection term.

number

[0099] The second term on the right-hand side of equation (43) can be rewritten as follows for any scalar quantity A:

number

[0100] The second term on the right-hand side of equation (43) can be rewritten as follows for any vector quantity A:

number

[0101] By rewriting equations (17) to (22) using equation (43), the governing equations to be solved by the SPH method are written as follows:

number

number

number

number

number

number

[0102] Equations (47) to (51) are governing equations obtained by modifying the Lagrangian derivative terms included in equations (18) to (22) to include an artificial advection term. By solving the governing equations shown in equations (46) to (51), the physical quantities change in a direction that corrects the distortion of the physical quantities caused by adding the artificial advection term. This makes it possible to reduce the distortion of the physical quantities caused by the division and merging of the computational particles 20 and the action of the artificial advection term.

[0103] [When to evaluate whether spatial resolution is sufficient or not] Next, a simulation method according to another embodiment in which the timing for evaluating whether the spatial resolution is excessive or insufficient is optimized will be described with reference to Fig. 12. Fig. 12 is a flowchart showing the procedure of the simulation method according to this embodiment. The following description will focus on the differences from the flowchart shown in Fig. 2.

[0104] In fluid analysis techniques using particle methods, analysis using a neighbor particle list in which, for each computational particle 20, the numbers of other computational particles 20 existing in the vicinity of that computational particle 20 are stored is widely used. The process of updating this neighbor particle list is executed at regular intervals. As shown in FIG. 12, for each time step, it is determined whether or not to update the neighbor particle list (step SA10). In the time step in which the neighbor particle list is to be updated, the neighbor particle list is updated (step SA11), and in other time steps, the neighbor particle list is not updated.

[0105] In the method shown in Fig. 2, it is determined whether or not to evaluate the spatial resolution for each time step (step SA4). In contrast, in the embodiment shown in Fig. 12, it is determined whether or not to evaluate the spatial resolution for each time step before updating the neighbor particle list (step SA4) only at the time step when the neighbor particle list is updated.

[0106] The interval for updating the neighbor particle list is generally the number of time steps N calculated by Eq. (42). r In step SA4, the number of time steps N from the time when the spatial resolution was most recently evaluated is shorter. r If this time has elapsed, it may be determined that the spatial resolution should be evaluated for excess or deficiency.

[0107] The neighbor particle list needs to be updated after the computational particles 20 are split or combined. Therefore, by splitting or combining the computational particles 20 at the same time as updating the neighbor particle list, the computational efficiency can be improved.

[0108] [Actual simulation results] Next, with reference to Figures 13 to 15, the results of a simulation performed using the simulation method shown in Figure 12 will be described. In this simulation, the governing equations of Equations (46) to (51) were used to analyze the flow between parallel plates.

[0109] Figure 13 is a schematic diagram of the analysis model. The analysis space 10 is a rectangular parallelepiped. The dimension in the flow direction (left-right direction in Figure 13) is 0.025 m, the distance between a pair of parallel wall surfaces 11 (height dimension in Figure 13) is 0.1 m, and the dimension in the width direction perpendicular to the flow is 0.025 m. Periodic boundary conditions are imposed as boundary conditions in the flow direction and width direction. A no-slip condition is imposed as a boundary condition for the wall surfaces 11.

[0110] The initial kernel width of the computational particles 20 was set to 0.005 m, and multiple computational particles 20 were placed at equal intervals in the analysis space 10. Air was assumed to be the working fluid, and an initial temperature value of 293 K and an initial pressure value of 101.3 kPa were given to all computational particles 20 in the analysis space 10. As the initial velocity condition, an initial velocity value was given to each computational particle 20 according to the distance from the wall 11 so that the average flow velocity was a parabolic distribution of 170 m / s. The viscosity coefficient was calculated using Sutherland's law, and a value 100 times this was given. During the analysis, a pressure gradient was applied in the flow direction so that the average velocity of the flow field was 170 m / s.

[0111] The time step size of the analysis is 7.28×10 -7 The time was set to 0 s. The flow field was allowed to evolve over time, with the initial state set to 0.1 s, and the velocity distribution value at time 0.1 s was output. For comparison, an analysis was also performed in which the kernel width of all calculation particles20 was fixed to 0.0005 m, in accordance with the maximum spatial resolution required for the analysis.

[0112] 14A and 14B are schematic diagrams showing the distribution of the computational particles 20 at the start and end of the calculation, respectively. At the start of the calculation, multiple computational particles 20 with a kernel width of 0.005 m are aligned at equal intervals. At the end of the calculation, it can be seen that the closer to the wall surface 11 the region, the smaller the computational particles 20 are distributed. This is because the velocity gradient is large near the wall surface 11, requiring higher spatial resolution, and so the computational particles 20 are divided and their size becomes smaller during the calculation. At the end of the calculation, the number of computational particles 20 was approximately 12,000. In a comparative example in which the kernel width of the computational particles 20 was fixed at 0.0005 m, the number of computational particles 20 was approximately 500,000.

[0113] Figure 15 is a graph showing the velocity distribution at the end of the calculation and the velocity distribution in the initial state. The horizontal axis of Figure 15 represents the flow velocity in units of m / s, and the vertical axis represents the distance from one wall surface 11 in units of m. The thin solid line in Figure 15 represents the velocity distribution in the initial state, and the thick solid line and dashed line represent the velocity distribution at the end of the calculation obtained by the simulation methods of this example and the comparative example, respectively.

[0114] It can be seen that the velocity distribution obtained by the simulation method of this embodiment is almost the same as the velocity distribution obtained by the simulation method of the comparative example. In other words, it was confirmed that this embodiment can reduce the number of calculation particles 20 to about 1 / 40 of that in the comparative example, and can perform analysis with sufficiently high resolution over the entire flow field.

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

[0116] 10 Analysis space 11 Wall 20, 20a, 20a1~20a8, 20ab, 20b, 20c, 20d, 20e, 20s calculation particles 30 Input section 31 Processing section 32 Output section 33 Storage section

Claims

1. A simulation method comprising: representing a fluid using 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; and evolving the physical quantity assigned to each of the plurality of computational particles and positions of the plurality of computational particles over time by solving a governing equation, evaluating, for each of the plurality of computational particles, whether or not a spatial resolution is sufficient or insufficient, depending on the size of the plurality of computational particles, according to the state of the flow field represented by the plurality of computational particles; A simulation method that adjusts spatial resolution according to the results of an evaluation of whether the spatial resolution is excessive or insufficient.

2. The fluid being analyzed is in contact with the wall, The evaluation of whether the spatial resolution is excessive or insufficient is performed based on the Reynolds number of the flow field or the distance from the wall surface at each position of the plurality of calculation particles; 2. The simulation method according to claim 1, wherein the spatial resolution is adjusted by dividing or combining the plurality of calculation particles.

3. adding an artificial advection term to the governing equation that acts in a direction to reduce distortion in the distribution of the plurality of computational particles caused by dividing or merging the plurality of computational particles; The simulation method according to claim 2 , wherein the governing equation to which the artificial advection term is added is solved during time evolution.

4. 4. The simulation method according to claim 3, wherein the governing equation includes a Lagrangian derivative term, and the Lagrangian derivative term is modified to take into account the artificial advection term, and the governing equation is solved, thereby changing the physical quantity in a direction that corrects distortion of the physical quantity for each of the plurality of calculation particles caused by adding the artificial advection term.

5. A simulation device that analyzes fluid flow using a particle method, an input section for inputting simulation conditions; a processing unit that uses a particle method to evolve the positions of each of a plurality of particles over time based on the simulation conditions input to the input unit; Equipped with The processing unit a step of representing a fluid using 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, 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; a step of evaluating whether a spatial resolution, which depends on the size of the plurality of calculation particles, is excessive or insufficient for each of the plurality of calculation particles according to a state of a flow field represented by the plurality of calculation particles; The procedure for adjusting the spatial resolution according to the results of the evaluation of whether the spatial resolution is excessive or insufficient. A simulation device that performs the above.

6. A program that causes a computer to execute a simulation that analyzes a fluid flow using a particle method, a step of representing a fluid using 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, 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; a step of evaluating whether a spatial resolution, which depends on the size of the plurality of calculation particles, is excessive or insufficient for each of the plurality of calculation particles according to a state of a flow field represented by the plurality of calculation particles; The procedure for adjusting the spatial resolution according to the results of the evaluation of whether the spatial resolution is excessive or insufficient. A program that executes the following.

Images