Simulation method, simulation device, and program

Abstract

To provide a simulation method that can analyze a flow of fluid contacting with a wall surface without arranging a wall particle for reproducing a wall surface.SOLUTION: A simulation method is configured to: express flowing fluid contacting with a wall surface by a plurality of particles; determine an inter-particle wall surface interaction between the plurality of particles and the wall surface, and an inter-particle interaction among the plurality of particles; and solve an equation of motion controlling a motion of the plurality of particles as for each of the plurality of particles, and thereby make positions of the plurality of particles and velocities thereof temporally evolve. Upon solving the equation of motion, regarding the particle in which a distance to the wall surface is equal to or less than a first distance set as a simulation condition of the plurality of particles, the simulation method is configured to act attenuation force receiving from the wall surface and random force in accordance with a temperature of the wall surface, besides force by the inter-particle interaction and inter-particle wall surface interaction, and make the position of the particle and the velocity thereof temporally evolve.SELECTED DRAWING: Figure 4

Description

[Technical Field]

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

[0002] In a conventional method for analyzing the flow of a fluid in contact with a wall surface using molecular dynamics, the wall surface is represented by multiple wall particles, and the interaction between the wall particles and the fluid particles is taken into consideration (see Patent Document 1). By controlling the temperature of the wall particles, the temperature of the fluid particles near the wall surface is reproduced. [Prior art documents] [Patent documents]

[0003] [Patent Document 1] Japanese Patent Application Laid-Open No. 2007-219831 Summary of the Invention [Problem to be solved by the invention]

[0004] Conventional simulation methods require the placement of wall particles in addition to fluid particles, which increases the number of particles to be calculated. In addition, the spatial resolution of the wall shape is limited by the particle size of the wall particles.

[0005] An object of the present invention is to provide a simulation method, a simulation device, and a program that can analyze the flow of a fluid that comes into contact with a wall surface without arranging wall particles to reproduce the wall surface. [Means for solving the problem]

[0006] According to one aspect of the present invention, The fluid flowing in contact with the wall is represented by multiple particles. determining particle-wall interactions between the plurality of particles and the wall and particle-particle interactions among the plurality of particles; a simulation method for evolving positions and velocities of the plurality of particles over time by solving an equation of motion that governs the motion of the plurality of particles, for each of the plurality of particles, the method comprising: When solving the equations of motion, A simulation method is provided in which, for particles among the plurality of particles whose distance to the wall surface is equal to or shorter than a first distance set as a simulation condition, the position and velocity of the particles are caused to evolve over time by applying, in addition to the forces due to the particle-particle interactions and the particle-wall interactions, a damping force received from the wall surface and a random force according to the temperature of the wall surface.

[0007] According to another aspect of the present invention, A simulation device for analyzing a fluid flow along a wall surface, comprising: an input section for inputting simulation conditions; a processing unit that analyzes the flow of the fluid based on the simulation conditions input to the input unit; an output unit that outputs the analysis result by the processing unit; Equipped with the processing unit represents the fluid with the plurality of particles based on the simulation conditions input to the input unit; evolving the positions and velocities of the particles over time by solving an equation of motion that governs the motion of the particles for each of the particles; When solving the equations of motion, A simulation device is provided in which, for particles among the plurality of particles whose distance to the wall surface is equal to or shorter than a first distance set in the simulation conditions, the forces due to the particle-particle interactions and particle-wall interactions set in the simulation conditions, the damping force received from the wall surface, and a random force according to the temperature of the wall surface are applied to cause the positions and velocities of the particles to evolve over time.

[0008] According to yet another aspect of the present invention, A program that causes a computer to execute a procedure for analyzing a fluid flow that flows along a wall surface, A procedure for obtaining simulation conditions; a step of analyzing the flow of the fluid based on the acquired simulation conditions; Execute The step of analyzing the fluid flow includes: expressing the fluid with the plurality of particles based on the acquired simulation conditions; a step of evolving the positions and velocities of the plurality of particles over time by solving an equation of motion that governs the motion of the plurality of particles for each of the plurality of particles; Including, When solving the equations of motion, A program is provided that causes the position and velocity of particles to evolve over time by applying forces due to particle-particle interactions and particle-wall interactions set in the simulation conditions, damping forces received from the wall surface, and random forces according to the temperature of the wall surface to particles among the plurality of particles whose distance to the wall surface is equal to or shorter than a first distance set in the simulation conditions. [Effects of the Invention]

[0009] Since wall particles that replicate wall surfaces are not placed, the increase in the number of particles to be calculated is suppressed, reducing the calculation load. By applying a random force to the fluid particles according to the damping force from the wall surface and the temperature of the wall surface, non-slip conditions on the wall surface are replicated, and the temperature of the fluid particles can be maintained at the temperature of the wall surface. [Brief explanation of the drawings]

