SIMULATION METHOD, SIMULATION DEVICE, AND NON-TRANSITORY COMPUTER READABLE MEDIUM STORING PROGRAM

A simulation method in which a fluid flowing in contact with a wall surface is represented by a plurality of particles, particle-wall surface interaction and interparticle interaction are determined, and an equation of motion governing motion of the plurality of particles is solved for each of the plurality of particles to develop positions and velocities of the plurality of particles over time includes causing, in a case where the equation of motion is solved, attenuation force received from the wall surface and random force according to a temperature of the wall surface, in addition to force due to the interparticle interaction and the particle-wall surface interaction, to act on a particle, among the plurality of particles, whose distance to the wall surface is equal to or less than a first distance set in a simulation condition to develop a position and a velocity of a particle over time.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority to Japanese Patent Application No. 2022-087797, filed on May 30, 2022, which is incorporated by reference herein in its entirety.

BACKGROUND

Technical Field

A certain embodiment of the present invention relates to a simulation method, a simulation device, and a non-transitory computer readable medium storing a program.

Description of Related Art

In a method in the related art of analyzing a flow of a fluid in contact with a wall surface by using a molecular dynamics method, the wall surface is represented by a plurality of wall particles, and an interaction between the wall particle and a fluid particle is considered (refer to the related art). With control of a temperature of the wall particle, a temperature of the fluid particle near the wall surface is reproduced.

SUMMARY

According to one aspect of the present invention, there is provided a simulation method in which a fluid flowing in contact with a wall surface is represented by a plurality of particles, particle-wall surface interaction between the plurality of particles and the wall surface and interparticle interaction between the plurality of particles are determined, and an equation of motion governing motion of the plurality of particles is solved for each of the plurality of particles to develop positions and velocities of the plurality of particles over time, which includes causing, in a case where the equation of motion is solved, attenuation force received from the wall surface and random force according to a temperature of the wall surface, in addition to force due to the interparticle interaction and the particle-wall surface interaction, to act on a particle, among the plurality of particles, whose distance to the wall surface is equal to or less than a first distance set in a simulation condition to develop a position and a velocity of a particle over time.

According to another aspect of the present invention, there is provided a simulation device that analyzes a flow of a fluid flowing along a wall surface, which includes an input unit that receives a simulation condition, a processing unit that analyzes the flow of the fluid based on the simulation condition input to the input unit, and an output unit that outputs an analysis result obtained by the processing unit, in which the processing unit represents the fluid with a plurality of particles based on the simulation condition input to the input unit, solves an equation of motion governing motion of the plurality of particles for each of the plurality of particles to develop positions and velocities of the plurality of particles over time, and causes, in a case where the equation of motion is solved, force due to interparticle interaction and particle-wall surface interaction set in the simulation condition, attenuation force received from the wall surface, and random force according to a temperature of the wall surface, to act on a particle, among the plurality of particles, whose distance to the wall surface is equal to or less than a first distance set in the simulation condition to develop a position and a velocity of a particle over time.

According to still another aspect of the present invention, there is provided a non-transitory computer readable medium storing a program that causes a computer to execute a procedure of analyzing a flow of a fluid flowing along a wall surface, which includes a procedure of acquiring a simulation condition, and a procedure of analyzing the flow of the fluid based on the acquired simulation condition, in which the procedure of analyzing the flow of the fluid includes a procedure of representing the fluid with a plurality of particles based on the acquired simulation condition, and a procedure of solving an equation of motion governing motion of the plurality of particles for each of the plurality of particles to develop positions and velocities of the plurality of particles over time, and in a case where the equation of motion is solved, force due to interparticle interaction and particle-wall surface interaction set in the simulation condition, attenuation force received from the wall surface, and random force according to a temperature of the wall surface are caused to act on a particle, among the plurality of particles, whose distance to the wall surface is equal to or less than a first distance set in the simulation condition to develop a position and a velocity of a particle over time.

DETAILED DESCRIPTION

In a simulation method in the related art, since it is necessary to dispose the wall particle in addition to the fluid particle, the number of particles to be calculated increases. Further, a spatial resolution for a wall surface shape is limited by a particle size of the wall particle.

