Python/PIVの検証用に管内の粒子流れ動画を作ってみた

粒子をトラッキングして流れの速度場を計測するPIVというのを学びたく、まずはサンプルの管内流れ動画をPythonで作ってみました。これで良いのかどうかはわかりませんが、Pythonによる流れ動画の作り方を備忘録として残します。

こんにちは。wat(@watlablog)です。本日はPIVコードのテストに使う粒子流れ動画を生成する案を記録します！

目次（項目クリックでジャンプできます）

PIVコード検証用のサンプル動画

PIVとは？

PIV(Particle Image Velocimetry)とは、粒子画像流速測定法を意味し、流れ場にトレーサ粒子を配置させ、粒子の速度を計測する技術です。

トレーサ粒子は液体であればまさに粒子、気体であれば人体に無害な煙状の物体を使うことが一般的です。

以下にPIVで行われている事をざっくり図で示します。
PIVは高速度カメラ等で動画（複数の画像）を撮影します。各画像間で何秒経っているかは撮影時の条件から自明なので、各画像内のトレーサ粒子の移動量を知る事ができれば、速度を算出する事が可能です。

簡単そうに見えますが、実際の撮影画像から精度良くトレーサ粒子の移動量を抽出するためには様々なノウハウが必要です。餅は餅屋という事で以下にPIV専門の会社様のリンクを記載しておきます。

外部リンク：カトウ光研株式会社：PIVとは

流体分野の実験について

機械工学の世界では航空機、自動車、その他産業機械は屋外や水中等で高速で移動したり、様々な流体の中で使用される機械を扱います。これら流体力を受ける機械は、場合によっては本来の性能を出せないことが多々あります。

そのため製品形状に対して数値シミュレーションを使って設計を行ったりするのですが、流体の解析は適切な壁の設定、乱流モデル、境界条件、メッシュ形状…と様々な要因で結果が変わるので、検証目的としても実験による流れの把握は重要です。

また、流体シミュレーションを精度良くやろうと思った時には、場合によっては億単位のセル数を必要とし、スパコンで何日も計算することもあります。実機があるなら実験の方がはやいというケースもあるのがこの分野の特徴ですね。

このように、計算における科学的発展が目覚ましい現代においても、流体分野は（他の分野もそうですが）理論流体力学や数値流体力学と共に、実験流体力学も重要な役割を担っています。

サンプル動画を作ってみたきっかけ

この記事を書いた理由は、勉強のために自分でもPIVのコードを作って使ってみたかったからです。そのためにはコードを試す撮影画像が必要なのですが…個人で実験設備を揃えるのは至難のワザなので、粒子の流れ画像をなんとなく作ってしまおうと思いました。

この記事で作成するサンプル動画

この記事では以下のようなGIF動画を作成します。

黒い点の一つ一つがトレーサ粒子を模擬しており、初期配置は完全ランダムにしてみました。
また、粒子の大きさも時刻によってある範囲で変わるようにしており、数値調整次第では途中で消えたり出てきたりといった挙動もできるようにしてあります。

図の\(y\)軸中心の流速が速く、\(y=\pm 2\)では遅いという管内流れの特徴も模擬しています。

流体は粘性を持ちますので、壁と流体の摩擦から壁付近では流速が遅くなります。実際は管径や粘性、レイノルズ数によってこの速度分布が変わりますが、今回のプログラムではここまで表現していません。ご興味のある方はハーゲンポアズイユ流れとかでGoogle検索してみて下さい。

以下に速度分布を同時プロットした動画を示しますが、今回はscipyのnorm.pdf(確率密度関数)を使って速度分布を模擬しました。実際は管のシミュレーションとかして求めるのが良いのでしょうけど、まだ未熟なためそれは後日。

norm.pdfは「Pythonでガウス分布を持つノイズの作り方と調整方法」でも使っていますので、是非参考にしてみて下さい。

Pythonで粒子流れ動画を作るコード例

全コード

全コードを以下に示します。GIF動画を作るcreate_gif関数と粒子画像を作るgenerate_particle_movie関数で構成されています。詳細はコード内コメントに記載しましたのでご確認下さい。

generate_particle_movie関数はwx（横幅）, wy（高さ）, sf（粒子密度のスケールファクター）, frame（フレーム数）, center（流れ中心）, sigma（速度分布の標準偏差）, vsf（速度のスケールファクター）を引数にしています。

是非色々な流れ場を作って遊んでみて下さい。

import numpy as np
from scipy.stats import norm
from PIL import Image
import os
import glob
from matplotlib import pyplot as plt

# GIFアニメーションを作成
def create_gif():
    path_list = sorted(glob.glob(os.path.join(*['out', '*'])))  # ファイルパスをソートしてリストする
    imgs = []                                                   # 画像をappendするための空配列を定義

    # ファイルのフルパスからファイル名と拡張子を抽出
    for i in range(len(path_list)):
        img = Image.open(path_list[i])                          # 画像ファイルを1つずつ開く
        imgs.append(img)                                        # 画像をappendで配列に格納していく

    # appendした画像配列をGIFにする
    imgs[0].save('fluid.gif',
                 save_all=True, append_images=imgs[1:], optimize=False, duration=100, loop=0)

def generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=0, sigma=0.8, vsf=0.1):
    # 粒子を初期配置を設定(ランダム）
    X = []
    Y = []

    for i in range(500):
        x = np.random.uniform(-wx * frame, wx, int(50 * frame * sf))      # x軸を作成
        y = np.random.uniform(-wy * 10, wy * 10, int(50 * frame * sf))    # y軸を作成
        X.append(x)
        Y.append(y)

    # 各フレームで速度を足しながら画像を作っていくループ
    for k in range(frame):
        print(k)                                                          # 現在のフレームをコンソールに表示
        X += vsf * norm.pdf(Y, center, sigma)                             # 速度分布を定義して粒子のx座標を変更

        # ここからグラフ描画-------------------------------------
        # フォントの種類とサイズを設定する。
        plt.rcParams['font.size'] = 14
        plt.rcParams['font.family'] = 'Times New Roman'

        # 目盛を内側にする。
        plt.rcParams['xtick.direction'] = 'in'
        plt.rcParams['ytick.direction'] = 'in'

        # グラフの上下左右に目盛線を付ける。
        fig = plt.figure(figsize=(wx,wy))
        ax1 = fig.add_subplot(111)
        ax1.yaxis.set_ticks_position('both')
        ax1.xaxis.set_ticks_position('both')

        # 軸のラベルを設定する。
        ax1.set_xlabel('x')
        ax1.set_ylabel('y')

        # スケールの設定をする。
        ax1.set_xticks(np.arange(0, wx * 10, wx / 10))
        ax1.set_xlim(0, wx)
        ax1.set_ylim(-wy / 2, wy / 2)

        # データプロット
        for j in range(len(X)):
            size =np.random.randint(1,20)
            ax1.scatter(X[j], Y[j], label='particle', s=size, color='black')

        #速度分布確認用のプロット
        #v = 2*norm.pdf(sorted(Y[0]), 0, 0.8)
        #ax1.plot(v, sorted(Y[0]), color='red', lw=6)

        # レイアウト設定
        fig.tight_layout()

        # グラフを表示する。
        #plt.show()

        # outフォルダが無い時に新規作成
        if os.path.exists('out'):
            pass
        else:
            os.mkdir('out')

        # 画像保存パスを準備
        path = os.path.join(*['out', str("{:05}".format(k)) + '.png'])

        # 画像を保存する
        plt.savefig(path)

        # GIFアニメーションを作成する関数を実行する
        create_gif()
        plt.close()
        # ---------------------------------------------------

generate_particle_movie(wx=10, wy=4, sf=1, frame=50)

import numpy as np

from scipy.stats import norm

from PIL import Image

import os

import glob

from matplotlib import pyplot as plt

# GIFアニメーションを作成

def create_gif():

path_list = sorted(glob.glob(os.path.join(*['out', '*']))) # ファイルパスをソートしてリストする

imgs = [] # 画像をappendするための空配列を定義

# ファイルのフルパスからファイル名と拡張子を抽出

for i in range(len(path_list)):

img = Image.open(path_list[i]) # 画像ファイルを1つずつ開く

imgs.append(img) # 画像をappendで配列に格納していく

# appendした画像配列をGIFにする

imgs[0].save('fluid.gif',

save_all=True, append_images=imgs[1:], optimize=False, duration=100, loop=0)

def generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=0, sigma=0.8, vsf=0.1):

# 粒子を初期配置を設定(ランダム）

X = []

Y = []

for i in range(500):

x = np.random.uniform(-wx * frame, wx, int(50 * frame * sf)) # x軸を作成

y = np.random.uniform(-wy * 10, wy * 10, int(50 * frame * sf)) # y軸を作成

X.append(x)

Y.append(y)

# 各フレームで速度を足しながら画像を作っていくループ

for k in range(frame):

print(k) # 現在のフレームをコンソールに表示

X += vsf * norm.pdf(Y, center, sigma) # 速度分布を定義して粒子のx座標を変更

# ここからグラフ描画-------------------------------------

# フォントの種類とサイズを設定する。

plt.rcParams['font.size'] = 14

plt.rcParams['font.family'] = 'Times New Roman'

# 目盛を内側にする。

plt.rcParams['xtick.direction'] = 'in'

plt.rcParams['ytick.direction'] = 'in'

# グラフの上下左右に目盛線を付ける。

fig = plt.figure(figsize=(wx,wy))

ax1 = fig.add_subplot(111)

ax1.yaxis.set_ticks_position('both')

ax1.xaxis.set_ticks_position('both')

# 軸のラベルを設定する。

ax1.set_xlabel('x')

ax1.set_ylabel('y')

# スケールの設定をする。

ax1.set_xticks(np.arange(0, wx * 10, wx / 10))