[0010] [Figure 1] 1A and 1B are cross-sectional views schematically showing an example of an analytical model to be analyzed by a simulation method according to an embodiment. [Figure 2] Figure 2A is a schematic diagram showing the forces acting on a particle located at a distance from a wall surface farther than a first distance, and Figure 2B is a schematic diagram showing the forces acting on a particle located at a distance from a wall surface shorter than the first distance. [Figure 3]FIG. 3 is a block diagram of the simulation device according to this embodiment. [Figure 4] FIG. 4 is a flowchart showing the procedure of the simulation method according to this embodiment. [Figure 5] FIG. 5A is a schematic diagram of an analytical model in which the fluid to be analyzed is represented by multiple particles. FIG. 5B is a schematic diagram showing a particle system in which isotropic renormalization has been performed on the particle system shown in FIG. 5A. FIG. 5C is a schematic diagram showing a particle system in which renormalization has been performed on the particle system shown in FIG. 5B in the x and z directions without performing renormalization in the y direction. [Figure 6] FIG. 6 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model without renormalization. [Figure 7] FIG. 7 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed once in the x, y, and z directions. [Figure 8] FIG. 8 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which one renormalization was performed in the x-direction. [Figure 9] FIG. 9 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed once in the y-direction. [Figure 10] FIG. 10 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which one renormalization was performed in the z-direction. [Figure 11] FIG. 11 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed twice in the x, y, and z directions. [Figure 12] FIG. 12 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed twice in the x-direction. [Figure 13] FIG. 13 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed twice in the y-direction. [Figure 14]FIG. 14 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed twice in the z-direction. [Figure 15] FIG. 15 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed three times in the x, y, and z directions. [Figure 16] FIG. 16 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed three times in the x-direction. [Figure 17] FIG. 17 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed three times in the y-direction. [Figure 18] FIG. 18 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed three times in the z-direction. [Figure 19] FIG. 19 is a graph showing the distribution of the x-direction flow velocity in the y-direction calculated from the analysis results using an analytical model in which renormalization was performed three times in the x, y, and z directions. [Figure 20] FIG. 20 is a graph showing the distribution of flow velocity in the y direction in the x direction calculated from the results of analysis performed by the method according to the comparative example using an analytical model in which renormalization was performed three times in the x direction. [Figure 21] FIG. 21 is a graph showing the distribution of flow velocity in the y direction in the x direction calculated from the results of analysis performed by the method according to the comparative example using an analytical model in which renormalization was performed three times in the y direction. [Figure 22] FIG. 22 is a graph showing the distribution of flow velocity in the y direction in the x direction calculated from the results of analysis performed by the method according to the comparative example using an analytical model in which renormalization was performed three times in the z direction. [Figure 23] 23A and 23B are cross-sectional views of an analytical model that is the subject of analysis in a simulation method according to another embodiment. DETAILED DESCRIPTION OF THE INVENTION

[0011] A simulation method, a simulation device, and a program according to an embodiment will be described with reference to FIGS. 1A to 4. FIG.

[0012] 1A and 1B are cross-sectional views schematically illustrating an example of an analytical model to be analyzed by a simulation method according to an embodiment. FIG. 1A is a cross-sectional view taken along dashed dotted line 1A-1A in FIG. 1B, and FIG. 1B is a cross-sectional view taken along dashed dotted line 1B-1B in FIG. 1A. A fluid flows within an analytical space 30 defined by a wall surface 40. The wall surface 40 is composed of, for example, a pair of surfaces arranged parallel to each other. An xyz Cartesian coordinate system is defined in which a plane parallel to the wall surface 40 is defined as the xz plane and the direction of fluid flow is defined as the x direction.

[0013] The dimensions of the analysis space 30 between the wall surfaces 40 in the x, y, and z directions are denoted as Lx, Ly, and Lz, respectively. Periodic boundary conditions are applied to the boundaries perpendicular to the z direction and the boundaries perpendicular to the x direction. The fluid in the analysis space 30 is represented by a plurality of particles 31. In this embodiment, the behavior of the plurality of particles 31 is analyzed using molecular dynamics or renormalization group molecular dynamics. Specifically, the equations of motion governing the motion of the plurality of particles 31 are numerically solved to evolve the positions and velocities of the particles 31 over time.

[0014] As the interaction potential U acting between the particles 31, for example, the following Lennard-Jones potential is used.

number

