有限要素法の基礎②：形状関数、要素剛性行列、全体行列の組立て【MATLABコード付き】

はじめに

　前の記事では、一次元弾性棒の支配方程式を導き、近似解、重み付き残差法、ガラーキン法までを説明しました。しかし、その段階では近似解を領域全体にわたる基底関数で表しており、まだ「領域を小さな要素に分けて計算する」という有限要素法らしい操作は現れていません。

　この記事では、ガラーキン法を要素ごとの計算へ分解します。まず、支配方程式を部分積分して弱形式を導きます。次に、節点に結び付いた折れ線状の基底関数を導入し、各要素の要素剛性行列と定数ベクトルを求めます。最後に、それらを全体行列へ組み立てる操作と、節点における力の釣合いとの関係を確認します。

境界値問題

　前の記事と同じく、\(0\le x\le 1\) の領域で、次の境界値問題を考えます。

\begin{align} -\frac{d^2u(x)}{dx^2}=f(x) \qquad (0<x<1) \tag{1} \end{align}

\begin{align} u(0)=0, \qquad u(1)=0 \tag{2} \end{align}

　式\((1)\)は領域内部のすべての点で成立すべき微分方程式であり、強形式と呼ばれます。\(u(x)\) は未知関数、\(f(x)\) は既知の分布荷重に相当する関数です。式\((2)\)は、領域の両端で変位を0に固定する境界条件です。

　有限要素法では、式\((1)\)を各点で直接満たす近似解を探すのではなく、まず領域全体の積分で表された式へ変換します。その理由は、要素ごとに折れ線状の近似解を使うと、その2階微分を要素境界で通常の関数として扱いにくいからです。

強形式から弱形式を導く

試験関数を掛けて領域全体で評価する

　式\((1)\)が領域全体としてどの程度満たされているかを調べるため、任意の試験関数 \(v(x)\) を掛けて積分します。

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

　ここで試験関数 \(v(x)\) は、残差のどの成分を調べるかを決める重みです。両端では変位 \(u\) が指定されているため、近似解を両端で変化させないように、試験関数にも

\begin{align} v(0)=0, \qquad v(1)=0 \tag{4} \end{align}

　を課します。

部分積分によって微分の階数を下げる

　式\((3)\)の第1項には \(u\) の2階微分が含まれています。そこで、部分積分

\begin{align} \int_0^1 v(x)\frac{d^2u}{dx^2}dx = \left[ v(x)\frac{du}{dx} \right]_0^1 – \int_0^1 \frac{dv}{dx}\frac{du}{dx}dx \tag{5} \end{align}

　を用います。式\((5)\)を式\((3)\)へ代入すると、

\begin{align} \int_0^1 \frac{dv}{dx}\frac{du}{dx}dx – \left[ v\frac{du}{dx} \right]_0^1 = \int_0^1 v f\,dx \tag{6} \end{align}

となります。さらに、式\((4)\)より境界項は0になるため、最終的に

\begin{align} \int_0^1 \frac{dv}{dx}\frac{du}{dx}dx = \int_0^1 v f\,dx \tag{7} \end{align}

を得ます。これが今回用いる弱形式です。

　強形式の式\((1)\)には \(u\) の2階微分が必要でしたが、弱形式の式\((7)\)に必要なのは \(u\) と \(v\) の1階微分だけです。弱形式は、近似解に必要な滑らかさを1段階下げます。

　ここまでで、支配方程式を弱形式へ変換できました。しかし、式\((7)\)の未知関数 \(u(x)\) は、まだ領域内で連続的に変化する関数のままです。このままでは、領域内の無数の位置における値を求めなければならず、有限個の未知数を使う数値計算にはなっていません。

　そこで有限要素法では、領域を小さな区間に分け、区切り位置に節点を置きます。そして、未知関数そのものを直接求める代わりに、節点における値だけを未知数として求めます。隣り合う節点値の間を最も単純な直線で結べば、要素内では直線、領域全体では折れ線状の近似関数が得られます。

　ところが、この折れ線は節点で曲がっています。強形式が要求する2階微分を節点で計算できないため、「このような関数を近似解として使ってよいのか」という疑問が生じます。ここで、先ほど導いた弱形式が必要になります。次節では、まず有限要素法が折れ線近似を導入する理由を整理し、その折れ線を弱形式なら扱えることを確認します。

