尊敬する「宇宙い入ったカマキリ」さんの
ブログをなぞり、Pythonを使い数値流体
解析を行います。
カマキリさんのURLはこちら↓↓↓
https://takun-physics.net/10186/
#6は「2次元ポアソン方程式をPythonで実装する」
に挑戦します。
①ポアソン方程式
$\frac{\partial^2 p}{\partial x^2} + \frac{\partial^2 p}{\partial y^2} = b(x,y)$
物理現象の例：点電価と電位の関係
$\frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = - \frac{\rho ( \vec{r})}{ \epsilon}$
$\phi $：電位（V）
$\rho ( \vec{r})$：点電荷（C）
②ポアソン方程式の離散化
離散化方法として、空間に対しては２次精度中心差分を用いる。
※2次精度中心差分については、こちらを参考にして下さい。
「位置xを位置番号i、位置yを位置番号j、時間tを計算のステップ数n」で
置き換えて整理する。
$\frac{p^{n}_{i+1,j} -2p^{n}_{i,j} + p^{n}_{i-1,j}}{\Delta x^2} + \frac{p^{n}_{i,j+1} -2p^{n}_{i,j} + p^{n}_{i,j-1}}{\Delta y^2} = b^{n}_{i,j}$
$p^{n}_{i,j}$で整理すると、
$p^{n}_{i,j}=\frac{(p^{n}_{i+1,j} + p^{n}_{i-1,j})\Delta y^2 + (p^{n}_{i,j+1} + p^{n}_{i,j-1})\Delta x^2 - b^{n}_{i,j} \Delta x^2 \Delta y^2}{2(\Delta x^2 + \Delta y^2)}・・・(1)$
③境界条件の設定
$0 \le x \le 2 $、$0 \le y \le 2 $のx-y平面を考える。
・x-y平面の縁において$p=0$、つまり、$p_{0,j},p_{i,0},p_{nx,j},p_{i,ny}=0$
とする。ただし、$0 \le i \le nx $、$0 \le j \le ny $
・$i = \frac{1}{4} nx, j = \frac{1}{4} nx$ において$b_{i,j} = 100$、
$i = \frac{3}{4} nx, j = \frac{3}{4} nx$ において$b_{i,j} = -100$ とする。
④ポアソン分布を$python$で実装する。
・境界条件
---
>>> nx = 50
>>> ny = 50
>>> nt = 100
>>> xmin = 0
>>> xmax = 2
>>> ymin = 0
>>> ymax = 2
>>>
>>> dx = (xmax - xmin) / (nx -1)
>>> print(dx)
0.04081632653061224
>>> type(dx)
<class 'float'>
>>> dy = (ymax - ymin) / (ny -1)
>>> print(dy)
0.04081632653061224
---
・初期状態
---
>>> import numpy as np
>>> p = np.zeros([ny,nx])
>>> print(p)
[[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
...
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]]
>>> type(p)
<class 'numpy.ndarray'>
>>> pd = np.zeros([ny,nx])
>>> print(pd)
[[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
...
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]]
>>> b = np.zeros([ny,nx])
>>> print(b)
[[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
...
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]]
>>> x = np.linspace(xmin, xmax, nx)
>>> print(x)
[0. 0.04081633 0.08163265 0.12244898 0.16326531 0.20408163
0.24489796 0.28571429 0.32653061 0.36734694 0.40816327 0.44897959
0.48979592 0.53061224 0.57142857 0.6122449 0.65306122 0.69387755
0.73469388 0.7755102 0.81632653 0.85714286 0.89795918 0.93877551
0.97959184 1.02040816 1.06122449 1.10204082 1.14285714 1.18367347
1.2244898 1.26530612 1.30612245 1.34693878 1.3877551 1.42857143
1.46938776 1.51020408 1.55102041 1.59183673 1.63265306 1.67346939
1.71428571 1.75510204 1.79591837 1.83673469 1.87755102 1.91836735
1.95918367 2. ]
>>> y = np.linspace(ymin, ymax, ny)
>>> print(y)
[0. 0.04081633 0.08163265 0.12244898 0.16326531 0.20408163
0.24489796 0.28571429 0.32653061 0.36734694 0.40816327 0.44897959
0.48979592 0.53061224 0.57142857 0.6122449 0.65306122 0.69387755
0.73469388 0.7755102 0.81632653 0.85714286 0.89795918 0.93877551
0.97959184 1.02040816 1.06122449 1.10204082 1.14285714 1.18367347
1.2244898 1.26530612 1.30612245 1.34693878 1.3877551 1.42857143
1.46938776 1.51020408 1.55102041 1.59183673 1.63265306 1.67346939
1.71428571 1.75510204 1.79591837 1.83673469 1.87755102 1.91836735
1.95918367 2. ]
---
*numpy.zeros(配列の形式)：要素が全て０の配列を生成
>>> b[int(ny /4), int(nx /4)] = 100
>>> b[int(3 * ny /4), int(3 * nx /4)] = -100
>>> print(b)
[[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
...
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]
[0. 0. 0. ... 0. 0. 0.]]
>>> print(b[50/4,50/4])
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
IndexError: only integers, slices (`:`), ellipsis (`...`), numpy.newaxis (`None`) and integer or boolean arrays are valid indices
>>> print(b[12,12])
100.0
>>> print(b[36,36])
0.0
>>> print(b[37,37])
-100.0
---
・(1)式を実装
$p^{n}_{i,j}=\frac{(p^{n}_{i+1,j} + p^{n}_{i-1,j})\Delta y^2 + (p^{n}_{i,j+1} + p^{n}_{i,j-1})\Delta x^2 - b^{n}_{i,j} \Delta x^2 \Delta y^2}{2(\Delta x^2 + \Delta y^2)}・・・(1)$
---
>>> for it in range(nt):
... pd = p.copy()
... p[1:-1,1:-1] = (((pd[1:-1,2:] + pd[1:-1, :-2]) * dy**2 + (pd[2:, 1:-1] + pd[:-2, 1:-1]) * dx**2 - b[1:-1, 1:-1] * dx**2 * dy**2) / (2 * (dx**2 + dy**2)))
...
>>> p[0, :] = 0
>>> p[ny-1, :] = 0
>>> p[:, 0] = 0
>>> p[:, nx-1] = 0
---
*q[1:-1,1:-1]：1番目から末尾の-1番目の行、1番目から末尾の-1番目の列を取り出す
*q[1:-1,2:]：1番目から末尾の-1番目の行、2番目以降の全ての列を取り出す
*q[1:-1,:-2]：1番目から末尾の-1番目の行、末尾の-2番目以前の列を取り出す
*pd = p.copy()：pのデータをpdにコピー(値渡し)する
**値渡しと参照渡し(pd = p)：参照渡しをすると、pの値が変わるとpdの値も連動してかわる
・p(x,y)をmatplotlibを使って可視化する
---
>>> from matplotlib import pyplot as plt
>>> from matplotlib import cm
>>> from mpl_toolkits.mplot3d import Axes3D
>>> def plot2D(x, y, p):
... fig = plt.figure(figsize=(11, 7), dpi=100)
... ax = fig.gca(projection = '3d')
... X, Y = np.meshgrid(x, y)
... surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)
... ax.view_init(30, 225)
... ax.set_xlabel('x')
... ax.set_ylabel('y')
... ax.set_zlabel('p')
...
>>> plot2D(x, y, p)
>>> plt.show()
---
・ｐの圧力分布
・b(x,y)をmatplotlibを使って可視化する
---
>>> def plot2D(x, y, b):
... fig = plt.figure(figsize=(11, 7), dpi=100)
... ax = fig.gca(projection = '3d')
... X, Y = np.meshgrid(x, y)
... surf = ax.plot_surface(X, Y, b[:], rstride=1, cstride=1, cmap=cm.viridis,linewidth=0, antialiased=False)
... ax.view_init(30, 225)
... ax.set_xlabel('x')
... ax.set_ylabel('y')
... ax.set_zlabel('b')
...
>>> plot2D(x, y, b)
>>> plt.show()
---
コメント
コメントを投稿