It is desirable to provide a simulation method, a simulation device, and a non-transitory computer readable medium storing a program capable of analyzing a flow of a fluid in contact with a wall surface without disposing a wall particle for reproducing the wall surface.

Since the wall surface particle that reproduces the wall surface is not disposed, it is possible to suppress an increase in the number of particles to be calculated and reduce a calculation load. With the attenuation force received from the wall surface and the random force according to the temperature of the wall surface being caused to act on the particle of the fluid, it is possible to reproduce the non-slip condition on the wall surface and maintain the temperature of the fluid particle at the temperature of the wall surface.

A simulation method, a simulation device, and a non-transitory computer readable medium storing a program according to one embodiment will be described with reference to drawings fromFIGS.1A to4.

FIGS.1A and1Bare cross-sectional views schematically showing an example of an analysis model to be analyzed in the simulation method according to one embodiment.FIG.1Ais a cross-sectional view taken along a line1A-1A of a single point chain line ofFIG.1B, andFIG.1Bis a cross-sectional view taken along a line1B-1B of a single point chain line ofFIG.1A. A fluid flows in an analysis space30defined by a wall surface40. The wall surface40is configured of, for example, a pair of surfaces disposed parallel to each other. An xyz orthogonal coordinate system is defined in which a plane parallel to the wall surface40is set as an xz plane and a flow direction of the fluid is set as an x direction.

Sizes of the analysis space30between the wall surfaces40in x, y, and z directions are respectively marked as Lx, Ly, and Lz. A periodic boundary condition is applied to a boundary perpendicular to the z direction and a boundary perpendicular to the x direction. The fluid in the analysis space30is represented by a plurality of particles31. In the embodiment, behavior of the plurality of particles31is analyzed by using a molecular dynamics method or a renormalization group molecular dynamics method. Specifically, an equation of motion governing motion of the plurality of particles31is numerically solved, and a position and a velocity of the particle31are developed over time.

For example, the following Leonard-Jones potential is used as interaction potential U acting between the particles31.

Here, r is a distance between particles, and ε and σ are fitting parameters having dimensions of energy and length, respectively. σ may be referred to as a collision diameter of the particles.

FIG.2Ais a schematic diagram showing force acting on the particle31at a position where a distance Liwfrom the wall surface40is farther than a first distance L1. Here, the subscript i means an i-th particle31in a case where the plurality of particles31are serially numbered.

The following equation can be applied as the equation of motion governing the motion of the particle31at a position where the distance Liwfrom the wall surface40is farther than the first distance L1.

Here, m is mass of the particle31, and viis a velocity vector of the i-th particle31. Fijis force that the i-th particle31receives from a j-th particle31, and can be calculated from the Leonard-Jones potential shown in equation (1). Fextis external force such as gravity that the particle31receives. A symbol on the right side of equation (2) means a sum of the particles31around the i-th particle31.

With representation of the wall surface40with a plurality of wall particles and consideration of interaction between the particle31and the wall particle, a non-slip condition in which a flow velocity becomes zero on the wall surface is reproduced. In the present embodiment, the non-slip condition is reproduced without representing the wall surface with the wall particle.

FIG.2Bis a schematic diagram showing force acting on the particle31(hereinafter may be simply referred to as particle31near the wall surface) at a position where the distance Liwfrom the wall surface40is equal to or less than the first distance L1. In a case where an attenuation term is added to the equation of motion governing the particle31near the wall surface in order to reproduce the non-slip condition, a temperature of the particle31near the wall surface decreases with the passage of time. In the present embodiment, the temperature of the particle31near the wall surface is controlled by using a Langevin method.

In a case where the temperature control by the Langevin method is applied, the equation of motion governing the motion of the particle31near the wall surface is represented by the following equation.

Here, Fiwis force in a normal direction that the i-th particle31receives from the wall surface40. viand vware the velocity vector of the i-th particle31and the velocity vector of the wall surface40to which the i-th particle31approaches, respectively. A fourth term on the right side of equation (3) corresponds to attenuation force (viscous resistance force). The magnitude of the attenuation force is proportional to a relative velocity of the particle31to the wall surface40, and a direction of the attenuation force is opposite to the relative velocity. γ is referred to as an attenuation coefficient. With this attenuation force, the non-slip condition on the wall surface40can be reproduced.

