有限要素法(FEM)の概略をメモ

重み付き残差法の大まかな流れをメモする。
ここでは
$$ \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f=0 $$を有限要素法で解くことを考える
これは、任意の重み関数$ w $をかけて領域積分した値が0に等しいという以下の式に等価である。
$$ \int_{\Omega}w\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
このまま($w$が任意)では難しいので$w=w_i$　$i=1,2\cdots ,n $について満たす解を求める
領域を少領域に分け、それぞれの領域で（例えば）1次式で解を
$$ u\approx \alpha_1+\alpha_2x+\alpha_3y $$と近似する。ただし、領域ごとに係数は異なる。このままでは扱いにくいので領域に接する節点における物理量で補間して解を表すことを考える。
小領域が三角形とすると、小領域$ e $、接点$j=(x_j,y_j)$における物理量は
$$ u^e_j=\alpha_1+\alpha_2x^e_j+\alpha_3y^e_j $$
なので、
$$
\begin{pmatrix}
u^e_1 \\
u^e_2 \\
u^e_3 \\
\end{pmatrix}
=
\begin{pmatrix}
1&x^e_1&y^e_1 \\
1&x^e_2&y^e_2 \\
1&x^e_3&y^e_3 \\
\end{pmatrix}
\begin{pmatrix}
\alpha_1\\
\alpha_2\\
\alpha_3\\
\end{pmatrix}
$$
となる。
これを解くと、$\alpha_i$が求めることができ、解（の近似）が求められる。つまり、$u^e_j$が求まれば、解が求まることがわかる。
扱いやすいように形を整えるために、
$$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $$を用いて以下のように表すことが出来る。ただし、$a^e_j,b^e_j,c^e_j$は$x^e_i,y^e_i$を用いて表され、既知である。
$$ u\approx=\alpha_1+\alpha_2x+\alpha_3y=\sum_{j=1}^3N^e_j(x,y)u^e_j$$
$N^e_j(x,y)$は、内挿関数や形状関数と呼ばれ、節点の物理量で補間するための重み（割合）であると考えることが出来る。
接点における重みを考えると、
$$ N^e_j(x_i,y_i)=\delta_{ij} $$
であり、任意の点$(x,y) $における解において各接点の割合の和は１になることから、
$$\sum_{j=1}^3N^e_j(x,y)=1$$
となる。
ここで一度原点に戻ると、少領域$e$で
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
となる$u$を求めたかったわけだが、$u=\sum_{j=1}^3N^e_j(x,y)u^e_j$を代入すると
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial x^2}+\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial y^2}+f \right)dS=0 $$
これは以下のように変形できる。
$$
\sum_{j=1}^3 \int_{\Omega_e}w_i\left(\frac{\partial^2 N^e_j}{\partial x^2}+\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j+\int_{\Omega_e}w_ifdS=0
$$

ここからガラーキン法(Galerkin法)を使う。
$$ w_i=-N^e_i(x,y) $$とする。ただし、小領域外では0であり、マイナスは式変形しやすくするために便宜上つけたものである。
すると、
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial x^2}+-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j=\int_{\Omega_e}N^e_ifdS
$$
これを式1とする。

ここで、次に行う式変形のために2次元のガウスの発散定理を思い出す。
$$ \int \mbox{div}{\bf{F}} dS=\int {\bf{F}}\cdot {\bf{n}}dC $$
$ {\bf{F}}=\left(N_i\frac{\partial N_j}{\partial x},N_i\frac{\partial N_j}{\partial y}\right)$を代入すると
$$\int\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+N_i\frac{\partial^2 N_j}{\partial x^2}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+N_i\frac{\partial^2 N_j}{\partial y^2}dS=\int\left(N_i\frac{\partial N_j}{\partial x}n_x+N_i\frac{\partial N_j}{\partial y}n_y\right)dC$$
これを用いると式1は（面倒なので$e$省略）
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dSu^e_j=\int_{\Omega_e}N_ifdS+\int_{\partial\Omega_e}N_i\frac{\partial u}{\partial {\bf{n}}}dC
$$
ただし最後の項は法線方向微分である。
$$ k^e_{ij}=\int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dS$$とすると、$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $より、
$$ k^e_{ij}=\int_{\Omega_e}\left(b_ib_j+c_ic_j\right)dS$$なので計算できる。
$$ \frac{\partial u}{\partial {\bf{n}}} $$は境界条件を設定するのに使える。
$K_e=(k^e_{ij})$とすると、$K_e$は対称行列であり要素剛性行列と呼ばれる。
各要素の$K_e$を上手く合わせれば1次連立方程式に帰着され、$ u^e_j $を求めることが出来る。

風上差分と数値拡散

