定常な移流拡散に対する有限要素法

初めに

もう当面，P1のみ考えることにした（cf. 1D Poisson, 2D Poisson）．

参考:

• 有限要素法による流れのシミュレーション
• 続・有限要素法による流れのシミュレーション
• 第3版 有限要素法による流れのシミュレーション

最新版を読めば事足りることはなくて，古いものにしか書いていないことも（逆も）沢山あり，全部手元に置いておきたい．

MAIN

定常な移流拡散方程式

以下である．

\begin{alignat}{2}
    \boldsymbol{\beta} \cdot \nabla u
    - \nabla \cdot \alpha \nabla u
    &= f
    \quad &&\text{in} \quad \Omega
    \\
    u
    &= g_D
    \quad &&\text{on} \quad \Gamma_D
    \\
    \boldsymbol{n} \cdot \alpha \nabla u
    &= g_N
    \quad &&\text{on} \quad \Gamma_N
\end{alignat}

ここで，$\boldsymbol{\beta}$は移流速度，$\alpha$は拡散率．

弱形式

適当な$v \in H_{\Gamma_D, 0}^1$を取り，

\begin{align}
    \int_{\Omega} \left(
        \boldsymbol{\beta} \cdot \nabla u
        - \nabla \cdot \alpha \nabla u
    \right) v
    \ d\boldsymbol{x}
    &= \int_{\Omega}
        f v
    \ d\boldsymbol{x}
    \\
    \int_{\Omega} \left(
        \boldsymbol{\beta} \cdot \nabla u
    \right) v
    \ d\boldsymbol{x}
    + \int_{\Omega}
        \alpha \nabla u \cdot \nabla v
    \ d\boldsymbol{x}
    &= \int_{\Omega}
        f v
    \ d\boldsymbol{x}
    + \int_{\Gamma_N}
        g_N v
    \ d\boldsymbol{s}
\end{align}

を考える（Gauss-Greenの式を使った）．

Galerkin

以降，簡単のため，1次元のDirichlet問題とし，$f=0$とする．また，表記の簡略化のため，暗に$u = u_h$などとする．

$\bar{\Omega} = \bigcup_k e_k$を集めて，

\begin{align}
    \sum_{k=1}^{K} \int_{e_k}
        \beta \frac{du}{dx} v
    \ dx
    + \sum_{k=1}^{K} \int_{e_k}
        \alpha \frac{du}{dx} \frac{dv}{dx}
    \ dx
    &= 0
\end{align}

を考える．要素$e_k = [x_{k-1}, x_{k}]$内でP1補間を行えば（要素長: $h_k$），

\begin{align}
    u (x)
    &= \sum_{j=1}^{2} N_j (x) u_j
    \\
    v (x)
    &= \sum_{i=1}^{2} N_i (x) v_i
    \\
    N_1 (x)
    &= \frac{x_{k} - x}{h_k}
    \\
    N_2 (x)
    &= \frac{x - x_{k-1}}{h_k}
\end{align}

より，

\begin{align}
    \int_{e_k}
        \beta \frac{du}{dx} v
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{dv}{dx}
    \ dx
    &= 0
    \\
    \int_{e_k}
        \beta \frac{d (N_j u_j)}{dx} (N_i v_i)
    \ dx
    + \int_{e_k}
        \alpha \frac{d (N_j u_j)}{dx} \frac{d (N_i v_i)}{dx}
    \ dx
    &= 0
    \\
    v_i \left[
    \int_{e_k}
        \beta \frac{d N_j}{d x} N_i
    \ dx \right]
    u_j
    + v_i \left[
    \int_{e_k}
        \alpha \frac{d N_j}{dx} \frac{d N_i}{dx}
    \ dx \right]
    u_j
    &= 0
\end{align}

$v$の任意性から，結局，

\begin{align}
    \boldsymbol{A}_k \boldsymbol{u}_k 
    + \boldsymbol{D}_k \boldsymbol{u}_k 
    &= 0
    \\
    \left( \boldsymbol{A}_k + \boldsymbol{D}_k \right) \boldsymbol{u}_k 
    &= 0
    \\
    \therefore
    \boldsymbol{K}_k \boldsymbol{u}_k 
    &= 0
    \quad (\boldsymbol{K}_k = \boldsymbol{A}_k + \boldsymbol{D}_k)
\end{align}

但し，

\begin{align}
    \boldsymbol{A}_k
    &= \left[
    \int_{e_k}
        \beta \frac{d N_j}{d x} N_i
    \ dx
    \right]
    = \beta \begin{bmatrix}
        -1/2 & 1/2 \\
        -1/2 & 1/2
    \end{bmatrix}
    \\
    \boldsymbol{D}_k
    &= \left[
    \int_{e_k}
        \alpha \frac{d N_j}{dx} \frac{d N_i}{dx}
    \ dx \right]
    = \alpha \begin{bmatrix}
         1/h_k & -1/h_k \\
        -1/h_k &  1/h_k
    \end{bmatrix}