Riis random force applied to maintain the temperature of the particle31near the wall surface40at a wall surface temperature. The mean of the magnitude of the random force Riis zero and the standard deviation σsthereof follows normal distribution represented by the following equation.

Here, kBis a Boltzmann constant, Twis a set temperature of the wall surface40, and Δtsis a time step width in a case where the equation of motion is numerically solved. A direction of the random force Riis random.

In the temperature control by the Langevin method, the strength of the control changes depending on the attenuation coefficient γ. In a case where the attenuation coefficient γ is set to be small, the temperature control is weakened and temperature followability of the particle31deteriorates. On the contrary, in a case where the attenuation coefficient γ is set to be large, a deviation from the set temperature Twmay occur or the calculation may fail due to a problem of balance with the time step width Δts.

The attenuation coefficient γ is represented by, for example, the following equation.

Here, ζ is an attenuation ratio, and k is a spring constant. For example, 0.707 can be used as the attenuation ratio ζ. The spring constant k can be obtained from a coefficient of a quadratic term in a case where the Leonard-Jones potential (equation (1)) is subjected to Taylor expansion near a saddle point, for example. In this case, the spring constant k can be calculated by the following equation.

For example, 5.23×10−13kg/s can be used as the attenuation coefficient γ.

Next, force Fiwin the normal direction that the i-th particle31receives from the wall surface40will be described. In a case where the distance Liwfrom the particle31to the wall surface40is equal to or less than the first distance L1, a virtual particle31V is disposed at a plane-symmetrical position with respect to the wall surface40. For example, ½ of the collision diameter σ of equation (1) is employed as the first distance L1. InFIG.2B, a circular diameter representing the particle31does not represent the collision diameter σ of the particle31.

The force Fiw, from the virtual particle31V, based on interparticle interaction potential defined by equation (1) is caused to act on the i-th particle31near the wall surface. That is, in a case where the particle31approach the wall surface40to a distance that is ½ or less of the collision diameter σ, the particle31receives repulsive force in the normal direction from the wall surface40.

Next, the simulation device according to the present embodiment will be described with reference toFIG.3.FIG.3is a block diagram of the simulation device according to the present embodiment. The simulation device according to the embodiment includes an input unit50, a processing unit51, an output unit52, and a storage unit53. A simulation condition and the like are input from the input unit50to the processing unit51. Further, various commands and the like are input from an operator to the input unit50. The input unit50is configured of, for example, a communication device, a removable media reading device, or a keyboard.

The processing unit51performs the simulation using the molecular dynamics method or the renormalization group molecular dynamics method (hereinafter simply referred to as molecular dynamics method) based on the input simulation condition and command. Further, a simulation result is output to the output unit52. The simulation result includes, for example, information representing a state of the particle31(FIGS.1A and1B) reproducing the fluid, which is a simulation object, a spatial change and temporal change in a physical quantity of the fluid, and the like. The processing unit51includes, for example, a central processing unit (CPU) of a computer. A program for causing the computer to execute the simulation by the molecular dynamics method is stored in the storage unit53. The output unit52includes a communication device, a removable media writing device, a display, or the like.

Next, the simulation method according to the present embodiment will be described with reference toFIG.4.FIG.4is a flowchart showing a procedure of the simulation method according to the present embodiment. Each procedure shown in the flowchart is performed by the processing unit51executing the program stored in the storage unit53.

First, the processing unit51acquires the simulation condition input to the input unit50(step S1). The simulation condition includes information that defines the fluid to be analyzed, information that defines a shape of the wall surface40(FIGS.1A and1B), information that represents the fluid with the plurality of particles31(FIGS.1A and1B), information that defines particle-wall surface interaction between the plurality of particles31and the wall surface40, information that defines interparticle interaction between the plurality of particles31, information that defines a temperature condition, information that defines the attenuation force that the particle31near the wall surface receives from the wall surface40, information that defines external force acting on the particle, the time step width, an analysis end condition, and the like.

