(57)【特許請求の範囲】
【請求項１】
粒子法を利用して、第１方向の長さに対する第２方向の長さの比が所定値以上の領域における流体の挙動を解析するための流動解析方法であって、
前記流体の第１方向の長さを、粒子法を利用した流動解析方法における１つの粒子により表現し、前記流体の第２方向の長さを当該第２方向の解像度に応じて複数の前記粒子により表現し、
前記粒子を前記流体の鉛直方向の長さを高さとする円柱により表現するとともに、前記円柱の体積を一定に維持した状態で前記円柱の直径が変化するようにし、前記流体の第１方向の長さの変化に応じて前記粒子の第１方向の長さを変化させることにより、第２方向における流体の挙動を解析することを特徴とする流動解析方法。
【請求項２】
請求項１に記載の流動解析方法において、前記第１方向の長さの変化に応じて粒子径を変化させることを特徴とする流動解析方法。
【請求項３】
請求項２に記載の流動解析方法において、粘性項を陰的に解くことにより流体の挙動を解析することを特徴とする流動解析方法。
【請求項４】
請求項３に記載の流動解析方法において、前記粒子が配置されていない前記流体の第１方向の位置に仮想的な粒子を配置し、当該仮想的な粒子と流体を表現する粒子の相互作用を考慮して流体の挙動を解析することを特徴とする流動解析方法。
【請求項５】
請求項４に記載の流動解析方法において、前記第１方向における前記流体の粒子と仮想的な粒子のレイノルズ数が一致するように動粘性係数を補正することを特徴とする流動解析方法。
【請求項６】
粒子法を利用して、第１方向の長さに対する第２方向の長さの比が所定値以上の領域における流体の挙動を解析するための流動解析装置であって、
前記流体の第１方向の長さを、粒子法を利用した流動解析方法における１つの粒子により表現し、前記流体の第２方向の長さを当該第２方向の解像度に応じて複数の前記粒子により表現する初期データ設定手段と、
前記粒子を前記流体の鉛直方向の長さを高さとする円柱により表現するとともに、前記円柱の体積を一定に維持した状態で前記円柱の直径が変化するようにし、前記流体の第１方向の長さの変化に応じて前記粒子の第１方向の長さを変化させるデータ更新手段とを備えることを特徴とする流動解析装置。
【請求項７】
粒子法を利用して、第１方向の長さに対する第２方向の長さの比が所定値以上の領域における流体の挙動を解析するための流動解析プログラムであって、
前記流体の第１方向の長さを、粒子法を利用した流動解析方法における１つの粒子により表現する処理と、前記流体の第２方向の長さを当該第２方向の解像度に応じて複数の前記粒子により表現する処理と、前記粒子を前記流体の鉛直方向の長さを高さとする円柱により表現するとともに、前記円柱の体積を一定に維持した状態で前記円柱の直径が変化するようにし、前記流体の第１方向の長さの変化に応じて前記粒子の第１方向の長さを変化させる処理とをコンピュータに実行させることを特徴とする流動解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、非常に薄い平面領域等のアスペクト比が極端な領域における流体の挙動解析に適用して好適な、流動解析方法、流動解析装置、及び流動解析プログラムに関する。
【背景技術】
【０００２】
従来より、流動解析方法の一つとして、均一な大きさの複数の粒子により流体を表現し、複数の粒子の挙動をコンピュータを用いて計算することにより流体の挙動を解析する方法、いわゆる粒子法を利用した流動解析方法が知られている（特許文献１−３参照）。この粒子法を利用した流動解析方法は、空気や水等の自由表面を有する物体、特に大きく変形したり、分裂，合体したりする物体の挙動解析に適しており、原子力分野を始めとしてプラント，航空，船舶，環境等の様々な技術分野において広く利用されている。
【特許文献１】特開平７−３３４４８４号公報
【特許文献２】特開２０００−３３９２９２号公報
【特許文献３】特開２００２−１３７２７２号公報
【発明の開示】
【発明が解決しようとする課題】
【０００３】
従来の粒子法を利用した流動解析方法では、解析領域及び解析時間内において粒子の大きさ、換言すれば空間解像度を変化させることはできない。このため、例えば水平方向の長さに対し鉛直方向の長さが非常に短い領域等、アスペクト比が極端な領域における流体の挙動を解析する場合、長さが短い方向の空間解像度を考慮して粒子の大きさを設定し、設定した大きさの粒子を解析領域内に配置しなければならない。このような背景から、従来の粒子法を利用した流動解析方法によりアスペクト比が極端な領域における流体の挙動を解析する場合には、解析領域内に非常に多くの粒子（計算点）が配置されるために、コンピュータの処理能力の制約上、現実的な時間で流動解析を行うことが困難になる。
【０００４】
本発明は、このような課題を解決するためになされたものであり、その目的は、アスペクト比が極端な領域における流体の挙動解析に要する時間を短縮可能な流動解析方法、流動解析装置、及び流動解析プログラムを提供することにある。
【課題を解決するための手段】
【０００５】
本発明に係る流動解析方法、流動解析装置、及び流動解析プログラムは、粒子法を利用して、第１方向の長さに対する第２方向の長さの比が所定値以上の領域における流体の挙動を解析する際、流体の第１方向の長さを、粒子法を利用した流動解析方法における１つの粒子により表現し、流体の第２方向の長さを第２方向の解像度に応じて複数の粒子により表現し、粒子を流体の鉛直方向の長さを高さとする円柱により表現するとともに、円柱の体積を一定に維持した状態で円柱の直径が変化するようにし、流体の第１方向の長さの変化に応じて粒子の第１方向の長さを変化させることにより、第２方向における流体の挙動を解析する。
【発明の効果】
【０００６】
本発明に係る流動解析方法、流動解析装置、及び流動解析プログラムによれば、解析領域内に配置される要素数を大幅に削減することができるので、アスペクト比が極端な領域における流体の挙動解析に要する時間を大幅に短縮することができる。
【発明を実施するための最良の形態】
【０００７】
以下、図面を参照して、本発明の実施形態となる流動解析装置の構成及びその動作について説明する。
【０００８】
〔流動解析装置の構成〕
本発明の実施形態となる流動解析装置１は、パーソナルコンピュータ，ワークステーション，汎用コンピュータ等の公知の情報処理装置により構成され、図１に示すように、ＲＡＭ（Random Access Memory）２，ＲＯＭ（Read Only Memory）３，初期設定データベース（ＤＢ）４，及びＣＰＵ（Central Processing Unit）５を主な構成要素として備える。ＲＡＭ２は、ＣＰＵ５が実行する流動解析装置１の起動プログラムや流動解析プログラム等のコンピュータプログラム及び各種データを一時的に格納するワーキングエリアを提供する。ＲＯＭ３は、コンピュータプログラム及びこのコンピュータプログラムを実行するために必要な各種データを記憶する。なおＲＯＭ３は、磁気的又は光学的な記憶媒体若しくは半導体メモリ等のＣＰＵ５により読み取り可能な記憶媒体を含む構成を有し、ＲＯＭ３に記憶されているコンピュータプログラムや各種データの一部又は全部は電気通信回線を介してダウンロードされるように構成してもよい。
【０００９】
初期設定ＤＢ４は、粒子法により流動解析を実行するために必要な各種パラメータを格納する。具体的には、このパラメータには、粒子法における粒子ｉの粒子径ｒiの初期値ｒiini，位置（重心の座標），速度，温度，圧力，流体の物性（例えば粘度，密度，熱伝導率，比熱、潜熱等），時間刻み幅Δｔ，計算終了時間Ｔ等の初期条件に関する情報が含まれる。ＣＰＵ５は、ＲＯＭ３に記憶されているコンピュータプログラムをＲＡＭ２内にロードし、読み出されたコンピュータプログラムを実行することにより、流動解析装置１全体の処理動作を制御する。具体的には、ＣＰＵ５は、ＲＯＭ３に記憶された流動解析プログラム及び流動解析プログラムを実行するために必要なデータをＲＡＭ２へロードし、流動解析プログラムに従って粒子法を利用した流動解析を実行する。本実施形態では、ＣＰＵ５は、ＲＯＭ３に記憶されている流動解析プログラムを実行することにより、外力項計算部１１，粘性項計算部１２，粒子移動追跡部１３，圧力項計算部１４，粒子データ更新部１５，及び終了判定部１６として機能する。各部の機能の詳細については後述する。
【００１０】
〔流動解析方法〕
このような構成を有する流動解析装置１は、以下に示す流動解析処理を実行することによりアスペクト比が極端な領域における流体の挙動を現実的な時間で解析する。以下、図２に示すフローチャートを参照して、この流動解析処理を実行する際の流動解析装置１の動作について詳しく説明する。なお以下の説明では、アスペクト比が極端な領域における流体の挙動として、水平方向の長さに対して高さ方向の長さが非常に小さい領域における流体の挙動、具体的には水平方向の長さと高さ方向の長さの比が１００対１以上の領域における流体の挙動を想定する。このような流体の挙動の具体例としては、２枚のガラス板の間に数点塗布された樹脂材料の２枚のガラス板を押し狭めた際における挙動を例示できる。樹脂材料の挙動を解析することにより、樹脂材料の塗布位置及び塗布量を最適化することができる。
【００１１】
図２に示すフローチャートは、流動解析装置１に対しオペレータが流動解析処理の開始を指示したタイミングで開始となり、流動解析処理はステップＳ１の処理に進む。
【００１２】
ステップＳ１の処理では、ＣＰＵ５が、初期設定ＤＢ４からＲＡＭ２内に解析対象となる流体の物性，時間刻み幅ΔＴ，計算終了時間Ｔ等の計算条件パラメータ、及び粒子の座標位置，速度，温度等の初期条件パラメータを読み出し、読み出されたパラメータに基づいて流体を複数の粒子により表現する。従来の粒子法を利用した流動解析手法では、図３に示すように高さ方向の解像度に合わせた大きさの粒子を解析領域全体に配置するが、本流動解析処理では、図４に示すように高さ方向の粒子を平均化して１要素（１粒子）で表現するように初期条件パラメータを設定する。
【００１３】
より具体的には、本流動解析処理では、図４に示すように各粒子は円柱形状を有する。このような粒子形状によれば、円柱の直径、すなわち粒子径を水平方向の解像度に合わせて調整することにより、水平方向の長さに対して高さ方向の長さが非常に小さい解析領域内に配置される粒子数を大幅に少なくし、現実的な時間で流動解析を行うことができる。つまり本流動解析処理では、高さ方向の流体の動きを平均化し、水平方向の流体の動きを重点的に解析する。
【００１４】
なお高さ方向を単に１要素で表現しただけでは、高さ方向の流体の動きを完全に無視して水平方向の流体の動きのみを解析する単純な２次元計算と同じになり、高さ方向の変化を考慮することができなくなる。そこで本流動解析処理では、図５（ａ），（ｂ）に示すように、高さの変化に応じて円柱の体積を一定に維持した状態で円柱の直径が変化するように初期条件パラメータを設定する。すなわち初期状態における粒子の半径及び高さをそれぞれｒ０，ｈ０とした時、任意の時刻における粒子の半径ｒ及び高さｈが以下の数式１を満たすように初期条件パラメータを設定する。これにより、ステップＳ１の処理は完了し、流動解析処理はステップＳ２の処理に進む。
【数１】
0005242224.tif 000002
【００１５】
ステップＳ２の処理では、ＣＰＵ５が、高さ方向の座標変化に応じて前述の数式１を満足するように円柱の半径ｒを更新することにより、各粒子の水平方向の位置変化を計算する。なおこのステップＳ２の処理は１回目の計算では省略される。これにより、ステップＳ２の処理は完了し、流動解析処理はステップＳ３の処理に進む。
【００１６】
ステップＳ３の処理では、外力項計算部１１が、重力等の各粒子に直接的に作用する外力項を計算する。粒子法における外力項の計算方法は本願発明の出願時点で既に公知であるので詳細な説明は省略する。これにより、ステップＳ３の処理は完了し、流動解析処理はステップＳ４の処理に進む。
【００１７】
ステップＳ４の処理では、粘性項計算部１２が、粘性項を計算し、粒子ｉが有する速度分布を近傍の粒子ｊに分配することにより、粘性による粒子の速度の拡散を計算する。具体的には、粘性項は以下の数式２により表される。数式２中、パラメータｕ，ｔ，υはそれぞれ速度ベクトル，計算時刻，動粘性係数を示す。この数式２を粒子法で離散化すると、粘性項は陽的には数式３，陰的には数式４のように表される。数式２，３中、パラメータｄ，λ，ｎ０，ｗｉｊｋはそれぞれ次元数，拡散を解析解に一致させるための係数，基準粒子密度，粒子ｉ，ｊ間の重み，及びタイムステップを示す。
【数２】
0005242224.tif 000003
【数３】
0005242224.tif 000004
【数４】
0005242224.tif 000005
【００１８】
粒子ｉの速度変化Δｕは、以下の数式５により表されるので、これを数式４に代入することにより、以下の数式６に示す連立一次方程式が導かれる。本実施形態では、この数式６を解くことにより粘性項を陰的に解く。このように粘性項を陰的に解くことにより、タイムステップｋを大きくすることができるので、発散することなく計算を安定的に行うことができる。
【数５】
0005242224.tif 000006
【数６】
0005242224.tif 000007
【００１９】
なお既述の通り、本流動解析処理は、図６に示すように粒子の形状を薄い円柱形状と見立てて解析を行うものであるので、図７に示すように高さ方向に仮想的な粒子Ｐを配置して計算する場合、粒子の直径分の隙間がある領域を計算することになってしまう。すなわち、粒子の直径をＬとした場合、高さ方向に仮想的な粒子Ｐを配置したとしても幅Ｌの領域を解析することになるが、実際の幅は幅Ｌよりも小さいために物理的に正確な解析を行うことはできない。そこで本流動解析処理では、幅Ｌの解析対象と実際の現象のレイノルズ数を合わせることにより物理的に正確な計算を行うことを可能にする。
【００２０】
具体的には、実際の現象のレイノルズ数Ｒｅは以下の数式７により表される。数式７中パラメータＵは特性速度を示す。従って実際のレイノルズ数Ｒｅｒ，解析解のレイノルズ数Ｒｅｓをそれぞれ以下の数式８，９のように表すと、これらが等しい場合、以下の数式１０が導出される。従って本流動解析処理では、動粘性係数υを以下の数式１０のように補正することにより、高さ方向の粘性を物理的に正確に計算することができる。これにより、このステップＳ４の処理は完了し、流動解析処理はステップＳ５の処理に進む。
【数７】
0005242224.tif 000008
【数８】
0005242224.tif 000009
【数９】
0005242224.tif 000010
【数１０】
0005242224.tif 000011
【００２１】
ステップＳ５の処理では、粒子移動追跡部１３が、ステップＳ３及びステップＳ４の計算結果に基づいて粒子の速度を仮更新する。これにより、ステップＳ５の処理は完了し、流動解析処理はステップＳ６の処理に進む。
【００２２】
ステップＳ６の処理では、粒子移動追跡部１３が、ステップＳ５の処理により仮更新された速度による時間刻み幅ΔＴの間での粒子の移動を追跡し、追跡結果に基づいて粒子の座標位置を仮更新する。これにより、ステップＳ６の処理は完了し、流動解析処理はステップＳ７の処理に進む。
【００２３】
ステップＳ７の処理では、圧力項計算部１４が、ステップＳ６の処理により算出された座標位置における各粒子の粒子数密度を算出し、算出された粒子数密度が所定値であるか否かを判別する。判別の結果、粒子数密度が所定値でない場合、圧力項計算部１４は、粒子数密度が所定値になるように粒子の圧力を修正することにより粒子の座標位置及び流速を修正する。そして圧力計算部１４は、修正された圧力によって生じる圧力項を算出する。なお粒子法における圧力項の計算方法は本願発明の出願時点で既に公知であるので詳細な説明は省略する。これにより、ステップＳ７の処理は完了し、流動解析処理はステップＳ８の処理に進む。
【００２４】
ステップＳ８の処理では、粒子データ更新部１５が、ステップＳ７の処理の算出結果に基づいてＲＡＭ２内に格納されている粒子の速度を更新する。これにより、ステップＳ８の処理は完了し、流動解析処理はステップＳ９の処理に進む。
【００２５】
ステップＳ９の処理では、粒子データ更新部１５が、ステップＳ７の処理の算出結果に基づいてＲＡＭ２内に格納されている粒子の座標位置を更新する。これにより、ステップＳ９の処理は完了し、流動解析処理はステップＳ１０の処理に進む。
【００２６】
ステップＳ１０の処理では、終了判定部１６が、解析時間ｔを時間刻み幅ΔＴだけインクリメントする。これにより、ステップＳ１０の処理は完了し、流動解析処理はステップＳ１１の処理に進む。
【００２７】
ステップＳ１１の処理では、終了判定部１６が、解析時間ｔが計算終了時刻Ｔであるか否かを判別する。判別の結果、解析時間ｔが計算終了時刻Ｔである場合、終了判定部１０は、ＲＡＭ２内に格納されている粒子データをＲＯＭ３の所定ファイルにコピーし、一連の流動解析処理を終了する。この一連の流動解析処理によれば、例えば図８（ａ）〜（ｅ）に示すように、２枚のガラス板の間に数点塗布された樹脂材料の２枚のガラス板を押し狭めた際における挙動を時系列で解析することができる。一方、解析時間ｔが計算終了時刻Ｔでない場合には、終了判定部１０は流動解析処理をステップＳ２の処理に戻す。
【００２８】
以上の説明から明らかなように、本発明の実施形態となる流動解析装置１は、高さ方向の長さに対する水平方向の長さの比が所定値以上の領域における流体の挙動を解析する際、流体の高さ方向の長さを１つの粒子により表現し、流体の水平方向の長さを第２方向の解像度に応じて複数の粒子により表現し、流体の高さ方向の長さの変化に応じて粒子の水平方向の長さを変化させることにより、水平方向における流体の挙動を解析する。そしてこのような流動解析装置１によれば、解析領域内に配置される粒子数を大幅に削減することができるので、アスペクト比が極端な領域における流体の挙動解析に要する時間を大幅に短縮することができる。
【００２９】
また本発明の実施形態となる流動解析装置１は、粘性項を陰的に解くことにより流体の挙動を解析するので、タイムステップｋを大きくすることにより、発散することなく計算を安定的に行うことができる。また本発明の実施形態となる流動解析装置１は、粒子の上下方向に仮想的な粒子Ｐを配置し、流体を構成する粒子と仮想的な粒子Ｐのレイノルズ数が一致するように動粘性係数を補正するので、高さ方向の粘性を物理的に正確に計算することができる。
【００３０】
以上、本発明者によってなされた発明を適用した実施の形態について説明したが、この実施形態による本発明の開示の一部をなす記述及び図面により本発明は限定されることはない。例えば、高さ方向の長さは強制的に与えることもできるが，力の釣り合いを考慮して算出してもよい。すなわち、粒子の上方向の壁に外力Ｆｅと流体の反発力Ｆｆが加わっているとすると、壁の運動方程式は以下の数式１１のように表される。数式１１中、Ｍは壁の質量，ａは壁の加速度を示す。従って、この数式１１により壁の動きを計算することにより高さ方向の長さを求めることができる。なお，流体の反発力Ｆｆは以下の数式１２により表される。数式１２中、Ｐは圧力，ｎは高さ方向の単位ベクトルを示す。
【数１１】
0005242224.tif 000012
【数１２】
0005242224.tif 000013
【００３１】
また本実施形態は本発明を粒子法を利用した流動解析方法に適用したものであるが、本発明は解析領域をメッシュにより表現する差分法や有限体積法等の粒子法以外の解析法を利用した流動解析方法にも適用することができる。具体的には本発明では、アスペクト比の極端な方向をｚ方向とするとｚ方向は1要素で表現される。その際、ｚ方向の高さをｈ、その他の方向の面積をＳとすると、要素の体積Ｖは以下の数式１３のように表される。また要素内にある質量をＭとすると、要素の密度ρは以下の数式１４のように表される。この時、高さｈの時間変化量を∂ｈ／∂ｔとすると、密度ρの時間変化量は以下の数式１５のように表される。これは高さｈが変化したために起こる密度変化である。
【数１３】
0005242224.tif 000014
【数１４】
0005242224.tif 000015
【数１５】
0005242224.tif 000016
【００３２】
本来の連続式は以下に示す数式１６で与えられる。しかしながら本発明では、ｚ方向の密度変化を特殊な扱い方で解析するので、連続式は以下の数式１７で与えられる。従って数式１７に示す連続式を採用することにより、本発明を粒子法以外の解析法を利用した流動解析方法にも適用することができる。
【数１６】
0005242224.tif 000017
【数１７】
0005242224.tif 000018
【００３３】
また本実施形態における流動解析装置１は、モデム等の通信用部材を備えてインターネットを含む電気通信回線と接続可能なように構成してもよい。この場合には、本実施形態における流動解析プログラムは、電気通信回線からダウンロード等によって、ＲＡＭ２又はＲＯＭ３内に格納させてもよい。なお、この場合における通信ネットワークからダウンロードするためのダウンロードプログラムは、予めＲＯＭ３内に格納されているか、別の記憶媒体からＲＯＭ３内にインストールされるものとする。このように、本実施形態に基づいて当業者等によりなされる他の実施の形態、実施例及び運用技術等は全て本発明の範疇に含まれる。
【産業上の利用可能性】
【００３４】
本発明は、２枚のガラス板間に塗布された樹脂材料が押し広げられる際の挙動、アスペクト比が極端な領域における流体の挙動解析に適用できる。
【図面の簡単な説明】
【００３５】
【図１】本発明の実施形態となる流動解析装置の構成を示すブロック図である。
【図２】本発明の実施形態となる流動解析処理の流れを示すフローチャート図である。
【図３】本願発明の流動解析処理における粒子の配置例を示す模式図である。
【図４】本願発明の流動解析処理における粒子の配置例を示す模式図である。
【図５】高さの変化に応じて体積を一定に維持した状態で円柱の直径が変化する様子を示す図である。
【図６】本願発明の流動解析処理における粒子の配置例を示す模式図である。
【図７】図６に示す粒子の上下方向に仮想的な粒子を配置した状態を示す図である。
【図８】２枚のガラス板の間に数点塗布された樹脂材料の２枚のガラス板を押し狭めた際における挙動を本発明により解析した例を示す図である。
【符号の説明】
【００３６】
１：流動解析装置
２：ＲＡＭ
３：ＲＯＭ
４：初期設定ＤＢ
５：ＣＰＵ
１１：外力項計算部
１２：粘性項計算部
１３：粒子移動追跡部
１４：圧力項計算部
１５：粒子データ更新部
１６：終了判定部
【図１】
0005242224.tif 000019
【図２】
0005242224.tif 000020
【図３】
0005242224.tif 000021
【図４】
0005242224.tif 000022
【図５】
0005242224.tif 000023
【図６】
0005242224.tif 000024
【図７】
0005242224.tif 000025
【図８】
0005242224.tif 000026