概要

境界埋め込み法の1種であるVolume Penalization法（Angot P., Bruneau C.-H., Fabrie P., Numerische Mathmatik, 81(1999) pp. 497-520)をOpenFOAMの非圧縮ソルバのpisoFoamに実装。
NACA0012の周りの流れの計算をテストで実施。

Volume Penalization法（VP法）

非圧縮の場合だけ説明します。VP法では物体を浸透率のある多孔質媒体で模擬します。下記のNavier-Stokes方程式の最後の項の外力項（penalization項）が物体による効果を表します。χはマスク関数（物体で1，流体で0），ηが浸透率，usが物体境界上の速度です。すべりなしではus=0です。（uは速度ベクトルだと思ってください）

\frac{\partial u}{\partial t} + (u\cdot \nabla)u = - \nabla p + \nu \nabla^2 u + \frac{\chi}{\eta}(u-u_s)

OpenFOAMへの実装

ベースはpisoFoamを使いました。OpenFOAMは9です。

マスク関数χは新しい物理量を作って，場の量として表します。
浸透率ηを物性値として，設定するようにします。

変更点だけ書きます。ソルバ名はMake/filesで変更してください。

createFields.Hにχとηの定義を追加。

Info<< "Reading field chi\n" << endl;
volScalarField chi
(
        IOobject
        (
            "chi",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
);

 IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
        )
    );

    dimensionedScalar eta
    (
        transportProperties.lookup("eta")
    );

UEqn.Hを変更

fvVectorMatrix UEqn
(
    fvm::ddt(U) + fvm::div(phi, U)
  + MRF.DDt(U)
  + turbulence->divDevSigma(U)
  + fvm::Sp(chi/eta, U)
 ==
    fvModels.source(U)
);

テストケース（NACA0012周りの流れ）

テストケースとして，NACA0012翼周りの流れを計算します。NACA0012の形状は代数的に与えられるらしいです。ここを参考にさせてもらいました。https://cattech-lab.com/science-tools/naca0012-validation/

流れの条件も参考にさせてもらいました。流入速度は44m/s，動粘性係数は$1\times10^{-5}$m2/s，コード長1mにします。

まず，χ (chi)のファイルを0ディレクトリに作ります。（注：今回周期境界は使いませんが，以前使った名残で残ってます）

0/chi

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      Q;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    ".*"
    {
        type            cyclic;
    }


}

// ************************************************************************* //

U, pはこんな感じ。

0/U

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    inlet
    {
        type            fixedValue;
        value           uniform (44 0 0);
    }

    outlet
    {
        type            zeroGradient;
    }

    ".*"
    {
        type            cyclic;

    }

//    front
//    {
//        type            empty;
//    }
}

// ************************************************************************* //

0/p

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 2 -2 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }

    ".*"
    {
        type           cyclic;
        
    }


//    front
//    {
//        type            empty;
//    }
}

// ************************************************************************* //

transportPropertiesを修正します。etaの値は大きいと物体に流体が結構はいるので，かなり小さめの値にします。

constant/transportProperties

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  9
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

transportModel  Newtonian;

nu              [0 2 -1 0 0 0 0] 1e-5; // Re_blockage = 100

eta           eta  [0 0  1 0 0 0 0] 1e-8;

// ************************************************************************* //

blockMeshDictはこんなかんじ。翼はchiの分布で与えるのでただの四角い領域です。

system/blockMeshDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
    (-1.5 -1.5 0)   
    (3 -1.5 0)
    (3 1.5 0)
    (-1.5 1.5 0)
    (-1.5 -1.5 3)   
    (3 -1.5 3)
    (3 1.5 3)
    (-1.5 1.5 3)

);

blocks
(
    hex (0 1 2 3 4 5 6 7) (900 600 1)
    simpleGrading (1 1 1)

);

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 3 7 4)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (1 2 6 5)
            
             (0 1 5 4)
             (3 2 6 7)

        );
    }
    bottom
    {
        type empty;
        
        faces
        (
            (0 1 2 3)
            (4 5 6 7)
            
        );
    }

);

mergePatchPairs
(
);

// ************************************************************************* //

NACA0012の形状をχの分布に与える

NACA0012の形状は

y=\pm 0.6(0.2969\sqrt{x}-0.126x-0.3516x^2+0.2843x^3-0.1015x^4)

らしいです。この境界で囲まれた領域のchiの値を1にします。ここではswak4Foamの機能の1つであるfunkySetFieldsを使います。このままだと迎え角が0度しかできないので，座標の回転を計算して，迎え角を与えます。

元の座標$(x,y)$を角度$\theta$回転して$(x',y')$とすると，

x' = (\cos \theta) x - (\sin \theta) y

y'= (\sin \theta) x + (\cos \theta) y

これつかって，迎え角aのときのNACA0012の形状を表す関数は以下（上側だけかきます）

(\sin a) x + (\cos a)y = 0.6(0.2969\sqrt{(\cos a)x - (\sin a)y}-0.126((\cos a)x - (\sin a)y)-0.3516((\cos a)x - (\sin a)y)^2+0.2843((\cos a)x - (\sin a)y)^3 - 0.1015((\cos a)x - (\sin a)y) ^4

迎え角10度にしたときのfunkySetFieldsDictが以下。ルートのなかが負にならないようにしてます。翼がある領域だけ設定するための条件設定してます。

system/funkySetFieldsDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/

FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      funkySetFieldsDict;
}

expressions
(
    naca0012
    {
        field chi;
        expression "1.0";
        
        condition "(0.984*pos().x-0.173*pos().y) >= 0.0 && (0.984*pos().x-0.173*pos().y) <= 1.0 && -0.6*( 0.2969*sqrt(mag(0.984*pos().x-0.173*pos().y)) - 0.126*(0.984*pos().x-0.173*pos().y) -0.3516*pow(0.984*pos().x-0.173*pos().y, 2) +0.2843*pow(0.984*pos().x-0.173*pos().y, 3) -0.1015*pow(0.984*pos().x-0.173*pos().y, 4)) <=0.173*pos().x+ 0.984*pos().y && 0.6*( 0.2969*sqrt(mag(0.984*pos().x-0.173*pos().y)) - 0.126*(0.984*pos().x-0.173*pos().y) -0.3516*pow(0.984*pos().x-0.173*pos().y, 2) +0.2843*pow(0.984*pos().x-0.173*pos().y, 3) -0.1015*pow(0.984*pos().x-0.173*pos().y, 4)) >=0.173*pos().x+ 0.984*pos().y";
        keepPatches 1;
    }
);

これつかってfunkySetFieldsでchiの値に分布を与えます。

速度分布

迎え角10度の結果です。流れが剥離してます。拡大するとわかりますが，表面のでこぼこが結構見えます。流れが剥離しがちなので，Cd値とかはあってないと思います。もっとメッシュ細かくすると改善すると思います。

迎え角4度。少し剥離してます。