有限要素法の基礎①：一次元弾性体、近似解、ガラーキン法【MATLABコード付き】

はじめに

　有限要素法の定式化は、最初から要素剛性行列を置くだけでは全体像をつかみにくいものです。そこで本記事では、まず一次元弾性棒の釣合い条件、適合条件、構成式からナビアの式を導きます。その後で、得られたナビアの式を一自由度ばねの \(F=kx\) と対応づけ、有限要素法の全体方程式へのつながりを示します。さらに、近似解、重み付き残差法、ガラーキン法へ進みます。

有限要素法とは

　最初に理解すべきは、有限要素法は微分方程式を解く道具であるということです。有限要素法は釣合い式や浸透流の連続の式を解くのが直接の目的ではありません。もちろん、解くべき微分方程式が、たまたま釣合い式や連続の式かもしれませんが、あくまでも直接に扱うのは微分方程式です。いわば、数学の道具です。

　有限要素法と聞いて、始めから離散的なトラス構造を思い浮かべるのは適切な連想ではありません。したがって、有限要素法が相手とする微分方程式をまず知らなければなりません。しかし、微分方程式ならなんでもよいというわけではありません。なぜなら、勝手に作った微分方程式にはどのような境界条件のもとで一つの解が存在するかがわかっていませんし、また、それに有限要素法を適用した場合、どのような定式化がふさわしいのかも明確でないからです。

　逆にいえば、有限要素法が最も得意とする微分方程式があります。それが楕円型、あるいはそれに近い型の微分方程式です。物理現象でいえば、変位を未知関数とする弾性体の釣合い式や、全水頭を未知関数とする浸透の連続の式です。ここでは、このような微分方程式の誘導方法や性質について述べたあと、いかにしてそれらに有限要素法を適用するのかを述べます。

　有限要素法の基礎は意外と分かりやすい概念に立脚しています。重要なキーワードは、弱形式、ガラーキン法、形状関数（基底関数）の三つです。これらのうち、弱形式はしばしば「仮想仕事の式」と呼ばれているものであり、有限要素法による解の意味を考える場合に、極めて大きな役割を果たしています。

一次元弾性体の基礎式

棒の基礎式を順に導く

　まず、ナビアの式を先取りせず、棒の力学から順に導きます。物体の変形を考察する静力学の解が満足しなければならないのは、以下の三つの関係式と境界条件です。

釣合い条件
適合条件
構成式
境界条件

　このうち、釣合い条件、適合条件、構成式の関係を先に図1へ示します。適合条件と構成式を組み合わせると応力を変位で表せます。さらに、その結果を釣合い条件へ代入すると、変位だけを未知量とするナビアの式が得られます。以下では、この関係を一つずつ導出します。

*図1 釣合い条件、適合条件、構成式からナビアの式を導く関係です。適合条件と構成式から応力を変位で表し、それを釣合い条件へ代入します。*

釣合い条件：棒の各位置における力のつり合い

　いま、図2のように、長さが \(l\) で両端を固定された断面積 \(A\) の棒が \(x\) 軸に沿っておかれています。棒の左端は \(x=0\) に、また、右端は \(x=l\) にあります。この棒には \(x\) 軸方向に単位体積当たり物体力 \(f(x)\) が作用します。棒の断面内に作用する引張り応力を \(\sigma\) とします。特に座標 \(x\) における \(\sigma\) の値を関数 \(\sigma(x)\) のように表示します。

図2 一次元弾性棒と微小区間の釣合い。（a）は両端を固定された断面積 \(A\) の棒に物体力 \(f(x)\) が作用する問題、（b）は棒から切り出した微小区間に作用する応力と物体力を示します。

\begin{align} (\sigma+d\sigma)A-\sigma A+fA\,dx=0 \tag{1} \end{align}

　あるいは、

\begin{align} \frac{d\sigma}{dx}+f=0 \tag{2} \end{align}

　のようになります。これは、応力で表現した釣合い条件です。

適合条件

　次に、変位とひずみの関係を述べます。一般に力学を考えるとき、力の釣合い条件ばかりが注目され、ややもすれば忘れがちなのが、この変位とひずみの関係です。その一つの理由は、力の釣合いが直感的に把握できることに比べて、変位–ひずみ関係がやや数学的あるいは幾何学的であるからとも考えられます。

　図3のように、直線棒の軸に沿って \(x\) 軸を設けたとき、座標 \(x\) の点において \(x\) 軸方向に \(u(x)\) だけ変位したとします。\(x\)（棒片の尻尾）と \(x+dx\)（棒片の頭）の範囲を取り上げます。変形前には棒片の長さは \(dx\) です。変形に伴って \(x+dx\) においては \(u(x+dx)\) だけ \(x\) 軸方向に移動し、\(x\) においては \(u(x)\) だけ同方向に移動します。したがって、その差である \(u(x+dx)-u(x)\) だけ棒片は伸びることになります。

図3 変位とひずみの関係。変形前に長さ \(dx\) であった棒片の両端が、それぞれ \(u(x)\)、\(u(x+dx)\) だけ移動する様子を示します。

　それを \(dx\) で割ることによって、その区間における平均的な伸び率が求められます。これがひずみ \(\varepsilon\) であり、

\begin{align} \varepsilon=\frac{du}{dx} \tag{3} \end{align}

となります。この式は変位–ひずみ関係ですが、ここではこれを適合条件と呼ぶことにします。

構成式

　いま、簡単のため棒は線形弾性体からなるとし、そのヤング（Young）率を \(E\) とします。\(x\) 軸方向の引張りひずみ \(\varepsilon\) と応力 \(\sigma\) とは、

\begin{align} \sigma=E\varepsilon \tag{4} \end{align}

　なる関係があります。これは最も簡単な応力–ひずみ関係ですが、一般的には構成式と呼ばれるものです。材料の特性を記述する式です。この弾性モデルがすべてのモデルの基本となっています。また、有限要素法を理解するうえではこれで十分です。

　以上で、境界条件以外の三つの条件を簡単な例を用いて説明しました。すると以下に示すように変位 \(u\) で表現した釣合い式を導くことができます。まず、構成式\((4)\)に適合条件式\((3)\)を代入すると、

