2019-06-29

Fluent: journal file を複数に分割して include するには→ file/read-journal を使う

ANSYS Fluent programming

まえがき

これ、できないと思いこんでたけど、よく見たら User's Guide に書いてあった。.set ファイルに頼らずに、なるべく journal file でやろうとすると、ファイルがどんどん長くなっていて、困っていたので、これができると非常に嬉しい。

方法

要は以下のコマンドで他の journal file(s) を read するというだけ：

file/read-journal "<path_to_another_jorunal_file>" ""

最後に2つのダブルクォーテーションを置くことを忘れないように（忘れると、なぜか次のコマンドがファイル名として認識されてエラーになる）。

ファイルパスは、同一ディレクトリで良ければ、ふつうに、

file/read-journal "./sub.jou" ""

って感じでOK.

もし複数の journal files を連続して読み込みたいなら、

file/read-journal "<path_to_1st_jorunal_file>" "<path_to_2nd_jorunal_file>" ""

みたいにいける。まぁ行数節約しなくていいなら複数行に分けて書いたほうが見やすい気もするが。

具体例

こんなのやろうって人は自分でいくらでも思いつくとは思うけど、たとえば、smoothing の方法を切り替えやすくなる。

つまり、以下のように、smoothing_common.jou, smoothing_spring.jou, smoothing_diffusion.jou, smoothing_linelast.jou というファイルを作っておいて、

smoothing_common.jou

;; ****************************************
;; COMMON settings for SPRING, DIFFUSION, or LIN-ELASTIC method (see below for method-specific settings)
;;; Hidden settings (not accessible from GUI)
;;;; Skewness smoothing on all deforming boundaries [yes]
;;;;   ? Enable/disable skewness smoothing for all deforming dynamic boundary zones. If disabled, only the deforming dynamic boundary zones are smoothed which have smoothing explicitly enabled or use local face remeshing.
define/dynamic-mesh/controls/smoothing-parameters/skew-smooth-all-deforming-boundaries?
yes
;;;; Skewness threshold for cell [0.9]
;;;;   ? Set the cell skewness threshold above which cells will be smoothed using the skewness method.
define/dynamic-mesh/controls/smoothing-parameters/skew-smooth-cell-skew-max
0.7
;;;; Skewness threshold for face [0.7]
;;;;   ? Set the face skewness threshold above which deforming boundary faces will be smoothed using the skewness method.
define/dynamic-mesh/controls/smoothing-parameters/skew-smooth-face-skew-max
0.5
;;;; Skewness-based smoothing number of iterations [4]
define/dynamic-mesh/controls/smoothing-parameters/skew-smooth-niter
4

smoothing_spring.jou

;; ****************************************
;; SPRING method
;;    As boundary layer deformation smoothing can only be used with SPRING, this might be the one we should use.
;;; Smoothing method [spring] (spring, diffusion, linelast)
define/dynamic-mesh/controls/smoothing-parameters/smoothing-method
spring
;;; Normal settings (accessible from GUI)
;;;; <SPRING> Max iterations [20]
define/dynamic-mesh/controls/smoothing-parameters/max-iter
100
;;;; <SPRING> Spring constant factor [1]
define/dynamic-mesh/controls/smoothing-parameters/constant-factor
1
;;;; <SPRING> Convergence Tolorance [0.001]
define/dynamic-mesh/controls/smoothing-parameters/convergence-tolerance
0.0001
;;;; <SPRING> Laplace Node Relaxation factor [1]
define/dynamic-mesh/controls/smoothing-parameters/laplace-node-relaxation
1
;;;; <SPRING> Enable/disable spring-based smoothing for all cell shapes [no]
define/dynamic-mesh/controls/smoothing-parameters/spring-on-all-elements?
no
;;;; <SPRING> Enable/disable spring-based smoothing for tri/tet elements in mixed element zones [no]
define/dynamic-mesh/controls/smoothing-parameters/spring-on-simplex-elements?
no
;;; Hidden settings (not accessible from GUI)
;;;; <SPRING> Set the spring boundary node relaxation factor [1]
;;;;   ? The boundary node relaxation is used by spring smoothing. The boundary node relaxation allows you to relax the update of the node positions at deforming boundaries. A value of 0 prevents deforming boundary nodes from moving and a value of 1 indicates no under-relaxation.
define/dynamic-mesh/controls/smoothing-parameters/bnd-node-relaxation
1
;;;; <SPRING> Set the stiffness factor for springs connected to boundary nodes [1]
define/dynamic-mesh/controls/smoothing-parameters/bnd-stiffness-factor
1