The processing unit51disposes the plurality of particles31in the analysis space30(FIGS.1A and1B) based on the input simulation condition. After that, procedures of step S4, step S5, and step S6are repeated for all the particles31(step S3). Hereinafter, the repetition processing in step S3will be described focusing on the i-th particle31.

First, determination is made whether or not the distance Liwfrom the i-th particle31to the wall surface40is equal to or less than the first distance L1(step S4). In a case where the distance Liwfrom the i-th particle31to the wall surface40is longer than the first distance L1, the equation of motion (equation (2)) governing the motion of the particle31is numerically solved in consideration of the force Fijdue to the interparticle interaction shown inFIG.2Aand the external force Fext(step S5). In a case where the distance Liwfrom the i-th particle31to the wall surface40is equal to or less than the first distance L1, the equation of motion (equation (3)) governing the motion of the particle31is numerically solved in consideration of the force Fijdue to the interparticle interaction shown inFIG.2B, the force Fiwdue to the particle-wall surface interaction, attenuation force −γ(vi−vw) according to the relative velocity between the particle31and the wall surface, the random force Riaccording to the set temperature of the wall surface40, and the external force Fext(step S6).

The equation of motion is solved for all the particles31and then the position and velocity of each of the particles31are developed over time (step S7). The repetition processing in step S3and step S7are repeated until the analysis ends (step S8). The analysis end condition is provided, for example, by the simulation condition acquired in step S1. In a case where the analysis ends, an analysis result is output to the output unit52(step S9).

Next, the excellent effects of the embodiments shown inFIGS.1A to4will be described. In the above embodiment, the analysis is performed without representing the wall surface40(FIGS.1A and1B) with the plurality of wall particles. Therefore, the number of particles to be analyzed is reduced. As a result, a calculation load can be reduced.

The non-slip condition on a surface of the wall surface40can be reproduced by causing the attenuation force −γ(vi−vw) from the wall surface40to act on the particles31near the wall surface (step S6,FIG.2B). Further, the temperature of the particles31near the wall surface can be maintained at the set temperature by causing the random force Riaccording to the temperature of the wall surface40to act on the particles31near the wall surface (step S6,FIG.2B).

In the above embodiment, the flow of the fluid flowing through the two wall surfaces disposed in parallel is analyzed. However, even in a case where the wall surface has another shape, the flow of the fluid can be analyzed by applying the simulation method and the simulation device according to the above embodiment.

Next, a simulation method, a simulation device, and a non-transitory computer readable medium storing a program according to another embodiment will be described with reference toFIGS.5A to5C. In the embodiment described with reference toFIGS.1A to4, the particle31is not renormalized. In the embodiment described below, the particle31is renormalized to reduce the number of particles, and the equations of motion shown in equations (3) and (4) are applied to a particle system after the renormalization.

First, a renormalization method applied in the present embodiment will be described.

FIG.5Ais a schematic diagram of an analysis model in which the fluid to be analyzed is represented by the plurality of particles31. The fluid to be analyzed is accommodated in a space sandwiched between a pair of wall surfaces40disposed in parallel. The fluid is represented by the plurality of particles31. Each shape of the particle31can be considered to correspond to an equipotential surface of the interaction potential generated by the particle31. In general, the equipotential surface of the interaction potential has a spherical shape. InFIG.5A, each of the particles31is represented by the equipotential surface of a certain size. An xyz orthogonal coordinate system is defined in which a direction in which the wall surfaces40are separated is set as the y direction. A size in the y direction of the space in which the fluid is accommodated is sufficiently smaller than sizes in the x direction and the z direction.

FIG.5Bis a schematic diagram showing a particle system in which the particle system shown inFIG.5Ais isotropically renormalized. The number of particles is reduced by the renormalization. The interaction potential is stretched isotropically in the x, y, and z directions. For this reason, inFIG.5B, each of particles32after the renormalization is represented by a spherical surface larger than the particle31shown inFIG.5A.

FIG.5Cis a schematic diagram showing a particle system in which the particle system shown inFIG.5Bis not renormalized in the y direction and is further renormalized in the x direction and the z direction. The number of particles is further reduced by the renormalization. Further, the interaction potential is stretched in the x direction and the z direction, and is not stretched in the y direction. For this reason, inFIG.5C, each of particles33after the renormalization is represented by a flat rotating ellipse with the y direction as a short axis direction.

