Pythonを用いて格子ボルツマン法流体の数値計算

格子ボルツマン法って何でしょうか？

僕は詳しくありません（笑）

格子ボルツマン法が何かはさておき、シミュレーションしてから理解するという所謂「形から入る」ということでPythonを用いてシミュレーションのソースコードを走らせてみました。

こちらのURLを参考にしています。
格子ボルツマン法のソースコード入手先

目的
1. Pythonのインストール
2. 外部ライブラリのインストール
  - numpyのインストール
  - matplotlibのインストール
Pythonコードの修正（Python2系→Python3系）
- 修正その1
  - 修正その2

格子ボルツマン法のPythonコードを実行

結果

まとめ

目的

とりあえずPythonを用いて格子ボルツマン法で流体の数値計算を実行してみます。

>>こちらにドキュメントもありますので詳しく知りたい方はどうぞ。（英文だが）

やることは、

流れ

格子ボルツマン法のソースコード入手先
よりPythonのソースコードを入手。

↓

Pythonのソースコード（Python2系）をPyhon3系に修正

↓

外部モジュール（numpy、matplotlib)のインストール

↓

Pythonの実行

環境

Windows10
Python3.6.4
(コマンドプロンプトからpythonを起動)

Pythonのインストール

Pythonはこちらからインストールできます。
https://qiita.com/taiponrock/items/f574dd2cddf8851fb02c

外部ライブラリのインストール

本記事のソースコードをPythonで実行するには以下の２つの外部モジュール(numpy、matplotlib)がインストールされている必要があります。

以下にインストールするにあたってわかりやすいサイトを添付しておきますので、これを参考に外部モジュールをインストールすることでコマンドプロンプトからPythonの計算が実行できます。

numpyのインストール

https://qiita.com/kenichi-t/items/7319b3876ae7f0d3817c

matplotlibのインストール

https://qiita.com/Morio/items/d75159bac916174e7654

Pythonコードの修正（Python2系→Python3系）

Pythonのソースコードを走らせてみたいと思うのですが、http://physics.weber.edu/schroeder/fluids/のPythonコードは2系で書かれています。

最近のPythonはバージョンアップされているため、例えばPython3.6.4(Python3系)で実行するためには、少し修正が必要となります。

Python2系を扱ったことがないので、詳しくは知りませんが以下の修正を行うことでPython3系でも正常に計算実行できます。

修正その1

【修正前】

barrier[(height/2)-8:(height/2)+8, (height/2)] = True

1	barrier[(height/2)-8:(height/2)+8, (height/2)] = True

これでは配列部分が実数になるので、以下のように整数型になるように修正する必要があります。

【修正後】

barrier[int(height/2)-8:int(height/2)+8, int(height/2)] = True

1	barrier[int(height/2)-8:int(height/2)+8, int(height/2)] = True

修正その2

【修正前】

print "%1.1f" % (100/(endTime-startTime)), 'frames per second'

1	print "%1.1f" % (100/(endTime-startTime)), 'frames per second'

これはpython2系の書き方である。python3系では以下ののように（）でくくる必要があります。

【修正後】

print ("%1.1f" % (100/(endTime-startTime)), 'frames per second')

1	print ("%1.1f" % (100/(endTime-startTime)), 'frames per second')

修正が終わったので、コマンドプロンプトを起動して、以下をタイプして計算実行してみましょう。
※ファイルの保存の際は、文字コードはUTF-8を選択

格子ボルツマン法のPythonコードを実行

コマンドプロンプトを立ち上げて以下を打ち込みます。

python Lattice_Boltzman.py

1	python Lattice_Boltzman.py

サンプルコードは以下に示します。（余計なコメントが多かったので、シンプルにしました。）

#ファイル名：Lattice_Boltzman.py
import numpy
import time
import matplotlib.pyplot
import matplotlib.animation
%matplotlib nbagg
#↑jupyter notebookを使っている方はこれを書いてください。

# Define constants:
height = 80 # lattice dimensions
width = 200
viscosity = 0.02 # fluid viscosity
omega = 1 / (3*viscosity + 0.5) # "relaxation" parameter
u0 = 0.1 # initial and in-flow speed
four9ths = 4.0/9.0 # abbreviations for lattice-Boltzmann weight factors
one9th = 1.0/9.0
one36th = 1.0/36.0
performanceData = False # set to True if performance data is desired

