Simulation device and program

Abstract

To provide a simulation device capable of reducing the calculation load while maintaining the required spatial resolution. [Solution] A simulation device is provided that uses a particle method to analyze the time change of at least one physical quantity that is the object of calculation. A processing unit performs calculations based on simulation conditions input to an input unit. The processing unit has a function of representing an object for which the physical quantity is defined using a plurality of calculation particles with different shapes and sizes, a function of assigning the physical quantity to the plurality of calculation particles, and a function of solving a governing equation that governs the change in the physical quantity for the plurality of calculation particles to evolve the physical quantity over time.

Description

[Technical Field]

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

[0002] In fluid analysis using particle methods such as the SPH (Smoothed Particle Hydrodynamics) method, the spatial resolution of the analysis is determined by the diameter of the calculation particles. The diameter of the calculation particles 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, near the wall of an object, the length scale of fluid motion is represented by the boundary layer thickness, and high spatial resolution is required. On the other hand, in regions away from the wall, the length scale of fluid motion is often larger than the boundary layer thickness. For this reason, analysis can often be performed with lower spatial resolution in regions away from the wall than in regions near the wall.

[0003] Also, in heat transfer analysis, high spatial resolution is required in the direction of a large temperature gradient, whereas analysis can be performed with low spatial resolution in the direction of a small temperature gradient.

[0004] Since numerical analysis must be performed with limited computational resources, it is desirable to ensure the necessary spatial resolution by setting the diameter of the computational particles small in areas where high spatial resolution is required, and to perform analysis by setting the diameter of the computational particles large and reducing the number of computational particles in areas where high spatial resolution is not required.In fluid analysis technology using particle methods, a technology is known that uses physical quantities of the flow field to evaluate whether the spatial resolution at the position of each computational particle is excessive or insufficient, and dynamically divides or combines the computational particles, thereby automatically setting an appropriate spatial resolution according to the flow field (Patent Document 1). [Prior art documents] [Patent documents]

[0005] [Patent Document 1] Japanese Patent Application Publication No. 2024-25512 Summary of the Invention [Problem to be solved by the invention]

[0006] There is a demand for a technique that can further reduce the calculation load compared to the conventional technique of setting an appropriate spatial resolution. An object of the present invention is to provide a simulation device and a program that can reduce the calculation load while maintaining the required spatial resolution. [Means for solving the problem]

[0007] According to one aspect of the present invention, A simulation device that analyzes time variation of at least one physical quantity that is a calculation target using a particle method, a processing unit that performs calculations based on simulation conditions input to the input unit; Equipped with The processing unit A function of representing an object for which the physical quantity is defined by a plurality of calculation particles having different shapes and sizes; a function of assigning the physical quantity to the plurality of computational particles; a function of solving a governing equation governing a change in the physical quantity for the plurality of computational particles and evolving the physical quantity over time; A simulation device is provided having:

[0008] According to another aspect of the present invention, A program that causes a computer to perform a function of analyzing the time change of at least one physical quantity that is a calculation target using a particle method, A function to acquire simulation conditions, A function of representing an object for which the physical quantity is defined by a plurality of calculation particles having different shapes and sizes; a function of assigning the physical quantity to the plurality of computational particles; a function of solving a governing equation governing a change in the physical quantity for the plurality of computational particles and evolving the physical quantity over time; A program for realizing this on a computer is provided. [Effects of the Invention]

[0009] By using multiple computational particles with different shapes and sizes, it is possible to flexibly respond to cases where the required spatial resolution differs depending on the direction. For example, by reducing the size of the computational particles in a direction where high spatial resolution is required, it is possible to perform analysis with high spatial resolution. Furthermore, by relatively increasing the size of the computational particles in a direction where high spatial resolution is not required, it is possible to suppress an increase in the number of computational particles and reduce the computational load. [Brief explanation of the drawings]

[0010] [Figure 1] 1A and 1B are schematic diagrams showing a portion of a wall surface 20 of a channel and a plurality of computational particles 10 positioned near the wall surface 20. FIG. [Figure 2] FIG. 2 is a schematic diagram for explaining the parameters defined for an ellipsoidal computational particle 10. In FIG. [Figure 3] FIG. 3A is a schematic diagram of a wall surface 20 and the s-th computational particle 10 in the vicinity thereof, and FIG. 3B is a schematic diagram of the computational particle 10 as viewed from the direction perpendicular to the wall surface 20. As shown in FIG. [Figure 4] FIG. 4 is a schematic diagram of a computational particle 10 for explaining a method for evaluating the spatial resolution at the position of the computational particle 10. In FIG. [Figure 5] FIG. 5 is a schematic diagram of a computational particle 10 for explaining a method for evaluating the spatial resolution at the position of the computational particle 10. In FIG. [Figure 6] FIG. 6 is a schematic diagram showing an example in which one computational particle 10 is divided into elliptical computational particles 10 with different aspect ratios. [Figure 7] FIG. 7A is a schematic diagram showing an example of a plurality of computational particles before merging, and FIG. 7B is a schematic diagram showing an example of a plurality of computational particles after merging. [Figure 8] FIG. 8 is a block diagram of a simulation device according to the first embodiment. [Figure 9]FIG. 9 is a flowchart showing the processing procedure executed by the processing unit 31 of the simulation device according to the first embodiment. [Figure 10] FIG. 10 is a schematic diagram of the analysis space when the simulation was actually performed. [Figure 11] 11A and 11B are schematic diagrams showing the distribution of calculation particles near the wall surface 20 at the start of the analysis (time t=0 seconds) and at the end of the analysis (time t=0.1 seconds), respectively. [Figure 12] FIG. 12 is a graph showing the velocity distribution in the width direction at the end of the analysis. [Figure 13] 13A and 13B are schematic diagrams showing an example of the distribution of a plurality of computational particles 10 when an analysis is performed using the simulation devices according to the comparative example and the second embodiment, respectively. [Figure 14] FIG. 14 is a flowchart showing the procedure of a simulation performed by the simulation device according to the second embodiment. [Figure 15] FIG. 15 is a perspective view of a heat transfer member to be simulated and a schematic diagram showing the shape of calculation particles. [Figure 16] FIG. 16 is a graph showing a comparison of the simulation results with the theoretical analysis results. DETAILED DESCRIPTION OF THE INVENTION

