JP2011248878A

JP2011248878A - 流体構造連成シミュレーション方法、装置及びプログラム

Info

Publication number: JP2011248878A
Application number: JP2011107604A
Authority: JP
Inventors: Takashi Yamazaki; 崇史山▲崎▼; Yoshimasa Kadooka; 良昌門岡; Toshiaki Hisada; 俊明久田; Seiryo Sugiura; 清了杉浦; Takumi Washio; 巧鷲尾; Junichi Okada; 純一岡田
Original assignee: Fujitsu Ltd; University of Tokyo NUC
Current assignee: Fujitsu Ltd; University of Tokyo NUC
Priority date: 2010-05-24
Filing date: 2011-05-12
Publication date: 2011-12-08
Anticipated expiration: 2031-05-12
Also published as: US20110288834A1; EP2390802A2; JP5737621B2; US10114911B2; EP2390802A3

Abstract

【課題】流体構造連成シミュレーション方法、装置及びプログラムにおいて、流体と構造間のメッシュの整合性に依存せず、正確なシミュレーションを行うことを目的とする。
【解決手段】流体構造連成シミュレーションにおいて、メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成し、物理現象をシミュレートする主要時間発展ループをコンピュータにより実行し、主要時間発展ループは、構造の境界上に構造の変位に伴って移動する相互作用仲介要素を配置し、相互作用仲介要素内に、流体の圧力と速度及び構造の変位と相互作用する、流体の圧力及び速度の補正函数を定義し、補正函数に基づいて流体と構造でメッシュが不整合の状態でシミュレーションを実行するように構成する。
【選択図】図１２

Description

本発明は、流体構造連成シミュレーション方法、装置及びプログラムに関する。本発明は、流体構造連成シミュレーションプログラムを格納したコンピュータ読み取り可能な記憶媒体にも関する。

流体構造連成シミュレーションとは、流体構造連成（FSI: Fluid Structure Interaction）解析を用いた数値解析シミュレーションのことを言う。

流体と構造間の自由境界（又は、移動境界）は普遍的に存在し、近年の計算機の性能向上と有限要素法等の数値シミュレーション手法の高度化に伴って、基礎科学においてのみならずＣＡＥ（Computer Aided Engineering）、生体力学等の応用的分野においても流体構造連成シミュレーションが行なわれるようになった。これにより、より現実的な問題設定下における現象の理解が可能となってきている。流体構造連成問題に対するシミュレーションには、境界追跡（IT:Interface-Tracking）型のＡＬＥ（Arbitrary Lagrangian-Eulerian）法が用いられることが多い（例えば、非特許文献４、５）。

図１は、ＩＴ型のＡＬＥ法を用いた数値シミュレーションを説明する図である。ＡＬＥ法のアルゴリズム又はプログラムは、図１に示す計算機システム１により実行される。計算機システム１は、入力装置２、計算機３、及び出力装置４を有する。入力装置２は、流体と構造の物理定数（例えば、粘性係数や密度、ヤング率（Young's Modulus）等）と初期条件、計算領域（又は、定義域）やメッシュ情報を計算機３に入力する。出力装置４は、計算機３で計算された各離散時刻ステップでの流体の圧力、速度や構造の変位等の自由度の数値情報、可視化情報等を出力する。

ＡＬＥ法のアルゴリズム又はプログラムは、入力装置２からのメッシュ情報を読み取るメッシュ情報読み取り部６、コネクティビティ（Connectivity）と呼ばれる計算領域を離散化した節点間のグラフ情報（即ち、有限要素を構成する節点リスト）を構成するグラフ情報構成部７、初期条件を設定する初期条件設定部８、及びＮタイムステップ（Ｎ Time Step）を繰り返す主要時間発展ループ部９を含む。

図２は、ＩＴ型のＡＬＥ法の主要時間発展ループ部９の構成の一例を説明する図である。主要時間発展ループ部９は、ＦＳＩ部１１とメッシュ制御部１２を有する。ＦＳＩ部１１は、物理現象をシミュレートする主体であり、時間の離散化手法に対応して構造の要素毎の小マトリクス及び等価節点力を構成するマトリクス作成ルーチン１１１と、時間の離散化手法に対応して流体の要素毎の小マトリクス及び等価節点力を構成するマトリクス作成ルーチン１１２と、構造と流体の境界条件を適用する境界条件適用ルーチン１１３を有する。ＦＳＩ部１１は更に、これらのルーチン１１１〜１１３から得られる要素毎の小行列からの寄与を再構成して全自由度に対する行列方程式を作成するマトリクスアセンブラルーチン（Matrix Assembler Routine）１１４、ＦＳＩ部１１内でニュートン・ラフソン（Newton-Raphson）ループを実行する毎に流体と構造の従属変数（圧力、速度、変位等）の変化量を求解して評価する線形ソルバ（Solver）１１５、及びソルバ１１５が求めた従属変数の変化量とその評価を例えば表示して出力する出力部１１６を有する。

構造の変位をｄで表すと、｜ｄ｜≠０であるため、ＦＳＩ部１１内でNewton-Raphsonループを実行する毎に流体と構造の境界で節点の不整合が生じる。主要時間発展ループ部９内のメッシュ制御部１２は、このメッシュの不整合を回復させる役割を担うものであり、数値的な求解を可能にするための補助的な部分を形成する。具体的には、メッシュ制御部１２は、構造の変位を境界条件とし、流体領域のメッシュを弾性体として扱って変形させることで、メッシュの整合性を回復させる。このため、メッシュ制御部１２は、構造の要素毎の小マトリクスを作成するマトリクス作成ルーチン１２１、構造の変位の境界条件を適用する境界条件適用ルーチン１２２、全自由度に対する行列方程式を作成するマトリクスアセンブラルーチン１２４、及びメッシュ制御部１２内でループを実行する毎にメッシュの整合性を評価するソルバ１２５を有する。メッシュ制御シミュレーションの対象を非線型弾性体として扱う場合には、メッシュ制御部１２内でNewton-Raphsonループが繰り返し実行される。一方、シミュレーションの対象を線型弾性体として扱えるのであれば、メッシュ制御部１２内でNewton-Raphsonループが１回実行されると処理はこのループから抜ける。

このように、ＡＬＥ法ではメッシュ制御部１１によりメッシュ制御を行う必要があるが、メッシュ制御は、特に大変形問題において破綻する可能性がある。メッシュ制御が破綻すると、シミュレーションの実行が止まるので、求めるべき近似解を得ることができない。又、構造同士が接触する問題であれば、メッシュが必ず破綻する（接触により体積がゼロとなる流体メッシュが存在する）のでメッシュ制御の適用が制限される。更に、メッシュ制御部１１が実行するべき計算量は比較的多いので、メッシュ制御部１１（即ち、計算機３）への計算負荷は比較的大きい。従って、従来のＩＴ型のＡＬＥ法を用いた流体構造連成シミュレーション方法は、流体と構造間のメッシュの整合性に依存し、大変形問題や構造同士が接触する問題を扱うことができず、メッシュ制御の計算量が比較的多いので計算機への計算負荷も比較的大きい。

一方、メッシュ制御を行なわずにＩＴ型のＡＬＥ法を用いた流体構造連成シミュレーションを行なう方法も提案されている（例えば、特許文献１、特許文献２、非特許文献３、及び非特許文献６）。しかし、このような提案方法は、流体と構造間のメッシュの整合性に依存するＩＴ型の手法であるため、例えば構造が複雑な幾何形状の場合は流体のメッシュが比較的細かくなり、計算量が比較的多くなるので、計算機への計算負荷も比較的大きい。

ラグランジュ（Lagrange）未定乗数法を用いて、流体と構造の境界付近での流体の非圧縮条件を制約条件として、流体と構造の運動を連成する手法も提案されている（例えば、非特許文献３）。しかし、この手法では、粘性による境界面での粘着条件が無視される。このように、メッシュの整合性に依存しない流体構造連成シミュレーション方法では、正確なシミュレーションを行うことは難しい。

特開平９−２４５０８０号公報特開２０００−３５２５４５号公報

久田俊明，「非線形有限要素法のためのテンソル解析の基礎」，丸善，東京，1992 久田俊明，野口裕之，「非線形有限要素法の基礎と応用」，丸善，東京，1995 久田俊明，鷲尾巧，「心臓弁の流体構造連成シミュレーション法に関する数理的考察」，応用数理，Vol.16,No.2,36-50,2006 登坂宣好，宮田秀明，島崎洋治，野村卓史，志村正幸，畑中勝守，数値流体力学編集委員会編，「数値流体力学シリーズ4移動境界流れ解析」，東京大学出版会，東京，1995 Q. Zhang and T. Hisada, "Analysis of Fluid-Structure Interaction Problems with Structural Buckling and Large Domain Changes by ALE Finite Element Method", Comput. Methods Appl. Mech. Engrg., Vol.190,pp.6341-6357,2001 J. de Hart, G. W. M. Peters, P. J. G. Schreurs, F. P. T. Baaijens, "A three-dimensional computational analysis of fluid-structure interaction in the aortic valve", J. Biomech., Vol.36, No.1,103-112,2003