# Initialize all the arrays to steady rightward flow:
n0 = four9ths * (numpy.ones((height,width)) - 1.5*u0**2) # particle densities along 9 directions
nN = one9th * (numpy.ones((height,width)) - 1.5*u0**2)
nS = one9th * (numpy.ones((height,width)) - 1.5*u0**2)
nE = one9th * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2)
nW = one9th * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2)
nNE = one36th * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2)
nSE = one36th * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2)
nNW = one36th * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2)
nSW = one36th * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2)
rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW # macroscopic density
ux = (nE + nNE + nSE - nW - nNW - nSW) / rho # macroscopic x velocity
uy = (nN + nNE + nNW - nS - nSE - nSW) / rho # macroscopic y velocity

# Initialize barriers:
barrier = numpy.zeros((height,width), bool) # True wherever there's a barrier
barrier[int(height/2)-8:int(height/2)+8, int(height/2)] = True # simple linear barrier
barrierN = numpy.roll(barrier, 1, axis=0) # sites just north of barriers
barrierS = numpy.roll(barrier, -1, axis=0) # sites just south of barriers
barrierE = numpy.roll(barrier, 1, axis=1) # etc.
barrierW = numpy.roll(barrier, -1, axis=1)
barrierNE = numpy.roll(barrierN, 1, axis=1)
barrierNW = numpy.roll(barrierN, -1, axis=1)
barrierSE = numpy.roll(barrierS, 1, axis=1)
barrierSW = numpy.roll(barrierS, -1, axis=1)

# Move all particles by one step along their directions of motion (pbc):
def stream():
    global nN, nS, nE, nW, nNE, nNW, nSE, nSW
    nN = numpy.roll(nN, 1, axis=0) # axis 0 is north-south; + direction is north
    nNE = numpy.roll(nNE, 1, axis=0)
    nNW = numpy.roll(nNW, 1, axis=0)
    nS = numpy.roll(nS, -1, axis=0)
    nSE = numpy.roll(nSE, -1, axis=0)
    nSW = numpy.roll(nSW, -1, axis=0)
    nE = numpy.roll(nE, 1, axis=1) # axis 1 is east-west; + direction is east
    nNE = numpy.roll(nNE, 1, axis=1)
    nSE = numpy.roll(nSE, 1, axis=1)
    nW = numpy.roll(nW, -1, axis=1)
    nNW = numpy.roll(nNW, -1, axis=1)
    nSW = numpy.roll(nSW, -1, axis=1)
# Use tricky boolean arrays to handle barrier collisions (bounce-back):
    nN[barrierN] = nS[barrier]
    nS[barrierS] = nN[barrier]
    nE[barrierE] = nW[barrier]
    nW[barrierW] = nE[barrier]
    nNE[barrierNE] = nSW[barrier]
    nNW[barrierNW] = nSE[barrier]
    nSE[barrierSE] = nNW[barrier]
    nSW[barrierSW] = nNE[barrier]

# Collide particles within each cell to redistribute velocities (could be optimized a little more):
def collide():
    global rho, ux, uy, n0, nN, nS, nE, nW, nNE, nNW, nSE, nSW
    rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW
    ux = (nE + nNE + nSE - nW - nNW - nSW) / rho
    uy = (nN + nNE + nNW - nS - nSE - nSW) / rho
    ux2 = ux * ux # pre-compute terms used repeatedly...
    uy2 = uy * uy
    u2 = ux2 + uy2
    omu215 = 1 - 1.5*u2 # "one minus u2 times 1.5"
    uxuy = ux * uy
    n0 = (1-omega)*n0 + omega * four9ths * rho * omu215
    nN = (1-omega)*nN + omega * one9th * rho * (omu215 + 3*uy + 4.5*uy2)
    nS = (1-omega)*nS + omega * one9th * rho * (omu215 - 3*uy + 4.5*uy2)
    nE = (1-omega)*nE + omega * one9th * rho * (omu215 + 3*ux + 4.5*ux2)
    nW = (1-omega)*nW + omega * one9th * rho * (omu215 - 3*ux + 4.5*ux2)
    nNE = (1-omega)*nNE + omega * one36th * rho * (omu215 + 3*(ux+uy) + 4.5*(u2+2*uxuy))
    nNW = (1-omega)*nNW + omega * one36th * rho * (omu215 + 3*(-ux+uy) + 4.5*(u2-2*uxuy))
    nSE = (1-omega)*nSE + omega * one36th * rho * (omu215 + 3*(ux-uy) + 4.5*(u2-2*uxuy))
    nSW = (1-omega)*nSW + omega * one36th * rho * (omu215 + 3*(-ux-uy) + 4.5*(u2+2*uxuy))