[0015] FIG. 2A shows the distance L from the wall surface 40. iw 1 is a schematic diagram showing a force acting on a particle 31 located farther than a first distance L1. Here, the subscript i indicates the i-th particle 31 when the particles 31 are numbered consecutively.

[0016] Distance L from wall 40iw The following equation can be applied as the equation of motion that governs the motion of the particle 31 at a position farther than the first distance L1.

number

[0017] By representing the wall surface 40 with a plurality of wall particles and considering the interaction between the particles 31 and the wall particles, a non-slip condition where the flow velocity becomes zero at the wall surface is reproduced. In this example, the wall surface is not represented by wall particles, and the non-slip condition is reproduced.

[0018] FIG. 2B shows the distance L from the wall surface 40. iw 1 is a schematic diagram showing the force acting on a particle 31 (hereinafter, simply referred to as a particle 31 near the wall) located at a position equal to or less than a first distance L1. If a damping term is added to the equation of motion governing the particle 31 near the wall in order to reproduce the non-slip condition, the temperature of the particle 31 near the wall will decrease over time. In this example, the temperature of the particle 31 near the wall is controlled using the Langevin method.

[0019] When temperature control by the Langevin method is applied, the equation of motion that governs the motion of the particles 31 near the wall surface is expressed by the following formula.

number

[0020] R i is a random force applied to maintain the temperature of the particles 31 near the wall surface 40 at the wall surface temperature. i The magnitude of has a mean of zero and a standard deviation of σ s follows the normal distribution expressed by the following formula:

number

[0021] In temperature control by the Langevin method, the strength of control changes depending on the attenuation coefficient γ. If the attenuation coefficient γ is set to a small value, the temperature control becomes weaker, and the temperature tracking of the particles 31 becomes poor. Conversely, if the attenuation coefficient γ is set to a large value, the time step width Δt s In consideration of this, the set temperature T w Deviations may occur or the calculation may fail.

[0022] The attenuation coefficient γ is expressed by the following formula, for example.

number

number

[0023] Next, the normal force F acting on the i-th particle 31 from the wall surface 40 is iw The distance L from the particle 31 to the wall surface 40 will be explained. iw becomes equal to or less than the first distance L1, a virtual particle 31V is placed at a position symmetrical with respect to the wall surface 40. As the first distance L1, for example, 1 / 2 of the collision diameter σ in equation (1) is adopted. Note that in FIG. 2B, the diameter of the circle representing the particle 31 does not represent the collision diameter σ of the particle 31.

[0024] The i-th particle 31 near the wall is subjected to a force F based on the interparticle interaction potential defined by equation (1) from the virtual particle 31V. iw That is, when the particle 31 approaches the wall surface 40 to a distance equal to or less than half of the collision diameter σ, the particle 31 receives a repulsive force from the wall surface 40 in the normal direction.

[0025] Next, the simulation device according to this embodiment will be described with reference to Fig. 3. Fig. 3 is a block diagram of the simulation device according to this embodiment. The simulation device according to this embodiment includes an input unit 50, a processing unit 51, an output unit 52, and a storage unit 53. Simulation conditions and the like are input from the input unit 50 to the processing unit 51. Furthermore, various instructions (commands) and the like are input from an operator to the input unit 50. The input unit 50 is composed of, for example, a communication device, a removable media reader, a keyboard, and the like.

[0026] The processing unit 51 performs a simulation using molecular dynamics or renormalization group molecular dynamics (hereinafter simply referred to as molecular dynamics) based on the input simulation conditions and commands. Furthermore, the processing unit 51 outputs the simulation results to the output unit 52. The simulation results include, for example, information representing the state of particles 31 (FIGS. 1A and 1B) that reproduce the fluid that is the simulation object, and spatial and temporal changes in the physical quantities of the fluid. The processing unit 51 includes, for example, a central processing unit (CPU) of a computer. A program for causing the computer to execute a simulation using the molecular dynamics method is stored in the storage unit 53. The output unit 52 includes a communication device, a removable media writing device, a display, etc.

[0027] Next, the simulation method according to this embodiment will be described with reference to Fig. 4. Fig. 4 is a flowchart showing the procedure of the simulation method according to this embodiment. Each procedure shown in the flowchart is performed by the processing unit 51 executing a program stored in the storage unit 53.

[0028] First, the processing unit 51 acquires the simulation conditions input to the input unit 50 (step S1). The simulation conditions include information defining the fluid to be analyzed, information defining the shape of the wall surface 40 (FIGS. 1A and 1B), information for representing the fluid with a plurality of particles 31 (FIGS. 1A and 1B), information defining particle-wall interactions between the plurality of particles 31 and the wall surface 40, information defining particle-particle interactions between the plurality of particles 31, information defining temperature conditions, information defining the damping force that particles 31 near the wall surface receive from the wall surface 40, information defining external forces acting on the particles, a time step size, an analysis termination condition, etc.