[0011] 1A and 1B, an overview of a simulation performed by a simulation device according to a first embodiment will be described. The simulation device according to the first embodiment represents a fluid using a plurality of computational particles, assigns a physical quantity to each of the plurality of computational particles, places the plurality of computational particles in an analysis space that models a flow path, and analyzes the behavior of the plurality of computational particles.

[0012] 1A and 1B are schematic diagrams showing a part of a wall surface 20 of a flow channel and a plurality of computational particles 10 located near the wall surface 20. Fig. 1A shows the computational particles 10 used in analysis by a simulation device according to a comparative example, and Fig. 1B shows the computational particles 10 used in analysis by a simulation device according to a first embodiment.

[0013] Generally, when analyzing a flow with a high Reynolds number, high spatial resolution is required near the wall surface 20. To obtain high spatial resolution, the diameter of the calculation particles 10 near the wall surface 20 must be small, as shown in Figure 1A.

[0014] However, the spatial resolution in the direction parallel to the wall surface 20 does not need to be as high as the spatial resolution in the direction perpendicular to the wall surface 20. For this reason, in the example shown in FIG. 1A, the spatial resolution in the direction parallel to the wall surface 20 is excessive.

[0015] 1B, by making the computational particles 10 flattened spheroids and aligning their short axes parallel to the direction perpendicular to the wall surface 20, the spatial resolution in the direction perpendicular to the wall surface 20 is increased, and excessive spatial resolution in the parallel direction is avoided. As shown in FIG. 1B, by making the computational particles 10 flattened spheroids, the number of computational particles 10 can be reduced, and the computational load can be reduced.

[0016] In the first embodiment described below, a flow field is analyzed using a plurality of computational particles 10 with different shapes and sizes, as shown in FIG. 1B.

[0017] [Overview of particle methods] Next, we will briefly explain the general SPH 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 uses a discretization model in which the superposition of the differential coefficients of the kernel function is the spatial derivative.

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

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

number

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

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

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

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

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

number

[0025] [Governing equations of flow] One way to analyze high Reynolds number compressible flows with limited computational resources is to use Reynolds-Averaged Navier-Stokes (RANS) analysis. This is a technique 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.

[0026] As a formulation method for RANS analysis, for example, the Shear Stress Transport (SST) model (FR Menter, "Two-equation eddy-viscosity turbulence models for engineering applications", AIAA Journal Vil.32, No.8, pp.1598-1605 (1994), https: / / doi.org / 10.2514 / 3.12149) can be used.

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

number

number

number

number

number

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

number

number

number

[0029] Other parameters σ k , σ ω , σ ω2 , β * are constants used in the SST model. γ, β, F l is a variable calculated using k and ω. Since the SST model is widely used, we will not explain these constants here.

[0030] By discretizing the spatial differential terms in equations (8) to (12) using equation (4) and solving them numerically, compressible flows at high Reynolds numbers can be analyzed using the SPH method. The governing equations discretized by the SPH method are shown below.

number

number

number

number

number

number

number

number

[0031] [Method using non-spherical computational particles] Next, a method using non-spherical computational particles 10 such as spheroids or superellipsoids as shown in FIG. 1B will be described. When calculating the kernel function and its derivative shown in Equation (2) and Equation (6), the absolute value of the distance vector between two calculation particles is used. In other words, the values ​​of the kernel function and its derivative do not depend on the direction. This is equivalent to using spherical calculation particles.

[0032] To use non-spherical computational particles, a direction vector of the computational particle is defined for each computational particle, and the kernel function and its derivatives are calculated using the vector. Next, a method for calculating the kernel function and its derivatives when ellipsoidal computational particles are used in two-dimensional space will be described with reference to FIG.

[0033] FIG. 2 is a schematic diagram for explaining parameters defined for an elliptical computational particle 10. The central position vector of the s-th computational particle 10 is expressed as x s The length of the major axis of the elliptical calculation particle 10 is denoted as h a,s and the length of the minor axis is h b,s In this case, the kernel function of the computational particle 10 with respect to an arbitrary position vector x and the derivative of the kernel function can be calculated as follows:

number

number

number

