OpenFOAM : Cavity flow チュートリアルをしっかりやってみます。その１
以下のチュートリアルケースをしっかりとトレースしてみました。
このチュートリアルでは、以下の形状で、上面以外は、壁に囲まれた空間です。上部壁は、x 方向に 1m/sで動いている。以下の図は、サイトから引用しています
初期は層流状態からスタートし、均一なメッシュで、icoFoam（非定常層流解析）ソルバーで、層流、等温、非圧縮流れとして解析します。
準備
設定ファイル、計算結果を保存するフォルダを作成、その中にチュートリアルケースをコピーします。以下のコマンドを実行します。。
cp -r $FOAM_TUTORIALS/incompressible/icoFoam/cavity/ ./フォルダ内には、
0 constant system の３つのフォルダがコピーされる。各フォルダに設定ファイルがあります。フォルダの構造は以下の用になる様です。
User Guideから引用しています。
メッシュの準備
OpenFOAMは、常に、３次元のデカルト座標系で、すべての形状は、３次元で認識されます。
z軸に垂直な面（0123 と 4567）に対して、empty境界条件を指定することで、２次元の流れをシミュレーションすることができます。
ここでは、メッシュのサイズは、z軸方向の幅を１単位とします。
z方向は、0.1の厚みとします。均一なメッシュを 20 x 20 で設定します。
system/blockMeshDict ファイルに生成したい形状、メッシュの定義を記します。今回は、コピーしてきたものをそのまま利用します。中身の確認だけ行います。
１ー７行目は、ヘッダー情報です。
８−１４行目は、FoamFileで、削除しても問題ない。（正確には理解できていない。）
17行目の scale で実際の問題の寸法の次元に調整することができる。ここでは、0.1倍としている。
３１−３４行目：blocks は、メッシュの基本サイズになる。
hex (0 1 2 3 4 5 6 6) で六面体、
(20 20 1) で、x方向に２０刻み、y方向に２０刻み、z方向は、1刻み
simpleGrading(1 1 1) で、各方向に等サイズにメッシュを切る
４０ー６９行目：boudnaryで境界条件を設定しています。
moivngWall {} の箇所は、type wall で壁と設定、face( (3 7 6 2) )で上面を指定している。
fixedWalls {}の箇所でも type wall で壁と設定、３面を指定している。
上面のみ移動することにしたいので、２つのグループに分けてwallを定義している。
frontAndBackは、Z軸に垂直な面で、type emptyとすることで、２次元のシミュレーションに相当する。（摩擦とか影響がないとしているのだろうと理解。）
movingWallなどの名前の付け方は、わかりやすいネーミングで良いと理解している。type と face の定義は、OpenFOAMのルールに従う。
1/*--------------------------------*- C++ -*----------------------------------*\
2| ========= | |
3| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
4| \\ / O peration | Version: v2006 |
5| \\ / A nd | Website: www.openfoam.com |
6| \\/ M anipulation | |
7\*---------------------------------------------------------------------------*/
8FoamFile
9{
10 version 2.0;
11 format ascii;
12 class dictionary;
13 object blockMeshDict;
14}
15// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
16
17scale 0.1;
18
19vertices
20(
21 (0 0 0)
22 (1 0 0)
23 (1 1 0)
24 (0 1 0)
25 (0 0 0.1)
26 (1 0 0.1)
27 (1 1 0.1)
28 (0 1 0.1)
29);
30
31blocks
32(
33 hex (0 1 2 3 4 5 6 7) (20 20 1) simpleGrading (1 1 1)
34);
35
36edges
37(
38);
39
40boundary
41(
42 movingWall
43 {
44 type wall;
45 faces
46 (
47 (3 7 6 2)
48 );
49 }
50 fixedWalls
51 {
52 type wall;
53 faces
54 (
55 (0 4 7 3)
56 (2 6 5 1)
57 (1 5 4 0)
58 );
59 }
60 frontAndBack
61 {
62 type empty;
63 faces
64 (
65 (0 3 2 1)
66 (4 5 6 7)
67 );
68 }
69);
70
71mergePatchPairs
72(
73);
74
75// ************************************************************************* //チュートリアルケースをコピーしたフォルダで
blockMeshを実行します。以下の応答が返ってきます。
/*---------------------------------------------------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2106 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
Build : _c15bfde3cb-20210624 OPENFOAM=2106
Arch : "LSB;label=32;scalar=64"
Exec : blockMesh
Date : Jul 25 2021
Time : 18:10:37
Host : penguin
PID : 3343
I/O : uncollated
Case : /home/$USER/ofV2106/cavity
nProcs : 1
trapFpe: Floating point exception trapping enabled (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster (fileModificationSkew 5, maxFileModificationPolls 20)
allowSystemOperations : Allowing user-supplied system call operations
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time
Creating block mesh from "system/blockMeshDict"
Creating block edges
No non-planar block faces defined
Creating topology blocks
Creating topology patches
Creating block mesh topology
Check topology
Basic statistics
Number of internal faces : 0
Number of boundary faces : 6
Number of defined boundary faces : 6
Number of undefined boundary faces : 0
Checking patch -> block consistency
Creating block offsets
Creating merge list (topological search)...
Deleting polyMesh directory "constant/polyMesh"
Creating polyMesh from blockMesh
Creating patches
Creating cells
Creating points with scale 0.1
Block 0 cell size :
i : 0.005 .. 0.005
j : 0.005 .. 0.005
k : 0.01 .. 0.01
There are no merge patch pairs
Writing polyMesh with 0 cellZones
----------------
Mesh Information
----------------
boundingBox: (0 0 0) (0.1 0.1 0.01)
nPoints: 882
nCells: 400
nFaces: 1640
nInternalFaces: 760
----------------
Patches
----------------
patch 0 (start: 760 size: 20) name: movingWall
patch 1 (start: 780 size: 60) name: fixedWalls
patch 2 (start: 840 size: 800) name: frontAndBack
End次に、ParaViewは独立してインストールして利用することにしています。
touch post.foamと同じフォルダ内で入力して、post.foamというファイルを作成します。
ParaViewを起動して、post.foamを読み込むと以下の画面になります。Applyをクリックします。
以下のSolidな四角が表示されます。黄色で囲ったSurfaceをWireframeに変更します。
以下の用にメッシュが表示されます。x, y方向に20個、z方向に1個分のメッシュが確認できます。
blockMeshDictファイル内の blocks 内の (20 20 1)を変更、blockMeshを実行するとその数値にあわせた刻み数に変更されます。
境界条件、初期条件
シミュレーション時刻 t = 0 secでの初期条件を設定します。０というフォルダがあり、ここに初期、境界条件を設定します。中には、p、Uというファイルがあり、それぞれ圧力と、速度を定義します。
以下は、内容です。ヘッダー部分は割愛しています。
17dimensions [0 2 -2 0 0 0 0];
18
19internalField uniform 0;
20
21boundaryField
22{
23 movingWall
24 {
25 type zeroGradient;
26 }
27
28 fixedWalls
29 {
30 type zeroGradient;
31 }
32
33 frontAndBack
34 {
35 type empty;
36 }
37}
38
39// ************************************************************************* //17行目：dimensionsは、このファイルの数値の次元を定義しています。次元の定義方法は、以下のテーブルのNo.で次元がわりあてられています。
[0 2 -2 0 0 0 0] は、Length(m)^2 x Time(s)^(-2) = m2/s2 になります。
圧力の次元なので、N / m2 = kg m/s^2 / m2 = kg / m / s^2 で、
[1 -1 -2 0 0 0 0]; となりそうですが、非圧縮性流体ソルバーでは、
の方程式を解いているとのことです。密度の変化がないため、両辺を密度で割った次元になっています。
圧力は、 kg / m / s^2 / (kg/m^3) = m^2 / s^2 の次元になっているとのことです。
テーブルは、以下のリンクより引用。
１９行目：internalFieldは、内部の状態を定義する。uniformは一様な値になる。 nonuniformの場合は、フィールドのすべての場所に値を設定する必要がある。多くのケースでは、一様に0の設定から始めるので問題ないと考えられるとのことです。（動的に変わる値は計算される。）
２１−３７行目：boundaryFieldは、境界条件の定義です。fixedWallsと、movingWallsの壁からの圧力勾配は、ゼロ、z軸に垂直な面は、emptyで2Dシミュレーションを実行することを意味しています。
同じく、Uファイル、速度の境界、初期条件について確認する。
17dimensions [0 1 -1 0 0 0 0];
18
19internalField uniform (0 0 0);
20
21boundaryField
21{
22 movingWall
23 {
24 type fixedValue;
25 value uniform (1 0 0);
26 }
27 fixedWalls
28 {
29 type noSlip;
30 }
31 frontAndBack
32 {
33 type empty;
34 }
35}
36
37// ************************************************************************* //１７行目：dimension [0 1 -1 0 0 0 0]; は、 m s^-1 の次元になります。（速度の通りになる。）
１８行目：internalField、内部空間の初期の速度は、一様に0とする
２１−３５行目：boundaryField、境界条件は、
movingWall（上面）は、固定値（fixedValue）で、x方向に 1m/sで一定で移動（ uniform(1 0 0) ）するとしています。
fixedWall（壁）は、noSlip で壁上では流体は動かない。速度が0となります。壁上で速度を持つときは、slipとします。（速度勾配はゼロになる）
frontAndBack（z軸に垂直な面）は、emptyで2Dシミュレーションを実行することを意味しています。
物性値
icoFoam（非定常層流解析）ソルバーで必要になる物性は、動粘度のみになります。
constant/transportPropertiesに、物性定義があります。
動粘度に対応するキーワードは、nuになります。このnuの値を定義します。
コピーした値をそのまま利用します。
17
18nu 0.01;
19
20
21// ************************************************************************* //代表長さ、0.1 m （溝の幅）、代表速度は、1m/sとして、動粘度が0.01[m2/s] とするとレイノルズ数は、10と判断されます。
コントロール（シミュレーション制御）
シミュレーションする時間、計算値の出力間隔などを定義します。
system/controlDict ファイルになります。
２２行目、startTimeは、0からにしておきます。。この数値から初期値のフォルダを参照することになります。（計算結果があるときに途中から続きを見たいときには、先に行った計算結果の時間をスタートタイムとすれば続きから計算がされるだろうと理解。）
２０行目、startFromは、startTimeと同じにしておくため、startTimeと記述する。
２６行目：endTimeは、シミュレーションの完了時間になります。一般論として、層流状態では、空間の10倍の距離を流れれば、定常に達していると考えられています。移動が1 m/sで、距離が0.1mであるので、 1s で十分だろうと判断されます。（チュートリアルの説明分では、0.5sで定常に達していたと書かれています。）
２４行目：stopAtで、シミュレーションの終了時刻の指定をすることができます。endTimeとおなじになるように、endTimeと記述しています。
２８行目：deltaTは、非定常計算の時刻の刻み幅になります。icoFoamソルバーでの計算の精度、安定性を考慮して決定する必要があります。チュートリアルでは、0.005 s となっています。
精度と計算の安定性の確認としては、クーラン数が１以下になっていることを確認します。
|U|速度の絶対値は、セル内での速度の大きさ、ここでは壁が1 m/sで動いているので、1 m/s、
δxは、計算空間のセルのサイズで、0.1m の幅を20刻みにしているので、0.005 m になります。
δt は、ここでの指定刻み時間で、クーラン数条件を満たすように指定します。
となり、0.005 sが指定されています。
３０行目：timeSpteと、３２行目：writeIntervalは、２０としているので、0.005 s x 20 = 0.1 sとなり、この値がtimeStepに自動で代入されます。この結果、0.1 sごとに結果を出力するという指定になります。
OpenFOAMは、計算を実行するとこのtimeStepごとにフォルダを作成し、その時刻のシミュレーション結果を保存していきます。
17
18application icoFoam;
19
20startFrom startTime;
21
22startTime 0;
23
24stopAt endTime;
25
26endTime 0.5;
27
28deltaT 0.005;
29
30writeControl timeStep;
31
32writeInterval 20;
33
34purgeWrite 0;
35
36writeFormat ascii;
37
38writePrecision 6;
39
40writeCompression off;
41
42timeFormat general;
43
44timePrecision 6;
45
46runTimeModifiable true;
47
48
49// ************************************************************************* //離散化と線形ソルバーの設定
system/fvSchemesファイル内で、有限体積法の離散化方法を指定します。
system/fvSolutionファイル内で、離散化された線形代数方程式の解法（行列計算）を指定します。
指定できるオプションは、User Guideに詳しくあるようです。（まだ詳しく理解できていないところです。）
計算の実行
<case>フォルダ（contantフォルダと、systemフォルダがある場所）で、
icoFoamを実行すると計算状況が表示されます。
描画
メッシュの確認のときに作った、post.foamをparaviewで開いて、表示してみます。
ファイルを指定して、OKボタンを押します。
Applyボタンを押します。
初期では、圧力が描画されいます。
赤丸のところを、圧力から速度に変えてみました。青矢印で示しているボタンをおしてみます。
上の画面が開くので、Jetを選択してみます。（Applyボタンをおして、Closeします。）
速度のベクトル表示を入れたいと思います。下の図にしめしたメニューにあるボタン（Glyph）を押します。
Pipeline BrowserにGlyph1が出てきます。それを選択すると下のPropertiesタブ内で設定ができます。向き（Orientation）は、速度U、Scale Arrayは、No Scale Arrayとしています。Scale Factorを初期値が0.01だったので、半分の0.005としました。Applyボタンを押します。
メニューバーになる、赤丸で囲んだ箇所を速度のU、Surfaceとすると、ベクトルの強度を配色で表せます。配色は先程、速度にあてた配色になります。
上の図では、post.foamを選んで、コンター図を消して、Edgeだけ表示にしています。下の図を参照してください。
ストリームラインの描画は、メニュー > Filters > CosmoTools > Stream Tracerを選択します。
Pipeline7 Browser上にStreamTracer1がでてくるので、Applyボタンをして、表示します。他の描画は消しています。（post.foamなど、名前の前の目玉アイコンをクリックして目玉が閉じた状態にしています。）
Stream Lineのシードポイントは、ジオメトリの中心を垂直に通る線源に沿って指定しています。Resolutionを100程度にしています。（ラインの数を適当に調整しています。）
Applyボタンを押すと、以下の図が表示されます。
（上の図は、配色は、デフォルトの青ー赤のものからな、予め定義されているJet配所似変更しています。）
まとめ
openfoam.comのHPのチュートリアル、2.1 Lid-driven cavity flow にしたがって、インストールしたopenfoamのtutorialからcavityケースのコピー、実行まで行った。
結果は、paraviewで参照した。
速度について、コンター図とベクトル図、ストリームラインまで書いてみた。
所感
このチュートリアルにはつづきがあり、メッシュをより細かくきったものや、等間隔メッシュから、勾配が急になりそうな箇所だけメッシュを細かくするなどの説明がある。
続きを実行して、書いていきたい思う。
よろしければサポートをお願いします。講習会、有料情報の取得などにあてたいと考えています。やっぱり、単純にうれしいです。