2017年12月8日金曜日

Python + Matplotlib.pyplot でカルマン渦(気流)を連続的な流線で描く


結果図はこちら



This page shows how to plot air flow past a cylinder with continuous streamline or how to plot vecter field with continuous streamline.
The air flow i.e. vecter field is calculated using COMSOL Multiphysics.
Sorry for inconvinience but I cannot give you the data used to draw the figures shown in this page.
このページではベクトル場である円柱後方の気流、つまりカルマン渦を連続的な流線で描きます。
ベクトル場は COMSOL Multiphysics を用いて計算しました。
その気流のベクトル場は、申し訳ありませんがアップロードの仕方がわかりません。
This page is based on follwing web sites:
Flow Past a Cylinder - COMSOL APPLICATION GALLERY - URL: https://www.comsol.com/model/flow-past-a-cylinder-97
and
Visualizing streamlines - Number Crunch -
URL: https://www.numbercrunch.de/blog/2013/05/visualizing-streamlines/
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import pandas as pd
import scipy as sp
import numpy as np
from matplotlib import cm
from scipy import interpolate
from scipy.integrate import ode
In [2]:
# load the data including pressure and vecter fields of each point
# calculated and exported by COMSOL Multiphysics
df = pd.read_csv("./cylinder_flow.txt",delim_whitespace=True,skiprows=9,
           names=["x","y","p","u","v"])
In [3]:
df.head(10)
Out[3]:
x y p u v
0 0.000367 0.000067 0.478731 0.000997 -3.151339e-09
1 0.001100 0.000067 0.478685 0.000997 -9.454016e-09
2 0.001833 0.000067 0.478638 0.000997 -1.575669e-08
3 0.002567 0.000067 0.478592 0.000997 -2.205937e-08
4 0.003300 0.000067 0.478546 0.000996 -2.836205e-08
5 0.004033 0.000067 0.478500 0.000996 -3.466472e-08
6 0.004767 0.000067 0.478454 0.000996 -4.096740e-08
7 0.005500 0.000067 0.478412 0.000996 -1.122428e-07
8 0.006233 0.000067 0.478375 0.000997 -2.769915e-07
9 0.006967 0.000067 0.478339 0.000998 -4.405037e-07
In [4]:
# interpolate.interp2d cannot treat Nan.
# fill Nan with sufficiently small number
df = df.fillna(10*10**(-15))
In [5]:
# the lines and rows of the input array
xlen = int(df["x"].shape[0]**0.5)
ylen = int(df["x"].shape[0]**0.5)
In [6]:
# make unique array from mesh-gridded input array
x_row = df["x"][: xlen]
y_row = df["y"][::ylen]
In [7]:
# generate interpolate function
fps = interpolate.interp2d(x_row,y_row,df["p"].values.reshape(xlen,ylen))
fus = interpolate.interp2d(x_row,y_row,df["u"].values.reshape(xlen,ylen))
fvs = interpolate.interp2d(x_row,y_row,df["v"].values.reshape(xlen,ylen))
In [8]:
def datalim(array,v0,v1):
    array = np.array(array)
    array[np.where(array<v0)] = v0
    array[np.where(array>v1)] = v1
    return array
In [9]:
# limit the range of the data
c0,c1 = -1,1
XX,YY = np.meshgrid(x_row,y_row)
limited = df["p"].values.reshape(xlen,ylen)
limited = np.array(limited)
limited = datalim(limited,c0*0.9999,c1*0.9999)

# plot the pressure data
fig = plt.figure(facecolor="w",figsize=(10,2))
plt.contourf(XX,YY,limited,100,cmap=cm.jet)
plt.clim(c0,c1)

# draw the cylinder
circle = plt.Circle((0.2,0.2),0.05,facecolor="w",edgecolor="k")
plt.axes().add_artist(circle)

# set aspect ratio to equal
plt.axes().set_aspect("equal","datalim")
plt.xlim(0,2.2)
plt.ylim(0,0.4)

# adjust colorbar
fig_coord = [0.91,0.15,0.01,0.7]
cbar_ax = fig.add_axes(fig_coord)
cbar = plt.colorbar(cax=cbar_ax)
cbar.set_ticks(np.linspace(c0,c1,5))
cbar.set_ticklabels(np.linspace(c0,c1,5))
ticklabs = cbar.ax.get_yticklabels()
cbar.ax.set_yticklabels(ticklabs,ha='right')
cbar.ax.yaxis.set_tick_params(pad=22)  # your number may vary
cbar.set_label("Pressure [Pa]")
In [10]:
def Vect_dir(t,p,fx,fy):
    ex = fx(p[0],p[1])
    ey = fy(p[0],p[1])
    n = (ex**2+ey**2)**0.5
    if n == 0:
        raise ValueError("hoge"+str(p)+str(n))
        print("hoge")
    return [ex/n,ey/n]
In [11]:
# starting point of the continuous stream line.
places=[[0.01,y] for y in np.arange(0.01,0.4,0.02)]
In [12]:
R=0.001
dt=0.8*R

# calculation area
x0, x1=0, 2.2
y0, y1=0, 0.4
xs,ys = [],[]

# initial setup of ordinary differential equation
r=ode(Vect_dir)
r.set_integrator('vode')
r.set_f_params(fus,fvs)

# loop over field lines starting in different directions 
for p in places:
    x=[p[0]]
    y=[p[1]]
    r.set_initial_value([x[0], y[0]], 0)
    while r.successful():
        r.integrate(r.t+dt)
        x.append(r.y[0])
        y.append(r.y[1])
        # check if field line left drwaing area or ends in some charge
        if (not (x0<r.y[0] and r.y[0]<x1)) or (not (y0<r.y[1] and r.y[1]<y1)):
            break
    xs.append(x)
    ys.append(y)
In [13]:
# limit the range of the data
c0,c1 = -1,1
XX,YY = np.meshgrid(x_row,y_row)
limited = df["p"].values.reshape(xlen,ylen)
limited = np.array(limited)
limited = datalim(limited,c0*0.9999,c1*0.9999)

# plot the pressure data
fig = plt.figure(facecolor="w",figsize=(10,2))
plt.contourf(XX,YY,limited,100,cmap=cm.jet)
plt.clim(c0,c1)

# draw continuous stream lines
for x,y in zip(xs,ys):
    plt.plot(x,y,color="k")

# draw the cylinder
circle = plt.Circle((0.2,0.2),0.05,facecolor="w",edgecolor="k")
plt.axes().add_artist(circle)

# set aspect ratio to equal
plt.axes().set_aspect("equal","datalim")
plt.xlim(0,2.2)
plt.ylim(0,0.4)

# adjust colorbar
fig_coord = [0.91,0.15,0.01,0.7]
cbar_ax = fig.add_axes(fig_coord)
cbar = plt.colorbar(cax=cbar_ax)
cbar.set_ticks(np.linspace(c0,c1,5))
cbar.set_ticklabels(np.linspace(c0,c1,5))
ticklabs = cbar.ax.get_yticklabels()
cbar.ax.set_yticklabels(ticklabs,ha='right')
cbar.ax.yaxis.set_tick_params(pad=22)  # your number may vary
cbar.set_label("Pressure [Pa]")
plt.savefig("cylinder_flow.png",bbox_inches="tight",pad_inches=0.02,dpi=250)
ラベル: , ,