[0034] [How to set the direction vector of the calculation particle] Next, referring to FIGS. 3A and 3B, the unit direction vector s of the non-spherical calculation particle s 3A is a schematic diagram of a wall surface 20 and the s-th computational particle 10 in its vicinity. The unit normal vector at the point on the wall surface 20 where the distance from the center of the computational particle 10 to the wall surface 20 is the shortest is defined as n s The velocity vector of the calculation particle 10 is denoted as v s In general fluid analysis, high spatial resolution is required in the direction perpendicular to the wall surface 20.

[0035] Therefore, when using an elliptical computational particle 10, it is effective to set the minor axis in a direction perpendicular to the wall surface 20. In this case, the major axis of the computational particle 10 is the unit normal vector n s The direction of this long axis is the unit direction vector s of the computational particle 10. s It can be defined as follows.

[0036] In three-dimensional space, the unit normal vector n s In the plane perpendicular to s The direction of the

[0037] 3B is a schematic diagram of the computational particle 10 as viewed from the direction perpendicular to the wall surface 20. As an example, the unit direction vector s s Calculate the direction of particle 10's velocity vector v s Specifically, the unit direction vector s of the computational particle 10 can be set to be aligned with the direction of sis the unit normal vector n of the nearest wall 20 s and a velocity vector v assigned to the calculation particle 10. s That is, when viewed from the vertical direction of the wall surface 20, the unit direction vector s s The direction of the velocity vector v s overlaps with the direction of.

[0038] [Method for evaluating spatial resolution in the position of calculated particles] Next, a method for evaluating the spatial resolution at the position of a computational particle will be described with reference to Fig. 4 and Fig. 5. Fig. 4 and Fig. 5 are schematic diagrams of a computational particle 10 for explaining a method for evaluating the spatial resolution at the position of the computational particle 10. Hereinafter, an example will be described in which the computational particle 10 is assumed to be elliptical, and the spatial resolution in the direction of the minor axis of the computational particle 10 is evaluated.

[0039] The spatial resolution at the position of the sth calculation particle 10 is determined by the kernel width h b,s and the shear strength γ in the direction of the minor axis at the position of the calculated particle 10 b,s Reynolds number Re defined by the dot γ,s It can be evaluated using the Reynolds number Re γ,s is defined by the following formula:

number

[0040] In addition, for the computational particle 10 in the vicinity of the wall surface 20 shown in FIG. 5, the dimensionless distance y s + The dimensionless distance y s + is defined by the following formula:

number

number

number

[0041] For the calculation particle 10 present in the vicinity of the wall surface 20, the wall friction stress can be calculated approximately as follows using the velocity gradient at the position of the s-th calculation particle 10.

number

[0042] At any time during the analysis, the value of Equation (27) or Equation (28) for each computational particle 10 can be calculated to evaluate whether the spatial resolution of each computational particle 10 is sufficient or not. For example, Re γ,s is greater than a predetermined first decision threshold, or y s + If Re is greater than a predetermined second threshold, it can be determined that the spatial resolution is insufficient. For example, 10 can be used as the first threshold and 3 can be used as the second threshold. γ,s is smaller than a predetermined third decision threshold, or y s + is smaller than a predetermined fourth decision threshold, it can be determined that the spatial resolution is excessive. For example, 2.5 can be adopted as the third decision threshold and 1 can be adopted as the fourth decision threshold. By making the third decision threshold smaller than the first decision threshold and the fourth decision threshold smaller than the second decision threshold, it is possible to prevent the computational particles 10 from repeatedly dividing and merging excessively.

[0043] [Calculation particle division method] Next, a method for dividing computational particles determined to have insufficient spatial resolution will be described with reference to Fig. 6. When dividing the computational particles 10, the physical quantities must be preserved. Furthermore, the shape, direction, number, and arrangement of the computational particles 10 after division are arbitrary.

[0044] FIG. 6 is a schematic diagram showing an example of dividing one computational particle 10A into elliptical computational particles 10 with different aspect ratios. When dividing a circular computational particle 10A with a diameter D, it can be divided into two elliptical computational particles 10B having a minor axis parallel to the direction in which spatial resolution is insufficient. The length of the minor axis is set to 1 / 2 of the diameter D of the computational particle 10 before division, and the length of the major axis is set to be equal to the diameter D of the computational particle 10A before division. By this operation, one circular computational particle 10A is divided into two elliptical computational particles 10B with an aspect ratio of 2.

[0045] When further dividing an elliptical computational particle 10B with an aspect ratio of 2, it can be divided into two elliptical computational particles 10C having a minor axis parallel to the direction in which spatial resolution is lacking (the direction of the minor axis). The length of the minor axis of the divided computational particle 10C is set to 1 / 2 of the length D / 2 of the minor axis of the computational particle 10B before division, and the length of the major axis is set to be equal to the length D of the major axis of the computational particle 10B before division. By this operation, the elliptical computational particle 10B with an aspect ratio of 2 is divided into two elliptical computational particles 10C with an aspect ratio of 4.

[0046] By repeating this operation, the elliptical computational particle 10C with an aspect ratio of 4 is divided into two elliptical computational particles 10D with an aspect ratio of 8. The length of the minor axis of this computational particle 10D is D / 8. By using computational particles with a large aspect ratio, the effect of reducing the number of computational particles can be enhanced.

[0047] If a computational particle 10D with an aspect ratio of 8 is judged to have insufficient resolution, repeating the same procedure will result in the aspect ratio of the computational particle after division becoming 16. If the aspect ratio of the computational particle becomes too large, the analysis will become numerically unstable.

