【Python】Pythonで流体のトレーサーを可視化する方法｜工場プログラマーのモダン開発プラクティス

どうも、こんにちは。最近ようやく生活が落ち着いて、ブログ記事を書けるようになったコンです。

本記事では、numpy + matplotlib だけを使って、最小の「1個動く点」から始めて、段階的にトレーサーアニメーションを組み立てていく流れで解説します。

最終的に、こんな感じで白い丸が動く GIF を作れるようになります。

トレーサーはパラパラ漫画
Step 1：まずは1個の点を動かす
Step 2：複数の点に増やす
Step 3：速度場に沿って動かす（流体っぽくする）
Step 4：GIF として書き出す
まとめ

トレーサーはパラパラ漫画

「トレーサー」と聞くと何か特別な仕組みを想像するかもしれませんが、中身はとても単純で、要するにパラパラ漫画です。

まず、各粒子の (x, y) 座標を NumPy 配列で持っておき、毎フレームその数字を少しずつ書き換える。次に、更新した座標を matplotlib の scatter に「この位置に白い丸を描いてね」と渡すと、白い丸が並んだ静止画（1フレーム）が1枚出来上がります。同じことを 100 〜 200 枚ぶん繰り返して、最後に PillowWriter で連結すれば 1 個の GIF になり、人間の目には「丸が動いている」ように見える ── やっているのはこれだけです。

つまり、いわゆるトレーサーは 「座標の数字を動かす」「静止画を描く」「パラパラ漫画として連結する」 の3要素しかありません。以降のセクションでは、この3つを順番に組み立てていきます。

Step 1：まずは1個の点を動かす

最初の一歩として、白い丸を1個だけ、画面上を直線で動かすサンプルから始めます。matplotlib に同梱の FuncAnimation という機能を使うと、こんなに短く書けます。

import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
fig, ax = plt.subplots(figsize=(5, 5))
ax.set_xlim(-2, 2); ax.set_ylim(-2, 2); ax.set_aspect("equal")
sc = ax.scatter([0], [0], s=80, c="white", edgecolors="black", linewidths=0.6)
def update(frame):
    x = -1.8 + frame * 0.02
    sc.set_offsets([[x, 0]])
    return sc,
anim = FuncAnimation(fig, update, frames=200, interval=33, blit=True)
plt.show()

たったこれだけで、左から右に丸がスーッと動く動画が表示されます。

何をしているかというと：

scatter で1個の点（座標は [0, 0]）をプロット
update(frame) 関数を定義：フレーム番号から現在の x 座標を計算して、set_offsets で点の座標を上書き
FuncAnimation に「この update を frames 回（=200回）呼んで」と頼む

FuncAnimation は 「update() を frames 回呼んでくれるループの皮」と思っておけば OK です。中で何が動いているかというと、frame=0, 1, 2, ..., 199 と順番に update を呼び出して、各フレームの figure を画像として束ねている、それだけ。

もう1つのキー要素が scatter.set_offsets() です。これは 「scatter の中の点の座標を上書きする」関数で、「scatter プロットを毎回作り直す必要はなく、座標だけ更新すればいい」というのがポイント。これで matplotlib は描画を最小限の差分で済ませられます。

ちなみに blit=True は「変更があった部分だけ再描画する」最適化フラグで、これを有効にすると体感がだいぶ滑らかになります。逆に blit=False にするとフレーム全体を毎回描き直すので重い。特別な理由がなければ、アニメ系は基本 blit=True で良いと思います。

Step 2：複数の点に増やす

1個動かせれば、複数も同じ要領で動かせます。NumPy 配列で「N個の点の座標」を持って、set_offsets に N×2 の配列を渡すだけ。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
N = 100
rng = np.random.RandomState(0)
xs = rng.uniform(-2, 2, N)
ys = rng.uniform(-2, 2, N)
vxs = rng.uniform(-0.03, 0.03, N)  # 各点ランダムな速度
vys = rng.uniform(-0.03, 0.03, N)
fig, ax = plt.subplots(figsize=(5, 5))
ax.set_xlim(-2, 2); ax.set_ylim(-2, 2); ax.set_aspect("equal")
sc = ax.scatter(xs, ys, s=20, c="white", edgecolors="black", linewidths=0.4)
def update(frame):
    global xs, ys
    xs = xs + vxs
    ys = ys + vys
    sc.set_offsets(np.column_stack([xs, ys]))
    return sc,
anim = FuncAnimation(fig, update, frames=200, interval=33, blit=True)
plt.show()

これで100個の白い丸が、それぞれランダムな方向に等速でゆっくり流れていく動きが出ます。np.column_stack([xs, ys]) で N×2 の配列にしているのが set_offsets に渡す形式の決まり、というだけで、構造は Step 1 とまったく同じです。粒子数が 1 個から 100 個になっても、コードはほぼ同じ長さで済みます。

逆に言うと、ここから先で「速度をどう決めるか」だけが本質的な工夫ポイントになります。

Step 3：速度場に沿って動かす（流体っぽくする）

各点が独立にランダムな速度を持つだけだと「流体」っぽくはなりません。流体っぽくするには 「位置に応じて決まる速度」、つまり速度場に沿って動かす必要があります。