[0029] The processing unit 51 places a plurality of particles 31 in the analysis space 30 (FIGS. 1A and 1B) based on the input simulation conditions. Thereafter, the procedures of steps S4, S5, and S6 are repeated for all particles 31 (step S3). Hereinafter, the repetitive processing of step S3 will be described, focusing on the i-th particle 31.

[0030] First, the distance L from the i-th particle 31 to the wall 40 iw It is determined whether the distance L from the i-th particle 31 to the wall surface 40 is equal to or less than the first distance L1 (step S4). iw is longer than the first distance L1, the force F due to the interparticle interaction shown in Figure 2A ij and external force F ext Taking into consideration the above, the equation of motion (Equation (2)) governing the motion of the particle 31 is numerically solved (step S5). iw When the distance is equal to or less than the first distance L1, the force F due to the interparticle interaction shown in Figure 2B ij , the force due to particle-wall interaction F iw , damping force -γ(v i -v w ), random force R according to the set temperature of wall surface 40 i , and external force F ext Taking this into consideration, the equation of motion (Equation (3)) that governs the motion of the particle 31 is numerically solved (step S6).

[0031] After solving the equations of motion for all particles 31, the positions and velocities of the particles 31 are evolved over time (step S7). The repetitive processing of step S3 and step S7 are repeated until the analysis is completed (step S8). The condition for completing the analysis is given, for example, by the simulation conditions acquired in step S1. When the analysis is completed, the analysis results are output to the output unit 52 (step S9).

[0032] Next, the excellent effects of the embodiment shown in Figures 1A to 4 will be described. In the above embodiment, the wall surface 40 (Figures 1A and 1B) is analyzed without being represented by multiple wall particles. This reduces the number of particles to be analyzed. As a result, the calculation load can be reduced.

[0033] In addition, the particle 31 near the wall is subjected to a damping force -γ(v i -v w) (step S6, FIG. 2B), it is possible to reproduce the non-slip condition on the surface of the wall surface 40. Furthermore, a random force R i By applying the force (step S6, FIG. 2B), the temperature of the particles 31 near the wall surface can be maintained at the set temperature.

[0034] In the above embodiment, the flow of a fluid flowing between two parallel wall surfaces is analyzed. However, even when the wall surfaces have other shapes, the simulation method and simulation device according to the above embodiment can be applied to analyze the flow of a fluid.

[0035] Next, a simulation method, a simulation device, and a program according to another embodiment will be described with reference to Figures 5A to 5C. In the embodiment described with reference to Figures 1A to 4, renormalization of particles 31 is not performed. In the embodiment described below, renormalization of particles 31 is performed to reduce the number of particles, and the equations of motion shown in formulas (3) and (4) are applied to the particle system after renormalization.

[0036] First, the renormalization technique applied in this embodiment will be described. FIG. 5A is a schematic diagram of an analytical model in which the fluid to be analyzed is represented by a plurality of particles 31. The fluid to be analyzed is contained in a space sandwiched between a pair of parallel wall surfaces 40. The fluid is represented by a plurality of particles 31. The shape of each particle 31 can be considered to correspond to an equipotential surface of the interaction potential generated by the particle 31. Generally, the shape of the equipotential surface of the interaction potential is spherical. In FIG. 5A, each particle 31 is represented by an equipotential surface of a certain size. An xyz Cartesian coordinate system is defined in which the direction in which the wall surfaces 40 are separated is the y direction. The dimension in the y direction of the space in which the fluid is contained is sufficiently smaller than the dimensions in the x and z directions.

[0037] FIG. 5B is a schematic diagram showing a particle system in which isotropic renormalization has been performed on the particle system shown in FIG. 5A. Renormalization reduces the number of particles. The interaction potential is isotropically stretched in the x, y, and z directions. For this reason, in FIG. 5B, each of the particles 32 after renormalization is represented by a spherical surface larger than the particle 31 shown in FIG. 5A.

[0038] FIG. 5C is a schematic diagram showing the particle system shown in FIG. 5B, where renormalization is not performed in the y direction, but is further performed in the x and z directions. Renormalization further reduces the number of particles. Furthermore, the interaction potential is elongated in the x and z directions, but not in the y direction. For this reason, in FIG. 5C, each particle 33 after renormalization is represented as a flattened spheroid with the y direction as its minor axis.

[0039] [Isotropic renormalization] Next, the renormalization transformation rule when performing isotropic renormalization on a particle system and the energy of the particle system will be described. The renormalization transformation rule is expressed by the following equation.