[0048] When dividing a computational particle 10D with an aspect ratio of 8 in the minor axis direction, the computational particle 10D is first divided in the major axis direction. This results in a computational particle 10E with a minor axis length of D / 8, a major axis length of D / 2, and an aspect ratio of 4. By dividing this computational particle 10E in the minor axis direction, four computational particles 10F with a minor axis length of D / 16, a major axis length of D / 2, and an aspect ratio of 8 are obtained. By this procedure, one computational particle 10D with an aspect ratio of 8 is divided into four computational particles 10F with a minor axis length of 1 / 2 without increasing the aspect ratio.

[0049] In this way, if the aspect ratio of the computational particle after division becomes too large, the length of the minor axis of the computational particle can be shortened sufficiently to obtain the desired spatial resolution while maintaining the aspect ratio at a certain value or less by dividing the computational particle also in the major axis direction. For example, when dividing an elliptical computational particle with an aspect ratio of 8 or more in the minor axis direction, it is recommended to divide the computational particle using the above procedure while maintaining the aspect ratio at 8 or less.

[0050] [Method of combining computational particles] Next, a method for merging computational particles will be described with reference to FIGS. 7A and 7B. By combining two calculation particles that are judged to have excessive spatial resolution, it is possible to save the calculation resources required for the analysis.

[0051] 7A is a schematic diagram showing an example of a plurality of computational particles before merging, and FIG. 7B is a schematic diagram showing an example of a plurality of computational particles after merging. In FIG. 7A, computational particles determined to have excessive spatial resolution are hatched. The other computational particle (hereinafter referred to as the nearest computational particle) that exists closest to each of the plurality of computational particles is determined. The arrows in FIG. 7A point from each computational particle to the nearest computational particle.

[0052] A pair of computational particles having a relationship of nearest neighbors to each other (hereinafter referred to as a nearest pair) is extracted. In Fig. 7A, computational particles 10a and 10b have a relationship of nearest neighbors, and computational particles 10d and 10e have a relationship of nearest neighbors.

[0053] When two computational particles having a nearest pair relationship are both determined to have excessive spatial resolution, the two computational particles of the nearest pair are merged. In the example shown in FIG. 7A, computational particles 10a and 10b are merged to become one computational particle 10ab as shown in FIG. 7B. Computational particles 10d and 10e have a nearest pair relationship, but since one computational particle 10e is not determined to have excessive spatial resolution, the computational particles 10d and 10e are not merged. Computational particle 10c is determined to have excessive spatial resolution, but is not merged because it does not form a nearest pair.

[0054] The method for determining the position and shape of the computational particles after merging is arbitrary. The position of the computational particle after merging can be determined, for example, using a weighted average value based on the mass of the two computational particles before merging. The size of the computational particle after merging can be determined so that the mass and volume do not change before and after merging. The aspect ratio of the computational particle after merging can be given, for example, a weighted average value based on the mass of the aspect ratio before merging. The minor axis direction of the computational particle after merging is determined by the unit normal vector n of the nearest wall surface 20 (Figure 3A). s The velocity vector v can be made parallel to s The magnitude and direction of (Figure 3A) can be determined so that the momentum vector does not change before and after merging. The unit direction vector s of the computational particle after merging s(Fig. 3A) is the unit normal vector n s and velocity vector v s It can be determined from the following.

[0055] [When to evaluate whether spatial resolution is sufficient or not] Next, the timing for evaluating whether the spatial resolution is excessive or insufficient will be described. If the computational particles are repeatedly divided and merged in a short number of calculation steps, the computational particles filling the analysis space may overlap, or there may be areas where no computational particles exist, resulting in an inconsistency in the distribution of the computational particles and making the calculation numerically unstable. To avoid this, the spatial resolution should be evaluated, and the process of dividing and merging the computational particles should be performed at intervals of at least a certain number of calculation steps. The number of calculation steps in this case, N r is, for example, the maximum kernel width h max Using the fluid's representative sound speed c0 and time step size Δt, it can be given as follows:

number

[0056] In fluid analysis methods using particle methods, a neighbor particle list is used, which stores the serial numbers of other calculation particles that exist in the vicinity of each calculation particle. It is efficient to perform the evaluation process of excess or deficiency in spatial resolution and the division and merging process of calculation particles at the timing of updating this neighbor particle list. In practice, for example, when the number of calculation steps N has elapsed since the previous evaluation process of excess or deficiency in spatial resolution was performed, r After the above has passed and the neighbor particle list is updated, the process of evaluating whether the spatial resolution is excessive or insufficient and the process of dividing and combining the calculation particles can be performed.

[0057] 8 is a block diagram of a simulation device according to the first embodiment. The simulation device according to the first 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.

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

[0059] Next, the processing procedure executed by the processing unit 31 (FIG. 8) of the simulation device according to the first embodiment will be described with reference to Fig. 9. Fig. 9 is a flowchart showing the processing procedure executed by the processing unit 31 of the simulation device according to the first embodiment.

