【OpenFOAM】移流方程式で離散化スキームの勉強をする

こんにちは（@t_kun_kamakiri）

今回はOpenFOAMを使って移流方程式を解きたいと思います。その際に離散化スキームをどのように選択すると解の振る舞いが変わってくるのかを見ていきます。

本記事の内容

OpenFOAMの移流方程式を使って中心差分と風上差分の結果比較を行います。

以前にも移流方程式を差分法を使ってFortranで解いたことがありますが、移流方程式は式が簡単で解析解があるにも関わらず数値的に解こうとすると、不安定だったり拡散してしまったり難しい問題があります。

１次元の移流方程式

\begin{align*}\frac{\partial T(x,t)}{\partial t}+u\frac{\partial T(x,t)}{\partial x}=0\tag{1}\end{align*}

※$u$：一定速度

1次元の移流方程式は、下記の絵のように何かの物理量を形を変えずに一定速度$u$で運ぶ方程式を意味しています。

横軸は空間$x$ですが、縦軸は「密度」や「温度」だったり色々考えることができます。本記事では、温度$T$として扱うことにします。

Fortranを使った差分法の結果は以下のようになりました。

この結果を見ると中心差分は振動してしまい、前進差分と後退差分は流れの向きによっては安定だったり不安定だったりします。上の例のように右側に進行する($u>0$)場合は、安定に進みますが解が拡散してしまっています。
このように、流体の第二項の対流項$u\frac{\partial T}{\partial x}$の微分の離散化方法によっては、数値解が異なることがあります。異なると言っても、通常は流体のナビエストークスには流体粘性ので不連続解を分散させて解が振動しないようにしますし、不連続な解が出ることも稀なのでどのような離散化手法をとってもあまり解が変わりません。ただ、衝撃波が出ると速度や圧力の不連続解が生じ、さらに粘性項より対流項が大きくなるので上記に近い方程式になります。離散化スキームを以下のようにざっくり分けて特徴をまとめておきます。

離散化スキームざっくりまとめ

1次精度：風上差分(upwind)→こちらは安定だが解が拡散してしまう
2次精度：中心差分、QUICK→こちらは不連続解をシャープにとらえようとするが、数値振動が起こりやすい
TVDスキーム：minmod、vanLeer、vanAlbada、superBeeが有名。
色々な種類がありますが、不連続部分で振動を抑えるようにしながら、解はシャープにとらえようとするいい感じのスキーム

離散化スキームの種類はむちゃくちゃありますので、まずはこれくらいざっくりわけて特徴を理解しておくとよいでしょう。

あわせて読みたい

対流項の発散スキームの比較
こちらはOpenFOAMの離散化スキームの比較として一番参考になる記事でしょう。
Divergence scheme example
公式サイトにも離散化スキームに対する検討がされています。メッシュと離散化スキームの選択次第で結果がずいぶんと変わるのがわかるでしょう。
OpenFOAM6.2 Numerical schemes
OpenFOAM Basic Training
第1回オープンCAE講習会 OpenFOAM初中級講習B
upwind_difference.pdf (xrea.com)
こちらのpdfはスキームについてよくまとめられています。

こんな方が対象

OpenFOAMの移流方程式を使って離散化スキームの基礎的な内容を理解したい方
OpenFOAMのチュートリアルから独自に検証するケースファイルを作りたい方

離散化スキームについて勉強したい方の参考になれば幸いです。

チュートリアルをコピーする
メッシュ生成
境界条件の設定
初期状態の設定
物性値の設定
離散化スキームの設定
計算制御の設定
行列ソルバの設定
解析実行（中心差分）
解析実行（風上差分）
まとめ
参考書

チュートリアルをコピーする

スカラー輸送方程式のチュートリアルをコピーします。

cp -r $FOAM_TUTORIALS/basic/scalarTransportFoam/pitzDaily 001_scalarTransportFoam

1	cp -r $FOAM_TUTORIALS/basic/scalarTransportFoam/pitzDaily 001_scalarTransportFoam