\end{align}

を集めて，

\begin{align}
    \left( \boldsymbol{A} + \boldsymbol{D} \right) \boldsymbol{u}
    &= 0
    \\
    \boldsymbol{K} \boldsymbol{u}
    &= 0
    \quad (\boldsymbol{K} = \boldsymbol{A} + \boldsymbol{D})
\end{align}

を解くことである．

いま，$x \in [0, 1]$において，$u (0) = 0, \ u (1) = 1$という境界条件を与え，

\begin{align}
    u
    &= \frac{\exp (\text{Pe} \cdot x) - 1}{\exp (\text{Pe}) - 1}
\end{align}

という厳密解を与える．但し，$\text{Pe} = |\beta| L / \alpha$はPeclet数．また，要素Peclet数を$\text{Pe}_k = (|\beta| h_k) / (2 \alpha)$と定めた上で，移流速度を$\beta = 1$に固定し，拡散率と解像度を変動させたとき，Galerkin型有限要素解は以下のようになっている．

$\text{Pe}=25$	$\text{Pe}=50$	$\text{Pe}=100$

Peclet数は解自体の様子を定めているが，有限要素解の視点からは然して重要でない．どちらかと言うと，要素Peclet数が1を超えるか否かによって有限要素解の安定性が様変わりしている（何だか，まるでCourant数のような？）．要素長を次々に小さくしてゆけば，$\text{Pe}_k < 1$を満たすことは可能であろうが，計算資源の制約から非現実的である．有限差分法と同様に，移流に対して最小限の粘性を加える安定化を求める．

Upwind Galerkin

ところで，先の$\boldsymbol{K}_k \boldsymbol{u}_k = 0$は，要素$e_k$での式であった．

$e_k$と$e_{k+1}$に跨って調べてみると，

\begin{align}
    \boldsymbol{K}_k \boldsymbol{u}_k
    &= 0
    \\
    \boldsymbol{K}_{k+1} \boldsymbol{u}_{k+1}
    &= 0
\end{align}

があって，真ん中の行を取り出すと，

\begin{align}
    \boldsymbol{K}_{k}^{21} u_{k-1}
    + (\boldsymbol{K}_{k}^{22} + \boldsymbol{K}_{k+1}^{11}) u_{k}
    + \boldsymbol{K}_{k+1}^{12} u_{k+1}
    &= 0
    \\
    \beta \frac{u_{k+1} - u_{k-1}}{2 h_k}
    - \alpha \frac{u_{k+1} - 2 u_k + u_{k-1}}{h_k^2}
    &= 0
\end{align}

が見える．中心差分は非常に不安定で，有限差分法では，風上化を施して数値的安定性を得たのであった（cf. 有限差分法）．例えば，1次精度風上差分では，

\begin{align}
    \alpha
    \to
    \alpha + \bar{\alpha};
    \;
    \bar{\alpha}
    &= \frac{|\beta| h_k}{2}
\end{align}

とした．これを引き継ぎ，false viscosityを入れよう，というのがupwind Galerkinである．が，当然受け入れられないよね．

=== 独り言ここから ===
false viscosityという言葉は非常に便利だと思っている．この$\bar{\alpha}$は，numerical viscとかartificial viscとか呼ばれるようだが，私は，numerical viscとartificial viscは異なると思っている．私の中の区別としては，有限差分法のfalse viscはnumericalだが，upwind Galerkinのfalse viscはartificialだと考える．だって有限差分法では移流項を数値的に離散化した結果"現れてしまった"ものだったのに対し，upwind Galerkinでは元の問題を"人間が積極的に"書き換えている．両者は形式的に同じに見えるが，登場するまでの意味が違うように感じる．なんて今の時代言っても，「どうでもいい」で一蹴されるんだろうけど．
=== 独り言ここまで ===

Upwind Petrov-Galerkin

Upwind Galerkinは，

\begin{align}
    \beta \frac{u_{k+1} - u_{k-1}}{2 h_k}
    - (\alpha + \bar{\alpha}) \frac{u_{k+1} - 2 u_k + u_{k-1}}{h_k^2}
    &= 0
\end{align}

であった．議論を逆向きに辿れば，これは，

\begin{align}
    \boldsymbol{A}_k \boldsymbol{u}_k 
    + \tilde{\boldsymbol{D}}_k \boldsymbol{u}_k 
    &= 0
    \\
    \boldsymbol{A}_k
    &= \left[
    \int_{e_k}
        \beta \frac{d N_j}{d x} N_i
    \ dx
    \right]
    \\
    \tilde{\boldsymbol{D}}_k
    &= \left[
    \int_{e_k}
        (\alpha + \bar{\alpha}) \frac{d N_j}{dx} \frac{d N_i}{dx}
    \ dx \right]
\end{align}

