連結浮体の波浪動揺解析方法に関する調査

背景

連結浮体の揺動という課題に対して、OpenFOAMでそのものズバリの解析例は記憶に無かったが、数年前に拙宅HP記事のAllrun(v1806標準チュートリアル)やってみた、の中で紹介した、

  • multiphase/overInterDyMFoam/boatAndPropeller

というチュートリアルに近い(rudderの部分をもう一つの船体(boat)に置き換えれば良さそう)」と考えられた。しかしこのチュートリアルは何故かv2012まで同梱されていたのが、現在では無くなってしまっている(⇒Allrunリスト)。しかも詳しく調べると、v1806の初出に比べ、v1812以降は形状も少し異なってしまっているという怪しい(?)チュートリアルであった。

それはともかく、このあたりを足掛かりに、最新のDEXCSを使って、どこまで出来るのか調べてみようとなった。

また、公開可能案件としてのご依頼であり、オーバーセットメッシュを使った解析であってもDEXCSを使うと、チャチャッとケースのセットアップができる好例…として公開できるかもしれないという目論見もあった。

仕様

  1. 正弦波の波周期・波長を任意に設定可能な境界条件
  2. 任意の数の浮体を,任意の指定位置で2点支持連結.
  3. 規則波中を一定の平均速度で前進する.(平均前進速度周りで,波の位相に応じた前後揺れが起きる)
  4. 向かい波(追い波)のみを対象.運動自由度は3自由度(surge, pitch, heave).
  5. 連結部に作用する荷重を出力できる.

予備調査

まずは、上述の標準チュートリアルケースを最新DEXCS(OpenFOAM-v2406)で動かすことが出来るかどうか? であったが、当然ながらすんなりとは動いてくれなかった。

しかしながら、動かなかった原因は、dynamicMrshDictの文法問題であり、OpenFOAMのバージョンによってパラメタオブジェクトの構造が異なっていたのが主な原因であった。古いパラメタファイル(dynamicMeshDict)を新しい書式に変更するだけで、すんなり動いてくれた。

因みに新しい書式で記述された例を理解するには、

  • multiphase/overInterDyMFoam/RAS/rigidBodyHull

が役に立ってくれた。

また、このチュートリアルは、前述のboatAndPropellerに比べて、rudderが無いとうだけで boat や propeller が実物に近い形状になっていた。rudderもpropellerもboatに対する連結体としては同じオブジェクトとみなせるので、結局のところ、boatAndPropeller が無くなったのは、これが rigidBodyHull に変わったから、と見て良さそうであった(rigidBodyHullはv2106〜)。

モデルの説明

モデルは当初、標準チュートリアルのboatAndPropellerをベースに動作確認してきたが、最終的にわかり易くするため極力シンプルな数字になるよう構成し直したもので公開することとした。

FreeCADで作成したCADモデルを示しておく。まずは全体外観である。

下の図も含めて寸法がmm(ミリメートル)で表示されているが、OpenFOAMで計算される際にはメートルで解釈される。プールは、断面が1メートル四方で5メートル長さの流路ということである。

奥の方にある赤い物体が連結された浮遊体で、以下に詳細図を示す。

長さ 1メートル、幅 0.15メートル、高さ 0.1メートルの直方体が、0.1メートルの隙間で連結されている。COMというのは重心位置である。

ケースファイル一式の説明と実行方法

ケースファイル一式は github で公開しておいた。ダウンロードして解凍すると、以下のように展開されるはずである。

ケースファイル一式は3つのフォルダに区分けして収納してある。それぞれのフォルダ名から推察していただけると思うが、ケースファイルの用途が異なっている。

Caseフォルダに収納したAllrunを実行すれば、全ケースファイルを実行、Allclean にて全初期化するようになっている。但し、OpenFOAM-v2406を実行できる環境(端末)にて。OpenFOAMの他のバージョンについては未検証。

Caseフォルダ下の個々のケースファイルにも、Allrun、Allclean が収納してあり、個別に実行することも可能。但し、これらを実行するには、Meshフォルダ下で必要なメッシュファイル(必要なメッシュファイルは、同梱のAllrun,preを見ればわかるようになっている)が作成済であることが前提である。

Meshフォルダ