弱形式から有限要素近似へ

まず節点間を直線で表す

　有限要素法では、連続な未知関数 \(u(x)\) を有限個の節点値で表します。たとえば領域の途中に節点を置き、隣り合う節点値を直線で結べば、各要素内の関数は2個の端点値だけで決まります。直線の傾きも要素内で一定になるため、変位からひずみに相当する1階微分を簡単に計算できます。

　複数の要素を並べると、隣接要素は同じ節点値を共有します。そのため、近似関数そのものは節点で途切れず、領域全体で連続になります。一方、要素ごとの直線の傾きは一般に異なるため、節点では折れ曲がります。ここで使う折れ線は、厳密解ではなく近似解です。

　次の「有限要素法で使う特殊な基底関数」の節では、この直線補間を形状関数によって具体的に組み立てます。その準備として、まず単純な折れ線を微分し、強形式では扱えない一方で、弱形式では積分できることを確認します。

折れ線は1階微分できても、節点では2階微分できない

　中央で折れ曲がる、次の連続な関数を考えます。

\begin{align} \tilde{u}(x) = \begin{cases} x, & 0\leq x\leq 1/2,\\ 1-x, & 1/2<x\leq 1 \end{cases} \tag{8} \end{align}

　式\((8)\)は \(x=1/2\) でつながっているため、関数そのものは連続です。これを区間ごとに1回微分すると、

\begin{align} \frac{d\tilde{u}}{dx} = \begin{cases} 1, & 0<x<1/2,\\ -1, & 1/2<x<1 \end{cases} \tag{9} \end{align}

となります。各区間では傾きを計算できますが、\(x=1/2\) では傾きが \(1\) から \(-1\) へ跳ぶため、1階微分は不連続です。さらに区間内でもう1回微分すれば、

\begin{align} \frac{d^2\tilde{u}}{dx^2}=0 \qquad \left(x\neq \frac{1}{2}\right) \tag{10} \end{align}

　ですが、折れ曲がっている \(x=1/2\) では2階微分を通常の値として定義できません。したがって、式\((8)\)を強形式の式\((1)\)へ直接代入しても、節点では左辺を計算できません。また、要素内部では式\((10)\)の左辺が0になるため、\(f(x)\neq 0\) の支配方程式を各点で満たすこともできません。

弱形式では要素ごとの傾きを積分できる

　一方、弱形式の式\((7)\)に現れるのは、\(u\) の2階微分ではなく1階微分です。式\((9)\)を弱形式の左辺へ入れ、折れ曲がる位置で積分区間を分けると、

\begin{align} \begin{aligned} \int_0^1 \frac{dv}{dx}\frac{d\tilde{u}}{dx}\,dx &= \int_0^{1/2}\frac{dv}{dx}\,dx – \int_{1/2}^{1}\frac{dv}{dx}\,dx \\ &= \left\{v(1/2)-v(0)\right\} – \left\{v(1)-v(1/2)\right\} \end{aligned} \tag{11} \end{align}

となります。試験関数の境界条件 \(v(0)=v(1)=0\) を使えば、式\((11)\)は、

\begin{align} \int_0^1 \frac{dv}{dx}\frac{d\tilde{u}}{dx}\,dx =2v(1/2) \tag{12} \end{align}

となります。この計算には \(d^2\tilde{u}/dx^2\) が一度も現れません。必要なのは、左右の区間で定義できる1階微分だけです。したがって、傾きが節点で跳んでいても、要素ごとに積分して足し合わせることができます。

　ここで重要なのは、節点の折れ曲がりを無視しているわけではないことです。式\((12)\)には折れ曲がる位置の試験関数値 \(v(1/2)\) が残っており、左右の傾きの違いが積分結果へ反映されています。有限要素法で要素ごとの式を全体方程式へ組み立てるときにも、要素境界の寄与は共有節点へ集められます。

「滑らかさを1段階下げる」の意味

　強形式では、領域内の各点で \(d^2u/dx^2\) を計算できるだけの滑らかさが必要です。これに対し、弱形式では \(u\) が連続であり、要素ごとに \(du/dx\) を計算できれば積分を評価できます。つまり、要求される微分の階数が2階から1階へ下がります。これが、滑らかさを1段階下げるという意味です。