「$FOAM_TUTORIALS」は環境変数で「/opt/OpenFOAM/OpenFOAM-v2012/tutorials」と割り当てられています（個人の環境に依る）。

ファイル構成は以下のようになっています。

.
├── 0
│   ├── T
│   └── U
├── constant
│   └── transportProperties
└── system
    ├── blockMeshDict
    ├── controlDict
    ├── fvSchemes
    └── fvSolution

├── 0

│ ├── T

│ └── U

├── constant

│ └── transportProperties

└── system

├── blockMeshDict

├── controlDict

├── fvSchemes

└── fvSolution

scalarTransportFoamは以下の偏微分方程式を解く非定常のソルバーです。

\begin{align*}
\frac{\partial T}{\partial t}+\nabla\cdot(\boldsymbol{u}T)-\nabla\cdot (D_{T}\nabla T)=\boldsymbol{S}_{T}\tag{2}
\end{align*}

あわせて読みたい

scalarTransportFoam

本記事では(2)を少し簡単にするために、

生成項を0：$\boldsymbol{S}_{T}=0$
拡散項を0：$D_{T}=0$
流速は一定：$\boldsymbol{u}=(0.1, 0, 0)$

と設定するようにします。

では、フォルダを移動します。

cd 001_scalarTransportFoam/

1	cd 001_scalarTransportFoam/

メッシュ生成

メッシュ生成はblockMeshを使います。

「system/blockMeshDict」ファイルを以下のように編集します。

system/blockMeshDict

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

scale  1.0;
xmax 1;
ymax 0.1;
zmax 0.1;

vertices        
(
   ( 0.0   0.0   0.0)
   ( $xmax 0.0   0.0)
   ( $xmax 0.0   $zmax)
   ( 0.0   0.0   $zmax)
   ( 0.0   $ymax 0.0)
   ( $xmax $ymax 0.0)
   ( $xmax $ymax $zmax)
   ( 0.0   $ymax $zmax)
);

blocks          
(
    hex (0 1 5 4 3 2 6 7)   (100 1 1) simpleGrading (1 1 1)
);

edges           
(
);

boundary         
(
    XMin
    {
        type patch;
        faces
        (
            (0 4 7 3)
        );
    }
    XMax
    {
        type patch;
        faces
        (
            (1 5 6 2)
        );
    }
    ZMin
    {
        type empty;
        faces
        (
            (0 1 5 4)
        );
    }
    ZMax
    {
        type empty;
        faces
        (
            (3 2 6 7)
        );
    }
    YMin
    {
        type wall;
        faces
        (
            (0 1 2 3)
        );
    }
    YMax
    {
        type wall;
        faces
        (
            (4 5 6 7)
        );
    }
);

mergePatchPairs
(
);

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

100

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class dictionary;

object blockMeshDict;

}

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

scale 1.0;

xmax 1;

ymax 0.1;

zmax 0.1;

vertices

(

( 0.0 0.0 0.0)

( $xmax 0.0 0.0)

( $xmax 0.0 $zmax)

( 0.0 0.0 $zmax)

( 0.0 $ymax 0.0)

( $xmax $ymax 0.0)

( $xmax $ymax $zmax)

( 0.0 $ymax $zmax)

);

blocks

(

hex (0 1 5 4 3 2 6 7) (100 1 1) simpleGrading (1 1 1)

);

edges

(

);

boundary

(

XMin

{

type patch;

faces

(

(0 4 7 3)

);

}

XMax

{

type patch;

faces

(

(1 5 6 2)

);

}

ZMin

{

type empty;

faces

(

(0 1 5 4)

);

}

ZMax

{

type empty;

faces

(

(3 2 6 7)

);

}

YMin

{

type wall;

faces

(

(0 1 2 3)

);

}

YMax

{

type wall;

faces

(

(4 5 6 7)

);

}

);

mergePatchPairs

(

);

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

以下のコマンドを実行してメッシュを作ります。

blockMesh

blockMesh