[0060] First, the processing unit 31 acquires the simulation conditions input by the user through the input unit 30 (step S1). The simulation conditions include information defining the fluid to be simulated, initial conditions, boundary conditions, an evaluation index for evaluating whether the spatial resolution is excessive or insufficient, the shape of the calculation particles, the maximum aspect ratio of the calculation particles, a division method for the calculation particles, and information for specifying the maximum kernel width. The information defining the fluid to be simulated includes physical property values ​​such as the density and viscosity of the fluid. The initial conditions include a method for specifying the initial physical property values ​​of the fluid, an initial value for the kernel width of the calculation particles, and information for arranging multiple calculation particles in the analysis space. The boundary conditions include information specifying the shape and size of the analysis space.

[0061] Next, the processing unit 31 places a plurality of computational particles in the analysis space based on the input simulation conditions (step S2). Next, the governing equations (16) to (21) for each of the computational particles are numerically solved to update the physical quantities and positions of each of the computational particles (step S3). Furthermore, the unit direction vector s s (Figures 2, 3A, and 3B) are updated (step S4).

[0062] Next, the processing unit 31 determines whether to update the neighbor particle list (step S5). For example, if the cumulative movement distance of at least one calculation particle from the position of each calculation particle when the neighbor particle list was most recently updated to the current time exceeds a predetermined value, it is determined that the neighbor particle list should be updated. If the cumulative movement distances of all calculation particles are equal to or less than the predetermined value, it is determined that the neighbor particle list should not be updated.

[0063] When it is determined that the neighbor particle list should be updated, the number of calculation steps N is determined based on the previous evaluation of the spatial resolution. r It is determined whether or not a predetermined number of calculation steps N (equation (32)) has elapsed since the previous evaluation of the spatial resolution. r If the time has not elapsed, the neighbor particle list is updated without evaluating whether the spatial resolution is excessive or insufficient (step S11).

[0064] In step S6, a predetermined number of calculation steps N r If it is determined that the time has elapsed, the spatial resolution of each calculation particle is evaluated (step S7). This evaluation is performed using equation (27) or equation (28).

[0065] Next, the computational particles determined to have insufficient spatial resolution are divided (step S8). The computational particles are divided using the procedure described with reference to FIG. 6, for example. After that, the computational particles determined to have excessive spatial resolution are merged (step S9). The computational particles are merged using the procedure described with reference to FIGS. 7A and 7B. After that, the shape and posture of the divided or merged computational particles are updated (step S10). Specifically, h in Equation (24) a,s , h b,s , and the unit direction vector s s Then, the neighbor particle list is updated (step S11).

[0066] If it is determined in step S5 that the neighbor particle list is not to be updated, or after the neighbor particle list is updated in step S11, it is determined whether or not the calculation is completed (step S12). If the calculation is to be continued, the procedure from step S3 is repeated. If the calculation is to be completed, the calculation result is output to the output unit 32 (FIG. 8) (step S13).

[0067] Next, with reference to FIGS. 10 to 12, the results of a simulation of a flow between parallel plates performed using the simulation device according to the first embodiment will be described.

[0068] Figure 10 is a schematic diagram of the analysis space. The analysis space is a plane with a length L = 0.025 m and a width W = 0.1 m. Both sides in the width direction are in contact with the wall surface 20. A fluid flows in the length direction. Periodic boundary conditions are applied to both ends in the length direction (the inlet and outlet ends of the fluid). A no-slip condition is applied as the boundary condition for the wall surface 20. Air is assumed as the working fluid, and an initial temperature of 293 K and an initial pressure of 101.3 kPa are given to all calculation particles. As an initial velocity condition, a velocity according to the distance from the wall surface 20 is given to the calculation particles so that the average flow velocity becomes a parabolic distribution of 170 m / s.

[0069] The viscosity coefficient was calculated using Sutherland's law, and 10 times the value was given to the viscosity coefficient. 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. The time step size of the analysis was 7.28 × 10 -7 The time was set to 0 seconds. The flow field was allowed to evolve over time, with time 0 seconds as the initial state, and the velocity distribution at time 0.1 seconds was output. The shape of the calculation particles was set to elliptical, with a maximum aspect ratio of 8. For comparison, an analysis was performed with the shape of the calculation particles limited to circles.

[0070] 11A and 11B are schematic diagrams showing the distribution of computational particles near the wall surface 20 at the start of the analysis (time t = 0 seconds) and the end of the analysis (time t = 0.1 seconds), respectively. At the start of the analysis, computational particles with a kernel width of 0.005 m are regularly distributed in the length and width directions. At the end of the analysis, it can be seen that the shorter the length of the minor axis of the computational particles closer to the wall surface 20, the larger the velocity gradient is near the wall surface 20, requiring higher spatial resolution. Therefore, the computational particles are divided in the direction perpendicular to the wall surface 20 during the analysis, and the dimensions in the direction perpendicular to the wall surface 20 become smaller.

[0071] In the comparative example in which the shape of the calculation particles is limited to a circle, not only the dimension in the direction perpendicular to the wall surface 20 but also the dimension in the direction parallel to the wall surface 20 becomes smaller when the calculation particles are divided. Therefore, when the spatial resolution in the direction perpendicular to the wall surface 20 is set to be the same in the comparative example and the first embodiment, the number of calculation particles is greater in the comparative example than in the first embodiment. In an actual simulation, when analysis was performed using the simulation device in the first embodiment, the number of calculation particles at the end of the analysis was about 0.24 times that in the comparative example.