smoothing_spring.jou

;; ****************************************
;; DIFFUSION method
define/dynamic-mesh/controls/smoothing-parameters/smoothing-method
diffusion
;;; Normal settings (accessible from GUI)
;;;; <DIFFUSION> Specify whether the diffusion coefficient is based on the boundary distance or the cell volume [boundary-distance] (boundary-distance, cell-volume)
define/dynamic-mesh/controls/smoothing-parameters/diffusion-coeff-function
boundary-distance
;;;; <DIFFUSION> Set the diffusion coefficient parameter used for diffusion-based smoothing [0] (0.0 -- 3.0)
define/dynamic-mesh/controls/smoothing-parameters/diffusion-coeff-parameter
1.5
;;; Hidden settings (not accessible from GUI)
;;;; <DIFFUSION> Set the numerical method used for diffusion-based smoothing [no] (no=FEM; yes=FVM)
;;;;   ? Answering yes at the prompt changes the diffusion-based smoothing method to the cell based finite volume approach that was the default in releases prior to Fluent 15.0. Answering no at the prompt changes the diffusion-based smoothing method to the default node-based finite element method.
define/dynamic-mesh/controls/smoothing-parameters/diffusion-fvm?
no
;;;; <DIFFUSION> Set the method used to evaluate the boundary distance for the diffusion coefficient calculation [no] (no=standard wall distance; yes=generalized boundary distance)
define/dynamic-mesh/controls/smoothing-parameters/boundary-distance-method
no
;;;; <DIFFUSION/LIN-ELAST> maximum number of iterations for mesh smoothing (FEM) [30]
define/dynamic-mesh/controls/smoothing-parameters/max-iter
30
;;;; <DIFFUSION/LIN-ELAST> Set the relative convergence tolerance for mesh smoothing (FEM) [1e-10]
define/dynamic-mesh/controls/smoothing-parameters/relative-convergence-tolerance
1e-15
;;;; <DIFFUSION/LIN-ELAST> Enable smoothing from reference position [no]
;;;;   ? Enables/disables smoothing from a reference position. Such smoothing may produce greater mesh quality consistency for cases with periodic or quasi-periodic motion, and is only available when the smoothing method is based on diffusion or the linearly elastic solid model.
define/dynamic-mesh/controls/smoothing-parameters/smooth-from-reference-position?
no
;;;; <DIFFUSION/LIN-ELAST> Set the verbosity for mesh smoothing (FEM) [0] (0 -- 1)
;;;;   ? Setting this to 1 will cause smoothing residuals to be printed to the text console. The default value of 0 suppresses this output.
define/dynamic-mesh/controls/smoothing-parameters/verbosity
1

smoothing_linelast.jou

;; ****************************************
;; LIN-ELAST method (linearly elastic solid based)
define/dynamic-mesh/controls/smoothing-parameters/smoothing-method
linelast
;;; Normal setting (accessible from GUI)
;;;; Poission ratio [0.45] (-1.0 -- 0.5)
define/dynamic-mesh/controls/smoothing-parameters/poisson-ratio
0.45
;;; Hidden settings (not accessible from GUI)
;;;; <DIFFUSION/LIN-ELAST> maximum number of iterations for mesh smoothing (FEM) [30]
define/dynamic-mesh/controls/smoothing-parameters/max-iter
30
;;;; <DIFFUSION/LIN-ELAST> Set the relative convergence tolerance for mesh smoothing (FEM) [1e-10]
define/dynamic-mesh/controls/smoothing-parameters/relative-convergence-tolerance
1e-4
;;;; <DIFFUSION/LIN-ELAST> Enable smoothing from reference position [no]
;;;;   ? Enables/disables smoothing from a reference position. Such smoothing may produce greater mesh quality consistency for cases with periodic or quasi-periodic motion, and is only available when the smoothing method is based on diffusion or the linearly elastic solid model.
define/dynamic-mesh/controls/smoothing-parameters/smooth-from-reference-position?
no
;;;; <DIFFUSION/LIN-ELAST> Set the verbosity for mesh smoothing (FEM) [0] (0 -- 1)
;;;;   ? Setting this to 1 will cause smoothing residuals to be printed to the text console. The default value of 0 suppresses this output.
define/dynamic-mesh/controls/smoothing-parameters/verbosity
1