\begin{align} \sigma=E\frac{du}{dx} \tag{5} \end{align}

　のように応力を変位で表すことができます。さらに、この結果を釣合い式\((2)\)に代入すると、変位で表した釣合い式、

\begin{align} E\frac{d^2u}{dx^2}+f=0 \tag{6} \end{align}

　が導かれます。このように変位で表した釣合い式をナビア（Navier）の式といいます。ナビアの式は、釣合い条件、適合条件、構成式を組み合わせた結果です。

ナビアの式と単純ばねの式\(F=kx\) の類似性

　ここで、得られた式\((6)\)を、一つのばねに対する \(F=kx\) と比較します。ばねの変位を \(q\)、ばね定数を \(k\)、外力を \(F_{\mathrm{ext}}\) とすると、ばねが変位 \(q\) に応じて生む復元力の大きさは \(kq\) です。静止している状態では、外力と復元力がつり合うため、

\begin{align} F_{\mathrm{ext}}-kq=0 \qquad\Longleftrightarrow\qquad kq=F_{\mathrm{ext}} \tag{7} \end{align}

となります。普段用いる \(F=kx\) は、内力を考えていない式ではありません。右辺の外力 \(F_{\mathrm{ext}}\) に対して、左辺の \(kq\) がばね内部に生じる復元力を表しており、式\((7)\)にはすでに「外力と内力のつり合い」が含まれています。

　棒でも基本的な考え方は同じです。ただし、ばねでは一つの変位 \(q\) を求めればよかったのに対し、棒では位置ごとに変位が異なるため、未知量は変位関数 \(u(x)\) になります。そこで棒から長さ \(dx\) の小さな区間を取り出します。

　この区間の伸びを \(\Delta u\) とすると、ひずみは \(\Delta u/dx\)、応力は \(E\Delta u/dx\)、軸力は応力に断面積 \(A\) を掛けた \(EA\Delta u/dx\) です。したがって、「軸力＝剛性×伸び」の形で見れば、この区間の軸方向剛性は \(EA/dx\) となります。実際の区間の伸びは \(u(x+dx)-u(x)\) なので、この区間が生む軸力 \(N(x)\) は、ばねの \(kq\) と同じ形で、

\begin{align} N(x) = \frac{EA}{dx} \left\{ u(x+dx)-u(x) \right\} \longrightarrow EA\frac{du}{dx} \tag{8} \end{align}

　と表されます。つまり、ばねの復元力 \(kq\) に対応する量が、棒では軸力 \(N=EA du/dx\) です。

　ここで大切なのは、棒の軸力が位置によって変わり得ることです。棒の各部分に単位体積当たりの物体力 \(f(x)\) が作用する場合、微小区間の左面と右面に作用する軸力の差が、その区間に作用する外力 \(Af(x) dx\) とつり合います。したがって、連続体における式\((7)\)に相当する式は、

\begin{align} \frac{dN}{dx}+Af(x)=0 \tag{9} \end{align}

です。式\((8)\)を式\((9)\)へ代入すると、

\begin{align} \frac{d}{dx} \left( EA\frac{du}{dx} \right) +Af(x) =0 \tag{10} \end{align}

となります。式\((10)\)は一次元弾性棒のナビアの式の一般形です。\(E\) と \(A\) が一定なら、式\((10)\)を断面積 \(A\) で割ることにより、すでに導出した式\((6)\)と一致します。この比較により、ナビアの式が \(F=kx\) と異なる発想ではなく、位置ごとに変位が異なる連続体について同じ力のつり合いを表した式であることが分かります。

　さらに先の見通しを示すと、有限要素法では棒を有限個の要素と節点に分け、連続的な変位関数 \(u(x)\) を有限個の節点変位ベクトル \(\mathbf{u}\) で近似します。その結果、式\((10)\)は、

\begin{align} \mathbf{K}\mathbf{u}=\mathbf{F} \tag{11} \end{align}

という連立方程式になります。これは式\((7)\)の \(kq=F_{\mathrm{ext}}\) を多自由度に拡張したものです。

　棒の軸に沿う弾性変形問題の支配方程式は2階の微分方程式となります。これは、代表的な一次元の楕円型微分方程式です。式\((6)\)の一般解には二つの積分定数が含まれますが、これを決定するのに必要なものが、以下に述べる境界条件です。

境界条件

　図2に示したように、ここで考える直線の棒は両端で移動しないように \(x\) 軸方向に固定されているとします。すなわち、

\begin{align} u(x=0)=0 \tag{12} \end{align}

\begin{align} u(x=l)=0 \tag{13} \end{align}

です。このように変位、すなわち未知量そのものを指定する境界条件を変位境界条件といいます。あるいはディリクレ（Dirichlet）条件といいます。

　ここで使う「斉次（せいじ）」と「非斉次（ひせいじ）」は、境界で指定する値が0か、0でないかを区別する数学上の用語です。材料が均一であるかどうかを表す言葉ではありません。一般に、未知関数 \(u\) に作用する境界条件を \(B[u]\) と書いたとき、\(B[u]=0\) の形を斉次境界条件、\(B[u]=g\) で \(g\neq0\) の形を非斉次境界条件と呼びます。

　変位を指定するディリクレ条件では、たとえば \(u(0)=0\) は斉次条件です。これは棒の端を基準位置から動かないように固定することを表します。一方、\(u(0)=u_0\) で \(u_0\neq0\) と指定すれば非斉次条件です。これは棒の端をあらかじめ \(u_0\) だけ強制的に移動させることを表します。

　応力や変位勾配を指定するノイマン条件でも考え方は同じです。右端で \(du/dx=0\) とする条件は、線形弾性棒では応力が0、すなわち自由端であることに対応する斉次条件です。これに対して、右端へ0でない応力や軸力を与える条件は非斉次ノイマン条件です。

　斉次条件が数値解析で扱いやすい理由は、斉次条件を満たす関数どうしを足したり定数倍したりしても、境界値が0のまま保たれるからです。後で近似解を基底関数の和として表すときには、この性質を利用します。非斉次条件の場合は、まず指定された0でない境界値を表す既知関数を用意し、その上に斉次条件を満たす未知の近似成分を重ねます。