Next, a renormalization transformation rule in a case where the isotropic renormalization is performed on the particle system and energy of the particle system will be described. The renormalization transformation rule is represented by the following equation.

Here, N is the number of particles, m is particle mass, r is a position vector indicating a particle position, V is a volume of a flow field of the fluid, and T is a temperature of the particle system. A parameter with a subscript R is represented as a parameter after the renormalization. λ is a parameter (renormalization factor) representing a degree of renormalization, and λ is a real number larger than 1. For example, in a case where the renormalization factor λ is represented by the following equation with n as an integer of 1 or more, n is referred to as the number of times of renormalization.

The interaction potential u(r) between the particles is provided with a distance between the particles as r. In a case where the distance r approaches infinity and u(r) approaches zero sufficiently quickly (for example, in case of Leonard-Jones potential), the renormalization transformation rule of the interaction potential u(r) is represented by the following equation.

In a case where the renormalization processing is performed by using the renormalization transformation rule of the above equations (7) and (9), macroscopic physical quantities such as energy and pressure of the particle system do not change before and after the renormalization. Hereinafter, it will be demonstrated that the energy does not change before and after the renormalization.

A distribution function ZNof the particle system is represented by the following equation.

Here, h is the Planck's constant, kBis the Boltzmann's constant, r is the position vector of the particle, p is momentum of the particle, and rijis a distance from an i-th particle to a j-th particle. The sigma in a first term on the right side of a third equation in equation (10) means that all of N particles are summed, and the sigma in a second term means that all of particle pairs are summed.

A portion ZN:kof a motion term of the distribution function ZNof equation (10) is represented by the following equation.

Therefore, kinetic energy Ekis represented by the following equation.

In a case where the renormalization of the renormalization factor λ is performed, both N and β on the right side of equation (12) become 1/λ3. Therefore, the kinetic energy Ekdoes not change before and after the renormalization.

A portion ZN:intof an interaction term of the distribution function ZNof equation (10) is represented by the following equation.

In a case where the distance rijapproaches infinity, u(rij) approaches zero sufficiently quickly. Therefore, the interaction potential u(r) in a case of r>rccan be approximated as follows, using a certain cutoff distance rc.

A case is considered in which the cutoff distance rcis sufficiently smaller than lengths Lx, Ly, and Lzof the space accommodating the fluid in the x, y, and z directions and particle number density N/V is sufficiently small. In this case, the probability that another particle is present in a range where the distance from each particle is equal to or less than the cutoff distance rcis extremely small. In particular, the probability that two or more other particles are present in the range where the distance from each particle is equal to or less than the cutoff distance rcmay be regarded as zero. Therefore, the multiple integral (equation (13)) of the distribution function ZN:intcan be approximated as follows.

Therefore, interaction energy Eintcan be approximated by the following equation.

The total energy E of the particle system is defined by the following equation.

In a case where the cutoff distance rcis sufficiently smaller than Lx, Ly, and Lz, and the interparticle distance r is larger than the cutoff distance rc, u(r)=0 can be approximated. Therefore, the volume integral of equation (16) can be approximated as follows.

Interaction energy Eint,Rafter the renormalization transformation is approximated by the following equation.

There is a need to satisfy a condition that a cutoff distance rcRafter the renormalization transformation is sufficiently smaller than Lx, Ly, and Lzin order to establish the approximation of equation (18). From equation (9), the cutoff distance rcRafter the renormalization transformation is A times the cutoff distance rcbefore the renormalization transformation. In a case where the renormalization factor λ is increased, the cutoff distance rcRbecomes longer. For this reason, an upper limit of the renormalization factor λ is constrained by the condition that the cutoff distance rcRafter the renormalization transformation is sufficiently smaller than Lx, Ly, and Lz.

In a case where r is variable-converted to λr in equation (19), equation (19) becomes the same as the right side of equation (16). Therefore, Eint,R=Eintis established, and the interaction energy Eintis also unchanged before and after the renormalization.

