【OpenFOAM】スロッシングの実験比較（設定ファイル解説あり）

こんにちは（@t_kun_kamakiri）

本記事ではOpenFOAMを用いたスロッシングタンクの流体解析を設定の解説を行います。

本記事の内容

OpenFOAMでスロッシングタンクの解析設定を解説
スロッシングタンクの実験データとの比較を行う

スロッシングタンクは、タンク内に入れた液体が周期的な振動によりうねる現象のことです。構造物には構造物特有の固有振動を持っているため、その周期と重なったときに共振を起こし小さな力でも問題になることがあります。

こんな方を対象

OpenFOAMでダイナミックメッシュをやってみたい
流体解析を実験データと比較したい

ベンチマークテストの試験データとOpenFOAMとの比較も行います。

OpenFOAMで再現するのはこちらの実験です。
タンクが振動して水がバシャバシャ暴れていますね。

CAEと実機との比較検討は本記事では十分に行うことができませんが、解析設定から実験データとの比較まで一連の流れを体験することができます。

解析対象

まずは解析対象を確認します。
実験の元データはこちらになります。

自身が解析するフォルダにとりあえず保存しておきましょう。
フォルダの中に試験の解説が書かれているpdfがあります。

試験の寸法はこちらにきっちり書かれています。
Sensor1~6は実験での圧力センサーの位置です。
OpenFOAMでも同じ位置にprobeを設けて圧力の時刻歴変化を出力するようにします。
試験ではz方向の厚みを31mm,62mm,124mmの3パターン行っていますが、データとして62mmのみを保存しているとのことです。

振動周期は以下のように書かれています。

$T_{1}$が振動の周期[s]です。
OpenFOAMでこれ通りに与えると、そもそも試験での振動と上記の数式による振動の周期が微妙にずれていました。
ですので、ここでは実験で得られた試験データの振動周期をOpenFOAMにも与えるようにします。試験条件が違うとなると実機とCAEのコリレーションの意味がなくなるので。

H[mm]の高さが2パターンありますが、これは液体の初期水位です。
今回示すのはH=93mmですので、93mだけ見ておけば良いです。
また、液体はwater（水）です。

物性値は以下のようになっています。
Water(水)による実験データと比較するのでwaterだけを見ておけば良いです。

以上で解析対象の説明が終わりました。
以下でOpenFOAMの解析設定の解説を行います。

チュートリアルから解析ケースをコピー

OpenFOAMにスロッシング用のチュートリアルがあるのでそれをコピーします。

cp -r $FOAM_TUTORIALS/multiphase/interFoam/laminar/sloshingTank3D6DoF 001_run

1	cp -r $FOAM_TUTORIALS/multiphase/interFoam/laminar/sloshingTank3D6DoF 001_run

※「$FOAM_TUTORIALS＝/usr/lib/openfoam/openfoam2212/tutorials」

ソルバはinterFoamを使います。
気体と液体を体積分率$\alpha$で区別するVOF法による解析ソルバです。

先ほど入手した試験データは001_runと同じ階層に保存しておきます。

001_run
expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10

9449af_05f4aa0eaa1d4164a1024d07c562592aは元々のフォルダ名ですので、この辺は好みで名前やフォルダ構成を変えてください。

メッシュ作成

形状は直方体なのでblockMeshで簡単に設定することができます。
形状作成とメッシュ作成、境界面の名前はblockMeshで行います。

system/blockMeshDict

scale   0.001;

xMin -450;
xMax  450;
yMin    0;
yMax  508;
zMin    0;
zMax   62;

dx 10;
dy 10;
dz 5;

xl #calc "$xMax-($xMin)";
yl #calc "$yMax-($yMin)";
zl #calc "$zMax-($zMin)";

nx #calc "$xl/$dx";
ny #calc "$yl/$dy";
nz #calc "$zl/$dz";

vertices
(
    ($xMin $yMin $zMin)
    ($xMax $yMin $zMin)
    ($xMax $yMax $zMin)
    ($xMin $yMax $zMin)
    ($xMin $yMin $zMax)
    ($xMax $yMin $zMax)
    ($xMax $yMax $zMax)
    ($xMin $yMax $zMax)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    walls
    {
        type wall;
        faces
        (
            (3 7 6 2)
            (0 4 7 3)
            (2 6 5 1)
            (1 5 4 0)
            (0 3 2 1)
            (4 5 6 7)
        );
    }
);