number

number

[0040] The interaction potential u(r) between particles is given, where r is the distance between particles. When u(r) approaches zero sufficiently quickly as the distance r approaches infinity (for example, in the case of the Lennard-Jones potential), the renormalization transformation law of the interaction potential u(r) is expressed by the following equation:

number

[0041] When renormalization processing is performed using the renormalization transformation rules of the above equations (7) and (9), the macroscopic physical quantities of the particle system, such as energy and pressure, remain unchanged before and after renormalization. Below, we will prove that energy remains unchanged before and after renormalization.

[0042] Partition function Z of a particle system N is expressed by the following formula:

number

[0043] The partition function Z in equation (10) N The part of the kinetic term Z N:k is expressed by the following formula:

number

[0044] Therefore, the kinetic energy E k is expressed by the following formula:

number

[0045] The partition function Z in equation (10) N The interaction term part Z N:int is expressed by the following formula:

number

number

[0046] Cutoff distance r c is the length L in the x, y, and z directions of the space containing the fluid. x , L y , L z Consider the case where the particle number density N / V is sufficiently small. In this case, the distance from each particle is the cutoff distance r c The probability that other particles exist within the range below is extremely small. In particular, the distance from each particle is the cutoff distance r c The probability that two or more other particles exist in the range below can be considered zero. Therefore, the partition function Z N:int The multiple integral (equation (13)) of can be approximated as follows:

[0047]

number

[0048] Therefore, the interaction energy E int is expressed as follows: approximation It is possible.

number

[0049] The total energy E of a particle system is defined by the following equation:

number

number

[0050] Interaction energy E after renormalization transformation int,R is approximated by the following formula:

number

[0051] In order for the approximation of Eq. (18) to hold, the cutoff distance r cR L x , L y , L z From equation (9), the cutoff distance r after the renormalization transformation is cR is the cutoff distance r before the renormalization transformation c When the renormalization factor λ is increased, the cutoff distance r cR Therefore, the upper limit of the renormalization factor λ is the cutoff distance r after the renormalization transformation. cR L x , L y , L z is constrained by the condition that it is sufficiently smaller than

[0052] When r is changed to λr in equation (19), equation (19) becomes the same as the right-hand side of equation (16). int,R =E int holds, and the interaction energy E int is also invariant before and after renormalization.

[0053] In this way, when isotropic renormalization is performed using the renormalization transformation rule shown in equation (7), the kinetic energy and interaction energy of the system remain unchanged before and after renormalization. Therefore, the energy of the entire system shown in equation (17) also remains unchanged before and after renormalization.

[0054] [Two-way renormalization] Next, the x-direction dimension of the flow field, L x and the dimension L in the z direction z Compared with the cutoff distance r c is sufficiently short, but the dimension L in the y direction y For example, this condition is met when the flow field has a thin plate-like shape with the y direction as the thickness direction.

[0055] When the shape of the flow field satisfies these conditions, renormalization is not performed in the y direction, but only in the x and z directions. The renormalization transformation law is expressed by the following equation.

number

[0056] The interaction potential follows the transformation rule:

number

[0057] Even when the renormalization transformation of equation (20) is performed, the kinetic energy expressed by equation (12) remains unchanged before and after the renormalization.

[0058] The volume integral in equation (16) can be approximated as follows:

number

[0059] Interaction energy E after renormalization transformation int,R is approximated by the following formula:

number

[0060] In this way, by performing anisotropic renormalization based on the renormalization transformation rule for two directions shown in equation (20), the energy of the particle system can be made unchanged before and after the renormalization.

[0061] [One-way renormalization] Next, the x-direction dimension of the flow field, L x Compared with the cutoff distance r c is sufficiently short, but the dimension L in the y direction y and the dimension L in the z direction z Compared with the cutoff distance r c For example, this condition is satisfied when the flow field has a long and thin cylindrical shape with the x direction as its length.

[0062] When the shape of the flow field satisfies these conditions, renormalization is not performed in the y and z directions, but only in the x direction. The renormalization transformation law is expressed by the following equation.

number

[0063] The interaction potential follows the transformation rule:

number

[0064] Even when the renormalization transformation of equation (24) is performed, the kinetic energy expressed by equation (12) remains unchanged before and after the renormalization.

[0065] The volume integral in equation (16) can be approximated as follows:

number

[0066] Interaction energy E after renormalization transformation int,R is approximated by the following formula:

number

[0067] In this way, by performing anisotropic renormalization based on the renormalization transformation rule for one direction shown in equation (24), the energy of the particle system can be made unchanged before and after the renormalization.

