Draw flow past cylinder with continuous stream line using Python Matplotlib.pyplot


The result is:


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.
This page is based on follwing web sites:



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.000167 0.000167 -0.115589 0.000997 2.946692e-08
1 0.000500 0.000167 -0.115595 0.000997 8.839726e-08
2 0.000833 0.000167 -0.115601 0.000997 1.473229e-07
3 0.001167 0.000167 -0.115608 0.000997 2.062439e-07
4 0.001500 0.000167 -0.115614 0.000996 2.651603e-07
5 0.001833 0.000167 -0.115620 0.000996 3.240720e-07
6 0.002167 0.000167 -0.115627 0.000996 3.829790e-07
7 0.002500 0.000167 -0.115633 0.000996 4.418814e-07
8 0.002833 0.000167 -0.115639 0.000995 5.007791e-07
9 0.003167 0.000167 -0.115646 0.000995 5.596721e-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=(4,4.1))
plt.contourf(XX,YY,limited,100,cmap=cm.jet)
plt.clim(c0,c1)

# draw the cylinder
circle = plt.Circle((0.1,0.5),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,1)
plt.ylim(0,1)

# adjust colorbar
fig_coord = [0.93,0.15,0.02,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.show()

c:\users\ejipc\miniconda3\envs\py36\lib\site-packages\matplotlib\cbook\deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.
  warnings.warn(message, mplDeprecation, stacklevel=1)
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,1,0.02)]
In [12]:
R=0.001
dt=0.8*R

# calculation area
x0, x1=0, 1
y0, y1=0, 1
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=(4,4.1))
ax = fig.add_subplot(111)
cax = ax.contourf(XX,YY,limited,100,cmap=cm.jet,vmin=c0,vmax=c1)
#ax.set_clim(c0,c1)
# draw continuous stream lines
for x,y in zip(xs,ys):
    ax.plot(x,y,color="k")

# draw the cylinder
circle = plt.Circle((0.1,0.5),0.05,facecolor="w",edgecolor="k")
ax.add_artist(circle)

# set aspect ratio to equal
ax.set_aspect("equal","datalim")
ax.set_xlim(0,1)
ax.set_ylim(0,1)

# adjust colorbar
fig_coord = [0.93,0.15,0.02,0.7]
cbar_ax = fig.add_axes(fig_coord)
cbar = fig.colorbar(cax,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]")

# save the figure
plt.savefig("cylinder_flow.png",bbox_inches="tight",pad_inches=0.02,dpi=250)

Author: Date:
Labels: , , , ,