　ただし、折れ線が強形式を各点で厳密に満たすようになるわけではありません。弱形式では、支配方程式を領域全体の積分として評価します。そのため、折れ線を有限要素近似の候補として使い、試験関数で選んだ各方向について残差が残らないように節点値を決められます。十分に滑らかな厳密解に対しては強形式と弱形式は同じ物理を表しており、弱形式は支配方程式を別の意味へ変更するものではありません。

有限要素法で使う特殊な基底関数

領域を要素に分割する

　領域 \(0\le x\le1\) を、図1のように三つの要素 \(e_1,e_2,e_3\) に分割します。要素境界に置かれた四つの節点座標を \(x_1,x_2,x_3,x_4\) とし、各節点における近似解の値を \(a_1,a_2,a_3,a_4\) とします。

図1 三要素に分割した一次元領域と、節点2、3に対応する全体基底関数です。下段は、節点値 \(a_2=3\)、\(a_3=2\) を与えたときの有限要素近似を示します。

　式\((2)\)より、両端の節点値は

\begin{align} a_1=0, \qquad a_4=0 \tag{13} \end{align}

です。したがって、未知量は内部節点の \(a_2\) と \(a_3\) だけです。

節点に結び付いた帽子型関数

　節点2に対応する基底関数 \(g_2(x)\) は、節点2で1、隣接節点1、3で0となり、それ以外の要素では0となる折れ線関数です。同様に、\(g_3(x)\) は節点3で1となります。一般に、節点 \(i\) に対応する基底関数は、

\begin{align} g_i(x_j)= \begin{cases} 1 & (i=j),\\ 0 & (i\ne j) \end{cases} \tag{14} \end{align}

という性質を持ちます。この性質のため、近似解

\begin{align} \tilde{u}(x) = a_2g_2(x)+a_3g_3(x) \tag{15} \end{align}

　を節点 \(x_2\)、\(x_3\) で評価すると、それぞれ

\begin{align} \tilde{u}(x_2)=a_2, \qquad \tilde{u}(x_3)=a_3 \tag{16} \end{align}

となります。つまり、未知係数 \(a_i\) は単なる展開係数ではなく、節点における近似解の値そのものです。

ガラーキン法を連立方程式にする

　ここで、弱形式に導入した試験関数 \(v(x)\) と、近似解を作る基底関数 \(g_i(x)\) を結び付けます。ガラーキン法では、節点 \(i\) に対応する方程式を作るとき、試験関数を同じ節点の基底関数に一致させ、\(v_i(x)=g_i(x)\) と置きます。本例の未知節点は2、3なので、式\((7)\)へ \(v_2=g_2\)、\(v_3=g_3\) を順に代入します。さらに、近似解には式\((15)\)を代入すると、

\begin{align} \begin{bmatrix} \displaystyle\int_0^1 g_2’g_2’\,dx & \displaystyle\int_0^1 g_2’g_3’\,dx\\[2mm] \displaystyle\int_0^1 g_3’g_2’\,dx & \displaystyle\int_0^1 g_3’g_3’\,dx \end{bmatrix} \begin{Bmatrix} a_2\\a_3 \end{Bmatrix} = \begin{Bmatrix} \displaystyle\int_0^1 g_2 f\,dx\\[2mm] \displaystyle\int_0^1 g_3 f\,dx \end{Bmatrix} \tag{17} \end{align}

を得ます。ここでは \(g_i’=dg_i/dx\) と書きました。式\((17)\)の第1行は \(v_2=g_2\) を選んだときの方程式、第2行は \(v_3=g_3\) を選んだときの方程式です。したがって、試験関数 \(v\) が途中で消えたのではありません。試験関数として選んだ \(g_2,g_3\) が、行列の各行を作る関数として残っています。

　式\((17)\)を

\begin{align} \mathbf{K}\mathbf{a}=\mathbf{b} \tag{18} \end{align}

と表します。\(\mathbf{K}\) は全体剛性行列、\(\mathbf{a}\) は未知節点値ベクトル、\(\mathbf{b}\) は分布荷重から作られる全体定数ベクトルです。

　式\((17)\)を領域全体で一度に積分しても解は求められます。しかし有限要素法では、基底関数が各節点の周囲にしか値を持たない性質を利用し、積分を要素ごとに分けます。これによって、同じ計算手順をすべての要素へ繰り返し適用できるようになります。

一要素の計算を完成させる

局所形状関数