[0068] [Generalization of anisotropic renormalization] The renormalization transformation rule when renormalization is performed in two directions and not in the remaining one is shown in equation (20), and the renormalization transformation rule when renormalization is not performed in two directions and renormalization is performed in only the remaining one is shown in equation (24). Next, we will explain the case where the degree of renormalization is determined for each of the three directions and anisotropic renormalization is performed.

[0069] Renormalization factors in the x, y, and z directions , respectively. λx , λ y , and λ z The renormalization transformation law in this case is expressed as follows:

number

[0070] The interaction potential follows the transformation rule:

number

[0071] Interaction potential u after renormalization transformation R (r) Cutoff distance r in x, y, z directions cxR , r cyR , r czR is expressed by the following formula:

number

[0072] In this case, it is preferable to satisfy the following condition so that an approximation similar to that of equations (18), (22), and (26) holds.

number

[0073] Next, in order to apply the equations of motion shown in Eqs. (3) and (4) to the particle system after renormalization, we must define a transformation law for the attenuation coefficient γ. The attenuation coefficient after renormalization is defined as γ R When performing isotropic renormalization, applying the following transformation rule will result in similar results between the particle system before and after renormalization.

number

[0074] When performing anisotropic renormalization, it is possible to apply the following transformation rule from equation (32).

number

[0075] However, the attenuation coefficient γ obtained by using the transformation rule of Eq. (33) for the renormalized particle system R When a simulation was performed using the equation (2), similar results were not obtained between the particle system before and after renormalization. Through various analyses by the inventors of the present application, it was found that similar results can be obtained before and after renormalization by applying the following transformation rule when the fluid flow direction is the x direction.

number

[0076] Next, the excellent effects of this embodiment will be described. In this embodiment, by performing renormalization, the number of particles can be reduced and the calculation load can be reduced. Furthermore, by performing anisotropic renormalization according to the shape of the analysis space, the number of particles can be further reduced and the calculation load can be reduced.

[0077] Next, the results of actual simulations will be explained with reference to Figs. 6 to 22. The analysis model of Poiseuille flow shown in Figs. 1A and 1B was analyzed using a method without renormalization, a method with isotropic renormalization, and a method with anisotropic renormalization. Periodic boundary conditions were applied to the x and z directions. Dimension L x and in the z direction Dimension L z Both are set to 7.20 nm, and the y direction Dimension L y was set to 6.29 nm.

[0078] An analytical model was created to keep the dimensionless length the same even after renormalization. Water (H2O) was assumed as the fluid to be analyzed, and the fitting parameters ε and σ in equation (1) were set to 404.5K and 0.264nm, respectively. In addition, the external force F in equations (2) and (3) was set to extAn external force in the x direction was applied as follows. The magnitude of the external force was set so that the maximum flow velocity was approximately 100 m / s. The temperature of the wall surface 40 was set equal to the initial value of the temperature of the particle system.

[0079] Analysis was continued until the flow reached a steady state, and the flow velocity in the x direction was calculated from the moving speed of the particles 31 in the steady state. Figures 6 to 22 are graphs showing the distribution of the x direction flow velocity calculated from the analysis results in the y direction. The horizontal axis represents the position in the y direction, i.e., the distance from one wall surface 40, in units of [nm], and the vertical axis represents the x direction component of the flow velocity in units of [m / s]. The solid line in each graph indicates the theoretical value, and the circle symbol indicates the flow velocity calculated from the analysis results. The n attached to each graph x , n y , n z represent the renormalization times in the x, y, and z directions, respectively. That is, the renormalization times are defined by the following equations:

number

[0080] Renormalization count n x , n y , n z , or the renormalization factor λ x , λ y , λ z The renormalization conditions specified by the above are given as simulation conditions acquired in step S1 of FIG. 4, for example.

[0081] Figure 6 shows the analytical results without renormalization. The damping coefficient γ is set to 5.23×10 -13The average error of the analytical results relative to the theoretical value was approximately 3.8%, indicating that the analytical results agree well with the theoretical value. Figure 7 shows the analytical results when renormalization was performed isotropically once each. Figures 8, 9, and 10 show the analytical results when renormalization was performed once in the x-direction, y-direction, and z-direction, respectively, and no renormalization was performed in the other directions. Note that the renormalization transformation rule of Equation (34) was used in Figures 7 to 10. The average errors of the analytical results shown in Figures 7, 8, 9, and 10 relative to the theoretical value were approximately 2.5%, 2.2%, 7.5%, and 3.2%, respectively, indicating that the analytical results agree well with the theoretical value.