　一方、問題によっては変位ではなく、境界に作用する応力、すなわち未知量の1階微分に相当する量を指定するような境界条件を与えることもあります。これを応力境界条件、あるいはノイマン（Neumann）条件といいます。したがって、たとえば右端が自由端であるというのは、

\begin{align} \frac{du}{dx}(x=l)=0 \tag{14} \end{align}

という斉次応力境界条件ということができます。

　与えられた領域のすべての境界において、変位あるいは応力境界条件のいずれかを指定しなければ問題の解は一意的に定まりません。ただし、同一の境界点において変位境界条件と応力境界条件を同時に与えることは不可能です。すなわち、境界で変位を指定すれば、そこでの応力は問題を解くことによって自動的に定まり、また、応力を指定すれば、そこでの変位は問題を解くことによって自動的に与えられるからです。

　一部が変位境界条件で、残りが応力境界条件であってもよいです。なお、応力境界条件ばかりを与える場合には、全体で釣合い条件を満足するような境界応力の設定が必要です。たとえば、ここで取り上げている例題のように、物体力の作用している棒の問題においては、物体力の総和に釣合うだけの力を両端の境界力として与えなければなりません。これを満たす場合においても、棒を左右に平行移動することは可能なので、絶対的な変位の値は不定となります。

　式\((6)\)の物体力を \(f(x)=f_0\) のように一定の大きさとすると、この微分方程式は容易に解けて、

\begin{align} u(x) = -\frac{f_0}{2E} \left( x^2+C_1x+C_2 \right) \tag{15} \end{align}

となります。ここに \(C_1,C_2\) は未知の積分定数です。図2の境界条件に合わせて積分定数を決定すれば、

\begin{align} u(x)=\frac{f_0}{2E}x(l-x) \tag{16} \end{align}

となります。

この解が表す変形を直感的に読む

　ここで一定としている \(f_0\) は、棒の右端だけに加える集中荷重ではありません。図2（a）の橙色の矢印で示したように、棒を構成する各部分に同じ向きで作用する、単位体積当たりの一様な物体力です。重力を横向きにしたような荷重をイメージすると、端部荷重との違いをつかみやすくなります。さらに、この棒は左端だけでなく右端も固定されているため、両端の変位はともに0です。

　式\((16)\)は下に凸ではなく上に凸の放物線であり、両端で0、中央で最大になります。この形を理解するには、変位そのものだけでなく、その傾きであるひずみと応力を見ると明確です。式\((16)\)を \(x\) で微分し、式\((5)\)を用いると、

\begin{align} \varepsilon(x) = \frac{du}{dx} = \frac{f_0}{2E}(l-2x), \qquad \sigma(x) = E\varepsilon(x) = \frac{f_0}{2}(l-2x) \tag{17} \end{align}

となります。左半分の \(0<x<l/2\) では \(du/dx>0\) なので、右へ進むほど変位が大きくなり、棒は引張り状態です。反対に、右半分の \(l/2<x<l\) では \(du/dx<0\) なので、右端の固定点へ近づくほど変位が小さくなり、棒は圧縮状態です。中央では左右の状態が切り替わり、

\begin{align} \left.\frac{du}{dx}\right|_{x=l/2}=0, \qquad u_{\max} = u\left(\frac{l}{2}\right) = \frac{f_0l^2}{8E} \tag{18} \end{align}

となります。つまり、中央で変位が大きいのは、中央に特別に大きな力が加わっているからではありません。左端から中央までは正のひずみが積み重なって変位が増え、中央から右端までは負のひずみが積み重なって変位が0へ戻るためです。変位が最大になる中央では変位の傾きが0なので、その位置のひずみと応力も0です。図4の右図は、この応力が左端から右端へ直線的に変化する様子を示しています。

　なお、「左端を固定し、右端に一定の集中荷重 \(P\) を加える棒」を考える場合は別の問題です。その場合、物体力がなければ棒の軸力と応力は一定となり、変位は固定端の0から荷重端へ向かって直線的に増加します。今回の放物線状の変位は、両端固定と一様な物体力を組み合わせた結果です。

図4 一様な物体力を受ける両端固定棒の解。左図は式\((16)\)の厳密解と、重み付き残差法およびガラーキン法などから得られる近似解の比較です。右図は式\((17)\)の応力分布であり、左半分が引張り、右半分が圧縮、中央ではひずみと応力が0になることを示します。

近似解の意味

　一般に微分方程式と適切な境界条件が与えられたとき、それに一切の近似や仮定を設けず、主として積分演算を用いながら求められる解を解析解といいます。解析解のなかには、フーリエ級数や多項式展開のように無限項の和として表現されるものもあります。ここで扱った例題では、弾性係数 \(E\) や物体力 \(f\) が初等関数である場合には、式\((16)\)のように解析解を求めることができます。

　しかし一次元問題でも、少し複雑な条件を伴ったりする場合や、二次元あるいは三次元の問題のほとんどの場合には、まず解析解は求められません。むしろ、求められるのが例外的です。そのような場合に近似的に解を求めるのですが、そのとき、数値的に求められる近似解を数値解といい、それを求める方法を数値解法といいます。

　さて、はたして近似解の「近似」とはどのような意味なのでしょうか。じつは、それにはいろいろな答え方があります。さらにいえば、厳密解ですらその定義の仕方によっては、解が存在しなかったり、あるいは存在しても一つでないこともあります。たとえば、弾性係数 \(E\) や物体力 \(f\) が特異な関数である場合には、そのような問題が生じます。むしろ逆に、「どのように解析解を定義すれば、その解析解がただ一つ存在するか」という問いかけもありえます。

　まして近似解には、いろいろな定義があっても不思議ではなく、また、それにより近似解が異なることも当然です。ところで一般に解析解（厳密解）には種々の属性があります。たとえば、

境界条件を満足する。
必要な階数だけ微分が可能である。
支配方程式をいたるところ満足する。

　などの性質を有しています。これらはすべて、解析解であるための必要条件です。つまり、解析解であるならばこれらすべてを満足するはずです。