　図2のように、左節点 \(p\)、右節点 \(q\)、長さ \(l_e=x_q-x_p\) の要素 \(e\) を取り出します。前節の \(g_i(x)\) は領域全体で定義された全体基底関数です。この全体基底関数を一つの要素 \(e\) の内部だけに制限したものを、局所形状関数 \(N_p(x)\)、\(N_q(x)\) と呼びます。

図2 一要素内の形状関数とその微分です。線形要素では形状関数の微分が要素内で一定になります。

　形状関数は、

\begin{align} \begin{aligned} N_p(x)&=\frac{x_q-x}{l_e}, \qquad N_q(x)=\frac{x-x_p}{l_e},\\ \left.g_p\right|_e&=N_p, \qquad \left.g_q\right|_e=N_q \end{aligned} \tag{19} \end{align}

です。記号 \(\left.g_p\right|_e\) は「全体基底関数 \(g_p\) を要素 \(e\) の区間だけで見る」という意味です。つまり、\(g_p\) と \(N_p\) は別の関数ではなく、同じ基底関数を全体から見るか、一要素の中だけで見るかという違いです。この要素上では、節点 \(p,q\) 以外に対応する全体基底関数は0です。

　式\((15)\)の近似解を要素 \(e\) の区間 \(x_p\le x\le x_q\) に限定し、式\((19)\)の対応 \(g_p|_e=N_p\)、\(g_q|_e=N_q\) を使うと、要素内の近似解は、

\begin{align} \tilde{u}^{\,e}(x) = N_p(x)a_p+N_q(x)a_q = \underbrace{ \begin{bmatrix} N_p(x) & N_q(x) \end{bmatrix} }_{\mathbf{N}^{e}(x)} \underbrace{ \begin{Bmatrix} a_p\\a_q \end{Bmatrix} }_{\mathbf{a}^{e}} \tag{20} \end{align}

　と表されます。式\((20)\)は、要素内部の値を両端の節点値から線形補間する式です。

近似解の微分と \(B\) 行列

　式\((20)\)を \(x\) で微分すると、

\begin{align} \frac{d\tilde{u}^{\,e}}{dx} = \underbrace{ \begin{bmatrix} \displaystyle\frac{dN_p}{dx} & \displaystyle\frac{dN_q}{dx} \end{bmatrix} }_{\mathbf{B}^{e}} \mathbf{a}^{e} \tag{21} \end{align}

となります。ここで \(\mathbf{B}^{e}\) は、節点値から要素内の勾配を求める行列です。式\((19)\)を微分すれば、

\begin{align} \mathbf{B}^{e} = \frac{1}{l_e} \begin{bmatrix} -1 & 1 \end{bmatrix} \tag{22} \end{align}

です。線形要素では \(\mathbf{B}^{e}\) が要素内で一定になるため、要素内の勾配も一定になります。

要素剛性行列

　式\((7)\)の左辺を要素 \(e\) の区間 \(x_p\le x\le x_q\) に限定します。近似解の微分には式\((21)\)を代入します。

　ここで、これまでの関数の置き換えを確認します。ガラーキン法では全体領域で \(v_i=g_i\) と置きました。さらに要素内では、式\((19)\)のように \(g_p|_e=N_p\)、\(g_q|_e=N_q\) と表します。したがって、要素内の試験関数は局所形状関数の線形結合として \(v^e=c_pN_p+c_qN_q=\mathbf{N}^e\mathbf{c}^e\) と書け、その微分は \(dv^e/dx=\mathbf{B}^e\mathbf{c}^e\) となります。ここで \(\mathbf{c}^e\) は、要素内で試験する形状関数を選ぶ任意係数ベクトルです。これらを式\((7)\)の左辺へ代入すると、

\begin{align} \int_{x_p}^{x_q} \frac{dv^e}{dx} \frac{d\tilde{u}^{\,e}}{dx}\,dx = \left(\mathbf{c}^{e}\right)^{\mathrm{T}} \left[ \int_{x_p}^{x_q} \left(\mathbf{B}^{e}\right)^{\mathrm{T}}\mathbf{B}^{e}\,dx \right] \mathbf{a}^{e} \tag{23} \end{align}

となります。この式では、試験関数 \(v\) は \(\mathbf{c}^e\) と \(\mathbf{B}^e\) の形で残っています。そこで、式\((23)\)の角括弧内を要素剛性行列