mergePatchPairs
(
);

scale 0.001;

xMin -450;

xMax 450;

yMin 0;

yMax 508;

zMin 0;

zMax 62;

dx 10;

dy 10;

dz 5;

xl #calc "$xMax-($xMin)";

yl #calc "$yMax-($yMin)";

zl #calc "$zMax-($zMin)";

nx #calc "$xl/$dx";

ny #calc "$yl/$dy";

nz #calc "$zl/$dz";

vertices

(

($xMin $yMin $zMin)

($xMax $yMin $zMin)

($xMax $yMax $zMin)

($xMin $yMax $zMin)

($xMin $yMin $zMax)

($xMax $yMin $zMax)

($xMax $yMax $zMax)

($xMin $yMax $zMax)

);

blocks

(

hex (0 1 2 3 4 5 6 7) ($nx $ny $nz) simpleGrading (1 1 1)

);

edges

(

);

boundary

(

walls

{

type wall;

faces

(

(3 7 6 2)

(0 4 7 3)

(2 6 5 1)

(1 5 4 0)

(0 3 2 1)

(4 5 6 7)

);

}

);

mergePatchPairs

(

);

設定が終わりましたら、blockMeshを実行します。

blockMesh

blockMesh

結果をParaViewで確認しましょう。

形状はいたってシンプルですね。

境界条件の設定

乱流モデルを使わない解析ですので、設定する物理量はU(流速)、p_rgh(圧力)、体積分率alpha.waterです。
圧力は計算の便宜上$p-\rho \boldsymbol{g}\cdot h$という値で行っている点に注意です。

いつでも初期化できるように0.origフォルダ内のファイルを編集します。

0.orig/U

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

internalField   uniform (0 0 0);

boundaryField
{
    walls
    {
        type            movingWallVelocity;
        value           uniform (0 0 0);
    }
}

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

internalField uniform (0 0 0);

boundaryField

