JP2017219378A

JP2017219378A - シミュレーション装置、シミュレーション装置の制御方法、および制御プログラム

Info

Publication number: JP2017219378A
Application number: JP2016112882A
Authority: JP
Inventors: 島田　直樹; Naoki Shimada; 直樹島田; 早紀仙田; Saki Senda
Original assignee: Sumitomo Chemical Co Ltd
Current assignee: Sumitomo Chemical Co Ltd
Priority date: 2016-06-06
Filing date: 2016-06-06
Publication date: 2017-12-14
Anticipated expiration: 2036-06-06
Also published as: JP6576303B2

Abstract

【課題】水系環境における数値シミュレーションの演算を発散させることなく高精度に実行する。【解決手段】三次元空間に設定した計算格子の少なくとも１つの座標スケールが他の座標スケールの少なくとも１つの２倍以上であり、基礎式から導出した偏微分方程式を各座標スケールに対応して分解し、さらに直交分解した方程式を用いるものであり、上記少なくとも１つの座標スケールに対応した方程式を用いて、圧力および三方向成分の速度を算出する長波長成分算出部と、上記他の座標スケールに対応した方程式を用いて、圧力および上記他の座標スケールに相当する方向成分の速度を算出する短波長成分算出部と、上記長波長成分算出部、および上記短波長成分算出部の算出結果を足し合わせる合成部と、を備える。【選択図】図１

Description

本発明は、非圧縮性流体あるいは擬圧縮性のニュートン流体の流れをシミュレーションするシミュレーション装置等に関する。

近年の生態系の保護に対する関心の高まり、人類の安全性確保に対する意識の向上に伴い、河川、湖、海洋等における水系環境の調査研究の重要性が増している。水系環境において、各種生態系に関与する原因要素を直接特定することは困難であり、水系環境に与える影響が大きい因子である、化学物質の濃度、温度、水況等の時空間変化を観測するか、または当該変化を試験にて推定することが行われている。

また、水況、すなわち化学物質を輸送する流動は、水系環境に与える影響が極めて大きいため、数日から数年に渡って多点で実測されている事例も存在する。しかし、多額の費用と労力とを必要とするため、実測データが十分に集まっているとは言い難い。

一方で、流体の数値シミュレーションをベースに、各種モデルを加えて水系環境を補完し予測する試みも行われている。

流体の数値シミュレーションとして、例えば、非特許文献１には、静水圧近似の代表的な海洋モデルが開示されている。また、非特許文献２には、流体解析モデルが開示されている。

