概要
- 境界埋め込み法の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度。少し剥離してます。