[0072] Figure 12 is a graph showing the velocity distribution in the width direction at the end of the analysis. The horizontal axis represents the position in the width direction in units of mm, and the vertical axis represents the flow velocity in units of m / s. The solid line in the graph shows the results of analysis performed using the simulation device according to the first embodiment, and the dashed line shows the results of analysis performed using the method according to the comparative example in which the calculation particles are limited to a circular shape. The velocity distributions of both cases are nearly identical. In other words, when analysis is performed using the simulation device according to the first embodiment, the number of calculation particles can be reduced by approximately 0.24 times compared to the comparative example in which the calculation particles are limited to a circular shape, and analysis results of equivalent accuracy can be obtained.

[0073] Next, the excellent effects of the first embodiment will be described. In the first embodiment, the shape of the calculation particles is changed depending on the direction of the required spatial resolution. Specifically, the size of the calculation particles in the direction requiring high spatial resolution is made smaller than the size of the calculation particles in the direction not requiring high spatial resolution. This makes it possible to perform analysis at the desired spatial resolution while suppressing an increase in the number of calculation particles.

[0074] Furthermore, changes in the shape and size of the calculation particles are performed automatically during the analysis. This means that the user does not need to have the skill to set the size of the calculation particles based on an appropriate spatial resolution. If the user sets the size of the calculation particles based on the spatial resolution and performs an analysis, and then finds that the set size of the calculation particles is too large, the user must reset the size of the calculation particles and redo the analysis. In the first embodiment, this eliminates the need to reset the dimensions of the calculation particles and redo the analysis.

[0075] In the simulation described with reference to FIGS. 10 to 12, an example was handled in which a fluid is sandwiched between a pair of wall surfaces 20, but the simulation device according to the first embodiment can also handle other fluid flows. For example, it can also handle so-called external flows that are not surrounded by wall surfaces, such as flows passing through a car body or flows around a baseball. In particular, it is possible to efficiently analyze fluid flows in which the required spatial resolution varies depending on the location. For example, it is possible to efficiently analyze cases in which the fluid flow is affected by the surface of an object.

[0076] Next, a simulation device according to a second embodiment will be described with reference to Figures 13A to 16. Below, a description of the configuration common to the simulation device according to the first embodiment described with reference to Figures 1A to 12 will be omitted.

[0077] The simulation device according to the first embodiment analyzes fluid flow. In contrast, the simulation device according to the second embodiment analyzes heat flow (heat transfer analysis) using the SPH method. It is preferable to set the size and shape of the calculation particles appropriately depending on the heat transfer phenomenon to be analyzed. High spatial resolution is required in areas and directions where the temperature gradient is large. On the other hand, analysis can be performed with low spatial resolution in areas and directions where the temperature gradient is small.

[0078] 13A and 13B are schematic diagrams showing an example of the distribution of a plurality of computational particles 10 when performing an analysis using the simulation devices according to the comparative example and the second embodiment, respectively. A plurality of computational particles 10 are distributed within an analysis space (a member formed of a heat-conductive material (hereinafter referred to as a heat-conductive member)) surrounded by a wall surface 20.

[0079] In the comparative example shown in FIG. 13A, spherical computational particles 10 of the same size are distributed. The size of the computational particles 10 needs to be set according to the largest temperature gradient that occurs in the heat transfer member to be analyzed. Generally, a region where a large temperature gradient occurs is a partial region of the heat transfer member. In other regions, the temperature gradient is small. In the comparative example, computational particles of the same size as those in the region with a large temperature gradient are placed in the region with a small temperature gradient, so the analysis is performed with excessive spatial resolution.

[0080] In a second example shown in FIG. 13B , the multiple computational particles 10 include computational particles 10 of different shapes or sizes. For example, the computational particles 10 are spherical, spheroidal (prolate spheroid and oblate spheroid), or superellipsoidal. A plurality of small spherical computational particles 10A are arranged near the wall surface 20. A spheroidal computational particle 10B is arranged in an area slightly far from the wall surface 20. The major axis of the computational particle 10B is longer than the diameter of the computational particle 10A and is parallel to the nearest wall surface 20. The minor axis of the computational particle 10B is approximately equal to or slightly longer than the diameter of the computational particle 10A and is approximately perpendicular to the nearest wall surface 20. A plurality of spherical computational particles 10C larger than the computational particle 10A are arranged in an area further away from the wall surface 20.

[0081] As an example, consider the case where a heat source exists on the wall surface 20, and a heat transfer phenomenon in which the temperature in the analysis space changes due to heat given to each wall surface 20 is analyzed. The largest temperature gradient occurs on the wall surface 20, and the direction of the temperature gradient is perpendicular to the wall surface 20. In the second embodiment, the desired spatial resolution is ensured by making the computational particles 10A near the wall surface 20 small.

[0082] The minor axes of the calculation particles 10B in the region slightly far from the wall surface 20 are perpendicular to the wall surface 20, so that the spatial resolution in the direction of the large temperature gradient can be satisfied. For the direction parallel to the wall surface 20, where high spatial resolution is not required, the number of calculation particles 10 is reduced by increasing the dimensions of the calculation particles 10B. Furthermore, for the region where high spatial resolution is not required, the diameter of the calculation particles 10C is increased, so that the number of calculation particles 10 can be reduced.