George L. Mellor、"USERS GUIDE for A THREE-DIMENSIONAL PRIMITIVE EQUATION NUMERICAL OCEAN MODEL"2002、[online]、2002年10月、Princeton University [平成27年12月1日検索] <URL:(http://www.ccpo.odu.edu/POMWEB/UG.10-2002.pdf> 小林敏雄[ほか]編「数値流体力学ハンドブック」丸善、2003年3月、pp.19-65 棚橋隆彦著者「ＣＦＤの基礎理論」アイピーシー、pp.201-220、1990 Hirt, C.W., Nichols, B.D., Romero, N.C. (1975): SOLA-A Numerical Solution Algorithm for Transient Fluid Flow ,Los Alamos Scientific Laboratory, LA-5852, pp. 1-50

しかしながら、上述した従来技術では、河川、湖、海洋等の水系環境における数値シミュレーションの課題を解決することができない。河川、湖、海洋等の水系環境では、数値シミュレーションに用いる計算格子の鉛直方向のスケールが、水平方向のスケールと比較して１／１０〜１／１０００のオーダで小さくなる。非特許文献１の方法では鉛直方向の運動量が考慮されておらず、これが駆動力となる流れについては解の精度が低下するおそれがある。一方、非特許文献２の方法では、一般に解析法に関する安定性の限界を逸脱することから、シミュレーションの演算が収束せず、数値解が得られないおそれがある。

本発明は、前記の問題点に鑑みてなされたものであり、その目的は、水系環境における数値シミュレーションの演算を発散させることなく高精度に実行できるシミュレーション装置等を実現することにある。

上記の課題を解決するために、本発明に係るシミュレーション装置は、三次元空間を運動する非圧縮性流体あるいは擬圧縮性のニュートン流体のシミュレーション装置であって、三次元空間に設定した計算格子の少なくとも１つの座標スケールが他の座標スケールの少なくとも１つの２倍以上であり、基礎式から導出した偏微分方程式を各座標スケールに対応して分解し、さらに直交分解した方程式を用いるものであり、上記少なくとも１つの座標スケールに対応した方程式を用いて、圧力および三方向成分の速度を算出する長波長成分算出部と、上記他の座標スケールに対応した方程式を用いて、圧力および小さな座標スケールに相当する方向成分の速度を算出する短波長成分算出部と、上記長波長成分算出部、および上記短波長成分算出部の算出結果を足し合わせる合成部と、を備える構成である。

また、上記の課題を解決するために、本発明に係るシミュレーション装置の制御方法は、三次元空間を運動する非圧縮性流体あるいは擬圧縮性のニュートン流体のシミュレーション装置の制御方法であって、三次元空間に設定した計算格子の少なくとも１つの座標スケールが他の座標スケールの少なくとも１つの２倍以上であり、基礎式から導出した偏微分方程式を各座標スケールに対応して分解し、さらに直交分解した方程式を用いるものであり、上記少なくとも１つの座標スケールに対応した方程式を用いて、圧力および三方向成分の速度を算出する長波長成分算出ステップと、上記他の座標スケールに対応した方程式を用いて、圧力および上記他の座標スケールに相当する方向成分の速度を算出する短波長成分算出ステップと、上記長波長成分算出ステップ、および上記短波長成分算出ステップにおける算出結果を足し合わせる合成ステップと、を含む方法である。

上記の構成によれば、運動方程式が偏微分方程式の楕円性を持つ場合において、該偏微分方程式を各座標スケールに応じた方程式に分解し、さらに直交分解して得られた方程式を複数個用いるので、発散させることなく分解したスケールごとに、流体の計算を行うことができる。すなわち、水系環境における数値シミュレーションの演算を発散させることなく高精度に実行できる。

なお、非圧縮性流体とは、圧力に依存して密度が変化しない流体である。また、擬圧縮性のニュートン流体とは、流体の密度を圧力以外の物理因子によって変化させることができる流体を示し、例えば、温度や化学種濃度に依存して密度を変えられる流体である。

本発明に係るシミュレーション装置では、現在の時刻の液面の位置を、次の時刻の圧力境界条件とすることが好ましい。

本発明に係るシミュレーション装置では、上記少なくとも１つの座標スケールが、上記他の座標スケールの１０の２乗倍より大きいものであってもよい。

上記の構成によれば、海洋計算、静水圧近似の計算に適用することができる。

本発明の各態様に係るシミュレーション装置は、コンピュータによって実現してもよく、この場合には、コンピュータを前記シミュレーション装置が備える各部（ソフトウェア要素）として動作させることにより前記シミュレーション装置をコンピュータにて実現させるシミュレーション装置の制御プログラム、およびそれを記録したコンピュータ読取り可能な記録媒体も、本発明の範疇に入る。

本発明によれば、運動方程式が偏微分方程式の楕円性を持つ場合において、該偏微分方程式を各座標スケールに応じた方程式に分解し、さらに直交分解して得られた方程式を複数個用いるので、発散させることなく分解したスケールごとに、流体の計算を行うことができるという効果を奏する。すなわち、本発明によれば、水系環境における数値シミュレーションの演算を発散させることなく高精度に実行することができ、計算負荷の面において実用性を得ることが困難な計算も可能になる。

シミュレーション装置の要部構成を示すブロック図である。スタガード配置の格子例を示す図である。自由表面における圧力境界条件を説明するための図である。シミュレーション装置におけるシミュレーションの流れを示すフローチャートである。計算対象とした矩形三次元領域を示す図である。計算条件を示す図である。計算結果を示す図である。海洋に用いた場合の計算例を説明するための図であり、瀬戸内海の深さ情報を示す図である。海洋に用いた場合の計算例を説明するための図であり、瀬戸内海の潮流の時間と速度との関係を示す図である。海洋に用いた場合の計算例の結果を示す図である。

本実施形態に係るシミュレーション装置１は、所定の式を予め記憶させておくことで、シミュレーションを行うことができる。

そこでまず、本実施形態に係るシミュレーション装置１が実行するシミュレーション方法で用いる式の導出方法を説明する。

なお、以下では、三次元空間に設定した計算格子のうち、鉛直方向の座標スケールが、２つの水平方向の座標スケールの２倍以上である場合の式の導出方法について説明する。

三次元直交座標系における非圧縮流体の数値計算では、安定性を確保するため、以下の式（１）（Ｃｏｕｒａｎｔ条件）を満たす必要がある。

ここで、ｘおよびｙは水平方向座標を示し、ｚは鉛直方向座標を示す。また、Δｘ、Δｙ、Δｚはそれぞれの方向の座標スケールを示し、Δｔは時間刻み幅を示す。また、ｕ、ｖ、ｗは各方向の速度を示す。また、下付添字の「ｍａｘ」は計算領域内における最大値を表している。

ここで、河川、湖、海洋等の水系環境では、数値シミュレーションに用いる計算格子の鉛直方向のスケールが、水平方向のスケールと比較して１／１０〜１／１０００のオーダで小さくなる。この点に基づき、Δｘ≒Δｙ＞＞Δｚとした場合、式（１）より、Δｔはｗ_ｍａｘ／Δｚにほぼ依存することになる。なお、Δｘ≒Δｙ≒Δｚとした場合は、式（１）を満たせば数値解を得ることができるが、水平方向に膨大な数の計算格子を要することになる。

そこで、以下の静水圧近似モデルが提案されている。上記のような計算格子スケールの制約から、シミュレーションの演算が収束せず、発散してしまうという問題を回避するため、まず、圧力を静水圧近似した以下の式を用いる。

ここで、ｐは圧力、ｐ_０は大気圧、ｈは水面高さ、ρは密度、ｇは重力加速度を示す。この近似式を用いることにより、Ｃｏｕｒａｎｔ条件は次式に緩和される。

ただし、前記の式（３）では、鉛直方向の動的挙動を計算することができない。近年のシミュレーションでは、陸地近海や湾内において数百ｍ以下の分解能での精度が求められており、水平方向速度の水深依存性と鉛直方向流の影響を考慮することは不可欠である。したがって、前記の式（３）は必ずしも妥当とは言えない。

そこで、本実施形態では、上記のような計算格子スケールの制約から、シミュレーションの演算が収束せず、発散してしまうという問題を回避するため、以下の方法で方程式を導出する。具体的には、まず、鉛直方向に対し準動的な運動方程式を導出し、次に、Ｐｒｏｊｅｃｔｉｏｎ法に適用するために離散化法を適用する。そして、自由表面を組み込むように拡張する。

〔式の導出〕
シミュレーション対象となる流体は、流体として相変化がない擬圧縮性のニュートン流体を仮定する。

まず、カーテシアン三次元直交座標系に対して、以下のように、基礎式として、質量保存式（４）および運動量保存式（５）（６）（７）を与える。

ここで、ｔは時間、Ｄ／Ｄｔは物質微分、μは粘度である。Δｘ≒Δｙ＞＞Δｚを考慮すれば、ｚ方向の式（７）は次式（８）のように近似できる。

次に、鉛直方向の波長を考慮して、変数ｗおよび変数ｐをそれぞれ以下の通り長波長成分と短波長成分との２つに分解する。なお、数式中に示したアルファベット文字の上にハット（^）が付いた文字を、本明細書中では、アルファベット文字の横にハットを付して表記する。以下、他の文字についても同様である。

ここで、上記の式（４）〜式（８）におけるｕ、ｖ、およびｗのうちのｗ＾は、波長がΔｘ≒Δｙ＞＞Δｚのスケールで捕捉できる波長成分（以下、長波長成分と称することもある。）の各方向の速度を意味し、ｐのうちのｐ＾は、長波長成分の圧力を意味する。また、ｗのうちのｗ’は、Δｚのスケールで捕捉できる波長成分（以下、短波長成分と称することもある。）のｚ方向の速度を意味し、ｐのうちのｐ’は、短波長成分の圧力を意味する。式（９）、および式（１０）を用いて、上記の式（４）〜式（６）、および式（８）を書き換えると、次の式（１１）、式（１２）、式（１３）、および式（１４）がそれぞれ得られる。

次に、数値解法としてＰｒｏｊｅｃｔｉｏｎ法（棚橋隆彦著者「ＣＦＤの基礎理論」アイピーシー、pp.201-220、1990を参照）を適用し、Ｅｕｌｅｒ陽解法（小林敏雄[ほか]編「数値流体力学ハンドブック」丸善、2003年3月、pp.19−65を参照）を用いて、上記の式（１２）、式（１３）、および式（１４）を時間的に離散化すると、次の式（１５）、式（１６）、および式（１７）を得ることができる。

ここで、上付添字の「ｎ＋１」は次の時刻（ｎ＋１）Δｔにおける値であることを示し、Ｇ_ｘ、Ｇ_ｙ、Ｇ_ｚは、現在時刻ｎから得られた既知の値をまとめた変数を示す。

（１．速度予測子）
上記の式（１５）、式（１６）、および式（１７）は、ｐ＾^ｎ＋１が未知であるため、圧力項を除いて速度予測子

を計算する。そのため、速度予測子はそれぞれ、次の式（１８）、式（１９）および式（２０）のように記載できる。

まず、ｕ、およびｖについて導出を進める。上記の式（１８）、および式（１９）を、上記の式（１５）、および式（１６）に代入することで、次の式（２１）、および式（２２）が得られる。

ここから、次の式（２３）、および式（２４）がそれぞれ得られる。

一方、ｗについては、上記式（２０）を、上記式（１７）に代入することで、次の式（２５）が得られる。

ここから、次の式（２６）が得られる。

次に、上記式（１１）を、時間的に離散化することで、次の式（２７）が得られる。

（２．圧力の長波長成分の式）
上記式（２７）について、Δｘ≒Δｙのスケールで捕捉できる波長成分を考慮して、ｗ’＝ｐ’＝０とし、直交分解すると、ｐ＾^ｎ＋１に関する以下のポアソン方程式（２８）（圧力の長波長成分の式）を得ることできる。

（３．速度の長波長成分の式）
また同様に、上記の式（２１）、式（２２）、および式（２５）について、Δｘ≒Δｙのスケールで捕捉できる波長成分を考慮して、ｗ’＝ｐ’＝０とすると、以下の式（２９）〜式（３１）で示す通り、ｕ、ｖ、ｗ＾の時間更新の式（速度の長波長成分の式）を得ることできる。

（４．圧力の短波長成分の式）
一方、上記式（２７）を、Δｚのスケールで捕捉できる波長成分に関して、直交分解して解くと、ｐ’^ｎ＋１に関する以下のポアソン方程式（３２）（圧力の短波長成分の式）を得ることできる。

（５．速度の短波長成分の式）
また同様に、上記式（２６）を、Δｚのスケールで捕捉できる波長成分に関して解くと、以下の式（３３）で示す通り、ｗ’の時間更新の式（速度の短波長成分の式）を得ることできる。

ここで、ｗ’の計算は、安定性の問題からΔｔの制約が必要である。そのため、Δｔの代わりに、例えば、次の式（３４）、および式（３５）で算出したΔｔ_ｚを用いるとよい。
Δｔ_ｚ＝Δｔ／Ｅ_ｚ（３４）

ここで、ＩＮＴは小数点以下を切り捨てて整数化することを意味している。

〔ポアソン方程式の数値解法の例〕
（１．圧力の長波長成分の式）
圧力の長波長成分の式である、上記のポアソン方程式（２８）を数値的に解く方法を述べる。

通常の数値流体力学の手法に倣い、変数はスタガード配置を用いる。図２に示すように、スカラーを格子中心、ベクトルを格子の境界面に配置する。図２は、スタガード配置の格子例を示す図である。なお、図２では、見易さを考慮して、二次元で描いている。これによって、変数の非物理的な振動を抑えられる。ｘ、ｙ、ｚ座標軸の格子位置をそれぞれｉ、ｊ、ｋとおくと、ポアソン方程式（２８）の第１項、および第２項（＝Ｓijkとおく）は、以下の式（３６）で示すように離散化できる。

また同様に、ポアソン方程式（２８）の第３項（圧力項）は、以下の式（３７）で示すように離散化できる。なお、以下では、記載を容易にするため、時刻をあらわす上付き添え字を省略している。

上記式（２８）に上記の式（３６）および式（３７）を代入すると、ｐに関する演算ができる。その演算法としては大きく分けて直接法と反復法とがあるが、ここでは、反復法の一例であるＧａｕｓｓ−Ｓｅｉｄｅｌ法（小林敏雄[ほか]編「数値流体力学ハンドブック」丸善、2003年3月、ｐ．19−65を参照）により以下の式（３８）を導出する。

ここで、Ｃ_Ｃ、Ｃ_Ｒ、Ｃ_Ｌ、Ｃ_Ｔ、Ｃ_Ｂ、Ｃ_Ｆ、Ｃ_Ｏは、上記式（３７）に含まれる添え字が同一のｐの係数である。上記式（３７）を用いて、ｉ、ｊ、ｋに関するループを組み、収束するまで繰り返し計算を行うことにより、圧力の長波長成分ｐ＾に関する収束値を得る。

（２．圧力の短波長成分の式）
同様に、圧力の短波長成分の式である、上記のポアソン方程式（３２）も、次式（３９）で示すように離散化できる。

上記式（３９）を用いて、圧力の長波長成分ｐ＾と同様に、収束するまで繰り返し計算を行うことにより、圧力の短波長成分ｐ’に関する収束値を得ることができる。

〔圧力境界条件〕
次に、液面の位置を考慮した圧力境界条件について説明する。

本実施形態のように、Δｘ≒Δｙ＞＞Δｚの計算格子を使用する場合、シミュレーションの不安定性や数値誤差を防止できるため、液面の位置を考慮することが好ましい。自由表面を水平方向に対する一価関数として計算することにより、液面の位置を考慮することができる。自由表面を水平方向に対する一価関数として計算する方法としては、例えば、ＳＯＬＡ−Ｓｕｒｆ（SOLution Algorithm-Surf）法等が挙げられる。

そこで、本実施形態では、ＳＯＬＡ−Ｓｕｒｆ法（Hirt, C.W., Nichols, B.D., Romero, N.C. (1975): SOLA-A Numerical Solution Algorithm for Transient Fluid Flow ,Los Alamos Scientific Laboratory, LA-5852, pp. 1-50を参照）に倣い、自由表面を水平方向に対する一価関数ｈ（海底からの液面の位置）として計算する。ｈは、次式（４０）を用いて得ることができる。なお、自由表面を再現できる手法であればよく、VOF（Volume of Fluid）法、ALE（Arbitrary Lagrangian-Eulerian）法、BFC（Boundary Fitted Coordinate）法（いずれも小林敏雄[ほか]編「数値流体力学ハンドブック」丸善、2003年3月、pp.19−65を参照）などの方法も適用可能である。

液面の影響を考慮するため、圧力境界条件を課す必要がある。図３は、自由表面における圧力境界条件を説明するための図である。

ｘ方向およびｚ方向に関する格子位置をそれぞれｉ、ｋと番号付けする。位置（ｉ，ｋ）における圧力をｐ_ｉｋと表記する。スタガード配置に倣い、ｐ_ｉｋは格子中心位置で定義するものとする。

一例として、ｐ_ｉｋの位置とｐ_ｉｋ＋１の位置との間に液面がある場合を考える。セル幅Δｚに対する、液面位置とｐ_ｉｋの位置との距離の比をηとする。液面位置で大気側の圧力がｐ_ａｔｍとなる条件を付与し、線形内挿関係を考慮すれば、次式（４１）を得ることができる。

上記式（４１）で得られる圧力を、圧力の式（３８）、および式（３９）の境界条件として付与する。

Ｐｒｏｊｅｃｔｉｏｎ法は、圧力を直接、反復計算するので、液面の位置を考慮した条件を直接反映できるのに加え、ポアソン方程式の解法を自由に選択することができる。そのため、本実施形態のように、数値解法としてＰｒｏｊｅｃｔｉｏｎ法を採用することにより、液面の位置を圧力境界条件への組み込むことも容易となる。また他の数値解法としてＳＯＬＡ法、ＳＭＡＣ（Simplified Maker and Cell）法等が挙げられる。中でもＳＯＬＡ法は速度と圧力とを同時反復する方法であるため、液面の位置を圧力境界条件として組み込み易い。

本実施形態に係るシミュレーション装置１は、後述するとおり、予め上記で算出された式を記憶させておくことでシミュレーションを行うことができる。

本実施形態に係るシミュレーション装置１の装置構成および処理フローを、以下、詳細に説明する。

〔シミュレーション装置１の構成〕
まず、図１を参照して、シミュレーション装置１の要部構成について説明する。図１は、シミュレーション装置１の要部構成を示すブロック図である。図１に示すように、シミュレーション装置１には、入力装置２と表示装置３とが接続されている。

入力装置２は、シミュレーション装置１を操作するための装置である。入力装置２は、シミュレーション装置１を操作するための入力操作を受け付け、受け付けた入力操作の内容を示す操作データをシミュレーション装置１に送信できるものであればよく、例えばマウス、キーボード等を適用することもできる。

表示装置３は、シミュレーション装置１が出力する表示用データを用いて画像を表示する装置である。表示装置３は、画像を表示できるものであればよく、液晶表示装置、有機ＥＬ（Electro Luminescence）表示装置等を適用することもできる。表示装置３には、シミュレーション装置１を操作するための操作画面、シミュレーション結果等が表示される。

シミュレーション装置１は、制御部１０、および記憶部２０を備えている。

制御部１０は、シミュレーション装置１の動作を統括して制御するものである。制御部１０は、記憶部２０に格納されているプログラムを、例えばＲＡＭ（Random Access Memory）等に読み出して実行することによって上記の制御を行う。

記憶部２０は、シミュレーション装置１が動作するために必要なデータやプログラム等を格納している。また、シミュレーションの過程で算出される途中結果、シミュレーションの最終結果なども記憶部２０に格納される。ここで、記憶部２０には、シミュレーションの過程で用いられる、流体の物性を示すパラメータ（密度、粘度、熱容量、熱伝導度、等）、初期値、境界条件（壁における圧力、潮の流れ、等）、変数等や、計算命令順序等があらかじめ登録されており、シミュレーションの開始時または参照時に読み込ませるものとする。

制御部１０には、シミュレーション算出部３０、シミュレーション制御部５０、および表示制御部６０が含まれ、シミュレーション算出部３０には、速度予測子算出部３１、圧力長波長成分算出部（長波長成分算出部）３２、圧力短波長成分算出部（短波長成分算出部）３３、圧力合成部（合成部）３４、速度長波長成分算出部（長波長成分算出部）３５、速度短波長成分算出部（短波長成分算出部）３６、速度合成部（合成部）３７、および圧力境界条件算出部３８が含まれている。

速度予測子算出部３１は、ｘ，ｙ，ｚ座標軸の各計算格子位置（ｉ、ｊ、ｋ）について、ｘ、ｙ、ｚ方向の速度予測子を算出するものである。具体的には、計算格子毎に、上述した式（１８）、式（１９）、および式（２０）を用いて算出する。

圧力長波長成分算出部３２は、ｘ，ｙ，ｚ座標軸の各計算格子位置（ｉ、ｊ、ｋ）について、圧力の長波長成分を算出するものである。具体的には、計算格子毎に、上述した式（３８）を用いて、収束するまで繰り返し演算を行い、収束値を圧力の長波長成分の算出結果として圧力合成部３４に通知する。

圧力短波長成分算出部３３は、ｘ，ｙ，ｚ座標軸の各計算格子位置（ｉ、ｊ、ｋ）について、圧力の短波長成分を算出するものである。具体的には、計算格子毎に、上述した式（３９）を用いて、収束するまで繰り返し演算を行い、収束値を圧力の短波長成分の算出結果として圧力合成部３４に通知する。

圧力合成部３４は、ｘ，ｙ，ｚ座標軸の各計算格子位置（ｉ、ｊ、ｋ）について、圧力長波長成分算出部３２が算出した圧力の長波長成分と、圧力短波長成分算出部３３が算出した圧力の短波長成分とを足し合わせ、合計値を当該計算格子の圧力としてシミュレーション制御部５０に通知する。

速度長波長成分算出部３５は、ｘ，ｙ，ｚ座標軸の各計算格子位置（ｉ、ｊ、ｋ）について、速度の長波長成分を算出するものである。具体的には、計算格子毎に、上述した式（２９）、式（３０）、および式（３１）を用いて、ｘ、ｙ、ｚ方向の速度の長波長成分を算出する。そして、算出結果のうち、ｘ、ｙ方向の速度の長波長成分については、当該計算格子のｘ，ｙ方向の速度としてシミュレーション制御部５０に通知する。また、ｚ方向の速度の長波長成分については、速度合成部３７に通知する。

速度短波長成分算出部３６は、ｘ，ｙ，ｚ座標軸の各計算格子位置（ｉ、ｊ、ｋ）について、ｚ方向の速度の短波長成分を算出するものである。具体的には、計算格子毎に、上述した式（３３）を用いてｚ方向の速度の短波長成分を算出し、算出結果を速度合成部３７に通知する。

速度合成部３７は、ｘ，ｙ，ｚ座標軸の各計算格子位置（ｉ、ｊ、ｋ）について、速度長波長成分算出部３５が算出したｚ方向の速度の長波長成分と、速度短波長成分算出部３６が算出したｚ方向の速度の短波長成分とを足し合わせ、合計値を当該計算格子のｚ方向の速度としてシミュレーション制御部５０に通知する。

圧力境界条件算出部３８では、各現在時刻ｔについて、上記式（４０）を用いて、次の時刻ｔ＋Δｔのシミュレーションの計算に用いる境界条件として、現在時刻の液面位置ｈを算出する。さらに、上記式（４１）を用いて、圧力境界条件を算出する。

シミュレーション制御部５０は、シミュレーション算出部３０を用いて、各時刻、各計算格子位置についての計算を行うことにより、流体のシミュレーションを実行し、実行結果を表示制御部６０に送信する。なお、シミュレーション制御部５０は、シミュレーションの実行結果を記憶部２０に格納してもよい。

より詳細には、シミュレーション制御部５０は、シミュレーション算出部３０による、ある計算格子位置における圧力および速度の算出が終了すると、対象の計算格子位置を更新して、更新後の計算格子位置における圧力および速度の算出をシミュレーション算出部３０に実行させる。これを、全ての計算格子位置（ｉ、ｊ、ｋの全ての組み合わせ）について計算が終了するまで、繰り返しシミュレーション算出部３０に実行させる。次に、シミュレーション制御部５０は、時刻を更新して、上述した、全ての計算格子位置における圧力および速度の計算をシミュレーション算出部３０に実行させる。そして、シミュレーション制御部５０は、以上の処理を、シミュレーション終了時刻まで時刻を更新しながら、繰り返す。

表示制御部６０は、シミュレーション制御部５０が実行したシミュレーション結果を表示装置３に表示させる。また、表示制御部６０は、表示させる表示データを記憶部２０から読み出し、表示装置３にこれを表示させる。すなわち、表示制御部６０は、上述したシミュレーション装置１によるシミュレーション結果を表示装置３に表示させたり、シミュレーション装置１を操作するための操作画面データ等を記憶部２０から読み出し、表示装置３に表示させる。

〔シミュレーション装置１における処理の流れ〕
次に、図４を参照して、シミュレーション装置１におけるシミュレーションの流れについて説明する。図４は、シミュレーション装置１におけるシミュレーションの流れを示すフローチャートである。

図４に示すように、まず、速度予測子算出部３１が、上記の式（１８）、式（１９）、および式（２０）を用いて速度予測子を算出する（Ｓ１０１）。

次に、圧力長波長成分算出部３２が、Gauss-Seidel法に倣い、上記式（３８）を用いて圧力の長波長成分ｐ＾^ｎ＋１を算出する（Ｓ１０２、長波長成分算出ステップ）。また、これと並行して、圧力短波長成分算出部３３が、上記式（３９）を用いて圧力の短波長成分ｐ’^ｎ＋１を算出する（Ｓ１０３、短波長成分算出ステップ）。ここで、上記式（４１）で得られる圧力は、上記式（３８）、および式（３９）の境界条件として与えられる。次に、速度長波長成分算出部３５が、ステップＳ１０２で算出した圧力の長波長成分ｐ＾^ｎ＋１を上記の式（２９）、式（３０）、および式（３１）に代入して速度の長波長成分ｕ^ｎ＋１、ｖ^ｎ＋１、ｗ＾^ｎ＋１を算出する（Ｓ１０５、長波長成分算出ステップ）。また、これと並行して速度短波長成分算出部３６が、ステップＳ１０３で算出した圧力の短波長成分ｐ’^ｎ＋１を、上記式（３３）に代入して速度の短波長成分ｗ’^ｎ＋１を算出する（Ｓ１０６、短波長成分算出ステップ）。

次に、圧力合成部３４は、圧力長波長成分算出部３２が算出した圧力の長波長成分ｐ＾^ｎ＋１と圧力短波長成分算出部３３が算出した圧力の短波長成分ｐ’^ｎ＋１とから、ｐ＾^ｎ＋１＋ｐ’^ｎ＋１を算出する（Ｓ１０４、合成ステップ）。

また、速度合成部３７は、速度長波長成分算出部３５が算出した速度の長波長成分ｕ^ｎ＋１、ｖ^ｎ＋１、ｗ＾^ｎ＋１うち、ｚ方向の長波長成分ｗ＾^ｎ＋１と、速度短波長成分算出部３６が算出したｚ方向の速度の短波長成分ｗ’^ｎ＋１とから、ｗ＾^ｎ＋１＋ｗ’^ｎ＋１を算出する（Ｓ１０７、合成ステップ）。

そして、シミュレーション制御部５０は、全ての計算格子位置について、ステップＳ１０１〜Ｓ１０７の算出が終了したか否かを判定し（Ｓ１０８）、終了していなければ（Ｓ１０８でＮＯ）、計算格子位置を更新して（Ｓ１１１）、ステップＳ１０１に戻る。

一方、ステップＳ１０８において、全ての計算格子位置について、ステップＳ１０１〜Ｓ１０７の算出が終了したと判定した場合（Ｓ１０８でＹＥＳ）、シミュレーション制御部５０は、位置を初期化して（Ｓ１０９）、シミュレーションの終了時刻まで、計算が完了しているか否かを判定する（Ｓ１１０）。終了時刻まで計算が完了していない場合（Ｓ１１０でＮＯ）、シミュレーション制御部５０は、圧力境界条件算出部３８に現在時刻における液面位置ｈを算出させて（Ｓ１１２）、時刻を更新する（Ｓ１１３）。そして、ステップＳ１０１に戻る。

一方、ステップＳ１１０において、シミュレーションの終了時刻まで、計算が完了していると判定した場合（Ｓ１１０でＹＥＳ）、シミュレーションンの計算を終了する。

なお、現在時刻とは、ステップＳ１０１〜Ｓ１０８、Ｓ１１１を実行する時刻である。そして、ステップＳ１０８で全ての位置について計算済みとなり（Ｓ１０８でＹＥＳ）、ステップＳ１０９で位置を初期化したとき、次の時刻（シミュレーションを実行する次の時刻）に進めることになる。

ここで、図４では、ステップＳ１０２とＳ１０３との演算を並列に実行するものとして説明した。すなわち、ステップＳ１０２とＳ１０３とは、ステップＳ１０１で速度予測子が算出された後は、互いに依存性のない独立して演算が可能なステップである。もちろん、ステップＳ１０２とＳ１０３とを、１ステップずつ演算することも可能であり、その場合、どのステップから実行しても構わない。

上記のように、シミュレーション装置１におけるシミュレーション方法は、長波長の方程式と短波長の方程式とを別途準備して計算するアプローチである。また、上記式（３２）（すなわち、上記式（３９））は、短いスケールの成分に限定して構成しているため、後述する実施例（計算事例）に示すように、収束に必要な計算回数は極端に増加しない。

ここで、シミュレーション装置１で用いる計算格子のサイズは、鉛直方向のスケール（Δｚ）が、水平方向のスケール（Δｘ＝Δｙ）と比較して、１／２より小さいものとする。ただし、鉛直方向のスケール（Δｚ）が、水平方向のスケール（Δｘ＝Δｙ）と同程度であっても、シミュレーションを実行することは可能である。なお、格子サイズ比と効果との関係については、後述する。

なお、シミュレーション装置１は、Δｘ≒Δｙ＞＞Δｚであることを前提とするので、ｕおよびｖの対流項の計算に、三次精度ＣＩＰ（Constrained Interpolation Profile Scheme）法を適用して、空間分解能の向上を図ることができる。また、ｗの対流項の計算に、一次精度風上差分法を適用して、安定性をさらに向上させることができる。さらに、計算ステップの中で計算負荷が最も高いポアソン方程式はＲｅｄ−ＢｌａｃｋＳＯＲ法を用いた並列計算によって高速化することができる。

〔実施例〕
次に、図５〜１０を参照して、前述した方法で実施したシミュレーションの結果を説明する。

〔計算格子のアスペクト比の影響〕
まず、図５〜７を参照して、本実施形態に係る方法で、計算格子のアスペクト比を変更しても、発散させることなく計算結果を得られ、安定性を保持できることを実証する。図５は、計算対象とした矩形三次元領域を示す図である。また、図６は、計算条件を示す図である。また、図７は、計算結果を示す図である。

ここでは、前述したように、図５に示す矩形三次元領域を計算対象とした。図５に示す矩形三次元領域は、Ｌ_ｘ＝Ｌ_ｙ＝２０ｍ、Ｌ_ｚ＝５ｍであり、計算領域のアスペクト比（Ｌ_ｘ／Ｌ_ｚ）は４となっている。水平方向の計算格子サイズはΔｘ＝Δｙ＝１ｍとしている。流れが三次元運動を起こすように、入口境界および出口境界（いずれも２ｍ×２．５ｍの面）を対向する面の対角に位置するように配置し、それ以外はすべり壁としている。入口境界および出口境界は、ｙｚ面と同じ法線を持ち、速度ｕを、ｕ＝１．０ｍ／ｓで与えた。

この設定によって，入口境界では、ｘ方向のみに運動量を持つことになるが、圧力を介してｙ方向およびｚ方向にも運動量が伝えられる。

初期条件として、全方向速度を０ｍ／ｓとした。また、時間刻み幅Δｔは０．１秒とし、滞留時間に相当する４００秒後までの過渡計算を実施した。また、鉛直方向の計算格子サイズΔｚは、図６に示すように、１ｍ、０．１ｍ、０．０１ｍの３種類に変えて、速度場に及ぼす影響を比較した。

計算結果を図７に示す。図７の（ａ）は、実行ナンバー１の結果を示し、（ｂ）は、実行ナンバー２の結果を示し、（ｃ）は、実行ナンバー３の結果を示す。また、図７では、計算開始から４００秒後の、ｚ＝２．５ｍの位置におけるｘｙ断面上の速度場を示している。図７を参照すると、全ての実行結果において、壁面付近で圧力を介して速度の向きが変化している傾向が一致していることが分かる。また、パーソナルコンピュータ（CPU:Intel Xeon 3.4Ghz,Memory:32GB）を用いて計算に要した時間は、最も負荷の高い実行ナンバー３で約７０時間であった。なお、従来（小林敏雄[ほか]編「数値流体力学ハンドブック」丸善、2003年3月、pp.19−65に記載の方法）の計算手法を用いた場合、実行ナンバー１以外は、時間ステップ１００以下で発散し、計算結果を得ることができなかった。

〔海洋に用いた場合の計算例〕
次に、図８〜１０を参照して、本実施形態を海洋に用いた場合の計算例を説明する。図８は、海洋に用いた場合の計算例を説明するための図であり、瀬戸内海の深さ情報を示す図である。図９は、海洋に用いた場合の計算例を説明するための図であり、瀬戸内海の潮流の時間と速度との関係を示す図である。図１０は、海洋に用いた場合の計算例の結果を示す図である。

本計算例では、海洋の事例として瀬戸内海の潮流計算をおこなった。図８に示す、香川県、愛媛県、広島県に接する瀬戸内海を計算対象とした。海底形状は、海上保安庁日本海洋データセンターの水深データから入力した。

計算格子はΔｘ＝Δｙ＝５００ｍ、Δｚ＝５ｍのサイズであり、総格子数は、２５５０００である。瀬戸内海の東端を入口境界、西端を出口境界として設定した。境界条件である入口流速として、図９に示す観測期間（1992年6月2日4:00〜1992年6月17日3:00）のデータを基に求めた時間平均値４．９７ｍ／ｓを与えた。

計算結果を図１０に示す。図８に示したように、計算対象となる領域は、多島を含み水深が５０ｍ以下の領域が多いにもかかわらず、図１０に示すように、本実施形態では、この領域の複雑な潮流を適切に計算できていることが分かる。

〔まとめ〕
以上のように、本実施形態では、鉛直方向に対し準動的な運動方程式を導出し、Ｐｒｏｊｅｃｔｉｏｎ法に適用するために離散化している。さらに、自由表面を境界条件に組み込んで計算方法を拡張している。これにより、Δｘ≒Δｙ＞＞Δｚの計算格子を用いた計算を行った結果、以下の結論を得ている。
（１）格子サイズのアスペクト比が１から離れた場合の計算であっても、本実施形態では、計算可能となった。
（２）例えば瀬戸内海のような、多島を含み、水深が数ｍの領域を多く含んだ海域を対象としても、本実施形態では、複雑な潮流を計算できた。
（３）格子サイズ比と本実施形態の効果との関係を具体的に示せば以下の通りである。格子サイズ比（Δｘ／Δｚ）が「１００以上」の場合、従来の方法では、演算が収束しない可能性があるが、本実施形態では上述したように、適切に演算を実行し、シミュレーションを行うことができる。なお、多くの海洋計算や静水圧近似の計算では格子サイズ比（Δｘ／Δｚ）が１００以上となる可能性が高いため、本実施形態による方法が有利となる。

なお、本実施形態は本発明の範囲を限定するものではなく、本発明の範囲内で種々の変更が可能であり、例えば、以下のように構成することができる。

本発明に係る計算方法は、三次元空間を運動する非圧縮性流体あるいは擬圧縮性のニュートン流体の計算方法であって、必要となる分解能になるように空間領域を分割した場合に、その三次元領域のうち少なくとも一つの座標スケールが、他の座標スケールの２倍以上、望ましくは１０の２乗倍より大きい要素を対象とし、その運動方程式が偏微分方程式の楕円性を持つ場合において、該偏微分方程式を各座標スケールに応じた方程式に分解し、さらに直交分解することで圧力について得られた方程式を複数個用いることで、発散させることなく分解したスケールごとに計算できる方法であってもよい。なお，擬圧縮性のニュートン流体とは、流体の密度が圧力以外の物理因子によって変化できる流体を指し、例えば温度や化学種濃度には依存して密度を変えられる。

さらに、本発明に係る計算方法は、液面などあらわすための圧力に関する境界条件を有した状態で安定な水象を計算する方法であってもよい。

本実施形態では、Δｘ≒Δｙ＞＞Δｚの場合について説明した。しかし、本発明は、少なくとも１つの座標スケールが他の座標スケールの少なくとも１つの２倍以上であればよく、例えば、Δｘ＞＞Δｙ≒Δｚ、Δｘ＜＜Δｙ≒Δｚ、Δｘ＞＞Δｙ＞＞Δｚのケース等にも適用することが可能である。すなわち、座標方向に分解する格子のスケールが異なる場合、各スケールに応じて方程式を分解すれば、水系環境における数値シミュレーションの演算を発散させることなく高精度に実行でき、安定性を向上させた計算が可能になる。

また本発明は、極座標等の直交系でも適用可能である。さらに、例えばBFC（Boundary Fitted Coordinate）法等の座標変換を施す方法に適用することで、非ユークリッド座標系（小林敏雄[ほか]編「数値流体力学ハンドブック」丸善、2003年3月、pp.19−65等を参照）にも適用できる。

〔ソフトウェアによる実現例〕
シミュレーション装置１の制御ブロック（特に、シミュレーション算出部３０、シミュレーション制御部５０、および表示制御部６０）は、集積回路（ＩＣチップ）等に形成された論理回路（ハードウェア）によって実現してもよいし、ＣＰＵ（Central Processing Unit）を用いてソフトウェアによって実現してもよい。

後者の場合、シミュレーション装置１は、各機能を実現するソフトウェアであるプログラムの命令を実行するＣＰＵ、上記プログラムおよび各種データがコンピュータ（またはＣＰＵ）で読み取り可能に記録されたＲＯＭ（Read Only Memory）または記憶装置（これらを「記録媒体」と称する）、上記プログラムを展開するＲＡＭ（Random Access Memory）などを備えている。そして、コンピュータ（またはＣＰＵ）が上記プログラムを上記記録媒体から読み取って実行することにより、本発明の目的が達成される。上記記録媒体としては、「一時的でない有形の媒体」、例えば、テープ、ディスク、カード、半導体メモリ、プログラマブルな論理回路などを用いることができる。また、上記プログラムは、該プログラムを伝送可能な任意の伝送媒体（通信ネットワークや放送波等）を介して上記コンピュータに供給されてもよい。なお、本発明は、上記プログラムが電子的な伝送によって具現化された、搬送波に埋め込まれたデータ信号の形態でも実現され得る。

本発明は上述した各実施形態に限定されるものではなく、請求項に示した範囲で種々の変更が可能であり、異なる実施形態にそれぞれ開示された技術的手段を適宜組み合わせて得られる実施形態についても本発明の技術的範囲に含まれる。さらに、各実施形態にそれぞれ開示された技術的手段を組み合わせることにより、新しい技術的特徴を形成することができる。

１シミュレーション装置
２入力装置
３表示装置
１０制御部
２０記憶部
３０シミュレーション算出部
３１速度予測子算出部
３２圧力長波長成分算出部（長波長成分算出部）
３３圧力短波長成分算出部（短波長成分算出部）
３４圧力合成部（合成部）
３５速度長波長成分算出部（長波長成分算出部）
３６速度短波長成分算出部（短波長成分算出部）
３７速度合成部（合成部）
３８圧力境界条件算出部
５０シミュレーション制御部
６０表示制御部

Claims

三次元空間を運動する非圧縮性流体あるいは擬圧縮性のニュートン流体のシミュレーション装置であって、
三次元空間に設定した計算格子の少なくとも１つの座標スケールが他の座標スケールの少なくとも１つの２倍以上であり、基礎式から導出した偏微分方程式を各座標スケールに対応して分解し、さらに直交分解した方程式を用いるものであり、
上記少なくとも１つの座標スケールに対応した方程式を用いて、圧力および三方向成分の速度を算出する長波長成分算出部と、
上記他の座標スケールに対応した方程式を用いて、圧力および上記他の座標スケールに相当する方向成分の速度を算出する短波長成分算出部と、
上記長波長成分算出部、および上記短波長成分算出部の算出結果を足し合わせる合成部と、を備えることを特徴とするシミュレーション装置。
現在の時刻の液面の位置を、次の時刻の圧力境界条件とすることを特徴とする請求項１に記載のシミュレーション装置。
上記少なくとも１つの座標スケールが、上記他の座標スケールの１０の２乗倍より大きいことを特徴とする請求項１または２に記載のシミュレーション装置。
三次元空間を運動する非圧縮性流体あるいは擬圧縮性のニュートン流体のシミュレーション装置の制御方法であって、
三次元空間に設定した計算格子の少なくとも１つの座標スケールが他の座標スケールの少なくとも１つの２倍以上であり、基礎式から導出した偏微分方程式を各座標スケールに対応して分解し、さらに直交分解した方程式を用いるものであり、
上記少なくとも１つの座標スケールに対応した方程式を用いて、圧力および三方向成分の速度を算出する長波長成分算出ステップと、
上記他の座標スケールに対応した方程式を用いて、圧力および上記他の座標スケールに相当する方向成分の速度を算出する短波長成分算出ステップと、
上記長波長成分算出ステップ、および上記短波長成分算出ステップにおける算出結果を足し合わせる合成ステップと、を含むことを特徴とするシミュレーション装置の制御方法。
請求項１に記載のシミュレーション装置としてコンピュータを機能させるための制御プログラムであって、上記長波長成分算出部、上記短波長成分算出部、および上記合成部としてコンピュータを機能させるための制御プログラム。