WO2014045416A1 - シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 - Google Patents

シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 Download PDF

Info

Publication number
WO2014045416A1
WO2014045416A1 PCT/JP2012/074295 JP2012074295W WO2014045416A1 WO 2014045416 A1 WO2014045416 A1 WO 2014045416A1 JP 2012074295 W JP2012074295 W JP 2012074295W WO 2014045416 A1 WO2014045416 A1 WO 2014045416A1
Authority
WO
WIPO (PCT)
Prior art keywords
fluid
interface
simulation
surface tension
calculated
Prior art date
Application number
PCT/JP2012/074295
Other languages
English (en)
French (fr)
Inventor
風間正喜
Original Assignee
富士通株式会社
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by 富士通株式会社 filed Critical 富士通株式会社
Priority to PCT/JP2012/074295 priority Critical patent/WO2014045416A1/ja
Priority to JP2014536514A priority patent/JP5892257B2/ja
Publication of WO2014045416A1 publication Critical patent/WO2014045416A1/ja
Priority to US14/621,896 priority patent/US10031984B2/en

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the present invention relates to a simulation program, a simulation method, and a simulation apparatus.
  • the handling is complicated when a moving boundary such as a free surface or a fluid / structure interaction problem occurs such as a free surface. In many cases, it is difficult to create a program.
  • Patent Document 1 is a technique of constructing a potential from the interparticle distance and applying an attractive force between the particles. This is a technique based on the image that surface tension is caused by non-uniformity of intermolecular forces on the surface.
  • Non-Patent Document 1 Non-Patent Document 2
  • Non-Patent Document 3 a model expressing wettability is also introduced.
  • a model for expressing the effect of surface tension that affects only the surface by the particle method is implemented in such a way that particles other than the surface are also affected, and non-physical behavior occurs. In particular, it becomes difficult to control the dynamic movement of the wetting angle.
  • Non-Patent Document 1 For example, in the calculation of the surface tension by the methods of Patent Document 1, Non-Patent Document 1, Non-Patent Document 2, and Non-Patent Document 3, not only the surface but also the internal particles are affected.
  • Non-Patent Document 1, Non-Patent Document 2, and Non-Patent Document 3 are methods for predicting particles on the surface from information on surrounding particles and introducing surface tension.
  • the particles inside the fluid are also determined to be surface particles, and surface tension is applied. For this reason, the contact surface cannot be captured, and handling such as when exercising while having a wetting angle becomes difficult.
  • the distribution of particles (Where i represents the particle number), in Non-Patent Document 1 and Non-Patent Document 2, it is determined that the particle i is a surface particle, and the surface tension is calculated by the following steps 1 to 3. Run with.
  • Non-Patent Document 1 and Non-Patent Document 2 use the following spline functions.
  • Equation 1 is the center of gravity of the particle within the radius h centered on the particle i, weighted by a kernel function, as shown in FIG.
  • the surface tension is applied by considering the particles with non-uniform distribution of the surrounding particles as the boundary particles. Exceeding "0.15 times the average particle spacing" The surface tension will be applied even if it has.
  • the present invention aims to suppress simulation results from exhibiting non-physical behavior.
  • the present invention employs the following configuration in order to solve the above problems. That is, according to one aspect of the present invention, in a simulation program, a simulation method, or a simulation apparatus for simulating the surface tension of a fluid, the fluid is collected as a group of particles based on boundary conditions and initial conditions input to a computer. Calculate the interface of the expressed fluid model, calculate the surface energy of the calculated interface, calculate the surface tension of the interface based on the calculated surface energy, and based on the calculated surface tension A simulation program, a simulation method, or a simulation apparatus is provided that outputs the fluid state at predetermined time intervals.
  • the calculation of the surface energy of the interface is expressed as the sum of the first term for calculating the surface energy of the interface with air and the second term for calculating the surface energy of the interface other than air. It is desirable to use a calculation formula.
  • the calculation of the surface energy at the interface with the air changes the surface tension coefficient of the first term and the surface tension coefficient of the second term.
  • the calculation of the surface energy at the interface with the air adds energy for correcting the nonuniformity of the surface distribution of the particles.
  • an appropriate simulation result can be output without showing non-physical behavior.
  • the present invention relates to a simulation program, a simulation method, and a simulation apparatus that are executed by a computer, capturing a fluid to be simulated as a collection of particles, and based on the particle distribution (point distribution) based on the convex hull construction method.
  • the interface is obtained from the particles at the interface constituting the boundary between the fluid to be simulated and the gas or solid other than the fluid.
  • the surface tension term is calculated by expressing the surface energy using the particles constituting the interface and calculating the first variation of the surface energy.
  • each of a liquid phase that is in a liquid state, a gas phase that is in a gaseous state, and a solid phase that is in a solid state that is in contact with the fluid that is the target of the particle method simulation via an interface is called “other phase”.
  • FIG. 5 is a diagram showing a configuration example of a simulation apparatus to which the present invention is applied.
  • the simulation apparatus 500 includes a processing unit 501, a storage unit 502, and an output unit 503, performs numerical calculation by a particle method based on an initial condition 511, and outputs a simulation result 512.
  • the storage unit 502 stores information on each calculation formula for executing a simulation program to which the present invention is applied.
  • the processing unit 501 executes a simulation process to which the present invention is applied, which will be described later as the first to third embodiments.
  • the output unit 503 outputs the simulation result 512 executed by the processing unit 501.
  • a two-dimensional interface extraction method will be described as the first embodiment.
  • the two-dimensional interface extraction in the first embodiment is realized by first obtaining a closed curve including two-dimensional particle groups and then obtaining a recess on the closed curve.
  • Two-dimensional particle group think of.
  • i represents the particle number.
  • FIG. 6 is a diagram for explaining a method for obtaining a closed curve including a point group.
  • a method for obtaining a closed curve in the first to fourth procedures using the Gift Wrapping method which is one of convex hull construction methods will be described.
  • Point sequence obtained by such a procedure As shown in FIG. 7, a closed curve that becomes a convex hull can be obtained.
  • FIG. 8 is a flowchart showing a flow of processing for obtaining a closed curve including a point group.
  • step S81 particle input data is acquired. This step S81 corresponds to the first procedure.
  • step S82 initial particles and reference line segments that are the second procedure are determined.
  • step S83 the angle formed between the reference line segment and the particle is calculated for all particles.
  • step S84 the particle having the smallest angle is defined as the boundary particle, and the reference component is updated in step S85.
  • steps S84 and S85 correspond to the third procedure and the fourth procedure.
  • FIG. 11 is a flowchart showing a flow of processing for obtaining a recess.
  • step S111 surface particle data by the convex hull construction method described with reference to FIGS. 6 to 8 is acquired, and in step S112, data of the surface particle i and the next surface particle i + 1 are acquired.
  • step S113 the distance between the particle i and the particle i + 1 is calculated. This step S113 corresponds to the fifth procedure.
  • step S114 it is determined whether the distance between the surface particles is longer than a predetermined threshold. If long (step S114: Yes), candidate k is selected from other particles in step S115, and k is set as i + 1, and the processing in step S113 and subsequent steps is repeated.
  • the flow from Step S114 (Yes) and Step S115 to Step S113 corresponds to the sixth to eighth procedures.
  • or step S114 are performed with respect to all the boundary particles (i). This corresponds to the ninth procedure.
  • a point sequence is used as the position of the surface particle.
  • the surface is composed of i, i + 1 particles as shown in FIG.
  • the surface is composed of triangular elements between three points.
  • the surface energy of a fluid can be defined as follows.
  • the integration region on the right side of the above (Equation 4) is the entire boundary surface of the fluid, Takes a value of 1 when the fluid is in contact with the air and 0 when the fluid is not touched.
  • Respectively, as shown in FIG. 13, are the partial surface tension coefficient in contact with air and the surface tension coefficient in the part in contact with the solid.
  • the point of the present invention is how to express the above (Formula 4) by the particle method.
  • the point sequence obtained by the interface extraction method Is used to determine the surface energy as follows.
  • Is a three-dimensional point Is a function that represents the distance between a solid and a solid. For example, if there is a solid wall in the plane of y 0 It is.
  • a function form obtained by performing some kind of interpolation (such as linear interpolation or interpolation by a cubic function) is used.
  • Integration range I a line segment as shown in FIG. Taken up.
  • the integral calculation can be numerically integrated using a Gaussian integration method or the like.
  • the sum (j) of the first term on the right side of (Formula 5) is taken between particles (between j and j + 1 particles), and the sum (k) of the second term is taken for all particles.
  • the first term on the right side is the one obtained by discretizing the surface energy of (Equation 4), while the second term on the right side is the potential energy representing the boundary condition that the particle does not pass through the wall.
  • the energy in the case of being in contact with the solid is lower in energy, so that a force to adhere to the solid is generated.
  • the particle method since the condition of not passing through the solid is expressed as a potential near the wall (Equation 6), the particle does not come into full contact with the solid. Therefore, in the second embodiment, as shown in FIG. 17, energy is set as shown in (Equation 5), and the distance from the solid is In the following cases, the repulsive force is received and the distance from the solid is more than In the case of less than, the adhesive force was made to work.
  • FIG. 18 is a flowchart showing the flow of calculation of the surface tension term.
  • step S181 data of a fluid to be simulated is acquired, and in step S182, the position of the fluid particle is updated by time dt / 2 using the current velocity.
  • step S183 the force applied to the fluid particles is calculated from the physical model.
  • step S184 fluid surface particle data is acquired according to the first embodiment described above.
  • step S185 surface energy is obtained as described above for all boundary particles, and the force applied to the surface particles is calculated.
  • step S186 the velocity of the fluid particles is updated in step S186, and the position of the fluid particles is updated for the time dt / 2 using the updated velocity in step S187.
  • the above steps S182 to S187 are repeatedly executed. Further, in place of the above energy, the following volume energy and contact energy can be introduced in order to suppress surface particle non-uniformity and fluid volume change.
  • Is the volume of the fluid Can be calculated.
  • Is the outward normal unit vector Is a constant
  • Is the mass of the surface particles s i.
  • the correction terms based on the energy in (Equation 8) and (Equation 9) can be calculated as follows. Is the initial volume and is a constant.
  • the correction term by volume energy is the volume as shown in FIG.
  • the correction term based on the contact energy has the effect of suppressing the non-uniformity of the surface particles as shown in FIG.
  • FIG. 21 is a diagram illustrating an example of a two-dimensional fluid.
  • Example 11 is a law of conservation of mass
  • Example 12 is a law of conservation of momentum
  • (Expression 13) is an equation of state. Is the viscous stress tensor, which determines the viscosity coefficient of the fluid (Constant) It is.
  • Equation 11 represents the effect that the density increases when there is a velocity field where fluids gather, and conversely the density decreases when the velocity field moves away from the fluid.
  • the first term on the right side of (Equation 12) is a pressure gradient term, which represents the effect that a force is generated from a portion where the fluid is high pressure to a portion where the pressure is low.
  • the second term on the right side is a viscous stress term that represents the effect of braking the flow.
  • the third term on the right side is the gravity term. Is the acceleration of gravity, Is a unit vector in the y direction.
  • the subscript represents the particle number. That is, the position vector, velocity vector, density, and pressure of the ath particle are It is. Is the mass of the b-th particle. Is a constant introduced by the parameter to calculate the viscosity term.
  • Equation 15 is obtained by discretizing the equation of motion of the above (Equation 12) by the particle method, the second term on the right side represents the pressure gradient term, and the third term represents the viscous stress term.
  • Is Y-direction component of Is the surface energy and has the same format as in the above (formula 5), and adopts the following discrete format.
  • the contact energy is as follows.
  • volume of fluid Is Can be calculated. here Is a vector A matrix to rotate.
  • I a two-dimensional renormalization matrix
  • each slope is calculated as follows.
  • FIG. 22 and 23 express wettability with a solid surface by applying surface tension by the above-described method.
  • FIG. so When so, Indicates the time.
  • the contact angle difference can be expressed by changing the surface tension coefficient.
  • FIG. 24 shows the case where the correction term calculated by (Equation 19) and (Equation 20) is not entered
  • FIG. 25 shows the case where the correction term is entered.
  • the simulation parameters are as follows.
  • Gravity acceleration is 9.8
  • the number of particles was 900, and the initial shape was made into a square by placing particles on a lattice having a length of 0.00167 [m].
  • the present invention is applicable to the calculation of fluid motion in which the surface tension is effective in the calculation of the particle method. This is particularly effective for molten metal pouring and simulation of resin pouring.
  • the simulation apparatus of FIG. 5 can be realized using, for example, an information processing apparatus (computer) as shown in FIG. 26 includes a CPU (Central Processing Unit) 2601, a memory 2602, an input device 2603, an output device 2604, an external storage device 2605, a medium drive device 2606, and a network connection device 2607. These are connected to each other by a bus 2608.
  • a CPU Central Processing Unit
  • the memory 2602 is a semiconductor memory such as a ROM (Read Only Memory), a RAM (Random Access Memory), a flash memory, and the like, and stores programs and data used for simulation processing.
  • the CPU 2601 performs the above-described simulation process by executing a program using the memory 2602.
  • the memory 2602 can also be used as the storage unit 502 in FIG.
  • the input device 2603 is, for example, a keyboard, a pointing device, etc., and is used for inputting instructions and information from the operator.
  • the output device 2604 is, for example, a display device, a printer, a speaker, or the like, and is used to output an inquiry to the operator or a processing result.
  • the output device 2604 can also be used as the output unit 503 in FIG.
  • the external storage device 2605 is, for example, a magnetic disk device, an optical disk device, a magneto-optical disk device, a tape device, or the like.
  • the external storage device 2605 includes a hard disk drive.
  • the information processing apparatus can store the program and data in the external storage device 2605 and load them into the memory 2602 for use.
  • the medium driving device 2606 drives the portable recording medium 2609 and accesses the recorded contents.
  • the portable recording medium 2609 is a memory device, a flexible disk, an optical disk, a magneto-optical disk, or the like.
  • This portable recording medium 2609 has Compact Disk Read Only. Memory (CD-ROM), Digital Versatile Disk (DVD), Universal Serial Bus (USB) memory, etc. are also included.
  • the operator can store programs and data in the portable recording medium 2609 and load them into the memory 2602 for use.
  • the computer-readable recording medium for storing the program and data used for the simulation processing is a physical (non-temporary) such as the memory 2602, the external storage device 2605, and the portable recording medium 2609. )
  • a recording medium is included.
  • the network connection device 2607 is a communication interface that is connected to the communication network 2610 and performs data conversion accompanying communication.
  • the information processing apparatus can receive a program and data from an external apparatus via the network connection apparatus 2607, and can use them by loading them into the memory 2602.
  • the network connection device 2607 can also be used as the output unit 503 in FIG.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

