JP6675785B2 - Simulation method, program, and simulation device by molecular dynamics method - Google Patents
Simulation method, program, and simulation device by molecular dynamics method Download PDFInfo
- Publication number
- JP6675785B2 JP6675785B2 JP2016119871A JP2016119871A JP6675785B2 JP 6675785 B2 JP6675785 B2 JP 6675785B2 JP 2016119871 A JP2016119871 A JP 2016119871A JP 2016119871 A JP2016119871 A JP 2016119871A JP 6675785 B2 JP6675785 B2 JP 6675785B2
- Authority
- JP
- Japan
- Prior art keywords
- particle
- particles
- temperature
- simulation
- velocity
- 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.)
- Active
Links
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)
Description
本発明は、分子動力学法によるシミュレーション方法、プログラム、及びシミュレーション装置に関する。 The present invention relates to a simulation method, a program, and a simulation device based on a molecular dynamics method.
分子動力学を用いたシミュレーションが行われている。分子動力学では、シミュレーション対象の系を構成する粒子の運動方程式が数値的に解かれる。分子動力学を用いたシミュレーションを行う場合、シミュレーション対象の系の温度をある目標温度に設定したい場合がある。 Simulations using molecular dynamics have been performed. In molecular dynamics, the equations of motion of the particles making up the system to be simulated are solved numerically. 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 a temperature setting is changed in a simulation based on molecular dynamics. Is disclosed. In this method, the velocities of the respective particles obtained in the middle of the simulation by the molecular dynamics method are updated based on the current system temperature and the 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.
シミュレーション対象の系に流れ場が存在する場合、従来の温度制御方法では、流れ場の流速も温度制御にともなって変化してしまう。このため、流れ場が存在する系のシミュレーション時に従来の温度制御方法を適用すると、シミュレーション結果はシミュレーション対象の系の各粒子の挙動を反映したものにならない。従って、従来の温度制御方法は、流れ場が存在する系のシミュレーションに適用することができない。 If a flow field exists in the system to be simulated, the flow rate of the flow field changes with the temperature control in the conventional temperature control method. Therefore, if a conventional temperature control method is applied during a simulation of a system having a flow field, the simulation result does not reflect the behavior of each particle of the system to be simulated. Therefore, the conventional temperature control method cannot be applied to the 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 of a system having a flow field by a molecular dynamics method.
本発明の一観点によると、
流れ場を形成する複数の粒子を含む粒子系を分子動力学法によりシミュレーションする方法であって、
着目している粒子の速度と、着目している粒子の位置における流れ場の流速との差に基づく運動エネルギ、及び着目している粒子の周囲に存在する他の複数の粒子の速度と、前記流速との差に基づく運動エネルギの平均を算出し、
算出された平均値に基づいて、粒子の各々の位置における温度を算出し、
粒子ごとに、各粒子の位置における算出された温度と、前記粒子系の目標温度とに基づいて、温度制御を行うことにより、前記粒子系の温度を目標温度に近づけるシミュレーション方法が提供される。
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,
The velocity of the particle of interest, the kinetic energy based on the difference between the flow velocity of the flow field at the position of the particle of interest, and the velocities of other particles present around the particle of interest, Calculate the average of kinetic energy based on the difference with the flow velocity,
Based on the calculated average value, calculate the temperature at each position of the particles,
A simulation method for controlling the temperature of the particle system close to the target temperature by performing temperature control based on the calculated temperature at the position of each particle and the target temperature of the particle system for each particle is provided.
本発明の他の観点によると、上記シミュレーション方法をコンピュータに実行させるプログラム、及び上記シミュレーション方法を実行するシミュレーション装置が提供される。 According to another aspect of the present invention, there are provided a program for causing a computer to execute the above-mentioned simulation method, and a simulation apparatus for executing the above-mentioned simulation method.
粒子の速度と、粒子の位置における流れ場の流速との差に基づいて温度制御を行うことにより、流れ場の流速が温度制御の影響を受けにくくなる。これにより、シミュレーション対象の粒子系の流れ場の挙動がシミュレーション結果に反映されることとなる。 By performing the temperature control based on the difference between the velocity of the particles and the flow velocity of the flow field at the position of the particles, the flow velocity of the flow field 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. In the simulation according to the embodiment, the molecular dynamics method is applied. In the molecular dynamics method, a system to be simulated is represented by a plurality of particles, and a behavior of the particles, for example, a change in position and a velocity are obtained by numerically solving an 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 the 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 an execution interval of temperature control are input by an operation of an operator. The physical quantities defining the particle system include, for example, the mass of the particles, the potential of the particles, and the like. The initial state of the particles includes, for example, the position of the particles, the speed of the particles, the temperature of the particle system, and the like. The temperature of the particle system is reflected in the velocity distribution of the plurality of particles contained in the particle system. The constraint conditions include a boundary condition, a flow velocity of the flow field, a target temperature of the particle system, 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. Among a 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 execution interval of the temperature control.
次に、ステップ１１において、現時点で実行するタイムステップが、温度制御を実行するタイムステップに該当するか否かを判定する。例えば、温度制御の実行間隔の設定値で示された回数のタイムステップごとに、温度制御が実行される。 Next, in step 11, it is determined whether or not the time step to be executed at the present time corresponds to a time step to execute temperature control. For example, the temperature control is executed for each time step indicated by the set value of the execution interval of the temperature control.
現時点で実行するタイムステップが温度制御を実行するタイムステップに該当しない場合、ステップ１２において、通常の１タイムステップの動解析を実行する。現時点で実行するタイムステップが温度制御を実行するタイムステップに該当する場合、ステップ１３において、１タイムステップの動解析及び温度制御を実行する。 When the time step to be executed at the present time does not correspond to the time step to execute the temperature control, in step 12, a normal dynamic analysis of one time step is executed. When the time step to be executed at the present time corresponds to the time step to execute temperature control, in step 13, dynamic analysis and temperature control of one time step are executed.
ステップ１２及びステップ１３の後、ステップ１４においてシミュレーションを終了するか否かを判定する。シミュレーションを終了しない場合、ステップ１１に戻ってステップ１１からステップ１４までの処理を繰り返す。シミュレーションを終了する場合には、ステップ１５においてシミュレーション結果を出力する。 After Steps 12 and 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 repeats the processing from step 11 to step 14. When ending the simulation, a simulation result is output in step 15.
次に、ステップ１３の処理について詳細に説明する。
まず、分子動力学法を用いて、１タイムステップの動解析を実行する。１タイムステップの時間刻み幅は、予め設定されている。この動解析により、シミュレーション対象の粒子系の各粒子の位置及び速度が更新される。
Next, the process of step 13 will be described in detail.
First, dynamic analysis of one time step is executed using a molecular dynamics method. The time interval 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 for one time step, the flow velocity of the flow field at the position of each particle is calculated for all the particles. The flow velocity of the flow field can be represented by an average value of the velocities of the particle of interest and other particles existing around the particle of interest. As the plurality of particles serving as the basis for calculating the flow velocity, for example, particles existing inside a partial space (hereinafter, referred to as a flow velocity calculation unit space) including the particle of interest can be adopted. The flow velocity calculation unit space may be, for example, a sphere having a radius of about three to four times the average distance between the particles in the equilibrium state with the target particle as the center.
次に、図２Ａを参照して、流れ場の流速を算出する具体例について説明する。
図２Ａに、シミュレーション対象の粒子系の各粒子の位置及び速度を表した模式図を示す。図２Ａにおいて、各粒子を丸記号で表し、各粒子の速度ベクトルを矢印で表す。各粒子の速度ベクトルは、流れ場の下流側の方（図２Ａにおいて右方）を向く傾向を示す。
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 velocity of each particle of 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 tends to point downstream (rightward in FIG. 2A) of the flow field.
着目する粒子Ｐｉを含む流速算出単位空間３０を定義することができる。着目する粒子Ｐｉの速度をｖｉと表し、流速算出単位空間３０内の他の粒子Ｐｊの速度をｖｊと表す。流速算出単位空間３０内の着目する粒子Ｐｉ以外の粒子の個数をＮＮと表す。着目する粒子Ｐｉの位置における流れ場の流速Ｖｉは、以下の式で表される。
流速Ｖｉを算出した後、各粒子Ｐｉ、Ｐｊの速度ｖｉ、ｖｊと、着目する粒子Ｐｉの位置における流れ場の流速Ｖｉとの差に基づいて、着目する粒子Ｐｉの位置における温度を算出する。 After calculating the flow velocity V i , the particle P i of interest is determined based on the difference between the velocity v i , v j of each particle P i , P j and the flow velocity V i of the flow field at the position of the particle P i. The temperature at the position is calculated.
図２Ｂに、シミュレーション対象の粒子系の各粒子Ｐｉ、Ｐｊの位置、及び速度ｖｉ、ｖｊと流速Ｖｉとの差を表した模式図を示す。各粒子に付された矢印が、その粒子の速度ｖｉ、ｖｊと流速Ｖｉとの差を表している。粒子の速度ｖｉ、ｖｊと流速Ｖｉとの差を流速算出単位空間３０内の粒子について平均した速度は０になる。 In Figure 2B, shows each particle P i to be simulated in the particle system, the position of P j, and the speed v i, a schematic diagram showing the difference between 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 . Velocity v i of the particle, v j and average velocity for particles of the flow velocity calculating unit space 30 the difference between the flow velocity V i becomes zero.
着目する粒子Ｐｉの位置における温度をＴｉ、ボルツマン定数をｋＢ、粒子Ｐｉ、Ｐｊの質量をそれぞれｍｉ、ｍｊと表すと、以下の式（２）が成り立つ。式（２）を用いて、着目する粒子Ｐｉの位置における温度Ｔｉを算出することができる。
式（２）の右辺は、着目している粒子Ｐｉの速度ｖｉと流速Ｖｉとの差に基づく運動エネルギ、及び着目している粒子Ｐｉの周囲に存在する他の複数の粒子Ｐｊの速度ｖｊと流速Ｖｉとの差に基づく運動エネルギの平均値を表している。この平均値に基づいて、着目する粒子Ｐｉの位置における温度Ｔｉが算出される。 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 mean value of the kinetic energy based on the difference between the speed v j and the flow velocity V i of j. Based on the average value, the temperature Ti at the position of the particle P i 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 one method of adjusting the speed, a speed scaling method can be applied. In the velocity scaling method, if the velocity after the velocity adjustment of the particle P i of interest is represented as v im , the velocity v im is represented by the following equation.
Here, a is a parameter that determines the strength of the scaling, and is a positive real number of 1 or less. s i is referred to as a scaling factor. Tc is a target temperature.
粒子の速度を調整することにより、１タイムステップの動解析及び温度制御の処理が終了する。 By adjusting the speed of the particles, the processing of dynamic analysis and temperature control in one time step is completed.
次に、分子動力学法によるシミュレーションにおいて、流れ場の流速を考慮しないで温度制御を適用する方法と比較して、上記実施例の効果について説明する。流れ場の流速を考慮しないで温度制御を行う場合には、例えば、着目する粒子Ｐｉの位置における温度Ｔｉは上述の式（２）に代えて、以下の式（４）を用いて算出される。
さらに、温度制御による速度調整は、式（３）に代えて、以下の式（５）が適用される。
上述の式（４）の速度ｖｉ、ｖｊには、図２Ａに示したように流れ場の流速Ｖｉが含まれている。流速Ｖｉを含んだ速度ｖｉ、ｖｊに基づいて温度制御が行われるため、温度制御によって流れ場の流速Ｖｉが乱されてしまう。上記実施例による方法では、式（３）に示したように、粒子の速度ｖｉ、ｖｊと、粒子Ｐｉの位置における流れ場の流速Ｖｉとの差に基づいて温度制御が行われ、温度制御実施後の速度に流速Ｖｉが加算される。このため、温度制御による流速Ｖｉの乱れを抑制することができる。 The speed v i, v j of the aforementioned formula (4) includes a flow velocity V i of the flow field, as shown in Figure 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 embodiment, as shown in equation (3), the speed v i of the particle, v and j, the temperature control based on the difference between the velocity V i of the flow field at the position of the particle P i done , the flow velocity V i is added to the speed after the temperature control embodiment. Therefore, it is possible to suppress the disturbance of the flow velocity V i by the temperature control.
上記実施例の効果を確認するために、同一のシミュレーションモデルを用いて、温度制御を行わない方法、流速を考慮しないで温度制御を行う方法、及び実施例による方法の３通りの方法でシミュレーションを行った。次に、これらのシミュレーション結果について説明する。 In order to confirm the effects of the above embodiment, simulations were performed using the same simulation model by three methods: a method not performing temperature control, a method performing temperature control without considering the flow velocity, and a method according to the embodiment. went. Next, the results of these simulations 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 determined by simulation. An xy rectangular 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 diameter D of the cylinder was 500 × 10 −10 m. The cylinder is located at the center of the square area with respect to the y direction. The distance L2 from the negative side (x = 0) in the x direction of the square area 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. As the interaction potential between particles, only the term of the repulsive force of the Leonard Jones potential was used. Specifically, the interaction potential φ (r) is defined by the following equation.
Here, ε = 1.65 × 10 −21 J and σ = 3.41 × 10 −10 m. The viscosity of the fluid of this particle system is 3.64 × 10 −4 Pa · s, and the density is 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 arranged on the cylinder surface were subjected to temperature control by the Langevin method. The temperature of the cylinder was 88.85K.
温度制御を行う場合の目標温度Ｔｃとして、温度制御を行わない方法でシミュレーションしたときの定常状態到達時の平均温度である２３５．１５Ｋを用いた。 As the target temperature Tc in the case of performing the temperature control, 235.15K, which is the average temperature at the time of reaching the steady state when simulating by the method without performing the temperature control, was used.
次に、境界条件について説明する。右端（ｘ＝Ｌ１）と左端（ｘ＝０）とで特殊な周期境界条件が課されている。具体的には、粒子が右端の境界を通過するとき、ｘ方向の流速６００ｍ／ｓを中心とし、温度を８８．８５Ｋとするボルツマン分布に基づいて、ランダムに新しい速度が割り当てられる。この粒子が、左端の境界から正方形の領域内に流入する。下端（ｙ＝０）と上端（ｙ＝Ｌ１）の境界は、ｘ方向の流速が６００ｍ／ｓ、温度が８８．８５Ｋの熱浴として働く。下端及び上端の境界に到達した粒子は、ｘ方向の流速６００ｍ／ｓだけシフトしたハーフマクスウェル分布に基づいて、ランダムに新しい速度が割り当てられる。 Next, the boundary condition will be described. Special periodic boundary conditions are imposed at the right end (x = L1) and the left end (x = 0). Specifically, when a particle passes the rightmost boundary, a new velocity is randomly assigned based on a Boltzmann distribution centered on a flow velocity of 600 m / s in the x direction and a temperature of 88.85 K. The particles flow into the square area from the leftmost boundary. The boundary between the lower end (y = 0) and the upper end (y = L1) functions as a heat bath having a flow velocity in the x direction of 600 m / s and a temperature of 88.85K. Particles reaching the lower and upper boundaries are randomly assigned new velocities based on a half-Maxwell distribution shifted by a flow velocity of 600 m / s in the x-direction.
図４Ａ、図４Ｂ、図５Ａ、図５Ｂ、図６Ａ、及び図６Ｂにシミュレーショ結果を示す。図４Ａ及び図４Ｂは、温度制御を行わない方法を用いてシミュレーションを行った結果を示し、図５Ａ及び図５Ｂは、流速を考慮しないで温度制御を行う方法を用いてシミュレーションを行った結果を示し、図６Ａ及び図６Ｂは、実施例による方法でシミュレーションを行った結果を示す。 4A, 4B, 5A, 5B, 6A, and 6B show simulation results. 4A and 4B show the results of a simulation performed using a method that does not perform temperature control, and FIGS. 5A and 5B show the results of a simulation performed using a method that performs temperature control without considering the flow velocity. 6A and 6B show the results of a simulation performed by the method according to the embodiment.
図４Ａ、図５Ａ、及び図６Ａは、流体に生じた渦を可視化した図である。渦を可視化する方法は、例えば特許第５８８８９２１号公報に開示されている。図４Ｂ、図５Ｂ、及び図６Ｂは、温度分布を濃淡で表す図である。色の相対的に濃い領域が相対的に高温であることを表している。図４Ａ〜図６Ｂまでの各々において、左から右に並ぶ複数の図は、時間の経過に対応している。 4A, 5A, and 6A are diagrams visualizing a vortex generated in the fluid. A method for visualizing a vortex is disclosed in, for example, Japanese Patent No. 5888921. FIG. 4B, FIG. 5B, and FIG. 6B are diagrams showing the temperature distribution by shading. A relatively dark region indicates a relatively high temperature. In each of FIGS. 4A to 6B, a plurality of figures 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 vortex is generated behind the cylinder. Further, as shown in FIG. 4B, it can be seen that the temperature of the fluid passing in the vicinity of the cylinder has increased relatively greatly, and a remarkable temperature distribution has occurred. The average temperature of the fluid when the steady state was reached was 235.15K. Thus, the remarkable temperature distribution occurred because temperature control was not performed in the simulation.
流速を考慮しないで温度制御を行った場合には、図５Ｂに示されているように、流体の温度は、図４Ｂの場合に比べて均一になっている。これは、温度制御が行われたことの効果である。ところが、図５Ａに示されているように、カルマン渦が確認できない。これでは、流体の挙動が適切にシミュレーションされたとはいえない。これは、温度制御時に流速を考慮していないため、温度制御が流速自体に影響を与えてしまったためである。 When the temperature control is performed without considering the flow velocity, as shown in FIG. 5B, the temperature of the fluid is more uniform than in the case of FIG. 4B. This is an effect of performing the temperature control. However, as shown in FIG. 5A, Karman vortices cannot be confirmed. This does not mean that the behavior of the fluid has been properly simulated. This is because the flow rate was not taken into account during the temperature control, and the temperature control affected the flow rate itself.
図６Ａに示されているように、実施例による方法で温度制御を行った場合には、カルマン渦が確認される。さらに、図６Ｂに示されているように、流体の温度がほぼ均一になっている。実施例による方法を採用することにより、流れ場を大きく乱すことなく、かつ適切な温度制御を行うことが可能であることが確認された。 As shown in FIG. 6A, when temperature control is performed by the method according to the embodiment, Karman vortices are confirmed. Further, as shown in FIG. 6B, the temperature of the fluid is substantially uniform. It was confirmed that by employing the method according to the example, it was possible to perform appropriate temperature control without significantly disturbing the flow field.
次に、上記実施例の変形例について説明する。上記実施例では、式（３）を用いて粒子の速度ｖｉを調整することにより温度制御を行った。その他に、速度に比例する摩擦力を与えて温度制御を行ってもよいし、Ｎｏｓｅ−Ｈｏｏｖｅｒ法を用いて温度制御を行ってもよい。 Next, a modification of the above embodiment will be described. In the above embodiment, the temperature was controlled by adjusting the velocity v i of the particles using the equation (3). In addition, temperature control may be performed by applying a frictional force proportional to the speed, or temperature control may be performed using the Nose-Hoover method.
速度に比例する摩擦力を与えて温度制御を行う場合には、ステップ１３（図１）において動解析を行う際に、通常の運動方程式に代えて、以下の式（７）を用いることができる。
ここで、ｒｉは着目する粒子Ｐｉの位置、ｆｉは粒子Ｐｉに作用する力を表す。γはスケーリングの強さを決めるパラメータである。
When the temperature control is performed by giving a frictional force proportional to the speed, the following equation (7) can be used instead of the normal equation of motion when performing the dynamic analysis in step 13 (FIG. 1). .
Here, r i is the position of the particle P i of interest, f i represents the force acting on the particles 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 potential between particles, and g represents the degree of freedom in space. Q is a parameter that determines the strength of the scaling.
分子動力学法によるシミュレーション方法において、計算対象の粒子数を減らすために、種々の方法が提案されている。例えば、流路の代表長さを短くする方法、繰り込みの手法を用いて粒子数を減らす方法（例えば、特許第５２４１４６８号公報）等により粒子数を減らすことができる。 In the simulation method by the molecular dynamics method, various methods have been proposed to reduce the number of particles to be calculated. For example, the number of particles can be reduced by a method of reducing the representative length of the flow path, a method of reducing the number of particles by using a renormalization method (for example, Japanese Patent No. 5241468).
シミュレーション対象である元粒子系の粒子数を減らすとき、シミュレーション対象である元粒子系と、粒子数を減らした粒子系とで、レイノルズ数が同一になるように粒子に与えられる物理量を対応付けることが好ましい。レイノルズ数を同一にすることにより、元粒子系と粒子数を減らした粒子系との間で、流れの相似性を確保することができる。 When reducing the number of particles in the source particle system to be simulated, it is possible to associate the physical quantities given to the particles so that the Reynolds number is the same between the source particle system as the simulation target and the particle system with the reduced number of particles. preferable. By making the Reynolds numbers 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 channel. For example, when the representative length L is reduced to 1 / x times by reducing the flow path in order to reduce the number of particles, the velocity v of the particles in the particle system having the reduced number of particles may be increased by x times.
繰り込みの手法を用いて粒子数を減らす場合には、例えば、繰り込みによって粒子系の密度ρ、代表長さＬが変化せず、粘度μがλ倍になる場合がある。ここで、λは繰り込み因子であり、繰り込み回数ｎを用いてλ＝２ｎと表される。この場合には、レイノルズ数Ｒｅを不変とするために、速度ｖをλ倍にすればよい。 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 may not be changed by the renormalization, and the viscosity μ may be increased by λ. Here, λ is a renormalization factor, and is expressed as λ = 2 n using the number of renormalizations n. In this case, in order to make the Reynolds number Re invariable, the speed v may be multiplied by λ.
図７に、実施例によるシミュレーション装置のブロック図を示す。このシミュレーション装置として、シミュレーションプログラムを実行するコンピュータを用いることができる。このコンピュータは、中央処理装置（ＣＰＵ）４０、メインメモリ４１、補助記憶装置４２、入力装置４３、及び出力装置４４を含む。補助記憶装置４２は、シミュレーションプログラムが格納されている記録媒体を含む。この記録媒体は、補助記憶装置に内蔵されたものでもよく、補助記憶装置に対して着脱可能なリムーバブル媒体でもよい。このシミュレーションプログラムがメインメモリ４１にロードされ、ＣＰＵ４０によって実行される。 FIG. 7 shows a block diagram of a simulation device according to the embodiment. A computer that executes a simulation program can be used as the simulation device. This 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. This recording medium may be built in the auxiliary storage device, or may be a removable medium that is removable 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), the output device 44 displays an image visualizing 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.
図７に示したシミュレーション装置により、上記実施例によるシミュレーションを行うことができる。 The simulation according to the above embodiment can be performed by the simulation apparatus shown in FIG.
上述の実施例及び変形例は例示であり、異なる実施例及び変形例で示した構成の部分的な置換または組み合わせが可能であることは言うまでもない。複数の実施例及び変形例の同様の構成による同様の作用効果については実施例及び変形例ごとには逐次言及しない。さらに、本発明は上述の実施例及び変形例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。 The above-described embodiments and modified examples are exemplifications, and it is needless to say that the configurations shown in different embodiments and modified examples can be partially replaced or combined. The same operation and effect of the same configuration of a plurality of embodiments and modified examples will not be sequentially described for each embodiment and modified example. Furthermore, the present invention is not limited to the above-described embodiments and modifications. For example, it will be apparent to those skilled in the art that various modifications, improvements, combinations, and the like can be made.
１０〜１５ ステップ
３０ 流速算出単位空間
４０ 中央処理装置（ＣＰＵ）
４１ メインメモリ
４２ 補助記憶装置
４３ 入力装置
４４ 出力装置
10 to 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 (5)
着目している粒子の速度と、着目している粒子の位置における流れ場の流速との差に基づく運動エネルギ、及び着目している粒子の周囲に存在する他の複数の粒子の速度と、前記流速との差に基づく運動エネルギの平均を算出し、
算出された平均値に基づいて、粒子の各々の位置における温度を算出し、
粒子ごとに、各粒子の位置における算出された温度と、前記粒子系の目標温度とに基づいて、温度制御を行うことにより、前記粒子系の温度を目標温度に近づけるシミュレーション方法。 A method of simulating a particle system including a plurality of particles forming a flow field by a molecular dynamics method,
The velocity of the particle of interest, the kinetic energy based on the difference between the flow velocity of the flow field at the position of the particle of interest, and the velocities of other particles present around the particle of interest, Calculate the average of kinetic energy based on the difference with the flow velocity,
Based on the calculated average value, calculate the temperature at each position of the particles,
A simulation method for bringing the temperature of the particle system closer to the target temperature by performing temperature control for each particle based on the calculated temperature at the position of each particle and the target temperature of the particle system.
シミュレーショによって計算される対象の前記粒子系は、シミュレーションを行う対象である元粒子系に対して繰り込み処理を行って粒子数を減らしたものであり、
前記元粒子系の粒子に定義される物理量と、くりこまれた前記粒子系の粒子に定義される物理量とは、前記元粒子系とくりこまれた前記粒子系とでレイノルズ数が同一になるように対応付けられている請求項１または２に記載のシミュレーション方法。 further,
The particle system to be calculated by the simulation is a reduction in the number of particles by performing a renormalization process on the original particle system to be simulated,
The physical quantity defined for the particles of the original particle system and the physical quantity defined for the renormalized particles of the particle system are such that the Reynolds number is the same in the original particle system and the renormalized particle system. simulation method according to claim 1 or 2 is associated.
着目している粒子の速度と、着目している粒子の位置における流れ場の流速との差に基づく運動エネルギ、及び着目している粒子の周囲に存在する他の複数の粒子の速度と、前記流速との差に基づく運動エネルギの平均を算出する手順と、
算出された平均値に基づいて、粒子の各々の位置における温度を算出する手順と、
粒子ごとに、各粒子の位置における温度が前記粒子系の目標温度に近づくように温度制御を行うことにより、前記粒子系の温度を目標温度に近づける手順と
を実行するシミュレーション装置。
When determining the position and velocity of each of the particles of a particle system including a plurality of particles forming a flow field by a molecular dynamics method at a certain time step,
The velocity of the particle of interest, the kinetic energy based on the difference between the flow velocity of the flow field at the position of the particle of interest, and the velocities of other particles present around the particle of interest, Calculating an average of the kinetic energy based on the difference with the flow velocity;
Based on the calculated average value, a procedure for calculating the temperature at each position of the particles,
A simulation apparatus for performing a procedure of controlling the temperature of the particle system to a target temperature by performing temperature control such that a temperature at a position of each particle approaches a target temperature of the particle system for each particle .
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 JP2017224197A (en) | 2017-12-21 |
| JP6675785B2 true 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) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| 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 |
-
2016
- 2016-06-16 JP JP2016119871A patent/JP6675785B2/en active Active
Also Published As
| Publication number | Publication date |
|---|---|
| JP2017224197A (en) | 2017-12-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Wan et al. | Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method | |
| CN109643334B (en) | Reverse simulation of multiple fibers | |
| KR100984048B1 (en) | An effective method to solve rigid body interactions in particle based fluid simulations | |
| JP5441422B2 (en) | Simulation method and program | |
| JP6675785B2 (en) | Simulation method, program, and simulation device by molecular dynamics method | |
| JP2009069929A (en) | Method for constructing surface of fluid simulation based on particle method, program and storage medium with its program stored | |
| JP5872324B2 (en) | Mesh generator | |
| JP5892257B2 (en) | Simulation program, simulation method, and simulation apparatus | |
| JP2019070896A (en) | Simulation method, simulation device, and program | |
| KR102459848B1 (en) | Modeling method and modeling apparatus of target object at high speed based on particles | |
| JP7246636B2 (en) | Information processing device, particle simulator system, and particle simulator method | |
| CN109434840A (en) | A kind of robot free path generation method based on spline curve | |
| JP5888921B2 (en) | Simulation method and analysis apparatus | |
| JP7192603B2 (en) | Program editing device, program editing method, and program editing program | |
| WO2015174483A1 (en) | Three-dimensional model creation device, three-dimensional model creation method, and program | |
| JP7086484B2 (en) | Simulation method, simulation equipment and program | |
| EP2800019A1 (en) | Analysis device and simulation method | |
| WO2022239441A1 (en) | Analysis device and program | |
| JP2002092639A (en) | Method and device for forming animation representing particle behavior | |
| JP5122410B2 (en) | Design verification display device, design verification display method, and design verification display computer program | |
| JP5720551B2 (en) | Simulation program, simulation method, and simulation apparatus | |
| JP6869621B2 (en) | Simulation method and simulation equipment | |
| JPH0561961A (en) | System and device for interactive processing of computer animation | |
| JP2007323585A (en) | Polygon data division method and device | |
| Kumar et al. | Improved Methodology for Mass Conservation in Sharp Interface Immersed Boundary Method for Moving Boundaries |
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 | Written amendment |
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 |