いわて駐在研究日誌2。

NEVER STAND BEHIND ME

OpenFOAMで磁界解析

 先日に引き続き、とある研究ネタ(?)の仕込み用に電磁界解析をやる必要があり、オープンソースな解析ソフトについて、いろいろ調べてみた。2次元解析に限れば、オープンソースでないものであれば有限要素法ベースのSATEとかあるし、オープンソースであればfemmなどが公開されており、自由に使えるようだ。ちょっと使ってみたが、正直な所、とっつきにくい。昔のいわゆる「研究室内製ソフトウェアにGUIつけてみました」といった感じである。

 そういえばOFにもあったよなーということで、今回は、OF(-dev)にある電磁界解析について調査してみた。

 

1.OF-devのソルバー

solvers/electromagneticsを確認したところ、以下の3つが存在

  • electrostaticFoam: 静電界ソルバー
  • magneticFoam:静磁界ソルバー(永久磁石による磁界解析)

  • mhdFoam:磁性流体ソルバー( incompressible, laminar flow of a conducting fluid under the influence of a magnetic field.)

自分の用途としてはまずは静磁界解析ができれば良いので、magneticFoamのチュートリアルをみた所、なぜか存在せず、以下のCFD onlineのスレッドで公開されてたチュートリアルを実行してみる。

magneticFoam -- CFD Online Discussion Forums

 1)ダウンロード+展開

wget https://www.cfd-online.com/Forums/attachments/openfoam/10560d1324730520-magneticfoam-magsandbox.tar.gz

mv 10560d1324730520-magneticfoam-magsandbox.tar.gz magsandbox.tar.gz

tar zxvf magsandbox.tar.gz

cd magsandbox

 2)実行

./Allrun →ポテンシャルΨ、B、Hの計算

 

u-blox M8P評価ボード

とある目的で、u-blox M8Pを搭載した評価ボードキットC94-M8P-4-11(日本向け)を購入。これは、基地局およびローバー(移動局)機能を持ったNEO-M8P-2が2セット入っていてアクティブアンテナおよび局間通信のためのUHF無線もついており、いろいろ遊べそうな評価ボードです。

C94-M8P | u-blox

キットの内容

  • 2 × アプリケーション・ボード(どちらもNEO‑M8P‑2)
  • 2 × 外部UHFアンテナ(Japan 920MHz)
  • 2 × 外部アクティブGNSSアンテナ
  • 2 × マイクロUSBケーブル

設定ソフトウェアの導入と接続

①u-center for windowsのインストール

モジュールの設定には、ドライバも兼ねたu-center for windowsをu-bloxのサイトから落としてきてインストールしておく。現在のバージョンはv8.27で、このほかにAndroid版もあるようだが、Linux版はなさそうなのでWindowsマシン必須のよう。

https://www.u-blox.com/en/evaluation-software-and-tools

インストーラが起動してWizardが立ち上がるが、Windows10では、driverはUse Windows Serial Driverを選択しておく。

 

②接続

付属のmicroUSB-Bのケーブルで、評価ボードとPCを接続する。ここでは、PC1台に2つのボードを接続する。

 

Windows10のデバイスマネージャで確認すると、自動的にセンサー>u-blox GNSS location sensorが2つ認識される...が、これでは仮想COMポートを使って通信できないという謎仕様のため、それぞれデバイスマネージャから>ドライバソフトウェアの更新>コンピュータを参照・・・>デバイス一覧から選択>互換性のあるハードウェアを表示>USB Serial device を選んで、ドライバを変更する必要がある。

バイスマネージャでポート(COMとLPT)にUSBシリアルデバイスが2つ出てくればOK(今回はCOM3,4で2つ認識されたので、COM3を基地局、COM4を移動局とする)

セットアップ

u-bloxのサイトから、ユーザーズガイドなどのPDFをダウンロードできる。

Product resources | u-blox

ここでは、Setup Guideを見ながら確認していく(以下の設定については後日更新予定)。

 

③u-centerからの基準基地局"Base station"の設定

u-centerを起動し、上述のように基地局とするCOM3につながったボードをメニュー>Reciever>Port>COM3を選んで接続する。とりあえずbaude rateはdefaultの9600bpsで。接続が確立すれば、緑のアイコンになるはず。

u-centerメニュー>View>Configuration Viewより以下を変更する。

(1)UBX-CFG-TMODE3 ・・・Time Mode 3

Mode: 0-Disabled/1-Survey-in/2-Fixed Modeから選択する。

基地局として用いる場合は、1-Survey-in/2-Fixed Modeから選択する。

〇前者の場合

Minimum Observation Time =300 s

Position Accuracy = 2 m

などとすると、10-15分で1m程度の精度でfixするようだ。

UBX-NAV-SVINのmessageでも確認できる。

〇後者の場合

Use Lat/Lon/Alt Position = ON(基地局の座標を緯度・経度・高度で指定)