　近似解というのは、これらのうち幾つかを弱くした条件を満足するものですが、その弱くさせる方法はもちろん一意的ではありません。しかし多くの場合、境界条件のうちディリクレ条件は常に要求します。したがって、「支配方程式など無視して、ディリクレ条件のみを満たすものをすべて解と定義する」というのが最も広い近似解の定義と思われます。

　いうまでもなく、このように乱暴な宣言をすれば解は容易に求められますが、それにどのような意味があるかは、それを宣言する人の主観的な価値判断によります。逆に狭い範囲で定義された近似解は求めにくいですが、それだけ多くの価値を有しています。詳しい議論はさておき、数値計算で求めた解がどのような意味での解であるかは、ある程度理解しておく必要があります。

　たとえば、有限要素法による解は、ディリクレ条件を満たすものの、必要な階数の微分可能性はまったく満足しませんし、支配方程式も弱い意味でのみ満足するような解です。

近似解の原型

　微分方程式の解は一つのなめらかな関数です。関数を正確に表現するためには、無限の多くの点におけるその値を決定しなければなりません。それは実際上、数値計算では不可能であるので、未知関数を有限の未知数で表現することを考えます。

　簡単のため、棒の長さ \(l=1\)、そして式\((6)\)で \(E=1\) と仮定した微分方程式、

\begin{align} \frac{d^2u(x)}{dx^2}+f(x)=0 \tag{19} \end{align}

　と斉次ディリクレ境界条件、

\begin{align} u(x=0)=0 \tag{20} \end{align}

\begin{align} u(x=1)=0 \tag{21} \end{align}

　を取りあげます。近似関数の一つの自然な表現は未知関数 \(u(x)\) を以下のように有限個の異なる既知関数 \(g_i(x)\) の線形和で表す方法です。

\begin{align} \tilde{u}(x) = \sum_{i=1}^{N}a_i g_i(x) \tag{22} \end{align}

　ここに、\(g_i(x)\ (i=1,2,\ldots,N)\) は、解 \(u(x)\) と同じ斉次境界条件、

\begin{align} g_i(0)=g_i(1)=0 \tag{23} \end{align}

　を満たす \(N\) 個の既知のなめらかな関数です。この \(g_i(x)\) の選び方により解の近似の度合が異なるかもしれませんが、基本的には自由に選ぶことができます。また、個々の \(g_i(x)\) が斉次の境界条件を満足することから、いかなる未知係数 \(a_i\) に対しても、近似関数 \(\tilde{u}(x)\) はやはり斉次境界条件を満たします。

　もしも、非斉次のディリクレ境界条件が与えられている場合には、その非斉次の境界条件を満たす適当な既知関数 \(g(x)\) を式\((22)\)の右辺に加えて、

\begin{align} \tilde{u}(x) = g(x) + \sum_{i=1}^{N}a_i g_i(x) \tag{24} \end{align}

　とすれば、やはり、いかなる未知係数 \(a_i\) に対しても、近似関数 \(\tilde{u}(x)\) はその非斉次境界条件を満たします。ノイマン条件の場合も同様にすればよいです。

　未知関数 \(u(x)\) の近似である式\((22)\)や式\((24)\)で重要なことは、未知関数 \(u(x)\) 自体が \(N\) 個の未知数に置き換わっていることです。だから、いったんこのような近似式を用いれば、あとは何らかの方法で \(N\) 個の方程式を見つければよいことになります。以下で、その代表的な方法を述べます。なお、ここでは \(g_i(x)\) のことを基底関数と呼ぶことにします。

重み付き残差法

　式\((22)\)の \(\tilde{u}(x)\) が厳密解であると仮定すると、このとき任意の関数 \(v(x)\) に対して、

\begin{align} \int_0^1 \left( \frac{d^2\tilde{u}(x)}{dx^2} +f(x) \right) v(x)\,dx =0 \tag{25} \end{align}

　すなわち、

\begin{align} \int_0^1 \left( \sum_{i=1}^{N} a_i\frac{d^2g_i(x)}{dx^2} +f(x) \right) v(x)\,dx =0 \tag{26} \end{align}

　が成立しなければなりません。すなわち、解ならば各点ごとに式\((19)\)を満たすことから、当然、上式の積分値は0でなければなりません。つまり、この式は \(\tilde{u}(x)\) が厳密解であるための一つの必要条件です。

　ただ、近似解は各点ごとで式\((19)\)を満たすのは困難であるので、式\((25)\)や式\((26)\)のように、重み付きの平均的な意味で支配方程式を満たすことで我慢をします。いわば、方程式の誤差（残差）を重み付きの意味で0にしようとする方法で、これを重み付き残差法といいます。

　もちろん、なるだけ多くの数の \(v(x)\) について上式を満足すれば、より精度のよい近似解といえますし、わずかの数の \(v(x)\) についてのみしか成立しないのであれば精度が悪いといえます。このような意味で関数 \(v(x)\) のことを試験関数と呼ぶことがあります。

　ところでいま、近似解 \(\tilde{u}(x)\) を決定するには、\(N\) 個の未知係数を決定すればよいのであるから、いろいろな \(v(x)\) を \(N\) 個、\(v_1(x),v_2(x),\ldots,v_N(x)\) 選び、式\((26)\)に代入して \(N\) 個の方程式を作成すれば \(a_i\) を決定する以下のような連立方程式が得られます。

\begin{align} \begin{bmatrix} \displaystyle\int_0^1 v_1(x)g_1”(x)\,dx & \cdots & \displaystyle\int_0^1 v_1(x)g_N”(x)\,dx \\ \displaystyle\int_0^1 v_2(x)g_1”(x)\,dx & \cdots & \displaystyle\int_0^1 v_2(x)g_N”(x)\,dx \\ \vdots & \ddots & \vdots \\ \displaystyle\int_0^1 v_N(x)g_1”(x)\,dx & \cdots & \displaystyle\int_0^1 v_N(x)g_N”(x)\,dx \end{bmatrix} \begin{Bmatrix} a_1\\a_2\\\vdots\\a_N \end{Bmatrix} = \begin{Bmatrix} \displaystyle-\int_0^1 v_1(x)f(x)\,dx\\ \displaystyle-\int_0^1 v_2(x)f(x)\,dx\\ \vdots\\ \displaystyle-\int_0^1 v_N(x)f(x)\,dx \end{Bmatrix} \tag{27} \end{align}