　流体の表面張力をシミュレーションするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置において、コンピュータに、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、前記算出された界面の表面エネルギーを算出させ、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置が提供される。

Description

シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
　本発明は、シミュレーションプログラム、シミュレーション方法及びシミュレーション装置に関する。
　近年の計算機パワーの向上に伴い、シミュレーションの手法も徐々に発展してきた。その結果として、様々な応用分野にシミュレーションが用いられるようになってきた。
　流体や弾性体等の連続体の問題を解く数値計算手法としては、格子をベースにして微分方程式の近似解を求解する有限差分法、有限要素法、又は有限体積法等が多く用いられてきた。近年では数値計算をＣＡＥ（Computer Aided Engineering）などの応用分野で活用するため、これらの数値計算手法も発展し、流体と構造物が相互作用する問題が解かれるようになってきた。
　しかしながら、有限要素法又は有限体積法等の格子を用いる手法では、自由表面などの界面の存在する問題や、流体・構造連成問題などの移動境界が発生する場合には、取り扱いが複雑なため、プログラム作成が困難である場合が多い。
　これに対して格子を用いないＭＰＳ（Moving Particle Semi-implicit）法、ＳＰＨ（Smoothed Particle Hydrodynamics）法等の粒子法では、移動境界の取り扱いに特別な処置を必要としない。それゆえ、近年粒子法は広く用いられるようになっている。
　自由表面などの、変形し移動する境界を容易に扱いたいがために開発された粒子法であるが、離散化方法から、連続体が単なる粒子群となってしまう場合には、図１に示すように、連続体の境界面がどこかはっきりしなくなる。そのため、粒子法では表面張力などの境界を陽に扱う必要がある問題に対しては、統一的手法は出来ておらず、個別に独自開発された技術が、いつくかあるという状態である（例えば、特許文献１、非特許文献１、非特許文献２、非特許文献３参照。）。
　例えば、特許文献１は、図２に示すように、粒子間距離からポテンシャルを構成し、粒子間に引力を働かせる手法である。表面張力は表面での分子間力の不均一性から生じるものである、という描像を元に作られた手法である。
　また、非特許文献１、非特許文献２、非特許文献３の手法は、近傍粒子の数が減ったものを境界にある粒子と見なして、表面張力を境界の粒子に加える手法である。
　また、非特許文献２、非特許文献３では、濡れ性を表現するモデルも導入されている。
特開２００８－１１１６７５号公報
T. Hongo, M. Shigeta, S.  Izawa, and Y. Fukunishi, "３次元非圧縮ＳＰＨ法における気液界面に作用する表面張力モデル", 第２３回数値流体力学シンポジウム, A8-5 (2009) M. Agawa, M. Shigeta, S.  Izawa, and Y. Fukunishi, "斜面上を流下する液体の非圧縮ＳＰＨシミュレーション", 第２３回数値流体力学シンポジウム, A9-4 (2009) K. Nomura, S. Koshizuka, Y. Oka and H. Obata, "Numerical Analysis of Droplet Breakup Behavior using Particle Method", Journal of Nuclear Science and Technology, Vol. 38, No. 12, pp. 1057-1064(2001)
　しかしながら、上述のような従来の技術においては、下記のような問題点がある。
　すなわち、表面のみに影響する表面張力の効果を、粒子法で表すためのモデルが、表面以外にある粒子も影響を受けるような形で実装されており、非物理的な挙動が発生する。特に濡れ角の動的な運動を制御することが難しくなってしまう。
　例えば、特許文献１、非特許文献１、非特許文献２、非特許文献３の手法による表面張力の計算では、表面だけではなく、内部の粒子にも影響が及んでしまう。
　特に、特許文献１の手法では、全粒子に対して、引力を働かせているため、内部の粒子にも表面張力による力がかかることになる。
　また、非特許文献１、非特許文献２、非特許文献３の手法は、周りの粒子の情報から、表面にある粒子を予想し、表面張力を導入する方法であるが、計算の中で、流体内部の粒子も、表面の粒子であると判定されてしまい、表面張力が加えられてしまうことが起こる。そのため、接触面を捕らえることが出来なくなり、濡れ角を持ちながら運動する場合などの扱いが困難となる。
　例えば、粒子の分布を
Figure JPOXMLDOC01-appb-M000001
 （ここで、ｉは粒子の番号を表す）とすると、非特許文献１、非特許文献２では粒子ｉが表面の粒子であると、判定して表面張力の計算を、以下の手順１乃至３で実行する。