風上差分については感覚的に流れの上流の情報を使うと説明されることが多いので数学的な説明をメモする
まず、拡散方程式について軽くメモする。
拡散方程式は、以下の2条件を満たす
・濃度が高い方から低い方へ流れる
・流れの速さは濃度の勾配に比例する
$ x-\Delta x,x,x+\Delta x $における濃度$ f(x) $を考えて、$ x $に流れ込む量と流れ出す量を足し合わせると拡散方程式
$$\frac{\partial f}{\partial t}=D\frac{\partial^2 f}{\partial x^2}$$
が得られる　ただし、$D$は正のときに拡散をするが、負だと濃度勾配をさらに大きくすることになる。詳しくは負の拡散係数とかでググれば出てくる

移流方程式
$$ \frac{\partial f}{\partial t}+u\frac{\partial f}{\partial x}=0$$
を数値的に解くことを考える
テイラー展開をすると
$$ f_{i-1}=f_i-\left(\frac{\partial f}{\partial x}\right)_i\Delta x+\frac{(\Delta x)^2}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+\mathcal{O}(\Delta x^3)$$
$$ f_{i+1}=f_i+\left(\frac{\partial f}{\partial x}\right)_i\Delta x+\frac{(\Delta x)^2}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+\mathcal{O}(\Delta x^3)$$
よって
\begin{equation} u\left(\frac{\partial f}{\partial x}\right)_i=u\frac{f_i-f_{i-1}}{\Delta x}+u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)\end{equation}

\begin{equation} u\left(\frac{\partial f}{\partial x}\right)_i=u\frac{f_{i+1}-f_{i}}{\Delta x}-u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)\end{equation}
上二式の二階偏微分の項が$ u $の正負によって拡散項になったり負の拡散項になることが本質である。$u$が正の時は上の1式が風上差分となるがその時にちゃんと拡散項になっていて負の拡散項ではないことを確かめる。

1式を変形すると
$$ u\frac{f_i-f_{i-1}}{\Delta x}=u\left(\frac{\partial f}{\partial x}\right)_i-u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)  $$
移流方程式に1式を用いて移流方程式を離散化すると
$$ \frac{\partial f}{\partial t}+u\frac{f_i-f_{i-1}}{\Delta x}=0$$
を解くことになるが、これは上式を代入すると
$$ \frac{\partial f}{\partial t}+u\frac{\partial f}{\partial x}=u\frac{\Delta x}{2}\left(\frac{\partial^2 f}{\partial x^2}\right)_i+u\mathcal{O}(\Delta x^2)$$
これは、右辺が拡散方程式と同じで拡散項になっていることがわかる。
(2)式を用いて同様のことをした場合(風上差分ではないほうの差分)右辺が負の拡散項になるので濃度勾配に対して物理量が更に濃度勾配が大きくなる作用を持つ。そのため、安定性が悪い。

ラグランジュ補間による差分法の2次精度微分

ラグランジュ補間による差分法の2次精度微分


関数$ f $の数値微分$f'(x_1)$を行う方法のメモ．

$$ (x_0,f_0),(x_1,f_1),(x_2,f_2)$$
が分かるとする($x_0<x_1<x_2$, $ f_0,f_1,f_2$は定数)

上記3点を通る2次関数は
$$g= f_0\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}+f_1\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}+f_2\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}$$
である．代入すると自明である．
$$\Delta x_0=x_1-x_0, \Delta x_1=x_2-x_1$$
とおく．
$g$を微分して$x=x_1$を代入すると
$$ f'(x_1)\approx g'(x_1)$$であり，
$$g'(x_1)=-f_0\frac{\Delta x_1}{\Delta x_0(\Delta x_0+\Delta x_1)}+f_1\frac{(\Delta x_1-\Delta x_0)}{\Delta x_0\Delta x_1} +f_2\frac{\Delta x_0}{\Delta x_1(\Delta x_0+\Delta x_1)}$$
\begin{equation}=\cfrac{f_2(\Delta x_0)^2-f_0(\Delta x_1)^2+f_1((\Delta x_1)^2-(\Delta x_0)^2)}{\Delta x_0 \Delta x_1(\Delta x_0+\Delta x_1)}\end{equation}
となる　これが2次精度の公式である．