上記の如く、従来のＩＴ型のＡＬＥ法を用いた流体構造連成シミュレーション方法は、流体と構造間のメッシュの整合性に依存し、メッシュ制御の計算量が比較的多いので計算機への計算負荷も比較的大きい。一方、ＩＣ（Interface-Capturing）型のメッシュの整合性に依存しない従来の流体構造連成シミュレーション方法では、正確なシミュレーションを行うことは難しい。このように、従来の流体構造連成シミュレーション方法では、流体と構造間のメッシュの整合性に依存せずに正確なシミュレーションを行うことは難しいという問題があった。

そこで、本発明は、流体と構造間のメッシュの整合性に依存せず、正確なシミュレーションを行うことのできる流体構造連成シミュレーション方法、装置及びプログラムを提供することを目的とする。

本発明の一観点によれば、メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成処理と、物理現象をシミュレートする主要時間発展ループ処理をコンピュータにより実行し、前記主要時間発展ループ処理は、前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する流体構造連成シミュレーション方法が提供される。

本発明の一観点によれば、メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成部と、物理現象をシミュレートする主要時間発展ループ部を備え、前記主要時間発展ループ部は、前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する流体構造連成シミュレーション装置が提供される。

本発明の一観点によれば、コンピュータに流体構造連成シミュレーションを実行させるプログラムであって、メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成手順と、物理現象をシミュレートする主要時間発展ループ手順を前記コンピュータにより実行させ、前記主要時間発展ループ手順は、前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行するプログラムが提供される。

開示の流体構造連成シミュレーション方法、装置及びプログラムによれば、流体と構造間のメッシュの整合性に依存せず、正確なシミュレーションを行うことができる。

ＩＴ型のＡＬＥ法を用いた数値シミュレーションを説明する図である。ＩＴ型のＡＬＥ法の主要時間発展ループ部の構成の一例を説明する図である。本発明の一実施例における流体構造連成シミュレーション装置の一例を示すブロック図である。曲面Σに適合するメッシュを用いた場合と、実施例における解析による、定常問題（ポワッソン方程式）の数値解を示す図である。曲面Σに適合するメッシュを用いた場合と、実施例における解析による、定常問題（粘性オペレータ）の数値解を示す図である。ＡＬＥ法によるＦＳＩシミュレーション結果の一例を示す図である。実施例によるＦＳＩシミュレーション結果の一例を示す図である。非圧縮条件に基づく拘束条件のみを適用した場合のシミュレーション結果の一例を示す図である。ＡＬＥ法、本実施例、及び非圧縮条件に基づく拘束条件のみを適用した場合の夫々のシミュレーションにおける変位を示す図である。本実施例における圧力と速度の基底函数を説明する図である。２次元三角形要素を例として、各有限要素毎に持つ自由度の独立性を説明する図である。流体と構造の連成を説明する図である。実施例におけるＩＣ型の数値シミュレーションを説明する図である。ＩＣ型の実施例の主要時間発展ループ部の構成の一例を説明する図である。第１のサブパートの構成の一例を示す図である。第１及び第２のサブパートの動作を説明するフローチャートである。図１６の処理Ａを説明する図である。図１６の処理Ｂを説明する図である。第２のサブパートの構成の一例を示す図である。第３のサブパートの構成の一例を示す図である。第３及び第４のサブパートの動作を説明するフローチャートである。図２１の処理ａを説明する図である。図２１の処理ｂを説明する図である。第４のサブパートの構成の一例を示す図である。運動する構造の変位に追随して移動するＩＭＥのイメージを示す図である。圧力の補正函数の基底函数を説明する図である。半開球と極座標の定義を説明する図である。直円柱と極座標の定義を説明する図である。速度の補正函数の基底函数を説明する図である。流体と構造の相互作用を説明する図である。流体と構造の相互作用を説明する図である。相互作用境界Σ上での粘着条件のラグランジュ未定乗数による結果のイメージを示す図である。

開示の流体構造連成シミュレーション方法、装置及びプログラムでは、メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成し、物理現象をシミュレートする主要時間発展ループを実行する。主要時間発展ループは、構造の境界上に構造の変位に伴って移動する相互作用仲介要素を配置し、相互作用仲介要素内に、流体の速度及び構造の変位と相互作用する、流体の圧力及び速度の補正函数を定義し、補正函数に流体と構造の相互作用を担わせて、流体と構造のメッシュが不整合の状態でシミュレーションを実行する。以下に、開示の流体構造連成シミュレーション方法、装置及びプログラムの各実施例を図面と共に説明する。

図３は、本発明の一実施例における流体構造連成シミュレーション装置の一例を示すブロック図である。本実施例では、流体構造連成シミュレーション装置は、計算機システム（又は、コンピュータシステム）により形成される。図３において、計算機システム２０は、ＣＰＵ（Central Processing Unit）２１、ネットワークのインタフェース（Ｉ／Ｆ）２２、記憶装置２３、キーボードやマウス等を含む入力装置２４、及び表示装置２５がバス２６を介して接続された構造を有する。ＣＰＵ２１は、計算機システム２０全体の制御を司ると共に、流体構造連成シミュレーションプログラムを実行することで、流体構造連成シミュレーション機能を実現する。Ｉ／Ｆ２２は、計算機システム２０から送信されるコマンド（又は、要求）やデータをネットワーク（図示せず）へ出力すると共に、ネットワークから受信したコマンドやデータを計算機システム２０に入力する機能を有する。記憶装置２３は、ＣＰＵ２１が実行する流体構造連成シミュレーションプログラムを含む各種プログラム、メッシュで表される流体と構造のメッシュ情報、ＣＰＵ２１が実行する演算の中間データを含む各種データ等を格納する。入力装置２４は、ＣＰＵ２１に各種コマンドやデータを入力するのに用いられる。表示装置２５は、ユーザに操作を促す画面やＣＰＵ２１が実行したシミュレーション結果や評価結果等を表示するのに用いられる。

記憶装置２３は、コンピュータ読み取り可能な記憶媒体により形成可能である。コンピュータ読み取り可能な記憶媒体は、一例として、磁気記録媒体、光記録媒体、光磁気記録媒体、ディスクを記録媒体として用いるディスク装置、ＲＯＭ（Read Only Memory）等を含む半導体記憶装置等を含む。

先ず、本実施例において用いる流体構造連成解析手法について、図４〜図１０と共に説明する。

上記の方程式(８)から、Lagrange未定乗数に対しては、Δuが曲面Σ周辺の節点上に置かれたポテンシャルのソースに対応しており、このΔuを調節することで曲面Σでの境界条件（５）を満足させるものと捉えることができる。

図６は、上記の問題に対するＡＬＥ法によるシミュレーション結果（数値解ｖ）の一例を示す図であり、図７は、本実施例によるシミュレーション結果（数値解ｖ）の一例を示す図である。図８は、比較のために例えば上記非特許文献３の非圧縮条件に基づく拘束条件（３０）のみを適用した場合のシミュレーション結果（数値解ｖ）の一例を示す図である。図６〜図８では、いずれも上端と接触して紙面に垂直におかれた暖簾)即ち、曲面Σ）に流れが当たり、後流に渦を発生させながら曲面Σを押して変位させ、曲面Σの変位が一定になった時刻で観測されるスナップショットを示し、定義域の奥行き方向中心断面での速度が描かれている。尚、図６〜図８描画条件はいずれも同じである。

図７のシミュレーション結果は、図８のシミュレーション結果と比較して曲面Σの表面付近での接線方向の速度が小さくなり、粘性による粘着条件が満たされていることが確認された。図７のシミュレーション結果では、曲面Σの後流での速度は比較的強めであるが、流れの向きや曲面Σの変位も図６のシミュレーション結果に近付いて図８と比較すると改善されることが確認された。

図９は、ＡＬＥ法を適用した場合、本実施例を適用した場合、及び例えば上記非特許文献３のような非圧縮条件に基づく拘束条件（３０）のみを適用した場合の夫々のシミュレーションにおける曲面Σの下端中央の変位(displacement)ｄのｘ成分の時間変化を示す図である。図９中、縦軸は変位ｄのｘ成分、横軸は時間ｔを表す。図９において、Ｉ_ＡはＡＬＥ法を適用した場合の変位、Ｉ_Ｂは本実施例を適用した場合の変位、Ｉ_Ｃは非圧縮条件に基づく拘束条件（３０）のみを適用した場合の変位を示す。尚、図６〜図９のシミュレーション結果の計算条件はρ＝１．０×１０^３，μ＝３．８×１０^−３，ρ_ｓ＝２００，ｇ＝０とし、境界条件は左端ｘ＝０においてのみ一様にＶ＝^ｔ（０．３，０，０）を与えた。