（手順１）
　以下の量を計算する。
Figure JPOXMLDOC01-appb-M000002
ここで、右辺の
Figure JPOXMLDOC01-appb-I000003
はカーネル関数（重み関数と呼ばれることもある）で、非特許文献１、非特許文献２では以下のスプライン関数を用いている。
Figure JPOXMLDOC01-appb-M000004
ここで、
Figure JPOXMLDOC01-appb-I000005
は粒子間の影響半径で、初期状態の平均粒子間隔の２倍から３倍程度が良く用いられる。
Figure JPOXMLDOC01-appb-I000006
はカーネル関数の全空間積分量が１になるように調整された値で、２次元の場合は
Figure JPOXMLDOC01-appb-M000007
 、３次元の場合は
Figure JPOXMLDOC01-appb-M000008
 と決められる。
　上記（式１）は、図３に示すように、カーネル関数による重みを付けた、粒子iを中心とした半径ｈ内の粒子の重心である。
（手順２）
　（式１）の重心と粒子iの距離を計算する。
Figure JPOXMLDOC01-appb-M000009
（手順３）
Figure JPOXMLDOC01-appb-I000010
の値が平均粒子間隔の所定倍数、例えば０．１５倍以上であれば粒子ｉを表面の粒子と見なして、表面張力を加える。
　これらの手順１乃至３では、図４に示すように、周りの粒子分布が不均一になった粒子を境界の粒子と見なして表面張力を加えているが、内部の粒子が計算中にたまたま閾値“平均粒子間隔の０．１５倍”を超える
Figure JPOXMLDOC01-appb-I000011
を持ってしまった場合にも表面張力を適用してしまうことになる。
　１つの側面では、本発明は、シミュレーション結果が非物理的な挙動を示すことを抑制することを目的とする。
　本発明は、上記課題を解決するため、下記のような構成を採用した。
　すなわち、本発明の一態様によれば、流体の表面張力をシミュレーションするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置において、コンピュータに、入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、前記算出された界面の表面エネルギーを算出させ、前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム、シミュレーション方法又はシミュレーション装置が提供される。
　また、本発明は、前記界面の表面エネルギーの算出が、空気との界面の表面エネルギーを算出する第１項と前記空気以外との界面の表面エネルギーを算出する第２項との和で表現した計算式を用いることが望ましい。
　また、本発明は、前記空気との界面の表面エネルギーの算出が、前記第１項の表面張力係数と前記第２項の表面張力係数とを変えることが望ましい。
　また、本発明は、前記空気との界面の表面エネルギーの算出が、前記粒子の表面分布の不均一性を補正するためのエネルギーを加算することが望ましい。
　また、本発明は、前記空気との界面の表面エネルギーの算出が、前記流体の体積を補正するためのエネルギーを加算することが望ましい。
　開示されたシミュレーションプログラム、シミュレーション方法及びシミュレーション装置によれば、非物理的な挙動を示すことなく、適切なシミュレーション結果を出力することができる。
粒子法の問題点を説明するための図である。 粒子間の引力を説明するための図である。 重みを付けた粒子の重心を示す図である。 境界粒子の重心を示す図である。 本発明を適用したシミュレーション装置の構成例を示す図である。 点群を包括する閉曲線を求める方法を説明するための図である。 凸包となる閉曲線を示す図である。 点群を包括する閉曲線を求める処理の流れを示すフローチャートである。 凹部を求める方法を説明するための図である。 凹部を求める方法により修正された閉曲線を示す図である。 凹部を求める処理の流れを示すフローチャートである。 流体の表面粒子の点列を示す図である。 表面張力係数を説明するための図である。 表面エネルギーのΧ（カイ）を幅εで滑らかに近似した関数を説明するための図である。 線分上での積分を説明するための図である。 接触核の表現を説明するための図である。 固体からの距離とエネルギーの関係を示す図である。 表面張力項の計算の流れを示すフローチャートである。 体積を戻す効果が示された図である。 表面粒子の負均一性を抑える効果が示された図である。 ２次元流体の例を示す図である。 表面張力を加えて固体面との濡れ性を表現した図（その１）である。 表面張力を加えて固体面との濡れ性を表現した図（その２）である。 補正項を入れなかった場合の例を示す図である。 補正項を入れた場合の例を示す図である。 情報処理装置の構成図である。
　以下、本発明の実施の形態について、図面を参照しながら詳細に説明する。
　本発明は、コンピュータに実行させるシミュレーションプログラム、シミュレーション方法及びシミュレーション装置であって、シミュレーションの対象である流体を粒子の集まりとして捕らえ、その粒子の分布（点の分布）から、凸包構成法を基にした計算幾何手法により、当該シミュレーションの対象である流体と当該流体以外の気体又は固体との境界を構成する界面にある粒子から界面を求める。そして、界面を構成する粒子を用いて表面エネルギーを表現し、表面エネルギーの第１変分を計算することにより、表面張力項の計算を行う。
　ここで、粒子法シミュレーションの対象である流体と界面を介して接する、液体状態にある相である液相、気体状態にある相である気相、固体状態にある相である固相の各々を「他相」という。
　また、表面エネルギーを構成する際に、表面張力の大きさを、各面毎に変えることで、流体の濡れ性や、付着現象も表現できるモデルを採用し、差分化する。
　図５は、本発明を適用したシミュレーション装置の構成例を示す図である。
　図５において、シミュレーション装置５００は、処理部５０１、格納部５０２及び出力部５０３を備え、初期条件５１１に基づいて粒子法による数値計算を行って、シミュレーション結果５１２を出力する。
　格納部５０２は、本発明を適用したシミュレーションプログラムを実行するための各計算式の情報を格納する。
　処理部５０１は、第１の実施の形態乃至第３の実施の形態として後述する、本発明を適用したシミュレーション処理を実行する。
　出力部５０３は、処理部５０１が実行したシミュレーション結果５１２を出力する。
（第１の実施の形態）
　第１の実施の形態として、２次元の界面抽出法について説明する。
　本第１の実施の形態における２次元の界面抽出は、まず、２次元上の粒子群を包括する閉曲線を求め、次に、閉曲線上の凹部を求めることで実現する。
　まず、閉曲線の求め方を説明する。
　２次元上の粒子群（点群）
Figure JPOXMLDOC01-appb-I000012
を考える。ここで、ｉは粒子の番号を表す。
　このような点群を包括する閉曲線を求める方法を説明する。
　図６は、点群を包括する閉曲線を求める方法を説明するための図である。
　本第１の実施の形態では、凸包構成法の１つであるGift Wrapping法を用いて、第１乃至第４の手順で閉曲線を求める方法を説明する。
（第１の手順）
　粒子群
Figure JPOXMLDOC01-appb-I000013
の中から最小のｘ座標の値をもつ粒子を
Figure JPOXMLDOC01-appb-I000014
とする（図６（Ａ））。
（第２の手順）
　参考となる線分
Figure JPOXMLDOC01-appb-M000015
 と境界粒子以外の全ての粒子
