|TOP Page|

FEM for Poisson's Equation

ポアソン方程式の有限要素法による数値解法







Last Update:2009年05月16日


概要

ポアソン方程式の弱形式化

ポアソン方程式(Poisson's Equation)を解く問題を次のように定義する。

問題設定
\{ \quad \nabla^2 \phi = f \qquad ( \quad in \quad V \quad ) \\ \quad \phi = g_1 \qquad (\quad on \quad S_1 \quad) \\ \quad \frac{\partial \phi}{\partial n } = g_2 \qquad (\quad on \quad S_2 \quad)

境界条件は値が決定されている場合のDirichlet条件と，法線方向の微分値が決定されている場合のNeumann条件の２種類に分けられる．これらは，第一種境界条件，第二種境界条件または，基本境界条件，自然境界条件と呼ばれる場合もある．境界はDirichlet条件かNeumann条件のどちらかが設定されていなくては，解が定まらない．また，Dirichlet条件とNeumann条件が与えられる境界は重複してはならない．つまり，

S_1 \cup S_2 = S \qquad S_1 \cap S_2 = \emptyset

以下強形式で書かれたポアソン方程式を弱形式に変形する.\delta\phiを

\delta\phi = 0 \qquad ( \quad on \quad S_1 \quad )

を満たす任意の関数とする.上式に\delta\phiをかけて，V上で積分すると

\int_V \nabla^2 \phi \delta \phi dV=\int_V f \delta\phi dV

となる. ここで，微分の連鎖則より

\frac{\partial}{\partial x_k}(\frac{\partial\phi}{\partial x_k}\delta\phi) = \frac{\partial}{\partial x_k}(\frac{\partial\phi}{\partial x_k})\delta\phi + \frac{\partial\phi}{\partial x_k}\frac{\partial\delta\phi}{\partial x_k}

が成り立つ．これをV内で積分すると

\int_V \nabla \cdot ( \nabla \phi \delta \phi ) dV = \int_V \nabla^2 \phi \delta \phi dV + \int_V (\nabla \phi) \cdot (\nabla \delta \phi) dV

第２項に前式を代入する．また，右辺と左辺を入れ替えると，

\int_V (\nabla \phi) \cdot (\nabla \delta \phi) dV = \int_V \nabla \cdot ( \nabla \phi \delta \phi ) dV + \int_V f \delta\phi dV

右辺第一項について考える。この項にGaussの発散定理を用いると

\int_V \nabla \cdot ( \nabla \phi \delta \phi ) dV = \int_S n \cdot \nabla \phi \delta \phi dS = \int_S \frac{ \partial \phi }{ \partial n } \delta \phi dS = \int_{S_1} \frac{ \partial \phi }{ \partial n } \delta \phi dS + \int_{S_2} g_2 \delta \phi dS\\ \qquad \qquad = \int_{S_2} g_2 \delta \phi dS

となる。ここでS_1上では\delta \phiは０であることを用いた

これを上式に代入すると弱形式で書かれたポアソン方程式は以下の通りになる。

弱形式ポアソン方程式
\{\int_V (\nabla \phi) \cdot (\nabla \delta \phi) dV = \int_V f \delta\phi dV + \int_{S_2} g_2 \delta \phi dS \qquad \forall\delta\phi\ni \{ \delta\phi \mid \delta\phi=0 \quad (on \quad S_1) \} \\ \phi = g_1 \qquad (\quad on \quad S_1 \quad)

弱形式ラプラス方程式の有限要素法離散化

弱形式で書かれたラプラス方程式は以下の通り

\{ \int_V (\nabla \phi) \cdot (\nabla \delta \phi) dV = \int_V f \delta \phi dV + \int_{S_2} g_2 \delta \phi dS \qquad \forall\delta\phi\ni \{ \delta\phi \mid \delta\phi=0 \quad (on \quad S_1) \} \\ \phi = g_1 \qquad (\quad on \quad S_1 \quad)

上式を成分で書くと

\int_V\frac{\partial\delta\phi}{\partial x_k} \frac{\partial\phi}{\partial x_k}dV = \int_Vf\delta\phi dV + \int_{S_2} g_2 \delta \phi dS

のようになる。ここではkについて総和規約を適応している。上式は積分で書かれている。積分とはつまり和であるので積分領域を分割してた中で積分を実行し、それを足し合わせて上式を計算してもよい。

そこでVを要素分割した領域をV_e、S_2を要素分割した領域を{S_2}_eとすると定式は次のように書ける。

\sum_e^{nV_e}\int_V_e\frac{\partial\delta\phi}{\partial x_k} \frac{\partial\phi}{\partial x_k}dV = \sum_e^{nV_e}\int_V_ef\delta\phi dV + \sum_e^{nS_2_e}\int_{{S_2}_e} g_2 \delta \phi dS

但し、nV_e、nS_2_eはそれぞれVとS_2_eの要素分割数である。

V_e内で関数\phi,\delta\phiは次のような形状関数Nによって離散化されているとする。
但し、固定境界条件が設定されている面S_1上の節点\bar{n}では\phi^{\bar{n}}=g_1、\delta\phi^{\bar{n}}=0である。

\phi = N^n\phi^n、\delta\phi = N^m\delta\phi^m

同様にS_2_e内で関数\phi,\delta\phiは次のような形状関数\bar{N}によって離散化されているとする。

\phi = \bar{N}^n\phi^n、\delta\phi = \bar{N}^m\delta\phi^m

形状関数で補間された\phi,\delta\phiを上式に代入して

\sum_{e}^{nV_e}\int_V_e\frac{\partial N^m}{\partial x_k}\frac{\partial N^n}{\partial x_k}\phi^n \delta\phi^m dV_e = \sum_{e}^{nV_e}\int_{V_e} f N^m \delta\phi^m dV + \sum_{e}^{nS_2_e}\int_{S_2_e} g_2 \bar{N}^m \delta\phi^m dS

ここで\phi^n,\delta\phi^mは節点の値であり，関数ではないので積分の外に出すことができる。それぞれ積分の外に出した後にまとめると次のようになる。

\delta\phi^m \( \sum_{e}^{nV_e}\int_V_e\frac{\partial N^m}{\partial x_k}\frac{\partial N^n}{\partial x_k}dV \) \phi^n = \delta\phi^m \( \sum_{e}^{nV_e}\int_{V_e} f N^m dV + \sum_{e}^{nS_2_e}\int_{S_2} g_2 \bar{N}^m dS \)

次のように行列表記できる。

離散化されたポアソン方程式
\delta\phi^mA^{m,n}\phi^n=\delta\phi^mb^m

A^{m,n} = \sum_{e}^{nV_e}A_e^{a,b}
(但し、要素内節点番号a,bは全体節点番号m,n)
要素剛性行列
A_e^{a,b}=\int_V_e\frac{\partial N^a}{\partial x_k}\frac{\partial N^b}{\partial x_k}dV
b^{m} = \sum_{e}^{nV_e}\int_{V_e} f N^m dV + \sum_{e}^{nS_2_e}\int_{S_2} g_2 \bar{N}^m dS

上式をマトリックス表記すると次のようになる。
\delta\b{\phi}^T\b{A}\b{\phi}=\delta\b{\phi}^T\b{b}

境界条件処理

このままでは連立一次方程式として解くことができない。
固定境界条件が設定されている面S_1上の節点\bar{n}では\phi^{\bar{n}}=g_1、\delta\phi^{\bar{n}}=0はGivenであった。面S_1上の\phiの値を\bar{\phi}とする。また、未知数のある面S_1上以外の点における\phi,\delta\phiの値を\hat{\phi},\delta\hat{\phi}とすると上の方程式は次のようになる。

\(\begin{array}\delta\hat{\phi},0\end{array}\)\b{A}\(\begin{array}\hat{\phi}\\ \bar{\phi} \end{array}\) = \(\begin{array}\delta\hat{\phi},0\end{array}\)\b{b}

これは展開して次のように書くことができる。

\delta\hat{\phi}\hat{\b{A}}\hat{\phi}=\delta\hat{\phi}\(\hat{\b{b}}-\tilde{\b{A}}\bar{\phi}\)

但し、\hat{\b{A}}は未知自由度について\b{A}の成分を取り出してきた行列である。また\tilde{\b{A}}は行が未知自由度、列が既知自由度である\b{A}の成分を取り出してきた行列である、さらに\hat{\b{b}は未知自由度について\b{b}の成分を取り出してきたベクトルである。

ここで任意の\hat{\phi}の値について上式が成り立たなければいけなかったことに注意して

\hat{\b{A}}\hat{\phi}=\hat{\b{b}}-\tilde{\b{A}}\bar{\phi}

のように書ける。この形にすることによって、計算して解を求めることが可能である。


残差と増分について解くとより簡単な式が得られる。上式を増分について解いてみる。\hat{\phi}=\hat{\phi}_0+\Delta\hat{\phi}を代入すると

\hat{\b{A}}(\hat{\phi}_0+\Delta\hat{\phi})=\hat{\b{b}}-\tilde{\b{A}}\bar{\phi}

これを変形すると、

\hat{\b{A}}\Delta\hat{\phi}=\hat{\b{b}}-\tilde{\b{A}}\bar{\phi} - \hat{\b{A}}\hat{\phi}_0=\hat{\b{r}}_0

ただし、\hat{\b{r}}_0は残差\b{r}_0=\b{b}-\b{A}\b{\phi}_0の未知自由度のみを取り出したものである。\b{A}\b{\phi}_0を計算することは、\tilde{A}\bar{\phi}を計算するよりも計算負荷が大きくなる可能性がある。しかし、非線形方程式を解く場合は通常残差と増分に対して解くので残差はどのみち計算しなければならず余分な手間とはならない。非線形方程式を解く場合はこの種類の境界条件処理をよく使う。

要素での離散化例

単純な形状の場合は数値積分を用いなくても解析的に積分を実行することができる．この場合について，簡単に纏めておく．

シンプレックス(三角形,４面体)一次要素

三角形や四面体の一次の内挿関数は線形である。節点aに対する内挿関数をL^aのように書く。

線形であるので空間微分は定数。つまり要素内で\frac{\partial L^a}{\partial x^i}は一定である。これから剛性行列を求める際の被積分も定数となるので、この定数に要素の面積や体積を掛けることで積分を実行することができる。これ用いると次のように要素剛性行列を書くことができる。

A_e^{a,b}=\int_V_e\frac{\partial N^a}{\partial x_k}\frac{\partial N^b}{\partial x_k}dV=\Delta\frac{\partial L^a}{\partial x_k}\frac{\partial L^b}{\partial x_k}

但し、\Deltaは２次元の場合は三角形の面積、３次元の場合は四面体の体積である。

ここで、

area=\Delta、dldx[a][i] = \frac{\partial L^a}{\partial x_i}

emat[a][b] = A_e^{a,b}

のように変数が対応している場合、プログラムの例は次のとおり

// 要素剛性行列を求める
for(unsigned int inoel=0;inoel<nnoel;inoel++){
for(unsigned int jnoel=0;jnoel<nnoel;jnoel++){
      double dtmp1 = 0.0;
      for(unsigned int idim=0;idim<ndim;idim++){
         dtmp1 += dldx[inoel][idim]*dldx[jnoel][idim];
      }
      emat[inoel][jnoel] = area*dtmp1;
}
}

直交四角形要素

\{\begin{array}{l}N^1=(1-\frac{x}{h_x})(1-\frac{y}{h_y})\\ N^2=\frac{x}{h_x}(1-\frac{y}{h_y})\\ N^3=\frac{x}{h_x}\frac{y}{h_y}\\ N^4=(1-\frac{x}{h_x})\frac{y}{h_y}\end{array}

これらのｘ微分、ｙ微分を計算すると、

\{\begin{array}{ll}\frac{\partial N^1}{\partial x}=-\frac{1}{h_x}(1-\frac{y}{h_y}) & \frac{\partial N^1}{\partial y}=-(1-\frac{x}{h_x})\frac{1}{h_y}\\ \frac{\partial N^2}{\partial x}=\frac{1}{h_x}(1-\frac{y}{h_y}) & \frac{\partial N^2}{\partial y}=-\frac{x}{h_x}\frac{1}{h_y}\\ \frac{\partial N^3}{\partial x}=\frac{1}{h_x}\frac{y}{h_y} & \frac{\partial N^3}{\partial y}=\frac{x}{h_x}\frac{1}{h_y}\\ \frac{\partial N^4}{\partial x}=-\frac{1}{h_x}\frac{y}{h_y} & \frac{\partial N^4}{\partial y}=(1-\frac{x}{h_x})\frac{1}{h_y}\end{array}

積分を行うと次のとおり、

\int_{V_e}\frac{\partial N^1}{\partial x}\frac{\partial N^1}{\partial x}dV = \int_{V_e}\frac{\partial N^2}{\partial x}\frac{\partial N^2}{\partial x}dV = \int_{V_e}\frac{\partial N^3}{\partial x}\frac{\partial N^3}{\partial x}dV = \int_{V_e}\frac{\partial N^4}{\partial x}\frac{\partial N^4}{\partial x}dV = \frac{1}{3}\frac{h_y}{h_x}\\ \int_{V_e}\frac{\partial N^1}{\partial x}\frac{\partial N^2}{\partial x}dV = \int_{V_e}\frac{\partial N^3}{\partial x}\frac{\partial N^4}{\partial x}dV = -\frac{1}{3}\frac{h_y}{h_x}\\ \int_{V_e}\frac{\partial N^1}{\partial x}\frac{\partial N^3}{\partial x}dV = \int_{V_e}\frac{\partial N^2}{\partial x}\frac{\partial N^4}{\partial x}dV = -\frac{1}{6}\frac{h_y}{h_x}\\ \int_{V_e}\frac{\partial N^1}{\partial x}\frac{\partial N^4}{\partial x}dV = \int_{V_e}\frac{\partial N^2}{\partial x}\frac{\partial N^3}{\partial x}dV = \frac{1}{6}\frac{h_y}{h_x}


\int_{V_e}\frac{\partial N^1}{\partial y}\frac{\partial N^1}{\partial y}dV = \int_{V_e}\frac{\partial N^2}{\partial y}\frac{\partial N^2}{\partial y}dV = \int_{V_e}\frac{\partial N^3}{\partial y}\frac{\partial N^3}{\partial y}dV = \int_{V_e}\frac{\partial N^4}{\partial y}\frac{\partial N^4}{\partial y}dV = \frac{1}{3}\frac{h_x}{h_y}\\ \int_{V_e}\frac{\partial N^1}{\partial y}\frac{\partial N^2}{\partial y}dV = \int_{V_e}\frac{\partial N^3}{\partial y}\frac{\partial N^4}{\partial y}dV = \frac{1}{6}\frac{h_x}{h_y}\\ \int_{V_e}\frac{\partial N^1}{\partial y}\frac{\partial N^3}{\partial y}dV = \int_{V_e}\frac{\partial N^2}{\partial y}\frac{\partial N^4}{\partial y}dV = -\frac{1}{6}\frac{h_x}{h_y}\\ \int_{V_e}\frac{\partial N^1}{\partial y}\frac{\partial N^4}{\partial y}dV = \int_{V_e}\frac{\partial N^2}{\partial y}\frac{\partial N^3}{\partial y}dV = -\frac{1}{3}\frac{h_x}{h_y}

よって、要素剛性行列は次のとおり

A^{a,b}=\int_V_e\frac{\partial N^a}{\partial x}\frac{\partial N^b}{\partial x}+\frac{\partial N^a}{\partial y}\frac{\partial N^b}{\partial y}dV =\frac{1}{6h_yh_x}\[\begin{array}{llll} 2(h_y^2+h_x^2) & -2h_y^2+h_x^2& -h_y^2-h_x^2 & h_y^2-2h_x^2\\ -2h_y^2+h_x^2 & 2(h_y^2+h_x^2) & h_y^2-2h_x^2 & -h_y^2-h_x^2\\ -h_y^2-h_x^2 & h_y^2-2h_x^2 & 2(h_y^2+h_x^2) & -2h_y^2+h_x^2\\ h_y^2-2h_x^2 & -h_y^2-h_x^2 & -2h_y^2+h_x^2 & 2(h_y^2+h_x^2)\end{array}\]

直交正方形要素

ここで、h_x=h_yの時は特別に次のようになる。

A^{a,b}=\[\begin{array}{cccc}\frac{2}{3}&-\frac{1}{6}&-\frac{1}{3}&-\frac{1}{6}\\ -\frac{1}{6}&\frac{2}{3}&-\frac{1}{6}&-\frac{1}{3}\\ -\frac{1}{3}&-\frac{1}{6}&\frac{2}{3}&-\frac{1}{6}\\ -\frac{1}{6}&-\frac{1}{3}&-\frac{1}{6}&\frac{2}{3}\\ \end{array}\]

各点におけるステンシルは次のとおり

\frac{1}{3}\[\begin{array}{ccc}-1&-1&-1\\-1&8&-1\\-1&-1&-1\end{array}\]

サンプルプログラム（直交四角形要素）

次のサンプルプログラムは大きさ１×１の矩形領域における次のようなポアソン方程式(Poisson's Equation)を解く問題である。

-{\nabla}^2\phi=1 \qquad ( \quad in \quad \Omega = \{(x,y):0.0<x,y< 1.0 \} \quad )

\phi=0 \qquad ( \quad on \quad \partial\Omega \quad )

結果の表示

プログラムを実行すると'solution.dat'という結果ファイルが生成される。

'solution.dat'はgnuplot用のデータ形式をしている

windows版のgnuplotの導入は以下のページが詳しい

初歩 gnuplot入門
http://auemath.aichi-edu.ac.jp/~khotta/ghost/gnuplot.html
物理の鍵しっぽ(Windows版のGnuplotの初期設定)
http://hooktail.org/computer/index.php?Windows%C8%C7Gnuplot%A4%CE%BD%E9%B4%FC%C0%DF%C4%EA



gnuplotを起動した後、プログラムを実行した'solution.dat'のあるディレクトリに移動し、次のようなコマンドを打つと結果が表示される。

gnuplot> set style data line
gnuplot> splot 'solution.dat'

軸対称問題

軸対称な問題を考えてみよう．つまり解析領域や境界条件や物性に周方向の変化が無い場合を考える．この場合解も軸対称と考えられ，厚さ方向の変化を無視した場合と同様に３次元の問題を２次元的に取り扱うことができる．

円柱座標におけるラプラス方程式

円柱座標系を考えた場合，対象となる場は軸までの距離とｚ座標の変数として表される．ここでは簡単のため，回転軸がｚ軸と一致していると考える．この時，場は次のように表される．

\phi(x,y,z) = \phi(r,z)\;\;\; r=\sqrt{x^2+y^2}

これらの微分を計算すると以下のとおりになる．

\{\begin{array}{l} \frac{\partial r}{\partial x} = \frac{x}{\sqrt{x^2+y^2}}=\frac{x}{r}\;\;\;\;\; & \frac{\partial^2 r}{\partial x^2} = \frac{\partial}{\partial x}(\frac{x}{r}) = \frac{1}{r} - \frac{x^2}{r^3}\\ \frac{\partial r}{\partial y} = \frac{y}{\sqrt{x^2+y^2}}=\frac{y}{r}\;\;\;\;\; & \frac{\partial^2 r}{\partial y^2} = \frac{\partial}{\partial y}(\frac{y}{r}) = \frac{1}{r} - \frac{y^2}{r^3}\\ \frac{\partial^2 r}{\partial x\partial y} = -\frac{xy}{r^3}\end{array}

これらを用いて

\frac{\partial^2 \phi}{\partial x^2} = \frac{\partial}{\partial x}\(\frac{\partial\phi}{\partial x}\) = \frac{\partial}{\partial x}\{\(\frac{\partial\phi}{\partial r}\)\(\frac{\partial r}{\partial x}\)\} = \(\frac{\partial^2\phi}{\partial r^2}\)\(\frac{\partial r}{\partial x}\)^2 + \(\frac{\partial\phi}{\partial r}\)\(\frac{\partial^2 r}{\partial x^2}\)\\\;\;\;\; =\(\frac{\partial^2\phi}{\partial r^2}\)\(\frac{x}{r}\)^2 + \(\frac{\partial\phi}{\partial r}\)\(\frac{1}{r} - \frac{x^2}{r^3}\)

\frac{\partial^2 \phi}{\partial y^2} = \(\frac{\partial^2\phi}{\partial r^2}\)\(\frac{y}{r}\)^2 + \(\frac{\partial\phi}{\partial r}\)\(\frac{1}{r} - \frac{y^2}{r^3}\)

これらを用いて直交座標におけるラプラス演算子を計算すると以下のとおりになる．

\(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}\)\phi(r,z) = \(\frac{\partial^2\phi}{\partial r^2}\)\(\frac{x^2+y^2}{r^2}\) + \(\frac{\partial\phi}{\partial r}\)\(\frac{2}{r} - \frac{x^2+y^2}{r^3}\) + \frac{\partial^2 \phi}{\partial z^2}\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \frac{\partial^2\phi}{\partial r^2} + \frac{1}{r}\(\frac{\partial\phi}{\partial r}\) + \frac{\partial^2 \phi}{\partial z^2}

よって，軸対称問題の円筒座標系におけるラプラス方程式は次のように書ける

軸対称問題の円筒座標系におけるラプラス方程式
\frac{\partial^2\phi}{\partial r^2} + \frac{1}{r}\(\frac{\partial\phi}{\partial r}\) + \frac{\partial^2 \phi}{\partial z^2} = f(r,z)

fは単位体積あたりのソースを表している．このfも軸対称であることに注意されたい．

弱形式化

軸対称に固定境界条件とし自然境界条件が定義されているとする．

固定境界上で０となる，軸対称な関数\delta\phi(r,z)を上式に掛け，領域全体で積分する．

\int_V\delta\phi (\nabla^2\phi+f) dxdydz = 0\;\; \forall \delta\phi

\{\begin{array}{l}x\rightarrow rsin\theta\\ y\rightarrow rcos\theta\\ z\rightarrow z\end{array}

この場合のヤコビアンの絶対値は

\det J(r,\theta,z)=\det \[\begin{array} \sin\theta & r \cos\theta & 0 \\ \cos\theta & -r\sin\theta & 0 \\ 0 & 0 & 1 \end{array}\] = r

となるので，dxdydz = r drd\theta dzとなる．これを上に代入すると，

\int_V\delta\phi (\nabla^2\phi+f) r drd\theta dz = 0\;\; \forall \delta\phi\\\Leftrightarrow \int_\theta\int_z\int_r \delta\phi\(\frac{\partial^2\phi}{\partial r^2} + \frac{1}{r}\(\frac{\partial\phi}{\partial r}\) + \frac{\partial^2 \phi}{\partial z^2} + f(r,z)\) r drdzd\theta = 0

ここで，\theta方向の変化が無いことを考慮して，\thetaの積分を除外できる．これをガウスの発散定理を用いて弱形式に変換する

\int_z\int_r \delta\phi\(\frac{\partial^2\phi}{\partial r^2} + \frac{1}{r}\(\frac{\partial\phi}{\partial r}\) + \frac{\partial^2 \phi}{\partial z^2} + f(r,z)\) r drdz = 0\;\;\forall \delta\phi\\ \Leftrightarrow \int_z\int_r \(\frac{\partial}{\partial r}\(\delta\phi r\frac{\partial\phi}{\partial r}\) - r\frac{\partial\delta\phi}{\partial r}\frac{\partial \phi}{\partial r} + \frac{\partial}{\partial z}\(\delta\phi r\frac{\partial \phi}{\partial z}\) - r\frac{\partial\delta\phi}{\partial z}\frac{\partial\phi}{\partial z} + \delta\phi rf(r,z)\) drdz = 0\;\;\forall \delta\phi\\\Leftrightarrow \int_z\int_r r\frac{\partial\delta\phi}{\partial r}\frac{\partial \phi}{\partial r} + r\frac{\partial\delta\phi}{\partial z}\frac{\partial\phi}{\partial z} drdz = \int_z\int_r \frac{\partial}{\partial r}\(\delta\phi r\frac{\partial\phi}{\partial r}\) + \frac{\partial}{\partial z}\(\delta\phi r\frac{\partial \phi}{\partial z}\) + \delta\phi rf(r,z) drdz\;\;\forall \delta\phi\\\Leftrightarrow \int_z\int_r r\frac{\partial\delta\phi}{\partial r}\frac{\partial \phi}{\partial r} + r\frac{\partial\delta\phi}{\partial z}\frac{\partial\phi}{\partial z} drdz = \int_\Gamma n_r\(\delta\phi r\frac{\partial\phi}{\partial r}\) + n_z\(\delta\phi r\frac{\partial \phi}{\partial z}\)ds + \int_z\int_r \delta\phi rf(r,z) drdz\;\;\forall \delta\phi\\\Leftrightarrow \int_z\int_r r\frac{\partial\delta\phi}{\partial r}\frac{\partial \phi}{\partial r} + r\frac{\partial\delta\phi}{\partial z}\frac{\partial\phi}{\partial z} drdz = \int_\Gamma \delta\phi r\frac{\partial\phi}{\partial n}ds + \int_z\int_r\delta\phi rf(r,z) drdz\;\;\forall \delta\phi\\

断面の２次元的な領域を\Omega，\frac{\partial \phi}{\partial n}=tとすると次のようになる．

弱形式化された円柱座標でのラプラス方程式
\int_\Omega r\frac{\partial\delta\phi}{\partial r}\frac{\partial \phi}{\partial r} + r\frac{\partial\delta\phi}{\partial z}\frac{\partial\phi}{\partial z} drdz = \int_\Gamma \delta\phi t rds + \int_\Omega\delta\phi f(r,z) rdrdz\;\;\forall \delta\phi\\

有限要素法離散化

上の弱形式化では軸対称な問題を断面の２次元問題に帰着させた．さて断面の領域がメッシュ分割されているとしてみよう．積分は積分領域を要素単位で分けて足し合わせて計算するとして上の式は以下のようになる．

\sum_e\int_{\Omega_e} r\frac{\partial\delta\phi}{\partial r}\frac{\partial \phi}{\partial r} + r\frac{\partial\delta\phi}{\partial z}\frac{\partial\phi}{\partial z} drdz = \sum_e\int_{\Gamma_e} \delta\phi trds + \sum_e\int_{\Omega_e}\delta\phi f(r,z) rdrdz\;\;\forall \delta\phi\\

ガラーキン法を用いた離散化をするとして，要素の中で\delta\phiと\phiが２次元の断面上で補間関数N(r,z)を使って次のように要素分割されているとする.

\phi = N^n\phi^n,\;\;\;\; \delta\phi = N^m\phi^m

要素の中での番号m,nは全体での番号a,bに対応しているとする．また，自然境界条件が課されている表面の要素は形状関数\bar{N}で離散化されているとする．これを上に代入すると以下のようになる．

\sum\int_{\Omega_e} r\frac{\partial N^m\delta\phi^m}{\partial r}\frac{\partial N^n\phi^n}{\partial r} + r\frac{\partial N^m\delta\phi^m}{\partial z}\frac{\partial N^n\phi^n}{\partial z} drdz = \sum_e\int_{\Gamma_e} \bar{N}^m\delta\phi^m trds + \sum_e\int_{\Omega_e} N^m\delta\phi^m f(r,z) rdrdz\;\;\forall \delta\phi\\ \Leftright \sum_e\[\delta\phi^m\int_\Omega r\frac{\partial N^m}{\partial r}\frac{\partial N^n}{\partial r} + r\frac{\partial N^m}{\partial z}\frac{\partial N^n}{\partial z} drdz\phi^n \] = \sum_e\[\delta\phi^m \int_{\Gamma_e} \bar{N}^mtrds\] + \sum_e\[\delta\phi^m\int_{\Omega_e} N^m f(r,z) rdrdz\]\;\;\forall \delta\phi\\ \Leftright \sum_e\[\delta\phi^m A_e^{mn}\phi^n \] = \sum_e\[\delta\phi^m {f_e}_{out}^m\] + \sum_e\[\delta\phi^mf_e_{in}^m\]\;\;\forall \delta\phi\\ \Leftright \delta\phi^a\[ \sum_e A_e^{mn} \]\phi^b = \delta\phi^a\[ \sum_e {f_e}_{out}^m\] + \delta\phi^a\[ \sum_e {f_e}_{in}^m\]\;\;\forall \delta\phi\\ \Leftright \delta\phi^aA^{ab}\phi^b = \delta\phi^a f_{out}^a + \delta\phi^a f_{in}^a\;\;\forall \delta\phi

但し，A_e^{mn}は要素剛性行列で，{f_e}^mは外力ベクトルであり，次のように書くことができる．

要素剛性行列
A_e^{mn}=\int_\Omega r\frac{\partial N^m}{\partial r}\frac{\partial N^n}{\partial r} + r\frac{\partial N^m}{\partial z}\frac{\partial N^n}{\partial z} drdz
外力ベクトル
{f_e}_{out}^m = \int_{\Gamma_e} \bar{N}^m\delta\phi^m trds
{f_e}_{in}^m = \int_{\Omega_e} N^m\delta\phi^m f(r,z) rdrdz

三角形一次要素（線形要素）における要素剛性行列

三角形一次要素のような単純な場合は解析的に要素剛性行列の積分を行うことができる．このような一次要素による線形な補間関数をNの変わりにLを使って表す．

軸までの距離rも三角形内で一次的に変化するから

r = L^p r^p

というように，節点の値の補間によって得ることができる．

A_e^{mn}=\int_\Omega L^pr^p\frac{\partial L^m}{\partial r}\frac{\partial L^n}{\partial r} + L^pr^p\frac{\partial L^m}{\partial z}\frac{\partial L^n}{\partial z} drdz\\\;\;\;\; = \frac{1}{3}(r^1+r^2+r^3)\{\frac{\partial L^m}{\partial r}\frac{\partial L^n}{\partial r} + \frac{\partial L^m}{\partial z}\frac{\partial L^n}{\partial z} \}\Delta

但し\Deltaは三角形の面積である．結局この場合，重心位置での距離に２次元三角形のラプラス方程式の要素剛性を掛けたものになることが分かる．

参考にしたもの

Books

有限要素法概説―理工学における基礎と応用 (FEM+BEM (3)) 菊池文雄 著
流れと熱伝導の有限要素法入門 矢川元基 著

Made by Nobuyuki UMETANI  梅谷　信行
n.umetani@gmail.com
Generated by ひとりWiki version 1.15.1