# Force steady rightward flow at ends (no need to set 0, N, and S components):
    nE[:,0] = one9th * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2)
    nW[:,0] = one9th * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2)
    nNE[:,0] = one36th * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2)
    nSE[:,0] = one36th * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2)
    nNW[:,0] = one36th * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2)
    nSW[:,0] = one36th * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2)

# Compute curl of the macroscopic velocity field:
def curl(ux, uy):
    return numpy.roll(uy,-1,axis=1) - numpy.roll(uy,1,axis=1) - numpy.roll(ux,-1,axis=0) + numpy.roll(ux,1,axis=0)

# Here comes the graphics and animation...
theFig = matplotlib.pyplot.figure(figsize=(8,3))
fluidImage = matplotlib.pyplot.imshow(curl(ux, uy), origin='lower', norm=matplotlib.pyplot.Normalize(-.1,.1), 
cmap=matplotlib.pyplot.get_cmap('jet'), interpolation='none')
# See http://www.loria.fr/~rougier/teaching/matplotlib/#colormaps for other cmap options
bImageArray = numpy.zeros((height, width, 4), numpy.uint8) # an RGBA image
bImageArray[barrier,3] = 255 # set alpha=255 only at barrier sites
barrierImage = matplotlib.pyplot.imshow(bImageArray, origin='lower', interpolation='none')

# Function called for each successive animation frame:
startTime = time.clock()
#frameList = open('frameList.txt','w') # file containing list of images (to make movie)
def nextFrame(arg): # (arg is the frame number, which we don't need)
    global startTime
    if performanceData and (arg%100 == 0) and (arg > 0):
        endTime = time.clock()
        print ("%1.1f" % (100/(endTime-startTime)), 'frames per second')
        startTime = endTime
#frameName = "frame%04d.png" % arg
#matplotlib.pyplot.savefig(frameName)
#frameList.write(frameName + '\n')
    for step in range(20): # adjust number of steps for smooth animation
        stream()
        collide()
    fluidImage.set_array(curl(ux, uy))
    return (fluidImage, barrierImage) # return the figure elements to redraw

animate = matplotlib.animation.FuncAnimation(theFig, nextFrame, interval=1, blit=True)
matplotlib.pyplot.show()

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

#ファイル名：Lattice_Boltzman.py

import numpy

import time

import matplotlib.pyplot

import matplotlib.animation

%matplotlib nbagg

#↑jupyter notebookを使っている方はこれを書いてください。

# Define constants:

height = 80 # lattice dimensions

width = 200

viscosity = 0.02 # fluid viscosity

omega = 1 / (3*viscosity + 0.5) # "relaxation" parameter

u0 = 0.1 # initial and in-flow speed

four9ths = 4.0/9.0 # abbreviations for lattice-Boltzmann weight factors

one9th = 1.0/9.0

one36th = 1.0/36.0

performanceData = False # set to True if performance data is desired

# Initialize all the arrays to steady rightward flow:

n0 = four9ths * (numpy.ones((height,width)) - 1.5*u0**2) # particle densities along 9 directions

nN = one9th * (numpy.ones((height,width)) - 1.5*u0**2)

nS = one9th * (numpy.ones((height,width)) - 1.5*u0**2)