Figure JPOXMLDOC01-appb-I000016
に対して、線分
Figure JPOXMLDOC01-appb-I000017
と
Figure JPOXMLDOC01-appb-M000018
 の成す角を測り、最小となるインデックスj2を探索し、
Figure JPOXMLDOC01-appb-M000019
 とする（図６（Ｂ））。
（第３の手順）
　ｋ＝２から参考となる線分を
Figure JPOXMLDOC01-appb-M000020
 として、
Figure JPOXMLDOC01-appb-M000021
 との成す角を計算し、最小となるインデックスj2を探索し、
Figure JPOXMLDOC01-appb-M000022
 とする。
（第４の手順）
　ｋを１インクリメントし、上記第３の手順を
Figure JPOXMLDOC01-appb-M000023
 となるまで繰り返す。
　このような手順により得られた点列
Figure JPOXMLDOC01-appb-I000024
を繋ぐと、図７に示すように、凸包となる閉曲線を求めることができる。
　図８は、点群を包括する閉曲線を求める処理の流れを示すフローチャートである。
　ステップＳ８１において、粒子の入力データを取得する。このステップＳ８１が上記第１の手順に相当する。
　ステップＳ８２において、上記第２の手順である初期粒子、参照線分を決定する。そして、ステップＳ８３において、全ての粒子に対して、参照線分と粒子間との成す角を計算する。これらステップＳ８２及びステップＳ８３が上記第２の手順に相当する。
　そして、全ての境界粒子に対して、ステップＳ８４において、成す角が最小となる粒子を境界粒子とし、ステップＳ８５において、参照成分を更新する。これらステップＳ８４及びステップＳ８５が上記第３の手順及び第４の手順に相当する。
　次に、上述のようにして求めた閉曲線（凸包）上の点列
Figure JPOXMLDOC01-appb-I000025
から、第５乃至第９の手順で凹部を求める求め方を説明する。
（第５の手順）
　粒子ｉと粒子ｉ＋１間の距離
Figure JPOXMLDOC01-appb-M000026
 を計算する。
（第６の手順）
　上記粒子間の距離
Figure JPOXMLDOC01-appb-I000027
が所定の閾値（影響半径ｈの数倍程度）よりも大きい場合（長すぎる場合）、図９に示すように、粒子群の中から線分
Figure JPOXMLDOC01-appb-M000028
 に射影を落とした距離（射影が線分上に無い場合は除く）が最小となる粒子