In this manner, in a case where the isotropic renormalization is performed using the renormalization transformation rule shown in equation (7), the kinetic energy and the interaction energy of the system do not change before and after the renormalization. Therefore, the energy of the entire system shown in equation (17) is also unchanged before and after the renormalization.

Renormalization in Two Directions

Next, a case will be described in which the cutoff distance rcis sufficiently shorter than the size Lxin the x direction and the size Lzin the z direction of the flow field, but is not sufficiently shorter than the size Lyin the y direction. For example, in a case where the flow field has a thin plate shape with the y direction as a thickness direction, the above condition is satisfied.

In a case where the shape of the flow field satisfies such a condition, the renormalization is not performed in the y direction, but is performed in only two directions of the x direction and the z direction. The renormalization transformation rule is represented by the following equation.

The interaction potential follows the following transformation rule.

Even in a case where the renormalization transformation of equation (20) is performed, the kinetic energy represented by equation (12) does not change before and after the renormalization.

The volume integral of equation (16) can be approximated as follows.

The interaction energy Eint,Rafter the renormalization transformation is approximated by the following equation.

Here, r-tilde is the same as r-tilde of equation (21). In a case where x is variable-converted to λx and z is variable-converted to λz in equation (23), the right side of equation (23) becomes the same as the right side of equation (16). Therefore, Eint,R=Eintis established, and the interaction energy Eintis also unchanged before and after the renormalization.

In this manner, with the anisotropic renormalization based on the renormalization transformation rule in the two directions shown in equation (20), the energy of the particle system can be unchanged before and after the renormalization.

Renormalization in One Direction

Next, a case will be described in which the cutoff distance rcis sufficiently shorter than the size Lxin the x direction of the flow field, but is not sufficiently shorter than the size Lyin the y direction and the size Lzin the z direction. For example, in a case where the flow field has an elongated cylindrical shape with the x direction as a length direction, the above condition is satisfied.

In a case where the shape of the flow field satisfies such a condition, the renormalization is not performed in the y direction and the z direction, but is performed only in the x direction. The renormalization transformation rule is represented by the following equation.

The interaction potential follows the following transformation rule.

Even in a case where the renormalization transformation of equation (24) is performed, the kinetic energy represented by equation (12) does not change before and after the renormalization.

The volume integral of equation (16) can be approximated as follows.

The interaction energy Eint,Rafter the renormalization transformation is approximated by the following equation.

Here, r-tilde is the same as the r-tilde of equation (25). In a case where x is variable-converted to λx in equation (27), the right side of equation (27) becomes the same as the right side of equation (16). Therefore, Eint,R=Eintis established, and the interaction energy Eintis also unchanged before and after the renormalization.

In this manner, with the anisotropic renormalization based on the renormalization transformation rule in one direction shown in equation (24), the energy of the particle system can be unchanged before and after the renormalization.

Generalization of Anisotropic Renormalization

Equation (20) shows the renormalization transformation rule in a case where the renormalization is performed in two directions and is not performed in a remaining one direction. Equation (24) shows the renormalization transformation rule in a case where the renormalization is not performed in two directions and is performed only in the remaining one direction. Next, a case will be described in which the degree of the renormalization is determined for each of the three directions and the anisotropic renormalization is performed.

The renormalization factors in the x direction, y direction, and z direction are respectively marked as λx, λy, and λz. The renormalization transformation rule in this case is represented by the following equation.

The interaction potential follows the following transformation rule.

Cutoff distances rcxR, rcyR, and rczRof the interaction potential uR(r) after the renormalization transformation in the x, y, and z directions are represented by the following equations.

In this case, it is preferable to satisfy the following condition in order to establish an approximation similar to the approximations of equations (18), (22), and (26).

Next, in order to apply the equations of motion shown in Equations (3) and (4) to the particle system after the renormalization, a transformation rule for the attenuation coefficient γ is required to be defined. The attenuation coefficient after the renormalization is marked as γR. In a case where the isotropic renormalization is performed, similar results can be obtained between the particle system before the renormalization and the particle system after the renormalization, in a case where the following transformation rule is applied.

In a case where the anisotropic renormalization is performed, the following transformation rule is considered to be applied, from equation (32).