このように1.0×0.1×0.1の計算領域を作りました。

境界条件の設定

今回はxの正の方向へ流れが一様に0.1m/sで流れるようにします。

	U	T
XMin	type fixedValue; value uniform (0.1 0 0); or type zeroGradient;	type zeroGradient;
XMax	type zeroGradient;	type zeroGradient;
YMin	type slip;	type slip;
YMax	type slip;	type slip;
ZMin	type empty;	type empty;
ZMax	type empty;	type empty;

Uがデフォルトでは既に速度の分布があるファイルになっているので、以下のように書き換えます。

0/U

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volVectorField;
    location    "288";
    object      U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

internalField   uniform (0.1 0 0);
boundaryField
{
    XMin
    {
        type fixedValue;
        value uniform (0.1 0 0);
    }

    XMax
    {
        type zeroGradient;
    }

    YMin
    {
        type slip;
    }

    YMax
    {
        type slip;
    }

    ZMin
    {
        type empty;
    }

    ZMax
    {
        type empty;
    }
}


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

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class volVectorField;

location "288";

object U;

}

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

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

internalField uniform (0.1 0 0);

boundaryField

{

XMin

{

type fixedValue;

value uniform (0.1 0 0);

}

XMax

{

type zeroGradient;

}

YMin

{

type slip;

}

YMax

{

type slip;

}

ZMin

{

type empty;

}

ZMax

{

type empty;

}

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

流速は一定：$\boldsymbol{u}=(0.1, 0, 0)$です。
これは時間が進行しても常に同じ値です。

次に温度境界条件を設定します。

0/T

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    location    "0";
    object      T;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions      [0 0 0 1 0 0 0];

internalField   uniform 0;

boundaryField
{
    XMin
    {
        type            zeroGradient;
    }
    XMax
    {
        type            zeroGradient;
    }
    ZMin
    {
        type            empty;
    }
    ZMax
    {
        type            empty;
    }
    YMin
    {
        type            slip;
    }
    YMax
    {
        type            slip;
    }
}
// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class volScalarField;

location "0";

object T;

}

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

dimensions [0 0 0 1 0 0 0];

internalField uniform 0;

boundaryField

{

XMin

{

type zeroGradient;

}

XMax

{

type zeroGradient;

}

ZMin

{

type empty;

}

ZMax

{

type empty;

}

YMin

{

type slip;

}

YMax

{

type slip;

}

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

初期状態の設定

「setFields」で温度分布の初期条件を作ります。

このように$x$方向0.1~0.2だけ$T=1$とし、その他を$T=0$となるように設定します。

設定ファイルはsetFieldsDictに書きますが、デフォルトではファイルが無いので適当なチュートリアルから持ってくる必要があります。
有名なのはinterFoamのdamBreakにsetFieldsDictがありますが、それを知らない場合は以下のコマンドでsetFieldsDictが使われているチュートリアルを探すことができます。

find $FOAM_TUTORIALS -name "setFieldsDict"

1	find $FOAM_TUTORIALS -name "setFieldsDict"

コマンドを実行すると以下のようにずら～っと候補が出てきます。

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/system/setFieldsDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sloshingTank2D/system/setFieldsDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/driftFluxFoam/RAS/mixerVessel2D/system/setFieldsDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/setFieldsDict
/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/setFieldsDict
（以下省略）

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/system/setFieldsDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sloshingTank2D/system/setFieldsDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/driftFluxFoam/RAS/mixerVessel2D/system/setFieldsDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleMultiphaseInterFoam/laminar/damBreak4phase/system/setFieldsDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/twoLiquidMixingFoam/lockExchange/system/setFieldsDict

/opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/interIsoFoam/damBreakWithObstacle/system/setFieldsDict

（以下省略）

適当に以下からコピーをしてきます。

cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/system/setFieldsDict system/

1	cp -r /opt/OpenFOAM/OpenFOAM-v2012/tutorials/multiphase/compressibleInterDyMFoam/laminar/sphereDrop/system/setFieldsDict system/

以下のように編集します。

system/setFieldsDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      setFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

defaultFieldValues
(
    volScalarFieldValue T 0.0
    //volVectorFieldValue U ( 0.1 0 0 )
);

regions
(
    boxToCell
    {
        box ( 0.1 -1 -1 ) ( 0.2 1 1 );
        fieldValues
        (
            volScalarFieldValue T 1
        );
    }
);


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

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class dictionary;

location "system";

object setFieldsDict;

}

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