Figure JPOXMLDOC01-appb-I000029
を新たな界面粒子の候補とする。
　（第７の手順）
Figure JPOXMLDOC01-appb-M000030
 及び
Figure JPOXMLDOC01-appb-M000031
 が両方とも
Figure JPOXMLDOC01-appb-M000032
 より小さければ、
Figure JPOXMLDOC01-appb-I000033
をｉとｉ＋１番目の界面粒子の間に挿入する。
　上述の大小関係が成立しない場合、後述の第９の手順に進む。
（第８の手順）
　線分
Figure JPOXMLDOC01-appb-M000034
 及び
Figure JPOXMLDOC01-appb-M000035
 に対して、上記第５乃至第７の手順を繰り返し、距離が閾値より小さくなるか、界面候補の粒子が見つからなくなった場合には、処理を終了する。
（第９の手順）
　上記第５乃至第８の手順を全てのiに対して実行する。
　このような手順により、図１０に示すような凹部を求めることができる。
　図１１は、凹部を求める処理の流れを示すフローチャートである。
　ステップＳ１１１において、図６乃至図８を用いて説明した凸包構成法による表面粒子データを取得し、ステップＳ１１２において、表面粒子ｉと次の表面粒子ｉ＋１のデータを取得する。
　ステップＳ１１３において、粒子ｉと粒子ｉ＋１の距離を計算する。このステップＳ１１３が上記第５の手順に相当する。
　ステップＳ１１４において、表面粒子間の距離が所定の閾値より長いか否か判断する。そして、長い場合（ステップＳ１１４：Ｙｅｓ）、ステップＳ１１５において、他の粒子の中から候補ｋを選択し、ｋを上記ｉ＋１として上記ステップＳ１１３以降の処理を繰り返す。これらステップＳ１１４（Ｙｅｓ）及びステップＳ１１５からステップＳ１１３へ戻る流れが上記第６の手順乃至第８の手順に相当する。
　そして、ステップＳ１１２乃至ステップＳ１１４を全ての境界粒子（ｉ）に対して実行する。これが上記第９の手順に相当する。
（第２の実施の形態）
　第２の実施の形態として、表面張力項の計算法について説明する。
　本第２の実施の形態における表面張力項の計算では、上述した第１の実施の形態により得られた閉曲面から表面張力を計算する。
　上述の第１の実施の形態により、表面の粒子の位置として点列
Figure JPOXMLDOC01-appb-I000036
が得られており、２次元の場合、図１２に示すように、表面はｉ、ｉ＋１の粒子間で構成されているとする。３次元の場合は３点間の三角形要素などで表面は構成されているものとする。
　流体の表面エネルギーは一般に以下のように定義できる。
Figure JPOXMLDOC01-appb-M000037
 　ここで、上記（式４）の右辺の積分領域は、流体の境界面全域であり、
Figure JPOXMLDOC01-appb-I000038
は流体が空気に触れている部分では１、それ以外の部分では０の値をとる。
Figure JPOXMLDOC01-appb-I000039
はそれぞれ、図１３に示すように、空気に触れている部分表面張力係数、固体に触れている部分の表面張力係数である。ここで、上記（式４）をどのように粒子法で表現するかが本発明のポイントとなる。
　本発明では、界面抽出法により得られた点列
Figure JPOXMLDOC01-appb-I000040
を用いて、以下のように表面エネルギーを求める。
Figure JPOXMLDOC01-appb-M000041
Figure JPOXMLDOC01-appb-M000042
 　ここで、
Figure JPOXMLDOC01-appb-I000043
は、図１４に示すように、（式４）の
Figure JPOXMLDOC01-appb-I000044
を幅
Figure JPOXMLDOC01-appb-I000045
の間で、０から１までの曲線で滑らかに近似した関数であり、
Figure JPOXMLDOC01-appb-I000046
は、３次元上の点
Figure JPOXMLDOC01-appb-I000047
と固体との距離を表す関数である。例えば、ｙ＝０の平面に固体壁がある場合は
Figure JPOXMLDOC01-appb-M000048
 である。
Figure JPOXMLDOC01-appb-I000049
は
Figure JPOXMLDOC01-appb-M000050
 のときは０、
Figure JPOXMLDOC01-appb-M000051
 のときは１の値を取り、
Figure JPOXMLDOC01-appb-M000052
 の範囲は何らかの補間（線形補間や３次関数による補間等）を行った関数形を用いる。
　積分範囲
Figure JPOXMLDOC01-appb-I000053
は、図１５に示すように、線分
Figure JPOXMLDOC01-appb-M000054
 の上で取られる。積分計算は、ガウス積分法などを用いて数値積分することができる。
　（式５）の右辺の第１項の和（ｊ）は、粒子間（ｊ，ｊ＋１の粒子間）で取られ、第２項の和（ｋ）は、全ての粒子で取るものである。右辺第１項は、（式４）の表面エネルギーを離散化したものに対して、右辺第２項は、粒子が壁を通り抜けないという境界条件を表したポテンシャルエネルギーであり、（式６）のような形式を用いると、粒子は固体との距離が０のときポテンシャルは無限大に発散するため、固体の壁を乗り越えて行く非物理的な挙動を抑えることができる。
　（式５）のエネルギーを用いて、表面の粒子
Figure JPOXMLDOC01-appb-I000055
にかかる力は以下のように計算できる。
Figure JPOXMLDOC01-appb-M000056
 なお、（式７）は差分近似して求めても良い。
　上記（式５）は接触している面が、固体か気体かによって、表面張力係数が変わっている。図１６に示すように、この差によって、接触面が変わる際の接触角が表現できる形式となっている。
Figure JPOXMLDOC01-appb-M000057
 の場合は固体と接触している場合の方がエネルギー的に低くなるため、固体に付着しようとする力が発生することになる。粒子法では、固体を通り抜けないという条件を壁近傍のポテンシャル（式６)として表すため、粒子が固体に完全に接触することはない。そのため、本第２の実施の形態では、図１７に示すように、（式５）のようにエネルギーを設定して、固体からの距離が
Figure JPOXMLDOC01-appb-I000058
以下の場合は反発力を受け、固体からの距離が
Figure JPOXMLDOC01-appb-I000059
以上
Figure JPOXMLDOC01-appb-M000060
 未満の場合は付着力が働くようにした。
　図１８は、表面張力項の計算の流れを示すフローチャートである。
　まず、ステップＳ１８１において、シミュレーションの対象とする流体のデータを取得し、ステップＳ１８２において、現在の速度を用いて流体粒子の位置を時間ｄｔ／２分だけ更新する。そして、ステップＳ１８３において、物理モデルにより流体粒子にかかる力を算出する。
　次に、ステップＳ１８４において、上述の第１の実施の形態により流体の表面粒子のデータを取得する。
　そして、ステップＳ１８５において、全ての境界粒子に対して、上述のようにして表面エネルギーを求め、表面粒子にかかる力を計算する。
　全ての境界粒子に対する計算が終了したら、ステップＳ１８６において、流体粒子の速度を更新し、ステップＳ１８７において、更新した速度を用いて流体粒子の位置を時間ｄｔ／２分だけ更新する。
　そして、上記ステップＳ１８２乃至ステップＳ１８７を繰り返し実行する。
　さらに、上記エネルギーに代え、表面粒子の不均一性と流体の体積変化を抑えるため以下の体積エネルギー及び接点エネルギーを導入することも可能である。
体積エネルギー：
Figure JPOXMLDOC01-appb-M000061
接点エネルギー：
Figure JPOXMLDOC01-appb-M000062
　ここで、
Figure JPOXMLDOC01-appb-I000063
は流体の体積であり、
Figure JPOXMLDOC01-appb-M000064
 と計算できる。ここで