このように、本実施例では、Newton流体の粘性オペレータに対してGreenの定理を適用した関係式に基づいて、メッシュの整合性に依存しない流体構造連成シミュレーションを行ない、シェル要素を用いた薄膜構造との相互作用を対象として、広く使われているＩＴ型のＡＬＥ法との比較検討を行なった。この比較検討の結果、本実施例を用いることで、構造境界での粘着条件を満たし、従来手法よりもＡＬＥ法によるシミュレーション結果に近い流れと構造の変位が再現できることが確認された。又、一般的な流体要素での基底函数による補間を用いたLagrange未定乗数法であれば、粘着条件によって生じてしまう流体の逆流が抑えられることも確認できた。

次に、上記実施例の流体構造連成シミュレーションを、図１１以降と共に、より詳細に説明する。上記の如く、本実施例では、流体と構造の相互作用についての境界捕捉（ＩＣ：Interface-Capturing）型数値シミュレーションを行う。これは、例えば風船や石鹸膜、暖簾や旗等の薄膜として扱える構造物と周囲の流体との相互作用、剛体として扱える構造や、弾性体として扱える心筋等の厚さを持つ構造と内外の流体との相互作用等を含む。尚、構造として、紐や梁と言った２次元的で表面を定義できない構造の有限要素は本実施例の対象としない。

本実施例の流体構造連成シミュレーション方法は、次のような手順（又は、ステップ）ＳＴ１〜ＳＴ３を含んでも良い。

手順ＳＴ１：流体領域と構造の相互作用境界は、構造の変位のために移動境界（又は、自由境界）となるが、この移動境界の表裏に運動する構造の変位に追随して移動する有限要素を定義する。即ち、この有限要素内でのみ非零の値を持つ流体の圧力及び速度の補正函数を定義する。

手順ＳＴ２：定義された補正函数が、流体と構造の相互作用を仲介しながら圧力と速度勾配の不連続性を表現できるように、値を持つ領域（サポート又は台）や、構造上での函数値、微分値、要素境界の函数値、微分値等について夫々に条件を課す。

手順ＳＴ３：流体と構造の相互作用を記述するために、従来の流体構造連成シミュレーションで用いられる弱形式とは異なる、構造上の境界積分が露に出現した弱形式を用いる。この弱形式は、Greenの定理に基づいており、流体の非圧縮条件と粘性オペレータの満たす関係式に由来する。

流体構造連成シミュレーションにおいて、運動する構造の変位に追随して移動する有限要素を定義し、非圧縮性や速度の連続性を満足しながら、流体と構造の境界上での粘着条件、構造の表裏での圧力と速度勾配の不連続性を表現し得る補正函数を導入することで、構造の変位ｄと流体の圧力ｐ、速度ｖの相互作用の仲介をさせる。以下の説明では、このように定義した有限要素を「相互作用仲介要素（ＩＭＥ：Interaction Mediating Element）」と呼ぶものとする。

図１１は、２次元三角形要素を例として、各有限要素毎に持つ自由度の独立性を説明する図である。図１１では、節点を○印で示し、流体メッシュ２９の要素毎に、相互作用を持つ節点のみが◎印で示され、流体と共有されない場合の構造の節点をハッチング付きの○印で示している。図１１中、（ａ）はＡＬＥ法の場合を示し、（ｂ）は本実施例の場合を示す。図１１（ａ）の場合、流体部を構成する三角形要素の辺と構造を構成する要素の辺が一致し、辺上の節点が両要素で共有され、この三角形要素単位で自由度間の相互作用が閉じている。このため、ＡＬＥ法の場合は、常に流体と構造のメッシュを整合させる必要がある。

一方、図１１（ｂ）の場合、流体部を構成する三角形要素とは独立して構造が存在する。この場合、各有限要素毎に持つ自由度の独立性（即ち、メッシュ整合性に対する非依存性）のため、流体と構造間の相互作用をもつ節点の集合は、構造の変位によって各時間ステップ（又は、離散時刻）毎に変化する。このため、流体と構造の両者が重なりを持つ部分をサーチして、流体と構造間の相互作用を持つ節点の集合を決定すれば良い。言い換えれば、このような相互作用を持つ節点の集合として要素を定義すれば良い。従って、本実施例の場合、流体と構造でメッシュが不整合の状態でシミュレーションを実行できる。

本実施例においてこのようにメッシュが不整合の状態でシミュレーションを実行可能にするのが、構造表面上に配置されている流体と構造の間の相互作用を仲介するＩＭＥである。流体と構造の連成には、Greenの定理に基づき、流体の非圧縮性と粘性オペレータが満たす関係式に依存する拘束条件を用いる。

図１２は、流体と構造の連成を説明する図である。図１２中、（ａ）はＡＬＥ法を用いた流体と構造の連成を示し、（ｂ）は本実施例による流体と構造の連成を示す。

図１２（ａ）に示すように、ＡＬＥ法を用いた連成処理３１０の場合、流体３１１と構造３１２は夫々の領域を重複なく分割する小領域毎に相互作用する節点の集合が閉じており（有限要素）、その節点の持つ自由度が相互作用し合う。節点の持つ自由度は、流体３１１ならば例えば圧力ｐ及び速度ｖであり、構造３１２ならば変位ｄ等である。

一方、図１２（ｂ）に示すように、本実施例を用いた連成処理３０の場合、流体３１の持つ自由度は直接に構造３２の自由度と相互作用しない。移動する構造３２上に相互作用の仲介となるＩＭＥ３３が定義され、このＩＭＥ３３内に流体の圧力と速度についての補正函数が定義される。つまり、構造３２の境界上に、構造３２の変位に伴って移動する有限要素であるＩＭＥ３３を配置し、このＩＭＥ３３内に流体３１の圧力及び速度に対する補正函数を定義する。これらの補正函数は、流体３１の大域的な速度及び構造３２の変位と相互作用する。又、これらの補正函数は、Greenの定理に基づいた、流体３１と構造３２の境界付近での流体３１の非圧縮条件、及び、Newton流体の粘性作用素の満たす関係を制約条件とするLagrange未定乗数である。このように、ＩＭＥ３３は構造３２の変位に追随して移動するため、構造３２の持つ自由度との相互作用ペアは変化しないが、時々刻々のＩＭＥ３３の移動に伴って、流体３１の持つ自由度との相互作用ペアは変化する。そこで、ＩＭＥ３３と重なりを持つバックグラウンドの流体要素をサーチし、その流体要素のメンバである流体節点が相互作用ペアとして決まる。ＩＣ型の数値シミュレーション、即ち、流体構造連成シミュレーションは、相互作用ペアに基づいて実行すれば良い。このようにして、流体３１の各自由度はＩＭＥ３３を通して移動する構造３２との相互作用境界を認識し、適切に相互作用することができるようになり、流体３１と構造３２のメッシュの整合性への依存性を断ち切ることができる。

図１３は、本実施例におけるＩＣ型の数値シミュレーションを説明する図である。ＩＣ型数値シミュレーションのアルゴリズム又はプログラムは、図１３に示す計算機システム４１により実行される。計算機システム４１は、入力装置４２、計算機４３、及び出力装置４４を有する。入力装置４２は、例えば図３の入力装置２４又はＩ／Ｆ２２で構成可能であり、流体と構造の物理定数（例えば、粘性や密度、ヤング率（Young's Modulus）等）と初期条件、計算領域（又は、定義域）やメッシュ情報を計算機４３に入力する。出力装置４４は、例えば図３の表示装置２５で構成可能であり、計算機４３で計算された各離散時刻ステップでの流体の圧力、速度や構造の変位等の自由度の数値情報、可視化情報等をシミュレーション結果として出力する。本実施例では、ＡＬＥ法を用いた場合と異なり、メッシュのグラフ（又は、コネクティビティ）を流体と構造の夫々について作成すると共に、構造が３次元的な分布を持つ場合にはこのグラフ作成時に流体と構造の相互作用境界を構成する２次元的な要素についてのグラフを別途作成する。尚、構造が３次元的な分布を持つ場合とは、シェルや膜要素といった２次元的な「面」の要素ではなく、体積を持つ要素である場合を指す。

ＩＣ型数値シミュレーションのアルゴリズム又はプログラムは、入力装置４２からのメッシュ情報を読み取るメッシュ情報読み取り処理を実行するメッシュ情報読み取り部４６、流体と構造の夫々についてコネクティビティと呼ばれる計算領域を離散化した節点間のグラフ情報（即ち、有限要素を構成する節点リスト）を構成するグラフ情報構成処理を実行するグラフ情報構成部４７、初期条件を設定する初期条件設定処理を実行する初期条件設定部４８、及びＮタイムステップ(N Time Step)が繰り返される主要時間発展ループ処理を実行する主要時間発展ループ部４９を含む。

