2020年4月11日土曜日

SU2のチュートリアル

はじめに

OpenFoamは何度か使ってことあるので、SU2も使ってみました。
SU2は非構造な計算格子データ(非構造格子)に対応したOSSの圧縮性流体ソルバーです。
https://su2code.github.io/
圧縮性流体ソルバーで広く使われる計算スキームは2000年ころにはだいたい固まってあまり変化しなくなりました。
2012年に登場したこのソフトウェアも割合オーソドックスなものを踏襲していますが、いくつか新しいものも取り入れられているようです。

Binaryファイルをダウンロードしてそのまま使うのも良いのだけど、ソースをこのページの内容に従って取得してビルドしてみました。

以下にソースファイルがあるものとします。
$HOME/SU2/
OSはUbuntu 16.04です。

SU2はver6以前とver7以降でビルドの仕方が少し変わったようです。

SU2-v6.0

いつものconfigureでビルドしていきます。
git clone -b v6.0.0 https://github.com/su2code/SU2.git
mv SU2 SU2-6.0.0
cd ./SU2-6.0.0
./bootstrap
./configure --prefix=$HOME/SU2/SU2-v7.0.3 --enable-mpi
make -j 4
make install
bin/以下に実行モジュールがインストールされます。

SU2-v7.0

Blackbird(v7.0)からmeson, ninjaを使うとのこと。
git clone -b v7.0.3 https://github.com/su2code/SU2.git
mv SU2 SU2-7.0.3
cd ./SU2-7.0.3
./meson.py build --prefix=$HOME/SU2/SU2-v7.0.3
./ninja -C build install
なお、このままではうまく行かなかったので、こちらの環境では以下変更を要しました。
/meson_scripts/extract_file.py
  with tarfile.open(fn) as z:
    z.extractall(str(outpath))
  with tarfile.open(str(fn)) as z:
    z.extractall(str(outpath))

チュートリアル

SU2-v7.0.3で非定常な翼型周りを対象にチュートリアルを行ってみます。
git clone https://github.com/su2code/Tutorials
Unsteady_NACA0012に移動。
cd Tutorials/compressible_flow/Unsteady_NACA0012
計算設定ファイルはunsteady_naca0012.cfg、計算格子ファイルはunsteady_naca0012_mesh.su2です。 計算実行は単列計算では
$HOME/SU2/SU2-v7.0.3/bin/SU2_CFD unsteady_NACA0012.cfg
並列計算(たとえば8並列)であれば、
mpirun -np 8 $HOME/SU2/SU2-v7.0.3/bin/SU2_CFD unsteady_NACA0012.cfg
デフォルトではVTKフォーマットで可視化ファイルが生成されます。Paraviewで可視化します。
後流の流れ場が不連続に見えます。
今回の格子は構造格子であればC型格子に相当するもので、翼後端から後流にかけてとても格子が密になっています。
後流を拡大してあげると分布は連続になっていて、おそらく収束が足りないのではないかと思います。

2020年3月18日水曜日

OpenFoam, CalculiX adapterのインストール

はじめに

流体構造連成解析(FSI)をpreCICEを用いて行いたい。
preCICEはインストール済みで、次はソルバーのアダプターを用意します。
ここではOpenFOAMアダプターとCalculiXアダプターでFSIを行ってみます。
なお、アダプターのBuilding(OF), Builiding(CalculiX)にすべて書いてあります。

OpenFOAMアダプター

最初にアダプターをDLする。
git clone https://github.com/precice/openfoam-adapter.git
次に該当のディレクトリに移動して
./Allwmake
とくにエラーもなく終了すればそのままOFカップリングのチュートリアルに進めばよいのですが、preCICEのバージョンによってはうまくいかないこともある。

CalculiXアダプター

アダプターをDLし、コンパイルします。
CalculiXディレクトリに移動して、
wget https://github.com/precice/calculix-adapter/archive/master.zip 
unzip master.zip 
cd calculix-adapter-master 

Makefileを編集する。SPOOLES, CCX, ARPACK, PRECICE_ROOT, YAMLのパスを正しく設定します。
以下でインストールします。
make
ここもとくに問題終了すればよいのですが、preCICEのバージョンでうまくいかないこともあります。

FSIチュートリアル

OpenFoam-CalculiXのチュートリアルを行ってみます。
基本的にここを参考にすればよいだけなのですが、まずは以下のチュートリアルケースを取得します。
https://github.com/precice/tutorials
このあと、対象のディレクトリに移動して、
cd tutorials/FSI/flap_perp/OpenFOAM-CalculiX
./Allrun