nE = one9th * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nW = one9th * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nNE = one36th * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nSE = one36th * (numpy.ones((height,width)) + 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nNW = one36th * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nSW = one36th * (numpy.ones((height,width)) - 3*u0 + 4.5*u0**2 - 1.5*u0**2)

rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW # macroscopic density

ux = (nE + nNE + nSE - nW - nNW - nSW) / rho # macroscopic x velocity

uy = (nN + nNE + nNW - nS - nSE - nSW) / rho # macroscopic y velocity

# Initialize barriers:

barrier = numpy.zeros((height,width), bool) # True wherever there's a barrier

barrier[int(height/2)-8:int(height/2)+8, int(height/2)] = True # simple linear barrier

barrierN = numpy.roll(barrier, 1, axis=0) # sites just north of barriers

barrierS = numpy.roll(barrier, -1, axis=0) # sites just south of barriers

barrierE = numpy.roll(barrier, 1, axis=1) # etc.

barrierW = numpy.roll(barrier, -1, axis=1)

barrierNE = numpy.roll(barrierN, 1, axis=1)

barrierNW = numpy.roll(barrierN, -1, axis=1)

barrierSE = numpy.roll(barrierS, 1, axis=1)

barrierSW = numpy.roll(barrierS, -1, axis=1)

# Move all particles by one step along their directions of motion (pbc):

def stream():

global nN, nS, nE, nW, nNE, nNW, nSE, nSW

nN = numpy.roll(nN, 1, axis=0) # axis 0 is north-south; + direction is north

nNE = numpy.roll(nNE, 1, axis=0)

nNW = numpy.roll(nNW, 1, axis=0)

nS = numpy.roll(nS, -1, axis=0)

nSE = numpy.roll(nSE, -1, axis=0)

nSW = numpy.roll(nSW, -1, axis=0)

nE = numpy.roll(nE, 1, axis=1) # axis 1 is east-west; + direction is east

nNE = numpy.roll(nNE, 1, axis=1)

nSE = numpy.roll(nSE, 1, axis=1)

nW = numpy.roll(nW, -1, axis=1)

nNW = numpy.roll(nNW, -1, axis=1)

nSW = numpy.roll(nSW, -1, axis=1)

# Use tricky boolean arrays to handle barrier collisions (bounce-back):

nN[barrierN] = nS[barrier]

nS[barrierS] = nN[barrier]

nE[barrierE] = nW[barrier]

nW[barrierW] = nE[barrier]

nNE[barrierNE] = nSW[barrier]

nNW[barrierNW] = nSE[barrier]

nSE[barrierSE] = nNW[barrier]

nSW[barrierSW] = nNE[barrier]

# Collide particles within each cell to redistribute velocities (could be optimized a little more):

def collide():

global rho, ux, uy, n0, nN, nS, nE, nW, nNE, nNW, nSE, nSW

rho = n0 + nN + nS + nE + nW + nNE + nSE + nNW + nSW

ux = (nE + nNE + nSE - nW - nNW - nSW) / rho

uy = (nN + nNE + nNW - nS - nSE - nSW) / rho

ux2 = ux * ux # pre-compute terms used repeatedly...

uy2 = uy * uy

u2 = ux2 + uy2

omu215 = 1 - 1.5*u2 # "one minus u2 times 1.5"

uxuy = ux * uy

n0 = (1-omega)*n0 + omega * four9ths * rho * omu215

nN = (1-omega)*nN + omega * one9th * rho * (omu215 + 3*uy + 4.5*uy2)

nS = (1-omega)*nS + omega * one9th * rho * (omu215 - 3*uy + 4.5*uy2)

nE = (1-omega)*nE + omega * one9th * rho * (omu215 + 3*ux + 4.5*ux2)

nW = (1-omega)*nW + omega * one9th * rho * (omu215 - 3*ux + 4.5*ux2)

nNE = (1-omega)*nNE + omega * one36th * rho * (omu215 + 3*(ux+uy) + 4.5*(u2+2*uxuy))

nNW = (1-omega)*nNW + omega * one36th * rho * (omu215 + 3*(-ux+uy) + 4.5*(u2-2*uxuy))

nSE = (1-omega)*nSE + omega * one36th * rho * (omu215 + 3*(ux-uy) + 4.5*(u2-2*uxuy))

nSW = (1-omega)*nSW + omega * one36th * rho * (omu215 + 3*(-ux-uy) + 4.5*(u2+2*uxuy))

# Force steady rightward flow at ends (no need to set 0, N, and S components):

nE[:,0] = one9th * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nW[:,0] = one9th * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nNE[:,0] = one36th * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nSE[:,0] = one36th * (1 + 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nNW[:,0] = one36th * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2)

nSW[:,0] = one36th * (1 - 3*u0 + 4.5*u0**2 - 1.5*u0**2)

# Compute curl of the macroscopic velocity field:

def curl(ux, uy):

return numpy.roll(uy,-1,axis=1) - numpy.roll(uy,1,axis=1) - numpy.roll(ux,-1,axis=0) + numpy.roll(ux,1,axis=0)

# Here comes the graphics and animation...

theFig = matplotlib.pyplot.figure(figsize=(8,3))

fluidImage = matplotlib.pyplot.imshow(curl(ux, uy), origin='lower', norm=matplotlib.pyplot.Normalize(-.1,.1),

cmap=matplotlib.pyplot.get_cmap('jet'), interpolation='none')

# See http://www.loria.fr/~rougier/teaching/matplotlib/#colormaps for other cmap options

bImageArray = numpy.zeros((height, width, 4), numpy.uint8) # an RGBA image

bImageArray[barrier,3] = 255 # set alpha=255 only at barrier sites

barrierImage = matplotlib.pyplot.imshow(bImageArray, origin='lower', interpolation='none')

# Function called for each successive animation frame:

startTime = time.clock()

#frameList = open('frameList.txt','w') # file containing list of images (to make movie)

def nextFrame(arg): # (arg is the frame number, which we don't need)

global startTime

if performanceData and (arg%100 == 0) and (arg > 0):

endTime = time.clock()

print ("%1.1f" % (100/(endTime-startTime)), 'frames per second')

startTime = endTime

#frameName = "frame%04d.png" % arg

#matplotlib.pyplot.savefig(frameName)

#frameList.write(frameName + '\n')

for step in range(20): # adjust number of steps for smooth animation

stream()

collide()

fluidImage.set_array(curl(ux, uy))

return (fluidImage, barrierImage) # return the figure elements to redraw

animate = matplotlib.animation.FuncAnimation(theFig, nextFrame, interval=1, blit=True)

matplotlib.pyplot.show()

結果

まとめ

雛形のコードを修正するだけですぐに実行できました。

格子ボルツマン法の理解はまだできていないので、ドキュメントを読んで別途投稿しようと思います。

速度空間という考え方がキーになると思われます。

また、もともとはセルオートマトン法という考え方から来ていると思われますが勉強しないとわかりません(- -)

セルオートマトン法―複雑系の自己組織化と超並列処理

加藤恭義

Amazon 楽天市場 Yahoo

Amazonの情報を掲載しています

目的

環境

Pythonのインストール

外部ライブラリのインストール

numpyのインストール

matplotlibのインストール

Pythonコードの修正（Python2系→Python3系）

修正その1

修正その2

格子ボルツマン法のPythonコードを実行

結果

まとめ

エンタルピーを用いてベルヌーイの定理を記述

流体力学 任意の面での応力ベクトルとコーシー応力の関係式【図でわかりやすく解説】

流体の抵抗が速度に比例する場合と速度の2乗に比例する場合【京都大学の物理、第Ⅰ問：力学編】

圧縮性流体のベルヌーイの定理

タンクからの噴出速度をベルヌーイの定理から計算

流体力学の非圧縮性と圧縮性

1次元のオイラーの運動方程式

流体力学 エネルギー保存則の導出

流体力学 循環の定義【渦ができる要因を理解する】

流体力学 ケルビンの循環定理【渦ができる要因を理解する】

等エントロピー流れでのベルヌーイの定理の計算

非圧縮性流体の基礎方程式

流体の密度伝搬は音速で伝わる。圧縮性の考慮は音速と流速の大小関係が重要

亜音速流れと超音速流れ

流体力学の力学的相似性

流体力学の運動量保存則の導出

流体力学 ラグランジュの渦定理【渦ができる要因を理解する】

流体力学の理解に必要な数学の知識をまとめた

流体力学 質量保存則（連続の式）の導出

流体力学 ヘルムホルツの渦定理【渦ができる要因を理解する】

COMMENT コメントをキャンセル

流体力学任意の面での応力ベクトルとコーシー応力の関係式【図でわかりやすく解説】

流体力学エネルギー保存則の導出

流体力学循環の定義【渦ができる要因を理解する】

流体力学ケルビンの循環定理【渦ができる要因を理解する】

流体力学ラグランジュの渦定理【渦ができる要因を理解する】

流体力学質量保存則（連続の式）の導出

流体力学ヘルムホルツの渦定理【渦ができる要因を理解する】