資料:FlexPDE が内蔵する有限要素法

2. Appendix: Principles of Finite Element Analysis (Gunnar Backstrom, Fields of Physics)

有限要素解析 (FEA: finite element analysis) の中心をなす考え方は、対象領域を簡単な形状をした小領域 (セル) に分割し、これら小領域上で同時に問題を解くというものです。FlexPDE では三角形あるいは、断面が三角形の角柱に分割します。境界に接するセルには、少し曲がった外面があります。従って、FlexPDE のスクリーンに表示されるセルは、対象平面上の角柱を投影したものです。

FlexPDE では、離散点 (ノード) つまり、三角形の頂点および頂点間の中点における従属変数の値を決めて PDE を解きます。従って、解プロセスで記録されるものは、これらノードの値のみとなります。三角形セルの他点での関数値および微分係数を得るには、ノード間での内挿をおこないます。


2D における内挿と微分

内挿アルゴリズムは FEA において非常に重要であり、これについて調べてみましょう。これと同時に、微分ルーティンのテストも行います。微分係数がどのように計算されるか予めわかっているわけではありませんが、結果が正確かどうかの判定は簡単に行えます。

簡単な 2変数関数を例にとり、微分係数も含めてプロットし、厳密解と比較します。関数を f (x) = sin (xy) とした場合の記述ファイルを以下に示します。

TITLE 'Interporation of sin(x*y) '
DEFINITIONS
      Lx=pi/2 Ly=pi/2
      f=sin(x*y) fx=dx(f) fx_ex=cos(x*y)*y
      fxy=dx(dy(f)) fxy_ex=sin(x*y)*x*y+cos(x*y)
BOUNDARIES
region 'domain' start(-Lx,-Ly)
      line to (Lx,-Ly) to (Lx, Ly) to (-Lx,Ly) to finish
PLOTS
      contour(f) surface(f) surface(fx)
      contour(fx-fx_ex) surface(fxy) contour(fxy-fxy_ex)
END

fx の誤差の等高線図を見ると、正確にゼロとなっています。これは、数学者が行うのと同じ解析手続きを用いて、数式処理による微分を内部的に行っていることを意味しています。この観察から、FlexPDE は DEFINITIONS で指定された微分、微分関数を厳密な手段によって取り扱っているころを示しています。言い換えれば、解析式を 1回微分しても誤差が発生していない、ということです。

次の図は、混合微分の面をプロットしたもので (surface(fxy))、乖離 fxy-fxy_ex のプロットに示されているように厳密解とは幾分異なっています。


解を得るための内挿

一旦解が得られますと、ノード単位に従属変数の値を格納するテーブルに入れられます。任意点の解を計算する場合には、FlexPDE は、このテーブルを元に二次アルゴリズムによる内挿を行います。この方法は簡単で、各三角形領域ごとに、以下の多項式を使用します。

P(x, y) = a0 + a1x + a2y + a3xy + a4x2 + a5y2

関数値は、3つの頂点、3つの中点、つまり 6点について既知であり、上式の係数の個数に合致します。係数 ai を決めるには、6つの線形方程式を連立して解けばいいことになります。

次に簡単なポアソン方程式を解き、FlexPDE がノードの値だけから、内挿して解を得る様子を観察することにします。任意の関数を選び、そのラプラシアンを取り、その結果の r (x, y) を用いて解として選択する関数を持つ PDE を作ります。方程式 ∇2U = r (x, y) はこのニーズに合致します。境界条件の値は、初期条件から直接得ることができます。記述ファイルは以下のようになります。

TITLE 'Poisson Equation '
SELECT errlim=1e-15 nodelimit=10
VARIABLES       u
DEFINITIONS
            Lx=1 Ly=1 u_ex=y^2*sin(x)
EQUATIONS
            del2(u)=-(y^2-2)*sin(x)
BOUNDARIES
region 'domain' start 'outer' (0,0) value(u)=u_ex
            line to (Lx,0) to (Lx,Ly) to (0,Ly) to finish
PLOTS
            grid(x,y) zoom(0,0,0.3,0.3)
            contour(u-u_ex) zoom(0,0,0.3,0.3) fixed range(-1e-8,1e-8)
END

ノード数を小さくしたいので、nodelimit の指定を行っています。更に fixed range を指定してプロットする関数値をゼロを中心に 1分間隔にするよう限定しています。

厳密解に対する誤差がここでの関心事です。以下のプロット (contour(u-u_ex)) を見ると、全ての等高線が重なり合い殆どゼロになっていることがわかります。等高線は目の粗い誤差曲面の内挿で生成されますので、ノード点を正確に通るとは限りません。下図と同領域のグリッド図を重ねあわしてみると、ゼロ線が頂点と中点を通過していることがわかります。