Meshフォルダ下のケースファイルは、パーツ毎にメッシュを作成するフォルダであり、同梱のFreeCADモデル(.FCStd)のメッシュコンテナを使って作成済のメッシュが同梱されている。DEXCS-WBが使える環境であれば、細分化パラメタを変更したり、モデルそのものを変更することも可能になっている。

各ケースファイルのフォルダ名だけを残して、市販ツールなど別の手段で作成したメッシュケースで置き換えて使っても良い。但し、パッチ名とパッチタイプが整合していることは必須条件である。

Caseフォルダ

Caseフォルダに収納されたケースファイルは、全体の目標仕様に対して個別の仕様技術を検証できるようにしたケースファイルとして収納したものである。

水面高さ、波の条件や、船体の剛体運動条件などは、パラメタ変更して計算可能になっているが、本ケースでは簡単にする為、乱流計算は想定していない。

Refフォルダ

Refフォルダは標準チュートリアルケースであるが、Allrunを実行した後、TreeFoamによる初期化を施したもの(計算結果は削除するが、メッシュは残った状態)である。上述のCaseフォルダは、これらのチュートリアルケースを雛形として修正変更したものである。具体的には、

  • Ref/rigidBodyHull/background ⇒ Case/1_hullInPool
  • Ref/stokesI ⇒ Case/2_wavePool

という関係になっている。これらをケース比較することにより、変更箇所を明示できるようにした。また、変更していない箇所については、筆者にも十分理解できていない点があることはお断りしておく。

プールに浮体を置いてみた

まずは、Ref/rigidBodyHull/backgroundのケースファイルを参考にして、浮体だけをプールに置くというケースファイルを作成した。

ここでは、オーバーセットメッシュを使ったケースファイルのセットアップと、dynamicMeshDict ファイルのセットアップ要領がポイントになっている。

オーバーセット計算のセットアップ要領

セットアップの要点は、Allrun.preに集約されている。