\begin{align} \mathbf{K}^{e} = \int_{x_p}^{x_q} \left(\mathbf{B}^{e}\right)^{\mathrm{T}}\mathbf{B}^{e}\,dx \tag{24} \end{align}

と定義します。式\((22)\)を代入し、積分過程を省略せずに計算すると、

\begin{align} \begin{aligned} \mathbf{K}^{e} &= \int_{x_p}^{x_q} \frac{1}{l_e^2} \begin{bmatrix} -1\\1 \end{bmatrix} \begin{bmatrix} -1&1 \end{bmatrix} dx\\ &= \frac{1}{l_e^2} \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix} \int_{x_p}^{x_q}dx\\ &= \frac{1}{l_e} \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix} \end{aligned} \tag{25} \end{align}

となります。要素が短いほど \(1/l_e\) が大きくなり、要素剛性は大きくなります。これは、同じ材料と断面を持つ短い棒ほど変形しにくいという、ばね剛性 \(EA/l_e\) の性質と対応しています。今回の無次元化された問題では \(EA=1\) に相当するため、式\((25)\)の形になります。

要素定数ベクトル

　弱形式の右辺を要素 \(e\) に限定すると、要素定数ベクトルは、

\begin{align} \mathbf{b}^{e} = \int_{x_p}^{x_q} \left(\mathbf{N}^{e}\right)^{\mathrm{T}}f(x)\,dx = \begin{Bmatrix} \displaystyle\int_{x_p}^{x_q}N_p(x)f(x)\,dx\\[2mm] \displaystyle\int_{x_p}^{x_q}N_q(x)f(x)\,dx \end{Bmatrix} \tag{26} \end{align}

です。これは、要素内に連続的に分布する \(f(x)\) を、両端節点へ等価な節点量として振り分ける操作です。\(f(x)=f_0\) が要素内で一定なら、

\begin{align} \mathbf{b}^{e} = \frac{f_0l_e}{2} \begin{Bmatrix} 1\\1 \end{Bmatrix} \tag{27} \end{align}

　となり、要素全体の荷重 \(f_0l_e\) が左右の節点へ半分ずつ配分されます。

　以上により、一要素について解くべき式は、

\begin{align} \mathbf{K}^{e}\mathbf{a}^{e}=\mathbf{b}^{e} \tag{28} \end{align}

　の形に整理されました。ただし、共有節点を持つ要素は互いに独立ではないため、式\((28)\)を要素ごとに別々に解くのではなく、次節で全体方程式へ組み立てます。

要素行列を全体行列へ組み立てる

局所節点番号を全体節点番号へ対応させる

　三要素の等分割では、各要素の節点対応は、

\begin{align} e_1:(1,2), \qquad e_2:(2,3), \qquad e_3:(3,4) \tag{29} \end{align}

です。要素 \(e_2\) の左上成分は全体行列の \((2,2)\) 成分へ、右下成分は \((3,3)\) 成分へ入ります。同じ全体成分へ複数要素の寄与が入る場合には、それらを加算します。

図3 三つの要素剛性行列を全体節点番号に対応する位置へ配置し、重なった成分を加算する組立て操作です。

　各要素長が \(l_e=1/3\) の場合、式\((25)\)より、

\begin{align} \mathbf{K}^{e_1} = \mathbf{K}^{e_2} = \mathbf{K}^{e_3} = 3 \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix} \tag{30} \end{align}

です。これらを式\((29)\)に従って組み立てると、

\begin{align} \mathbf{K} = 3 \begin{bmatrix} 1&-1&0&0\\ -1&2&-1&0\\ 0&-1&2&-1\\ 0&0&-1&1 \end{bmatrix} \tag{31} \end{align}

となります。対角成分には、その節点へ接続する要素の寄与が足し合わされています。たとえば節点2には要素 \(e_1\) と \(e_2\) が接続するため、\((2,2)\) 成分は \(3+3=6\) です。

　同様に、要素定数ベクトルも全体節点番号に従って加算します。\(f(x)=1\) の場合、各要素について

\begin{align} \mathbf{b}^{e} = \frac{1}{6} \begin{Bmatrix} 1\\1 \end{Bmatrix} \tag{32} \end{align}

　なので、全体定数ベクトルは、

\begin{align} \mathbf{b} = \frac{1}{6} \begin{Bmatrix} 1\\2\\2\\1 \end{Bmatrix} \tag{33} \end{align}