2020年2月25日火曜日

preCICEインストール(修正)

はじめに

preCICE(Precise Code Interaction Coupling Environment)を手元にインストールしましたので、わすれる前にメモを残しておきたい。これは複数の独立したソフトウェアを接続し、マルチフィジックスの実行を容易にしてくれるカップリングライブラリとのこと。

preCICE Wikiに従ってインストールしています。ただしMPIによる並列を使うソフトウェアの場合、一つ一つソースからインストールしていく必要があるようです。

サーバー

サーバーOS:Ubuntu Server 16.04.6 LTS

準備

最初に
sudo apt-get update
sudo apt-get install build-essential g++ python-dev autotools-dev libicu-dev build-essential libbz2-dev
sudo apt install build-essential libeigen3-dev libxml2-dev
sudo apt install gfortran

Boost, YAML, MPICH, PETSC、(OPENFOAM)のインストールが必要です。

boost (1.66): boost1.65.1以上を要するとのこと
wget https://dl.bintray.com/boostorg/release/1.66.0/source/boost_1_66_0.tar.gz
tar xzvf boost_1_66_0.tar.gz
cd boost_1_66_0/
./bootstrap.sh --prefix=/usr/local/
./b2 --help
./b2 -j8 variant=release link=shared threading=multi runtime-link=shared
sudo ./b2 -j8 install 

YAML(0.5.3)
wget https://github.com/jbeder/yaml-cpp/archive/yaml-cpp-0.5.3.zip
 unzip yaml-cpp-0.5.3.zip
 cd yaml-cpp-yaml-cpp-0.5.3
--静的ライブラリ
mkdir build
cd build
cmake ..
make
--動的ライブラリ
cd yaml-cpp-yaml-cpp-0.5.3
mkdir build2
cd build2
cmake -DBUILD_SHARED_LIBS=ON ..
make
.bsahrcに以下追加
export CPLUS_INCLUDE_PATH=$HOME/app/yaml-cpp-yaml-cpp-0.5.3/include:$CPLUS_INCLUDE_PATH
export LD_LIBRARY_PATH=$HOME/app/yaml-cpp-yaml-cpp-0.5.3/build2:$LD_LIBRARY_PATH

BoostとYAMLについて
Boostが1.67以上のときはYaml-ccp0.6以上, Boost1.66以下のときはYaml-ccp0.6より古いバージョンにしないといけない。

MPICH
適当なバージョンをDL(ここではpetscとして3.3, OpenFOAM-v1812として3.2。厳密にはバージョンは合わせるべき)
wget http://www.mpich.org/static/downloads/3.3.2/mpich-3.3.2.tar.gz
tar xvf  mpich-3.3.2.tar.gz
cd mpich-3.3.2
./configure --prefix=/usr/lib/mpich/mpich-3.3.2 |& tee log_c.txt
make 2>&1 | tee log_m.txt
sudo make install |& tee log_mi.txt
以下を.bashrcに追加
MPI_ROOT=/usr/lib/mpich/mpich-3.3.2
PATH=$MPI_ROOT/bin:$PATH
LD_INCLUDE_PATH=$MPI_ROOT/include:$LD_INCLUDE_PATH
LD_LIBRARY_PATH=$MPI_ROOT/lib:$LD_LIBRARY_PATH
export MPI_ROOT PATH LD_INCLUDE_PATH LD_LIBRARY_PATH

PETSC
まずは準備として
sudo apt install libopenblas-base
sudo apt install libopenblas-dev
git clone -b maint https://bitbucket.org/petsc/petsc petsc
cd petsc
./configure --with-debugging=0 --with-mpi-dir=/usr/lib/mpich/mpich-3.3.2
make PETSC_DIR=$HOME/app/petsc PETSC_ARCH=arch-linux2-c-opt all
.bashrcに以下:
export PETSC_DIR=$HOME/app/petsc
export PETSC_ARCH=arch-linux2-c-opt
export LD_LIBRARY_PATH=$PETSC_DIR/$PETSC_ARCH/lib:$LD_LIBRARY_PATH
export CPATH=$PETSC_DIR/include:$PETSC_DIR/$PETSC_ARCH/include:$CPATH
export LIBRARY_PATH=$PETSC_DIR/$PETSC_ARCH/lib:$LIBRARY_PATH
export PYTHONPATH=$PETSC_DIR/$PETSC_ARCH/lib