　なお、\(g”(x)\) は関数 \(g(x)\) の2階微分を指します。これを解けば式\((22)\)より \(\tilde{u}(x)\) が得られます。できれば \(N\) 個以上の \(v(x)\) に対しても、上式は満足すべきではありますが、一般にはそれは不可能です。つまり、本当の解 \(u(x)\) を有限個の自由度で近似していることの帰結です。

　したがって、よい近似解を得るには、正解をうまく近似するような \(g_i(x)\) を選ぶように努力しなければなりません。ここでは、\(g_i(x)\) として2階微分可能な関数を仮定しているから、近似解は境界条件のほかに、微分可能性も満たしています。ただ、支配方程式\((19)\)を厳密に満足しているのではなく、式\((25)\)あるいは式\((26)\)のように領域で平均的に満たしているにすぎません。したがって、支配方程式の満足度が弱められています。

例題1

　式\((22)\)において \(N=1\) として、

\begin{align} g_1(x)=\sin(\pi x) \tag{28} \end{align}

　を採用します。これは明らかに境界条件である式\((23)\)を満たします。任意に選ぶことのできる一つの試験関数を、

\begin{align} v(x)=1 \tag{29} \end{align}

　とすれば、未知数 \(a_1\) に関する方程式\((27)\)は、

\begin{align} -a_1\pi^2 \int_0^1\sin(\pi x)\,dx = -\int_0^1 f(x)\,dx \tag{30} \end{align}

となります。さらに \(f(x)=1\) とすれば、\(a_1=1/(2\pi)\) が得られます。よってこの場合の近似解は、

\begin{align} \tilde{u}(x) = \frac{1}{2\pi}\sin(\pi x) \tag{31} \end{align}

となります。式\((16)\)で表される厳密解と比較すると図4に示すようになります。

ガラーキン法

　重み付き残差法の基本式\((25)\)あるいは式\((26)\)において \(v_i(x)=g_i(x)\) とする方法をガラーキン法といいます。つまり、解 \(u(x)\) を近似するのに用いた関数、すなわち基底関数 \(g_i(x)\) そのものを試験関数とするような重み付き残差法です。式\((27)\)から明らかなように、解くべき連立方程式は以下のようになります。

\begin{align} \begin{bmatrix} \displaystyle\int_0^1 g_1(x)g_1”(x)\,dx & \cdots & \displaystyle\int_0^1 g_1(x)g_N”(x)\,dx \\ \displaystyle\int_0^1 g_2(x)g_1”(x)\,dx & \cdots & \displaystyle\int_0^1 g_2(x)g_N”(x)\,dx \\ \vdots & \ddots & \vdots \\ \displaystyle\int_0^1 g_N(x)g_1”(x)\,dx & \cdots & \displaystyle\int_0^1 g_N(x)g_N”(x)\,dx \end{bmatrix} \begin{Bmatrix} a_1\\a_2\\\vdots\\a_N \end{Bmatrix} = \begin{Bmatrix} \displaystyle-\int_0^1 g_1(x)f(x)\,dx\\ \displaystyle-\int_0^1 g_2(x)f(x)\,dx\\ \vdots\\ \displaystyle-\int_0^1 g_N(x)f(x)\,dx \end{Bmatrix} \tag{32} \end{align}

　こうして左辺の係数行列と右辺の定数ベクトルはすべて既知関数の積分として計算されるので、これより未知量である \(a_i\ (i=1,2,\ldots,N)\) が求められます。

例題2

　先の例題1をガラーキン法で解けば、

\begin{align} -a_1\pi^2 \int_0^1\sin^2(\pi x)\,dx = -\int_0^1\sin(\pi x)\,dx \tag{33} \end{align}

　より、\(a_1\) が得られ、近似解は、

\begin{align} \tilde{u}(x) = \frac{4}{\pi^3}\sin(\pi x) \tag{34} \end{align}

となります。この近似解を図4に示します。この問題の場合には、例題1における近似解よりもガラーキン法による近似解の方が精度がよいです。

なぜ試験関数を基底関数に一致させるのか

　試験関数を基底関数と一致させるという操作は、単に未知数と方程式の数をそろえるための便宜ではありません。その意味を理解するには、基底関数と試験関数がそれぞれ何を表しているかを分けて考える必要があります。

　まず、式\((22)\)の基底関数 \(g_i(x)\) は、近似解を変化させることのできる方向を表します。たとえば係数 \(a_i\) を \(\delta a_i\) だけ変えれば、近似解の変化は、

\begin{align} \delta\tilde{u}(x) = \delta a_i g_i(x) \tag{35} \end{align}

となります。したがって、\(g_1,g_2,\ldots,g_N\) は、有限個の自由度で近似解を修正できる \(N\) 個の方向と考えることができます。

　一方、試験関数 \(v_i(x)\) は解を表すための関数ではなく、近似解を支配方程式へ代入したときに生じる誤差を調べるための重みです。式\((19)\)に対する残差を、

\begin{align} R(x) = \frac{d^2\tilde{u}(x)}{dx^2} +f(x) \tag{36} \end{align}

　と定義すると、重み付き残差法は、

\begin{align} \int_0^1 v_i(x)R(x)\,dx=0 \qquad (i=1,2,\ldots,N) \tag{37} \end{align}

　を要求します。つまり、試験関数は残差のどの成分を検査するかを決める物差しです。

　ガラーキン法では、この試験関数として基底関数そのものを選びます。すなわち \(v_i(x)=g_i(x)\) とするので、式\((37)\)は、

\begin{align} \int_0^1 g_i(x)R(x)\,dx=0 \qquad (i=1,2,\ldots,N) \tag{38} \end{align}

となります。この式は、残差 \(R(x)\) がすべての基底関数に対して直交することを表します。言い換えると、近似解を修正できるどの方向にも、修正に使える残差成分を残さないように未知係数を決めています。