となります。

境界条件を反映する

　式\((13)\)の \(a_1=a_4=0\) を用いると、未知量は \(a_2,a_3\) だけになります。式\((31)\)と式\((33)\)から、内部節点に対する方程式は、

\begin{align} \begin{bmatrix} 6&-3\\ -3&6 \end{bmatrix} \begin{Bmatrix} a_2\\a_3 \end{Bmatrix} = \begin{Bmatrix} 1/3\\1/3 \end{Bmatrix} \tag{34} \end{align}

です。ここまで来れば、有限要素法の問題は通常の連立一次方程式を解く問題になっています。

組立ては節点力の釣合いでもある

　ここで、全体行列への組立ては、共有節点における節点力の釣合いを表します。全体行列の組立ては、行列成分を機械的に足す操作として理解できます。同時に、力学的には共有節点における力の釣合いを作る操作でもあります。

　要素 \(e\) の節点力ベクトルを、

\begin{align} \mathbf{F}^{e} = \mathbf{K}^{e}\mathbf{a}^{e}-\mathbf{b}^{e} \tag{35} \end{align}

と定義します。右辺第1項は節点変位によって要素内部に生じる抵抗、第2項は要素内の分布荷重に対応する等価節点量です。

図4 二つの隣接要素が共有する節点 \(q\) における節点力の釣合いです。要素ごとの節点力を全体節点 \(q\) へ集めると、その和が0になります。

　図4の節点 \(q\) では、左要素 \(e\) と右要素 \(f\) の寄与を足し合わせて、

\begin{align} F_q^{e}+F_q^{f}=0 \tag{36} \end{align}

　を満たします。これが全体方程式の節点 \(q\) に対応する1行です。したがって、要素行列の同じ全体行・列への加算は、その節点に接続するすべての要素の力を集め、釣合い式を作る操作と解釈できます。

数値例による確認

例1：一様荷重を受ける三等分要素

　まず、\(f(x)=1\) とし、領域を \(x=0,1/3,2/3,1\) で三等分します。解くべき方程式は式\((34)\)です。これを解くと、

\begin{align} \begin{Bmatrix} a_2\\a_3 \end{Bmatrix} = \begin{Bmatrix} 1/9\\1/9 \end{Bmatrix} \tag{37} \end{align}

となります。一方、式\((1)\)の厳密解は、

\begin{align} u(x)=\frac{1}{2}x(1-x) \tag{38} \end{align}

です。式\((38)\)を \(x=1/3,2/3\) で評価すると、どちらも \(1/9\) になります。したがって、この問題では有限要素解の節点値が厳密解と一致します。

図5 一様荷重 \(f(x)=1\) に対する厳密解と有限要素解の比較です。左は等分割、右は不等分割です。有限要素解は要素内では直線ですが、いずれも節点では厳密解と一致します。

　ただし、図5から分かるように、要素内部では有限要素解が直線、厳密解が放物線であるため、両者は一致しません。「節点値が厳密解と一致する」ことと、「領域内のすべての点で厳密解と一致する」ことは別です。

例2：要素長を不均一にする

　次に、節点を \(x=0,1/2,3/4,1\) に置きます。各要素長は \(1/2,1/4,1/4\) なので、式\((25)\)より、

\begin{align} \mathbf{K}^{e_1} = 2 \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix}, \qquad \mathbf{K}^{e_2} = \mathbf{K}^{e_3} = 4 \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix} \tag{39} \end{align}

となります。定数ベクトルは、

\begin{align} \mathbf{b}^{e_1} = \frac{1}{4} \begin{Bmatrix} 1\\1 \end{Bmatrix}, \qquad \mathbf{b}^{e_2} = \mathbf{b}^{e_3} = \frac{1}{8} \begin{Bmatrix} 1\\1 \end{Bmatrix} \tag{40} \end{align}

です。組立てと境界条件の反映を行うと、

\begin{align} \begin{bmatrix} 6&-4\\ -4&8 \end{bmatrix} \begin{Bmatrix} a_2\\a_3 \end{Bmatrix} = \begin{Bmatrix} 3/8\\1/4 \end{Bmatrix} \tag{41} \end{align}

　となり、

