1 連立方程式の設定

1.1 ラプラス方程式について

ラプラス方程式は物理学の分野ではよく出て来る式である。 ３変数$x,y,z$の関数 $\phi (x,y,z)$の場合

$$\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2} \right)\phi (x,y,z) = 0$$ (1)

をいう。これを $\Delta\phi = 0$と書くことも多い。 この方程式に従う関数の例は、静電場、重力ポテンシャル、定常状態の温度分布 等がある。

今は、２次元の場合のラプラス方程式を考える。２次元のラプラス方程式は以下のようになる。

$$\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)\phi (x,y) = 0$$ (2)

有限要素法ではこの微分方程式を変分形の方程式に直して計算していく。

1.2 領域の要素分割

次に示すラプラス方程式の変分形を、有限要素法で解く方法を導出する。

$$\frac{1}{2}\int\!\!\!\int_D\left\{\left(\frac{\partial u}{\partial x}\right)^2 +\left(\frac{\partial u}{\partial y}\right)^2\right\}dxdy$$ (3)

領域$D$の中で、ラプラス方程式を解くときに、まず領域を三角形要素に分割する。 要素には番号がふられているものとして、$i$番目の要素の領域は、 （便宜上）$(I)$と表すものとする。 このときの、三角形の頂点にあたる節点の値を求めようとする。 この節点には番号がふられているものとし、節点$i$の$u$の値を$u_i$とする。 節点での$u$の値 $u_1,u_2,\cdots ,u_n$をすべて計算することにより、 領域全体にわたり$u$の値を近似していく。

1.3 決定方程式

未知関数$u(x,y)$を近似するために、次のような一次関数の形状関数を導入する。 要素$(I)$の頂点が、節点$i,j,k$だったとするとき、形状関数 $\phi_{(I)i}$は、

$$\phi_{(I)i}= \begin{cases} 0 &(領域(I)以外)　\\ 節点iで1、節点j,kで0となる１次関数 &(領域(I)の中) \end{cases}$$ (4)

のようになり、具体的な $\phi_{(I)i}$を表すと、以下のようになる。

$$\phi_{(I)i} = \frac{(y_j - y_k)x + (x_k-x_j)y + (x_j y_k - x_k y_j)}{s_{I}}$$ (5)

ここで、$(x_i,y_i)$は節点$i$の座標、$s_{(I)}$は領域$(I)$の面積である。

一つの領域$(I)$の中の求める方程式は、次の形で近似される。

$$u = u_i\phi_{(I)i}+u_j\phi_{(I)j}+u_k\phi_{(I)k}$$ (6)

解くべき領域全体では、この方程式の足し合わせで近似される。 領域$D$全体では以下のようになる。

$$\begin{split} u &= \sum^N_{I=1}\left\{ u_{(I1)}\phi_{(I1)}+... ...\ &= \sum^N_{I=1}\sum^3_{t=1}u_{(It)}\phi_{(It)} \end{split}$$ (7)

この式の添え字の説明をする。 $u_{(I1)}$の意味は、要素番号$I$の三角形領域が持っている一つ目の節点での、$u$の値。 同様に、$u_{(I2)}$は、要素番号$I$の三角形領域が持っている二つ目の節点での、$u$の値。 また、 $\phi_{(I1)}$の意味は、要素番号$I$の三角形領域が持っている一つ目の節点が 節点$i$だったとすると、 $\phi_{(I1)}$は $\phi_{(I)i}$と等価である。

このような形状関数を使い、$J[u]$の変分を計算していく。

1.4 変分の計算

汎関数$J[u]$を決定関数を使って表すと、次のようになる。

$$\begin{split} J[u] &= \frac{1}{2}\int\!\!\!\int_D\left\{\lef... ...tial \phi_{(It)}}{\partial y}\right)^2\right\}dxdy \end{split}$$ (8)

汎関数$J[u]$が、関数 $J(u_1,u_2,\cdots ,u_n)$になったので、 $u_1,u_2,\cdots ,u_n$で偏微分していく。 ただし、$u_i$が固定点（ディレクレ条件）だった場合は除く。

$u_i$で偏微分した場合を計算していく。 この場合、変数として$u_i$を含む積分領域だけ抜き出して計算すればよい。 すなわち、以下のようになる。

$$\begin{split} \frac{\partial J[u]}{\partial u_i} &= \frac{... ...}}{\partial y}\right)^2\right\}dxdy \\ &+ \cdots \end{split}$$ (9)

ここでの、領域 $(I),(J),(K),\cdots\$は頂点に節点$i$をもつ三角形領域である。 つまり、節点$(It)$のどれか（節点 $(I1),(I2),(I3)$のどれか）は節点$i$である （節点$(Jt)$も同様）。

さらに計算していくと以下のようになる。

$$\begin{split} \frac{\partial J[u]}{\partial u_i} &= \frac{... ...}}{\partial y}\right)^2\right\}dxdy \\ &+ \cdots \end{split}$$ (10)

ここで、 $\frac{\partial}{\partial u_i} \left(\sum^3_{t=1}u_{(It)}\frac{\partial \phi_{(It)}}{\partial x}\right)^2$ を計算すると、以下のようになる。

$$\frac{\partial}{\partial u_i} \left(\sum^3_{t=1}u_{(It)} \frac{\... ...rtial \phi_{(It)}}{\partial x} \right) \frac{\partial \phi_{(I1)}}{\partial x}$$ (11)

ここで、 $u_i=u_{(I1)}$であるとしている。

これをもとに、

$$\begin{split} \frac{\partial J[u]}{\partial u_i} &= \int\!... ...i_{(K)i}}{\partial y} \right\}dxdy \\ &+ \cdots \end{split}$$ (12)

となる。ここで、 $(Jt),(Kt),\cdots\$のそれぞれにおいて、 ３つのうちどれかは節点$i$を示すものである。 さらに、積分の中には、変数として$x,y$が入っていないので、以下のようになる。

$$\begin{split} \frac{\partial J[u]}{\partial u_i} &= \left\... ... \right\} \int\!\!\!\int_{(K)} dxdy \\ &+ \cdots \end{split}$$ (13)

ここで、積分 $\int\!\!\!\int_{(I)} dxdy$はその領域の面積である。 これを$s_{(I)}$とすると、

$$\begin{split} \frac{\partial J[u]}{\partial u_i} &= \left\... ...(K)i}}{\partial y} \right\} s_{(K)} \\ &+ \cdots \end{split}$$ (14)

となる。

1.5 連立方程式

連立方程式は、 $\frac{\partial J[u]}{\partial u_i}=0$としたものを連立させる。 また、ディレクレ条件の場合はその点での値がわかっているので、$u_i=C_i$とする。

$$\begin{cases} \frac{\partial J[u]}{\partial u_1} = 0 \\ \f... ... u_{m+1} = C_{m+1} \\ \ \ \vdots \\ u_n = C_n \end{cases}$$ (15)

ここで、ディレクレ条件ではない節点を、節点$1〜m$として、 ディレクレ条件の節点を、節点$m+1〜n$としている。 これを行列 ${\bf Au = b}$の形で表すと、

$$\begin{pmatrix} a_{11} &a_{12} &\cdots &a_{1m} &a_{1m+1} &\cdots... ...egin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ C_{m+1} \\ \vdots \\ C_n \end{pmatrix}$$ (16)

これより、 $a_{11}〜a_{mn}$までを求めなくてはいけないことが分かる。

ホームページ: Yamamoto's laboratory
著者: 夏井拓也