However, in a case where the simulation is performed on the particle system subjected to the renormalization using an attenuation coefficient γRobtained by using the transformation rule of equation (33), similar results cannot be obtained between the particle systems before and after the renormalization. With various analyzes by the inventors of the present application, it has been found that similar results can be obtained before and after the renormalization in a case where the following transformation rule is applied with the flow direction of the fluid of the x direction.

Next, the excellent effects of the present embodiment will be described.

In the present embodiment, with the renormalization, it is possible to reduce the number of particles and thus the calculation load. Further, with the anisotropic renormalization according to the shape of the analysis space, it is possible to further reduce the number of particles and thus the calculation load.

Next, results of actual simulation will be described with reference toFIGS.6to22. The analysis model of the Poiseuille flow shown inFIGS.1A and1Bis analyzed by a method of not performing the renormalization, a method of performing the isotropic renormalization, and a method of performing the anisotropic renormalization. The periodic boundary condition is applied to the x direction and the z direction. Both the size Lxin the x direction and the size Lzin the z direction are set to 7.20 nm, and the size Lyin the y direction is set to 6.29 nm.

An analysis model is created such that dimensionless lengths are maintained the same even after the renormalization. Water (H2O) is assumed as the fluid to be analyzed, the fitting parameter c of equation (1) is set to 404.5K, and σ is set to 0.264 nm. Further, external force in the x direction is caused to act as the external force Fextof the equations (2) and (3). The magnitude of the external force is set such that a maximum value of the flow velocity is approximately 100 m/s. The temperature of the wall surface40is made equal to an initial value of the temperature of the particle system.

The analysis is performed until the flow reaches a steady state, and the flow velocity in the x direction is calculated from a movement speed of the particle31in the steady state.FIGS.6to22are graphs showing distribution of the flow velocity in the x direction calculated from analysis results in the y direction. The horizontal axis represents a position in the y direction, that is, a distance from one wall surface40in unit [nm], and the vertical axis represents an x direction component of the flow velocity in unit [m/s]. A solid line of each graph indicates a theoretical value, and a circle symbol indicates the flow velocity obtained from the analysis result. The nx, ny, and nzassigned to each graph represent the number of times of renormalization in the x direction, the y direction, and the z direction, respectively. That is, the number of times of renormalization is defined by the following equation.

The renormalization condition designated by the number of times of renormalization nx, ny, and nz, the renormalization factor λx, λy, and λz, or the like is provided, for example, by the simulation condition acquired in step S1ofFIG.4.

FIG.6shows an analysis result in a case where the renormalization is not performed. The attenuation coefficient γ is set to 5.23×10−13kg/s. An average error of the analysis result with respect to the theoretical value is about 3.8%, and it can be seen that the analysis result matches well the theoretical value.FIG.7shows an analysis result in a case where the renormalization is performed isotropically once each.FIG.8,FIG.9, andFIG.10show analysis results in a case where the renormalization is performed once in the x direction, the y direction, and the z direction, respectively, and the renormalization is not performed in the other directions. InFIGS.7to10, the renormalization transformation rule of equation (34) is used. The average errors of the analysis results shown inFIG.7,FIG.8,FIG.9, andFIG.10with respect to the theoretical values are respectively approximately 2.5%, approximately 2.2%, approximately 7.5%, and approximately 3.2%, and it can be seen that the analysis results match well the theoretical values.

FIG.11shows an analysis result in a case where the renormalization is performed isotropically twice each.FIG.12,FIG.13, andFIG.14show analysis results in a case where the renormalization is performed twice in the x direction, the y direction, and the z direction, respectively, and the renormalization is not performed in the other directions. InFIGS.11to14, the renormalization transformation rule of equation (34) is used. The average errors of the analysis results shown inFIGS.11,12,13, and14with respect to the theoretical values are respectively approximately 1.9%, approximately 4.2%, approximately 8.4%, and approximately 2.7%, and it can be seen that the analysis results match well the theoretical values.