OpenFOAM
MPIはデフォルトではsystemopenmpiになっている。 preciceはopenmpiはだめらしい(ポート?)ので、mpichに変更する。 ここではmpich-3.2としているが、3.3が良いかも
cd path/openfoam1812/OpenFOAM-v1812/etc
bashrcを編集:
#export WM_MPLIB=SYSTEMOPENMPI
export WM_MPLIB=MPICH
source bashrc
cd ../../ThirdParty-v1812
wget http://www.mpich.org/static/downloads/3.2.1/mpich-3.2.1.tar.gz
tar xfz mpich-3.2.1.tar.gz
mv mpich-3.2.1 mpich-3.2
./makeMPICH

preCICEインストール

まずはファイルの取得(xyzはバージョン、例えばx.y.z=1.6.1)
wget https://github.com/precice/precice/archive/v1.6.1.tar.gz
tar -xzvf v1.6.1.tar.gz
cd precice-1.6.1/
つぎに
mkdir build && cd build
cmake -DBUILD_SHARED_LIBS=ON -DCMAKE_BUILD_TYPE=RelWithDebInfo  ..
make -j $(nproc)
テスト
make test_base
make test
インストール
sudo make install
(CMake optionして指定しなければ/usr/localにインストールされる)

2019年4月10日水曜日

pimpleDyMFoamを用いた流体と剛体運動の連成解析


はじめに

運動している剛体周りの流体挙動とその連成を数値解析的手法で調べる。
OpenFoamのソルバーにpimpleDyMFoamがある。移動格子での流体解析が実装されており、これによって剛体周りの流体とその運動を記述できるかもとのこと。
ここでは1軸自由度の流体運動連成解析を考える、剛体運動は以下の運動方程式で記述する。
\begin{align} I\frac{\partial^2 \theta}{\partial t^2} = M \end{align}
流体によるモーメントから上記方程式に従って迎角が得られる、それを移動格子を通じて流体解析にフィードバックする。格子移動時の物理量内挿はALEでやってる?

以下、いろいろ適当な点/間違いあるかと思います。

計算対象

なんでも良いのだけどGrabCADから適当なumbrellaのSTLを見つけてくる(傘の運動に興味があった)
umbrella.stl

適当にリスケールしています、だいたい代表長0.8mくらい。

計算条件

ヴァージョン:OpenFoam v5.0
サーバーOS:Ubuntu Server 16.04.6 LTS
並列数:72

ベースにしたチュートリアルケース:/opt/openfoam/openfoam5/openfoam5/tutorials/incompressible/pimpleDyMFoam/wingMotion
(インストール場所/opt/openfoam/openfoam5です)

まず次のようなスクリプトを準備
./qfoam_header.sh
#!/bin/bash

ncpu=72
dir_pre_mesh=mesh
#dir_main_simple=simpleFoam
dir_main_dym=pimpleDyMFoam

source /opt/openfoam/openfoam5/openfoam5/etc/bashrc

# Source tutorial run functions
. $WM_PROJECT_DIR/bin/tools/RunFunctions

計算格子の作成

計算格子はsnappyHexMeshで作成する。
流体解析的な精度は今回は気にしないので格子は真面目に作らない。200万cells数くらいと思う。

以下を編集
{dir_pre_mesh}/system/controlDict
{dir_pre_mesh}/system/blockMeshDict
{dir_pre_mesh}/system/snappyHexMeshDict
{dir_pre_mesh}/system/surfaceFeatureExtractDict
{dir_pre_mesh}/system/decomposeParDict

#!/bin/bash

. qfoam_header.sh

# Make mesh
cd ${dir_pre_mesh}

cp -r 0.orig 0

surfaceFeatureExtract > log_surfaceFeatureExtract 2>&1
blockMesh > log_blockMesh 2>&1
decomposePar > log_decomposePar_1 2>&1
mpirun -np $ncpu snappyHexMesh -overwrite -parallel > log_snappyHexMesh 2>&1
#snappyHexMesh -overwrite > log_snappyHexMesh 2>&1

# For parallel process
reconstructParMesh -constant > log_reconstructParMesh 2>&1
rm -r processor*
decomposePar > log_decomposePar_2 2>&1

計算格子を作成した。

なぜか傘の柄の一部が消えてしまった。。

格子ファイルのコピー

メイン計算は別ディレクトリで行う。
#!/bin/bash

. qfoam_header.sh