を解いていることに対応しており，さらに辿れば，

\begin{align}
    \int_{e_k}
        \beta \frac{du}{dx} v
    \ dx
    + \int_{e_k}
        (\alpha + \bar{\alpha}) \frac{du}{dx} \frac{dv}{dx}
    \ dx
    &= 0
\end{align}

を解いている（false viscで元の系を書き換えたと解釈して良い）．整理すると，

\begin{align}
    (\text{LHS})
    &= \int_{e_k}
        \beta \frac{du}{dx} v
    \ dx
    + \int_{e_k}
        (\alpha + \bar{\alpha}) \frac{du}{dx} \frac{dv}{dx}
    \ dx
    \\
    &= \int_{e_k}
        \beta \frac{du}{dx} v
    \ dx
    + \int_{e_k}
        \bar{\alpha} \frac{du}{dx} \frac{dv}{dx}
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{dv}{dx}
    \ dx
    \\
    &= \int_{e_k}
        \beta \frac{du}{dx} \left( v + \frac{\bar{\alpha}}{\beta} \frac{dv}{dx} \right)
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{dv}{dx}
    \ dx
\end{align}

を解いていることになっている．即ち，移流についてのみテストの書き換えが起きている．人々はこれを不自然と言う．私には自然か不自然か分からないが，初めの弱形式化の手続きの中で，両辺に同一のテストを乗じたことを思い出せば不自然と言えるような気もする．

こういうことを思い出して，全ての項に同一の試験関数を採用するのだそうだ．即ち，（突然だがソースも導入して）

\begin{align}
    \int_{e_k}
        \beta \frac{du}{dx} \tilde{v}
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{d\tilde{v}}{dx}
    \ dx
    - \int_{e_k}
        f \tilde{v}
    \ dx
    &= 0
\end{align}

但し，

\begin{align}
    \tilde{v}
    &= v + \frac{\bar{\alpha}}{\beta} \frac{dv}{dx}
    \\
    &= v + \tau_k \beta \frac{dv}{dx}
\end{align}

試行と試験が異なるのでPetrov-Galerkin型の枠組みになっている（ソースもupwindしなければ上手くいかないことがLeonardから報告されており，現代ではもう知られた事実である）．$\tau_k$は，要素$e_k$毎に指定する安定化パラメータである．このPetrov-Galerkinの方針を，再度展開してゆく．

\begin{align}
    (\text{LHS})
    &= \int_{e_k}
        \beta \frac{du}{dx} \tilde{v}
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{d\tilde{v}}{dx}
    \ dx
    - \int_{e_k}
        f \tilde{v}
    \ dx
    \\
    &= \int_{e_k}
        \beta \frac{du}{dx} \left(
            v + \tau_k \beta \frac{dv}{dx}
        \right)
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{d}{dx} \left(
            v + \tau_k \beta \frac{dv}{dx}
        \right)
    \ dx
    - \int_{e_k}
        f \left(
            v + \tau_k \beta \frac{dv}{dx}
        \right)
    \ dx
    \\
    &= \int_{e_k}
        \beta \frac{du}{dx} v
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{dv}{dx}
    \ dx
    - \int_{e_k}
        f v
    \ dx
    \\
    &\quad+ \tau_k \int_{e_k}
        \left(
            \beta \frac{du}{dx}
            - \alpha \frac{d^2u}{dx^2}
            - f
        \right)
        \left(
            \beta \frac{dv}{dx}
        \right)    
    \ dx
\end{align}

結局，

\begin{align}
    &
    \int_{e_k}
        \beta \frac{du}{dx} v
    \ dx
    + \int_{e_k}
        \alpha \frac{du}{dx} \frac{dv}{dx}
    \ dx
    - \int_{e_k}
        f v
    \ dx
    \\
    &\quad+ \tau_k \int_{e_k}
        \left(
            \beta \frac{du}{dx}
            - \alpha \frac{d^2u}{dx^2}
            - f
        \right)
        \left(
            \beta \frac{dv}{dx}
        \right)    
    \ dx
    \\
    &\quad = 0
\end{align}

を解いているのだと気付く．$\tau_k$を乗じていない部分は通常のGalerkinになっており，$\tau_k$を乗じた部分は重み付き残差である．陽なfalse viscで安定化させるのではなく，要素内における方程式の残差が安定化させる働きを担う．

P1で具体的に書き下すと，P1なので安定化項の内側にある粘性項は消えてしまって，

\begin{align}
    \boldsymbol{A}_k \boldsymbol{u}_k 
    + \check{\boldsymbol{A}}_k \boldsymbol{u}_k 
    + \boldsymbol{D}_k \boldsymbol{u}_k 
    &= \boldsymbol{S}_k
    + \check{\boldsymbol{S}}_k
    \\
    \left( \boldsymbol{A}_k + \check{\boldsymbol{A}}_k + \boldsymbol{D}_k \right)
    \boldsymbol{u}_k 
    &= \boldsymbol{S}_k
    + \check{\boldsymbol{S}}_k
