OpenFOAMにおける多孔質体の設定（porous model / fvOption)

#はじめに
本記事は多孔質体を含めた領域の流体解析を行うことを目的としています。

[こちら](http://penguinitis.g1.xrea.com/study/OpenFOAM/porous/porous.html)に多孔質体ソルバー`porousSimpleFoam`の詳細がまとめられていますが、ここでは、pisoFoam等の基本ソルバーに多孔質体モデルを導入するやり方(fvOptionを利用）をまとめたいと思います。

# 環境
なお、本記事はFoundation版`OpenFOAM v8`、`pisoFoam/laminar/porousBlockage`を使用しています。

# Porous model
多孔質モデルは、基本、運動方程式のソース項 $S\_m$ に組み込まれます。

```math
\frac{\partial \rho U}{\partial t} + (U \cdot \nabla) U = \nabla \sigma + S_m
```


## Darcy-Forchheiman equation
Darcy-Forchheiman 多孔質モデルでは、ソース項は以下のように表されます。

```math
S_m = - \left( \mu D + \frac{\rho |U|}{2} f  \right) U 
```

ここで、D は Darcy coefficient [1/m$^2$]、f は Forchheimer coefficient [1/m]を表しています。

## Power law
べき乗則モデルでは、

```math
S_m = - \rho C_0 |U|^{(C_1 -1)/2}
```

のようにべき乗則で流動抵抗が表されます。ここで、$C_0$ は model linear coefficient、$C_1$ は model exponent coefficientであり、ユーザーが設定する必要があります。

## その他モデル
その他にも、以下のようなモデルが用意されています。

![Screen Shot 2021-01-05 at 23.53.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/227359/f253382a-2636-d47f-c3ab-b1e1555a83e5.png)

# 設定方法
ここでは、先に述べたとおり`pisoFoam/laminar/porousBlockage`チュートリアルを見ていきます。

fvOptionを使った多孔質体モデルの実装では、
- topoSetによる多孔質領域の設定
- 多孔質体の物性値の設定
の大きく２つを行う必要があります。

## topoSet

```c++:system/topoSetDict
    // porousBlockage
    {
        name    porousBlockageCellSet;
        type    cellSet;
        action  new;
        source  boxToCell;
        sourceInfo
        {
            box (-0.5 -0.5 -1) (0.5 0.5 1);
        }
    }
    {
        name    porousBlockage;
        type    cellZoneSet;
        action  new;
        source  setToCellZone;
        sourceInfo
        {
            set porousBlockageCellSet;
        }
    }
```

## fvOption

```c++:constant/fvOption
    type            explicitPorositySource;

    explicitPorositySourceCoeffs
    {
        selectionMode   cellZone;
        cellZone        porousBlockage;

        type            DarcyForchheimer;

        // D 100;  // Very little blockage
        // D 200;  // Some blockage but steady flow
        // D 500;  // Slight waviness in the far wake
        D 1000; // Fully shedding behavior

        d   ($D $D $D);
        f   (0 0 0);

        coordinateSystem
        {
            type    cartesian;
            origin  (0 0 0);
            coordinateRotation
            {
                type    axesRotation;
                e1  (1 0 0);
                e2  (0 1 0);
            }
        }
    }
```

このチュートリアルでは多孔質体が等方的な圧損を有している場合を表しています。

ハニカムなどの一方向のみ流れやすく、多方向には圧損が生じるような場合は、$d$、$f$の値を以下のように方向性を持たせて設定します。

```c++
    DarcyForchheimerCoeffs
    {
        d   d [0 -2 0 0 0 0 0] (7.65e7 3e10 3e10);
        f   f [0 -1 0 0 0 0 0] (255 1e5 1e5);

        coordinateSystem
        {
            type    cartesian;  // global co-ordinate system (redundant)
            origin  (0 0 0);    // redundant in this case
            coordinateRotation
            {
                type    axesRotation; // local Cartesian co-ordinates
                e1      (0.998 0.061 0);
                e2      (0.016 -1.000 0);
            }
        }
    }
```

なお、`coordinateSystem`は、多孔質体のlocal座標を設定するためのようです（どうしてlocal座標を設定するのか、その意図は指し図りきれていないです・・・）。local座標は３つのベクトル（$\bf{e_1}$, $\bf{e_2}$, $\bf{e_3}$）を使って設定しますが、入力ではうち２つを設定する必要があります。特にlocal座標を考える必要がない場合は、

```c++
coordinateRotation
{
    type    axesRotation; // local Cartesian co-ordinates
    e1      (1 0 0);
    e2      (0 1 0);
}
```
のようにx, yベクトルを設定します。なお、$\bf{e_3} = \bf{e_1} \times \bf{e_2}$ で自動的に規定されます。

- porousBlockageの解析結果例
![Screen Shot 2021-01-06 at 0.29.08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/227359/be6a6559-bc52-31cc-c559-1da40e49daca.png)


#さいごに
以上、駆け足でfvOptionを使った多孔質体の設定方法を見てきました。以前は、自らでソース項をUEqn.Hなどに追加する必要がありましたが、fvOptionを使えばそのようなソースからのビルドをする必要がありません。
まだ、説明を端折っているところもありますので、順次、加筆していきたいと思います。