【OpenFOAM 2006】ドライバ部を設けたバックステップ流れの計算

平行平板間の完全に発達した乱流場を計算し,その流れ場を流入部に与えたバックステップ流れを計算する

計算に使ったファイルは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. チャネル乱流を計算する
  2. 1.の計算結果をもとに,バックステップ流れの初期条件を作成する
  3. バックステップ流れを計算する

最終的には以下のような結果が得られるようになる

バックステップの上流側の四角で囲まれた領域がドライバ部であり,チャネル乱流の計算を行っている

最終時刻における速度場(上)および圧力場(下)
速度勾配テンソルの第二不変量の等値面:時刻0(上)および最終時刻(下)

それでは順番に説明していく

チャネル乱流の計算

まずはチャネル乱流の計算を行う

計算に用いたファイルは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のファイルを読み込む

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 volField
volFieldを書き出す

既存の初期条件ファイルの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で,inlettypemappedPatchに設定し,流れ方向に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

inlettype 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時間ほどかかったので,実行するときは覚悟するように

結果はこんな感じ

最終時刻における瞬時値:速度場(上),圧力場(下)
速度勾配テンソルの第二不変量の等値面:時刻0(上),最終時刻(下)
時間平均:速度場(上),圧力場(下)

まとめ

ドライバ部を設定したバックステップ流れの計算を行った

チャネル流れの計算結果をバックステップ流れの初期条件に設定するやり方は力任せなので,もっといい方法があるかもしれない

Sponsored Link

コメント

タイトルとURLをコピーしました