\end{align}

但し，

\begin{align}
    \check{\boldsymbol{A}}_k
    &= \left[ \tau_k \int_{e_k}
        \beta \frac{d N_j}{d x}
        \beta \frac{d N_i}{d x}
    \ dx
    \right]
    = \tau_k \beta^2 \begin{bmatrix}
         1/h_k & -1/h_k \\
        -1/h_k &  1/h_k
    \end{bmatrix}
    \\
    \boldsymbol{S}_k
    &= \left[ \int_{e_k}
        f N_i
    \ dx
    \right]
    = \begin{bmatrix}
         f h_k / 2 \\
         f h_k / 2
    \end{bmatrix}
    \\
    \check{\boldsymbol{S}}_k
    &= \left[ \tau_k \int_{e_k}
        f \beta \frac{d N_i}{d x}
    \ dx
    \right]
    = \tau_k \beta \begin{bmatrix}
         - f \\
           f
    \end{bmatrix}
\end{align}

である．$\check{\boldsymbol{A}}_k$は$\boldsymbol{D}_k$に良く似ている．安定化パラメータは，1次元の問題では

\begin{align}
    \tau_k
    &= \frac{h_k}{2 |\beta|} \xi (\text{Pe}_k) \\
    \xi (\text{Pe}_k)
    &= \coth (\text{Pe}_k) - \frac{1}{\text{Pe}_k} \\
\end{align}

が知られている．

さて，先と同じく，

\begin{align}
    u
    &= \frac{\exp (\text{Pe} \cdot x) - 1}{\exp (\text{Pe}) - 1}
\end{align}

が厳密解となる定常な移流拡散を解く．

$\text{Pe}=25$	$\text{Pe}=50$	$\text{Pe}=100$

驚くべき結果である．

Streamline Upwind/Petrov-Galerkin

（何故スラッシュなんだろうという疑問がずっとある．．．まあいいか）

1次元の問題を考えてきたが，多次元でも解けるようになりたい．取り敢えず2次元が解ければ満足である．混合境界値問題まで広げる．

\begin{alignat}{2}
    \boldsymbol{\beta} \cdot \nabla u
    - \nabla \cdot \alpha \nabla u
    &= f
    \quad &&\text{in} \quad \Omega
    \\
    u
    &= g_D
    \quad &&\text{on} \quad \Gamma_D
    \\
    \boldsymbol{n} \cdot \alpha \nabla u
    &= g_N
    \quad &&\text{on} \quad \Gamma_N
\end{alignat}

テスト$\tilde{v}$を以て，

\begin{align}
    \int_{\Omega} \left(
        \boldsymbol{\beta} \cdot \nabla u
    \right) \tilde{v}
    \ d\boldsymbol{x}
    + \int_{\Omega}
        \alpha \nabla u \cdot \nabla \tilde{v}
    \ d\boldsymbol{x}
    - \int_{\Omega}
        f \tilde{v}
    \ d\boldsymbol{x}
    &= \int_{\Gamma_N}
        g_N \tilde{v}
    \ d\boldsymbol{s}
\end{align}

と弱形式化する．$\tilde{v}$は，1次元の自然な拡張になっており，

\begin{align}
    (1\text{D})\quad
    \tilde{v}
    &= v + \tau_k \beta \frac{dv}{dx}
    \\
    (n\text{D})\quad
    \tilde{v}
    &= v + \tau_k \beta_i \frac{\partial v}{\partial x_i}
\end{align}

である．これを，Galerkin項と安定化項に分けるように整理する．

\begin{align}
    (\text{LHS})
    &= \int_{\Omega} \left(
        \beta_j \frac{\partial u}{\partial x_j}
    \right) \tilde{v}
    \ d\boldsymbol{x}
    + \int_{\Omega}
        \alpha \frac{\partial u}{\partial x_j} \frac{\partial \tilde{v}}{\partial x_j}
    \ d\boldsymbol{x}
    - \int_{\Omega}
        f \tilde{v}
    \ d\boldsymbol{x}
    \\
    &= \int_{\Omega} \left(
        \beta_j \frac{\partial u}{\partial x_j}
    \right) v
    \ d\boldsymbol{x}
    + \int_{\Omega}
        \alpha \frac{\partial u}{\partial x_j} \frac{\partial v}{\partial x_j}
    \ d\boldsymbol{x}
    - \int_{\Omega}
        f v
    \ d\boldsymbol{x}
    \\
    &\quad+ \sum_k \tau_k \int_{e_k} \left(
        \beta_j \frac{\partial u}{\partial x_j}
        - \alpha \frac{\partial^2 u}{\partial x_j \partial x_j}
        - f
    \right)
    \left(
        \beta_i \frac{\partial v}{\partial x_i}
    \right)
    \ d\boldsymbol{x}
    \\
    &\quad+ \tau_k \int_{\Gamma_N}
        \alpha \frac{\partial u}{\partial x_j} n_j
        \left(
            \beta_i \frac{\partial v}{\partial x_i}
        \right)
    \ d\boldsymbol{s}
    \\
    (\text{RHS})
    &= \int_{\Gamma_N}
        g_N \tilde{v}
    \ d\boldsymbol{s}
    \\
    &= \int_{\Gamma_N}
        g_N v
    \ d\boldsymbol{s}
    + \tau_k \int_{\Gamma_N}
        g_N \left( \beta_i \frac{\partial v}{\partial x_i} \right)
    \ d\boldsymbol{s}
