DynamicMeshを使った泳動作中の計算が発散することについて

[wd.ta...@gmail.com]

お世話になっております．

OpenFoam初心者です．

早速ですが，現在自分の研究で

人の泳動作中の流体シミュレーションを実施しようとしております．

使用ソフトはopenfoam6，流体は非圧縮性の非定常流体，乱流モデルはk-εモデルで計算しています．

シミュレーションの前に

・水泳選手の形状情報を，3次元人体計測装置を用いて取得

・水泳選手が水中で泳いでいる時の運動学的データをモーションキャプチャを用いて取得
これら取得したデータを用いて，シミュレーションを実施しようとしています．

OpenFoam演算の流れをまとめると

① blockMeshで，11.0×1.5×1.5(x, y z )のメッシュを作成

② 水泳選手の形状メッシュを，snappyHexMeshを使って，blockMeshで作成した計算格子の中に作成

③ topoSetDictとCreatPatchを用いて，水泳選手の形状メッシュを，下肢・上肢というようにセグメント毎に定義し直す（添付；図1）

④ pimpleFoamを用いて，水中で泳動作をしていない状態のシミュレーションを実施(添付；図2)

⑤ ④のlatestTimeでのデータを，初期値として泳動作中のシミュレーションを実施

　・移動メッシュには，dynamicMotionSolverFvMeshのdisplacementLaplacian，diffusivityには，inverseDistanceを使用

　・0ディレクトリに，PointDisplacementファイルを作成し，solidBodyMotionDisplacementのtabulated6DoFMotionを用いて，実際に測定した泳動作の運動学データを入力

　なお，運動学データは，10ms毎にしております．

以上の流れで行っています．

④までは計算，出力した抗力係数から計算し直した抗力等は，先行研究の結果と同程度の結果が得られています．

しかし，⑤を実施するとおそらく圧力だと思うのですが，計算開始後，0.03秒後には発散してしまい計算が止まってしまいます（添付；図3）．

fvsolutionで緩和係数を小さくするや，solverをGAMGからPCG/PBiCGStabに変えるなど，計算を安定させる方法を試してみましたがうまく計算できません．

発散した時の結果をparaFoamで見てみたところ，水泳選手の身体に作用している圧力がおかしいかと思いますが，解決策が見つかりません．

こちらの問題を解決する方法がありましたらご教授頂ければ幸いです．

なお，先行研究でOpenFoamを用いて人の泳動作中のシミュレーションは行われているので，計算はできるかと思います．

長文で失礼しましたが，もし足りない情報等もありましたら

ご教授ください．

よろしくお願いいたします．

下記，発散する前のエラーです．

-----------------------------------------------------

PIMPLE: Iteration 1

DICPBiCGStab:  Solving for cellDisplacementx, Initial residual = 1.23062e-05, Final residual = 9.41443e-06, No Iterations 1

DICPBiCGStab:  Solving for cellDisplacementy, Initial residual = 0.0287695, Final residual = 0.000273505, No Iterations 5

DICPBiCGStab:  Solving for cellDisplacementz, Initial residual = 9.33004e-06, Final residual = 9.33004e-06, No Iterations 0

DICPCG:  Solving for pcorr, Initial residual = 1, Final residual = 9.5428e-07, No Iterations 314

DICPCG:  Solving for pcorr, Initial residual = 0.10179, Final residual = 9.47478e-07, No Iterations 273

time step continuity errors : sum local = 4.06619e-15, global = 1.63796e-17, cumulative = -7.61575e-08

DILUPBiCGStab:  Solving for Ux, Initial residual = 0.0137191, Final residual = 4.79097e-05, No Iterations 2

DILUPBiCGStab:  Solving for Uy, Initial residual = 0.000233217, Final residual = 1.89075e-06, No Iterations 1

DILUPBiCGStab:  Solving for Uz, Initial residual = 6.94858e-05, Final residual = 7.60524e-06, No Iterations 1

DICPCG:  Solving for p, Initial residual = 0.999997, Final residual = 0.0919308, No Iterations 6

DICPCG:  Solving for p, Initial residual = 0.0345613, Final residual = 9.39654e-06, No Iterations 258

time step continuity errors : sum local = 1.31394e-13, global = -8.74637e-16, cumulative = -7.61575e-08

PIMPLE: Iteration 2

DILUPBiCGStab:  Solving for Ux, Initial residual = 0.905905, Final residual = 0.00498697, No Iterations 1000

