UDF: cell thread に含まれる face thread にアクセスする
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 ポインタのままで比較しても特に問題はなさそうだった。
…あれ？でもなんかこれあまり有意義ではないかも・・・
ペンギン遊泳・フリッパー関連の先行研究
論文リスト
- 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
- 計測はいくつかの手法で行っている。
- 時間計測
- 2台のSony 3210カメラ（情報が出てこないがおそらくアナログフィルム）により60 fpsで水路の中を泳ぐペンギンを撮影。2台使っているがこれはステレオ撮影のためではなく、単に広い領域を撮影するためのようだ。calibration とか three-dimension とかいう言葉は全く出てこないし、結果もない。また、水路の外から撮影しているにもかかわらず、ガラスと水による屈折の影響を全く考慮していない。
- 1台のカメラによる撮影から、体長で割った速度と、羽ばたき周波数を計測。
- 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
フリッパーバンドはコストを増大させる、という話。
Fluent UDF: CURRENT_TIME はコール時点の現在時刻を示すが、DEFINE_GRID_MOTION の引数として渡される time は「メッシュ移動後の時刻」を示している（ぽい）
タイトルのとおり。まだ確証はないんだけどテストした結果おそらくこれ。確かに各時刻でメッシュ移動は流れ計算より先に行われるので、「流れ場にとっての時刻」はメッシュ移動「後」になる、と思えば理屈はわかる。しかしいつものとおりマニュアルでこのへんが明確には説明されていないので非常に紛らわしい…。
Fluent UDF: NODE_POS_UPDATED で付与された MOVED フラグは NODE_POS_MOVABLE コマンドでリセットできる
問題編
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); を忘れてハマりがちなので注意したい（ハマった）。
Fluent で dynamic mesh の motion preview を TUI (journal) からやるには → solve/mesh-motion
タイトルで完結。コマンド名に preview って入ってないから TUI マニュアルで検索かけても出てこない…
ANSYS Meshing と Fluent (or Fluent Meshing) の表示する maximum skewness や minimum orthogonal quality が微妙に違う（原因不明）
再現手順
- ANSYS Meshing でメッシュを作成し、maximum skewness を確認する。今回は 0.944 くらいだった。まぁ < 0.95 だし、いけるかな、と。
- Fluent 用に .msh を export する（Workbench で Update とするのと同じことだが、自分はこちらを好む）。
- Fluent や Fluent Meshing（Fluent launcher で Meshing mode を ON）で読み込む。
- TUI の
/mesh/check-verbosityを 2 以上にする（デフォルトは 0）。 - Mesh quality をチェックする (TUI だと
mesh/quality)。すると Maximum Cell Equivolumen Skewness = 0.967 とか出てくる。なんだこれ。minimum orthogonal quality も違う。誤差ってレベルじゃねーぞ…
考えられる原因
おそらく
- 定義が微妙に違う
- 実は .msh に export するときに微妙に形状が変わっている
のどちらかだろう。後者だとさすがにひどいから前者だと思いたいが、もしそうだとしてもなんでだよ感がすごい。もともとは別の会社のアプリだったとかこっちは知ったことじゃない。統合して売ってるんだから整合性とってくれよ…