[0083] Next, the second embodiment will be described in detail. In heat transfer analysis, it is necessary to calculate the second derivative of temperature. In the SPH method, the second derivative with respect to any scalar quantity A can be calculated using the following formula:

number

[0084] [Governing equations for heat transfer phenomena] Next, we will explain the governing equations used to analyze the heat transfer phenomenon in solids. The heat transfer equations shown below are used to analyze the heat transfer phenomenon.

number

number

[0085] The kernel function and the kernel function derivatives for the non-spherical computational particles 10 may be calculated using the method described with reference to FIG.

[0086] Next, the procedure of the simulation performed by the simulation device according to the second embodiment will be described with reference to Fig. 14. Fig. 14 is a flowchart showing the procedure of the simulation performed by the simulation device according to the second embodiment.

[0087] First, the processing unit 31 (FIG. 8) acquires the simulation conditions input to the input unit 30 (FIG. 8) (step SA1). The simulation conditions include information defining the heat transfer member to be simulated, boundary conditions, and initial conditions. The information defining the heat transfer member includes physical property values ​​such as the density, thermal conductivity, and specific heat of the heat transfer member. The boundary conditions include information specifying the shape and size of the analysis space (heat transfer member) and information specifying the shape and size of the calculation particles. The initial conditions include a method for specifying the initial physical property values ​​of the heat transfer member, a kernel width value for the calculation particles, and information for arranging multiple calculation particles in the analysis space.

[0088] When the simulation conditions are acquired, the processing unit 31 places a plurality of calculation particles in the analysis space based on the acquired simulation conditions (step SA2). Thereafter, the governing equation (Equation (35)) is solved for each of the plurality of calculation particles, and the physical quantity (temperature T) assigned to the calculation particle is updated (step SA3). Step SA3 is repeatedly executed until a calculation end condition is satisfied (step SA4). When the calculation is completed, the calculation result is output (step SA5).

[0089] Next, the excellent effects of the second embodiment will be described. 13A and 13B, the number of computational particles 10 required for analysis can be reduced. In this case, by reducing the dimensions of the computational particles 10 in a direction where high spatial resolution is required, for example, a direction in which the temperature gradient is steep, and increasing the dimensions of the computational particles 10 in a direction perpendicular to that direction, sufficient spatial resolution can be ensured even if the number of computational particles 10 is reduced.

[0090] Next, the results of an actual simulation performed using the simulation device according to the second embodiment will be described with reference to FIGS.

[0091] FIG. 15 is a perspective view of the heat transfer member to be simulated and a schematic diagram showing the shape of the computational particle 10. The heat transfer member 25 to be simulated is a long rectangular parallelepiped with a length of 10 m, a height of 0.2 m, and a width (depth direction in FIG. 15) of 0.2 m. The physical properties of the heat transfer member 25 were those of iron. Specifically, the density was 8900 kg / m, the thermal conductivity coefficient was 372 W / (m·K), and the specific heat was 491 J / (kg·K). As boundary conditions, an isothermal wall condition with a temperature of 373.15 K was applied to one end face of the heat transfer member 25 (the left end face in FIG. 15, hereinafter referred to as the "heat source"), and an adiabatic wall condition was applied to the other face. As an initial condition, the temperature of the entire heat transfer member 25 was set to 273.15 K.

[0092] The computational particles 10 placed in the region up to 0.5 m from the heat source were spherical with a diameter of 0.01 m, and the computational particles 10 placed in the region 0.5 m to 1 m from the heat source were spherical with a diameter of 0.02 m. The computational particles 10 placed in the region 1 m to 2 m from the heat source were prolate spheroids with a major axis length of 0.04 m and a minor axis length of 0.02 m. The computational particles 10 placed in the region 2 m to 3.33 m from the heat source were prolate spheroids with a major axis length of 0.08 m and a minor axis length of 0.04 m. The computational particles 10 placed in the region 3.33 m to 5 m from the heat source were prolate spheroids with a major axis length of 0.12 m and a minor axis length of 0.04 m. The computational particles 10 placed in the region 5 m or more from the heat source were prolate spheroids with a major axis length of 0.16 m and a minor axis length of 0.04 m. In addition, the major axis of each of the computational particles 10 was set parallel to the longitudinal direction of the heat transfer member 25 .

[0093] Figure 16 is a graph comparing the simulation results with the theoretical analysis results. The horizontal axis represents the distance from the heat source in units of [m], and the vertical axis represents the temperature in units of [K]. The circle, triangle, square, and inverted triangle symbols in the graph represent the simulation results for temperature when the elapsed time from the start of the calculation is 0, 5053, 25266, and 39921 seconds, respectively. The solid line in the graph represents the theoretical analysis results for temperature when the elapsed time corresponds to the simulation results indicated by each symbol.

[0094] It can be seen that the simulation results closely match the theoretical analysis results within a distance of 2 m or less from the heat source. Even within a distance of 2 m or more from the heat source, the simulation results closely match the theoretical analysis results when the temperature gradient is gentle. When the distance from the heat source is far and the temperature gradient is relatively large, the difference between the simulation results and the theoretical analysis results becomes slightly larger. This is because the longitudinal dimension of the calculation particles was increased. To reduce the difference between the simulation results and the theoretical analysis results, the longitudinal dimension of the calculation particles can be reduced.