{

walls

{

type movingWallVelocity;

value uniform (0 0 0);

}

壁条件ではありますが振動をしているのでtypeをmovingWallVelocityとしています。

※movingWallVelocity：パッチでのフラックスが0となるよう、パッチでの値の法線方向を置き換える。

0.orig/p_rgh

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

internalField   uniform 0;

boundaryField
{
    walls
    {
        type            fixedFluxPressure;
    }
}

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

internalField uniform 0;

boundaryField

{

walls

{

type fixedFluxPressure;

}

fixedFluxPressureは境界のフラックスが速度の境界条件にマッチするように圧力勾配を調整するものです。

0.orig/alpha.water

dimensions      [0 0 0 0 0 0 0];

internalField   uniform 0;

boundaryField
{
    walls
    {
        type            zeroGradient;
    }
}

dimensions [0 0 0 0 0 0 0];

internalField uniform 0;

boundaryField

{

walls

{

type zeroGradient;

}

typeを法線方向の勾配0のzeroGradientとしておきます。

初期状態の設定

初期状態としては水を93mm浸した状態を作るためsetFieldsで場の設定を行います。

system/setFieldsDict

defaultFieldValues
(
    volScalarFieldValue alpha.water 0
);

regions
(
    boxToCell
    {
        box (-0.5 -0.1 -0.1) (0.5 0.093 0.1);
        fieldValues
        (
            volScalarFieldValue alpha.water 1
        );
    }
);

defaultFieldValues

(

volScalarFieldValue alpha.water 0

);

regions

(

boxToCell

{

box (-0.5 -0.1 -0.1) (0.5 0.093 0.1);

fieldValues

(

volScalarFieldValue alpha.water 1

);

}

);

では初期状態を確認しましょう。
まず、先ほど境界条件で0.origで設定したので0フォルダという名前にコピーします。

cp -r 0.orig 0

1	cp -r 0.orig 0

setFieldsにより0フォルダ内はalpha.waterの分布を持ったファイルができてしまうので、0.origを残しておくと後で境界条件を設定しなおすときに便利です。

では、setFieldsを実行します。

setFields

setFields

ParaViewで結果を見てみましょう。

ParaViewのデフォルトが0フォルダを読みこまない設定になっているので、「Skip Zero Time」のチェックを外して読み込んでください。
図のように93mmに体積分率$\alpha=1$（水）が表示されていればOKです。

移動メッシュの設定

移動メッシュはconstant /dynamicMeshDictで行います。

dynamicFvMesh   dynamicMotionSolverFvMesh;

motionSolver    solidBody;

solidBodyMotionFunction tabulated6DoFMotion;

CofG            (0 0 0);
timeDataFileName "<constant>/1DoF_exp.dat";

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolver solidBody;

solidBodyMotionFunction tabulated6DoFMotion;

CofG (0 0 0);

timeDataFileName "<constant>/1DoF_exp.dat";

solidBodyMotionFunctionは以下のように選択ができます。

SDA
axisRotationMotion
drivenLinearMotion
linearMotion multiMotion
oscillatingLinearMotion
oscillatingRotatingMotion
rotatingMotion
tabulated6DoFMotion

SDA

axisRotationMotion

drivenLinearMotion

linearMotion multiMotion

oscillatingLinearMotion

oscillatingRotatingMotion

rotatingMotion

tabulated6DoFMotion

3次元6自由度の振動を与えるためtabulated6DoFMotionとしています。
CofGは重心位置のことです。
「timeDataFileName “<constant>/1DoF_exp.dat”;」で振動のデータを与えています。

1DoF_exp.dat”には決まったフォームでデータを書かないといけないので、エクセルで処理します。

SPHERIC_TestCase10/data_files/lateral_oil_1x.txtをExcelにタブ区切りで貼り付ける
以下のように入力用のフォームにする
コピーしてconstant/1DoF_exp.datにコピーする

="( " & A3 &"( (0 0 0) ( 0 0 " & F3 & ") )) "）

1	="( " & A3 &"( (0 0 0) ( 0 0 " & F3 & ") )) "）

constant/1DoF_exp.dat

167001
(
(  0 ( (0 0 0)  (0   0  0.000000001143341) ) )
(  0.00005 ( (0 0 0)  (0   0  0.001194) ) )
(  0.0001 ( (0 0 0)  (0   0  0.001194) ) )
(  0.00015 ( (0 0 0)  (0   0  0.001194) ) )
(  0.0002 ( (0 0 0)  (0   0  0.001194) ) )
(  0.00025 ( (0 0 0)  (0   0  0.001194) ) )
(  0.0003 ( (0 0 0)  (0   0  0.001194) ) )
(  0.00035 ( (0 0 0)  (0   0  0.001194) ) )
(  0.0004 ( (0 0 0)  (0   0  0.001194) ) )

(省略)

167001

(

( 0 ( (0 0 0) (0 0 0.000000001143341) ) )

( 0.00005 ( (0 0 0) (0 0 0.001194) ) )

( 0.0001 ( (0 0 0) (0 0 0.001194) ) )

( 0.00015 ( (0 0 0) (0 0 0.001194) ) )

( 0.0002 ( (0 0 0) (0 0 0.001194) ) )

( 0.00025 ( (0 0 0) (0 0 0.001194) ) )

( 0.0003 ( (0 0 0) (0 0 0.001194) ) )

( 0.00035 ( (0 0 0) (0 0 0.001194) ) )

( 0.0004 ( (0 0 0) (0 0 0.001194) ) )

(省略)

※余談ですがチュートリアルには数式で上記の振動データを作るプログラムがありました。こちらの数式に従いプログラムを作成してタンクの振動のインプットデータとして使うことができます。試しにやってみましょう。

型定義などはOpenFOAM用のC++って感じですがわりと簡単に関数を作れます。
作ったプログラムが以下です。

gen6DoF/gen6DoF.C

#include "List.H"
#include "vector.H"
#include "Vector2D.H"
#include "Tuple2.H"
#include "OFstream.H"
#include "fvCFD.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
    // End-time of the table
    const scalar endTime = 8;

    // Number of entries in the table
    const label nTimes = 100;

    // Amplitude of the translation [m]
    const vector transAmp(0, 0, 0);

    // Frequency of the translation [rad/s]
    const vector transOmega(0, 0, 0);

    // Amplitude of the rotation [deg]
    const scalar rotAmpz = 4;
    const vector rotAmp(0, 0, rotAmpz);

    // Frequency of the rotation [rad/s]
    const scalar L = 0.9;
    const scalar H = 0.093;
    const scalar g = 9.81;
    const scalar T1 = constant::mathematical::twoPi*std::pow(std::pow(constant::mathematical::pi*g/L*std::tanh(constant::mathematical::pi*H/L), 0.5),-1);
    const scalar rotOmegaz = constant::mathematical::twoPi /T1;
    const vector rotOmega(0, 0, rotOmegaz);

    Info << "T1 = " << T1 << endl;

    List<Tuple2<scalar,  Vector2D<vector>>> timeValues(nTimes);

    forAll(timeValues, i)
    {
        scalar t = (endTime*i)/(nTimes - 1);
        timeValues[i].first() = t;

        timeValues[i].second()[0] = vector
        (
            transAmp.x()*Foam::sin(transOmega.x()*t),
            transAmp.y()*Foam::sin(transOmega.y()*t),
            transAmp.z()*Foam::sin(transOmega.z()*t)
        );

        timeValues[i].second()[1] = vector
        (
            rotAmp.x()*Foam::sin(rotOmega.x()*t),
            rotAmp.y()*Foam::sin(rotOmega.y()*t),
            rotAmp.z()*Foam::sin(rotOmega.z()*t)
        );
    }

    {
        OFstream dataFile("MyGen1DoF.dat");
        dataFile << timeValues << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

#include "List.H"

#include "vector.H"

#include "Vector2D.H"

#include "Tuple2.H"

#include "OFstream.H"

#include "fvCFD.H"

using namespace Foam;

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

// Main program:

int main(int argc, char *argv[])

{

// End-time of the table

const scalar endTime = 8;

// Number of entries in the table

const label nTimes = 100;

// Amplitude of the translation [m]

const vector transAmp(0, 0, 0);

// Frequency of the translation [rad/s]

const vector transOmega(0, 0, 0);

// Amplitude of the rotation [deg]

const scalar rotAmpz = 4;

const vector rotAmp(0, 0, rotAmpz);

// Frequency of the rotation [rad/s]

const scalar L = 0.9;

const scalar H = 0.093;

const scalar g = 9.81;

const scalar T1 = constant::mathematical::twoPi*std::pow(std::pow(constant::mathematical::pi*g/L*std::tanh(constant::mathematical::pi*H/L), 0.5),-1);

const scalar rotOmegaz = constant::mathematical::twoPi /T1;

const vector rotOmega(0, 0, rotOmegaz);

Info << "T1 = " << T1 << endl;

List<Tuple2<scalar, Vector2D<vector>>> timeValues(nTimes);

forAll(timeValues, i)

{

scalar t = (endTime*i)/(nTimes - 1);

timeValues[i].first() = t;

timeValues[i].second()[0] = vector

(

transAmp.x()*Foam::sin(transOmega.x()*t),

transAmp.y()*Foam::sin(transOmega.y()*t),

transAmp.z()*Foam::sin(transOmega.z()*t)

);

timeValues[i].second()[1] = vector

(

rotAmp.x()*Foam::sin(rotOmega.x()*t),

rotAmp.y()*Foam::sin(rotOmega.y()*t),

rotAmp.z()*Foam::sin(rotOmega.z()*t)

);

}

{

OFstream dataFile("MyGen1DoF.dat");

dataFile << timeValues << endl;

}

Info<< "End\n" << endl;

return 0;

}

こちらをコンパイルします。

wmake

wmake

同じフォルダに「MyGen1DoF.dat」が出力されて先ほどのconstant/1DoF_exp.datとして使うことができます。

こちらのインプットデータで解析をしましたが、試験とタンクの振動周期がずれてきます。

試験が理論通りにいかないことはよくあることです。今回は振動のインプットデータは試験と同じにしたいので試験データの振動をそのままOpenFOAMのデータとして使います。
もし、振動が理論通りになっているかどうかも含めて知りたい、もしくは試験データがまだないと言った場合は理論式を使わざるを得ないのです。

物性値

物性値はconstant/transportPropertiesで行います。
こちらは資料通りに設定します。

constant/transportProperties

phases          (water air);

water
{
    transportModel  Newtonian;
    nu              8.96e-07;
    rho             998.2;
}

air
{
    transportModel  Newtonian;
    nu              1.48e-05;
    rho             1;
}

sigma           0.0728;

phases (water air);

water

{

transportModel Newtonian;

nu 8.96e-07;

rho 998.2;

}

air

{

transportModel Newtonian;

nu 1.48e-05;

rho 1;

}

sigma 0.0728;

離散化スキーム

離散化スキームはデフォルトの設定のままでよいです。

system/fvSchemes

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    div(rhoPhi,U)   Gauss vanLeerV;
    div(phi,alpha)  Gauss vanLeer;
    div(phirb,alpha) Gauss vanLeer;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

ddtSchemes

{

default Euler;

}

gradSchemes

{

default Gauss linear;

}

divSchemes

{

div(rhoPhi,U) Gauss vanLeerV;

div(phi,alpha) Gauss vanLeer;

div(phirb,alpha) Gauss vanLeer;

div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;

}

laplacianSchemes

{

default Gauss linear corrected;

}

interpolationSchemes

{

default linear;

}

snGradSchemes

{

default corrected;

}

対流項は2次精度のTVDスキームvanLeerを使用しています。

代数ソルバ

代数ソルバもデフォルトのままでよいです。

system/fvSolution

solvers
{
    alpha.water
    {
        nAlphaCorr      1;
        nAlphaSubCycles 3;
        cAlpha          1;
    }

    "pcorr.*"
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-05;
            relTol          0;
            smoother        DICGaussSeidel;
            cacheAgglomeration no;
        }

        tolerance       1e-05;
        relTol          0;
        maxIter         100;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.01;
        smoother        DIC;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       2e-09;
            relTol          0;
            nVcycles        2;
            smoother        DICGaussSeidel;
            nPreSweeps      2;
        }

        tolerance       2e-09;
        relTol          0;
        maxIter         20;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-06;
        relTol          0;
        nSweeps         1;
    }
}

PIMPLE
{
    momentumPredictor no;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    correctPhi      no;

    pRefPoint       (0 0 0);
    pRefValue       1e5;
}

relaxationFactors
{
    equations
    {
        "U.*"           1;
    }
}

solvers

{

alpha.water

{

nAlphaCorr 1;

nAlphaSubCycles 3;

cAlpha 1;

}

"pcorr.*"

{

solver PCG;

preconditioner

{

preconditioner GAMG;

tolerance 1e-05;

relTol 0;

smoother DICGaussSeidel;

cacheAgglomeration no;

}

tolerance 1e-05;

relTol 0;

maxIter 100;

}

p_rgh

{

solver GAMG;

tolerance 1e-08;

relTol 0.01;

smoother DIC;

}

p_rghFinal

{

solver PCG;

preconditioner

{

preconditioner GAMG;

tolerance 2e-09;

relTol 0;

nVcycles 2;

smoother DICGaussSeidel;

nPreSweeps 2;

}

tolerance 2e-09;

relTol 0;

maxIter 20;

}

{

solver smoothSolver;

smoother GaussSeidel;

tolerance 1e-06;

relTol 0;

nSweeps 1;

}

PIMPLE

{

momentumPredictor no;

nCorrectors 2;

nNonOrthogonalCorrectors 0;

correctPhi no;

pRefPoint (0 0 0);

pRefValue 1e5;

}

relaxationFactors

{

equations

{

"U.*" 1;

}

以下で基準圧を指定しています。

pRefPoint (0 0 0);
pRefValue 1e5;

pRefValueの値を0Paとしていると実験と同じ基準圧になります。

計算制御

controlDictでは計算制御の設定と物理量の出力設定を行います。

system/controlDict

application     interFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         8;

deltaT          0.01;

writeControl    adjustable;

writeInterval   0.05;

purgeWrite      0;

writeFormat     binary;//ascii;

writePrecision  6;

writeCompression on;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           0.5;

maxAlphaCo      0.5;

maxDeltaT       1;

functions
{
    probes
    {
        type            probes;
        libs            (sampling);
        writeControl   timeStep;
        writeInterval  1;
        probeLocations
        (
            (-0.450  0.093 0.031) //sensor1
            (-0.450  0.299 0.031) //sensor2
            (-0.425  0.508 0.031) //sensor3
            (0.400   0.508 0.031) //sensor4
            (0.450   0.483 0.031) //sensor5
            (0.450   0.222 0.031) //sensor6
        );
        fixedLocations  false;
        fields
        (
            p
        );
    }
}
じ

application interFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 8;

deltaT 0.01;

writeControl adjustable;

writeInterval 0.05;

purgeWrite 0;

writeFormat binary;//ascii;

writePrecision 6;

writeCompression on;

timeFormat general;

timePrecision 6;

runTimeModifiable yes;

adjustTimeStep yes;

maxCo 0.5;

maxAlphaCo 0.5;

maxDeltaT 1;

functions

{

probes

{

type probes;

libs (sampling);

writeControl timeStep;

writeInterval 1;

probeLocations

(

(-0.450 0.093 0.031) //sensor1

(-0.450 0.299 0.031) //sensor2

(-0.425 0.508 0.031) //sensor3

(0.400 0.508 0.031) //sensor4

(0.450 0.483 0.031) //sensor5

(0.450 0.222 0.031) //sensor6

);

fixedLocations false;

fields

(

);

}

じ

時間刻みは「maxCo 0.5; 」のクーラン数0.5で指定しています。
ちなみに「deltaT 0.01;」は初期値にしか使っていないので（初期値にも使ってないかも）、こちらの値を小さくしても時間刻みは変わらないので注意です。
時間刻みを変えたい場合はクーラン数を小さくしましょう。

functionsで指定した位置での圧力を出力する設定を行っています。

測定点は移動する場合と移動しない場合を考える必要があります。
メッシュ移動するけど初期の固定位置で測定したい場合はyes(or true)
メッシュ移動すると移動した分だけ測定位置も変わる場合はno
今回は容器に圧力計測を取り付けているので、fixedLocationsはnoとする必要があります。

fixedLocations Do not recalculate cells if mesh moves no true

以上で解析の設定が終了です。
次に計算を実行します。

計算実行

interFoamコマンドを実行して計算を実行します。
※今回は並列計算無しで計算をさせますが、8秒前の解析で数十分かかるので、待ってられない方は並列計算をさせてください。

interFoam > log.interFoam &

1	interFoam > log.interFoam &

interFoam：ソルバの実行
log.interFoam：log.interFoamにログを出力
&：バックグラウンドで流す

ParaViewで確認します。CAEは水しぶきの水がぼやけてしまっていますね。

実験の動画も載せておきます。

実験と比較するとCAEは解像度から水しぶきがぼやけているため、実機との乖離があります。今回使っているVOFによる解法はメッシュが粗いと空気と水との境界で体積分率がなまってしまうという問題があるため、ここからCAEの精度を上げるためには改善が必要です。

メッシュを細かくして解像度を上げる
ただし計算負荷がかかる
界面を追跡しながらダイナミックメッシュを行う
毎時刻メッシュを切り直すの計算負荷がかかる
乱流モデルを設定する
細かい渦などを見ることを諦めて乱流状態を数理モデルで表現する
こちらのように
界面圧縮項を大きくする
interIsoFoamソルバを使う

interIsoFoam
気液混相流に対してisoAdvector位相分率ベースの界面捕捉アプローチを用いたinterFoamから派生したソルバーで、適応的再メッシュを含むメッシュの動きやメッシュトポロジー変更を任意に行うことができます。interIsoFoamを使うと界面のぼやけるのがinterFoamと比べるとかなり改善されています。

今度のCAE精度アップの改善活動は必要ですが、本記事では深追いはしません。

圧力データの比較

実機とCAEの圧力データの比較を行うために簡単にPythonスクリプトを書きました。
今回はsensor1のみをグラフ化してみます。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

plt.figure(figsize=(16,10))
#data_p = np.loadtxt(f'./expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10/data_files/lateral_water_1x.txt')
df_p1_exp = pd.read_table(f"../expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10/data_files/lateral_water_1x.txt").iloc[:,0:2]#, columns=['time','p1','p2','p3','p4','p5','p6'])  
df_p1_exp['Pressure[Pa]'] = df_p1_exp['Pressure[mbar]']*100

caedata_p = np.loadtxt(f'./postProcessing/probes/0/p')
df_p1_cae = pd.DataFrame(caedata_p, columns=['time','CoGp', 'p1','p2','p3','p4','p5','p6']) 
df_p1_cae['p1(Pa)abs'] = df_p1_cae['p1']-df_p1_cae['p1'][0]

plt.title('p1 pressure',fontsize=24)

plt.plot(df_p1_exp.iloc[:,0],df_p1_exp.iloc[:,2],linestyle='solid',color='black',label="Exp")
plt.plot(df_p1_cae.iloc[:,0],df_p1_cae.iloc[:,8],linestyle='solid',color='red',label="OpenFOAM")
plt.xlabel('time(sec)',fontsize=16)
plt.ylabel('Pressure(Pa)',fontsize=16)
plt.legend(fontsize=24)
plt.grid()
plt.xticks(list(np.arange(0,8,1)),fontsize=12)
plt.yticks(list(np.arange(-1000,9000,1000)),fontsize=12)

import pandas as pd

import numpy as np

import matplotlib.pyplot as plt

plt.figure(figsize=(16,10))

#data_p = np.loadtxt(f'./expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10/data_files/lateral_water_1x.txt')

df_p1_exp = pd.read_table(f"../expdata/9449af_05f4aa0eaa1d4164a1024d07c562592a/SPHERIC_TestCase10/data_files/lateral_water_1x.txt").iloc[:,0:2]#, columns=['time','p1','p2','p3','p4','p5','p6'])

df_p1_exp['Pressure[Pa]'] = df_p1_exp['Pressure[mbar]']*100

caedata_p = np.loadtxt(f'./postProcessing/probes/0/p')

df_p1_cae = pd.DataFrame(caedata_p, columns=['time','CoGp', 'p1','p2','p3','p4','p5','p6'])

df_p1_cae['p1(Pa)abs'] = df_p1_cae['p1']-df_p1_cae['p1'][0]

plt.title('p1 pressure',fontsize=24)

plt.plot(df_p1_exp.iloc[:,0],df_p1_exp.iloc[:,2],linestyle='solid',color='black',label="Exp")

plt.plot(df_p1_cae.iloc[:,0],df_p1_cae.iloc[:,8],linestyle='solid',color='red',label="OpenFOAM")

plt.xlabel('time(sec)',fontsize=16)

plt.ylabel('Pressure(Pa)',fontsize=16)

plt.legend(fontsize=24)

plt.grid()

plt.xticks(list(np.arange(0,8,1)),fontsize=12)

plt.yticks(list(np.arange(-1000,9000,1000)),fontsize=12)

おおまかな圧力の挙動は再現できていますね。

ちなみにポリヘドラルメッシュに変えてみましょう。

これまた何とも言えない結果ですね。
ずいぶんノイジーな感じです。

ビタビタに実験とCAEがあることは難しいので、これをどう定量的に合っていると表現するか、どこまでの実機とCAEのコリレーションレベルで妥協するのかは納期や目的にも依るでしょう。

実験データはかなりノイジーなデータなので、n増し（複数回試験をすること）しないとノイズなのか本当にそういう振動なのかが判断付かないというのもあります。このような衝撃のある圧力データなどは1回のデータで判断せずにn増し試験を追加するのが賢明です。

追試験1

どうもOpenFOAMの結果の圧力が1.5秒3.5秒・・・あたりで負圧へ振動しているのが気になっていたのですが、原因が基準圧の座標位置の問題でした。

system/fvSolution

（省略）

PIMPLE
{
    momentumPredictor no;
    nCorrectors     2;
    nNonOrthogonalCorrectors 0;
    correctPhi      no;

    pRefPoint       (-0.45 0.508 0.031);
    pRefValue       0;
}

（省略）

（省略）

PIMPLE

{

momentumPredictor no;

nCorrectors 2;

nNonOrthogonalCorrectors 0;

correctPhi no;

pRefPoint (-0.45 0.508 0.031);

pRefValue 0;

}

（省略）

以下で基準圧を指定しています。

pRefPoint (-0.45 0.508 0.031));
pRefValue 0;

pRefValueの値を0Paとしていると実験と同じ基準圧になります。

これで余計な振動がなくなりました。

追試験2

実験データの「Position_original [deg]」の値が重複しているので、重複を除いたデータにして再度解析を行いました。

重複データの除去は以下のように適当にPythonで行いました。

import pandas as pd

df_data = pd.read_table("../expdata/SPHERIC_TestCase10/data_files/lateral_water_1x.txt")
df_data_new = pd.DataFrame()
i_list = []
data_old = ""

for i, data in enumerate(df_data["Position_original [deg]"]):
    if data != data_old:
        i_list.append(i)

    data_old = data

df_data_new = df_data.iloc[i_list]
df_data_new.to_excel("lateral_water_1x_new.xlsx", index=False)