\end{align}

両辺にNeumann条件と風上重みの積が登場するのでそれぞれが打ち消して，結局，

\begin{align}
    &
    \int_{\Omega}
        \beta_j \frac{\partial u}{\partial x_j} v
    \ d\boldsymbol{x}
    + \int_{\Omega}
        \alpha \frac{\partial u}{\partial x_j} \frac{\partial v}{\partial x_j}
    \ d\boldsymbol{x}
    - \int_{\Omega}
        f v
    \ d\boldsymbol{x}
    \\
    &\quad+ \sum_k \tau_k \int_{e_k} \left(
        \beta_j \frac{\partial u}{\partial x_j}
        - \alpha \frac{\partial^2 u}{\partial x_j \partial x_j}
        - f
    \right)
    \left(
        \beta_i \frac{\partial v}{\partial x_i}
    \right)
    \ d\boldsymbol{x}
    \\
    &\quad= \int_{\Gamma_N}
        g_N v
    \ d\boldsymbol{s}
\end{align}

が解くべき式である．これまで通り，Galerkinの枠組みに安定化項が付け加えられた形になっている．

再度P1補間を考える．即ち，2次元で，

\begin{align}
    u (x, y)
    &= \sum_{q=1}^{3} N_q (x, y) u_q
    \quad ((x, y) \in e_k)
    \\
    v (x, y)
    &= \sum_{p=1}^{3} N_p (x, y) v_p
    \quad ((x, y) \in e_k)
\end{align}

を与える．$N_r \ (r = 1, 2, 3)$はP1要素の形状関数で，節点$V_r$で1, その他の節点で0となるように選び，

\begin{align}
    N_r (x, y)
    &= a_r + b_r x + c_r y
\end{align}

に対して，

\begin{align}
    a_p
    &= \frac{1}{2 |e_k|} \left( x_q y_r - x_r y_q \right)
    \\
    b_p
    &= \frac{1}{2 |e_k|} \left( y_q - y_r \right)
    \\
    c_p
    &= \frac{1}{2 |e_k|} \left( x_q - x_r \right)
\end{align}

但し，$(p,q,r) = (1,2,3), (2,3,1), (3,1,2)$で回る．また，$|e_k|$は要素$e_k$の面積であり，2D Poissonのときと同様に，相対位置ベクトルの外積を利用して

\begin{align}
    |e_k|
    &= \frac{1}{2} 
    \det
    \begin{bmatrix}
        x_2 - x_1 & x_3 - x_1 \\
        y_2 - y_1 & y_3 - y_1 \\
    \end{bmatrix}
    \\
    &= \frac{1}{2} \left(
        (x_2 - x_1) (y_3 - y_1) - (x_3 - x_1) (y_2 - y_1)
    \right)
\end{align}

で与える．$u, \ v$にP1補間を与えれば，やはり安定化項の中の粘性は消えてしまって，

\begin{align}
    \int_{e_k}
        \beta_j \frac{\partial u}{\partial x_j} v
    \ d\boldsymbol{x}
    &= v_p \left[ \int_{e_k}
        \beta_j \frac{\partial N_q}{\partial x_j} N_p
    \ d\boldsymbol{x} \right]
    u_q
    \\
    \int_{e_k}
        \alpha \frac{\partial u}{\partial x_j} \frac{\partial v}{\partial x_j}
    \ d\boldsymbol{x}
    &= v_p \left[ \int_{e_k}
        \alpha \frac{\partial N_q}{\partial x_j} \frac{\partial N_p}{\partial x_j}
    \ d\boldsymbol{x} \right] u_q
    \\
    \int_{e_k}
        f v
    \ d\boldsymbol{x}
    &= v_p \left[ \int_{e_k}
        f N_p
    \ d\boldsymbol{x} \right]
    \\
    \int_{e_k} \left(
        \beta_j \frac{\partial u}{\partial x_j}
    \right)
    \left(
        \beta_i \frac{\partial v}{\partial x_i}
    \right)
    \ d\boldsymbol{x}
    &= v_p \left[
        \int_{e_k} \left(
            \beta_j \frac{\partial N_q}{\partial x_j}
        \right)
        \left(
            \beta_i \frac{\partial N_p}{\partial x_i}
        \right)
    \ d\boldsymbol{x} \right] u_q
    \\
    \int_{e_k}
        f
    \left(
        \beta_i \frac{\partial v}{\partial x_i}
    \right)
    \ d\boldsymbol{x}
    &= v_p \left[ \int_{e_k}
        f
    \left(
        \beta_i \frac{\partial N_p}{\partial x_i}
    \right)
    \ d\boldsymbol{x} \right]
    \\
    \int_{\Gamma_N}
        g_N v
    \ d\boldsymbol{s}
    &= v_p \left[ \int_{\Gamma_N}
        g_N N_p
    \ d\boldsymbol{s} \right]
