平行平板間の完全に発達した乱流場を計算し，その流れ場を流入部に与えたバックステップ流れを計算する
計算に使ったファイルはGithubにアップしてある
≫mtkbirdman.com/OpenFOAM/
キーワード：ドライバ部，バックステップ流れ，チャネル乱流，LES，pitzDailyMapped
はじめに
OpenFOAMで，完全に発達したチャネル乱流をドライバ部としたバックステップ流れを計算する
チャネル乱流のチュートリアルはこちら
≫$FOAM_TUTORIALS/incompressible/pimpleFoam/LES/channel395
ドライバ部を設けたバックステップ流れのチュートリアルはこちら
≫$FOAM_TUTORIALS/incompressible/pisoFoam/LES/pitzDailyMapped
使用するLESモデルについてはこちら
≫【OpenFOAM v2006】非線形SGSモデルの実装
Pythonに関する情報はこちら
≫Pythonに関する情報 _ note.nkmk.me
概要
ドライバ部を設けたバックステップ流れは以下の手順で行う
- チャネル乱流を計算する
- 1.の計算結果をもとに，バックステップ流れの初期条件を作成する
- バックステップ流れを計算する
最終的には以下のような結果が得られるようになる
バックステップの上流側の四角で囲まれた領域がドライバ部であり，チャネル乱流の計算を行っている
それでは順番に説明していく
チャネル乱流の計算
まずはチャネル乱流の計算を行う
計算に用いたファイルはGithubを参照
≫mtkbirdman.com/OpenFOAM/BFS_DR/
注意すべき点は以下のとおり
- バックステップ流れと同じ条件設定(バルク速度，粘性係数，格子点)で行う
- 特に格子点については，
blockMeshDictで宣言する順番なども同じにする
blockMeshDictについては，まずバックステップ流れの格子を作ってから，ドライバ部以外のコードを消していくようにして作ればいいと思う
channel395の計算結果を格子点数の違う計算の初期条件として使いたい場合は，mapFieldを使えばよい
≫データのマッピング - PENGUINITIS
計算結果はこんな感じ
わかりやすいようにバックステップ流れの格子のfeature edgeも表示している
初期条件の作成
上で計算したチャネル乱流をもとに，バックステップ流れの初期条件を作成する
想定しているディレクトリ構成は以下のとおり
BFSのディレクトリにあるapplydriverchannel.pyを使って，BFS_DR/100の計算結果をもとにBFS/0/*の初期条件ファイルを作成する
$FOAM_USER_DIR
├── BFS/
│ ├── 0/
│ │ ├── U*
│ │ ├── k*
│ │ ├── nut*
│ │ └── p*
│ ├── constant/
│ │ ├── transportProperties*
│ │ └── turbulenceProperties*
│ ├── system/
│ ├── Allclean*
│ ├── Allrun*
│ ├── applydriverchannel.py*
│ └── BFS.foam*
└── BFS_DR/
├── 0/
├── 100/
│ ├── uniform/
│ ├── U*
│ ├── k*
│ ├── nut*
│ ├── :
│ └── p*
├── constant/
│ ├── polyMesh/
│ ├── transportProperties*
│ └── turbulenceProperties*
├── postProcessing/
├── system/
├── Allclean*
├── Allrun*
└── BFS_DR.foam*
applydriverchannel.pyの解説
ソースコード全体はGithubを参照
≫mtkbirdman.com/OpenFOAM/BFS/applydriverchannel.py
main
まずはメインの部分から
- ドライバ計算へのパス
sourceTimeと時刻0ディレクトリへのパスcaseTimeをキーボードから入力する - バックステップ流れのセル数をあらかじめ調べておき，
targetFieldSizeに設定する - 初期条件として必要な変数の
object(変数名)，class(ベクトル，スカラー，テンソル，対称テンソル)を設定する
import os
import sys
import numpy as np
import itertools
# input sourceTime and caseTime from keyboard
sourceTime = input("sourceTime: ")
caseTime = input("caseTime: ")
# set fieldname and fieldtype
fieldNames = ["U", "k", "nut", "p"]
fieldTypes = ["vector","scalar","scalar","scalar"]
targetFieldSize = 1132800
# write initial condition using driver calculation
for i in range(len(fieldNames)):
# read driver calculation file
sorceField = readVolField(caseTime=sourceTime,fieldName=fieldNames[i])
# copy driver calculation result to initial condition and set value 0 in region except driver's one
targetField = np.zeros((targetFieldSize,sorceField.shape[1]))
targetField[:sorceField.shape[0],:]=sorceField[:,:]
targetField[sorceField.shape[0]:,:] = 0.
# write volField file
volField = writeVolField(fieldType=fieldTypes[i], caseTime=caseTime, fieldName=fieldNames[i],volFieldValue=targetField)
utility
- リストの中から部分一致で検索してindexを返してくれる関数
inclusiveIndex inclusiveIndexは，該当するindexが1つならintを，複数あるならlistを返す
def inclusiveIndex(lst, purpose):
# search word "purpose" from "lst" and return the index
idxs = []
for idx, line in enumerate(lst):
if purpose in line:
idxs.append(idx)
if len(idxs) == 1:
idxs = idxs[0]
return idxs
volFieldの読み込み
volFieldのファイルを読み込む
- Pythonでファイルの読み込み、書き込み（作成・追記） _ note.nkmk.me
- Pythonで文字列の一部を削除（stripなど） _ note.nkmk.me
- NumPy配列ndarrayとPython標準のリストを相互に変換 _ note.nkmk.me
def readVolField(caseTime,fieldName):
# read volField file as list
PATH = os.path.join(caseTime, fieldName)
print(PATH)
with open(PATH, "r") as f:
volField = f.readlines()
# split line into values
volField = [line.strip().strip('(').strip(')').replace(' ',',').split(',') for line in volField]
# read data size written in next line of "internalField"
idx = inclusiveIndex(volField,"internalField")
datasize = int(volField[idx+1][0])
# trim lines except internalField and convert values from int to float
volField = np.array(volField[idx+3:datasize+idx+3],dtype='float64')
return volFieldvolFieldを書き出す
既存の初期条件ファイルのinternalFieldにドライバの計算結果を書き込む
internalFieldまでの行と，boudaryField以降の行はきちんと記述されているとする
def writeVolField(fieldType,caseTime,fieldName,volFieldValue):
PATH = os.path.join(caseTime, fieldName)
volField = []
# write header, footer and internalField as list
header, footer = writeHeaderAndFooter(fieldType=fieldType, caseTime=caseTime, fieldName=fieldName)
internalField = writeInternalField(volFieldValue=volFieldValue,fieldType=fieldType)
# concatenate header, footer and internalField
volField = header + internalField + footer
# write volField into raw file
print(PATH)
with open(PATH, mode='w') as f:
f.write('\n'.join(volField))
return volField
def writeHeaderAndFooter(fieldType, caseTime, fieldName):
# read volField file named as "fieldName"
PATH = os.path.join(caseTime,fieldName)
with open(PATH, "r") as f:
volField = f.readlines()
volField = [line.strip('\n') for line in volField]
# split volField into header and footer
header = volField[:inclusiveIndex(volField,"internalField") + 1]
footer = volField[inclusiveIndex(volField,"boundaryField") - 1:]
return header, footer
def writeInternalField(volFieldValue,fieldType):
internalField = ["{}".format(volFieldValue.shape[0])]
internalField.append("(")
if fieldType.upper() == "scalar".upper():
volFieldValue = [" ".join(map(str,line)) for line in volFieldValue.tolist()]
else:
volFieldValue = ["("+" ".join(map(str,line))+")" for line in volFieldValue.tolist()]
internalField = internalField + volFieldValue
internalField.append(");")
internalField.append("")
return internalField実行
以下のコマンドで実行できる
/BFS$ python applydriverchannel.py
sourceTime: ../BFS_DR/100
caseTime: ./0
../BFS_DR/100/U
./0/U
../BFS_DR/100/k
./0/k
../BFS_DR/100/nut
./0/nut
../BFS_DR/100/p
./0/p
/BFS$このプログラムによって，次のような初期条件が作成できる
ドライバ部以外の値が0に設定されていることが確認できる
バックステップ流れの計算
計算に使ったファイルはGithubを参照
≫mtkbirdman.com/OpenFOAM/BFS/
system/blockMeshDict
blockMeshDict内のboundaryで，inletのtypeをmappedPatchに設定し，流れ方向に6.4の位置の結果を流入条件として持ってきたいので，offset (6.4 0 0);と設定する
boundary
(
inlet
{
type mappedPatch;
offset (6.4 0 0);
sampleRegion region0;
sampleMode nearestCell;
samplePatch none;
faces ((0 10 19 9));
}
outlet
{
type patch;
faces ((4 14 15 5)(5 15 16 6));
}
bottomWall
{
type wall;
faces ((0 1 11 10)(1 2 12 11)(2 3 13 12)(3 4 14 13));
}
topWall
{
type wall;
faces ((9 8 18 19)(8 7 17 18)(7 6 16 17));
}
sides1
{
type cyclic;
neighbourPatch sides2;
faces ((0 1 8 9)(1 2 7 8)(2 3 4 5)(2 5 6 7));
}
sides2
{
type cyclic;
neighbourPatch sides1;
faces ((19 18 11 10)(18 17 12 11)(15 14 13 12)(17 16 15 12));
}
);0/U 0/k
inletをtype mapped;に設定する
Uは流入部で平均流速を1に固定したいので，setAverage true;，average (1 0 0);とする
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1706 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
location "0";
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField nonuniform List<vector>
1132800
(
(0.0271687 1.75494e-05 0.00252807)
:
(0.0 0.0 0.0)
);
boundaryField
{
inlet
{
type mapped;
value uniform (0 0 0);
interpolationScheme cell;
setAverage true;
average (1 0 0);
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (0 0 0);
}
"(topWall|bottomWall)"
{
type fixedValue;
value uniform (0 0 0);
}
"(sides1|sides2)"
{
type cyclic;
}
}
// ************************************************************************* ///*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v1706 |
| \\ / A nd | Web: www.OpenFOAM.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volScalarField;
location "0";
object k;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 2 -2 0 0 0 0];
internalField nonuniform List<scalar>
1132800
(
1.44104e-05
:
0.0
);
boundaryField
{
inlet
{
type mapped;
value uniform 0;
interpolationScheme cell;
setAverage false;
average 0;
}
outlet
{
type zeroGradient;
}
"(topWall|bottomWall)"
{
type fixedValue;
value uniform 0;
}
"(sides1|sides2)"
{
type cyclic;
}
}
// ************************************************************************* //結果
計算を実行する
/BFS$ ./Allrun
Running blockMesh on /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/BFS
Running decomposePar on /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/BFS
Running pimpleFoam (4 processes) on /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/BFS
Running reconstructPar on /home/mtkbirdman/OpenFOAM/mtkbirdman-v2006/BFS
/BFS$ノートPCだと5万ステップの計算に60時間ほどかかったので，実行するときは覚悟するように
結果はこんな感じ
まとめ
ドライバ部を設定したバックステップ流れの計算を行った
チャネル流れの計算結果をバックステップ流れの初期条件に設定するやり方は力任せなので，もっといい方法があるかもしれない
↓関連記事
コメント