#!/bin/sh
cd ${0%/*} || exit 1                        # Run from this directory
. $WM_PROJECT_DIR/bin/tools/RunFunctions    # Tutorial run functions

cp -r ../../Mesh/background/constant/polyMesh ./constant/
runApplication mergeMeshes . ../../Mesh/hull-1 -overwrite
runApplication checkMesh

restore0Dir

runApplication setFields
runApplication topoSet -dict system/topoSetDict.cHullProp

オーバーセット計算では、いわゆるバックグラウンドと称される全体メッシュ(5行目でコピーしたメッシュ)と、個別に詳細を調べたいオーバーセットメッシュとを、単純にmergeMeshesコマンドを使ってメッシュを合体させる(6行目)。但しそれだけでは不十分で、それぞれのメッシュ領域のセルに対して、zoneIDというスカラー番号を付与する必要がある。これを実行しているのが11行目のsetFieldsコマンドであるが、

system/setFieldsDict(17〜)
defaultFieldValues
(
    volScalarFieldValue zoneID 123
    volScalarFieldValue alpha.water 0
);

regions
(
    // Set cell values
    // (does zerogradient on boundaries)
    cellToCell
    {
        set region0;

        fieldValues
        (
            volScalarFieldValue zoneID 0
        );
    }

    cellToCell
    {
        set region1;

        fieldValues
        (
            volScalarFieldValue zoneID 1
        );
    }

// Set cell values
    // (does zerogradient on boundaries)
    boxToCell
    {
        box (-1000 -1000 -1000) (1000 1000 0.0);

        fieldValues
        (
            volScalarFieldValue alpha.water 1
        );
    }
);

において、バックグラウンドメッシュに対して zoneID 0(17行目)オーバーセットメッシュに対して zoneID 1 (27行目)を付与している。その際に、ベースとしたcellSet情報(13行目と23行目のregion#)はどうやって作ったか?となるが、これはAllrun.preの7行目checkMeshコマンドを実行したことにより作成されている。このあたり、標準チュートリアルケースではtopoSetコマンドを駆使してcellSet情報を作成しているが、そのようにせずともcheckMeshコマンドを使えば良いというのを今回発見した。

なお、オーバーセット計算の説明からは外れるが、31行目以降は、プールの水位を決めており、ボートが丁度半分沈んだ位置としてある。

また、オーバーセットメッシュ(Mesh/hull-1)については、DEXCS-WBを使って、以下のイメージ図に従って作成したケースファイルが同梱されている。かような単純形状であれば、blockMeshで作成することも可能であり、そうしても良かったが、そうすると船体(hullWall1)やオーバーセット領域(overset1)の形状や、メッシュ細分化方案の変更は容易でない。かような応用変更を前提としたケースファイルということである。

変更するに際しての留意点を以下に記しておく。

  • FreeCADモデルのコンポーネントツリー上に非表示オブジェクトがたくさんあるが、本モデルを形状変更する際に、全体のバックグラウンドモデル(inlet, outlet, など)との対比や干渉を調べることが出来るよう残しておいたものである。メッシュ作成時には、これらを必ず非表示状態に戻しておく。非表示パーツは削除しても問題は無い。
  • パーツの名前(hullWall1, overset1)は変更しない(変更しても良いが、計算ケースの境界条件名も変更する必要が生じて、面倒になる)。
  • パーツそのものを変更する場合、既存パーツを削除して同名の別形状パーツを新規作成するなり、インポートするなりで良いが、それだけでなくメッシュ細分化コンテナ(MeshRefinement, MeshRefinement001)で対象パーツを指定し直す必要がある。
  • 本例におけるメッシュ細分化コンテナの役割は境界タイプ(wallまたはoverset)を指定するのが主たる目的で、細分化レベルは「0」(細分化しない)を指定してあるが、数字を大きくしても良い。
  • ソルバーコンテナ(CfdSolver)は使用しない。

dynamicMesh計算のセットアップ要点

dynamicMesh計算は、constant/dynamicMeshDictファイルで動作が規定されるが、これの元になったのは標準チュートリアルケースrigidBodyHullで使われていた同名ファイルである。これをどう変更したのか少しでもわかり易くすべく、両者のDictファイルをDEXCSに同梱のKDiff3ツールを使って対比表示したイメージ図を使って説明する。左側が変更後、右側が変更前(標準チュートリアルケース)である。

まず、変更前では、solversブロック(19行目)の中に2つのsolver(21行目のVFと30行目のboat_propeller)を定義する構造となっていた。VFは、boat_propellerの動きに連動して(37, 38行目)バックグラウンドメッシュ(26行目の cellSet c0)全体を動かそうとするもので、今回の計算でこの機能は使っていない(後述するが、使おうとしたが使えなかった)。

ソルバーの名前は何でも良さそうなので、標準チュートリアルでboat_propeller(30行目)としていたところを、doubleHull(22行目)に変更した(といっても本ケースはsingleHullでしかないが…)。

標準チュートリアルケースのVFと連動させるパラメタ相当部分(37, 38行目)は不要になるので削除した(残しておいても問題は無さそうであった)。

左側の35, 36行目のaccelaration…パラメタの値が変更されているが、試行錯誤の結果がそのまま残してあるだけで、どちらが良いとかの意味は無い。違いもよくわかっていない。

左側の28行目

cellSet cHullProp;

にて、メッシュが動く領域を規定しているが、これはAllru.preの12行目を実行することで作成されている。因みに、system/topoSetDict.cHullProp は以下の内容で、

system/topoSetDict.cHullProp(17〜)
actions
(
    {
        name    cHullProp;
        type    cellSet;
        action  new;
        source  cellToCell;
        set     region1;
    }
);

cellSet の region1 を指定しているだけである。従って28行目を、

cellSet region1;

としても良かったが、後で複数の物体を繋ぐことを見越して、この行はそのまま残して使うこととした。

左側39行目(右側49行目)のbodies は剛体を定義するブロックで、複数の剛体を定義でき、右側標準チュートリアルではhull(51行目)とpropeller(83行目)を定義しているが、左側本ケースでは、hull-1(41行目)だけであり、右側hullブロックのパラメタ名をそのまま流用、パラメタの値を本ケースに適合した。

45行目のmassと48行目のinertiaはsurfaceInertiaというユーティリティコマンドを使って計算したものである。Mesh/hull-1 フォルダ下に、hullWall1のSTLデータも同梱してあるので、ここで以下のように入力する。

$ surfaceInertia -density 700 hullWall1.stl

-density は容易に推測できると思うが密度で、700(kg/m3)を指定したということは、水に浮かべた時に、7割方が水中に沈むということになる。以下のような結果が帰ってくるので、ここから数値を丸めて転記したものである。

(途中省略)

Entries for sixDoFRigidBodyDisplacement boundary condition:
mass            10.5000005737;
centreOfMass    (0.5 0 0);
momentOfInertia (0.0284375033791 0.883750048547 0.894687550448);
orientation     (1 -0 5.84081946706e-18 0 1 0 -5.84081946706e-18 -0 1);

Writing scaled principal axes at centre of mass of "hullWall1.stl" to "axes.obj"

End

なお、centerOfMassについては、上の数字とdynamicMeshDict中(47行目)の数字が違うが、これについては、KDiff3の比較イメージ図においてdynamicMeshDict中(50行目)表示が途中で切れているので、以下にこの部分だけ表示しておくが、

transform       (1 0 0 0 1 0 0 0 1)  (0.5 0 0);

この部分の数字(0.5 0 0 )と併せて解釈が必要である。

つまり、44行目で、

parent root;

とあり、全体座標系(root)でいうところの上記位置(0.5 0 0)が回転中心であることを指している。dynamicMeshDict中(47行目)のcenterOfMassについては、この回転中心からの相対座標で定義する。本例では中身が一様な物体を想定して回転中心と重心を合致させているということである。標準チュートリアルでは実際の船体を想定して回転中心からずらした値になっているようで、本例でも任意に変更して良い(補足3)。

55行目の joints ブロックも変更してあるが、今回の計算仕様として、

自由度の説明図

補足4

運動自由度は3自由度(surge, pitch, heave)ということであったから。標準チュートリアルでは、Sway(Py)も考えているということである。なお、Pxyz という書き方は許されても、Pxz という書き方は不可で、Px とPz に分けて定義する必要があった。

68行目の patches ブロックの変更内容は説明するまで無いであろう。

標準チュートリアルにおける83行目propeller を定義する部分は削除しているが、5_hullsInPool のケースで、hull-2に変更して使うことになる。

また、restraints (標準チュートリアル105行目)ブロックの中身(107〜118行)を空にしているが、この部分は船体に推力を付与しており、4_hullMoveInPool のケースで数値を適合させて使うことになる。

境界条件

TreeFoamのGridEditor を使って表示した境界条件を以下に示しておく。但し、Allrun もしくは、Allrun.pre を実行した後でないと表示できない点には注意されたい。

参考に標準チュートリアル(rigidBodyHull/background)でのGridEditor表示も以下に示しておく。但し、乱流変数については非表示としてある。

境界名は変わっているが、プールの外周部分を除いてほとんどの条件はコピペされたものになっている。

プールの外周部分について補足しておくと、back, base, front については簡単にする為対称境界とした。inlet, outlet については閉じた境界を想定しており、本来であれば、type wall とすべきであったが、他のケースでは全て何らかの流出を伴う条件で、type patch を使うことになる。層流計算の場合には type patch であっても問題無いので、他のケースと共用できるよう wallとはしていない。

その他諸々

KDiff3を使って、本ケースのケースファイル自体を、標準チュートリアル(rigidBodyHull/background)と比較したのが以下の図になる。左欄にツリー表示されたフォルダやファイル名に対して、四角の色付きマークが2つ並んでおり、左側が本ケースで右側が標準チュートリアルケースで、どちらも緑色であれば全く同一で、赤マークは内容が異なっており、黒マークはファイルそのものが欠落しているという意味である。

本ケースにおいて欠落した(黒マーク)ファイルについては、説明するまでないであろう。

これまでの説明の中で言及していないものについて、以下補足しておく。

system/forces, forces.dplt, moments.dplt

system/forces はファイルの名前から想像出来るかと思われる。浮体に加わる流体力を算出するファンクションブロックである。controlDictファイルから#include文で組み込まれる。

forces.dplt, moments.dplt は上記のポストプロセスデータDEXCSプロットツールを使ってグラフ表示させる為のパラメタファイルであり、それを使って作成したプロット図も以下に示しておく。

 

船体の密度設定(700kg/m3)からは7割方水面下になろうとする。水面の初期位置が船体高さの丁度半分の位置からの計算であるから、慣性により沈みすぎて、また浮かび上がるという、プカリプカリした状況を見てとれる。モーメントが増大するのは、プールの端面(inlet, otlet)での反射波の影響か。

system/fvschems

system/fvSchemsについては、標準チュートリアルのものをそのまま使うと、実行時に下記のエラーが表示されるので、変更が必要であるとわかる。

[0]
[0]
[0] --> FOAM FATAL IO ERROR: (openfoam-2406 patch=241212)
[0] Number of searchBoxDivisions 3 should  equal the number of zones 2
[0]
[0] file: system/fvSchemes/oversetInterpolation at line 67 to 77.
[0]
[0]     From virtual bool Foam::cellCellStencils::trackingInverseDistance::update()
[0]     in file cellCellStencil/trackingInverseDistance/trackingInverseDistanceCellCellStencil.C at line 593.
[0]
FOAM parallel run exiting
[0]

変更は、以下に示す1箇所(73行目)。//(コメントアウト)で残したのは、ケース 5_hullsInPool の連結浮体とした時に、復活させることになるからである。

プールで波を作ってみた

前のケース(1_hullInPool)のinlet, outlet は静止壁であったが、これを波が流出入する境界条件としたい。波が流出入する境界条件は、標準チュートリアルに多く存在するが、ベースソルバーはinterFoamにて動作するものであった(overInterDyMFoamで使用できるかどうかわからなかった)。

そこでまずは、今回の対象領域(プールだけで浮体の無い状態)で標準チュートリアルに準じた波浪状態を再現できるプールモデルを作成した。

DEXCS-WBを使って、今回のプールモデルを対象に、標準チュートリアルケース(interFoam/laminar/wave/stokesI)を雛形とした解析手順のイメージ図を以下に示しておく。

解析コンテナ(dexcsCfdAnalysis)のプロパティ画面で、Template Case を選択、右隅下の❶「…」部分をクリックすると、フォルダの選択画面が現れるので、標準チュートリアルケースの在所(/usr/lib/openfoam/…,/stokes または、 Ref/stokesI でも良い)を選択する。

メッシュコンテナ(CfdMesh)をダブルクリックすると「メッシュ作成タスク画面」に切り替わる、

このボタンを押すと、TemplateCaseのケースファイル一式(0, constant, system)をコピーし、メッシュ(cfMesh)作成用のfmsファイル、system/meshDict ファイルが作成される。

❸ケース作成に成功すると、このボタンが有効になるので、これをクリックしてメッシュ計算に進む。

このボタンを押すとファイルマネージャーが起動されるので、パラメタファイルの変更が必要な場合にはこれを使う。

本例の場合には、system/setFieldsDict ファイルの変更が必要。具体的には、水面位置を設定する箇所で、26行目の

box (0 0 0) (20.0 1.0 0.6);

を、以下のように変更する。

box (-1000 -1000 -1000) (1000 1000 0.0);

このボタンを押すと、タスク画面が閉じて、FreeCAD画面上段のDEXCSツールバー(本例における❼〜❾)が使えるようになる。

本モデルと標準チュートリアルでは、パッチ名が一致していないので、GridEditorを使って境界条件を整合させる。

具体的な整合(編集)方法は、次項で説明する。

OF-端末が起動したら、setField コマンドを実行する。

ソルバーコンテナのタスク画面が表示されるので、DEXCSチュートリアルでの説明通り、ボタンを順番に押していけば良い。

GridEditorを使って、境界条件を適合する作業イメージは以下に示しておく。

atomosphere, back, base, front の境界条件がブランクになっているところを、コピペ作業で埋めていくだけである。

Case/2_wavePool は以上の手順で作成したケースファイルから、不要ファイルを削除し、system/controlDict(計算出力パラメタを変更),  constant/waveProperty を適合し直したものである。

constant/waveProperties(17〜)
inlet
{
    alpha           alpha.water;

    waveModel       StokesI;

    nPaddle         1;

    waveHeight      0.2;

    waveAngle       0.0;

    rampTime        0.5;

    activeAbsorption yes;

    wavePeriod      2.0;
}

outlet
{
    alpha           alpha.water;

    waveModel       shallowWaterAbsorption;

    nPaddle         1;
}

(9行目)waveHeight, (13行目)rampTime, (17行目)wavePeriod といったあたりを変更した。

造波理論については浅学にて正しく説明できないが、英文解釈である程度推定できる。パラメタを変更してほぼ用語のイメージ通りに変化することは確認できた。

波のあるプールに浮体を浮かべてみた

前節にて、今回の対象とするプールで、波を発生する為の境界条件とパラメタファイル(constant/waveProperties)がわかったので、これをオーバーセット計算でも使えるのかどうかの確認である。

3_hullInWavePool は、1_hullInPool に対して、境界条件を以下のように変更、

constant/waveProperties を追加しただけである。但し、waveHeight は、0.2 ⇒ 0.05に変更した(オーバーセット領域がプールから大きくはみ出してしまうので)。

結果は、めでたしめでたし…であった。

浮体を動かしてみた

次に浮体を動かす方法である。当初の目標仕様は一定速度で動かすということであった。

これを実現するのに、当初は標準チュートリアルの rigidBodyHull にならって dynamicMeshDict のsolversブロックの中でメッシュそのものを動かせないかであった。つまり、標準チュートリアルでは、

solvers
{
    VF
    {
        motionSolverLibs (libfvMotionSolvers);
        motionSolver            solidBody;
        solidBodyMotionFunction drivenLinearMotion;
        cellSet                 c0;
        cOfGdisplacement        CofG;

    }

となっていた。この部分を以下のように改変する。

solvers
{
    VF
    {
        motionSolverLibs (libfvMotionSolvers);
        motionSolver            solidBody;
        solidBodyMotionFuncton linearMotion;
        velocity    (0.1 0 0);
    }

こうすると、確かにメッシュ全体が指定の速度で動くことは確認できた(linearMotionの由来は補足2を参照)。さらにcellSetを指定すれば、オーバーセット領域のみを動かすことも可能であった…が、しかしである。浮体の方がうまくついてきてくれない。浮体の回転中心が絶対座標系で固定されたままであるかの挙動であった。但しこの問題分析は後になってわかったもので、本調査は本記事で書いている順番通りにやったものではない。当初はこのメッシュを動かす方法も、回転中心や重心の指定方法の手探りと同時進行しながらやっていた為、これ以上の調査には相当な困難が予想された。

結局、この方法は諦めて、標準チュートリアルでやっているのと同じように、浮体に推力を付与して動かすこととした。具体的には、

constant/dynamicMeshDict(77〜)
restraints
{
    force
    {
        type        externalForce;
        body        hull-1;
        location    (0 -0.03 0);
        force       table
        (
            ( 0   (6 0 0))
            ( 2   (1.5 0 0))
            (20   (1.5 0 0))
        );
    }
}

とした。7行目 location(推力の作用点)は、絶対座標系で指定。force(推力)は時間で変化させる形式である。最初大きく加速して、あとは一定、というのが現実のアクセルワークであろう。推力の値そのものは、水の抵抗がわからなかったので試行錯誤(10秒間でプールから出ることのないように)で決めた。水の抵抗が無いとしたら、10秒間で6メートルほど進む事になるのだが、本ケースでは1.7メートルしか動いていない。どんなものであろう?

なお、定常速度になっているかどうかもわからない。これを実現するには、もっと長い流路が必要になりそうではある。

また、本計算の流出入境界(inlet, outlet)は静止壁条件に戻して実施している、というか 1_hullInPool のケースに対して上記変更を加えたものである。

浮体を繋いでみた

これも要因を単独で調べるべく、1_hullInPool のケースに対して、もうひとつの浮体を追加した。これも変更箇所をわかり易くするため。KDiff3によるケース比較図を示しておく。左側が本ケース(5_hullsInPool)、右側が 1_hullInPool である。

変更ファイルがたくさんあるが、ほとんど(0.orig下のフィールド変数ファイルやsystemフォルダ下のforces関連ファイル)はパーツが増えたことに伴う整合性を考えて機械的に変更しただけであるので説明は省く。

また、

  •  system/setFieldDict
  • system/topoSetDict.cHullProp
  • system/fvschemes
については、オーバーセット領域が一つ増えたことに伴う必須の変更(追加)項目があるが、これらも1_hullInPool の解説を思い起こせば容易に理解できるであろう。

 

dynamicMeshDict への以下の追加については補足しておく。

constant/dynamicMeshDict(76〜)
hull-2
{
    type            rigidBody;
    parent          hull-1;
    mass           10.5;

    centreOfMass     (-0.6 0.0 0);
    inertia         (0.0284 0 0 0.884 0 0.895);
    // Transformation tensor and centre of rotation
    transform       (1 0 0 0 1 0 0 0 1)  (-0.5 0.0 0);
    joint
    {
        type            composite;

        joints
        (
            {
                type Ry;
            }
        );
    }
    patches
    (
        hullWall2
    );
    innerDistance   100;
    outerDistance   200;
}

追加したhull-2の形状はhull-1と全く同一であるので、(5行目)mass (8行目)inertia は同一の値としている。

(10行目)transformの座標(-0.5 0 0)は (5行目)parent hull-1 の回転中心からの相対座標値であるので、hull-1 の後端面中心に接続することを想定している。(7行目) centerOfMass はhull-1で見たのと同じで、そこ(transform座標)からの相対座標値ということである。hull-1とhull-2の間隔が0.1なので、ここでも一様物体を想定していることになる。

(15行目)jointsの(18行目)type Ry;はピッチ運動のみ考慮するということである。

Allrun.pre

Allru.preについては、7行目でオーバーセット領域を追加するコマンドが増えることは当然として、6行目においてマージ対象のメッシュが変更されている点にも着目されたい。

計算原理的には、これまで使ってきたオーバーセットメッシュ(Mesh/hull-1)をそのまま使って計算できた可能性もあるが、残念ながら今回の連結構成ではメッシュを作り直す必要があった。というのも連結浮体間の隙間が小さく、オーバーセット計算の分解能が確保できなかったからである。逆に言えば、連結浮体間の間隔がもっと大きければ、これまでのメッシュでも計算可能であったという点は付け加えておく。

cellType

前項で、連結体での計算ではhull-1 のオーバーセットメッシュをそのまま使えなかったと記したが、オーバーセット計算では、計算が開始されると、最初にcellTypeというフィールド変数が自動的に作成されるので、これを調べると計算ができそうなのか、そうでないのかの判断が、ある程度可能になるので記しておく。

下の図は、連結体が無い場合の hull-1 のオーバーセットメッシュを対象に(zoneIdでThresholdフィルターをかける)、カット半断面を作成、cellType を表示したもので、参考の為、中央にhullWall1も表示させてある。

cellTypeは、「0」が計算されるセル、「1」がベースメッシュとの間で補間されるセル、「2」が計算しないセルであり、この場合は「2」は存在しない。

一方、このメッシュを使って連結されるオーバーセットメッシュが存在した場合に、上と同じ表示をさせると、以下のようになってしまう。

参考の為に hullWall2 も表示しておいたが、cellTypeのスケールが変更されている点にも注意されたい。hullWall1の左端面にcellType 「4」のセルが存在するのと、その先が全面「2」となっている点である。cellType 「4」でどのような処理になるか不明であるが、左の全面が「2」であるということは、hullWall1 周りの計算で、仮に計算できたとしても何を計算しているのか、わからないということになる。

そこで、Mesh/hull-1F では、以下のように変更した次第である。

具体的には、左端面近傍でメッシュを細分化しただけである。オーバーセットの全体領域も大きくなっているが、これは試行錯誤の過程が残っただけで、上の例と同じでも良かったと思う。

Mesh/hull-2F について

前項で、Mwsh/hull-1# について説明したが、連結される浮体についてはどうであろう?

実はこちらも、かなりの試行錯誤を経て、公開ケース(下図イメージ)に至っている。

連結隙間のメッシュを細分化し、加えてhull-1F のオーバーセットメッシュを包含する形でhull-2F のオーバーセットメッシュとした点である。これも細分化領域の大きさやオーバーセットの領域全体のサイズをもう少し小さくしても良いかと思っているが、そのあたりの適合にまでは時間をかけられなかった。

当初は、下図に示すように2つのオーバーセット領域を交差させるような形で検討を進めていた。

試行錯誤の過程では、浮体のサイズがキリの良い数字でなかったことも理由なのか、かような構成で計算できたこともあった。オーバーセット領域のサイズを変えているところがミソであった。サイズを同一にすると、気液界面のオーバーセットパッチが接するあたりでたちどころに発散していたが、オーバーセット領域のサイズを変えると計算できるようになった。

しかし最終的に本記事にとりまとめる際に、浮体のサイズを変更して検証しようとしたら、どうしてもうまくいかなかった。ある程度計算は進行するのであるが、気液界面のオーバーセットパッチが交差するあたりが発散起点となった。サイズが異なっていた時には、当該部分で擾乱は生じていたが発散にまで至らなかったということである。

流体力、連結力、 について

流体力計算は、以下のように計算対象を追加した。

system/forces(9〜)
forces1
{
    type          forces;

    libs          ("libforces.so");

    writeControl  timeStep;
    timeInterval  1;

    log           yes;

    patches       (hullWall1);
    rho           rhoInf;     // Indicates incompressible
    log           true;
    rhoInf        1000;          // Redundant for incompressible

    CofR          (0.5 0 0);    // Rotation around centre line of propeller
    pitchAxis     (0 1 0);
}

forces2
{
    type          forces;

    libs          ("libforces.so");

    writeControl  timeStep;
    timeInterval  1;

    log           yes;

    patches       (hullWall2);
    rho           rhoInf;     // Indicates incompressible
    log           true;
    rhoInf        1000;          // Redundant for incompressible

    CofR          (-1.1 0 0);    // Rotation around centre line of propeller
    pitchAxis     (0 1 0);
}

forces3
{
    type          forces;

    libs          ("libforces.so");

    writeControl  timeStep;
    timeInterval  1;

    log           yes;

    patches       (hullWall1 hullWall2);
    rho           rhoInf;     // Indicates incompressible
    log           true;
    rhoInf        1000;          // Redundant for incompressible

    CofR          (-0.5 0 0);    // Rotation around centre line of propeller
    pitchAxis     (0 1 0);
}

(1行目)forces1 のブロックは、単一の浮体で使っていたもので、これをべースに(21行目)forces2、(41行目)forces3のブロックへコピペして要所変更して使っている。forces2は、連結された浮体を想定、forces3は両者を併せて(57行目)CofR(回転中心)を連結箇所の座標値に設定したものである。

 

まとめ

DEXCS最新版(OF-v2406)のoversetInterDyMFoamを使用、単純形状サンプルにて以下確認し、主要パラメタの設定要点をとりまとめた。

  1. 任意の位置で連結された2つの浮体が流体力による3自由度(surge,pitch,heave)の剛体運動を伴う二相流計算ができた。
  2. 仕様1(正弦波)を含むかどうかは定かでないがOpenFOAMに標準のwaveタイプの境界条件は適用可能であった。
  3. 仕様3(等速運動)は、dynamicMeshMotionLib の適用を考えたが剛体運動が所望の挙動をしてくれなかった。バグがあるのか、そもそも使用法に問題があるのかどうかは不明。
  4. 上記代案として、浮体に推進力を付加させることで前進挙動を実現可能であった。推進力のタイムテーブルを適合すれば、準等速運動は可能。
  5. 仕様5(連結部に作用する力)はforcesポスト処理にて、連結支持部を作用点として算出した値で代用できると考えられるが、十分に検証できていない。
最後に、全仕様を組み込んだケース(Case/6_hullsMoveInWavePool)の結果アニメーションと、連結支持部における流体力プロットを以下に示しておく。

補足事項

補足1

テキストファイルの説明上の行番号は「1」から始まっているが、実際のテキストファイルの行番号は(〜#)の#(数字)を加えた行番号になる。

補足2

drivenLinearMotionで検索すると、ソースコードの在所が、src/meshTools/solidMotionFunctions/drivenLinearMotion であることがわかる。このあたりをファイルマネージャーで調べると、

となっているので、linearMotion が使えそうだとわかる。他にもある。

補足3

標準チュートリアルケース(rigidBodyHull)で使用されていた座標値の値は、あたかも重心も回転中心も絶対座標で定義してあるかの同一数値となっていた。

centreOfMass     (0.178 0 0.3323);
inertia         (564 0 0 8535 0 8535);
// Transformation tensor and centre of rotation
transform       (1 0 0 0 1 0 0 0 1)  (0.178 0 0.3323);

これを変えて挙動を調べれば検証できたかもしれないが、計算負荷も大きいので、1_hullInPool で以下変更して検証した。

centreOfMass     (0.5 0.0 0.0);
inertia         (0.0284 0 0 0.884 0 0.895);
// Transformation tensor and centre of rotation
transform       (1 0 0 0 1 0 0 0 1)  (0.5 0 0);

こうすると、重心位置は浮体の右端になっているはずである。結果のアニメーションを追加しておく。

補足4

補足5

Share

コメントする

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください