いわて駐在研究日誌2。

NEVER STAND BEHIND ME

includeAngleとresolveFeatureAngleについて

いまさらながらsHMねたです。

 

STLで出力した図形のresolveFeatureAngleをうまく設定するには、paraviewのFeature edges filterを使うと、角度毎に表示させて確認ができる。この値をsHMDictで指定すれば良い。

 

surfaceFeatureExtractDictに指定するincludeAngleは、

includeAngle⁼180 - resolveFeatureAngle

でOK。

あらかじめ境界メッシュの動きが分かっている場合のdynamicMesh

境界パッチの動きが時間や空間の関数として分かっている場合にどうするかを調べたい。

 

境界パッチ形状の移動に伴う流体メッシュ変形については、例によってPenguintsさんのサイトにまとめられていた。ただし、パッチ形状が変形するわけではなく、あくまで剛体運動する場合のようである。

PENGUINITIS - 移動メッシュ

 

constant/dynamicMeshDictには以下の例の通り記述する。dynamicFvMeshには、dynamicMotionSolverFvMeshを指定する。

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolverLibs ("libfvMotionSolvers.so");

solver velocityComponentLaplacian;

velocityComponentLaplacianCoeffs
{
component x;
diffusivity inverseDistance (move);
}

 ここで、solverの指定と、指定solver毎のCoeffsサブディクショナリの組み合わせでパッチ移動に伴うメッシュ変形計算を行う。solverの指定により、境界パッチの移動を記述するファイルが異なる(0/pointMotionUx やpointMotionU, pointDisplacementなど)。

 

4.xで選択可能なsolverは以下の7つで、motionSolverLibs ("libsixDoFRigidBodyMotion.so")と指定すれば、solver で sixDoFRigidBodyMotionも指定できるらしい。

  1. displacementComponentLaplacian
  2. displacementInterpolation
  3. displacementLaplacian
  4. displacementLayeredMotion
  5. displacementSBRStress
  6. velocityComponentLaplacian
  7. velocityLaplacian

この中で、パッチ移動の変位(displacement)を指定するのは、1~5、速度を指定するのは6,7ということである。Interpolationは、メッシュ移動前後の変形から補間するもので、Laplacianはメッシュ移動のラプラス方程式を解くので、より滑らかなメッシュになる。

 

今回の用途では、直接変位が計算できるので、移動パッチ"wing"に対するものであれば、以下のようになるのかな。

solver displacementLaplacian;

displacementLaplacianCoeffs
{
diffusivity quadratic inverseDistance ( wing );
}

 移動変位は 0/pointDisplacement で指定することになり、

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

dimensions [0 1 0 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
wing
{
type calculated;
value uniform (0 0 0);
}

[snip]
}
// ************************************************************************* //

 ここで指定するtypeには様々な変位指定が可能である。(追記予定)

mesh tutorial (OF4.x)

ということで、これまでmeshの剛体回転ばかりやってきたので、少し目先を変えてmesh morphing関係の調査。まずはチュートリアルから。

① tutorials/incompressible/pimpleDyMFoam/movingCone

円筒領域(くさび型で2D化)の中をconeが平行移動する。メッシュ変形は以下のdynamicMeshDictに記述の通り、移動する部分速度を指定して、それに応じてメッシュが拡大・縮小する。

メッシュ数の増減はないタイプ。

dynamicFvMesh   dynamicMotionSolverFvMesh;
motionSolverLibs ( "libfvMotionSolvers.so" );
solver          velocityComponentLaplacian;
velocityComponentLaplacianCoeffs
{
    component       x;
    diffusivity     directional (1 200 0);
}

移動する部分(patch "movingWall")の速度の指定は、0/pointMotionUx でx軸成分の速度を指定している。

