JP2017224197A - Molecular dynamics-based simulation method, program, and simulation device - Google Patents
Molecular dynamics-based simulation method, program, and simulation device Download PDFInfo
- Publication number
- JP2017224197A JP2017224197A JP2016119871A JP2016119871A JP2017224197A JP 2017224197 A JP2017224197 A JP 2017224197A JP 2016119871 A JP2016119871 A JP 2016119871A JP 2016119871 A JP2016119871 A JP 2016119871A JP 2017224197 A JP2017224197 A JP 2017224197A
- Authority
- JP
- Japan
- Prior art keywords
- particle
- velocity
- simulation
- temperature
- particles
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 69
- 238000004088 simulation Methods 0.000 title claims abstract description 60
- 238000000329 molecular dynamics simulation Methods 0.000 title claims abstract description 19
- 239000002245 particle Substances 0.000 claims abstract description 152
- 230000008569 process Effects 0.000 claims description 5
- 239000006185 dispersion Substances 0.000 claims description 2
- 238000012935 Averaging Methods 0.000 claims 1
- 239000012530 fluid Substances 0.000 description 13
- 238000010586 diagram Methods 0.000 description 12
- 238000004364 calculation method Methods 0.000 description 7
- 230000004048 modification Effects 0.000 description 7
- 238000012986 modification Methods 0.000 description 7
- 230000006399 behavior Effects 0.000 description 5
- 230000000694 effects Effects 0.000 description 4
- 230000008859 change Effects 0.000 description 2
- 230000003993 interaction Effects 0.000 description 2
- 239000000203 mixture Substances 0.000 description 1
- 230000000737 periodic effect Effects 0.000 description 1
- 238000004904 shortening Methods 0.000 description 1
- 230000000087 stabilizing effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
- G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
- G16C20/30—Prediction of properties of chemical compounds, compositions or mixtures
Landscapes
- Theoretical Computer Science (AREA)
- Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Life Sciences & Earth Sciences (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Health & Medical Sciences (AREA)
- General Health & Medical Sciences (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Physics & Mathematics (AREA)
- Chemical & Material Sciences (AREA)
- Crystallography & Structural Chemistry (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
Description
本発明は、分子動力学法によるシミュレーション方法、プログラム、及びシミュレーション装置に関する。 The present invention relates to a simulation method, a program, and a simulation apparatus based on a molecular dynamics method.
分子動力学を用いたシミュレーションが行われている。分子動力学では、シミュレーション対象の系を構成する粒子の運動方程式が数値的に解かれる。分子動力学を用いたシミュレーションを行う場合、シミュレーション対象の系の温度をある目標温度に設定したい場合がある。 Simulations using molecular dynamics have been carried out. In molecular dynamics, the equations of motion of the particles that make up the system to be simulated are numerically solved. When performing a simulation using molecular dynamics, it may be desired to set the temperature of the system to be simulated to a certain target temperature.
下記の特許文献１に、分子動力学によるシミュレーションにおいて、温度設定を変更した場合に、解析状態が迅速に安定化し、所望の設定温度での状態を正確に得ることができるようにする温度制御方法が開示されている。この方法では、分子動力学法によるシミュレーションの途中段階で求められた各粒子の速度を、現時点の系の温度及び目標温度に基づいて更新する。更新された速度に基づいて、系の温度を更新する。この操作により、系の温度を目標温度に近づける。 Patent Document 1 listed below discloses a temperature control method for quickly stabilizing an analysis state and accurately obtaining a state at a desired set temperature when the temperature setting is changed in a molecular dynamics simulation. Is disclosed. In this method, the velocity of each particle obtained in the middle of the simulation by the molecular dynamics method is updated based on the current system temperature and target temperature. Update the system temperature based on the updated speed. By this operation, the temperature of the system is brought close to the target temperature.
シミュレーション対象の系に流れ場が存在する場合、従来の温度制御方法では、流れ場の流速も温度制御にともなって変化してしまう。このため、流れ場が存在する系のシミュレーション時に従来の温度制御方法を適用すると、シミュレーション結果はシミュレーション対象の系の各粒子の挙動を反映したものにならない。従って、従来の温度制御方法は、流れ場が存在する系のシミュレーションに適用することができない。 When a flow field exists in the system to be simulated, the flow rate of the flow field also changes with temperature control in the conventional temperature control method. For this reason, if a conventional temperature control method is applied during simulation of a system in which a flow field exists, the simulation result does not reflect the behavior of each particle in the system to be simulated. Therefore, the conventional temperature control method cannot be applied to simulation of a system in which a flow field exists.
本発明の目的は、流れ場が存在する系の分子動力学法によるシミュレーションにおいて、温度制御を行うことが可能なシミュレーション方法を提供することである。 An object of the present invention is to provide a simulation method capable of performing temperature control in a simulation by a molecular dynamics method of a system in which a flow field exists.
本発明の一観点によると、
流れ場を形成する複数の粒子を含む粒子系を分子動力学法によりシミュレーションする方法であって、
前記粒子系の粒子の各々の位置及び速度を算出する際に、粒子ごとに、粒子の速度と、粒子の位置における流れ場の流速との差に基づいて温度制御を行うことにより、前記粒子系の温度を目標温度に近づけるシミュレーション方法が提供される。
According to one aspect of the invention,
A method of simulating a particle system including a plurality of particles forming a flow field by a molecular dynamics method,
When calculating the position and velocity of each particle of the particle system, the temperature of the particle system is controlled for each particle based on the difference between the velocity of the particle and the flow velocity of the flow field at the particle position. A simulation method is provided to bring the temperature of the target close to the target temperature.
本発明の他の観点によると、上記シミュレーション方法をコンピュータに実行させるプログラム、及び上記シミュレーション方法を実行するシミュレーション装置が提供される。 According to another aspect of the present invention, a program for causing a computer to execute the simulation method and a simulation apparatus for executing the simulation method are provided.
粒子の速度と、粒子の位置における流れ場の流速との差に基づいて温度制御を行うことにより、流れ場の流速が温度制御の影響を受けにくくなる。これにより、シミュレーション対象の粒子系の流れ場の挙動がシミュレーション結果に反映されることとなる。 By performing the temperature control based on the difference between the velocity of the particle and the flow field velocity at the particle position, the flow field velocity is less affected by the temperature control. Thereby, the behavior of the flow field of the particle system to be simulated is reflected in the simulation result.
図１、図２Ａ〜図２Ｂを参照して、実施例によるシミュレーション方法について説明する。実施例によるシミュレーションにおいて分子動力学法が適用される。分子動力学法では、シミュレーションの対象となる系を複数の粒子で表し、各粒子の運動方程式を数値的に解くことにより粒子の挙動、例えば位置の変化及び速度が求まる。実施例においては、シミュレーション対象である粒子系の複数の粒子によって流れ場が形成される。 A simulation method according to the embodiment will be described with reference to FIGS. 1 and 2A to 2B. The molecular dynamics method is applied in the simulation according to the embodiment. In the molecular dynamics method, a system to be simulated is represented by a plurality of particles, and the behavior of the particles, for example, change in position and velocity, can be obtained by numerically solving the equation of motion of each particle. In the embodiment, a flow field is formed by a plurality of particles of a particle system to be simulated.
図１に、実施例によるシミュレーション方法のフローチャートを示す。以下に説明するシミュレーションは、コンピュータがシミュレーションプログラムを実行することにより実現される。 FIG. 1 shows a flowchart of a simulation method according to the embodiment. The simulation described below is realized by a computer executing a simulation program.
ステップ１０において、シミュレーション対象の粒子系を定義する物理量、各粒子の初期状態、拘束条件、及び温度制御の実行間隔を設定する。例えば、オペレータの操作によって初期状態、拘束条件、及び温度制御の実行間隔が入力される。粒子系を定義する物理量には、例えば粒子の質量、粒子のポテンシャル等が含まれる。粒子の初期状態には、例えば粒子の位置、粒子の速度、粒子系の温度等が含まれる。粒子系の温度は、粒子系に含まれる複数の粒子の速度の分散に反映される。拘束条件には、境界条件、流れ場の流速、粒子系の目標温度等が含まれる。流れ場の流速は、粒子系に含まれる粒子の速度の平均値に反映される。分子動力学法を用いた動解析の複数のタイムステップのうち、温度制御の実行間隔の設定値で示されたタイムステップ数ごとに温度制御が実行される。 In step 10, a physical quantity that defines a particle system to be simulated, an initial state of each particle, a constraint condition, and a temperature control execution interval are set. For example, an initial state, a constraint condition, and a temperature control execution interval are input by an operator's operation. The physical quantity defining the particle system includes, for example, particle mass, particle potential, and the like. The initial state of the particle includes, for example, the position of the particle, the velocity of the particle, the temperature of the particle system, and the like. The temperature of the particle system is reflected in the dispersion of the velocities of a plurality of particles contained in the particle system. The constraint conditions include boundary conditions, flow field velocity, particle system target temperature, and the like. The flow velocity of the flow field is reflected in the average value of the velocity of the particles contained in the particle system. Of the plurality of time steps of the dynamic analysis using the molecular dynamics method, the temperature control is executed for each number of time steps indicated by the set value of the temperature control execution interval.
次に、ステップ１１において、現時点で実行するタイムステップが、温度制御を実行するタイムステップに該当するか否かを判定する。例えば、温度制御の実行間隔の設定値で示された回数のタイムステップごとに、温度制御が実行される。 Next, in step 11, it is determined whether or not the time step executed at the present time corresponds to a time step for executing temperature control. For example, the temperature control is executed every time step indicated by the set value of the temperature control execution interval.
現時点で実行するタイムステップが温度制御を実行するタイムステップに該当しない場合、ステップ１２において、通常の１タイムステップの動解析を実行する。現時点で実行するタイムステップが温度制御を実行するタイムステップに該当する場合、ステップ１３において、１タイムステップの動解析及び温度制御を実行する。 When the time step executed at the present time does not correspond to the time step for executing the temperature control, in Step 12, a normal dynamic analysis of one time step is executed. When the time step executed at the present time corresponds to the time step for executing the temperature control, the dynamic analysis and the temperature control of one time step are executed in step 13.
ステップ１２及びステップ１３の後、ステップ１４においてシミュレーションを終了するか否かを判定する。シミュレーションを終了しない場合、ステップ１１に戻ってステップ１１からステップ１４までの処理を繰り返す。シミュレーションを終了する場合には、ステップ１５においてシミュレーション結果を出力する。 After step 12 and step 13, it is determined in step 14 whether or not to end the simulation. If the simulation is not terminated, the process returns to step 11 and the processes from step 11 to step 14 are repeated. When ending the simulation, a simulation result is output in step 15.
次に、ステップ１３の処理について詳細に説明する。
まず、分子動力学法を用いて、１タイムステップの動解析を実行する。１タイムステップの時間刻み幅は、予め設定されている。この動解析により、シミュレーション対象の粒子系の各粒子の位置及び速度が更新される。
Next, the process of step 13 will be described in detail.
First, a dynamic analysis of one time step is executed using the molecular dynamics method. The time increment of one time step is set in advance. By this dynamic analysis, the position and velocity of each particle of the particle system to be simulated are updated.
１タイムステップの動解析を実行した後、すべての粒子について各粒子の位置における流れ場の流速を算出する。流れ場の流速は、着目する粒子、及び着目する粒子の周囲に存在する他の複数の粒子の速度の平均値で表すことができる。流速を算出する基礎となる複数の粒子として、例えば着目する粒子を含む一部の空間（以下、流速算出単位空間という。）の内側に存在する粒子を採用することができる。流速算出単位空間として、例えば、着目する粒子を中心として、平衡状態のときの粒子間の平均距離の３倍〜４倍程度の半径を持つ球体とすることができる。 After performing the dynamic analysis of one time step, the flow velocity of the flow field at the position of each particle is calculated for all particles. The flow velocity of the flow field can be represented by an average value of the velocities of the target particle and other particles existing around the target particle. As a plurality of particles serving as a basis for calculating the flow velocity, for example, particles existing inside a part of a space including the focused particle (hereinafter referred to as a flow velocity calculation unit space) can be employed. As the flow velocity calculation unit space, for example, a sphere having a radius of about 3 to 4 times the average distance between particles in the equilibrium state with the focused particle at the center can be used.
次に、図２Ａを参照して、流れ場の流速を算出する具体例について説明する。
図２Ａに、シミュレーション対象の粒子系の各粒子の位置及び速度を表した模式図を示す。図２Ａにおいて、各粒子を丸記号で表し、各粒子の速度ベクトルを矢印で表す。各粒子の速度ベクトルは、流れ場の下流側の方（図２Ａにおいて右方）を向く傾向を示す。
Next, a specific example of calculating the flow velocity of the flow field will be described with reference to FIG. 2A.
FIG. 2A is a schematic diagram showing the position and speed of each particle in the particle system to be simulated. In FIG. 2A, each particle is represented by a circle symbol, and the velocity vector of each particle is represented by an arrow. The velocity vector of each particle shows a tendency toward the downstream side of the flow field (right side in FIG. 2A).
着目する粒子Ｐｉを含む流速算出単位空間３０を定義することができる。着目する粒子Ｐｉの速度をｖｉと表し、流速算出単位空間３０内の他の粒子Ｐｊの速度をｖｊと表す。流速算出単位空間３０内の着目する粒子Ｐｉ以外の粒子の個数をＮＮと表す。着目する粒子Ｐｉの位置における流れ場の流速Ｖｉは、以下の式で表される。
流速Ｖｉを算出した後、各粒子Ｐｉ、Ｐｊの速度ｖｉ、ｖｊと、着目する粒子Ｐｉの位置における流れ場の流速Ｖｉとの差に基づいて、着目する粒子Ｐｉの位置における温度を算出する。 After calculating the flow velocity V i, each particle P i, P j velocity v i of, v j and, based on the difference between the velocity V i of the flow field at the position of the focused particle P i, the focused particle P i The temperature at the position is calculated.
図２Ｂに、シミュレーション対象の粒子系の各粒子Ｐｉ、Ｐｊの位置、及び速度ｖｉ、ｖｊと流速Ｖｉとの差を表した模式図を示す。各粒子に付された矢印が、その粒子の速度ｖｉ、ｖｊと流速Ｖｉとの差を表している。粒子の速度ｖｉ、ｖｊと流速Ｖｉとの差を流速算出単位空間３０内の粒子について平均した速度は０になる。 FIG. 2B is a schematic diagram showing the positions of the particles P i and P j of the particle system to be simulated and the difference between the velocities v i and v j and the flow velocity V i . The arrow attached to each particle represents the difference between the velocity v i , v j of the particle and the flow velocity V i . The average velocity of the particles in the flow velocity calculation unit space 30 is 0 as the difference between the particle velocity v i , v j and the flow velocity V i .
着目する粒子Ｐｉの位置における温度をＴｉ、ボルツマン定数をｋＢ、粒子Ｐｉ、Ｐｊの質量をそれぞれｍｉ、ｍｊと表すと、以下の式（２）が成り立つ。式（２）を用いて、着目する粒子Ｐｉの位置における温度Ｔｉを算出することができる。
式（２）の右辺は、着目している粒子Ｐｉの速度ｖｉと流速Ｖｉとの差に基づく運動エネルギ、及び着目している粒子Ｐｉの周囲に存在する他の複数の粒子Ｐｊの速度ｖｊと流速Ｖｉとの差に基づく運動エネルギの平均値を表している。この平均値に基づいて、着目する粒子Ｐｉの位置における温度Ｔｉが算出される。 Right side of the equation (2) is the kinetic energy based on the difference between the speed v i and the flow velocity V i of the particles P i of interest, and other plurality of particles P existing around the grains P i of interest It represents the average value of kinetic energy based on the difference between the velocity v j of j and the flow velocity V i . Based on this average value, the temperature Ti at the position of the particle Pi of interest is calculated.
温度Ｔｉを算出した後、各粒子の位置における温度が目標温度に近づくように、各粒子の速度を調整する。速度を調整する一手法として、速度スケーリング法を適用することができる。速度スケーリング法においては、着目する粒子Ｐｉの速度調整後の速度をｖｉｍと表すと、速度ｖｉｍは、以下の式で表される。
ここで、ａはスケーリングの強さを決めるパラメータであり、１以下の正の実数である。ｓｉは、スケーリングファクタと呼ばれる。Ｔｃは目標温度である。
After calculating the temperature T i, the temperature at the position of each particle so as to approach the target temperature, adjust the speed of each particle. As a method for adjusting the speed, a speed scaling method can be applied. In the speed scaling method, when the speed after speed adjustment of the particle P i of interest is expressed as v im , the speed v im is expressed by the following expression.
Here, a is a parameter that determines the strength of scaling, and is a positive real number of 1 or less. s i is called the scaling factor. T c is the target temperature.
粒子の速度を調整することにより、１タイムステップの動解析及び温度制御の処理が終了する。 By adjusting the velocity of the particles, the dynamic analysis and temperature control process in one time step is completed.
次に、分子動力学法によるシミュレーションにおいて、流れ場の流速を考慮しないで温度制御を適用する方法と比較して、上記実施例の効果について説明する。流れ場の流速を考慮しないで温度制御を行う場合には、例えば、着目する粒子Ｐｉの位置における温度Ｔｉは上述の式（２）に代えて、以下の式（４）を用いて算出される。
さらに、温度制御による速度調整は、式（３）に代えて、以下の式（５）が適用される。
上述の式（４）の速度ｖｉ、ｖｊには、図２Ａに示したように流れ場の流速Ｖｉが含まれている。流速Ｖｉを含んだ速度ｖｉ、ｖｊに基づいて温度制御が行われるため、温度制御によって流れ場の流速Ｖｉが乱されてしまう。上記実施例による方法では、式（３）に示したように、粒子の速度ｖｉ、ｖｊと、粒子Ｐｉの位置における流れ場の流速Ｖｉとの差に基づいて温度制御が行われ、温度制御実施後の速度に流速Ｖｉが加算される。このため、温度制御による流速Ｖｉの乱れを抑制することができる。 The speeds v i and v j in the above equation (4) include the flow field velocity V i as shown in FIG. 2A. Since the flow velocity V i laden velocity v i, the temperature control on the basis of v j is performed, resulting in disturbed flow velocity V i of the flow field by the temperature control. In the method according to the above embodiment, as shown in the equation (3), the temperature control is performed based on the difference between the particle velocities v i and v j and the flow velocity V i of the flow field at the position of the particle P i. The flow velocity V i is added to the speed after the temperature control is performed. Therefore, it is possible to suppress the disturbance of the flow velocity V i by the temperature control.
上記実施例の効果を確認するために、同一のシミュレーションモデルを用いて、温度制御を行わない方法、流速を考慮しないで温度制御を行う方法、及び実施例による方法の３通りの方法でシミュレーションを行った。次に、これらのシミュレーション結果について説明する。 In order to confirm the effect of the above-described embodiment, the simulation is performed by using the same simulation model in three ways: a method in which temperature control is not performed, a method in which temperature control is performed without considering the flow velocity, and a method according to the embodiment. went. Next, these simulation results will be described.
図３にシミュレーションモデルを示す。シミュレーションは２次元空間で行った。２次元の正方形の領域内に円形のシリンダが配置されている。この空間に流れる流体を構成する粒子の挙動をシミュレーションにより求めた。ｘｙ直交座標系を定義し、正方形の領域の１つの頂点を原点とし、正方形の各辺を、ｘ軸またはｙ軸に平行にした。正方形の一辺の長さＬ１を３９２２×１０−１０ｍとし、シリンダの直径Ｄを５００×１０−１０ｍとした。シリンダは、ｙ方向に関して正方形の領域の中心に配置されている。正方形の領域のｘ方向の負側の辺（ｘ＝０）からシリンダの中心までの距離Ｌ２をＬ１／４＝９８０．５×１０−１０ｍとした。 FIG. 3 shows a simulation model. The simulation was performed in a two-dimensional space. A circular cylinder is arranged in a two-dimensional square area. The behavior of the particles constituting the fluid flowing in this space was obtained by simulation. An xy orthogonal coordinate system was defined, with one vertex of the square area as the origin, and each side of the square parallel to the x-axis or y-axis. The length L1 of one side of the square was 3922 × 10 −10 m, and the cylinder diameter D was 500 × 10 −10 m. The cylinder is arranged at the center of a square area with respect to the y direction. The distance L2 from the negative side (x = 0) in the x direction of the square region to the center of the cylinder was L1 / 4 = 980.5 × 10 −10 m.
流体を構成する粒子の質量を６．６３×１０−２６ｋｇとした。正方形の領域内に配置される粒子の個数は約２６万個である。粒子間の相互作用ポテンシャルとして、レナードジョーンズポテンシャルの斥力の項のみを用いた。具体的には、相互作用ポテンシャルφ（ｒ）は、以下の式で定義される。
ここで、ε＝１．６５×１０−２１Ｊ、σ＝３．４１×１０−１０ｍである。この粒子系の流体の粘度は３．６４×１０−４Ｐａ・ｓであり、密度は１．１８×１０３ｋｇ／ｍ３である。
The mass of the particles constituting the fluid was 6.63 × 10 −26 kg. The number of particles arranged in the square area is about 260,000. Only the repulsive force term of the Leonard Jones potential was used as the interaction potential between particles. Specifically, the interaction potential φ (r) is defined by the following equation.
Here, ε = 1.65 × 10 −21 J and σ = 3.41 × 10 −10 m. This particle-based fluid has a viscosity of 3.64 × 10 −4 Pa · s and a density of 1.18 × 10 3 kg / m 3 .
さらに、シリンダは、円周上に相互にバネで接続された複数の粒子を配置することにより表現した。シリンダ表面に配置した粒子に対しては、ランジュバン法による温度制御を行った。シリンダの温度は８８．８５Ｋとした。 Furthermore, the cylinder was expressed by arranging a plurality of particles connected to each other by a spring on the circumference. The particles placed on the cylinder surface were temperature controlled by the Langevin method. The cylinder temperature was 88.85K.
温度制御を行う場合の目標温度Ｔｃとして、温度制御を行わない方法でシミュレーションしたときの定常状態到達時の平均温度である２３５．１５Ｋを用いた。 As the target temperature Tc in the case of performing temperature control, 235.15 K that is an average temperature when reaching a steady state when a simulation is performed by a method that does not perform temperature control is used.
次に、境界条件について説明する。右端（ｘ＝Ｌ１）と左端（ｘ＝０）とで特殊な周期境界条件が課されている。具体的には、粒子が右端の境界を通過するとき、ｘ方向の流速６００ｍ／ｓを中心とし、温度を８８．８５Ｋとするボルツマン分布に基づいて、ランダムに新しい速度が割り当てられる。この粒子が、左端の境界から正方形の領域内に流入する。下端（ｙ＝０）と上端（ｙ＝Ｌ１）の境界は、ｘ方向の流速が６００ｍ／ｓ、温度が８８．８５Ｋの熱浴として働く。下端及び上端の境界に到達した粒子は、ｘ方向の流速６００ｍ／ｓだけシフトしたハーフマクスウェル分布に基づいて、ランダムに新しい速度が割り当てられる。 Next, boundary conditions will be described. Special periodic boundary conditions are imposed at the right end (x = L1) and the left end (x = 0). Specifically, when the particle passes through the right end boundary, a new velocity is randomly assigned based on a Boltzmann distribution centered at a flow velocity of 600 m / s in the x direction and a temperature of 88.85K. The particles flow into the square area from the left edge boundary. The boundary between the lower end (y = 0) and the upper end (y = L1) acts as a heat bath having a flow velocity in the x direction of 600 m / s and a temperature of 88.85K. Particles that reach the lower and upper boundary are randomly assigned new velocities based on a half-Maxwell distribution shifted by a flow velocity in the x direction of 600 m / s.
図４Ａ、図４Ｂ、図５Ａ、図５Ｂ、図６Ａ、及び図６Ｂにシミュレーショ結果を示す。図４Ａ及び図４Ｂは、温度制御を行わない方法を用いてシミュレーションを行った結果を示し、図５Ａ及び図５Ｂは、流速を考慮しないで温度制御を行う方法を用いてシミュレーションを行った結果を示し、図６Ａ及び図６Ｂは、実施例による方法でシミュレーションを行った結果を示す。 4A, 4B, 5A, 5B, 6A, and 6B show the simulation results. 4A and 4B show the results of simulation using a method that does not perform temperature control, and FIGS. 5A and 5B show the results of simulation using a method that performs temperature control without considering flow velocity. FIG. 6A and FIG. 6B show the results of simulation by the method according to the embodiment.
図４Ａ、図５Ａ、及び図６Ａは、流体に生じた渦を可視化した図である。渦を可視化する方法は、例えば特許第５８８８９２１号公報に開示されている。図４Ｂ、図５Ｂ、及び図６Ｂは、温度分布を濃淡で表す図である。色の相対的に濃い領域が相対的に高温であることを表している。図４Ａ〜図６Ｂまでの各々において、左から右に並ぶ複数の図は、時間の経過に対応している。 4A, FIG. 5A, and FIG. 6A are views visualizing vortices generated in a fluid. A method for visualizing vortices is disclosed in, for example, Japanese Patent No. 5888921. 4B, FIG. 5B, and FIG. 6B are diagrams showing the temperature distribution in shades. A relatively dark area represents a relatively high temperature. In each of FIGS. 4A to 6B, a plurality of diagrams arranged from left to right correspond to the passage of time.
図４Ａに示されているように、温度制御を行わない場合には、シリンダの後方にカルマン渦が発生していることを確認できる。また、図４Ｂに示されているように、シリンダ近傍を通過する流体の温度が相対的に大きく上昇しおり、顕著な温度分布が発生していることがわかる。定常状態に達した時点における流体の平均温度は２３５．１５Ｋであった。このように、顕著な温度分布が発生したのは、シミュレーションにおいて温度制御を行っていないためである。 As shown in FIG. 4A, when temperature control is not performed, it can be confirmed that Karman vortices are generated behind the cylinder. Further, as shown in FIG. 4B, it can be seen that the temperature of the fluid passing through the vicinity of the cylinder is relatively large, and a remarkable temperature distribution is generated. When the steady state was reached, the average temperature of the fluid was 235.15K. Thus, the remarkable temperature distribution was generated because temperature control was not performed in the simulation.
流速を考慮しないで温度制御を行った場合には、図５Ｂに示されているように、流体の温度は、図４Ｂの場合に比べて均一になっている。これは、温度制御が行われたことの効果である。ところが、図５Ａに示されているように、カルマン渦が確認できない。これでは、流体の挙動が適切にシミュレーションされたとはいえない。これは、温度制御時に流速を考慮していないため、温度制御が流速自体に影響を与えてしまったためである。 When temperature control is performed without considering the flow rate, the temperature of the fluid is uniform as compared with the case of FIG. 4B as shown in FIG. 5B. This is an effect of the temperature control being performed. However, Karman vortices cannot be confirmed as shown in FIG. 5A. In this case, it cannot be said that the behavior of the fluid is appropriately simulated. This is because the temperature control has influenced the flow rate itself because the flow rate is not taken into account during temperature control.
図６Ａに示されているように、実施例による方法で温度制御を行った場合には、カルマン渦が確認される。さらに、図６Ｂに示されているように、流体の温度がほぼ均一になっている。実施例による方法を採用することにより、流れ場を大きく乱すことなく、かつ適切な温度制御を行うことが可能であることが確認された。 As shown in FIG. 6A, when temperature control is performed by the method according to the embodiment, Karman vortices are confirmed. Furthermore, as shown in FIG. 6B, the temperature of the fluid is substantially uniform. By adopting the method according to the embodiment, it was confirmed that appropriate temperature control can be performed without greatly disturbing the flow field.
次に、上記実施例の変形例について説明する。上記実施例では、式（３）を用いて粒子の速度ｖｉを調整することにより温度制御を行った。その他に、速度に比例する摩擦力を与えて温度制御を行ってもよいし、Ｎｏｓｅ−Ｈｏｏｖｅｒ法を用いて温度制御を行ってもよい。 Next, a modification of the above embodiment will be described. In the above example, the temperature was controlled by adjusting the particle velocity v i using the equation (3). In addition, the temperature control may be performed by applying a frictional force proportional to the speed, or the temperature control may be performed using a Nose-Hoover method.
速度に比例する摩擦力を与えて温度制御を行う場合には、ステップ１３（図１）において動解析を行う際に、通常の運動方程式に代えて、以下の式（７）を用いることができる。
ここで、ｒｉは着目する粒子Ｐｉの位置、ｆｉは粒子Ｐｉに作用する力を表す。γはスケーリングの強さを決めるパラメータである。
When temperature control is performed by applying a frictional force proportional to speed, the following equation (7) can be used instead of the normal equation of motion when performing dynamic analysis in step 13 (FIG. 1). .
Here, r i represents the position of the particle P i of interest, and f i represents the force acting on the particle P i . γ is a parameter that determines the strength of scaling.
Ｎｏｓｅ−Ｈｏｏｖｅｒ法を用いて温度制御を行う場合には、通常の運動方程式に代えて、以下の式（８）を用いることができる。
ここで、ｒｉは粒子Ｐｉの位置、Φは粒子間ポテンシャル、ｇは空間の自由度を表す。Ｑはスケーリングの強さを決めるパラメータである。
When performing temperature control using the Nose-Hoover method, the following equation (8) can be used instead of the normal equation of motion.
Here, r i represents the position of the particle P i , Φ represents the interparticle potential, and g represents the degree of freedom in space. Q is a parameter that determines the strength of scaling.
分子動力学法によるシミュレーション方法において、計算対象の粒子数を減らすために、種々の方法が提案されている。例えば、流路の代表長さを短くする方法、繰り込みの手法を用いて粒子数を減らす方法（例えば、特許第５２４１４６８号公報）等により粒子数を減らすことができる。 In the simulation method based on the molecular dynamics method, various methods have been proposed in order to reduce the number of particles to be calculated. For example, the number of particles can be reduced by a method of shortening the representative length of the channel, a method of reducing the number of particles using a retraction method (for example, Japanese Patent No. 5241468), or the like.
シミュレーション対象である元粒子系の粒子数を減らすとき、シミュレーション対象である元粒子系と、粒子数を減らした粒子系とで、レイノルズ数が同一になるように粒子に与えられる物理量を対応付けることが好ましい。レイノルズ数を同一にすることにより、元粒子系と粒子数を減らした粒子系との間で、流れの相似性を確保することができる。 When reducing the number of particles in the original particle system to be simulated, it is possible to associate physical quantities given to the particles so that the Reynolds number is the same in the original particle system to be simulated and the particle system in which the number of particles is reduced. preferable. By making the Reynolds number the same, flow similarity can be ensured between the original particle system and the particle system with a reduced number of particles.
レイノルズ数Ｒｅは、以下の式（９）で与えられる。
ここで、μは流体の粘度を表し、ρは流体の密度を表し、ｖは粒子の速度を表し、Ｌは流路の代表長さを表す。例えば、粒子数を減らすために流路を小さくして代表長さＬを１／ｘ倍にした場合、粒子数を減らした粒子系の粒子の速度ｖをｘ倍にすればよい。
The Reynolds number Re is given by the following equation (9).
Here, μ represents the viscosity of the fluid, ρ represents the density of the fluid, v represents the velocity of the particles, and L represents the representative length of the flow path. For example, when the flow path is reduced to reduce the number of particles and the representative length L is increased by 1 / x, the particle velocity v of the particle system with the reduced number of particles may be increased by x.
繰り込みの手法を用いて粒子数を減らす場合には、例えば、繰り込みによって粒子系の密度ρ、代表長さＬが変化せず、粘度μがλ倍になる場合がある。ここで、λは繰り込み因子であり、繰り込み回数ｎを用いてλ＝２ｎと表される。この場合には、レイノルズ数Ｒｅを不変とするために、速度ｖをλ倍にすればよい。 When the number of particles is reduced by using the renormalization technique, for example, the density ρ and the representative length L of the particle system do not change due to the renormalization, and the viscosity μ may be λ times. Here, λ is a renormalization factor and is expressed as λ = 2n using the number n of renormalizations. In this case, in order to make the Reynolds number Re unchanged, the speed v may be multiplied by λ.
図７に、実施例によるシミュレーション装置のブロック図を示す。このシミュレーション装置として、シミュレーションプログラムを実行するコンピュータを用いることができる。このコンピュータは、中央処理装置（ＣＰＵ）４０、メインメモリ４１、補助記憶装置４２、入力装置４３、及び出力装置４４を含む。補助記憶装置４２は、シミュレーションプログラムが格納されている記録媒体を含む。この記録媒体は、補助記憶装置に内蔵されたものでもよく、補助記憶装置に対して着脱可能なリムーバブル媒体でもよい。このシミュレーションプログラムがメインメモリ４１にロードされ、ＣＰＵ４０によって実行される。 FIG. 7 shows a block diagram of a simulation apparatus according to the embodiment. A computer that executes a simulation program can be used as the simulation apparatus. The computer includes a central processing unit (CPU) 40, a main memory 41, an auxiliary storage device 42, an input device 43, and an output device 44. The auxiliary storage device 42 includes a recording medium in which a simulation program is stored. The recording medium may be a built-in auxiliary storage device or a removable medium that is detachable from the auxiliary storage device. This simulation program is loaded into the main memory 41 and executed by the CPU 40.
入力装置４３からシミュレーションに必要な情報が入力される。例えば、ステップ１０（図１）で設定される種々の情報が入力される。シミュレーション結果が出力装置４４に出力される。例えば、ステップ１５（図１）において、図６Ａに示した流体に生じた渦を可視化した画像、図６Ｂに示した流体の温度分布を示す画像等が出力装置４４に表示される。 Information necessary for the simulation is input from the input device 43. For example, various information set in step 10 (FIG. 1) is input. The simulation result is output to the output device 44. For example, in step 15 (FIG. 1), an image that visualizes the vortex generated in the fluid shown in FIG. 6A, an image showing the temperature distribution of the fluid shown in FIG. 6B, and the like are displayed on the output device 44.
図７に示したシミュレーション装置により、上記実施例によるシミュレーションを行うことができる。 The simulation according to the above embodiment can be performed by the simulation apparatus shown in FIG.
上述の実施例及び変形例は例示であり、異なる実施例及び変形例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。複数の実施例及び変形例の同様の構成による同様の作用効果については実施例及び変形例ごとには逐次言及しない。さらに、本発明は上述の実施例及び変形例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。 It is needless to say that the above-described embodiments and modifications are examples, and partial replacement or combination of configurations shown in different embodiments and modifications is possible. About the same effect by the same composition of a plurality of examples and a modification, it does not mention sequentially for every example and a modification. Furthermore, the present invention is not limited to the above-described embodiments and modifications. It will be apparent to those skilled in the art that various modifications, improvements, combinations, and the like can be made.
１０〜１５ ステップ
３０ 流速算出単位空間
４０ 中央処理装置（ＣＰＵ）
４１ メインメモリ
４２ 補助記憶装置
４３ 入力装置
４４ 出力装置
10-15 Step 30 Flow velocity calculation unit space 40 Central processing unit (CPU)
41 Main memory 42 Auxiliary storage device 43 Input device 44 Output device
Claims (6)
粒子ごとに、粒子の速度と、粒子の位置における流れ場の流速との差に基づいて温度制御を行うことにより、前記粒子系の温度を目標温度に近づけるシミュレーション方法。 A method of simulating a particle system including a plurality of particles forming a flow field by a molecular dynamics method,
A simulation method for bringing the temperature of the particle system closer to a target temperature by performing temperature control for each particle based on the difference between the velocity of the particle and the flow velocity of the flow field at the particle position.
算出された平均値に基づいて、粒子の各々の位置における温度を算出する請求項１または２に記載のシミュレーション方法。 In addition, the kinetic energy based on the difference between the velocity and flow velocity of the particle of interest and the average of the kinetic energy based on the difference between the velocity and flow velocity of other particles around the particle of interest are calculated. And
The simulation method according to claim 1, wherein the temperature at each position of the particle is calculated based on the calculated average value.
シミュレーショによって計算される対象の前記粒子系は、シミュレーションを行う対象である元粒子系に対して繰り込み処理を行って粒子数を減らしたものであり、
前記元粒子系の粒子に定義される物理量と、くりこまれた前記粒子系の粒子に定義される物理量とは、前記元粒子系とくりこまれた前記粒子系とでレイノルズ数が同一になるように対応付けられている請求項１乃至３のいずれか１項に記載のシミュレーション方法。 further,
The target particle system to be calculated by simulation is a renormalization process performed on the original particle system to be simulated to reduce the number of particles.
The physical quantity defined for the original particle system particle and the physical quantity defined for the renormalized particle system particle have the same Reynolds number in the original particle system and the regenerated particle system. The simulation method according to claim 1, wherein the simulation method is associated.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2016119871A JP6675785B2 (en) | 2016-06-16 | 2016-06-16 | Simulation method, program, and simulation device by molecular dynamics method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2016119871A JP6675785B2 (en) | 2016-06-16 | 2016-06-16 | Simulation method, program, and simulation device by molecular dynamics method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| JP2017224197A true JP2017224197A (en) | 2017-12-21 |
| JP6675785B2 JP6675785B2 (en) | 2020-04-01 |
Family
ID=60686396
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2016119871A Active JP6675785B2 (en) | 2016-06-16 | 2016-06-16 | Simulation method, program, and simulation device by molecular dynamics method |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP6675785B2 (en) |
Cited By (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2020101894A (en) * | 2018-12-20 | 2020-07-02 | 住友重機械工業株式会社 | Simulation method, simulation device, and program |
| JP2020101387A (en) * | 2018-12-20 | 2020-07-02 | 住友重機械工業株式会社 | Simulation method, simulation device, and program |
| JP7538728B2 (en) | 2021-01-15 | 2024-08-22 | 住友重機械工業株式会社 | Analysis method, analysis device, and program |
-
2016
- 2016-06-16 JP JP2016119871A patent/JP6675785B2/en active Active
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP2020101894A (en) * | 2018-12-20 | 2020-07-02 | 住友重機械工業株式会社 | Simulation method, simulation device, and program |
| JP2020101387A (en) * | 2018-12-20 | 2020-07-02 | 住友重機械工業株式会社 | Simulation method, simulation device, and program |
| JP7086483B2 (en) | 2018-12-20 | 2022-06-20 | 住友重機械工業株式会社 | Simulation method, simulation equipment and program |
| JP7086484B2 (en) | 2018-12-20 | 2022-06-20 | 住友重機械工業株式会社 | Simulation method, simulation equipment and program |
| JP7538728B2 (en) | 2021-01-15 | 2024-08-22 | 住友重機械工業株式会社 | Analysis method, analysis device, and program |
Also Published As
| Publication number | Publication date |
|---|---|
| JP6675785B2 (en) | 2020-04-01 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Greifzu et al. | Assessment of particle-tracking models for dispersed particle-laden flows implemented in OpenFOAM and ANSYS FLUENT | |
| Wan et al. | Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method | |
| Gualtieri et al. | Anisotropic clustering of inertial particles in homogeneous shear flow | |
| Heinen et al. | Droplet coalescence by molecular dynamics and phase-field modeling | |
| CN105069184A (en) | Stirred tank reactor simulation method based on immersed boundary method | |
| Jin et al. | A unified moving grid gas-kinetic method in Eulerian space for viscous flow computation | |
| Luo et al. | An improved direct-forcing immersed boundary method with inward retraction of Lagrangian points for simulation of particle-laden flows | |
| JP6675785B2 (en) | Simulation method, program, and simulation device by molecular dynamics method | |
| Akkermans et al. | Overset DNS with application to sound source prediction | |
| JP5892257B2 (en) | Simulation program, simulation method, and simulation apparatus | |
| JP2019070896A (en) | Simulation method, simulation device, and program | |
| JP2020194285A (en) | Information processing device, particle simulator system, and particle simulator method | |
| Chen et al. | Study on the performances of supply air for uniform air supply square hood by numerical simulation | |
| JP5888921B2 (en) | Simulation method and analysis apparatus | |
| US12152978B2 (en) | Controlling a multiphase flow | |
| Giahi et al. | Fully resolved simulation of spherical and non-spherical particles in a turbulent channel flow | |
| Winter et al. | Reduced-order modeling of transonic buffet aerodynamics | |
| US20140303946A1 (en) | Analysis device and simulation method | |
| Martín et al. | Modeling double concentric jets using linear and non-linear approaches | |
| JP6091402B2 (en) | Analysis apparatus and analysis method | |
| Wu | Application of the hybrid Local Domain Free Discretization and Immersed Boundary Method (LDFD-IBM) to simulate moving boundary flow problems | |
| JP7086484B2 (en) | Simulation method, simulation equipment and program | |
| JP7086483B2 (en) | Simulation method, simulation equipment and program | |
| Kumar et al. | Improved Methodology for Mass Conservation in Sharp Interface Immersed Boundary Method for Moving Boundaries | |
| JP6869621B2 (en) | Simulation method and simulation equipment |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20181113 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20200114 |
|
| A521 | Request for written amendment filed |
Free format text: JAPANESE INTERMEDIATE CODE: A523 Effective date: 20200227 |
|
| TRDD | Decision of grant or rejection written | ||
| A01 | Written decision to grant a patent or to grant a registration (utility model) |
Free format text: JAPANESE INTERMEDIATE CODE: A01 Effective date: 20200310 |
|
| A61 | First payment of annual fees (during grant procedure) |
Free format text: JAPANESE INTERMEDIATE CODE: A61 Effective date: 20200310 |
|
| R150 | Certificate of patent or registration of utility model |
Ref document number: 6675785 Country of ref document: JP Free format text: JAPANESE INTERMEDIATE CODE: R150 |