次に2次精度であることを確認する．
テイラー展開をして
$$f_2=f_1+\Delta x_1(f'(x_1))+\frac{(\Delta x_1)^2}{2!}f''(x_1)+\mathcal{O}((\Delta x_1)^3)$$
$$f_0=f_1-\Delta x_0(f'(x_1))+\frac{(\Delta x_0)^2}{2!}f''(x_1)+\mathcal{O}((\Delta x_0)^3)$$

これらを(1)に代入する
(1)の分子は
$$ f_1(\Delta x_0)^2+\Delta x_1(\Delta x_0)^2(f'(x_1))+\frac{(\Delta x_1)^2 (\Delta x_0)^2}{2!}f''(x_1)+\mathcal{O}((\Delta x_1)^3(\Delta x_0)^2)-f_1(\Delta x_1)^2+\Delta x_0(\Delta x_1)^2(f'(x_1))-\frac{(\Delta x_0)^2(\Delta x_1)^2}{2!}f''(x_1)-\mathcal{O}((\Delta x_0)^3(\Delta x_1)^2)+f_1((\Delta x_1)^2-(\Delta x_0)^2)
$$
$$=\Delta x_0\Delta x_1(\Delta x_0+\Delta x_1)(f'(x_1))+\mathcal{O}\left((\Delta x_0)^2(\Delta x_1)^2(\Delta x_1+\Delta x_0)\right)$$
よって
$$(1)=f'(x_1)+\mathcal{O}\left(\Delta x_0\Delta x_1\right)$$

なので2次精度を満たす！

レイノルズ平均方程式の導出

乱流を扱うときに使うモデルです．
速度を平均と変動成分に分けます．
レイノルズ分解
$$ u_i=\bar{u}_i+u'_i $$
流れが定常→時間平均
一様な流れ→空間平均
一般→アンサンブル平均

アンサンブル平均
$$\bar{u}=\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{i=1}^Nu_i$$
次を満たすアンサブル平均をレイノルズ平均という
$$　\overline{u'}_i=0$$
$$\bar{\bar{u}}_i=\bar{u}_i$$
$$\overline{u'_i\bar{u}_j}=0$$
次が成りたつ
$$ \overline{\frac{\partial u_i}{\partial x_j}}=\frac{\partial \bar{u}_i}{\partial x_j}$$
$$ \overline{\frac{\partial u'_i}{\partial x_j}}=\frac{\partial \bar{u'}_i}{\partial x_j}=0$$
非圧縮流体の連続の式は$$ \frac{\partial u_i}{\partial x_i}=0$$
なのでレイノルズ平均を取ると$$\frac{\partial \overline{u}_i}{\partial x_i}=0$$
両辺引くと$$\frac{\partial u'_i}{\partial x_i}=0$$
ナビエ・ストークス方程式は
$$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\nu \frac{\partial^2 u_i}{\partial x_j\partial x_j} $$
$ u_i=\bar{u}_i+u'_i $を代入すると
$$\frac{\partial (\bar{u}_i+u'_i)}{\partial t}+(\overline{u}_j+u'_j)\frac{\partial (\bar{u}_i+u'_i)}{\partial x_j}=-\frac{1}{\rho}\frac{\partial p}{\partial x_i}+\nu \frac{\partial^2 (\bar{u}_i+u'_i)}{\partial x_j\partial x_j} $$
レイノルズ平均を取る．
$$\frac{\partial \bar{u}'_i}{\partial t}=0$$

随時更新

テンソルを用いたナビエストークス方程式の導出

テンソルを使ってナビエストークス方程式の導出手順についてメモします．

①
応力テンソル＝圧力＋粘性応力テンソル
応力テンソルの対角成分は面に垂直にかかる力なので応力テンソルのトレースの平均を圧力とする(外向き正)
つまり，応力テンソル＝平均応力テンソル＋偏差応力テンソルを考えている

$$
\mbox{応力テンソル　}\sigma= \left(
    \begin{array}{ccc}
      -p & 0 & 0 \\
      0 & -p & 0 \\
      0& 0 & -p
    \end{array}
  \right) +\mbox{粘性応力テンソル　}\tau
$$

②
ニュートン流体は偏差応力テンソル（粘性応力テンソル，せん断応力)が変形速度テンソル(ひずみ率テンソル)の線形結合で表せる　（剛体回転するときの応力は圧力のみ）

速度勾配テンソル${\bm{l}}$
$$ l_{ij}=\frac{\partial u_i}{\partial x_j} $$
直交分解する
$ {\bm{d}} $を変形速度テンソル(ひずみ速度)，${\bm{w}}$を回転速度テンソル（スピンテンソル）として$$ {\bm{d}}=\frac{{\bm{l}}+{\bm{l}}^t}{2} $$
成分は
$$ d_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) $$

となり$$ {\bm{w}}=\frac{{\bm{l}}-{\bm{l}}^t}{2} $$である
$ {\bm{d}} $は引き延ばしを，${\bm{w}}$ は回転を表す

流体は座標変換しても変わらないという等方性を仮定すると，変形速度テンソルの係数は4階の等方テンソルで表せる
4階の等方テンソル$ A $は次の形でかける
$$ A_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu_1\delta_{ik}\delta_{jl}+\mu_2\delta_{il}\delta_{jk} $$
変形速度テンソルは対称なので
$$\tau_{ij}=A_{ijkl}d_{kl}$$
$$=\lambda \delta_{ij}d_{kk}+(\mu_1+\mu_2) d_{ij}$$
$$=\lambda \delta_{ij}d_{kk}+2\mu d_{ij}$$
\begin{equation}=\lambda \delta_{ij}(d_{11}+d_{22}+d_{33})+2\mu d_{ij}\end{equation}
とかける
偏差応力テンソルのトレース＝０なので$$ \tau_{11}+\tau_{22}+\tau_{33}=0$$
から(1)をつかって
$$(3\lambda+2\mu)\mbox{div}{\bm{u}}=0  $$
よってストークスの仮定を使うときは$$\lambda=-\frac{2}{3}\mu$$とおく
以上より応力テンソルは
$$\sigma_{ij}=-p\delta_{ij}-\frac{2}{3}\mu\delta_{ij}\mbox{div}{\bm{u}}+\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$
非圧縮なら連続の式から$ div u=0 $なので
$$\sigma_{ij}=-p\delta_{ij}+\mu\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$となる

コーシーの運動方程式から(もしくはテンソルの力の向きを考えて)x成分を考えると
$$\rho \frac{Du_x}{Dt}=f_{x}+\frac{\partial \sigma_{xx}}{\partial x}+\frac{\partial \sigma_{xy}}{\partial y}+\frac{\partial \sigma_{xz}}{\partial z}$$
$$=f_{x}+\frac{\partial}{\partial x}(-p+2\mu\frac{\partial u_x}{\partial x}-\frac{2}{3}\mu\mbox{div}{\bm{u}})+\mu\frac{\partial}{\partial y}(\frac{\partial u_y}{\partial x}+\frac{\partial u_x}{\partial y})+\mu\frac{\partial}{\partial z}(\frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z})$$

$$ =f_{x}-\frac{\partial p}{\partial x}+\mu(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2})u_x+\frac{\mu}{3}\frac{\partial}{\partial x}\mbox{div}{\bm{u}}$$

よって導けた
非圧縮の場合はdiv u=0なので
$$\rho \frac{D{\bm{u}}}{Dt}={\bm{f}}-\nabla p+\mu \nabla^2 {\bm{u}} $$

圧縮性流体の運動方程式(運動量保存則)

（一般の）圧縮性流体の運動方程式は、テンソルが出てきて分かりにくいので、レイノルズの輸送定理を用いた具体的な導出をメモします。

レイノルズの輸送定理
$$ \int_{B_t}f({\bm{x}},t) dv $$
空間座標でtの変化とともに積分領域が変化する時に上の物質微分をしたいとすると、積分領域が変化するのでそのままでは上手くいかない。
一度物質座標に変換（変形勾配テンソルを使う）すれば、物質座標では積分領域が変化しない。なぜなら、tに合わせて空間座標における物質座標も変化するからである。（ゴムをのばすとゴム上の座標ものびるイメージ）
なので、物質座標で積分してから、また空間座標に戻せば良い。すると以下の結果が得られる。
$$
\frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv= \int_{B_t}(\dot{f}({\bm{x},t})+f({\bm{x},t})\mbox{div}{\bm{v}})dv
$$
$ f $の物質微分を考えて、
$$
\frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv= \int_{B_t}(\frac{\partial f}{\partial t}+{\bm{v}}\cdot\mbox{grad}f+f({\bm{x},t})\mbox{div}{\bm{v}})dv $$
よって
$$ \frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv=\int_{B_t}(\frac{\partial f}{\partial t}+\mbox{div}(f{\bm{v}})) dv
$$

運動量変化は以下のように表される
$$ \frac{D}{Dt}\int_{V}\rho{\bm{v}}dv $$
速さがベクトルだと扱いづらいので$ {\bm{e_k}} $成分について考える。レイノルズの輸送定理をつかうと
\begin{equation} \frac{D}{Dt}\int_{V}\rho{v_k}dv=\int_V(\frac{\partial \rho v_k}{\partial t}+\frac{\partial \rho v_k v_j}{\partial x_j})dv\end{equation}

また、$$ \nabla\cdot(\rho {\bm{v}}\otimes{\bm{v}})= ({\bm{e_i}}\frac{\partial }{\partial x_i})\cdot (\rho v_j v_k{\bm{e_j}}\otimes{\bm{e_k}})$$
$$=\frac{\partial \rho v_j v_k}{\partial x_j}{\bm{e}}_k$$
は(1)の右辺第二項と一致する
よって圧縮性流体の運動方程式の右辺は
$$
\frac{\partial\rho{\bm{u}}}{\partial t}+\nabla\cdot(\rho{\bm{v}}\otimes{\bm{v}})=\mbox{外部からの力積}
$$