\begin{align} \begin{Bmatrix} a_2\\a_3 \end{Bmatrix} = \begin{Bmatrix} 1/8\\3/32 \end{Bmatrix} \tag{42} \end{align}

を得ます。この値も、式\((38)\)を \(x=1/2,3/4\) で評価した値と一致します。要素長が異なっても、各要素について同じ式\((25)\)と式\((26)\)を使い、\(l_e\) と積分区間だけを変えれば計算できます。これは有限要素法をプログラム化しやすい重要な理由です。

例3：分布荷重が途中で不連続になる場合

　最後に、分布荷重を、

\begin{align} f(x)= \begin{cases} 1 & \left(0\le x<1/2\right),\\ 2 & \left(1/2\le x\le1\right) \end{cases} \tag{43} \end{align}

　とします。要素分割は例1と同じ \(x=0,1/3,2/3,1\) です。要素剛性行列は変わりませんが、要素定数ベクトルは式\((26)\)の積分から、

\begin{align} \mathbf{b}^{e_1} = \frac{1}{6} \begin{Bmatrix} 1\\1 \end{Bmatrix}, \qquad \mathbf{b}^{e_2} = \frac{1}{24} \begin{Bmatrix} 5\\7 \end{Bmatrix}, \qquad \mathbf{b}^{e_3} = \frac{1}{3} \begin{Bmatrix} 1\\1 \end{Bmatrix} \tag{44} \end{align}

となります。中央の要素 \(e_2\) の内部で \(f(x)\) が変化するため、左右の節点へ等分されない点が重要です。組立て後の方程式は、

\begin{align} \begin{bmatrix} 6&-3\\ -3&6 \end{bmatrix} \begin{Bmatrix} a_2\\a_3 \end{Bmatrix} = \begin{Bmatrix} 3/8\\5/8 \end{Bmatrix} \tag{45} \end{align}

　であり、その解は、

\begin{align} \begin{Bmatrix} a_2\\a_3 \end{Bmatrix} = \begin{Bmatrix} 11/72\\13/72 \end{Bmatrix} \tag{46} \end{align}

です。厳密解は、

\begin{align} u(x)= \begin{cases} \displaystyle-\frac{1}{2}x^2+\frac{5}{8}x & \left(0\le x\le1/2\right),\\[2mm] \displaystyle-x^2+\frac{9}{8}x-\frac{1}{8} & \left(1/2\le x\le1\right) \end{cases} \tag{47} \end{align}

となります。

図6 領域中央で不連続になる分布荷重と、その厳密解および有限要素解です。荷重が大きい右側では曲率が大きくなり、変位の最大位置も中央より右へ移動します。

　図6では、荷重が右半分で大きいため、厳密解の曲率 \(-u”=f\) も右半分で大きくなります。その結果、変位分布は左右非対称となり、最大変位の位置が中央より右へ移動します。一方、有限要素解は各要素内で直線ですが、節点では今回も厳密解と一致しています。

支配方程式から計算手順までを整理する

　ここまでに導いた有限要素法の流れを、図7にまとめます。試験関数 \(v\) はガラーキン法によって全体基底関数 \(g_i\) と一致させ、全体基底関数は各要素内で局所形状関数 \(N_i\) として扱います。その後、要素行列を全体行列へ組み立て、境界条件を反映して節点値を求めます。図の上半分は支配方程式を要素行列へ落とし込む定式化、下半分は得られた行列から節点値と要素内の近似解を求める計算手順です。

図7 支配方程式から弱形式、要素行列、全体方程式、解の評価へ進む一次元有限要素法の流れです。

おわりに

　この記事では、前の記事で学んだガラーキン法を、要素ごとの計算へ分解しました。弱形式を導くことで、2階微分を含む支配方程式を、1階微分だけを必要とする積分式へ変換しました。そのうえで、節点に結び付いた形状関数を使い、一要素の要素剛性行列と定数ベクトルを求めました。

　有限要素法の中心的な操作は、同じ形式で計算した要素行列を、接続関係に従って全体行列へ加算することです。この組立て操作は、行列の作成手順であると同時に、共有節点における力の釣合いを作る操作でもあります。最終的には、連続体の微分方程式が、節点値を未知量とする連立一次方程式へ置き換えられました。

　ここまでで、一次元有限要素法の最小構成が一通りそろいました。次の記事では、今回の一次元線形要素で得た考え方を土台として、より一般的な要素や数値積分へ進みます。