ノード値を求める

FEA においては、従属変数のノード値を知ることは、問題の解を得るのと同じことです。グリッド上の全ノード点において、従属変数が解っていれば、内挿を行って解領域の任意点での値を求めることができますし、多項式を微分すれば、1次、2次の微分係数を求めることもできます。ノード値を計算する方が、結果から内挿するより困難です。2次の一般的な線形 PDE は以下のように書けます。

ここで、fkxy の関数です。 u(x,y) が PDE の厳密解であれば、Eq は解領域 D 全面にわたってゼロとなります。この領域について積分すると、

しかし、この関係を満たす u(x,y) が必ず解とは限りません。というのは Eq は正であったり、負であったりして、積分内で相殺される場合もあるからです。一方、PDE の二乗をとると、

なる条件が得られます。この場合には、 u(x,y) は解になります。事実、この関係を利用して方程式を数値的に解くことができます。

所与の PDE について、 の上記積分はノード値の関数と見なすことができます。

勿論、数値計算の精度は限定されることを考えると、ノード値 vi の集合が積分 I (u1, u2, ..., un) を正確にゼロにすることは期待できず、積分を最小化することで満足しなくてはなりません。これは多変数関数の最小値を見つけることと等価であり、このような問題を解く標準的な方法があります。しかし、Eq の二乗をとると、被積分関数は複雑かつ非線形になってしまいます。

PDE 自身と境界条件の値が線形であれば、線形解析しか必要としない代替方法もあります。つまり、問題を上記のように単一の積分形式にしないで、上記と同じタイプの一群の方程式を導入します(j=1,2,...m)。

これらの方程式は二乗されていません。Wj(x,y) はウエイト関数で、全解領域 D 内部の特定の小領域 (セル) を強調することができます。ウエイト関数は自由に選択できますが、計算の観点からすれば、最も単純な方法は、ウエイトを 0 か 1 かにすることです。たとえば、最初の方程式では、三角形セルの 1つに 1のウエイトを用い、他全てのウエイトをゼロとします。要は試行関数 (最終的には解) が全てのセル上で、PDE を満足するようにするということです。

上述方法を適用する1つのやり方は、セル番号 j のウエイトを1、その他をゼロとするもので、その結果、解領域の各部分に関して m 個の方程式が得られます。この場合、セル数 (m) をノード数 (n) よりはるかに小さくする必要があります。

今、従属変数 (u) は1つだけで、境界条件全ては value で指定するものとします。従って、境界上の全ての値は既知であり、内部ノードの値だけを求めればよいことになります。その数は一般的いセル数より大きいので、方程式の数を増加させる柔軟な方法が必要になります。

ウエイト関数の選択については、幾つかの方法が知られています(5)。正しい方程式の数を得る一群のルールは以下の通りです。ノード j が中点である場合、この点を共有する2つのセルについて Wj=1 とし、それ以外の場合は Wj=0 とします。ノード j が頂点の場合には、この頂点を共有するセルについて Wj=1 とし、それ以外の場合は Wj=0 とします。点が境界上にある場合には、1つのセルだけを使用します。この手続きによって正しい方程式の数が得られ、全ての小領域の積分は異なるので、独立した方程式が得られます。

従属変数が複数個ある場合には、もっと高度なウエイト方式が必要になります。FlexPDE でどのような方法を用いているかは、明らかではありません。

PDE および境界条件が線形である場合、解析結果は線形方程式系になります。PDE もしくは境界条件が非線形である場合には、より一般的な代数方程式が得られます。FlexPDE は非線形系を反復法で解いており、一般的には計算に時間がかかり、予測できない場合もあります。


自然境界条件

これまで、全境界上で解を value で与えるものとしてきました。全境界あるいはその一部で自然条件である場合には、その情報は特別な方法で内蔵しておく必要があります。PDE の積分と場 F の外向き法線成分を関係づける定理は、

PDE を小領域上で積分すると、ラプラシアンの体積積分を自然境界条件が関与する面積分に置き換えることができます。セルの他 2 面は、隣接セルに合致すべき積分となります。

自然境界条件がある場合にはいつでも、PDE にあるファクタを乗じると体積積分を変更でき、境界上の面積分にできることを記憶しておくのが大切です。従って、ファクタ f を乗じると、これに対応する自然境界条件にも同じファクタが乗じられることになります。


  • (5) Bickford, W.B., A First Course in the Finite Element Method, Richard D. Irwin, 1990