# Copy the mesh
cd ${dir_main_dym}
cp -r ../${dir_pre_mesh}/processor* ./
cp -r ../${dir_pre_mesh}/0 ./
cp 0.orig/pointDisplacement ./0/
cp -r ../${dir_pre_mesh}/constant/polyMesh ./constant/
#mapFields ../${dir_pre_mesh} -parallelTarget -consistent -sourceTime latestTime -parallelSource > log_mapFields 2>&1
#
count=0
while [ $count -lt $ncpu ] ;
do
  mv processor$count/0/pointDisplacement.unmapped* processor$count/0/pointDisplacement
  count=`expr $count + 1`
done

物理モデル・計算モデル

分子粘性係数:1.e-5 Pa.s
乱流モデル:なし(層流計算)
一様流速度:10 m/s
一様流密度:1 kg/m3 (メイン計算では使わないけど空気力算出に必要)
Re数はだいたい1e6

以下を編集
{dir_main_dym}/constant/transportProperties
{dir_main_dym}/constant/turbulenceProperties
{dir_main_dym}/constant/dynamicMeshDict
慣性モーメントや重心、質量はすごく適当な値になっている。

流体+移動格子計算

次を編集
{dir_main_dym}/system/controlDict
{dir_main_dym}/system/decomposeParDictの"numberOfSubdomains 72;"も編集する必要があるはず。

pimpleDyMFoaMソルバーを実行する
#!/bin/bash

. qfoam_header.sh

# Solve flowfield
cd ${dir_main_dym}
mpirun -np $ncpu pimpleDyMFoam -parallel > log_pimpleDyMFoam 2>&1
しばらく待ちます・・5秒間の計算に実時間50時間程度かかった。

可視化

計算結果からVTKファイルを作成する
#!/bin/bash

. qfoam_header.sh

cd ${dir_main_dym}

#reconstructPar -latestTime > log_reconstructPar_post 2>&1
reconstructPar > log_reconstructPar_post 2>&1
foamToVTK > log_foamToVTK 2>&1

ParaView 5.6.0で可視化
・t=3.0sのときの速度ベクトル分布


空気力と迎角の時間履歴

ピッチングモーメント・迎角いずれも収束に向かっており、動的不安定性は出なかった。

おわりに

計算モデルとしてそれなりにロバストな印象はあるが、もう少し高速化できればいいかもしれない。
なお、いろいろ調べていくとやっぱり乱流モデルの影響は無視できないっぽい。

参考文献等

参考にさせていただきました:

{dir_main_dym}/constant/dynamicMeshDict


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

dynamicFvMesh      dynamicMotionSolverFvMesh;

motionSolverLibs ("libsixDoFRigidBodyMotion.so");

motionSolver    sixDoFRigidBodyMotion;

patches         (wall);
innerDistance   0.3;
outerDistance   5;

mass            15.6;
centreOfMass    (0.4 0 0);
momentOfInertia (1 0.5 0.5);
orientation
(
    1 0 0
    0 1 0
    0 0 1
);
angularMomentum (0 0 -0.1);
g               (0 0 0);
rho             rhoInf;
rhoInf          1;
report          on;

solver
{
    type symplectic;
}

constraints
{
    yLine
    {
        sixDoFRigidBodyMotionConstraint line;
        centreOfRotation    (0.4 0 0);
        direction           (0 1 0);
    }

    zLine
    {
        sixDoFRigidBodyMotionConstraint line;
        centreOfRotation    (0.4 0 0);
        direction           (0 0 1);
    }

    zAxis
    {
        sixDoFRigidBodyMotionConstraint axis;
        axis                (0 0 1);
    }
}

restraints
{
    verticalSpring
    {
        sixDoFRigidBodyMotionRestraint linearSpring;

        anchor          (0.4 0 0);
        refAttachmentPt (0.4 0 0);
        stiffness       4000;
        damping         2;
        restLength      0;
    }

    axialSpring
    {
        sixDoFRigidBodyMotionRestraint linearAxialAngularSpring;

        axis            (0 0 1);
        stiffness       0.0;
        damping         0.0;
        referenceOrientation $orientation;
    }
}


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

{dir_main_dym}/constant/turbulenceProperties

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

//simulationType  RAS;
simulationType laminar;

RAS
{
    RASModel        kOmegaSST;

    turbulence      on;

    printCoeffs     on;
}


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

{dir_main_dym}/constant/transportProperties

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

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

transportModel Newtonian;

nu              [0 2 -1 0 0 0 0] 1e-05;

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