図１４は、ＩＣ型の実施例の主要時間発展ループ部の構成の一例を説明する図である。主要時間発展ループ部４９は、ＦＳＩ部５１を有する。ＦＳＩ部５１は、物理現象をシミュレートする主体であり、時間の離散化手法に対応して構造の要素毎の小マトリクス及び等価節点力を構成するマトリクス作成ルーチン１５１と、時間の離散化手法に対応して流体の要素毎の小マトリクス及び等価節点力を構成するマトリクス作成ルーチン１５２と、ＩＭＥ毎に速度、圧力と補正函数の相互作用を表すマトリクスを構成するマトリクス作成ルーチン１５３、これらのルーチン１５１〜１５３から得られる要素毎の小行列からの寄与を再構成して全自由度に対する行列方程式を作成するマトリクスアセンブラルーチン（matrix assembler routine）１５５、ＦＳＩ部内でニュートン・ラフソン（Newton-Raphson）ループを実行する毎に流体と構造の従属変数（圧力、速度、変位等）の変化量を求解して評価する線形ソルバ１５６、及びソルバ１５６が求めた従属変数の変化量とその評価を例えば表示して出力する出力部１５７を有する。

マトリクス作成部１５３内の境界条件適用部は、補正函数と構造の変位の相互作用を表すモジュール６１〜６４で形成されている。これらのモジュール６１〜６４によって、運動する構造の変位に追随して移動するＩＭＥが生成され、流体と構造のメッシュの分離（即ち、メッシュ整合性の非依存）が実現される。以下の説明では、モジュール６１〜６４のことを、サブパート６１〜６４と呼ぶ。

第１のサブパート６１は、圧力の補正函数と流体の速度の相互作用マトリクスを計算する。第２のサブパート６２は、速度の補正函数と流体の速度の相互作用マトリクスを計算する。第３のサブパート６３は、圧力の補正函数と構造の変位速度の相互作用マトリクスを計算する。第４のサブパート６４は、速度の補正函数と構造の変位速度の相互作用マトリクスを計算する。以下の説明でも、流体領域（又は、積分領域）をΩで表し、その境界∂Ωのうち移動境界（又は、自由境界）となる流体と構造の相互作用境界（即ち、曲面）をΣで表すものとする。

図１５は、第１のサブパート６１の構成の一例を示す図である。サブパート６１には、ＩＭＥのサイズ（又は、空間スケール）κ、予め作成されている流体領域と構造領域のメッシュの夫々についてのグラフ（又は、コネクティビティ）、及び節点座標が入力される。又、このサブパート６１は、流体と構造の相互作用境界Σを構成する２次元的な要素毎に定義される、構造の変位に追随して移動する有限要素の夫々が、相互作用する流体節点の集合をサーチした結果を返すサーチ部６１１と、サーチ部６１１のサーチ結果に基づいて必要なメモリ領域を確保するアロケーション（allocation）部６１２、及び各節点の自由度への行列の寄与を返す流体寄与計算（fluid contribution computation）部６１３を有する。

本実施例では、流体のメッシュとは独立して構造の変位が決まり、ＩＭＥが相互作用する流体節点のリストは時間発展及びNewton-Raphsonのループの都度異なる。そこで、要素毎のマトリクスを作成するために、サーチを行ない相互作用する節点リストを更新する必要がある。サーチ部６１１は、このサーチを行うのに例えばセル・インデックス法（又は、バケット・ソート）を用いる。セル・インデックス法は定義域Ωを矩形セルに分割して、全ての要素をその体積中心の座標に従ってセルに対応させ、セル単位でサーチを行なうものである。このように、サーチ部６１１は、各ＩＭＥが相互作用する流体節点の数を決定する。具体的には、ＩＭＥ内部に離散的に配置された積分点を用いて積分を行なう場合であれば、夫々の積分点毎に、内包する流体要素をセル・インデックス法を用いてサーチし、その流体要素の構成メンバである速度節点をカウントする。大域的な流体要素毎にGauss積分を用いて積分を行なう場合には、ＩＭＥの体積中心を含むセルを中心とした周囲のセルに含まれる流体要素のGauss積分の積分点をＩＭＥが包含するかどうかを確認し、包含する場合はその流体要素の構成メンバである速度節点をカウントする。従ってサーチ部６１１は、相互作用境界Σを構成する各２次元的要素毎に定義されるＩＭＥ毎に、（ｓ１）相互作用境界と共に移動するＩＭＥ内の積分点座標を求め、（ｓ２）積分点を包含する流体要素をサーチし、（ｓ３）メッシュのグラフ情報に基づいて、サーチされた流体要素のメンバである流体節点（又は、自由度）を相互作用節点数として重複なくカウントする。或いは、（ｓ１’）相互作用境界と共に移動するＩＭＥの体積中心座標を求め、（ｓ２’）その体積中心から周囲の流体要素をサーチし、（ｓ３’）ＩＭＥが流体要素のGauss積分点を包含する場合は、流体要素のメンバである流体節点（又は、自由度）を相互作用節点数として重複なくカウントする。

アロケーション部６１２は、サーチ部６１１がカウントした相互作用節点数に基づいて、マトリクスへの寄与を記憶するためのメモリ領域を確保する。

流体寄与計算部６１３は、サーチ部６１１と同様の処理を行なうことで、ＩＭＥ毎に相互作用する流体節点のリストを作成し、同時に相互作用マトリクスへの寄与を計算する。従って、流体寄与計算部６１３は、ＩＭＥ内部に離散的に配置された積分点を用いて積分を行なう場合であれば、相互作用境界Σを構成する各２次元的要素毎に定義されるＩＭＥ毎に、（ｓ１１）相互作用境界と共に移動するＩＭＥ内の積分点座標を求め、（ｓ１２）積分点を包含する流体要素をサーチし、（ｓ１３）メッシュのグラフ情報に基づいて、サーチされた流体要素のメンバである流体節点（又は、自由度）をメモリに記憶し、（ｓ１４）積分点での基底函数値及び微分値を求め、（ｓ１５）要素内小マトリクスへの寄与を足し合わせる。或いは、大域的な流体要素毎にGauss積分を用いて積分を行なう場合には、（ｓ１１’）相互作用境界と共に移動するＩＭＥの体積中心座標を求め、（ｓ１２’）その体積中心から周囲の流体要素をサーチし、（ｓ１３’）ＩＭＥが流体要素のGauss積分点を包含する場合は、メッシュのグラフ情報に基づいて、流体要素のメンバである流体節点（又は、自由度）をメモリに記憶し、（ｓ１４’）積分点での基底函数値及び微分値を求め、（ｓ１５’）要素内小マトリクスへの寄与を足し合わせる。

このメモリ、及び後述する各メモリは、例えば図３に示す記憶装置２３により形成可能である。これにより、圧力の補正函数と大域的な流体の速度との相互作用を表す要素毎の小マトリクスが作成される。

図１６は、第１及び第２のサブパートの動作を説明するフローチャートである。又、図１７は図１６の処理Ａを説明する図であり、図１８は図１６の処理Ｂを説明する図である。図１６の処理は、ＩＭＥ内部に離散的に配置された積分点を用いて積分を行なう場合を例として、処理の流れをまとめたものである。尚、第１のサブパート６１の処理の場合、上記の手順（ｓ１）〜（ｓ３）では処理Ａ，Ｂは行なわれず、上記の手順（ｓ１１）〜（ｓ１５）で処理Ａ，Ｂが行なわれる。

