Go to Qiita Advent Calendar Top
LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@ShotaDeguchi(Shota DEGUCHI)

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

Last updated at Posted at 2025-09-30

初めに

時間を有限差分，空間をP1有限要素で離散化する．このとき，GLSはSUPGと同一．

Methods

非定常な移流拡散方程式

以下である．

\begin{alignat}{2}
    \partial_t u
    + \boldsymbol{\beta} \cdot \nabla u
    - \nabla \cdot \alpha \nabla u
    &= f
    \quad &&\text{in} \quad \Omega \times \mathcal{I}
    \\
    u
    &= g_D
    \quad &&\text{on} \quad \Gamma_D \times \mathcal{I}
    \\
    \boldsymbol{n} \cdot \alpha \nabla u
    &= g_N
    \quad &&\text{on} \quad \Gamma_N \times \mathcal{I}
    \\
    u
    &= u^{0}
    \quad &&\text{in} \quad \bar{\Omega} \times \lbrace 0 \rbrace
\end{alignat}

空間離散化

SUPG法に則って空間離散化を施し，

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

$u = N_{\beta} u_{\beta}, \ v = N_{\alpha} v_{\alpha}$とP1補間し，

\begin{align}
    \int_{e_k}
        v \frac{\partial u}{\partial t}
    \ d\boldsymbol{x}
    &= v_{\alpha} \left[
        \int_{e_k}
            N_{\alpha} 
            N_{\beta}
        \ d\boldsymbol{x}
    \right] \frac{\partial}{\partial t} u_{\beta}
    = v_{\alpha} \boldsymbol{M}_k \frac{\partial}{\partial t} u_{\beta}
    \\
    \int_{e_k}
        v \beta_j \frac{\partial u}{\partial x_j}
    \ d\boldsymbol{x}
    &= v_{\alpha} \left[
        \int_{e_k}
            \beta_j N_{\alpha} \frac{\partial N_{\beta}}{\partial x_j}
        \ d\boldsymbol{x}
    \right] u_{\beta}
    = v_{\alpha} \boldsymbol{A}_k u_{\beta}
    \\
    \int_{e_k}
        \alpha \frac{\partial v}{\partial x_j} \frac{\partial u}{\partial x_j}
    \ d\boldsymbol{x}
    &= v_{\alpha} \left[
        \int_{e_k}
            \alpha \frac{\partial N_{\alpha}}{\partial x_j} \frac{\partial N_{\beta}}{\partial x_j}
        \ d\boldsymbol{x}
    \right] u_{\beta}
    = v_{\alpha} \boldsymbol{D}_k u_{\beta}
    \\
    \int_{\Omega}
        v f
    \ d\boldsymbol{x}
    &= v_{\alpha} \left[
        \int_{e_k}
            N_{\alpha} f
        \ d\boldsymbol{x}
    \right]
    = v_{\alpha} \boldsymbol{S}_k
    \\
    \int_{e_k}
    \left(
        \beta_i \frac{\partial v}{\partial x_i}
    \right)
    \left(
        \frac{\partial u}{\partial t}
    \right)
    \ d\boldsymbol{x}
    &= v_{\alpha} \left[
        \tau_k
        \int_{e_k}
        \left(
            \beta_i \frac{\partial N_{\alpha}}{\partial x_i}
        \right)
        N_{\beta}
        \ d\boldsymbol{x}
    \right] \frac{\partial}{\partial t} u_{\beta}
    = v_{\alpha} \check{\boldsymbol{M}}_k \frac{\partial}{\partial t} u_{\beta}
    \\
    \int_{e_k}
    \left(
        \beta_i \frac{\partial v}{\partial x_i}
    \right)
    \left(
        \beta_j \frac{\partial u}{\partial x_j}
    \right)
    \ d\boldsymbol{x}
    &= v_{\alpha} \left[
        \tau_k
        \int_{e_k}
        \left(
            \beta_i \frac{\partial N_{\alpha}}{\partial x_i}
        \right)
        \left(
            \beta_j \frac{\partial N_{\beta}}{\partial x_j}
        \right)
        \ d\boldsymbol{x}
    \right] u_{\beta}
    = v_{\alpha} \check{\boldsymbol{A}}_k u_{\beta}
    \\
    \int_{e_k}
    \left(
        \beta_i \frac{\partial v}{\partial x_i}
    \right)
    f
    \ d\boldsymbol{x}
    &= v_{\alpha} \left[
        \tau_k
        \int_{e_k}
        \left(
            \beta_i \frac{\partial N_{\alpha}}{\partial x_i}
        \right)
        f
        \ d\boldsymbol{x}
    \right]
    = v_{\alpha} \check{\boldsymbol{S}}_k
    \\
    \int_{\Gamma_N}
        v g_N
    \ d\boldsymbol{s}
    &= v_{\alpha} \left[
        \int_{\Gamma_N}
            N_{\alpha} g_N
        \ d\boldsymbol{s}
    \right]
    = v_{\alpha} \boldsymbol{Q}_k