DILUPBiCGStab:  Solving for Uy, Initial residual = 0.467712, Final residual = 5.84424e+30, No Iterations 1000

DILUPBiCGStab:  Solving for Uz, Initial residual = 0.329059, Final residual = 1.01585e+18, No Iterations 1000

DICPCG:  Solving for p, Initial residual = 1, Final residual = 0.0970522, No Iterations 20

DICPCG:  Solving for p, Initial residual = 1.68339e-08, Final residual = 1.68339e-08, No Iterations 0

time step continuity errors : sum local = 1.68036e+25, global = -1.48677e+10, cumulative = -1.48677e+10

DILUPBiCGStab:  Solving for epsilon, Initial residual = 1, Final residual = 9.93296e-08, No Iterations 21

bounding epsilon, min: -1.32257e+56 max: 1.13488e+87 average: 1.89213e+84

DILUPBiCGStab:  Solving for k, Initial residual = 1, Final residual = 8.76068e-08, No Iterations 6

bounding k, min: -1.85603e-20 max: 1.6557e+88 average: 4.8212e+82

ExecutionTime = 4625.21 s  ClockTime = 7088 s

forces forces write:

    sum of forces:

        pressure : (-1.21769e+99 8.29774e+98 2.17187e+99)

        viscous  : (-4.90015e+100 5.42906e+100 7.04984e+100)

        porous   : (0 0 0)

    sum of moments:

        pressure : (-1.43268e+98 1.60044e+99 -6.96082e+98)

        viscous  : (4.69946e+99 5.00865e+100 -3.81125e+100)

        porous   : (0 0 0)

forceCoeffs forceCoeffs write:

    Cm    = 2.87446e+97

    Cd    = 6.0046e+97

    Cl    = 8.68902e+97

    Cl(f) = 7.21897e+97

    Cl(r) = 1.47005e+97

Courant Number mean: 1.37561e+26 max: 5.06915e+32

deltaT = 1.49848e-64

Time = 0.0305065

PIMPLE: Iteration 1

DICPBiCGStab:  Solving for cellDisplacementx, Initial residual = 3.28872e-05, Final residual = 9.38091e-06, No Iterations 1

DICPBiCGStab:  Solving for cellDisplacementy, Initial residual = 0.00339951, Final residual = 3.13542e-05, No Iterations 13

DICPBiCGStab:  Solving for cellDisplacementz, Initial residual = 9.10639e-06, Final residual = 9.10639e-06, No Iterations 0

DICPCG:  Solving for pcorr, Initial residual = 1, Final residual = 9.23793e-07, No Iterations 326

DICPCG:  Solving for pcorr, Initial residual = 0.0670329, Final residual = 9.95298e-07, No Iterations 267

time step continuity errors : sum local = 4.82617e-14, global = 2.29628e-16, cumulative = -1.48677e+10

[1] #0  Foam::error::printStack(Foam::Ostream&) at ??:?

[1] #1  Foam::sigFpe::sigHandler(int) at ??:?

[1] #2  ? in "/lib/x86_64-linux-gnu/libc.so.6"

[1] #3  Foam::DILUPreconditioner::calcReciprocalD(Foam::Field<double>&, Foam::lduMatrix const&) at ??:?

[1] #4  Foam::DILUPreconditioner::DILUPreconditioner(Foam::lduMatrix::solver const&, Foam::dictionary const&) at ??:?

[1] #5  Foam::lduMatrix::preconditioner::addasymMatrixConstructorToTable<Foam::DILUPreconditioner>::New(Foam::lduMatrix::solver const&, Foam::dictionary const&) at ??:?

[1] #6  Foam::lduMatrix::preconditioner::New(Foam::lduMatrix::solver const&, Foam::dictionary const&) at ??:?

[1] #7  Foam::PBiCGStab::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?

[1] #8  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"

[1] #9  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"

[1] #10  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"

[1] #11  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"

[1] #12  ? in "/opt/openfoam6/platforms/linux64GccDPInt32Opt/bin/pimpleFoam"

[df096d4daeb9:00668] *** Process received signal ***

[df096d4daeb9:00668] Signal: Floating point exception (8)

[df096d4daeb9:00668] Signal code:  (-6)

[df096d4daeb9:00668] Failing at address: 0x1f50000029c

[df096d4daeb9:00668] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x3ef20)[0x7f721d963f20]

[df096d4daeb9:00668] [ 1] /lib/x86_64-linux-gnu/libc.so.6(gsignal+0xc7)[0x7f721d963e97]