Figure JPOXMLDOC01-appb-I000065
は外向き法線単位ベクトルであり、
Figure JPOXMLDOC01-appb-I000066
は定数であり、
Figure JPOXMLDOC01-appb-I000067
は表面粒子ｓ，ｉの質量である。
（式８）、（式９）のエネルギーによる補正項は以下のように計算できる。
Figure JPOXMLDOC01-appb-M000068
Figure JPOXMLDOC01-appb-I000069
は初期体積で、定数である。
　体積エネルギーによる補正項は、図１９のように体積を
Figure JPOXMLDOC01-appb-I000070
に戻す効果があり、接点エネルギーによる補正項は、図２０のように表面粒子の不均一性を抑える効果がある。
（第３の実施の形態）
　第３の実施の形態として、ＳＰＨ法による非圧縮性粘性流体の運動を計算する際の表面張力導入方法について説明する。
　図２１は、２次元流体の例を示す図である。
　図２１に示すような状況下（ｙ＝０の平面の上におかれた流体、ｙ方向は鉛直上向きとする）における、非圧縮性粘性流体の運動方程式を考える。
Figure JPOXMLDOC01-appb-M000071
Figure JPOXMLDOC01-appb-M000072
Figure JPOXMLDOC01-appb-M000073
Figure JPOXMLDOC01-appb-I000074
は、それぞれ、流体の密度場、速度場、圧力場、音速である。
　（式１１）は質量保存則、（式１２）は運動量保存則、（式１３）は状態方程式である。
Figure JPOXMLDOC01-appb-I000075
は粘性応力テンソルで、流体の粘性係数を
Figure JPOXMLDOC01-appb-I000076
（定数）とすると、
Figure JPOXMLDOC01-appb-M000077
 である。
　（式１１）は流体が集まってくるような速度場がある場合には密度は上昇し、逆に流体が離れていくような速度場の場合には密度は低下する効果を表す。（式１２）の右辺第１項は圧力勾配項で、流体が圧力の大きい部分から圧力の小さい部分へ向かって力が発生する効果を表す。右辺第２項は粘性応力項で、流れにブレーキがかかるような効果をあらわす。右辺第３項は重力項である。
Figure JPOXMLDOC01-appb-I000078
は重力加速度、
Figure JPOXMLDOC01-appb-I000079
はｙ方向の単位ベクトルである。
　本第３の実施の形態では、上記（式１３）のような密度と圧力のみの関係を用いているが、一般の温度や内部エネルギー、エントロピー等を用いた状態方程式を用いても良い。
　また、空気との接触部には
Figure JPOXMLDOC01-appb-I000080
の表面張力係数（定数）、固体との接触部には
Figure JPOXMLDOC01-appb-I000081
の表面張力係数（定数）がかかっているものとする。
　（式１１）乃至（式１３）に示したＳＰＨ法を用いて離散化し、上述の第２の実施の形態による表面張力項を導入すると、以下のようになる。
Figure JPOXMLDOC01-appb-M000082
Figure JPOXMLDOC01-appb-M000083
Figure JPOXMLDOC01-appb-M000084
Figure JPOXMLDOC01-appb-M000085
　ここで、添え字は粒子の番号を表す。すなわちa番目の粒子の位置ベクトル、速度ベクトル、密度、圧力はそれぞれ
Figure JPOXMLDOC01-appb-I000086
である。
Figure JPOXMLDOC01-appb-I000087
はｂ番目の粒子の質量である。
Figure JPOXMLDOC01-appb-I000088
は粘性項を計算するために導入したパラメータで定数である。
　（式１５）は、上記（式１２）の運動方程式を粒子法により離散化したもので、右辺第２項は圧力勾配項、第３項は粘性応力項を表したものである。ここで、
Figure JPOXMLDOC01-appb-I000089
は
Figure JPOXMLDOC01-appb-I000090
のy方向成分、
Figure JPOXMLDOC01-appb-I000091
は表面エネルギーであり、上記（式５）と同様の形式となっており、具体的に離散化した以下の形式を採用する。
Figure JPOXMLDOC01-appb-M000092
　ここで
Figure JPOXMLDOC01-appb-I000093
はｙ＝０にある壁からの反発力を表すポテンシャルエネルギーであり、
Figure JPOXMLDOC01-appb-M000094
 を用いる。
Figure JPOXMLDOC01-appb-I000095
は定数である。
Figure JPOXMLDOC01-appb-I000096
は上述の第１の実施の形態により抽出された表面粒子であり、粒子群
Figure JPOXMLDOC01-appb-I000097
の何れか１つと対応している。
　上記（式１８）の微分
Figure JPOXMLDOC01-appb-M000098
 は以下のように計算する。
　仮に、
Figure JPOXMLDOC01-appb-I000099
が表面粒子
Figure JPOXMLDOC01-appb-I000100
と対応していない場合は、
Figure JPOXMLDOC01-appb-M000101
Figure JPOXMLDOC01-appb-M000102
 となり、仮に、
Figure JPOXMLDOC01-appb-I000103
が表面粒子
Figure JPOXMLDOC01-appb-I000104
と対応している場合は、
Figure JPOXMLDOC01-appb-M000105
Figure JPOXMLDOC01-appb-M000106
 を計算し、
Figure JPOXMLDOC01-appb-M000107
 が求められる。
　接点エネルギーは、２次元の場合以下のようになる。
Figure JPOXMLDOC01-appb-M000108
Figure JPOXMLDOC01-appb-M000109
 ここで
Figure JPOXMLDOC01-appb-I000110
，
Figure JPOXMLDOC01-appb-I000111
はそれぞれｘ方向単位ベクトル、ｙ方向単位ベクトルである。
　また、流体の体積
Figure JPOXMLDOC01-appb-I000112
は
Figure JPOXMLDOC01-appb-M000113
 と計算できる。ここで
Figure JPOXMLDOC01-appb-I000114
はベクトルを
Figure JPOXMLDOC01-appb-I000115
回転させる行列である。
　また、体積エネルギー
Figure JPOXMLDOC01-appb-I000116
は（式８）のように計算でき、
Figure JPOXMLDOC01-appb-M000117
 となる。
　さらに、
Figure JPOXMLDOC01-appb-I000118
は２次元の再標準化行列であり、
Figure JPOXMLDOC01-appb-M000119
Figure JPOXMLDOC01-appb-M000120
 である。
ここで、
Figure JPOXMLDOC01-appb-M000121
 、
Figure JPOXMLDOC01-appb-M000122
 、
Figure JPOXMLDOC01-appb-M000123
 、
Figure JPOXMLDOC01-appb-M000124
 である。