\end{align}

$v$の任意性から，

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

時間微分に関する行列（質量行列と呼ぶ）は，面積座標の便利公式から，

\begin{align}
    \boldsymbol{M}_k
    &= \left[
        \int_{e_k}
            N_{\alpha} 
            N_{\beta}
        \ d\boldsymbol{x}
    \right]
    = \frac{|e_k|}{12} \begin{bmatrix}
        2 & 1 & 1 \\
        1 & 2 & 1 \\
        1 & 1 & 2
    \end{bmatrix}
    \\
    \check{\boldsymbol{M}}_k
    &= \left[
        \tau_k
        \int_{e_k}
        \left(
            \beta_i \frac{\partial N_{\alpha}}{\partial x_i}
        \right)
        N_{\beta}
        \ d\boldsymbol{x}
    \right]
    = \frac{\tau_k |e_k|}{3} \left(
        \beta_1
        \begin{bmatrix}
            b_1 & b_1 & b_1 \\
            b_2 & b_2 & b_2 \\
            b_3 & b_3 & b_3
        \end{bmatrix}
        +
        \beta_2
        \begin{bmatrix}
            c_1 & c_1 & c_1 \\
            c_2 & c_2 & c_2 \\
            c_3 & c_3 & c_3
        \end{bmatrix}
    \right)
\end{align}

であるから，要素を集めて，

\begin{align}
    ( \boldsymbol{M} + \check{\boldsymbol{M}} ) \partial_t \boldsymbol{u}
    + ( \boldsymbol{A} + \check{\boldsymbol{A}} + \boldsymbol{D} ) \boldsymbol{u}
    &= \boldsymbol{S} + \check{\boldsymbol{S}} + \boldsymbol{Q}
\end{align}

を解く．

時間離散化

$\theta$法を考える（$0 \le \theta \le 1$）．但し，ソースと境界条件は時間的に変動しないとする．

\begin{align}
    ( \boldsymbol{M} + \check{\boldsymbol{M}} ) \frac{\boldsymbol{u}^{n+1} - \boldsymbol{u}^{n}}{\Delta t}
    + ( \boldsymbol{A} + \check{\boldsymbol{A}} + \boldsymbol{D} ) ( \theta \boldsymbol{u}^{n+1} + ( 1 - \theta) \boldsymbol{u}^{n} )
    &= \boldsymbol{S} + \check{\boldsymbol{S}} + \boldsymbol{Q}
    \\
    \left(
    \frac{1}{\Delta t} (\boldsymbol{M} + \check{\boldsymbol{M}})
    + \theta ( \boldsymbol{A} + \check{\boldsymbol{A}} + \boldsymbol{D} )
    \right)
    \boldsymbol{u}^{n+1}
    &= \boldsymbol{S} + \check{\boldsymbol{S}} + \boldsymbol{Q}
    + \left(
        \frac{1}{\Delta t} (\boldsymbol{M} + \check{\boldsymbol{M}})
        - ( 1 - \theta) ( \boldsymbol{A} + \check{\boldsymbol{A}} + \boldsymbol{D} )
    \right) \boldsymbol{u}^{n}
\end{align}

このとき，

  • $\theta = 0$: Explicit Euler
  • $\theta = 1$: Implicit Euler
  • $\theta = 1/2$: Crank-Nicolson

である．

$\theta = 0$を選ぶと，lumpingでもしないと$\theta = 0$にした意味を説明し辛い．が，いまのところlumpingは受け入れられない．

$\theta = 1/2$は時間方向に2次精度を与えるので嬉しいが，暫く，簡単のため，$\theta = 1$を選択しておく．

安定化

SUPGの安定化パラメータには，定常なときに使った

\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}

と，非定常な問題に拡張された

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

の，2つを試す（cf. 第3版 有限要素法による流れのシミュレーション）．それぞれ，v1, v2と呼ぶこととする．

Results

SUPGによる定常解を初期解に採用し，移流方向を$\pi / 2$だけ回転させる．いまのところ，shock-capturingには懐疑的であるので，入れていない．

$\text{Galerkin}$ $\text{SUPG-v1}$ $\text{SUPG-v2}$
galerkin.gif supg_v1.gif supg_v2.gif

順当な改善傾向．

終わりに

なんて面白いのだろう．

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?