main.jou

メインとなる journal では以下のように、必要ないものをコメントアウトすればよい：

<略>
;; Enable smoothing
define/dynamic-mesh/controls/smoothing?
yes
;
file/read-journal "./smoothing_common.jou" ""
file/read-journal "./smoothing_spring.jou" ""
; file/read-journal "./smoothing_diffusion.jou" ""
; file/read-journal "./smoothing_linelast.jou" ""
;
<略>

他にも、たとえば display/views/camera の設定をいろいろと保存しておいて、見たい角度などに応じて呼び出しを切り替える、といったこともできるだろう。

2019-06-26

UDF: cell thread に含まれる face thread にアクセスする

Fluent programming

UDF において、ある cell thread tc（たとえば翼を含む空間）に含まれる face thread tf（たとえば翼表面）を求めるのに、次のようなことができた：

begin_c_loop(c,tc);
c_face_loop(c,tc,i)
{
  if(C_FACE_THREAD(c,tc,i)==tf)
  {
    f=C_FACE(c,tc,i);
    f_node_loop(f,tf,n)
    {
       np=F_NODE(f,tf,n);
       if( NODE_POS_NEED_UPDATE(np) )
      {
         NODE_POS_UPDATED(np);
         do something;
      }
    }
  }
}
end_c_loop(c,tc);

if(C_FACE_THREAD(c,tc,i)==tf) のところ、UDFマニュアルのサンプルだと if( THREAD_ID(C_FACE_THREAD(c,tc,i)) == THREAD_ID(tf) ) のように書いてるんだけど、Thread ポインタのままで比較しても特に問題はなさそうだった。

…あれ？でもなんかこれあまり有意義ではないかも・・・

2019-06-25

Fluent UDF: マクロ末尾の _M1 は「ひとつ前の時刻」の意味だった

Fluent

UDF マニュアルより

The _M1 suffix can be applied to some of the cell variable macros in Table 3.8: Macros for Cell Flow Variables Defined in mem.h or sg_mem.h (p. 262) to allow access to the value of the variable at the previous time step (that is, $t - \Delta t$ ).

同様に、_M2 は「2つ前の時刻 ( $t - 2 \Delta t$ )」の意味だった。

2019-06-16

ペンギン遊泳・フリッパー関連の先行研究

swimming journal club biomechanics

論文リスト

Kinematics of swimming of penguins at the Detroit Zoo (1979)
The energetics of 'flying' and 'paddling' in water: locomotion in penguins and ducks (1985)
Penguin Swimming. I. Hydrodynamics (1988)
Penguin Swimming. II. Energetics and Behavior (1988)
Functional anatomy of the penguin flipper (1992)
The potential costs of flipper-bands to penguins (2002)
Lovvorn, J., Liggins, G. A., Borstad, M. H., Calisal, S. M. & Mikkelsen, J. 2001 Hydrodynamic drag of diving birds: effects of body size, body shape and feathers at steady speeds. J. Exp. Biol. 204, 1547–1557.
Shufeldt, 1901
Meister, 1962
(Dablow, 1925; Neu, 1931; Kooyman et al., 1971).

Kinematics of swimming of penguins at the Detroit Zoo