defaultFieldValues

(

volScalarFieldValue T 0.0

//volVectorFieldValue U ( 0.1 0 0 )

);

regions

(

boxToCell

{

box ( 0.1 -1 -1 ) ( 0.2 1 1 );

fieldValues

(

volScalarFieldValue T 1

);

}

);

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

デフォルトは「volScalarFieldValue T 0.0」とし、「box ( 0.1 -1 -1 ) ( 0.2 1 1 );」の領域だけを「volScalarFieldValue T 1」という値にします。

では、以下のコマンドで初期状態を設定します。

setFields

setFields

物性値の設定

今回は「拡散項を0：$D_{T}=0$」とします。物性値の設定は「constant/transportProperties」で行います。

constant/transportProperties

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

DT              0.0;


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

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class dictionary;

location "constant";

object transportProperties;

}

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

DT 0.0;

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

これにより、(2)のスカラー輸送方程式は移流方程式になります。

１次元の移流方程式

\begin{align*}\frac{\partial T(x,t)}{\partial t}+u\frac{\partial T(x,t)}{\partial x}=0\tag{1}\end{align*}

※$u$：一定速度

1次元の移流方程式は、下記の絵のように何かの物理量を形を変えずに一定速度$u$で運ぶ方程式を意味しています。

数値計算の結果が正しく解けていれば、上のように形を変えずに右側へ温度分布が進行していくはずです。

離散化スキームの設定

次に、離散化スキームを設定します。離散化スキームは「system/fvSchemes」で行います。

今回は移流方程式の$u\frac{\partial T(x,t)}{\partial x}$を移流項といいます。

この移流項の離散化方法によって大きく結果が変わることを確認します。

まずは、linear（中心差分）から見ていきます。

system/fvSchemes

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    default         none;
    // div(phi,T) Gauss linearUpwind grad(T); //線形風上差分
    div(phi,T) Gauss linear grad(T); //中心差分
    //div(phi,T) Gauss upwind grad(T); //風上差分
    // div(phi,T) Gauss QUICK grad(T); //QUICK
    // div(phi,T) Gauss vanLeer grad(T);
    // div(phi,T) Gauss limitedLinear 0.1 grad(T);
}

laplacianSchemes
{
    default         none;
    laplacian(DT,T) Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}


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

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class dictionary;

location "system";

object fvSchemes;

}

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

ddtSchemes

{

default Euler;

}

gradSchemes

{

default Gauss linear;

}

divSchemes

{

default none;

// div(phi,T) Gauss linearUpwind grad(T); //線形風上差分

div(phi,T) Gauss linear grad(T); //中心差分

//div(phi,T) Gauss upwind grad(T); //風上差分

// div(phi,T) Gauss QUICK grad(T); //QUICK

// div(phi,T) Gauss vanLeer grad(T);

// div(phi,T) Gauss limitedLinear 0.1 grad(T);

}

laplacianSchemes

{

default none;

laplacian(DT,T) Gauss linear corrected;

}

interpolationSchemes

{

default linear;

}

snGradSchemes

{

default corrected;

}

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

その他のスキームは別の記事で詳しく確認したいと思います。

計算制御の設定

解析時間や時間刻みは「system/controlDict」で設定を行います。

計算時間：5(秒）まで
時間刻み：0.001（秒）

system/controlDict

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     scalarTransportFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         5;

deltaT          0.001;

writeControl    timeStep;

writeInterval   100;

purgeWrite      0;

writeFormat     ascii;

writePrecision  6;

writeCompression off;

timeFormat      general;

timePrecision   6;

runTimeModifiable true;


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

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class dictionary;

location "system";

object controlDict;

}

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