ax1.set_xlim(0, wx)

ax1.set_ylim(-wy / 2, wy / 2)

# データプロット

for j in range(len(X)):

size =np.random.randint(1,20)

ax1.scatter(X[j], Y[j], label='particle', s=size, color='black')

#速度分布確認用のプロット

#v = 2*norm.pdf(sorted(Y[0]), 0, 0.8)

#ax1.plot(v, sorted(Y[0]), color='red', lw=6)

# レイアウト設定

fig.tight_layout()

# グラフを表示する。

#plt.show()

# outフォルダが無い時に新規作成

if os.path.exists('out'):

pass

else:

os.mkdir('out')

# 画像保存パスを準備

path = os.path.join(*['out', str("{:05}".format(k)) + '.png'])

# 画像を保存する

plt.savefig(path)

# GIFアニメーションを作成する関数を実行する

create_gif()

plt.close()

# ---------------------------------------------------

generate_particle_movie(wx=10, wy=4, sf=1, frame=50)

色々な流れ場を作ってみる

デフォルトの流れ

デフォルトは以下のコードです（関数実行時に何も指定しなくても動くようにしました）。

generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=0, sigma=0.8, vsf=0.1)

1	generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=0, sigma=0.8, vsf=0.1)

一様流

速度分布のsigmaとvsf(速度に対するスケールファクター)を変更して一様流を作ってみます。

generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=0, sigma=10, vsf=1)

1	generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=0, sigma=10, vsf=1)

グラフの上下端で速度変化が少なくなりました。流体シミュレーションで言えばスリップ境界になったようなイメージ。

センターオフセット流れ

流速の速い部分を\(y=1\)とオフセットさせたセンターオフセット流れ（ネーミングは一般的ではないと思いますが）を以下に示します。

コード変更点は以下。centerの値を0以外にします。ここでは流れの速さをわかりやすくするためにvsfも変えています。

generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=1, sigma=0.8, vsf=0.5)

1	generate_particle_movie(wx=10, wy=4, sf=1, frame=50, center=1, sigma=0.8, vsf=0.5)

軸を消して撮影画像のようにしてみる

画像解析で使うために、グラフの情報（軸の値、ラベル等）は要らないかも知れません。以下のコードのようにラベル定義をコメントアウトし、作ったtickの設定をFalseにする事で消す事ができます。

# 軸のラベルを設定する。
#ax1.set_xlabel('x')
#ax1.set_ylabel('y')

# 軸を消して撮影画像のようにする-----------
ax1.tick_params(labelbottom=False,
               labelleft=False,
               labelright=False,
               labeltop=False)

ax1.tick_params(bottom=False,
               left=False,
               right=False,
               top=False)
#------------------------------------

# 軸のラベルを設定する。

#ax1.set_xlabel('x')

#ax1.set_ylabel('y')

# 軸を消して撮影画像のようにする-----------

ax1.tick_params(labelbottom=False,

labelleft=False,

labelright=False,

labeltop=False)

ax1.tick_params(bottom=False,

left=False,

right=False,

top=False)

#------------------------------------

こんな感じ。

枠線すら要らない、という場合はもっと簡単で、.axis('off')を書けば良いだけです。

ax1.axis('off')

1	ax1.axis('off')

こんな感じ。これならそのまま解析できそうですね。

まとめ

PIVという流体の実験計測、分析技術をちょっと学んでみたく、まずはサンプルの流れ場の動画を作成する事をしてみました。

あまり需要の無い記事だとは思うので完全に自分用の備忘録のようなものです。

本記事ではPIVの概要と流れを模擬して動画を生成するPythonコードを紹介しました。

パラメータ変更で色々な流れ方をするようにできたと思います。

この動画で本当にPIV解析ができるのかどうかはまだ謎ですが、やってみて不具合があったらこのページを修正してみようと思います。

本当にPIV解析ができるようになるのか…は今後の勉強次第！
Twitterでも関連情報をつぶやいているので、wat(@watlablog)のフォローお待ちしています！

PIVコード検証用のサンプル動画

PIVとは？

流体分野の実験について

サンプル動画を作ってみたきっかけ

この記事で作成するサンプル動画

Pythonで粒子流れ動画を作るコード例

全コード

色々な流れ場を作ってみる

デフォルトの流れ

一様流

センターオフセット流れ

軸を消して撮影画像のようにしてみる

まとめ

コメント

コメントを残すコメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

PIVコード検証用のサンプル動画

PIVとは？

流体分野の実験について

サンプル動画を作ってみたきっかけ

この記事で作成するサンプル動画

Pythonで粒子流れ動画を作るコード例

全コード

色々な流れ場を作ってみる

デフォルトの流れ

一様流

センターオフセット流れ

軸を消して撮影画像のようにしてみる

まとめ

SNSでもご購読できます。

コメント

コメントを残す コメントをキャンセル

最近の投稿

アーカイブ

カテゴリー

コメントを残すコメントをキャンセル