計測はいくつかの手法で行っている。
1. 時間計測
2. 2台のSony 3210カメラ（情報が出てこないがおそらくアナログフィルム）により60 fpsで水路の中を泳ぐペンギンを撮影。2台使っているがこれはステレオ撮影のためではなく、単に広い領域を撮影するためのようだ。calibration とか three-dimension とかいう言葉は全く出てこないし、結果もない。また、水路の外から撮影しているにもかかわらず、ガラスと水による屈折の影響を全く考慮していない。
3. 1台のカメラによる撮影から、体長で割った速度と、羽ばたき周波数を計測。
4. 24 fps で？撮影。ちょっとここは自分のアナログ機器の知識不足で、何を言ってるのかよくわからない…。Fig. 2 はこの結果のようだが、羽ばたき1周期がたったの8フレーム・6フレームしかないので、24 fps というのは妥当に思える。
種は色々いるが、撮影できたのは主に African penguin（論文中での呼称は blackfoot penguin).
連続羽ばたきのみを解析の対象とし、間欠的な羽ばたき⇔滑空は解析から除いている。
翼の迎え角は当然計測できていない。というか羽ばたき振幅も画像の解像度不足で計測できていない。これらは率直に書いてあって、印象が良い。
やたらと stride length の利用を推しているが、この方法で整理することでなにが嬉しいのかわからない。
アデリーで、gliding 直前の downstroke duration (0.35 s) は連続羽ばたき中のそれ (0.54 s) よりも短い、つまりストロークが速い。これはちょっと興味深い。つまり、glide する前には1回だけ、ただし強めに羽ばたいておく、ということをやっている可能性がある。
エンペラーが姿勢を変えずに 2.18 m 滑空したケースがあり、ここから抗力係数を推算している。体重すら計測せずに文献値なのか…水族館あまり協力的じゃなかったのかな…。"This force is assumed to be exclusively friction drag (D)" は意味がわからない。別に pressure drag 成分を除く根拠はないだろうに。普通に total drag でいい。Discussion の最後の記述とか見ると流体力学わかってないという感じではないので、なにか勘違いしてたのかな…。
- また、結果としての抗力係数は 0.0021 から 0.003 とあるがさすがに小さすぎる。これは代表面積のとり方が表面積だからというだけでは説明がつかない。Lovvorn et al., 2001 によると little penguin （LIPEと略記）の模型・冷凍死骸とも、同等の Re（約100万）で最小でも 0.02 程度（代表面積を表面積にした場合。前面投影面積にすると1桁増える）。つまり、計測誤差が大きいというのもあり得るし、and/or 実際にはわずかに羽ばたくとか胴体を揺するとかして、このエンペラーはわずかでも推進力を生み出していた、ということもありそうに思える。
Neu (1931) では、水面付近を泳ぐペンギンの羽ばたき運動を調べていたが、これは今回の結果と異なる。水面付近ならではの制約があったのではと推測している。

The energetics of 'flying' and 'paddling' in water: locomotion in penguins and ducks

体重が同じくらいのアヒルとリトルペンギンを比較。アヒルは足漕ぎ遊泳 (paddling) の代表、ペンギンは翼での遊泳 (flying) の代表という設定。
メインは酸素（正確には体積流量 $\dot{V}_{O_2}$）の計測。
水面を遊泳中のペンギンを 200 fps のカメラ1台で横から撮影し、周波数を取得（アヒルは足漕ぎ、ペンギンは羽ばたきの周波数）。カメラは Milliken DBH-5 というアナログフィルムの高速度カメラ。
死骸を使って水面での抗力を計測（水中にすると不安定で計測不能だった）。

The potential costs of flipper-bands to penguins

フリッパーバンドはコストを増大させる、という話。

2019-06-06

Fluent UDF: CURRENT_TIME はコール時点の現在時刻を示すが、DEFINE_GRID_MOTION の引数として渡される time は「メッシュ移動後の時刻」を示している（ぽい）

ANSYS CFD Fluent troubleshooting

タイトルのとおり。まだ確証はないんだけどテストした結果おそらくこれ。確かに各時刻でメッシュ移動は流れ計算より先に行われるので、「流れ場にとっての時刻」はメッシュ移動「後」になる、と思えば理屈はわかる。しかしいつものとおりマニュアルでこのへんが明確には説明されていないので非常に紛らわしい…。

2019-06-06

Fluent UDF: NODE_POS_UPDATED で付与された MOVED フラグは NODE_POS_MOVABLE コマンドでリセットできる