FIG.15shows an analysis result in a case where the renormalization is performed isotropically three times each.FIG.16,FIG.17, andFIG.18show analysis results in a case where the renormalization is performed three times in the x direction, the y direction, and the z direction, respectively, and the renormalization is not performed in the other directions. InFIGS.15to18, the transformation rule of equation (34) is used. The average errors of the analysis results shown inFIGS.15,16,17, and18with respect to the theoretical values are respectively approximately 2.7%, approximately 4.4%, approximately 10.0%, and approximately 7.2%, and it can be seen that the analysis results match well the theoretical values.

FIGS.19to22are graphs showing results of the analysis performed by using the renormalization transformation rule shown in equation (33).FIG.19shows an analysis result in a case where the renormalization is performed isotropically three times each.FIG.20,FIG.21, andFIG.22show analysis results in a case where the renormalization is performed three times in the x direction, the y direction, and the z direction, respectively, and the renormalization is not performed in the other directions. The average errors of the analysis results shown inFIGS.19,20,21, and22with respect to the theoretical values are respectively approximately 4.1%, approximately 33.6%, approximately 86.9%, and approximately 71.1%. Although the renormalization conditions are the same inFIGS.15and19, there is a difference in the analysis results due to the influence of initial disposition of the particles31, a difference in random numbers in a case where the random force Riin equation (3) is generated, or the like. The difference therein is considered to be smaller in a case where the ensemble is taken for a long time.

It can be seen that in a case where the renormalization is performed isotropically, the analysis result matches well the theoretical value even in a case where the transformation rule of equation (33) is used. However, in a case where the anisotropic renormalization is performed, the analysis result deviates significantly from the theoretical value in a case where the transformation rule of equation (33) is applied. It can be seen that in a case where the anisotropic renormalization is performed, the transformation rule of equation (33) cannot be applied.

From the analysis result shown inFIG.6, it is confirmed that in a case where the renormalization is not performed, the flow of the fluid in contact with the wall surface40can be accurately analyzed by the simulation method according to the embodiment described with reference toFIGS.1A to4. From the analysis results shown inFIGS.7to18, it is confirmed that in a case where the isotropic or anisotropic renormalization is performed, the flow of the fluid in contact with the wall surface40can be accurately analyzed by using the transformation rule of equation (34). From the analysis result shown inFIG.19, it is confirmed that in a case where the isotropic renormalization is performed, the flow of the fluid in contact with the wall surface40can be accurately analyzed by using the transformation rule of equation (33).

Next, a simulation method and a simulation device according to still another embodiment will be described with reference toFIGS.23A and23B.

In the embodiment shown inFIG.5, the renormalization method has been described by taking the fluid flowing between two parallel wall surfaces40in one direction as an example. The transformation rule shown in equation (34) is also applicable to a rotating flow.

FIG.23AandFIG.23Bare cross-sectional views of an analysis model to be analyzed in the simulation method according to the present embodiment.FIG.23Ais a cross-sectional view taken along a line23A-23A of a single point chain line ofFIG.23B, andFIG.23Bis a cross-sectional view taken along a line23B-23B of a single point chain line ofFIG.23A.

Two cylindrical wall surfaces40A and40B having different diameters are concentrically disposed. A flow path around which the fluid revolves is formed between the wall surfaces40A and40B. An xyz orthogonal coordinate system is defined in which a direction parallel to a center axis of the cylindrical wall surfaces40A and40B is set as the y direction. Since the x direction and the z direction are equivalent to each other, the number of times of renormalization in the x direction and the number of times of renormalization in the z direction are set to be the same in a case where the renormalization is performed. In this case, the renormalization factor λxin the x direction and the renormalization factor λzin the z direction are equal to each other. In a case where the renormalization factor in the x direction and the z direction is marked as λxz, equation (34) is modified as follows.

With the use of the transformation rule of equation (36), it is possible to analyze the rotating fluid. The above transformation rule can be applied to, for example, an analysis of a flow of a fluid in a stirring tank, with the outer wall surface40B as a casing and the inner wall surface40A as a stirring blade.

It is needless to say that each of the above embodiments is an example, and partial substitutions or combinations of the configurations shown in different embodiments are possible. Similar action effects due to the configuration similar to the plurality of embodiments will not be sequentially referred to for each embodiment. Further, the present invention is not limited to the above embodiment. For example, it will be obvious to those skilled in the art that various changes, improvements, combinations, and the like are possible.