\end{align}

を集めて，

\begin{align}
    &v_p \left[ \int_{e_k}
        \beta_j \frac{\partial N_q}{\partial x_j} N_p
    \ d\boldsymbol{x} \right]
    u_q
    + v_p \left[ \int_{e_k}
        \alpha \frac{\partial N_q}{\partial x_j} \frac{\partial N_p}{\partial x_j}
    \ d\boldsymbol{x} \right] u_q
    + v_p \left[ \tau_k 
        \int_{e_k} \left(
            \beta_j \frac{\partial N_q}{\partial x_j}
        \right)
        \left(
            \beta_i \frac{\partial N_p}{\partial x_i}
        \right)
    \ d\boldsymbol{x} \right] u_q
    \\
    &= v_p \left[ \int_{e_k}
        f N_p
    \ d\boldsymbol{x} \right]
    + v_p \left[ \tau_k \int_{e_k}
        f
    \left(
        \beta_i \frac{\partial N_p}{\partial x_i}
    \right)
    \ d\boldsymbol{x} \right]
    + v_p \left[ \int_{\Gamma_N}
        g_N N_p
    \ d\boldsymbol{s} \right]
\end{align}

であり，やはり$v_p$の任意性から，

\begin{align}
    &
    \left[ \int_{e_k}
        \beta_j \frac{\partial N_q}{\partial x_j} N_p
    \ d\boldsymbol{x} \right]
    u_q
    + \left[ \int_{e_k}
        \alpha \frac{\partial N_q}{\partial x_j} \frac{\partial N_p}{\partial x_j}
    \ d\boldsymbol{x} \right] u_q
    + \left[ \tau_k 
        \int_{e_k} \left(
            \beta_j \frac{\partial N_q}{\partial x_j}
        \right)
        \left(
            \beta_i \frac{\partial N_p}{\partial x_i}
        \right)
    \ d\boldsymbol{x} \right] u_q
    \\
    &= \left[ \int_{e_k}
        f N_p
    \ d\boldsymbol{x} \right]
    + \left[ \tau_k \int_{e_k}
        f
    \left(
        \beta_i \frac{\partial N_p}{\partial x_i}
    \right)
    \ d\boldsymbol{x} \right]
    + \left[ \int_{\Gamma_N}
        g_N N_p
    \ d\boldsymbol{s} \right]
\end{align}

より，

\begin{align}
    \boldsymbol{A}_k \boldsymbol{u}_k 
    + \check{\boldsymbol{A}}_k \boldsymbol{u}_k 
    + \boldsymbol{D}_k \boldsymbol{u}_k 
    &= \boldsymbol{S}_k
    + \check{\boldsymbol{S}}_k
    + \boldsymbol{Q}_k
    \\
    \left( \boldsymbol{A}_k + \check{\boldsymbol{A}}_k + \boldsymbol{D}_k \right)
    \boldsymbol{u}_k 
    &= \boldsymbol{S}_k
    + \check{\boldsymbol{S}}_k
    + \boldsymbol{Q}_k
\end{align}

但し，