Fluent CFD ANSYS C/C++ troubleshooting

問題編

DEFINE_GRID_MOTION を使って板とかを曲げたりしたいとき、メッシュ内の各 node は複数の face にまたがっている可能性があるので、UDF マニュアルのサンプルにあるように、

#include "udf.h"

DEFINE_GRID_MOTION(my_dyn_mesh, domain, dt, time, dtime)
{
  Thread *tp = DT_THREAD(dt); // tp = thread pointer
  face_t f; // face index
  Node *np; // np = node pointer

  begin_f_loop(f,tp)
  {
    f_node_loop(f,tp,n)
    {
      np = F_NODE(f,tp,n);

      if ( NODE_POS_NEED_UPDATE(np) )
      {
        NODE_POS_UPDATED(np); // Update the flag
        
        do_something;
      }
    }
  }
  end_f_loop(f,tp);
}

というように、

if( NODE_POS_NEED_UPDATE(np) ) で「既に訪れたかどうか」のフラグを調べて、
NODE_POS_UPDATED(np); コマンドで「訪れた」という情報に書き換える。

しかし、「計算開始直後など、時間進行を伴わずに全ノードをループして何かをしたい」ということがある（自分はあった）。このとき、何も考えずにループ構造をコピペしてしまうと、1回目のループで既にフラグが「既に訪れた」ことになっているので、2回目のループでは「全ての」node をスキップしてしまう。つまり何も起きない。この問題にしばらく悩んでいた。つまり、

DEFINE_GRID_MOTION(my_dyn_mesh, domain, dt, time, dtime)
{
  Thread *tp = DT_THREAD(dt); // tp = thread pointer
  face_t f; // face index
  Node *np; // np = node pointer

  // First loop
  begin_f_loop(f,tp)
  {
    f_node_loop(f,tp,n)
    {
      np = F_NODE(f,tp,n);

      if ( NODE_POS_NEED_UPDATE(np) )
      {
        NODE_POS_UPDATED(np);
        
        do_something;
      }
    }
  }
  end_f_loop(f,tp);

  // Another loop
  begin_f_loop(f,tp)
  {
    f_node_loop(f,tp,n)
    {
      np = F_NODE(f,tp,n);

      if ( NODE_POS_NEED_UPDATE(np) )
      {
        NODE_POS_UPDATED(np);
        
        do_something_else; // This doesn't work!
      }
    }
  }
  end_f_loop(f,tp);
}

ということだ。…と思っていた。調査編でわかるように、実際には単に「訪れたかどうか」ではなかったのだが…。

調査編

まず、フラグの実体はなんなのか。適当なフォルダ（C:\Program Files\ANSYS Inc\v191\fluent\fluent19.1.0\ とか）において、Devas で NODE_POS_NEED_UPDATE や NODE_POS_UPDATED を検索すると、dynamesh_tools.h というヘッダファイルにヒットする。

/* node marks */
#define NODE_NULL       0x0000
#define NODE_FIXED      0x0001
#define NODE_MOVABLE    0x0002
#define NODE_CONTACT    0x0004
#define NODE_MOVED      0x0008
#define NODE_ADAPT      0x0010
#define NODE_DEFORMED   0x0020
＜中略＞
#define NODE_POS_MOVABLE(v) \
  (NODE_CLEAR_MARK(v, NODE_MOVED | NODE_FIXED | NODE_ADAPT), \
   NODE_SET_MARK(v, NODE_MOVABLE))
#define NODE_POS_NEED_UPDATE(v) (NODE_MARK(v) & NODE_MOVABLE)
#define NODE_POS_UPDATED(v) \
  (NODE_CLEAR_MARK(v, NODE_MOVABLE | NODE_FIXED | NODE_ADAPT), \
   NODE_SET_MARK(v, NODE_MOVED))
#define NODE_POS_UPDATED_P(v) (NODE_MARK(v) & NODE_MOVED)
＜略＞