application scalarTransportFoam;

startFrom startTime;

startTime 0;

stopAt endTime;

endTime 5;

deltaT 0.001;

writeControl timeStep;

writeInterval 100;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;

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

行列ソルバの設定

今回は特に設定を変えずにいきます。

system/fvSolution

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    T
    {
        solver          PBiCGStab;
        preconditioner  DILU;
        tolerance       1e-06;
        relTol          0;
    }
}

SIMPLE
{
    nNonOrthogonalCorrectors 0;
}


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

/*--------------------------------*- C++ -*----------------------------------*\

| ========= | |

| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |

| \\ / O peration | Version: v2012 |

| \\ / A nd | Website: www.openfoam.com |

| \\/ M anipulation | |

\*---------------------------------------------------------------------------*/

FoamFile

{

version 2.0;

format ascii;

class dictionary;

location "system";

object fvSolution;

}

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

solvers

{

solver PBiCGStab;

preconditioner DILU;

tolerance 1e-06;

relTol 0;

}

SIMPLE

{

nNonOrthogonalCorrectors 0;

}

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

解析実行（中心差分）

では計算を実行します。

scalarTransportFoam

1	scalarTransportFoam

計算が終わったらParaViewで確認しましょう。
中心差分だと温度分布の形を保ったまま右へ進行するどころか、振動してとても不安定な結果になってしまいました。

これはこちらの記事の内容と同じ結果になっていますね。

OpenFOAMは有限体積法でFortranでの結果は差分法ですが、移流項を中心差分で離散化するとどちらも振動する結果となっています。

ナビエストークス方程式には拡散項などがあるので振動も抑えられてしまいますが、発散項の中心差分での離散化は単独だと振動する解が出てきます。移流項を中心差分で解くと空間離散化の2次精度とは言え、振動する解となってしまうことは十分注意すべき点です。

解析実行（風上差分）

linear(中心差分で)は振動する解が得られたので今度は安定な差分法である風上差分を使います。以下「system/fvSchemes」でコメントアウトして風上差分を設定します。

system/fvSchemes

div(phi,T) Gauss upwind grad(T); //風上差分
// div(phi,T) Gauss linear grad(T); //中心差分

1 2	div(phi,T) Gauss upwind grad(T); //風上差分 // div(phi,T) Gauss linear grad(T); //中心差分

風上差分で拡散するやつ#移流方程式 pic.twitter.com/xFUtG5nKWR

— カマキリ🐲技術士勉強中 (@t_kun_kamakiri) February 1, 2023

Twitterにアップしたような拡散しながら右へ進行する解となりました。
こちらも記事の内容と同じ結果になっていますね。

まとめ

今回はスカラー輸送方程式のチュートリアルを変更して移流方程式にして、離散化スキームによる結果を見ました。

風上差分で拡散するやつ#移流方程式 pic.twitter.com/xFUtG5nKWR

— カマキリ🐲技術士勉強中 (@t_kun_kamakiri) February 1, 2023

微分を代数的に解くために離散化しているため、完ぺきな解法というものはなく一長一短ということですね。今回は中心差分と風上差分のみ紹介しましたが、その他に「線形風上差分」、「QUICK」、「TVD」など色々とあります。このようにスキームに対してどのような結果になるのかを知るためには離散化手法を知ることが大事ですが、何より計算で確かめてみることが良いでしょう。

有限体積法の「中心差分」による離散化はこちらにまとめています。1次元の熱拡散方程式ですが、参考になるでしょう。

あわせて読みたい

次回はスキームを色々変えて結果の違いを確かめたいと思います。複数ケースの計算を行うためスキームを変えて計算するプログラムを作って計算します。

「対流項の発散スキームの比較」にシェルスクリプトで自動で計算するプログラムがあります。シェルスクリプトはちょっと慣れないのでPythonで行いたいと思います。