[0095] The simulations shown in Figures 15 and 16 confirmed that for regions and directions with a gentle temperature gradient, analysis can be performed with high accuracy even if the dimensions of the calculation particles are made smaller than for regions or directions with a steep temperature gradient.

[0096] When spherical calculation particles with a diameter of 0.01 mm are placed throughout the entire analysis space, the number of calculation particles is approximately 400,000. In contrast, when the shape and size of the calculation particles are changed as shown in Figure 15, the number of calculation particles is 27,250. In this way, when performing a simulation using the simulation device according to the second embodiment, it is possible to reduce the number of calculation particles compared to when performing a simulation by placing only spherical calculation particles of the same size, and it is possible to perform a highly accurate simulation.

[0097] Next, a modification of the second embodiment will be described. In the second embodiment, in equation (34), it is assumed that the thermal conductivity k, density ρ, and specific heat C do not change over time. However, the thermal conductivity k, density ρ, and specific heat C may be made a function of temperature T, and may be changed over time.

[0098] Next, a simulation device according to a third embodiment will be described. In the first embodiment, a fluid flow is analyzed, and in the second embodiment, a heat transfer phenomenon is analyzed. A simulation device according to the third embodiment more generally analyzes the time change of at least one physical quantity that is the object of calculation using a particle method. In the first embodiment, the physical quantity that is the object of calculation is the position, velocity, density, pressure, etc. that are assigned to each calculation particle, and in the second embodiment, the physical quantity that is the object of calculation is the temperature that is assigned to each calculation particle.

[0099] In the third embodiment, when the physical quantity to be calculated evolves over time, a governing equation that governs the change in the physical quantity is solved. In the first embodiment, this governing equation is the governing equation for compressible flow shown in equations (8) to (12), and in the second embodiment, this governing equation is the heat transfer equation shown in equation (34).

[0100] As described above, an object to be analyzed can be represented by a plurality of calculation particles with different shapes and sizes, and the particle method can be used to numerically calculate the time change of physical quantities.

[0101] In the above embodiment, the SPH method has been described as an example of a particle method. However, the technique of changing the shape and size of calculation particles according to the required spatial resolution can also be applied to other particle methods, such as the Moving Particle Semi-implicit (MPS) method.

[0102] The above-described embodiments are merely examples, and the present invention is not limited to the above-described embodiments. For example, it will be obvious to those skilled in the art that various modifications, improvements, combinations, etc. are possible. [Explanation of symbols]

[0103] 10, 10A, 10B, 10C calculation particles 20 Wall 25 Heat transfer material 30 Input section 31 Processing section 32 Output section 33 Storage section

Claims

1. A simulation device that analyzes time variation of at least one physical quantity that is a calculation target using a particle method, a processing unit that performs calculations based on simulation conditions input to the input unit; Equipped with The processing unit A function of representing an object for which the physical quantity is defined by a plurality of calculation particles having different shapes and sizes; a function of assigning the physical quantity to the plurality of computational particles; a function of solving a governing equation governing a change in the physical quantity for the plurality of computational particles and evolving the physical quantity over time; A simulation device having:

2. the object is a fluid, The simulation device according to claim 1 , wherein the physical quantities include a position, a velocity, a density, and a pressure imparted to each of the plurality of computational particles.

3. The simulation device according to claim 2 , wherein the physical quantities include internal energy and temperature imparted to each of the plurality of computational particles.

4. 4. The simulation device according to claim 2, wherein the processing unit further has a function of determining whether or not to divide a computational particle into a plurality of computational particles based on at least one of the physical quantities assigned to each of the plurality of computational particles, and dividing the computational particle determined to be divided into the plurality of computational particles.

5. The simulation device according to claim 4 , wherein the processing unit changes the shape of the computational particle when dividing the computational particle determined to be divided into a plurality of computational particles.

6. 5. The simulation device according to claim 4, wherein the processing unit further has a function of determining whether or not to combine the plurality of computational particles into one computational particle based on at least one of the physical quantities assigned to each of the plurality of computational particles, and combining the plurality of computational particles determined to be combined into one computational particle.

7. The processing unit further has a function of calculating an index for evaluating spatial resolution for each of the plurality of calculation particles, and determining whether the spatial resolution is insufficient or excessive based on the index; In the determination of whether or not to divide the calculation particle, it is determined that the calculation particle has insufficient spatial resolution, and the calculation particle is to be divided; 7. The simulation device according to claim 6, wherein in the determination of whether or not to combine the computational particles, it is determined that two computational particles whose spatial resolution is determined to be excessive are to be combined.

8. At least some of the calculation particles are flat or prolate spheroids, 4. The simulation device according to claim 2, wherein the processing unit further has a function of changing the attitude of the computational particles based on at least one of the physical quantities assigned to the computational particles.

9. The simulation device according to claim 1 , wherein the object is a member formed of a heat-conductive material, and the physical quantity includes a temperature applied to the plurality of calculation particles.

10. A program that causes a computer to perform a function of analyzing the time change of at least one physical quantity that is a calculation target using a particle method, A function to acquire simulation conditions, A function of representing an object for which the physical quantity is defined by a plurality of calculation particles having different shapes and sizes; a function of assigning the physical quantity to the plurality of computational particles; a function of solving a governing equation governing a change in the physical quantity for the plurality of computational particles and evolving the physical quantity over time; A program that makes the computer realize this.

Images