Figure JPOXMLDOC01-appb-I000125
は粒子ａ、ｂ間での１次元リーマン問題を解いて求めた時空間の中間値である。具体的には以下のように決定する。
　粒子ａ、ｂに関して、以下の特性量を定める。
Figure JPOXMLDOC01-appb-M000126
Figure JPOXMLDOC01-appb-M000127
Figure JPOXMLDOC01-appb-M000128
Figure JPOXMLDOC01-appb-M000129
　さらに、各勾配を以下のよう計算する。
Figure JPOXMLDOC01-appb-M000130
Figure JPOXMLDOC01-appb-M000131
Figure JPOXMLDOC01-appb-M000132
Figure JPOXMLDOC01-appb-M000133
Figure JPOXMLDOC01-appb-M000134
Figure JPOXMLDOC01-appb-M000135
　ここで、上記（式２７）の速度（ｖ）の上付き添え字は、各成分を表す。
これらを用いて、
Figure JPOXMLDOC01-appb-I000136
は以下のように決定する。
Figure JPOXMLDOC01-appb-M000137
Figure JPOXMLDOC01-appb-M000138
Figure JPOXMLDOC01-appb-M000139
Figure JPOXMLDOC01-appb-M000140
Figure JPOXMLDOC01-appb-M000141
　図２２及び図２３は、上記の手法により表面張力を加えて固体面との濡れ性を表現したものである。図２２は、
Figure JPOXMLDOC01-appb-M000142
 で、
Figure JPOXMLDOC01-appb-M000143
 のときを示し、
Figure JPOXMLDOC01-appb-M000144
 で、
Figure JPOXMLDOC01-appb-M000145
 のときを示す。
　これらの図２２及び図２３は、静止しときの流体の形状である。それぞれ、表面張力係数を変化させることで、接触角の違いを表現することができている。
　また、図２４は、（式１９）及び（式２０）により計算される補正項を入れなかった場合を示し、図２５は、補正項を入れた場合を示す。
　図２４では不自然な穴ができ、計算が不安定になるのに対して、図２５では安定に計算できており、補正項の効果が現れている。
　シミュレーションのパラメータは以下の通りである。
　なお、鉛直下向き方向をｙ軸とし、ｙ＝０を壁面とした。重力加速度は９．８
Figure JPOXMLDOC01-appb-I000146
、粒子数は９００個、初期形状は長さ０．００１６７［ｍ］の格子上に粒子を置き正方形にした。
　ここで、粘性係数
Figure JPOXMLDOC01-appb-I000147
は０．０１
Figure JPOXMLDOC01-appb-I000148
、
Figure JPOXMLDOC01-appb-M000149
 、影響半径ｈは０．０００３３３［ｍ］、
Figure JPOXMLDOC01-appb-M000150
 、密度は１０００
Figure JPOXMLDOC01-appb-I000151
、表面張力は
Figure JPOXMLDOC01-appb-M000152
 、音速ｃは２３．７８
Figure JPOXMLDOC01-appb-I000153
、補正項のパラメータは
Figure JPOXMLDOC01-appb-M000154
 とした。
　上述したように、本発明は粒子法の計算で、表面張力が効果的な流体運動の計算に対して適用可能である。特に溶けた金属の流し込みや、樹脂の流し込みのシミュレーション等には効果的である。
　図５のシミュレーション装置は、例えば、図２６に示すような情報処理装置（コンピュータ）を用いて実現することが可能である。図２６の情報処理装置は、ＣＰＵ（Central Processing Unit）２６０１、メモリ２６０２、入力装置２６０３、出力装置２６０４、外部記憶装置２６０５、媒体駆動装置２６０６及びネットワーク接続装置２６０７を備える。これらはバス２６０８により互いに接続されている。
　メモリ２６０２は、例えば、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フラッシュメモリ等の半導体メモリであり、シミュレーション処理に用いられるプログラム及びデータを格納する。例えば、ＣＰＵ２６０１は、メモリ２６０２を利用してプログラムを実行することにより、上述のシミュレーション処理を行う。メモリ２６０２は、図５の格納部５０２としても使用できる。
　入力装置２６０３は、例えば、キーボード、ポインティングデバイス等であり、オペレータからの指示や情報の入力に用いられる。出力装置２６０４は、例えば、表示装置、プリンタ、スピーカ等であり、オペレータへの問い合わせや処理結果の出力に用いられる。出力装置２６０４は、図５の出力部５０３としても使用できる。
　外部記憶装置２６０５は、例えば、磁気ディスク装置、光ディスク装置、光磁気ディスク装置、テープ装置等である。この外部記憶装置２６０５には、ハードディスクドライブも含まれる。情報処理装置は、この外部記憶装置２６０５にプログラム及びデータを格納しておき、それらをメモリ２６０２にロードして使用することができる。
　媒体駆動装置２６０６は、可搬型記録媒体２６０９を駆動し、その記録内容にアクセスする。可搬型記録媒体２６０９は、メモリデバイス、フレキシブルディスク、光ディスク、光磁気ディスク等である。この可搬型記録媒体２６０９には、Compact Disk Read Only
 Memory （ＣＤ－ＲＯＭ）、Digital Versatile Disk（ＤＶＤ）、Universal Serial Bus（ＵＳＢ）メモリ等も含まれる。オペレータは、この可搬型記録媒体２６０９にプログラム及びデータを格納しておき、それらをメモリ２６０２にロードして使用することができる。
　このように、シミュレーション処理に用いられるプログラム及びデータを格納するコンピュータ読み取り可能な記録媒体には、メモリ２６０２、外部記憶装置２６０５、及び可搬型記録媒体２６０９のような、物理的な（非一時的な）記録媒体が含まれる。
　ネットワーク接続装置２６０７は、通信ネットワーク２６１０に接続され、通信に伴うデータ変換を行う通信インタフェースである。情報処理装置は、プログラム及びデータを外部の装置からネットワーク接続装置２６０７を介して受け取り、それらをメモリ２６０２にロードして使用することができる。ネットワーク接続装置２６０７は、図５の出力部５０３としても使用できる。
　開示の実施形態とその利点について詳しく説明したが、当業者は、特許請求の範囲に明確に記載した本発明の範囲から逸脱することなく、様々な変更、追加、省略をすることができるであろう。

Claims (7)

  1. 　流体の表面張力をシミュレーションするシミュレーションプログラムにおいて、
    　コンピュータに、
    　入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出させ、
    　前記算出された界面の表面エネルギーを算出させ、
    　前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出させ、
    　前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力させることを特徴とするシミュレーションプログラム。
  2. 　前記界面の表面エネルギーの算出は、
    　空気との界面の表面エネルギーを算出する第１項と前記空気以外との界面の表面エネルギーを算出する第２項との和で表現した計算式を用いることを特徴とする請求項１記載のシミュレーションプログラム。
  3. 　前記空気との界面の表面エネルギーの算出は、
    　前記第１項の表面張力係数と前記第２項の表面張力係数とを変えることを特徴とする請求項２記載のシミュレーションプログラム。
  4. 　前記空気との界面の表面エネルギーの算出は、前記粒子の表面分布の不均一性を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
  5. 　前記空気との界面の表面エネルギーの算出は、前記流体の体積を補正するためのエネルギーを加算することを特徴とする請求項３記載のシミュレーションプログラム。
  6. 　流体の表面張力をシミュレーションするシミュレーション方法において、
    　コンピュータが、
    　入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出し、
    　前記算出された界面の表面エネルギーを算出し、
    　前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出し、
    　前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力することを特徴とするシミュレーション方法。
  7. 　流体の表面張力をシミュレーションするシミュレーション装置において、
    　境界条件及び初期条件を含む諸条件を入力する入力部と、
    　前記入力された境界条件及び初期条件に基づいて、前記流体を粒子の集まりで表現した流体モデルの界面を算出する界面算出部と、
    　前記算出された界面の表面エネルギーを算出する表面エネルギー算出部と、
    　前記算出された表面エネルギーに基づいて、前記界面の表面張力を算出する表面張力算出部と、
    　前記算出された表面張力に基づいて、所定時間間隔毎の前記流体の状態を出力する出力部を備えることを特徴とするシミュレーション装置。