boundaryField
{
    movingWall
    {
        type            uniformFixedValue;
        uniformValue    constant 1;
    }
    farFieldMoving
    {
        type            slip;
    }
    fixedWall
    {
        type            uniformFixedValue;
        uniformValue    constant 0;
    }
    left
    {
        type            uniformFixedValue;
        uniformValue    constant 0;
    }

② tutorials/incompressible/pimpleDyMFoam/wingMotion

翼型(patch "wing")が6DOFベースで運動方程式に基づいた運動を行い、それに伴ってメッシュが変形する。同じフォルダに3つケースがあるが、sHMでメッシュを作り、simpleFoamで初期値を作成し、mapFieldでpimpleDyMFoamで実際にメッシュ移動を計算しているという流れ。

constant/dynamicMeshDictを見てみると、運動方程式を解くのに必要な質量、重心、慣性モーメント、重力などを指定しているほか、constraintsで拘束条件、restraintsで制約条件を追加しており、かなり複雑な移動をシミュレートできそう。このあたり、詳しい解説がどこかにないだろうか。

メッシュの変形結果については、各時刻フォルダにpointDisplacementが生成されるほか、polyMeshフォルダも作成される。

dynamicFvMesh dynamicMotionSolverFvMesh;

motionSolverLibs ("libsixDoFRigidBodyMotion.so");


solver            sixDoFRigidBodyMotion;

sixDoFRigidBodyMotionCoeffs
{
    patches         (wing);
    innerDistance   0.3;
    outerDistance   1;

    mass            22.9;
    centreOfMass    (0.4974612746 -0.01671895744 0.125);
    momentOfInertia (1.958864357 3.920839234 2.057121362);
    orientation
    (
        0.9953705935 0.09611129781 0
        -0.09611129781 0.9953705935 0
        0 0 1
    );
    angularMomentum (0 0 -2);
    g               (0 -9.81 0);
    rho             rhoInf;
    rhoInf          1;
    report          on;

    solver
    {
        type symplectic;
    }

    constraints
    {
        yLine
        {
            sixDoFRigidBodyMotionConstraint line;
            centreOfRotation    (0.25 0.007 0.125);
            direction           (0 1 0);
        }

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

    restraints
    {
        verticalSpring
        {
            sixDoFRigidBodyMotionRestraint linearSpring;

            anchor          (0.25 0.007 0.125);
            refAttachmentPt (0.25 0.007 0.125);
            stiffness       4000;
            damping         2;
            restLength      0;
        }

        axialSpring
        {
            sixDoFRigidBodyMotionRestraint linearAxialAngularSpring;

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

③tutorials/mesh/moveDynamicMesh/SnakeRiverCanyon

blockMeshで作成したベースメッシュのパッチminZを、moveDynamichMeshコマンドを用いてtriSurface/AcrossRiver.stl(地形)に合うように変形させるチュートリアル

ベースメッシュは、stl地形データよりもちょっとz方向上に生成されている。

motionSolverLibs ("libfvMotionSolvers.so");

solver displacementSBRStress;   //displacementLaplacian;
//solver velocityComponentLaplacian z;

displacementSBRStressCoeffs
{
    // diffusivity  uniform;
    // diffusivity  directional (1 200 0);
    // diffusivity  motionDirectional (1 1000 0);
    // diffusivity  file motionDiffusivity;
    diffusivity  quadratic inverseDistance 1(minZ);
}

 実際のメッシュ変形は、0/pointDisplacementに記載。

上面のmaxZパッチは移動なし、minX, minYなどの側面はそれぞれの法線方向ベクトルnに対し、fixedNormalSlipと指定して、nに直交方向する方向のみ許可している。

実際に変形するminZパッチ面については、指定速度で指定stlに向かって変形していく、すなわち時間変化する。プロジェクションの設定から単純にz軸方向のみのプロジェクションを行っていることがわかる。

boundaryField
{
    maxZ
    {
        type            fixedValue;
        value           uniform (0 0 0);
    }

    minZ
    {
        type            surfaceDisplacement;
        value           uniform (0 0 0);

        // Clip displacement to surface by max deltaT*velocity.
        velocity            (10 10 10);

        geometry
        {
            AcrossRiver.stl
            {
                type triSurfaceMesh;
            }
        };

        // Find projection with surface:
        //     fixedNormal : intersections along prespecified direction
        //     pointNormal : intersections along current pointNormal of patch
        //     nearest     : nearest point on surface
        // Other
        projectMode fixedNormal;

        // if fixedNormal : normal
        projectDirection (0 0 1);

        //- -1 or component to knock out before doing projection
        wedgePlane      -1;

        //- Points that should remain fixed
        //frozenPointsZone fixedPointsZone;
    }

[snip]

    maxX
    {
        type            fixedNormalSlip;
        n               (1 0 0);
    }

[snip]
}

 ④tutorials/mesh/refineMesh/refineFieldDirs

refineMeshコマンドの使い方らしいので、ちょっと違うかも。

OF-devのビルドエラー

学生に教える都合もあって、基本的に4.xの系列しか入れてないのだけど、岩手大の学生さんからの相談対応で、meshMorphing機能とかIB法関連の実装を調べるため、devとfoam-extend4.0を入れてみる。

 

extend4.0は問題なかったが、devがビルドエラーになり、/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/にあるtetPoints.Hがないと言われる。数日前のコミットでtetPoints.HがtetPointRef.Hに変わったらしい。

 

とりあえず、

cd ~/OpenFOAM/OpenFOAM-dev/src/OpenFOAM/meshes/primitiveShapes/tetrahedron

ln -sf tetPointRef.H tetPoints.H

 で逃げておけばいいのかな。

新年度です。

結局、動的ライブラリが放置状態なので今年度こそはと思う。

 

今年度の目標は、

校務

  • 授業は確実に準備をすすめること(前倒しでやること)
  • 卒研・特研は学生との作業時間を増やすこと
  • 地域連携を進め、共同研究3件を目指す(1件決定済み)

研究

  • マイクロORCのコア部分(タービン+ポンプ)を中心に開発をすすめる
  • マイクロ水平軸風車システム開発に区切りをつける(製品化の手前まで)
  • 最適化関係+OpenCAE(OpenFOAM/Salome_meca)のアウトプットを増やす
  • 生物規範工学の勉強をする
  • MAVなどの制御システムやIoT技術について勉強する
  • 論文書く

家庭&自己研鑽

  • 夏ごろに子供が2名体制になるので、子育ても頑張る
  • 英語でアウトプットできるようにTOEICの勉強を再開する+英語論文を書く

といったところでしょうか。

 

動的なsolidBodyMotionFvMeshのライブラリ作成(その2)

基礎にもどって、dynamicMeshの調査。

 

OF4.xでのdynamicメッシュは大きく分けると、

  • dynamicFvMesh ・・・Automatic mesh motion
  • topoChangerFvMesh・・・Topological mesh changes
に大別される(Ref.: A.O.Gonzalez,Mesh motion alternatives in OpenFOAM,2009)。
参考文献はちょっと古いが、現状の4.xでも、以下の通り。
< dynamicFvMesh >
  • dynamicInkJetFvMesh   relatively small mesh motion like ink-jet pumping
  • dynamicMotionSolverFvMesh  relatively small mesh motion
  • dynamicRefineFvMesh    refinement and corcening mesh
  • solidBodyMotionFvMesh    solid body motion
  • staticFvMesh    no motion (static)
  • (dynamicFvMesh)  just Abstract base class
< topoChangerFvMesh >
  • linearValveFvMesh  sliding meshes of rlative linear motions
  • linearValveLayersFvMesh   linearValveFvMesh+layer addtion/removal functions
  • mixerFvMesh   rotating slider mesh
  • movingConeTopoFvMesh  layer addtion/removal depending on cell layer thickness
  • rawTopoChangerFvMesh topoChangerFvMesh without any added functionality
  • (topoChangerFvMesh)  just Abstract base class
 
今回の自分の用途には、やはりsolidBodyMotionFvMesh が一番近い機能を持ってそうなので(topological change機能は不要)、これをベースに開発してみよう。

directDrivenMotionFvMeshの作成(OF 4.x)

参考とするsolidBodyMotionFvMeshは、dynamicFvMeshに継承されているので、solidBodyMotionFvMesht同じレベルでdirectDrivenMotionFvMeshというライブラリを作成してみる(Chen et al.と同じ名前だとまずいかなということで)。Motionを具体的に記述するmotionFunctionも作成してみよう。

 

(1)まずは、solidBodyMotionFvMeshをcopyしてフォルダ名やファイル名をrename

$ cp -r $FOAM_SRC/dynamicFvMesh/solidBodyMotionFvMesh/ .

$ mv solidBodyMotionFvMesh directDrivenMotionFvMesh

 

$ cd  directDrivenMotionFvMesh/

$ mv solidBodyMotionFvMesh.C directDrivenMotionFvMesh.C

$ mv solidBodyMotionFvMesh.H directDrivenMotionFvMesh.H

$ mv multiSolidBodyMotionFvMesh.C multiDirectDrivenMotionFvMesh.C

$ mv multiSolidBodyMotionFvMesh.H multiDirectDrivenMotionFvMesh.H

 

$ mv  solidBodyMotionFunctions directDrivenMotionFunctions

(2)sedコマンド(-iオプション)を利用したファイル中のclass名などのrename

$ sed -i s/solidBodyMotionFvMesh/directDrivenMotionFvMesh/g directDrivenMotionFvMesh.C

$ sed -i s/solidBodyMotionFvMesh/directDrivenMotionFvMesh/g directDrivenMotionFvMesh.H

 

$ sed -i s/multiSolidBodyMotionFvMesh/multiDirectDrivenMotionFvMesh/g multiDirectDrivenMotionFvMesh.C

$ sed -i s/multiSolidBodyMotionFvMesh/multiDirectDrivenMotionFvMesh/g multiDirectDrivenMotionFvMesh.H

$ sed -i s/solidBodyMotionFvMesh/directDrivenMotionFvMesh/g multiDirectDrivenMotionFvMesh.C

$ sed -i s/solidBodyMotionFvMesh/directDrivenMotionFvMesh/g multiDirectDrivenMotionFvMesh.H

 

(solidBodyMotionFunction関連)

$ sed -i s/SBMF/DDMF/g directDrivenMotionFvMesh.C

$ sed -i s/SBMF/DDMF/g directDrivenMotionFvMesh.H

$ sed -i s/SBMF/DDMF/g multiDirectDrivenMotionFvMesh.C

$ sed -i s/SBMF/DDMF/g multiDirectDrivenMotionFvMesh.H

 

$ sed -i s/solidBodyMotionFunction/directDrivenMotionFunction/g directDrivenMotionFvMesh.C

$ sed -i s/solidBodyMotionFunction/directDrivenMotionFunction/g directDrivenMotionFvMesh.H

$ sed -i s/solidBodyMotionFunction/directDrivenMotionFunction/g multiDirectDrivenMotionFvMesh.C

$ sed -i s/solidBodyMotionFunction/directDrivenMotionFunction/g multiDirectDrivenMotionFvMesh.H

※ これ以降、multiDirectDorivenMotion.C/Hについては詳述しません。後日、追記予定。

(3)pointPatchFields/derived/solidBodyMotionDisplacement関連のrname

$ pointPatchFields/derived

$ mv solidBodyMotionDisplacement/ directDrivenMotionDisplacement/

$ cd directDrivenMotionDisplacement/

$ mv solidBodyMotionDisplacementPointPatchVectorField.C directDrivenMotionDisplacementPointPatchVectorField.C

$ mv solidBodyMotionDisplacementPointPatchVectorField.H directDrivenMotionDisplacementPointPatchVectorField.H

 

$ sed -i s/solidBodyMotion/directDrivenMotion/g directDrivenMotionDisplacementPointPatchVectorField.C

$ sed -i s/solidBodyMotion/directDorivenMotion/g directDrivenMotionDisplacementPointPatchVectorField.H

$ sed -i s/SBMF/DDMF/g directDrivenMotionDisplacementPointPatchVectorField.C

$ sed -i s/SBMF/DDMF/g directDrivenMotionDisplacementPointPatchVectorField.H

(4)Makeフォルダの作成および修正

$ cp -r  $FOAM_SRC/dynamicFvMesh/Make .

Edit Make/files as:

directDrivenMotionFvMesh.C
multiDirectDrivenMotionFvMesh.C

directDrivenMotionFunctions/solidBodyMotionFunction/solidBodyMotionFunction.C
directDrivenMotionFunctions/solidBodyMotionFunction/solidBodyMotionFunctionNew.C

directDrivenMotionFunctions/directDrivenMotionFunction/directDrivenMotionFunction.C
directDrivenMotionFunctions/directDrivenMotionFunction/directDrivenMotionFunctionNew.C
directDrivenMotionFunctions/rotatingMotion/rotatingMotion.C

directDrivenMotionFunctions/axisRotationMotion/axisRotationMotion.C
directDrivenMotionFunctions/oscillatingRotatingMotion/oscillatingRotatingMotion.C

pointPatchFields/derived/directDrivenMotionDisplacement/directDrivenMotionDisplacementPointPatchVectorField.C

LIB = $(FOAM_USER_LIBBIN)/libdirectDrivenMotionFvMesh

(※ motionFunction関係は現時点で手つかずです。)

Edit Make/option as:

 EXE_INC = \
    -I$(LIB_SRC)/triSurface/lnInclude \
    -I$(LIB_SRC)/meshTools/lnInclude \
    -I$(LIB_SRC)/dynamicMesh/lnInclude \
    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
    -I$(LIB_SRC)/finiteVolume/lnInclude

LIB_LIBS = \
    -ltriSurface \
    -lmeshTools \
    -ldynamicMesh \
    -ldynamicFvMesh \
    -lfiniteVolume

 (5)directDrivenMotionFunctionの作成および修正

$ cd directDrivenMotionFunctions

(使わなそうなmotionFunctionをいくつか削除)

$ rm -rf SDA  linearMotion  multiMotion  oscillatingLinearMotion tabulated6DoFMotion

$ ls
axisRotationMotion  oscillatingRotatingMotion  rotatingMotion  solidBodyMotionFunction(solidBodyMotionFunctionをコピー&リネーム)

$ cp -r solidBodyMotionFunction directDrivenMotionFunction

$ cd directDrivenMotionFunction

$ mv solidBodyMotionFunction.C directDrivenMotionFunction.C

$ mv solidBodyMotionFunction.H directDrivenMotionFunction.H

$ mv solidBodyMotionFunctionNew.C directDrivenMotionFunctionNew.C

(directDrivenMotionFunction.Cなどの修正)

$ sed -i s/solidBodyMotion/directDrivenMotion/g directDrivenMotionFunction.C

$ sed -i s/solidBodyMotion/directDrivenMotion/g directDrivenMotionFunction.H

$ sed -i s/solidBodyMotion/directDrivenMotion/g directDrivenMotionFunctionNew.C

$ sed -i s/SBMF/DDMF/g directDrivenMotionFunction.C

$ sed -i s/SBMF/DDMF/g directDrivenMotionFunction.H

$ sed -i s/SBMF/DDMF/g directDrivenMotionFunctionNew.C

 (6)wmake libso

とりあえず、ここまでの変更を反映してライブラリが作れるかwmake libsoしてみる。

$ cd ../../   (ライブラリフォルダのトップに戻って)

$ wclean

$ wmake libso

[snip]

'/home/waku/OpenFOAM/waku-4.x/platforms/linux64GccDPInt32Opt/lib/libdirectDrivenMotionFvMesh.so' is up to date.

 以上の変更は、ライブラリやクラス名称を変更しただけなので、中身はオリジナルのもののまま。これをいろいろ変更していく。ちなみにライブラリのビルドは通るが、motionFunctionは多重定義となっているので、実行時にエラーになるので、こちらも名前を変える必要あり。

動的なsolidBodyMotionFvMeshのライブラリ作成

OF、あるいは、動的メッシュ機能をもつたいていのCFDソフトは、メッシュの移動・回転速度を指定して計算するのが一般的だろう。

 

風車の研究をしているので、ブレードトルク→回転の運動方程式→角速度→回転角のような動的なメッシュ移動ができないか試行錯誤しているが、自分のC++の知識ではOFのC++の高レベルの利用に追いついていけず、遅々として進んでいない(一応それらしい動きができるようにしてみたが、回転の運動方程式をfunctionObject機能を利用して外部で解いているので、次のステップとして、できれば全部OFの中でクローズしたいと考えている)。

 

もっと時間をかけてC++勉強しないとだめだなぁと思いつつ、いろいろ調べてたらそのようなライブラリを使用している論文を発見(Goong Chen et al, "OpenFOAM computation of interacting wind turbine flows and control(I): free rotating case", WPのpreprint?)。2014年かということで、自分がアイデア考えだしたころに論文になってたということになる。ちょっと残念。Chen先生は結構OF関連では、論文を出してるみたいだ。

 

彼らの論文によると、OF-extendを利用して(???明確な記述なし)、solidBodyMotionFvMeshからforceDrivenFvMeshというライブラリを作成し、その中で、forceDrivenMotionFunction( FDRotatingMoton, codedFDRotatingMotion)の回転計算functionを定義して利用しているようである。

 

ライブラリ自体のあまり詳しい情報は小論文にのっておらず、ネットでも検索でヒットしないので、情報を得るには直接コンタクトするか、CFDonlineに聞くしかないようだ。もう少し調査しつつ、正規版4.xでも作れないか作業してみたい。