Lat 基地局の緯度を設定(Google Mapなど参照)

Long 基地局の経度を設定

Alt 基地局の高度を設定

Accuracy 0.01m

として、あらかじめ分かっているデータを入力する。

(2) UBX-CFG-PRTでRTCM3形式データがUART1(無線に接続?)のみ出力される設定とする(移動局との通信データ)

Target = 1-UART1

Protocol in = none

Protocol out = 5-RTCM3

Baudrate   = 19200

(3) UBX-CFG-MSGで、1005/1077/1087/1230の4つのメッセージでUART1に出力されるように選択する

※Setup Guideで1230のメッセージで10が入力されているが1の間違いと思われる。

 

④u-centerからの移動局"Rover"の設定

Rover側のボードは無線でRTCM messageを受け取ると、自動的にRTKモードになるらしいので、radio linkの設定だけすれば良いらしい。

RTKLIBを用いなくてもよいということかな?

(1) UBX-CFG-PRTでRTCM3形式データがUART1(無線に接続?)のみ受信される設定とする(移動局との通信データ)

Target = 1-UART1

Protocol in = 5-RTCM3

Protocol out = none

Baudrate   = 19200 (基地局と合わせる)

RTKモードに切り替われば、Rover側のボードのGreenLEDが点灯するはず。点滅でFLOAT、点灯がFIXEDということらしい....一応GreenGreenLEDの点滅までは各h人できたが、点灯状態までいかない。

cfMesh + OF5.x

安心して使えないsHMをやめてcfMeshでメッシュ生成にトライしようとしたが、OF5.x/devだとビルドに失敗するので、ちょっと調査。以下、OF5.xでの検証だが、-devでもいけると思う。

 

まず、ググったところ、以下のとおり

Compiling cfMesh with OpenFOAM 4.x -- CFD Online Discussion Forums

で、cfMeshのソースを変更しないといけないようだ。変更点は3つ。

 

① metaDict_ not declared. Added "IOdictionary.H" to meshLibrary/utilities/meshes/polyMeshGen/polyMeshGen.H

#ifndef polyMeshGen_H
#define polyMeshGen_H

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

#include "polyMeshGenCells.H"
#include "dictionary.H"
#include "IOdictionary.H"

namespace Foam
{

 ② edgeMesh has been incorporated into meshTools in OF5 and dev. I did an in-place replacement of -ledgeMesh with nothing

cfMeshのフォルダに散らばって存在しているMake/optionsで-ledgeMeshを削除すれば良いということかな?  ということで、findとsed -i で対処。

find . -name "options" | xargs sed -i s/"-ledgeMesh"/""/g

 

③ Pstream::blocking and Pstream::scheduled have been replaced by Pstream::commsTypes::blocking and Pstream::commsTypes::scheduled - I replaced these

これらもいっぱいあるので、findとsed -iで対応。

 

ここまでやって、tetMeshOptimisation.C, triSurface2DCheck.C, quadricFittingI.H のビルドで、なんかよくわからないエラーが...。eigenVector関数の引数が違うとのこと。

Foam::vector Foam::eigenVector(const tensor&, Foam::scalar, const vector&, const vector&)

vector normal1 = eigenVector(nt, ev[1]);

 :

該当するソースを見ると、

nt:  symmTensor

ev:  eigenVectors of "nt" ・・・実際の引数はev[1]などの成分(スカラー)を渡している。

ということのようなので、OFのsrc/OpenFOAM/primitives/Tensor/tensor/tensor.Cで定義されているeigenValue()の実装が4.xまでと変わったせいのようだ。以下のように修正してみる。

const vector ev = eigenValues(nt);
const vector Ux(1, 0, 0), Uy(0, 1, 0), Uz(0, 0, 1);

//- make sure the point stays on the surface
vector disp = simplex.centrePoint() - points[nodeI];

if( mag(ev[2]) > (mag(ev[1]) + mag(ev[0])) )
{
//- ordinary surface vertex
# ifdef OpenCFDSpecific
vector normal = eigenVectors(nt, ev).z();
# else
vector normal = eigenVector(nt, ev[2], Ux, Uy);
# endif

normal /= (mag(normal)+VSMALL);
disp -= (disp & normal) * normal;
}
else if( mag(ev[1]) > 0.5 * (mag(ev[2]) + mag(ev[0])) )
{
//- this vertex is on an edge
# ifdef OpenCFDSpecific
vector normal1 = eigenVectors(nt, ev).y();
# else
vector normal1 = eigenVector(nt, ev[1], Uz, Ux);
# endif

normal1 /= (mag(normal1)+VSMALL);

# ifdef OpenCFDSpecific
vector normal2 = eigenVectors(nt, ev).z();
# else
vector normal2 = eigenVector(nt, ev[2], Ux, Uy);
# endif

bunnyOctreeのチュートリアルとかやってみたけど、ぱっと見だいじょうぶそう。

 

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

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