\begin{align}
    \boldsymbol{A}_k
    &= \left[ \int_{e_k}
        \beta_j \frac{\partial N_q}{\partial x_j} N_p
    \ d\boldsymbol{x} \right]
    \\
    &= \left[ \int_{e_k}
        \beta_1 \frac{\partial N_q}{\partial x_1} N_p
        + \beta_2 \frac{\partial N_q}{\partial x_2} N_p
    \ d\boldsymbol{x} \right]
    \\
    &= \frac{|e_k|}{3} \begin{bmatrix}
        \beta_1 b_1 + \beta_2 c_1 & \beta_1 b_2 + \beta_2 c_2 & \beta_1 b_3 + \beta_2 c_3 \\
        \beta_1 b_1 + \beta_2 c_1 & \beta_1 b_2 + \beta_2 c_2 & \beta_1 b_3 + \beta_2 c_3 \\
        \beta_1 b_1 + \beta_2 c_1 & \beta_1 b_2 + \beta_2 c_2 & \beta_1 b_3 + \beta_2 c_3
    \end{bmatrix}
    \\
    \check{\boldsymbol{A}}_k
    &= \left[ \tau_k 
        \int_{e_k} \left(
            \beta_j \frac{\partial N_q}{\partial x_j}
        \right)
        \left(
            \beta_i \frac{\partial N_p}{\partial x_i}
        \right)
    \ d\boldsymbol{x} \right]
    \\
    &= \left[ \tau_k 
        \int_{e_k} \left(
            \beta_1 \frac{\partial N_q}{\partial x_1}
            \beta_1 \frac{\partial N_p}{\partial x_1}
            +
            \beta_1 \frac{\partial N_q}{\partial x_1}
            \beta_2 \frac{\partial N_p}{\partial x_2}
            +
            \beta_2 \frac{\partial N_q}{\partial x_2}
            \beta_1 \frac{\partial N_p}{\partial x_1}
            +
            \beta_1 \frac{\partial N_q}{\partial x_2}
            \beta_2 \frac{\partial N_p}{\partial x_2}
        \right)
    \ d\boldsymbol{x} \right]
    \\
    &= \tau_k |e_k| \left(
    \beta_1 \beta_1 \begin{bmatrix}
        b_1 b_1 & b_1 b_2 & b_1 b_3 \\
        b_2 b_1 & b_2 b_2 & b_2 b_3 \\
        b_3 b_1 & b_3 b_2 & b_3 b_3
    \end{bmatrix}
    +
    \beta_1 \beta_2 \begin{bmatrix}
        b_1 c_1 & b_1 c_2 & b_1 c_3 \\
        b_2 c_1 & b_2 c_2 & b_2 c_3 \\
        b_3 c_1 & b_3 c_2 & b_3 c_3
    \end{bmatrix}
    +
    \beta_2 \beta_1 \begin{bmatrix}
        c_1 b_1 & c_1 b_2 & c_1 b_3 \\
        c_2 b_1 & c_2 b_2 & c_2 b_3 \\
        c_3 b_1 & c_3 b_2 & c_3 b_3
    \end{bmatrix}
    +
    \beta_2 \beta_2 \begin{bmatrix}
        c_1 c_1 & c_1 c_2 & c_1 c_3 \\
        c_2 c_1 & c_2 c_2 & c_2 c_3 \\
        c_3 c_1 & c_3 c_2 & c_3 c_3
    \end{bmatrix}
    \right)
    \\
    \boldsymbol{D}_k
    &= \left[ \int_{e_k}
        \alpha \frac{\partial N_q}{\partial x_j} \frac{\partial N_p}{\partial x_j}
    \ d\boldsymbol{x} \right]
    \\
    &= \left[ \int_{e_k}
        \alpha \frac{\partial N_q}{\partial x_1} \frac{\partial N_p}{\partial x_1}
        +
        \alpha \frac{\partial N_q}{\partial x_2} \frac{\partial N_p}{\partial x_2}
    \ d\boldsymbol{x} \right]
    \\
    &= \alpha |e_k| \begin{bmatrix}
        b_1 b_1 + c_1 c_1 & b_1 b_2 + c_1 c_2 & b_1 b_3 + c_1 c_3 \\
                          & b_2 b_2 + c_2 c_2 & b_2 b_3 + c_2 c_3 \\
        \text{sym.}       &                   & b_3 b_3 + c_3 c_3
    \end{bmatrix}
    \\
    \boldsymbol{S}_k
    &= \left[ \int_{e_k}
        f N_p
    \ d\boldsymbol{x} \right]
    \\
    &= \frac{f |e_k|}{3} \begin{bmatrix}
        1 \\ 1\\ 1
    \end{bmatrix}
    \\
    \check{\boldsymbol{S}}_k
    &= \left[ \tau_k \int_{e_k}
        f
    \left(
        \beta_i \frac{\partial N_p}{\partial x_i}
    \right)
    \ d\boldsymbol{x} \right]
    \\
    &= \tau_k |e_k| \left(
        \beta_1 f \begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix}
        +
        \beta_2 f \begin{bmatrix} c_1 \\ c_2 \\ c_3 \end{bmatrix}
    \right)
    \\
    \boldsymbol{Q}_k
    &= \left[ \int_{\Gamma_N}
        g_N N_p
    \ d\boldsymbol{s} \right]
    \\
    &= \frac{g_N h_k}{2} \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix}
\end{align}

である．各係数行列の計算には，面積座標の便利公式

\begin{align}
    \int_{e_k}
        L_1^p L_2^q L_3^r
    \ d\boldsymbol{x}
    &= \frac{2 \; |e_k| \; p! \; q! \; r!}{(2 + p + q + r)!}
\end{align}

を用いた．但し，$(L_1, L_2, L_3)$は面積座標であり，形状関数$(N_1, N_2, N_3)$に等しい．

なお，多次元における安定化パラメータ$\tau_k$は，