　たとえば二つの基底関数を用いて、

\begin{align} \tilde{u}(x) = a_1g_1(x)+a_2g_2(x) \tag{39} \end{align}

　と近似した場合、近似解を変更する方法は、\(a_1\) を変えて \(g_1\) の方向へ修正する方法と、\(a_2\) を変えて \(g_2\) の方向へ修正する方法の二つです。そこでガラーキン法は、

\begin{align} \begin{aligned} \int_0^1 g_1(x)R(x)\,dx &=0,\\ \int_0^1 g_2(x)R(x)\,dx &=0 \end{aligned} \tag{40} \end{align}

という二つの条件から \(a_1\) と \(a_2\) を決定します。1行目は \(g_1\) の方向に修正すべき残差成分がないことを、2行目は \(g_2\) の方向に修正すべき残差成分がないことを意味します。

　ここで、残差 \(R(x)\) そのものが領域内のすべての点で0になるわけではないことに注意が必要です。有限個の基底関数では表現できない残差成分は残る可能性があります。それでも、現在採用している基底関数で表現できる範囲では、どの方向へ係数を動かしても改善に使える残差成分が残っていない状態を選んでいます。

　したがって、試験関数を基底関数と一致させる真意は、近似解を動かせる方向と残差を検査する方向を一致させることにあります。これにより、有限個の基底関数で作られた近似空間の中で、支配方程式の残差を偏りなく検査した近似解が得られます。

MATLABコード

　今回用いたコードを以下に記載します。

MATLABコード

% make_blog_figures.m
%
% 1次元弾性棒の支配方程式と、重み付き残差法・ガラーキン法による
% 近似解を確認するためのブログ掲載用スクリプトです。
%
% このファイルを実行すると、次の図を同じフォルダへ保存します。
%   図2: 1次元弾性棒と微小区間の力のつり合い
%   図3: 変位差と軸ひずみの関係
%   図4: 厳密解、近似解および応力分布の比較
%
% スクリプト末尾の補助関数も含めて、この1ファイルだけで実行できます。

clear; close all; clc;
set(groot,'defaultFigureVisible','off');

fontName = 'Arial';
fontSize = 18;
lineWidth = 2.0;
blue = [0.12 0.36 0.72];
orange = [0.90 0.30 0.08];
green = [0.10 0.52 0.30];
purple = [0.48 0.24 0.62];
dark = [0.15 0.17 0.20];
lightBlue = [0.82 0.89 0.98];

%% 図2: 1次元弾性棒と微小区間の力のつり合い
fig = figure('Color','w','Position',[100 100 1400 560], ...
    'ToolBar','none','MenuBar','none');
tiledlayout(fig,1,2,'TileSpacing','compact','Padding','compact');

ax = nexttile;
hold(ax,'on'); axis(ax,[0 1 0 1]); axis(ax,'off');

% 両端の固定壁と棒を描きます。
rectangle(ax,'Position',[0.12 0.30 0.035 0.38], ...
    'FaceColor',dark,'EdgeColor','none');
rectangle(ax,'Position',[0.845 0.30 0.035 0.38], ...
    'FaceColor',dark,'EdgeColor','none');
rectangle(ax,'Position',[0.155 0.43 0.69 0.13], ...
    'FaceColor',lightBlue,'EdgeColor',dark,'LineWidth',1.8);

% 棒の軸方向に作用する物体力 f(x) を矢印で示します。
for xq = linspace(0.23,0.72,5)
    filledArrow(ax,xq,0.64,xq+0.075,0.64,orange,lineWidth,0.018,0.014);