[0082] Figure 11 shows the analytical results when renormalization was performed twice isotropically. Figures 12, 13, and 14 show the analytical results when renormalization was performed twice in the x, y, and z directions, respectively, and no renormalization was performed in other directions. Note that the renormalization transformation rule of equation (34) is used in Figures 11 to 14. The average errors of the analytical results shown in Figures 11, 12, 13, and 14 relative to the theoretical values ​​are approximately 1.9%, 4.2%, 8.4%, and 2.7%, respectively, and it can be seen that the analytical results agree well with the theoretical values.

[0083] Figure 15 shows the analytical results when renormalization was performed isotropically three times. Figures 16, 17, and 18 show the analytical results when renormalization was performed three times in the x, y, and z directions, respectively, and no renormalization was performed in other directions. Note that the transformation rule of equation (34) is used in Figures 15 to 18. The average errors of the analytical results shown in Figures 15, 16, 17, and 18 relative to the theoretical values ​​are approximately 2.7%, approximately 4.4%, approximately 10.0%, and approximately 7.2%, respectively, and it can be seen that the analytical results agree well with the theoretical values.

[0084] 19 to 22 are graphs showing the results of analysis using the renormalization transformation law shown in Equation (33). Figure 19 shows the analysis results when renormalization was performed isotropically three times. Figures 20, 21, and 22 show the analysis results when renormalization was performed three times in the x, y, and z directions, respectively, and no renormalization was performed in other directions. The average errors of the analysis results shown in Figures 19, 20, 21, and 22 relative to the theoretical values ​​are approximately 4.1%, 33.6%, 86.9%, and 71.1%, respectively. Note that although the renormalization conditions are the same in Figures 15 and 19, differences in the analysis results arise due to factors such as the initial arrangement of particles 31 and the random numbers used to generate the random force Ri in Equation (3). It is expected that the difference between the two will become smaller if the time ensemble is performed over a long period of time.

[0085] When renormalization is performed isotropically, the analytical results agree well with the theoretical values ​​even when the transformation rule of Equation (33) is used. However, when renormalization is performed anisotropically, when the transformation rule of Equation (33) is applied, the analytical results deviate significantly from the theoretical values. It can be seen that the transformation rule of Equation (33) cannot be applied when renormalization is performed anisotropically.

[0086] The analysis results shown in Fig. 6 confirm that, when renormalization is not performed, the flow of a fluid in contact with the wall surface 40 can be analyzed with high accuracy by the simulation method according to the embodiment described with reference to Figs. 1A to 4. The analysis results shown in Figs. 7 to 18 confirm that, when renormalization is performed isotropically or anisotropically, the flow of a fluid in contact with the wall surface 40 can be analyzed with high accuracy by using the transformation rule of Equation (34). The analysis results shown in Fig. 19 confirm that, when renormalization is performed isotropically, the flow of a fluid in contact with the wall surface 40 can be analyzed with high accuracy by using the transformation rule of Equation (33).

[0087] Next, a simulation method and a simulation device according to still another embodiment will be described with reference to FIGS. 23A and 23B.

[0088] 5, the renormalization method was explained using an example of a fluid flowing in one direction between two parallel wall surfaces 40. The transformation law shown in equation (34) can also be applied to a rotating flow.

[0089] 23A and 23B are cross-sectional views of an analytical model that is the subject of analysis by the simulation method according to this embodiment. Fig. 23A is a cross-sectional view taken along dashed line 23A-23A in Fig. 23B, and Fig. 23B is a cross-sectional view taken along dashed line 23B-23B in Fig. 23A.

[0090] Two cylindrical wall surfaces 40A and 40B with different diameters are arranged concentrically. A flow path through which the fluid circulates is formed between the wall surfaces 40A and 40B. An xyz Cartesian coordinate system is defined, with the direction parallel to the central axis of the cylindrical wall surfaces 40A and 40B being the y direction. Since the x and z directions are equivalent, the number of renormalizations in the x direction and the z direction are set to be the same when performing renormalization. In this case, the renormalization factor λ in the x direction is x and the renormalization factor λ in the z direction z is equal to the renormalization factor in the x and z directions, λ xz Then, equation (34) is transformed as follows:

number

[0091] Rotating fluids can be analyzed by using the transformation law of equation (36). For example, this transformation law can be applied to the analysis of fluid flow in a stirred tank, where the outer wall surface 40B is the casing and the inner wall surface 40A is the stirring blade.

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

[0093] 30 Analysis space 31 particles 31V Virtual Particle 32, 33 Renormalized particles 40, 40A, 40B Wall 50 Input section 51 Processing section 52 Output section 53 Memory section

Claims