PCT/JP2012/074295 2012-09-21 2012-09-21 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置 WO2014045416A1 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
PCT/JP2012/074295 WO2014045416A1 (ja) 2012-09-21 2012-09-21 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
JP2014536514A JP5892257B2 (ja) 2012-09-21 2012-09-21 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
US14/621,896 US10031984B2 (en) 2012-09-21 2015-02-13 Method and device for simulating surface tension

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2012/074295 WO2014045416A1 (ja) 2012-09-21 2012-09-21 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置

Related Child Applications (1)

Application Number Title Priority Date Filing Date
US14/621,896 Continuation US10031984B2 (en) 2012-09-21 2015-02-13 Method and device for simulating surface tension

Publications (1)

Publication Number Publication Date
WO2014045416A1 true WO2014045416A1 (ja) 2014-03-27

Family

ID=50340766

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2012/074295 WO2014045416A1 (ja) 2012-09-21 2012-09-21 シミュレーションプログラム、シミュレーション方法及びシミュレーション装置

Country Status (3)

Country Link
US (1) US10031984B2 (ja)
JP (1) JP5892257B2 (ja)
WO (1) WO2014045416A1 (ja)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109960841A (zh) * 2017-12-26 2019-07-02 中国科学院深圳先进技术研究院 一种流体表面张力的仿真方法、终端设备及存储介质

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10019827B2 (en) 2013-12-31 2018-07-10 Disney Enterprises, Inc. Material point method for simulation of granular materials
US10210287B2 (en) * 2013-12-31 2019-02-19 Disney Enterprises, Inc. Augmented material point method for simulating phase changes and varied materials
JP6547547B2 (ja) * 2015-09-25 2019-07-24 富士通株式会社 粒子シミュレーションプログラム、計算機資源配分方法、および粒子シミュレーション装置
US9740368B1 (en) 2016-08-10 2017-08-22 Quid, Inc. Positioning labels on graphical visualizations of graphs
US11243156B2 (en) 2020-07-08 2022-02-08 International Business Machines Corporation Determining interfacial tension for fluid-fluid-solid environments

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008111675A (ja) * 2006-10-27 2008-05-15 Mizuho Information & Research Institute Inc 流体解析装置、流体解析方法及び流体解析プログラム
JP2009129193A (ja) * 2007-11-22 2009-06-11 Univ Of Tokyo コンピュータによる流体解析方法
WO2010032656A1 (ja) * 2008-09-18 2010-03-25 国立大学法人京都大学 粒子法における界面粒子の判定方法および装置、ならびに界面粒子の判定用プログラム
JP2011123835A (ja) * 2009-12-14 2011-06-23 Nec Corp 境界粒子探索システム、境界粒子探索装置、制御方法、及びプログラム

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001156140A (ja) 1999-11-25 2001-06-08 Nec Corp 配線ｅｍシミュレーション方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008111675A (ja) * 2006-10-27 2008-05-15 Mizuho Information & Research Institute Inc 流体解析装置、流体解析方法及び流体解析プログラム
JP2009129193A (ja) * 2007-11-22 2009-06-11 Univ Of Tokyo コンピュータによる流体解析方法
WO2010032656A1 (ja) * 2008-09-18 2010-03-25 国立大学法人京都大学 粒子法における界面粒子の判定方法および装置、ならびに界面粒子の判定用プログラム
JP2011123835A (ja) * 2009-12-14 2011-06-23 Nec Corp 境界粒子探索システム、境界粒子探索装置、制御方法、及びプログラム

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
SHUNSUKE KOIZUMI: "Meniscus no Keisei o Tomonau Undo", 2007, pages 1 - 91, Retrieved from the Internet <URL:http://oacis.lib.kaiyodai.ac.jp/dspace/bitstream/123456789/614/1/kam571.pdf> [retrieved on 20121119] *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109960841A (zh) * 2017-12-26 2019-07-02 中国科学院深圳先进技术研究院 一种流体表面张力的仿真方法、终端设备及存储介质

Also Published As

Publication number Publication date
US10031984B2 (en) 2018-07-24
JPWO2014045416A1 (ja) 2016-08-18
US20150161305A1 (en) 2015-06-11
JP5892257B2 (ja) 2016-03-23

Similar Documents

Publication Publication Date Title
JP5892257B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
US9576087B2 (en) Lattice Boltzmann collision operators enforcing isotropy and Galilean invariance
JP5644872B2 (ja) シミュレーション装置、シミュレーション方法、及びプログラム
KR101328739B1 (ko) 형상 제어가 가능한 다상유체 시뮬레이션 장치 및 방법
JP6098190B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
CN109783935B (zh) 一种基于isph提高飞溅流体稳定性的实现方法
JP2012243288A (ja) シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
KR102181986B1 (ko) 더미 입자를 이용한 입자 기반의 유체 해석 시뮬레이션 방법 및 유체 해석 시뮬레이션 장치
JP6129193B2 (ja) 解析装置
JP5720551B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
JP5842992B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
WO2014045493A1 (ja) 解析方法および解析装置
JP6065616B2 (ja) シミュレーションプログラム、シミュレーション方法及びシミュレーション装置
JP6091402B2 (ja) 解析装置および解析方法
JPWO2012025995A1 (ja) 連続体運動解析プログラム、連続体運動解析方法及び連続体運動解析装置
JP6176029B2 (ja) シミュレーション装置、シミュレーションプログラム及びシミュレーション方法
JP2018185743A (ja) 流体シミュレーションプログラム、流体シミュレーション装置および流体シミュレーション方法
JP7395456B2 (ja) シミュレーション装置、及びプログラム
Kang et al. An accurate and efficient method for automatic deformation of unstructured polyhedral grids to simulate the flow induced by the motion of solid objects
JP6029457B2 (ja) 解析装置および解析方法
Mohaghegh et al. A simplified model for the normal collision of arbitrary shape particles in a viscous flow
Li An experimental and numerical study of normal particle collisions in a viscous liquid
Kumar et al. Improved Methodology for Mass Conservation in Sharp Interface Immersed Boundary Method for Moving Boundaries
Mandal et al. Haptic Rendering of Solid Object Submerged in Flowing Fluid with Environment Dependent Texture
Zhang et al. The Non-Newtonian Fluid Simulation Based on Predictive-Corrective Incompressible SPH

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 12885019

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2014536514

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 12885019

Country of ref document: EP

Kind code of ref document: A1