多孔質体の解析

2012年11月4日

はじめに

多孔質体 (porous media) を使った解析について。

使用バージョン

OpenFOAM 2.1.1

ファイル

多孔質体の解析

多孔質体の解析は、porousSimpleFoam で行える。多孔質体の設定は constant/porousZone で行う。設定できるタイプが 2 つある。

  • Darcy (Darcy-Forchheimer 則)
  • powerLaw (べき乗則)

多孔質体は、運動方程式のソース項に抵抗を与えることで表現する。Darcy-Forchheimer 則は次式で表される。

Si = -(μ D + 1/2 ρ |u| F) ui

ここで Si は運動方程式のソース項の i 成分、μ は粘性係数、ρ は密度、u は速度ベクトル、ui は u の i 成分である。D と F は係数で、ユーザーが指定することになる。

べき乗則の場合は次式になる。

Si = -ρ C0 |u|^((C1 - 1)/2)

C0, C1 が指定すべき係数である。

Darcy 則で指定する場合、constant/porousZone は次のようにする。

1
(
    porosity
    {
        coordinateSystem
        {
            e1  (1 0 0);
            e2  (0 1 0);
        }

        Darcy
        {
            d   d [0 -2 0 0 0 0 0] (2.5e6 -1000 -1000);
            f   f [0 -1 0 0 0 0 0] (0 0 0);
        }
    }
)

ここで "porosity" はセルゾーンの名前である。coordinateSystem は、多孔質体の方向を決めるための座標系の定義である。

べき乗則では次のようにする。

1
(
    porosity
    {
        powerLaw
        {
            C0  25000;
            C1  1;
        }
    }
)

多孔質体の解析の方法として、陽的なものと陰的なものがある。system/fvSolution の "SIMPLE" において "nUCorrectors" で反復回数を指定すると、陰的な計算が行われる。

SIMPLE
{
    nUCorrectors 2;
    ...

矩形領域 (0.1 m x 0.1 m、厚み 0.01 の 2D 領域) に左から右に向かって 1 m/s の流れを作る。真ん中に 0.04 m の幅の多孔質体の仕切りを作り、1000 Pa の圧力損失を与えることにする。動粘性係数 ν は 0.01 m2/s とする。

ここでは多孔質体のセルゾーンを setSet で作ることにする。次のような setBatch ファイルを作る。

cellSet porosity new boxToCell (0.03 0 0) (0.07 0.1 0.01)
cellZoneSet porosity new setToCellZone porosity

setSet の実行。

$ setSet -batch setBatch

Darcy 則は本来次式で表される。

Δp/L = -μ D u

Δp は圧力差、L は幅である。ここでは Δp = 1000、L = 0.04、u = 1、μ は ν と読み替えて μ = 0.01 として、D = 2.5e6 である。上記のように D と F を設定すると、下図のような圧力分布が得られる。

simpleFoam への多孔質体解析機能の追加

simpleFoam へ多孔質体解析機能を追加するための最低限の作業は、以下の通りである。

simpleFoam.C で porousZones.H をインクルード。

#include "porousZones.H"

createFields.H で porousZones のオブジェクトを定義。

porousZones pZones(mesh);

UEqn.H で方程式に抵抗を追加。

    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
    );

    UEqn().relax();
    
    pZones.addResistance(UEqn());

    solve(UEqn() == -fvc::grad(p));

Make/options で meshTools のインクルードパスとライブラリの指定をして wmake。