第１のサブパート６１の処理の場合、図１６のステップＳ１は、ポインタをpointer(1)=1に設定すると共に、図１７に示す処理Ａを実行する。処理Ａは、フラグndflag及びフラグmatimeがアロケート（allocate）されている場合にのみ実行される。フラグndflag，matimeは、ＩＭＥ内部に離散的に配置された積分点を用いて積分を行なう場合にアロケートされる。図１７に示すステップＳ１０１は、フラグndflagをndflag(1:pointer(#IME+1)-1)=0に設定し、第１のサブパート６１の場合はフラグmatimeをmatime(1:3,ndflag(1:pointer(#IME+1)-1)=0に設定し、後述する第２のサブパート６２の場合はフラグmatimeをmatime(1:3,1:3,ndflag(1:pointer(#IME+1)-1)=0に設定する。

図１６の説明に戻るに、ステップＳ２は、ＩＭＥのＩＤ（IDentification）ime_idをime_id=0に初期化する。ステップＳ３は、フラグflagをflag(1:nnode)=0に設定すると共に、リストlistをlist(1:nnode)=0に設定する。ステップＳ４は、ＩＭＥのＩＤであるime_idをインクリメントしてime_id=ime_id+1とし、ステップＳ５は、コネクティビティカウンタncountをncount=0に設定する。ステップＳ６は、三角形要素上のガウス点ＩＤであるigauss_triをigauss_tri=0に設定する。

ステップＳ７は、ガウス点ＩＤをインクリメントしてigauss_tri=igauss_tri+1とし、ステップＳ８は、ガウス点ＩＤ（又は、ガウス点座標）igauss_triと重みweight_triをメモリに記憶する。ステップＳ９は、ＩＭＥ内の離散積分点ＩＤ integをinteg=0に設定する。

ステップＳ１０は、離散積分点ＩＤをインクリメントしてinteg=integ+1とし、ステップＳ１１は、積分点座標tpointとガウス点座標gauss_triの相対関係をrel=tpoint-gauss_triに設定する。ステップＳ１２は、セル・インデックス法によって積分点座標tpointを包含する流体要素jelemをサーチする。ステップＳ１３は、流体要素jelem内の節点ＩＤ jnodeを0に初期化する。

ステップＳ１４は、節点ＩＤをインクリメントしてjnode=jnode+1とする。ステップＳ１５は、節点ＩＤjnodeに対するフラグflagがflag(jnode)==0であるか否かを判定し、判定結果が否（FALSE）であると処理はステップＳ１４へ戻り、判定結果が正（TRUE）であると処理はステップＳ１６へ進む。ステップＳ１６は、コネクティビティカウンタをインクリメントしてncount=ncount+1とすると共に、フラグ及びリストを夫々flag(jnode)=ncount、list(ncount)=jnodeに設定して図１８に示す処理Ｂを実行する。処理Ｂは、フラグndflag及びフラグmatimeがアロケート（allocate）されている場合にのみ実行される。図１８に示すステップＳ１６１は、流体要素jelemにおける、積分点での基底函数と基底函数の勾配を求める。ステップＳ１６２は、ＩＭＥにおける、積分点での基底函数と基底函数の勾配を求める。ステップＳ１６３は、基底函数とその勾配に積分点での重みを乗算して、弱形式における積分項の寄与contribを評価する。ステップＳ１６４は、インデックスをindex=pointer(ime_id)+ncountに設定し、フラグndflagをndflag(index)="global id of inode"に設定し、フラグmatimeをmatime(:,:,index)+=contribに設定する。

図１６の説明に戻るに、ステップＳ１７は、流体要素jelem内の節点ＩＤjnodeが節点数#nodesと等しくjnode==#nodesであるか否かを判定し、判定結果がFALSEであると処理はステップＳ１４へ戻り、判定結果がTRUEであると処理はステップＳ１８へ進む。ステップＳ１８は、ＩＭＥ内の離散積分点ＩＤintegが離散積分点数#integ_pointsと等しくinteg==#integ_pointsであるか否かを判定し、判定結果がFALSEであると処理はステップＳ１０へ戻り、判定結果がTRUEであると処理はステップＳ１９へ進む。ステップＳ１９は、ガウス点ＩＤigauss_triがガウス点数#gauss_triと等しくigauss_tri==#gauss_triであるか否かを判定し、判定結果がFALSEであると処理はステップＳ７へ戻り、判定結果がTRUEであると処理はステップＳ２０へ進む。

ステップＳ２０は、ポインタをpointer(ime_id+1)=pointer(ime_id)+ncountに設定し、ステップＳ２１は、フラグflagをflag(list(1:ncount))=0に設定する。ステップＳ２２は、ＩＭＥのＩＤime_idがＩＭＥ数#IMEと等しくime_id==#IMEであるか否かを判定し、判定結果がFALSEであると処理はステップＳ４へ戻り、判定結果がTRUEであると処理は終了する。

図１９は、第２のサブパート６２の構成の一例を示す図である。サブパート６２への入力は、サブパート６１への入力と同じである。このサブパート６２は、流体と構造の相互作用境界Σを構成する２次元的な要素毎に、構造の表裏に定義される、構造の変位に追随して移動する有限要素の夫々が、相互作用する流体節点の集合をサーチした結果を返すサーチ部６２１と、サーチ部６２１のサーチ結果に基づいて必要なメモリ領域を確保するアロケーション部６２２、及び各節点の自由度への行列の寄与を返す流体寄与計算部６２３を有する。

上記の如く、本実施例では、流体のメッシュとは独立して構造の変位が決まり、ＩＭＥが相互作用する流体節点のリストは、時間発展及びNewton-Raphsonのループの都度異なる。そこで、要素毎のマトリクスを作成するために、サーチを行ない相互作用する節点リストを更新する必要がある。サーチ部６２１は、このサーチを行うのにやはり、例えばセル・インデックス法を用いる。このように、サーチ部６２１は、各ＩＭＥが相互作用する流体節点の数を決定する。具体的には次のようになる。ＩＭＥ内部に離散的に配置された積分点を用いて積分を行なう場合は、夫々の積分点毎に、内包する流体要素をセル・インデックス法を用いてサーチし、その流体要素の構成メンバである速度節点をカウントする。大域的な流体要素毎にGauss積分を用いて積分を行なう場合には、ＩＭＥの体積中心を含むセルを中心とした周囲のセルに含まれる流体要素のGauss積分の積分点をＩＭＥが包含するかどうかを確認し、包含する場合はその流体要素の構成メンバである速度節点をカウントする。従ってサーチ部６２１は、相互作用境界Σを構成する各２次元的要素毎に定義されるＩＭＥ毎に、（ｓ２１）相互作用境界と共に移動するＩＭＥ内の積分点座標を求め、（ｓ２２）積分点を包含する流体要素をサーチし、（ｓ２３）メッシュのグラフ情報に基づいて、サーチされた流体要素のメンバである流体節点（又は、自由度）を相互作用節点数として重複なくカウントする。或いは、（ｓ２１’）相互作用境界と共に移動するＩＭＥの体積中心座標を求め、（ｓ２２’）その体積中心から周囲の流体要素をサーチし、（ｓ２３’）ＩＭＥが流体要素のGauss積分点毎を包含する場合は、流体要素のメンバである流体節点（又は、自由度）を相互作用節点数として重複なくカウントする。

アロケーション部６２２は、サーチ部６２１がカウントした相互作用節点数に基づいて、マトリクスへの寄与を記憶するためのメモリ領域を確保する。

流体寄与計算部６２３は、サーチ部６２１と同様の処理を行なうことで、ＩＭＥ毎に相互作用する流体節点のリストを作成し、同時に相互作用マトリクスへの寄与を計算する。従って、流体寄与計算部６２３は、ＩＭＥ内部に離散的に配置された積分点を用いて積分を行なう場合であれば、相互作用境界Σを構成する各２次元的要素毎に、相互作用境界Σの表裏に定義されるＩＭＥ毎に、（ｓ３１）相互作用境界と共に移動するＩＭＥ内の積分点座標を求め、（ｓ３２）積分点を包含する流体要素をサーチし、（ｓ３３）メッシュのグラフ情報に基づいて、サーチされた流体要素のメンバである流体節点（又は、自由度）をメモリに記憶し、（ｓ３４）積分点での基底函数値及び微分値を求め、（ｓ３５）要素内小マトリクスへの寄与を足し合わせる。或いは、（ｓ３１’）相互作用境界と共に移動するＩＭＥの体積中心座標を求め、（ｓ３２’）その体積中心から周囲の流体要素をサーチし、（ｓ３３’）ＩＭＥが流体要素のGauss積分点毎を包含する場合は、メッシュのグラフ情報に基づいて、流体要素のメンバである流体節点（又は、自由度）をメモリに記憶し、（ｓ３４’）積分点での基底函数値及び微分値を求め、（ｓ３５’）要素内小マトリクスへの寄与を足し合わせる。

第２のサブパート６２の処理も、上記第１のサブパート６１の処理の場合と同様に、図１６に示す処理の手順を実行すれば良い。尚、第２のサブパート６２の処理の場合、上記の手順（ｓ２１）〜（ｓ２３）では処理Ａ，Ｂは行なわれず、上記の手順（ｓ３１）〜（ｓ３５）で処理Ａ，Ｂが行なわれる。又、処理ＡのステップＳ１０１の処理が図１７と共に説明したように、第１のサブパート６１の処理の場合と一部異なる。

これにより、速度の補正函数と大域的な流体の速度との相互作用を表す要素毎の小マトリクスが作成される。

図２０は、第３のサブパート６３の構成の一例を示す図である。サブパート６３への入力は、サブパート６１への入力と同じである。このサブパート６３では、各要素と相互作用する構造要素は、要素が配置されている構造要素のみであるとして、サーチは行なわない。

サブパート６３は、必要なメモリ領域を確保するアロケーション部６３１、及び各節点の自由度への行列の寄与を返す構造寄与計算（structure contribution computation）部６３２を有する。構造寄与計算部６３２は、圧力の補正函数と構造の速度との相互作用を表すためにGreenの定理を用いて取り出した、相互作用境界Σ上の境界積分項の計算を担う。

アロケーション部６３１は、各ＩＭＥ毎に、相互作用を行なう構造の変位節点は、ＩＭＥが配置されている構造要素の節点だけであるとして、マトリクスへの寄与を記憶するために必要なメモリ領域を確保する。

構造寄与計算部６３２は、各ＩＭＥが配置されている構造要素の変位節点を相互作用する節点とし、相互作用マトリクスへの寄与を計算する。具体的には、相互作用境界Σを構成する各２次元的要素毎に定義されるＩＭＥ毎に、（ｓ４１）相互作用ペアである、ＩＭＥが配置されている構造要素のメンバである構造節点（又は、自由度）をメモリに記憶し、（ｓ４２）要素内小マトリクスへの寄与を足し合わせる。これにより、非圧縮条件が適用され、圧力の補正函数と構造の速度との相互作用を表す要素毎の小マトリクスが作成される。

図２１は、第３及び第４のサブパートの動作を説明するフローチャートである。又、図２２は図２１の処理ａを説明する図であり、図２３は図２１の処理ｂを説明する図である。

第３のサブパート６３の処理の場合、図２１のステップＳ３１は、ポインタをpointer(1)=1に設定すると共に図２２に示す処理ａを実行する。図２２に示すステップＳ３１１は、フラグndflagをndflag(1:(3×#IME))=0に設定し、第３のサブパート６３の処理の場合はフラグmatimeをmatime(1:3, ndflag(1:(3×#IME)))=0に設定し、後述する第４のサブパート６４の処理の場合はフラグmatimeをmatime(1:3,1:3, ndflag(1:(3×#IME)))=0に設定する。

図２１の説明に戻るに、ステップＳ３２は、ＩＭＥのＩＤime_idをime_id=0に初期化する。ステップＳ３３は、ＩＭＥのＩＤをインクリメントしてime_id=ime_id+1とし、ステップＳ３４は、裏表を含めた法線ベクトルvnを定義する。又、ステップＳ３５は、三角形の面積areaを定義する。ステップＳ３６は、ポインタをpointer(ime_id+1)=pointer(ime_id)+3に設定し、ステップＳ３７は、三角形要素の頂点ＩＤivertをivert=0に初期化する。

ステップＳ３８は、三角形要素の頂点ＩＤをインクリメントしてivert=ivert+1とする。ステップＳ３９は、インデックスindexをindex=pointer(ime_id)+ivertに設定する。ステップＳ４０は、図２３に示す処理ｂを実行する。図２３に示すステップＳ４０１は、法線ベクトルvnと、要素面積areaに注目して、弱形式における積分項の寄与contribを評価する。ステップＳ４０２は、フラグndflagをndflag(index)="global id of ivert"に設定すると共に、フラグmatimeをmatime(:,:,index)+=contribに設定する。

図２１の説明に戻るに、ステップＳ４１は、三角形要素の頂点ＩＤがivert==3であるか否かを判定し、判定結果がFALSEであると処理はステップＳ３８へ戻り、判定結果がTRUEであると処理はステップＳ４２へ進む。ステップＳ４２は、ＩＭＥのＩＤime_idがＩＭＥ数#IMEと等しくime_id==#IMEであるか否かを判定し、判定結果がFALSEであると処理はステップＳ３３へ戻り、判定結果がTRUEであると処理は終了する。

図２４は、第４のサブパート６４の構成の一例を示す図である。サブパート６４への入力は、サブパート６１への入力と同じである。このサブパート６４では、各要素と相互作用する構造要素は、要素が配置されている構造要素のみであるとして、サーチは行なわない。

サブパート６４は、必要なメモリ領域を確保するアロケーション部６４１、及び各節点の自由度への行列の寄与を返す構造寄与計算（structure contribution computation）部６４２を有する。構造寄与計算部６４２は、速度の補正函数と構造の速度との相互作用を表すためにGreenの定理を用いて取り出した、相互作用境界Σ上の境界積分項の計算を担う。

アロケーション部６４１は、各ＩＭＥ毎に、相互作用を行なう構造の変位節点は、ＩＭＥが配置されている構造要素の節点だけであるとして、マトリクスへの寄与を記憶するために必要なメモリ領域を確保する。

構造寄与計算部６４２は、各ＩＭＥが配置されている構造要素の変位節点を相互作用する節点とし、相互作用マトリクスへの寄与を計算する。具体的には、相互作用境界Σを構成する各２次元的要素毎に、相互作用境界Σの表裏に定義されるＩＭＥ毎に、（ｓ５１）相互作用ペアである、ＩＭＥが配置されている構造要素のメンバである構造節点（又は、自由度）をメモリに記憶し、（ｓ５２）要素内小マトリクスへの寄与を足し合わせる。これにより、粘着条件が適用され、速度の補正函数と構造の速度との相互作用を表す要素毎の小マトリクスが作成される。

第４のサブパート６４の処理も、上記第３のサブパート６３の処理の場合と同様に、図２１に示す処理の手順を実行すれば良い。又、処理ａのステップＳ３１１の処理が図２２と共に説明したように、第３のサブパート６３の処理の場合と一部異なる。

上記の如く、本実施例では、流体領域と構造領域の相互作用境界（即ち、移動境界、自由境界）となる表面（曲面）Σに配置する、構造の変位に伴って移動する有限要素（ＩＭＥ）を用いる。このＩＭＥに基づき、ＩＭＥ内でのみ非零の値を持ち、且つ、相互作用境界以外の要素境界では零となる圧力補正函数、及びＩＭＥ内でのみ非零の値を持ち、且つ、相互作用境界以外の要素境界では勾配が零となる速度の補正函数を用いる。又、流体と構造の相互作用を記述するために通常の流体、構造連成とは異なる弱形式の構成を用いる。

次に、ＩＭＥについてより詳細に説明する。移動境界（又は、自由境界）となる流体と構造の相互作用境界Σは、構造の有限要素分割によって２次元的に要素分割されている。この相互作用境界Σを構成する２次元的な有限要素の中心や頂点、積分点を代表点として圧力の補正函数の節点を配置する。この節点を中心として、相互作用境界Σの表裏にサイズκの球形又は直円柱形の要素を配置し、この要素内部に流体の圧力の補正函数を定義する。

同様にして、相互作用境界Σを構成する２次元的な有限要素の中心や頂点、積分点を代表点として速度の補正函数の節点を相互作用境界Σの表裏に夫々配置する。この節点を中心として、相互作用境界Σの表裏の夫々にサイズκの半球形又は直円柱形の要素を配置し、この要素内部に流体の速度の補正函数を定義する。

図２５は、流体と構造の相互作用境界Σ上に離散的に配置された、運動する構造の変位に追随して移動するＩＭＥのイメージを示す図である。図２５（ａ）は、圧力の補正函数に対するＩＭＥであり、表裏に跨って定義され、相互作用境界Σの表裏で要素が１つの節点に対応する。図２５（ｂ）は、速度の補正函数に対するＩＭＥであり、相互作用境界Σの裏面にも同様に要素が配置される。速度の補正函数に対するＩＭＥでは、相互作用境界Σの表裏の要素に夫々対応する節点が存在する。図２５は、一例として半球形の要素を示し、見やすいように本来よりも要素サイズが小さく描かれているが、ＩＭＥ同士が重なり合うように要素サイズを決めても良い。

運動する構造の変位に追随して移動するＩＭＥのサイズκは、流体と構造の相互作用境界Σ上でＩＭＥ同士が隙間なく重なり合う程度の大きさを採用可能である。例えば、相互作用境界Σの有限要素分割が三角形要素でなされており、その中心に１つずつＩＭＥを配置する場合には、相互作用境界Σを構成する有限要素のサイズの最大値Ｌと、１程度の実数αに対して、κ＝αＬを半径とする半球か、κを半径や高さとする直円柱を流体と構造の相互作用境界Σの表裏に配置すれば良い。

ところで、本実施例における流体構造連成シミュレーションは、流体と構造のメッシュが独立になることを目的の１つとしているが、メッシュの独立性に依らず、流体と構造間の相互作用を表現できなければ、流体と構造の整合性のある運動をシミュレートすることは不可能である。つまり、物理現象に基づいて、流体の物理量（圧力ｐ，速度ｖ）と構造の物理量（変位ｄ）を適切に関連付ける必要がある。そこで、以下の説明では、どのような物理現象が起きているかを確認しながら、定式化との関連を再度説明する。

先ず、どのような現象が起きているか、言い換えれば実現したい現象について簡単に説明する。

いずれにしても、構造の運動に影響を与えるトータルとしての応力は、構造表裏の圧力や速度勾配の不連続性によってもたらされる。これらの物理量の不連続性は、弁によって仕切られた二つのシリンダ内部の流体が弁を挟んで不連続に変化し得る状況を思い浮かべれば理解し易い。速度の連続性は、弁付近に真空が発生することがないことから明らかである。又、弁を挟む二つの部屋で全く異なる流体運動が生じている状況があり得ることを想像すれば、壁の前後で全く異なる応力が発生し、弁の両サイドで全く異なる圧力や速度勾配が生まれ、不連続になることが理解できる。

次に、構造が圧力に押され、圧力が不連続となることを如何に実現するかについて、図３０と共に説明する。図３０は、流体と構造の相互作用を説明する図である。図３０中、梨地付きの○印は高圧節点、白抜きの○印は低圧節点、ハッチング付きの○印は構造の節点を示し、矢印は変位を示す。又、「高圧」は高圧領域を示し、「低圧」は低圧領域を示す。図３０では、構造は流体の圧力に押されて変位し、流体の圧力は構造の前後で不連続となる様子を示す。

図３０に示したように、流体の圧力と構造の相互作用の結果、構造に阻まれて流体の圧力は構造の表裏で空間的に不連続になり、主として構造の表裏の圧力差によって構造は変位する。これは、メッシュに整合性がない状況下であれば、圧力に対して要素単位で連続な基底函数が用いられる通常の方法では実現不可能である（有限要素の基底函数は、通常要素内で圧力や速度勾配が連続函数となるように定義されており、要素内の不連続性を表現する能力はない）。そこで、相互作用境界Σを中心として、構造の運動に追随して移動する有限要素（ＩＭＥ）を定義し、その要素内部をサポートとし、構造の表裏で符合を変える圧力の補正函数を導入して流体の圧力を修飾することで流体圧力の不連続を実現する。又同時に、弱形式における流体の圧力勾配力項と非圧縮条件項の、圧力の補正函数と大域的な流体速度との相互作用項において、Greenの定理を用いて構造上の境界積分を出現させ、流体と構造の速度を一致させることで、変位を制御する。そのため、圧力補正函数は通常の圧力よりも高い正則性が要求される。

又、上記の如く、補正函数のサポートサイズκは、流体と構造の相互作用境界Σ上でＩＭＥ同士が隙間なく重なり合う程度の大きさを採用すれば良い。具体的には、例えばΣの有限要素分割が三角形要素でなされており、その中心に一つずつＩＭＥを配置する場合には、相互作用境界Σを構成する有限要素のサイズの最大値Ｌと、１程度の実数αに対して、κ＝αＬを半径とする半球か、κを半径や高さとする直円柱を流体と構造の相互作用境界Σの表裏に配置する。これらの有限要素内で圧力の補正函数を定義することで、いずれの補正函数もそのサポートがκ程度の長さスケールを持つことになる。

図３２は、例えば非特許文献６のような、相互作用境界Σ上での粘着条件のLagrange未定乗数による結果のイメージを示す図である。相互作用境界Σ上の一点（図３２中、白抜きの○印で示す）での粘着条件をLagrange未定乗数で実装すると、粘着境界条件を実現するために同一三角要素に含まれる速度節点では、構造の一方の端に存在する流体節点での速度（図３２中、太い実線の矢印で示す）からの寄与をキャンセルするように、構造のもう一方の端に存在する流体節点上では逆向きの速度（図３２中、破線の矢印で示す）が発生してしまうことになる。このような現象は、物理現象と整合するとは限らない。

又、速度の補正函数のサポートサイズについても、圧力の場合と同様に決めることができる。即ち、補正函数のサポートサイズκは、流体と構造の相互作用境界Σ上でＩＭＥ同士が隙間なく重なり合う程度の大きさとする。ここでもやはり、相互作用境界Σが巻き上がるように変形して、表裏が一意に決まらないような運動は想定していない。ｖは通常の有限要素法で用いられる流体の速度であり、要素内でその勾配も連続である。しかしこのように速度を補正することで、全体の速度Ｖとしては、相互作用境界Σ上で連続、且つ、勾配の不連続が実現できる。

次に、弱形式の定義について再び説明する。本実施例では、流体の自由度（又は、従属変数）を補正し、補正函数に構造との相互作用を担わせている。そこで、流体の支配方程式であるNavier-Stokes方程式と構造の運動方程式を基礎方程式とする。

上記の如く、流体構造連成シミュレーションにおいて、流体と構造の相互作用境界に、運動する構造の変位に追随して移動する有限要素を定義し、この要素内に圧力と速度の補正函数を配置し、境界積分を露出させた弱形式を用いて流体構造間の相互作用を記述することで、両領域のメッシュが独立した状態でシミュレーションを行うことが可能となる。これにより、メッシュの不整合問題が解決され、メッシュ制御が不要になる。

本実施例では、メッシュの整合性に依存しない流体構造連成シミュレーションを行うので、次のような特徴（ｋ１）〜（ｋ７）を実現できる。
（ｋ１）メッシュの破綻によってシミュレーションが止まることがない。つまり、構造が流体領域のメッシュ範囲を飛び出す運動をしない限り、メッシュが問題になることはない。
（ｋ２）構造同士が接触するような問題にも適用可能である。
（ｋ３）ＡＬＥ法と比べてメッシュ制御のための自由度が必要無くなるので、演算処理に必要となるメモリ領域を減少させることができる。厳密には、自由境界上での補正函数を表現するための自由度が必要になるが、３次元から２次元に定義域の次元が下がるため、ＡＬＥ法と比べてその自由度は著しく小さくできる。
（ｋ４）メッシュ制御として、弾性体問題を求解する必要がないので、その分処理の高速化が図れる。
（ｋ５）流体メッシュの作成が、ＡＬＥ法と比べてデリケートではないため、メッシュ作成に要する人的、経験的リソースを削減できる。
（ｋ６）ＩＴ型の手法であれば、メッシュの整合性のために、構造が複雑な幾何形状であれば流体部のメッシュも細かくなるが、本実施例では流体のシミュレーションが円滑に進むメッシュを準備するだけで良い。
（ｋ７）既存の流体コードに相互作用ルーチンを付け足すことで、流体構造連成シミュレーションを行なうことができる。その際、既存の流体コードとしては、有限要素法に限らず、差分法、有限体積法等任意の数値シミュレーション手法を用いることができ、流体の圧力と速度が与えられる離散点以外の任意の点での、それら物理量の補間方法が確立されていれば良い。

以上の実施例を含む実施形態に関し、更に以下の付記を開示する。
（付記１）
メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成処理と、
物理現象をシミュレートする主要時間発展ループ処理
をコンピュータにより実行し、
前記主要時間発展ループ処理は、
前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、
前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、
前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する、
流体構造連成シミュレーション方法。
（付記２）
前記主要時間発展ループ処理は、
前記相互作用仲介要素の移動に伴う前記流体の持つ自由度との相互作用ペアの変化に基づいて、前記相互作用仲介要素と重なりを持つバックグラウンドの流体要素をサーチし、
前記流体要素のメンバである流体節点を相互作用ペアとして決定し、
前記相互作用ペアに基づいて前記シミュレーションを実行する、
付記１記載の流体構造連成シミュレーション方法。
（付記３）
前記補正函数は、グリーンの定理に基づいた、前記流体と前記構造の境界付近での前記流体の非圧縮条件及びニュートン流体の粘性作用素の満たす関係を制約条件とするラグランジュ未定乗数である、付記１又は２記載の流体構造連成シミュレーション方法。
（付記４）
前記主要時間発展ループ処理は、
時間の離散化手法に対応して前記構造及び前記流体の夫々の要素毎のマトリクス及び等価節点力を構成する第１の処理と、
前記相互作用仲介要素毎に速度、圧力と補正函数の相互作用を表すマトリクスを構成する第２の処理と、
前記要素毎のマトリクス及び前記相互作用を表すマトリクスに基づいて、要素毎の行列からの寄与を再構成して全自由度に対する行列方程式を作成する第３の処理と、
前記主要時間発展ループ処理においてニュートン・ラフソンループを実行する毎に流体と構造の圧力、速度及び変位を含む従属変数の変化量を求解して評価する第４の処理と、
前記従属変数の変化量とその評価を出力する第５の処理を含む、
付記１乃至３のいずれか１項記載の流体構造連成シミュレーション方法。
（付記５）
前記第２の処理は、
前記圧力の前記補正函数と前記流体の速度の相互作用マトリクスを計算し、
前記速度の前記補正函数と前記流体の速度の相互作用マトリクスを計算し、
前記圧力の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算し、
前記速度の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算する、
付記４記載の流体構造連成シミュレーション方法。
（付記６）
前記補正函数は、前記相互作用仲介要素内でのみ非零の値を持ち、
前記主要時間発展ループ処理は、前記補正函数が、前記流体と構造の相互作用を仲介しながら圧力と速度勾配の不連続性を表現するように、値を持つ領域、前記構造上での函数値、微分値、要素境界の函数値、微分値について夫々に条件を課す、付記１乃至５のいずれか１項記載の流体構造連成シミュレーション方法。
（付記７）
メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成部と、
物理現象をシミュレートする主要時間発展ループ部を備え、
前記主要時間発展ループ部は、
前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、
前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、
前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する、
体構造連成シミュレーション装置。
（付記８）
前記主要時間発展ループ部は、
前記相互作用仲介要素の移動に伴う前記流体の持つ自由度との相互作用ペアの変化に基づいて、前記相互作用仲介要素と重なりを持つバックグラウンドの流体要素をサーチし、
前記流体要素のメンバである流体節点を相互作用ペアとして決定し、前記相互作用ペアに基づいて前記シミュレーションを実行する、
付記７記載の流体構造連成シミュレーション装置。
（付記９）
前記補正函数は、グリーンの定理に基づいた、前記流体と前記構造の境界付近での前記流体の非圧縮条件及びニュートン流体の粘性作用素の満たす関係を制約条件とするラグランジュ未定乗数である、付記７又は８記載の流体構造連成シミュレーション装置。
（付記１０）
前記主要時間発展ループ部は、
時間の離散化手法に対応して前記構造及び前記流体の夫々の要素毎のマトリクス及び等価節点力を構成する第１の部分と、
前記相互作用仲介要素毎に速度、圧力と補正函数の相互作用を表すマトリクスを構成する第２の部分と、
前記要素毎のマトリクス及び前記相互作用を表すマトリクスに基づいて、要素毎の行列からの寄与を再構成して全自由度に対する行列方程式を作成するマトリクスアセンブラと、
前記主要時間発展ループ処理においてニュートン・ラフソンループを実行する毎に流体と構造の圧力、速度及び変位を含む従属変数の変化量を求解して評価するソルバと、
前記従属変数の変化量とその評価を出力する出力部を有する、
付記７乃至９のいずれか１項記載の流体構造連成シミュレーション装置。
（付記１１）
前記第２の部分は、
前記圧力の前記補正函数と前記流体の速度の相互作用マトリクスを計算する第１のサブパートと、
前記速度の前記補正函数と前記流体の速度の相互作用マトリクスを計算する第２のサブパートと、
前記圧力の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算する第３のサブパートと、
前記速度の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算する第４のサブパートを有する、
付記１０記載の流体構造連成シミュレーション装置。
（付記１２）
前記補正函数は、前記相互作用仲介要素内でのみ非零の値を持ち、前記主要時間発展ループ処理は、前記補正函数が、前記流体と構造の相互作用を仲介しながら圧力と速度勾配の不連続性を表現するように、値を持つ領域、前記構造上での函数値、微分値、要素境界の函数値、微分値について夫々に条件を課す、付記７乃至１１のいずれか１項記載の流体構造連成シミュレーション装置。
（付記１３）
コンピュータに流体構造連成シミュレーションを実行させるプログラムであって、
メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成手順と、
物理現象をシミュレートする主要時間発展ループ手順
を前記コンピュータにより実行させ、
前記主要時間発展ループ手順は、
前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、
前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、
前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する、
プログラム。
（付記１４）
前記主要時間発展ループ手順は、
前記相互作用仲介要素の移動に伴う前記流体の持つ自由度との相互作用ペアの変化に基づいて、前記相互作用仲介要素と重なりを持つバックグラウンドの流体要素をサーチし、
前記流体要素のメンバである流体節点を相互作用ペアとして決定し、
前記相互作用ペアに基づいて前記シミュレーションを実行する、
付記１３記載のプログラム。
（付記１５）
前記補正函数は、グリーンの定理に基づいた、前記流体と前記構造の境界付近での前記流体の非圧縮条件及びニュートン流体の粘性作用素の満たす関係を制約条件とするラグランジュ未定乗数である、付記１３又は１４記載のプログラム。
（付記１６）
前記主要時間発展ループ手順は、
時間の離散化手法に対応して前記構造及び前記流体の夫々の要素毎のマトリクス及び等価節点力を構成する第１の処理と、
前記相互作用仲介要素毎に速度、圧力と補正函数の相互作用を表すマトリクスを構成する第２の処理と、
前記要素毎のマトリクス及び前記相互作用を表すマトリクスに基づいて、要素毎の行列からの寄与を再構成して全自由度に対する行列方程式を作成する第３の処理と、
前記主要時間発展ループ処理においてニュートン・ラフソンループを実行する毎に流体と構造の圧力、速度及び変位を含む従属変数の変化量を求解して評価する第４の処理と、
前記従属変数の変化量とその評価を出力する第５の処理
を実行する、
付記１３乃至１５のいずれか１項記載のプログラム。
（付記１７）
前記第２の処理は、
前記圧力の前記補正函数と前記流体の速度の相互作用マトリクスを計算し、
前記速度の前記補正函数と前記流体の速度の相互作用マトリクスを計算し、
前記圧力の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算し、
前記速度の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算する、
付記１６記載のプログラム。
（付記１８）
前記補正函数は、前記相互作用仲介要素内でのみ非零の値を持ち、前記主要時間発展ループ処理は、前記補正函数が、前記流体と構造の相互作用を仲介しながら圧力と速度勾配の不連続性を表現するように、値を持つ領域、前記構造上での函数値、微分値、要素境界の函数値、微分値について夫々に条件を課す、付記１３乃至１７のいずれか１項記載のプログラム。

以上、開示の流体構造連成シミュレーション方法、装置及びプログラムを実施例により説明したが、本発明は上記実施例に限定されるものではなく、本発明の範囲内で種々の変形及び改良が可能であることは言うまでもない。

２０計算機システム
２１ＣＰＵ
２２Ｉ／Ｆ
２３記憶装置
２４入力装置
２５表示装置
２６バス
３０連成処理
３１流体
３２構造
３３ＩＭＥ
４９主要時間発展ループ
６１〜６４サブパート

Claims

メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成処理と、
物理現象をシミュレートする主要時間発展ループ処理
をコンピュータにより実行し、
前記主要時間発展ループ処理は、
前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、
前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、
前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する、
流体構造連成シミュレーション方法。
コンピュータに流体構造連成シミュレーションを実行させるプログラムであって、
メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成手順と、
物理現象をシミュレートする主要時間発展ループ手順
を前記コンピュータにより実行させ、
前記主要時間発展ループ手順は、
前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、
前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、
前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する、
プログラム。
前記主要時間発展ループ手順は、
前記相互作用仲介要素の移動に伴う前記流体の持つ自由度との相互作用ペアの変化に基づいて、前記相互作用仲介要素と重なりを持つバックグラウンドの流体要素をサーチし、
前記流体要素のメンバである流体節点を相互作用ペアとして決定し、前記相互作用ペアに基づいて前記シミュレーションを実行する、
請求項２記載のプログラム。
前記補正函数は、グリーンの定理に基づいた、前記流体と前記構造の境界付近での前記流体の非圧縮条件及びニュートン流体の粘性作用素の満たす関係を制約条件とするラグランジュ未定乗数である、請求項２又は３記載のプログラム。
前記主要時間発展ループ手順は、
時間の離散化手法に対応して前記構造及び前記流体の夫々の要素毎のマトリクス及び等価節点力を構成する第１の処理と、
前記相互作用仲介要素毎に速度、圧力と補正函数の相互作用を表すマトリクスを構成する第２の処理と、
前記要素毎のマトリクス及び前記相互作用を表すマトリクスに基づいて、要素毎の行列からの寄与を再構成して全自由度に対する行列方程式を作成する第３の処理と、
前記主要時間発展ループ処理においてニュートン・ラフソンループを実行する毎に流体と構造の圧力、速度及び変位を含む従属変数の変化量を求解して評価する第４の処理と、
前記従属変数の変化量とその評価を出力する第５の処理
を実行する、
請求項２乃至４のいずれか１項記載のプログラム。
前記第２の処理は、
前記圧力の前記補正函数と前記流体の速度の相互作用マトリクスを計算し、
前記速度の前記補正函数と前記流体の速度の相互作用マトリクスを計算し、
前記圧力の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算し、
前記速度の前記補正函数と前記構造の変位速度の相互作用マトリクスを計算する、
請求項５記載のプログラム。
前記補正函数は、前記相互作用仲介要素内でのみ非零の値を持ち、
前記主要時間発展ループ処理は、前記補正函数が、前記流体と構造の相互作用を仲介しながら圧力と速度勾配の不連続性を表現するように、値を持つ領域、前記構造上での函数値、微分値、要素境界の函数値、微分値について夫々に条件を課す、請求項２乃至６のいずれか１項記載のプログラム。
メッシュで表される流体と構造の夫々について計算領域を離散化した節点間のグラフ情報を構成するグラフ情報構成部と、
物理現象をシミュレートする主要時間発展ループ部を備え、
前記主要時間発展ループ部は、
前記構造の境界上に前記構造の変位に伴って移動する相互作用仲介要素を配置し、
前記相互作用仲介要素内に、前記流体の圧力と速度及び前記構造の変位と相互作用する、前記流体の圧力及び速度の補正函数を定義し、
前記補正函数に基づいて、前記流体と前記構造でメッシュが不整合の状態でシミュレーションを実行する、
流体構造連成シミュレーション装置。