[df096d4daeb9:00668] [ 2] /lib/x86_64-linux-gnu/libc.so.6(+0x3ef20)[0x7f721d963f20]

[df096d4daeb9:00668] [ 3] /opt/openfoam6/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam18DILUPreconditioner15calcReciprocalDERNS_5FieldIdEERKNS_9lduMatrixE+0x90)[0x7f721ed42ec0]

[df096d4daeb9:00668] [ 4] /opt/openfoam6/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam18DILUPreconditionerC1ERKNS_9lduMatrix6solverERKNS_10dictionaryE+0x49)[0x7f721ed43009]

[df096d4daeb9:00668] [ 5] /opt/openfoam6/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam9lduMatrix14preconditioner31addasymMatrixConstructorToTableINS_18DILUPreconditionerEE3NewERKNS0_6solverERKNS_10dictionaryE+0x2e)[0x7f721ed4313e]

[df096d4daeb9:00668] [ 6] /opt/openfoam6/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZN4Foam9lduMatrix14preconditioner3NewERKNS0_6solverERKNS_10dictionaryE+0x23d)[0x7f721ed2e07d]

[df096d4daeb9:00668] [ 7] /opt/openfoam6/platforms/linux64GccDPInt32Opt/lib/libOpenFOAM.so(_ZNK4Foam9PBiCGStab5solveERNS_5FieldIdEERKS2_h+0x5d5)[0x7f721ed36e25]

[df096d4daeb9:00668] [ 8] pimpleFoam(+0x65239)[0x55c52e6ec239]

[df096d4daeb9:00668] [ 9] pimpleFoam(+0x7691c)[0x55c52e6fd91c]

[df096d4daeb9:00668] [10] pimpleFoam(+0x2c104)[0x55c52e6b3104]

[df096d4daeb9:00668] [11] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0xe7)[0x7f721d946b97]

[df096d4daeb9:00668] [12] pimpleFoam(+0x2dfea)[0x55c52e6b4fea]

[df096d4daeb9:00668] *** End of error message ***

--------------------------------------------------------------------------

mpirun noticed that process rank 1 with PID 0 on node df096d4daeb9 exited on signal 8 (Floating point exception).

--------------------------------------------------------------------------

----------------------------------------------------------------------------

　

[Inaba]

こんにちは。

稲葉と申します。

とても面白そうな内容ですね。

静止状態できちんと回っているようなので、移動でメッシュが歪んでいるのではないかと思いました。

⑤のケースをまるまるコピーして`moveDynamicMesh`を実行するとpimpleFoamなしで泳動のみをtimeStepごとに行ってくれます。

paraViewやlogファイルでちゃんと指定通りに移動しているか、その際にメッシュが反転したり重複したりしていないかを確認すると良いかと思います。

以上です。

ご参考まで。

稲葉

2020年8月20日木曜日 9:18:53 UTC+9 wd.ta...@gmail.com:

[wd.ta...@gmail.com]

稲葉様

田中と申します．

ご返事ありがとうございます．

"moveDynamicMesh"を実行したところ，入力した運動学データで選手の形状が動いているのは観察できました．

しかし，logとparaViewを確認すると，稲葉様がおっしゃる通り，

・メッシュの非直行性

・アスペクト比

・skewness

・zero or negative cell volume detected

とエラーが出ており，値もかなりおかしいことがわかりました．

paraviewでの確認で，特に胸あたりでメッシュがのびたり