\begin{align}
    \tau_k
    &= \frac{h_k}{2 |\boldsymbol{\beta}|} \xi (\text{Pe}_k)
    \\
    \xi (\text{Pe}_k)
    &= \coth (\text{Pe}_k) - \frac{1}{\text{Pe}_k}
    \\
    \text{Pe}_k
    &= \frac{|\boldsymbol{\beta}| h_k}{2 \alpha}
\end{align}

という，1次元の自然な拡張（cf. 続・有限要素法による流れのシミュレーション）と，

\begin{align}
    \tau_k
    &= \left(
        \left( \frac{2 |\boldsymbol{\beta}|}{h_k} \right)^2
        +
        \left( \frac{4 \alpha}{h_k^2} \right)^2
    \right)^{-1/2}
    \\
\end{align}

が知られている（cf. 第3版 有限要素法による流れのシミュレーション）．前者をSUPG-v1，後者をSUPG-v2と呼ぶこととする．

Shock-Capturing

SUPGは，衝撃面付近ではovershoot/undershootを起こすことが知られている．これに対して，

\begin{align}
    \sum_k \delta_k \int_{e_k}
        \frac{\partial u}{\partial x_j} \frac{\partial v}{\partial x_j}
    \ d\boldsymbol{x}
\end{align}

という，衝撃補足項を付け加えることで回避できる．但し，

\begin{align}
    \delta_k
    &= \frac{|\boldsymbol{\beta}| h_k}{2} \zeta (\text{Pe}_k)
    \\
    \zeta (\text{Pe}_k)
    &= \begin{cases}
        \text{Pe}_k / 3 & (\text{Pe}_k \le 3) \\
        1               & (\text{Pe}_k \gt 3)
    \end{cases}
\end{align}

は安定化パラメータ．

最後に，斜め移流の問題に対して，Galerkin，SUPG-v1，SUPG-v2を比較する．

$\text{Galerkin}$	$\text{SUPG-v1}$	$\text{SUPG-v2}$

SUPGの安定化パラメータは，この問題では殆ど同じ数字になっていた．なお，教科書と違う見た目なのは，恐らくメッシュが少し異なるからであり，要素の並び方を僅かに変えると，教科書と同じような結果が得られる（Galerkinに激しい振動）．

$\text{Galerkin}$	$\text{SUPG-v1}$	$\text{SUPG-v2}$

続いて，SUPG-v1，SUPG-v2にshock-captureを追加する．

$\text{SUPG-v1}$	$\text{SUPG-v1 + SC}$	$\text{SUPG-v2}$	$\text{SUPG-v2 + SC}$

あまりにも拡散が強い．何か間違えているに違いない．．．が，未だ良く分からない．

終わりに

なんて難しいのだろう．

Appendix

Gauss-Greenの式

Gaussの発散の定理

\begin{align}
    \int_{\Omega}
        \partial_j A_j
    \ d\boldsymbol{x}
    &= \int_{\Gamma}
        A_j n_j
    \ d\boldsymbol{s}
\end{align}

において，ヴェクトル$\boldsymbol{A}$を，スカラー$u, \ v$を以て$\boldsymbol{A} = v \nabla u$と表すと，

\begin{align}
    \partial_j A_j
    &= \partial_j (v \partial_j u)
    \\
    &= \partial_j v \partial_j u
        + v \partial_j \partial_j u
    \\
    A_j n_j
    &= (v \partial_j u) n_j
    \\
    &= v \partial_j u n_j
\end{align}

であるから，

\begin{align}
    \int_{\Omega}
        \left( 
            \partial_j v \partial_j u
            + v \partial_j \partial_j u
        \right)
    \ d\boldsymbol{x}
    &= \int_{\Gamma}
        v \partial_j u n_j
    \ d\boldsymbol{s}
\end{align}

が成り立つ．

弱形式化

弱形式化の手続きの中では，

\begin{align}
    - \int_{\Omega}
        v \partial_j \partial_j u
    \ d\boldsymbol{x}
    &= \int_{\Omega}
        \partial_j v \partial_j u
    \ d\boldsymbol{x}
    - \int_{\Gamma}
        v \partial_j u n_j
    \ d\boldsymbol{s}
\end{align}

を使った．

安定化項

安定化項を取り出すときには，

\begin{align}
    \int_{\Omega}
        \partial_j w \partial_j u
    \ d\boldsymbol{x}
    &= - \int_{\Omega}
        w \partial_j \partial_j u
    \ d\boldsymbol{x}
    + \int_{\Gamma}
        w \partial_j u n_j
    \ d\boldsymbol{s}
\end{align}

に対して，$w = \partial_i v$とすることで，

\begin{align}
    \int_{\Omega}
        \partial_j (\partial_i v) \partial_j u
    \ d\boldsymbol{x}
    &= - \int_{\Omega}
        (\partial_i v) \partial_j \partial_j u
    \ d\boldsymbol{x}
    + \int_{\Gamma}
        (\partial_i v) \partial_j u n_j
    \ d\boldsymbol{s}
\end{align}

を使った．