1. The fluid flowing in contact with the wall is represented by multiple particles. determining particle-wall interactions between the plurality of particles and the wall and particle-particle interactions among the plurality of particles; a simulation method for evolving positions and velocities of the plurality of particles over time by solving an equation of motion that governs the motion of the plurality of particles, for each of the plurality of particles, the method comprising: When solving the equations of motion, A simulation method in which, for particles among the plurality of particles whose distance to the wall surface is equal to or shorter than a first distance set as a simulation condition, the position and velocity of the particles are caused to evolve over time by applying, in addition to the forces due to the particle-particle interactions and the particle-wall interactions, a damping force received from the wall surface and a random force according to the temperature of the wall surface.

2. The fluid flows in one direction, When an xyz Cartesian coordinate system is defined in which the direction of the flow of the fluid is the x direction, the plurality of particles are renormalized in at least one direction of the x direction, the y direction, and the z direction, The number of renormalizations in the x, y, and z directions is n ｘ , n ｙ , n ｚ and marked it as Renormalization factor λ, which represents the degree of renormalization ｘ , λ ｙ , λ ｚ of, [Equation 1] and marked it as The damping coefficient before the renormalization of the damping force is γ, and the damping coefficient after the renormalization is γ Ｒ When denoted as [Equation 2] Applying the damping coefficient γ after renormalization Ｒ Calculate When solving the equation of motion, the damping coefficient γ after renormalization Ｒ The simulation method according to claim 1, wherein the following formula is used:

3. A simulation device for analyzing a fluid flow along a wall surface, comprising: an input section for inputting simulation conditions; a processing unit that analyzes the flow of the fluid based on the simulation conditions input to the input unit; an output unit that outputs the analysis result by the processing unit; Equipped with the processing unit represents the fluid with a plurality of particles based on the simulation conditions input to the input unit; evolving the positions and velocities of the particles over time by solving an equation of motion that governs the motion of the particles for each of the particles; When solving the equations of motion, A simulation device that applies forces due to particle-particle interactions and particle-wall interactions set in the simulation conditions, damping forces received from the wall surface, and random forces according to the temperature of the wall surface to particles among the plurality of particles whose distance to the wall surface is equal to or shorter than a first distance set in the simulation conditions, thereby evolving the positions and velocities of the particles over time.

4. the simulation conditions include renormalization conditions for renormalizing the plurality of particles, The fluid flows in one direction, When an xyz orthogonal coordinate system is defined in which the direction of the fluid flow is the x direction, the renormalization condition is the number of renormalizations in the x direction, y direction, and z direction, n ｘ , n ｙ , n ｚ containing information specifying Renormalization factor λ, which represents the degree of renormalization ｘ , λ ｙ , λ ｚ of, [Equation 3] and marked it as The damping coefficient before the renormalization of the damping force is γ, and the damping coefficient after the renormalization is γ Ｒ When the transformation rule is expressed as [Equation 4] Applying the damping coefficient γ after renormalization Ｒ Calculate When solving the equation of motion, the damping coefficient γ after renormalization Ｒ 4. The simulation device according to claim 3, wherein:

5. A program that causes a computer to execute a procedure for analyzing a fluid flow that flows along a wall surface, A procedure for obtaining simulation conditions; a step of analyzing the flow of the fluid based on the acquired simulation conditions; Execute The step of analyzing the fluid flow includes: expressing the fluid with a plurality of particles based on the acquired simulation conditions; a step of evolving the positions and velocities of the plurality of particles over time by solving an equation of motion that governs the motion of the plurality of particles for each of the plurality of particles; Including, When solving the equations of motion, A program that applies forces due to particle-particle interactions and particle-wall interactions set in the simulation conditions, damping forces received from the wall surface, and random forces according to the temperature of the wall surface to the particles among the plurality of particles whose distance to the wall surface is equal to or shorter than a first distance set in the simulation conditions, thereby evolving the positions and velocities of the particles over time.

6. the simulation conditions include renormalization conditions for renormalizing the plurality of particles, The fluid flows in one direction, When an xyz orthogonal coordinate system is defined in which the direction of the fluid flow is the x direction, the renormalization condition is the number of renormalizations in the x direction, y direction, and z direction, n ｘ , n ｙ , n ｚ containing information specifying Renormalization factor λ, which represents the degree of renormalization ｘ , λ ｙ , λ ｚ of, [Equation 5] and marked it as The damping coefficient before the renormalization of the damping force is γ, and the damping coefficient after the renormalization is γ Ｒ When marked as The step of analyzing the fluid flow includes: Transformation Law [Equation 6] Applying the damping coefficient γ after renormalization Ｒ Calculating When solving the equation of motion, the damping coefficient γ after renormalization Ｒ The program according to claim 5, wherein the program is

Images