例として、よくある単一渦の場 \(\mathbf{v}(x, y) = (-y, x)\) を使います。原点中心にぐるぐる回るやつです。式は短いですが、なぜぐるぐる回るかというと、各点 \((x, y)\) における速度ベクトル \((-y, x)\) が、原点と結ぶ動径ベクトル \((x, y)\) と垂直になっているから。垂直方向に進み続けるので、結果として円軌道を描きます。

動きを1個ずつ追えるように、まずは粒子5個の最小例で書きます。ここを増やすと粒子の密度が上がるだけで、コード構造は変わりません。

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
N = 5  # まずは5粒子で動きを追う。慣れたら数百〜数千に増やせばよい
rng = np.random.RandomState(0)
xs = rng.uniform(-1.5, 1.5, N)
ys = rng.uniform(-1.5, 1.5, N)
DT = 0.04  # 1ステップで動かす最大距離
fig, ax = plt.subplots(figsize=(5, 5))
ax.set_xlim(-2, 2); ax.set_ylim(-2, 2); ax.set_aspect("equal")
sc = ax.scatter(xs, ys, s=120, c="white", edgecolors="black", linewidths=0.8)
def update(frame):
    global xs, ys
    # 単一渦 v = (-y, x) に沿って動かす
    u = -ys
    v = xs
    # 1ステップで動く距離が DT を超えないようクリップ
    mag = np.sqrt(u ** 2 + v ** 2) + 1e-12
    scale = np.minimum(1.0, DT / mag)
    xs = xs + u * scale
    ys = ys + v * scale
    sc.set_offsets(np.column_stack([xs, ys]))
    return sc,
anim = FuncAnimation(fig, update, frames=200, interval=33, blit=True)
plt.show()

このコードが内部で実際にどんな数字の更新をしているのか、粒子0番の座標を最初の5ステップだけトレースしてみると、こんな感じです

frame	開始 (x, y)	速度 (u, v) = (-y, x)	終了 (x, y)
0	(0.1464, 0.4377)	(-0.4377, 0.1464)	(0.1085, 0.4504)
1	(0.1085, 0.4504)	(-0.4504, 0.1085)	(0.0696, 0.4597)
2	(0.0696, 0.4597)	(-0.4597, 0.0696)	(0.0301, 0.4657)
3	(0.0301, 0.4657)	(-0.4657, 0.0301)	(-0.0098, 0.4683)
4	(-0.0098, 0.4683)	(-0.4683, -0.0098)	(-0.0498, 0.4675)

毎ステップやっているのは 「現在の (x, y) から速度 (-y, x) を計算 → スケールを掛けて移動 → 新しい (x, y) を上書き」 だけ。x が徐々に減って 0 を跨ぎ、y は緩やかに増えていく ── 円軌道の左半分を反時計回りに動いている、というのが数字からも見て取れます。

クリップを入れる理由はちょっとした実用上のテクニックで、場によっては局所的に速度がめちゃくちゃ大きくなる場所（点電荷の中心付近や、強い渦の中心）があり、そのまま xs += u すると粒子が一気に画面の外へ吹っ飛ぶことがあるからです。DT で「1ステップの最大移動距離」を頭打ちにすると、強い場でも安全にアニメが回ります。

u, v の計算式を別の場（点電荷の電場、双極子場、CFD で吐き出した数値データなど）に差し替えれば、構造をいじらずそのまま応用できます。「動く粒子」のロジックは速度場と完全に分離されている、というのが綺麗な点ですね。

Step 4：GIF として書き出す

ここまでは plt.show() で対話的に表示していましたが、ブログに貼ったり共有したりするには GIF として書き出したいところです。matplotlib 同梱の PillowWriter を使い、Step 3 のコード末尾の plt.show() をこう置き換えるだけです。

from matplotlib.animation import FuncAnimation, PillowWriter
# (省略：xs, ys, sc, update の定義は Step 3 と同じ)
anim = FuncAnimation(fig, update, frames=200, interval=33, blit=True)
anim.save("tracer.gif", writer=PillowWriter(fps=30))
plt.close(fig)

anim.save() は update() を 200 回呼んで 200 枚のフレーム画像を作り、PillowWriter がそれを 1個の GIF ファイルとして連結 しているだけです。fps=30 にしているので、200 フレーム ÷ 30 fps ≒ 6.7 秒の GIF になります。

GIF にする前に、200枚作られるフレームのうち5枚だけを抜き出して並べてみるとこんな感じです。グリッドと目盛を入れて、各粒子の座標がだいたい読めるようにしてあります。先ほど座標トレース表で追っていた粒子0番が青丸です。

青丸の粒子0が frame 0 → 50 → 100 → 150 → 199 と進むにつれて、原点の周りをぐるりと半周以上回っているのが目で追えます。「動いている」ように見える GIF も、実態はこういう静止画 200 枚の連結であり、これを 30 fps で連続再生すれば人の目には「滑らかに動いている粒子」に見える、というだけのことです。

これで 「数字を動かす → scatter で描く → FuncAnimation でループ → PillowWriter で GIF 化」 が一通り通せるようになりました。あとは速度場と粒子数と見た目をいじれば、好きなトレーサー動画が作れます。