メッシュの重なりも生じていそうな部分もございます(図1，2）．

そもそも人の身体を剛体ではなく，皮膚の伸び縮み等があるものなので

こうゆう問題が生じるとは思うのですが・・・

僕の分野では人の身体を剛体としてみなして解析するので

皮膚の伸び縮み等を考慮しない場合が多いです．

これらを修正しなければならないと思ったのですが，どのように修正すればよろしいでしょうか．

入力する運動学データをメッシュが歪まないように加工するや，水泳選手の形状をいじるということもできなくはないですが

自分が解明したいこととずれてしまうので，極力避けたいと考えています．

一部分ですが，動かした時の選手の形状のメッシュと選手の周りのメッシュのスクリーンショット（図1-3），

"moveDynamicMesh"のlogを添付いたします．

よろしくお願いいたします．

2020年8月21日金曜日 5:42:28 UTC+9 Inaba:

[wd.ta...@gmail.com]

もしくは，メッシュをもう少し細かくした方がよろしいでしょうか．

今，メッシュの数は約67万セルです．

先行研究では，約100万セルと書かれていたので，この研究と比べると数が少ないですが，

細かくしすぎると，それはそれで静止状態でも計算が発散しやすくはなりました．

2020年8月21日金曜日 8:07:57 UTC+9 wd.ta...@gmail.com:

[Inaba]

田中様

こんにちは。

稲葉です。

返事が遅くなってしまったためすでに解決済みでしたら申し訳ないです。

このような伸縮する物体を計算したことがないので感想程度しか申し上げることができません。

① もしすべての部分が伸びているシーンがあるのであれば、その時の形状でメッシュを作成し、それをdynamicMeshで巻き戻し、それを初期条件とするのはどうでしょうか。

例えばひざと胸の部分でつの面が縦長になってしまっていますが、このシーンでメッシュを切ると初期状態ではそこは縮むだけですので計算できる可能性は高いかもしれません。

ただ今度は背中であったり首後であったりの状況が変わるので、最もバランスのいいタイミングがもしかしたらあるかもしれません。

② 先行研究と比較する場合、もし先行研究がsnappyHexMeshで切ったのでないのであれば、メッシュ数を参考にするのはよくないかもしれません。

snappyであっても、人体の周りのみ細分化されていたりなどもありますので、先行研究のメッシュ数に引っ張られずにいろいろと試してみる方がよいと感じました。

③ メッシュを断面で確認する際には

まずVTK polyhedraにチェックを入れた後に

Clipの場合：Crinkle clip

Sliceの場合：Crinkle Slice

で確認する方が、メッシュのセル一つ一つを確認できるのでお勧めです。

以上です。

これだというアドバイスができず心苦しい限りです。

頑張ってください。

稲葉

2020年8月21日金曜日 8:22:23 UTC+9 wd.ta...@gmail.com:

[wd.ta...@gmail.com]

稲葉様

ご連絡ありがとうございます．

いろいろ試してはいますが，まだ解決策が見つかっていない状態です．

たくさんの貴重なご意見ありがとうございます．

試してみます．

②ですが，snappyHexMeshを使っているようです．当方でも，メッシュ数を変えてみましたが，極端な設定にしない限りは静止状態でのメッシュの品質はほとんど同じでした．

また，blockmeshの分割数や形状の細分化レベルを上げすぎると，メッシュが作成されないこともありましたので，そこまで大きな修正はできないかもしれないと感じています．

現在，当方でも試しに，全身ではなく，例えば足だけ上半身だけといったセグメントにわけ，身体形状のメッシュの伸び等が含まれない形(セグメントを剛体とみなして計算)でシミュレーションをしてみてますが，

動きが最も大きい足だけのモデルだと，すぐに計算が発散してしまいます．

同じように，メッシュの歪み等を確認すると，足の形状メッシュではなく，足の周りの計算メッシュが歪むことで，計算が発散しているのでは？と感じます．

moveDynamicMeshをしメッシュの品質も確認すると，メッシュの非直行性など，全身で動かした時と似たような数字になっていました．

今，DynamicMeshDictの設定を，

motionSolverをdisplacementLaplacianでdiffusivityをinverseDistance(足のパッチ名)に設定しているのですが，

今回のような剛体移動の場合，うまく計算できそう・計算できたという設定の方法等ありませんでしょうか．

個人的には，diffusivityをdirectionalといった他のものに変える方法が一つかと思うのですが．．

ご教授頂ければ幸いです．

田中

2020年9月2日水曜日 21:42:26 UTC+9 Inaba:

[wd.ta...@gmail.com]

ちなみにですが，

計算時間は，1.3秒で

全身だと 0.03秒

足だけのセグメントだと 0.05秒

現在計算中の上半身だけのセグメントだと，0.31秒(まだ発散せず計算できています)

となっております．

2020年9月4日金曜日 11:06:17 UTC+9 wd.ta...@gmail.com:

[wd.ta...@gmail.com]

みなさま

ご無沙汰しております．

先日，こちらのフォームで質問させていただきました件ですが，

2016年ごろに，こちらの投稿で再メッシュ切りの議論があり，参考にさせていただきましたところ，うまく計算を行うことができました．

https://groups.google.com/g/openfoam/c/aJL0kP-rvyQ/m/uG24pQxUDgAJ

ちなみに，先行研究でもこちらの方法が使われているようでした．

まだ，データの精度は上げられそうですが，今行っている計算結果と先行研究の結果を定性的に比較したところ，正しく計算されていそうです．

行き詰まっていたところ，こちらのGoogle Groupに助けていただきました．

ありがとうございました．感謝しております．

また，二つほど質問がございます．

①ControlDictのFunction部分で，forces，forcesCoeffs，Q valueを出力しているようにしております．

これらのデータは，シミュレーション(pimpleFoam)をしながら計算されておりますが，

シミュレーションが全て終了した後，pimpleFoamを実行せずに，これらの値のみ計算し直すことは可能でしょうか．

②forcesとforcesCoeffsですが，力情報や係数情報を出したいpatchを選択する時に，

patchごと(例えば，私の場合だと人の足部，下腿，大腿・・・　それぞれのセグメントに作用している圧力やその力情報から計算される係数)を出すことは可能でしょうか．

今，patchの部分を

patches         (upperbody upchest chest upperwaist lowerwaist thigh shank foot);

のように記述しており．この設定であると人の全身に作用する圧力として算出されております．

以上，もし方法をご存知の方がいらっしゃいましたら，ご教授いただければ幸いです．

よろしくお願いいたします．

2020年9月4日金曜日 11:17:23 UTC+9 wd.ta...@gmail.com:

[penguinitis]

① 計算し終わったケースにおいて、controlDict で functions の設定をして

$ pimpleFoam -postProcess

のようにすれば、結果処理だけ実行できます。この場合は、-latestTime などの時間も可能です。

② controlDict の functions の設定は

functions

{

    forces

    {

        type forces;

        patches         (upperbody upchest chest upperwaist lowerwaist thigh shank foot);

        ...

    }

}

のような感じになっていると想像しますが、最初の forces は任意の名前なので、

これを変えて複数の forces の設定をすることが可能だったと思います。

functions

{

    forcesUpperbody

    {

        type forces;

        patches (upperbody);

        ...

    }

    forcesUpchest

    {

        type forces;

        patches (upchest);

        ...

    }

   ...

}

こんな感じで。

2020年10月6日火曜日 22:19:13 UTC+9 wd.ta...@gmail.com:

[wd.ta...@gmail.com]

ありがとうございます．

ご教示いただきました方法で計算ができ，自分が欲しかったデータを出力することができました．

ありがとうございました．

2020年10月7日水曜日 10:31:09 UTC+9 penguinitis:

[md]

お世話になっております。openfoamの初学者です。

私も現在 solidBodyMotionDisplacementのtabulated6DoFMotionを用いた解析を行おうと考えております。

運動をdatファイルでしていると思うのですが、水中での泳動作のような複雑な運動をどのように定義したのか教えていただきたいです。

よろしくお願いいたします。

2020年8月20日木曜日 9:18:53 UTC+9 wd.ta...@gmail.com:

お世話になっております．

[wd.ta...@gmail.com]

こんにちは。

①ヒトの体をtopoSetのコマンドを使って、大腿部・下腿部・足部・・・と身体のセグメントごとにパッチの定義を変更。

②それぞれの身体セグメントの運動データを、matlabを使ってdatファイルを作成。

　作成したdat fileは、各セグメントの位置と角度の時系列データ。

　中身は下記のような形になります。

　

　XXX（タイムステップ数）

　(

　(0.0 ((x y z) (x y z)))

    (0.1 ((x y z) (x y z)))

    ・

　・

　・

　・

　(1.0 ((x y z) (x y z)))

    )

　横から(時間 ((xyzの位置) (xyz軸周りの回転角度)))とデータが並んでいます。

③0ディレクトリに入れるpointDisplacementは、チュートリアル等を参考にして

　

　動かしたいパッチ名

　{

　　　type　　　　solidBodyMotionDisplacement;

　　　CofG　　(x y z); // 動かしたいパッチの回転中心の座標を定義

　　　timeDataFileName　　"FOAM_CASE/constant/XXXX.dat"; //②で作ったdat file名

　}

　

　これを①で定義したパッチ分記載する。

としました。

少々わかりづらく申し訳ありません。

　

2022年11月7日月曜日 15:32:45 UTC+9 md:

[md]

早速ご返信ありがとうございます。

なるほど。topoSetで分割することができるのですね。

大変参考になりました。ご丁寧にありがとうございます。

また何かわからないことがありましたらご教授いただけると幸いです。

よろしくお願いいたします。

2022年11月7日月曜日 15:50:56 UTC+9 wd.ta...@gmail.com: