これは数値計算Advent Calendar 2018の12月12日の記事ですが, 公開日は2018年12月26日です.
遅くなってしまい大変申し訳ありません.
みなさん, IGA(Isogeometric Analysis)はご存知でしょうか.
IGAは日本語の情報が少なく, たぶん知らない方も多いなので, 今回の記事ではざっくりと解説することを目標にします.
ざっくり, 本当にざっくり言えばIGAは「数値計算・ミーツ・形状表現」です.
ここでの数値計算とはFEM(Finite Element Method, 有限要素法)でも使われているGalerkin法のことで, 形状表現とはNURBS(Non-uniform Rational B-spline)のことです.
通常のFEMでは, 近似解を構成するために解析対象の形状をメッシュ分割します.
それに対してIGAでは三角形分割のようなメッシュ分割は行いません.
これを可能にしているのがNURBSで, 近似解構成のための基底関数が形状表現から自然に誘導される事を利用しています.
IGAで解ける問題は幾つかありますが, 本記事では2次元Poisson方程式を例として解説します.
領域がNURBS多様体で表されているとします.
次の2次元Poisson方程式を考えます.
ここで, は与えられた上の関数で, が求めるべき未知関数です.
はDirichlet境界, はNeumann境界でを満たします.
とくにであればLaplace方程式と呼ばれ, その解は調和関数と呼ばれます.
定常の熱伝導方程式と思えば, が温度, が上の熱源, が境界上の温度, は上での断熱条件に相当します.
IGAの手順を書けば次のようになります:
通常のFEMとは基底関数の作り方が違うだけなので, 計算手順も基本的に同じになります.
Galerkin法のために弱形式に書き換えます.
上でを取るテスト関数を掛けてで積分します.
ここでは外向き単位法線ベクトルで, 記号は次のように定めました.
つまり, 任意のテスト関数に対して
が充たされるようなを求めれば良いということになります.
これを弱形式と言います.
これに対して元の偏微分方程式は強形式と呼ばれ, 弱形式と(数学的な厳密さを言わなければ)等価です.
NURBSの詳細はNURBS多様体による形状表現pdfに書いたのでそちらを参照してください.
今回使う性質で重要なものは次の2点です.
領域の具体例として次の図のような円環領域の半分を考えましょう.
この形状は次のようにして作れます.
図には書いてませんが, 内半径1, 外半径2としています.
は矩形領域でです.
IGAにおける基底関数という用語には次の2つの意味があるため注意が必要です.
Galerkin法で使う基底関数はNURBS基底関数から誘導されるため実質的にこれらは同じで, 混乱の恐れはありません.
(はのために2つの添字が走りますが, Galerkin法での基底関数は添字一つで十分なので束ねてと書いています)
ただしNURBSによる形状表現が決まった段階で基底関数(と制御点)の数は決まってしまい, このままでは近似解を得るためには不十分という問題があります.
つまり近似解の精度向上には基底関数を増やす必要があり, ここで細分(refinement)と呼ばれる操作が役立ちます.
細分はパラメータ付けを変えることなく, 基底関数と制御点の数を増やす操作のことです.
(こちらも詳しくはNURBS多様体による形状表現pdfを参照して下さい)
に対して細分を行って制御点を増やした例が次の図です.
この細分の操作は, 有限要素法で言えばメッシュを細かくする操作に相当します.
理論的背景は他の文献を当ってもらうとして, ここでは解法のみ簡単に述べます.
をテスト関数とした弱形式は次のように表されるのでした.
NURBS基底関数を用いて, 厳密解をで近似します:
ここでは係数です.
はDirichlet条件から決まる添字の集合で, 次を満たします.
はを上で近似できるように係数を決定して構成されます. (の決定は難しくなく, とくに今回の例では定数のみを考えるのでより簡単です)
なので, あとはを求めればよいことになります.
も同様にのようにして離散化します.
これらを先の弱形式に代入すれば, を未知変数とする連立方程式()に帰着されます.
積分はNURBS多様体で入る座標で行うため, 次のようにして変数を取り替えます.
簡単のためにとし, 次のようにを定めます.
この条件で数値計算した結果が次のグラフです.
厳密解が調和関数になることを思い出せば, これはの実部と同じ形になることが分かります.
つまりは対数関数のグラフを原点中心に回転させたものに一致し, 次の図はその断面です.
簡単のためにとし, 次のようにを定めます.
この条件で数値計算した結果が次のグラフです.
厳密解が調和関数になることを思い出せば, これはの虚部と同じ形(常螺旋面)になることが分かります.
全てDirichlet条件としてとすれば次のようになります.
この条件で解いた結果が次のグラフです.
この計算結果では, 中央付近に不自然な山があるように見えます.
これは基底関数の数が十分でないために起こっており, もう少し細かい細分を取ると改善できます.
新しい基底関数で計算した結果が次の図です.
不自然な山は見える範囲では現れず, 近似解の精度が上がっている事が確認できます.
これまでの例と異なり, この条件に対する厳密解を陽に書き下すのは難しそうです.
前述のようにこれも物理的解釈が出来て, 「全体に熱源があって境界での温度がで固定された定常状態」というものです.
に対しても厳密解()が得られるような形状(単位円板)・境界条件として次を与えましょう.
この領域のNURBSによる形状表現は次のようになります.
基底関数は次のようになります.
この条件で解いた結果が次のグラフです.
厳密解と近似解の差を取れば次のようになります.
ただし差を大きく見せるため, 高さ方向に200倍して表示しています
細分を取って基底関数の数を増やすと次のようになります.
これらの基底関数を使って近似解を構成すると次のようになりました.
の見た目だけでは細分前後の変化が分かり難いですが, 同じように差を取ると誤差が減っていることが確認できます.
isoは「等しい」, geometricは「幾何的に」という意味なので, IGA(isogeometric analysis)は直訳すれば「幾何的に等しい解析」ということになります.
通常の有限要素法では, 領域を三角形分割する必要があるため, その意味で形状に対する近似が入ります.
これに対してIGAでは, 今回見たようにがNURBSで表現されていれば形状に対する近似無しに数値解を得ることが出来ます.
これが「幾何的に等しい解析」の意味する所です.
CADなどで使われる形状の多くはNURBSで表現可能であり, その意味でNURBSは非常に柔軟に様々な形状を表すことが可能です.
IGAではNURBSで表現可能な形状しか扱えませんが, これはデメリットでは無く, 寧ろCADで作成した形状からシームレスに数値計算を行えるためにメリットと考える方が適切です.
今回の数値計算のためのコードはJuliaで書きました.
全然綺麗には書けてませんが, 一応githubに上げておきます.
https://github.com/hyrodium/IGA-PoissonEquation
数学的にはGalerkin法を使っているため, 有限要素法に関する文献も有用です.
QiitaではなくHackMDで書いたため, コメントなどが出来ません.
質問・間違いなどありましたら筆者のTwitterまでご連絡ください.