end
text(ax,0.50,0.76,'f(x)','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize+2,'Interpreter','tex');
text(ax,0.20,0.495,'A','HorizontalAlignment','center', ...
    'VerticalAlignment','middle','FontName',fontName, ...
    'FontSize',fontSize+1,'FontWeight','bold');

% 座標軸、棒長さ、材料定数および境界条件を表示します。
filledArrow(ax,0.155,0.20,0.905,0.20,dark,1.5,0.018,0.012);
plot(ax,[0.155 0.155],[0.18 0.23],'Color',dark,'LineWidth',1.3);
plot(ax,[0.845 0.845],[0.18 0.23],'Color',dark,'LineWidth',1.3);
text(ax,0.155,0.12,'0','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize);
text(ax,0.845,0.12,'l','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'FontAngle','italic');
text(ax,0.93,0.20,'x','VerticalAlignment','middle', ...
    'FontName',fontName,'FontSize',fontSize,'FontAngle','italic');
title(ax,'(a) fixed elastic bar','FontName',fontName, ...
    'FontSize',fontSize,'FontWeight','normal');

ax = nexttile;
hold(ax,'on'); axis(ax,[0 1 0 1]); axis(ax,'off');

% 微小区間を取り出し、左右の応力と物体力を描きます。
rectangle(ax,'Position',[0.39 0.34 0.22 0.32], ...
    'FaceColor',[0.96 0.96 0.96],'EdgeColor',dark,'LineWidth',1.8);
filledArrow(ax,0.39,0.50,0.21,0.50,blue,lineWidth,0.026,0.018);
filledArrow(ax,0.61,0.50,0.79,0.50,blue,lineWidth,0.026,0.018);
filledArrow(ax,0.45,0.50,0.55,0.50,orange,lineWidth,0.022,0.016);
text(ax,0.50,0.57,'f','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize+1,'FontAngle','italic');
text(ax,0.23,0.40,'\sigma(x)','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize-1,'Interpreter','tex');
text(ax,0.77,0.40,'\sigma(x+dx)','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize-1,'Interpreter','tex');
text(ax,0.39,0.23,'x','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'FontAngle','italic');
text(ax,0.61,0.23,'x+dx','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'Interpreter','tex');
title(ax,'(b) equilibrium of a differential segment', ...
    'FontName',fontName,'FontSize',fontSize,'FontWeight','normal');

exportgraphics(fig,'fig02_fem_1d_bar_equilibrium.png','Resolution',400);
exportgraphics(fig,'fig02_fem_1d_bar_equilibrium.pdf','ContentType','vector');
close(fig);

%% 図3: 変位差と軸ひずみの関係
fig = figure('Color','w','Position',[100 100 1150 650], ...
    'ToolBar','none','MenuBar','none');
ax = axes(fig,'Position',[0.05 0.08 0.90 0.86]);
hold(ax,'on'); axis(ax,[0 1 0 1]); axis(ax,'off');

xL = 0.22; xR = 0.52;
xLd = 0.30; xRd = 0.67;
yBefore = 0.35; yAfter = 0.68; barH = 0.09;

rectangle(ax,'Position',[xL yBefore xR-xL barH], ...
    'FaceColor',[0.91 0.92 0.94],'EdgeColor',dark,'LineWidth',1.7);
rectangle(ax,'Position',[xLd yAfter xRd-xLd barH], ...
    'FaceColor',lightBlue,'EdgeColor',dark,'LineWidth',1.7);

% 変形前の座標軸と端点位置を読み取るための補助線を描きます。
filledArrow(ax,0.14,0.18,0.84,0.18,dark,1.5,0.022,0.014);
for xv = [xL xR]
    plot(ax,[xv xv],[0.18 yBefore],'--','Color',[0.45 0.45 0.45], ...
        'LineWidth',1.4);
end
for xv = [xLd xRd]
    plot(ax,[xv xv],[yAfter-0.04 yAfter+barH+0.02],':', ...
        'Color',blue,'LineWidth',1.4);
end

% 左端と右端の変位を矢印で示します。
filledArrow(ax,xL,0.60,xLd,0.60,orange,lineWidth,0.022,0.016);
filledArrow(ax,xR,0.53,xRd,0.53,orange,lineWidth,0.022,0.016);
text(ax,(xL+xLd)/2,0.625,'u(x)','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'Interpreter','tex');
text(ax,(xR+xRd)/2,0.555,'u(x+dx)','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'Interpreter','tex');

% 変形前の微小長さ dx を両矢印で示します。
doubleFilledArrow(ax,xL,0.10,xR,0.10,green,1.7,0.022,0.014);
text(ax,(xL+xR)/2,0.045,'dx','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'Interpreter','tex');

text(ax,xL,0.14,'x','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'FontAngle','italic');
text(ax,xR,0.14,'x+dx','HorizontalAlignment','center', ...
    'FontName',fontName,'FontSize',fontSize,'Interpreter','tex');
text(ax,0.72,yBefore+barH/2,'before deformation', ...
    'VerticalAlignment','middle','FontName',fontName,'FontSize',fontSize);
text(ax,0.72,yAfter+barH/2,'after deformation', ...
    'VerticalAlignment','middle','FontName',fontName,'FontSize',fontSize);
text(ax,0.88,0.18,'x','VerticalAlignment','middle', ...
    'FontName',fontName,'FontSize',fontSize,'FontAngle','italic');
title(ax,'displacement difference and axial strain', ...
    'FontName',fontName,'FontSize',fontSize+2,'FontWeight','normal');

exportgraphics(fig,'fig03_fem_displacement_strain_kinematics.png','Resolution',400);
exportgraphics(fig,'fig03_fem_displacement_strain_kinematics.pdf','ContentType','vector');
close(fig);

%% 図4: 厳密解と近似解の比較
% 次の境界値問題を考えます。
%
%   -d^2u/dx^2 = 1,   0 < x < 1
%   u(0) = u(1) = 0
%
% 厳密解は u(x)=x(1-x)/2 です。近似解には、両端で0となる
% 基底関数 g(x)=sin(pi*x) を使い、u_tilde(x)=a*g(x) と置きます。
x = linspace(0,1,1201);
uExact = 0.5*x.*(1-x);

g = sin(pi*x);

% 試験関数 v(x)=1 を用いる重み付き残差法では、
% integral_0^1 R(x) dx = 0 から a=1/(2*pi) を得ます。
aWeightedResidual = 1/(2*pi);
uWeightedResidual = aWeightedResidual*g;

% ガラーキン法では試験関数を基底関数と一致させ、v(x)=g(x) とします。
% integral_0^1 g(x)R(x) dx = 0 から a=4/pi^3 を得ます。
aGalerkin = 4/pi^3;
uGalerkin = aGalerkin*g;

% 比較用として、節点を x=0.5 に置いた区分線形近似も計算します。
gHalf = 2*x;
gHalf(x > 0.5) = 2*(1-x(x > 0.5));
uPiecewiseHalf = (1/8)*gHalf;

% 同様に、節点を x=1/3 に置いた区分線形近似を計算します。
alpha = 1/3;
gAlpha = x/alpha;
gAlpha(x > alpha) = (1-x(x > alpha))/(1-alpha);
uPiecewiseThird = alpha*(1-alpha)/2*gAlpha;

% 近似係数とガラーキン近似の最大誤差をコマンドウィンドウへ表示します。
fprintf('Weighted-residual coefficient = %.10f\n', aWeightedResidual);
fprintf('Galerkin coefficient          = %.10f\n', aGalerkin);
fprintf('Maximum Galerkin error        = %.3e\n', ...
    max(abs(uGalerkin-uExact)));

fig = figure('Color','w','Position',[100 100 1500 680], ...
    'ToolBar','none','MenuBar','none');
t = tiledlayout(fig,1,2,'TileSpacing','compact','Padding','compact');
ax = nexttile(t);
plot(ax,x,uExact,'-','Color',dark,'LineWidth',3.0); hold(ax,'on');
plot(ax,x,uWeightedResidual,'--','Color',orange,'LineWidth',2.4);
plot(ax,x,uGalerkin,'-.','Color',blue,'LineWidth',2.5);
plot(ax,x,uPiecewiseHalf,':','Color',green,'LineWidth',2.7);
plot(ax,x,uPiecewiseThird,'--','Color',purple,'LineWidth',2.2);
grid on; box on;
xlim([0 1]); ylim([0 0.20]);
xticks(0:0.2:1); yticks(0:0.025:0.20);
xlabel('x','FontName',fontName,'FontSize',fontSize);
ylabel('displacement  u(x)','FontName',fontName,'FontSize',fontSize);
title('exact and approximate displacement solutions', ...
    'FontName',fontName,'FontSize',fontSize+2,'FontWeight','normal');
lgd = legend(ax,{'Eq. (16): exact','Eq. (31): weighted residual', ...
    'Eq. (34): Galerkin','piecewise linear: node at x = 0.5', ...
    'piecewise linear: node at x = 1/3'}, ...
    'NumColumns',2, ...
    'FontName',fontName,'FontSize',fontSize-4);
lgd.Layout.Tile = 'south';
set(ax,'FontName',fontName,'FontSize',fontSize, ...
    'LineWidth',1.1,'Layer','top');

% 厳密解を微分すると、ひずみと応力は x に対して直線になります。
% ここでは分布形状を見るため、応力を f0*l で正規化して表示します。
ax = nexttile(t);
sigmaNormalized = 0.5-x;
plot(ax,x,sigmaNormalized,'-','Color',blue,'LineWidth',3.0); hold(ax,'on');
yline(ax,0,'-','Color',dark,'LineWidth',1.2);
xline(ax,0.5,':','Color',[0.35 0.35 0.35],'LineWidth',1.8);
plot(ax,0.5,0,'o','MarkerSize',8,'MarkerFaceColor',orange, ...
    'MarkerEdgeColor',dark,'LineWidth',1.2);
grid(ax,'on'); box(ax,'on');
xlim(ax,[0 1]); ylim(ax,[-0.58 0.58]);
xticks(ax,0:0.2:1);
xlabel(ax,'x/l','FontName',fontName,'FontSize',fontSize);
ylabel(ax,'normalized stress  \sigma/(f_0 l)', ...
    'FontName',fontName,'FontSize',fontSize,'Interpreter','tex');
title(ax,'stress and strain implied by the exact solution', ...
    'FontName',fontName,'FontSize',fontSize+2,'FontWeight','normal');
text(ax,0.28,0.48,'tension: du/dx > 0', ...
    'HorizontalAlignment','center','FontName',fontName, ...
    'FontSize',fontSize-1,'Color',blue);
text(ax,0.72,-0.48,'compression: du/dx < 0', ...
    'HorizontalAlignment','center','FontName',fontName, ...
    'FontSize',fontSize-1,'Color',blue);
text(ax,0.53,0.08,'du/dx = 0', ...
    'HorizontalAlignment','left','FontName',fontName, ...
    'FontSize',fontSize-1,'Color',dark);
set(ax,'FontName',fontName,'FontSize',fontSize, ...
    'LineWidth',1.1,'Layer','top');

exportgraphics(fig,'fig04_fem_exact_approximate_solutions.png','Resolution',400);
exportgraphics(fig,'fig04_fem_exact_approximate_solutions.pdf','ContentType','vector');
close(fig);

fprintf('\nSaved blog figures:\n');
fprintf('  fig02_fem_1d_bar_equilibrium.png\n');
fprintf('  fig03_fem_displacement_strain_kinematics.png\n');
fprintf('  fig04_fem_exact_approximate_solutions.png\n');

function filledArrow(ax,x0,y0,x1,y1,color,lineWidth,headLength,headWidth)
%FILLEDARROW 塗りつぶした三角形の矢じりを持つ片矢印を描きます。
dx = x1-x0;
dy = y1-y0;
arrowLength = hypot(dx,dy);
ux = dx/arrowLength;
uy = dy/arrowLength;
px = -uy;
py = ux;
baseX = x1-headLength*ux;
baseY = y1-headLength*uy;

plot(ax,[x0 baseX],[y0 baseY],'-','Color',color,'LineWidth',lineWidth);
patch(ax, ...
    [x1,baseX+headWidth*px/2,baseX-headWidth*px/2], ...
    [y1,baseY+headWidth*py/2,baseY-headWidth*py/2], ...
    color,'EdgeColor',color,'LineWidth',0.8);
end

function doubleFilledArrow(ax,x0,y0,x1,y1,color,lineWidth,headLength,headWidth)
%DOUBLEFILLEDARROW 両端の矢じりが同じ大きさの寸法線を描きます。
dx = x1-x0;
dy = y1-y0;
arrowLength = hypot(dx,dy);
ux = dx/arrowLength;
uy = dy/arrowLength;
px = -uy;
py = ux;

leftBaseX = x0+headLength*ux;
leftBaseY = y0+headLength*uy;
rightBaseX = x1-headLength*ux;
rightBaseY = y1-headLength*uy;
plot(ax,[leftBaseX rightBaseX],[leftBaseY rightBaseY], ...
    '-','Color',color,'LineWidth',lineWidth);
patch(ax, ...
    [x0,leftBaseX+headWidth*px/2,leftBaseX-headWidth*px/2], ...
    [y0,leftBaseY+headWidth*py/2,leftBaseY-headWidth*py/2], ...
    color,'EdgeColor',color,'LineWidth',0.8);
patch(ax, ...
    [x1,rightBaseX+headWidth*px/2,rightBaseX-headWidth*px/2], ...
    [y1,rightBaseY+headWidth*py/2,rightBaseY-headWidth*py/2], ...
    color,'EdgeColor',color,'LineWidth',0.8);
end

おわりに

　この記事では、一次元弾性棒を題材として、釣合い条件、適合条件、構成式から変位を未知関数とする微分方程式を導きました。また、その微分方程式を一自由度ばねの \(F=kx\) と対応づけることで、外力に対して変位を求めるという見慣れた考え方と、連続体の支配方程式との関係を確認しました。

　さらに、厳密解を求める代わりに、有限個の基底関数を組み合わせて近似解を表す考え方へ進みました。近似解を支配方程式へ代入すると残差が生じますが、重み付き残差法では、その残差を領域全体で評価して未知係数を決めます。なかでもガラーキン法は、近似解を表す基底関数と残差を検査する試験関数を一致させる方法です。これにより、近似解を修正できる各方向について、修正に使える残差成分が残らないように係数を決定できます。

　ここまでの議論により、微分方程式として表された連続体の問題を、有限個の未知係数に対する連立方程式へ置き換えるまでの基本的な流れが見えてきました。次の記事では、この考え方を土台として、有限要素法の具体的な定式化へ進みます。