とある。まず、フラグは TURE/FALSE のような2値ではなく、複数の値をとる。それがソース冒頭の NODE_NULL, NODE_FIXED, ... というやつだ。そしてこのフラグに関連するコマンドは「現在のフラグ状況を調べるコマンド」と「フラグをセットするコマンド」の2種類に分かれる。ここでは MOVABLE と MOVED 関連のみを示したが、実際には他のフラグに対応するコマンドもある。

現在のフラグ状況を調べるコマンド： - NODE_POS_NEED_UPDATE： MOVABLE かどうかを調べる - NODE_POS_UPDATED_P： MOVED かどうかを調べる

フラグをセットするコマンド： - NODE_POS_MOVABLE： MOVED 等のフラグを CLEAR_MARK して、MOVABLE をセットする - NODE_POS_UPDATED： MOVABLE 等のフラグを CLEAR_MARK して、MOVED をセットする

となっている。調べるコマンドは基本的に suffix として _P が付くのだが、なぜか MOVABLE だけは MOVABLE_P ではなく NEED_UPDATE となっていて一貫性に欠ける。

解決編

したがって、解決策としては、最初のループの後でフラグリセット用のループを回して、そこで NODE_POS_MOVABLE を実行する、ということになる。つまり、

DEFINE_GRID_MOTION(my_dyn_mesh, domain, dt, time, dtime)
{
  Thread *tp = DT_THREAD(dt); // tp = thread pointer
  face_t f; // face index
  Node *np; // np = node pointer

  // First loop
  begin_f_loop(f,tp)
  {
    f_node_loop(f,tp,n)
    {
      np = F_NODE(f,tp,n);

      if ( NODE_POS_NEED_UPDATE(np) ) // Check the MOVED flag
      {
        NODE_POS_UPDATED(np); // Set MOVED flag, clear any other flags
        
        do_something;
      }
    }
  }
  end_f_loop(f,tp);

  // Loop for resetting the MOVED flag
  begin_f_loop(f,tp)
  {
    f_node_loop(f,tp,n)
    {
      np = F_NODE(f,tp,n);
      NODE_POS_MOVABLE(np); // Clear MOVED flag (and other flags if any), set the MOVABLE flag
    }
  }
  end_f_loop(f,tp);

  // Another loop
  begin_f_loop(f,tp)
  {
    f_node_loop(f,tp,n)
    {
      np = F_NODE(f,tp,n);

      if ( NODE_POS_NEED_UPDATE(np) ) // Check the MOVED flag
      {
        NODE_POS_UPDATED(np); // Set MOVED flag, clear any other flags
        
        do_something_else; // This works!
      }
    }
  }
  end_f_loop(f,tp);
}

ということだ。

ところで、こうやってループを複数書いてると、end_f_loop(f,tp); を忘れてハマりがちなので注意したい（ハマった）。

2019-06-01

Fluent で dynamic mesh の motion preview を TUI (journal) からやるには → solve/mesh-motion

ANSYS Fluent CFD troubleshooting

タイトルで完結。コマンド名に preview って入ってないから TUI マニュアルで検索かけても出てこない…

dynamicsoar's log

主に研究関係のメモ

Fluent: journal file を複数に分割して include するには→ file/read-journal を使う

まえがき

方法

具体例

smoothing_common.jou

smoothing_spring.jou

smoothing_spring.jou

smoothing_linelast.jou

main.jou

UDF: cell thread に含まれる face thread にアクセスする

Fluent UDF: マクロ末尾の _M1 は「ひとつ前の時刻」の意味だった

ペンギン遊泳・フリッパー関連の先行研究

論文リスト

Kinematics of swimming of penguins at the Detroit Zoo

The energetics of 'flying' and 'paddling' in water: locomotion in penguins and ducks

The potential costs of flipper-bands to penguins

Fluent UDF: CURRENT_TIME はコール時点の現在時刻を示すが、DEFINE_GRID_MOTION の引数として渡される time は「メッシュ移動後の時刻」を示している（ぽい）

Fluent UDF: NODE_POS_UPDATED で付与された MOVED フラグは NODE_POS_MOVABLE コマンドでリセットできる

問題編

調査編

解決編

Fluent で dynamic mesh の motion preview を TUI (journal) からやるには → solve/mesh-motion