900 lines
31 KiB
Markdown
900 lines
31 KiB
Markdown
|
+++
|
|||
|
|
date = 2021-07-10T08:23:47Z
|
|||
|
|
description = ""
|
|||
|
|
draft = false
|
|||
|
|
slug = "jeonsanyuceyeoghag-cfd-with-python-navier-stokes-equation"
|
|||
|
|
title = "[전산유체역학] CFD with Python (Navier-Stokes Equation)"
|
|||
|
|
|
|||
|
|
+++
|
|||
|
|
|
|||
|
|
|
|||
|
|
## 1-D Linear Convection
|
|||
|
|
|
|||
|
|
1차원 선형 열전도 방정식은 가장 심플하면서도 가장 기초적인 방정식입니다.
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0 $$
|
|||
|
|
|
|||
|
|
이 식을 오일러 방정식으로 변환하여 수치해석적으로 해를 구할 수 있도록 변환을 해줍니다.
|
|||
|
|
|
|||
|
|
$$ u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n) $$
|
|||
|
|
|
|||
|
|
이제 이 오일러 방정식을 파이썬으로 구현해봅니다.
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot
|
|||
|
|
import time, sys
|
|||
|
|
%matplotlib inline
|
|||
|
|
|
|||
|
|
nx = 41 # try changing this number from 41 to 81 and Run All ... what happens?
|
|||
|
|
dx = 2 / (nx-1)
|
|||
|
|
nt = 25 #nt is the number of timesteps we want to calculate
|
|||
|
|
dt = .025 #dt is the amount of time each timestep covers (delta t)
|
|||
|
|
c = 1 #assume wavespeed of c = 1
|
|||
|
|
|
|||
|
|
u = numpy.ones(nx) #numpy function ones()
|
|||
|
|
u[int(.5 / dx):int(1 / dx + 1)] = 2 #setting u = 2 between 0.5 and 1 as per our I.C.s
|
|||
|
|
|
|||
|
|
un = numpy.ones(nx) #initialize a temporary array
|
|||
|
|
|
|||
|
|
for n in range(nt): #loop for values of n from 0 to nt, so it will run nt times
|
|||
|
|
un = u.copy() ##copy the existing values of u into un
|
|||
|
|
for i in range(1, nx): ## you can try commenting this line and...
|
|||
|
|
#for i in range(nx): ## ... uncommenting this line and see what happens!
|
|||
|
|
u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
|
|||
|
|
|
|||
|
|
pyplot.plot(numpy.linspace(0, 2, nx), u);
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/ZbM8j/btq9fWovXzY/D1HOqkCqgkw9YLDpyMFxb1/img.png" >}}
|
|||
|
|
|
|||
|
|
## 1-D Convection Equation (Non-Linear)
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = 0 $$
|
|||
|
|
|
|||
|
|
$$ u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy # we're importing numpy
|
|||
|
|
from matplotlib import pyplot # and our 2D plotting library
|
|||
|
|
%matplotlib inline
|
|||
|
|
|
|||
|
|
|
|||
|
|
nx = 41
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
nt = 20 #nt is the number of timesteps we want to calculate
|
|||
|
|
dt = .025 #dt is the amount of time each timestep covers (delta t)
|
|||
|
|
|
|||
|
|
u = numpy.ones(nx) #as before, we initialize u with every value equal to 1.
|
|||
|
|
u[int(.5 / dx) : int(1 / dx + 1)] = 2 #then set u = 2 between 0.5 and 1 as per our I.C.s
|
|||
|
|
|
|||
|
|
un = numpy.ones(nx) #initialize our placeholder array un, to hold the time-stepped solution
|
|||
|
|
|
|||
|
|
for n in range(nt): #iterate through time
|
|||
|
|
un = u.copy() ##copy the existing values of u into un
|
|||
|
|
for i in range(1, nx): ##now we'll iterate through the u array
|
|||
|
|
|
|||
|
|
###This is the line from Step 1, copied exactly. Edit it for our new equation.
|
|||
|
|
###then uncomment it and run the cell to evaluate Step 2
|
|||
|
|
|
|||
|
|
###u[i] = un[i] - c * dt / dx * (un[i] - un[i-1])
|
|||
|
|
|
|||
|
|
|
|||
|
|
pyplot.plot(numpy.linspace(0, 2, nx), u) ##Plot the results
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
## 1-D Diffusion Equation
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t}= \nu \frac{\partial^2 u}{\partial x^2} $$
|
|||
|
|
|
|||
|
|
$$ u_{i}^{n+1}=u_{i}^{n}+\frac{\nu\Delta t}{\Delta x^2}(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}) $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy #loading our favorite library
|
|||
|
|
from matplotlib import pyplot #and the useful plotting library
|
|||
|
|
%matplotlib inline
|
|||
|
|
|
|||
|
|
nx = 41
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
nt = 20 #the number of timesteps we want to calculate
|
|||
|
|
nu = 0.3 #the value of viscosity
|
|||
|
|
sigma = .2 #sigma is a parameter, we'll learn more about it later
|
|||
|
|
dt = sigma * dx**2 / nu #dt is defined using sigma ... more later!
|
|||
|
|
|
|||
|
|
|
|||
|
|
u = numpy.ones(nx) #a numpy array with nx elements all equal to 1.
|
|||
|
|
u[int(.5 / dx):int(1 / dx + 1)] = 2 #setting u = 2 between 0.5 and 1 as per our I.C.s
|
|||
|
|
|
|||
|
|
un = numpy.ones(nx) #our placeholder array, un, to advance the solution in time
|
|||
|
|
|
|||
|
|
for n in range(nt): #iterate through time
|
|||
|
|
un = u.copy() ##copy the existing values of u into un
|
|||
|
|
for i in range(1, nx - 1):
|
|||
|
|
u[i] = un[i] + nu * dt / dx**2 * (un[i+1] - 2 * un[i] + un[i-1])
|
|||
|
|
|
|||
|
|
pyplot.plot(numpy.linspace(0, 2, nx), u);
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/oaKgN/btq9iA57vBH/weCFyYoImjFkasFiDMir3k/img.png" >}}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
## Burger's Equation
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial ^2u}{\partial x^2} $$
|
|||
|
|
|
|||
|
|
$$ u_i^{n+1} = u_i^n - u_i^n \frac{\Delta t}{\Delta x} (u_i^n - u_{i-1}^n) + \nu \frac{\Delta t}{\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n) $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
import sympy
|
|||
|
|
from sympy import init_printing
|
|||
|
|
from matplotlib import pyplot
|
|||
|
|
from sympy.utilities.lambdify import lambdify
|
|||
|
|
|
|||
|
|
%matplotlib inline
|
|||
|
|
init_printing(use_latex=True)
|
|||
|
|
|
|||
|
|
x, nu, t = sympy.symbols('x nu t')
|
|||
|
|
phi = (sympy.exp(-(x - 4 * t)**2 / (4 * nu * (t + 1))) +
|
|||
|
|
sympy.exp(-(x - 4 * t - 2 * sympy.pi)**2 / (4 * nu * (t + 1))))
|
|||
|
|
|
|||
|
|
phiprime = phi.diff(x)
|
|||
|
|
|
|||
|
|
u = -2 * nu * (phiprime / phi) + 4
|
|||
|
|
ufunc = lambdify((t, x, nu), u)
|
|||
|
|
|
|||
|
|
###variable declarations
|
|||
|
|
nx = 101
|
|||
|
|
nt = 100
|
|||
|
|
dx = 2 * numpy.pi / (nx - 1)
|
|||
|
|
nu = .07
|
|||
|
|
dt = dx * nu
|
|||
|
|
|
|||
|
|
x = numpy.linspace(0, 2 * numpy.pi, nx)
|
|||
|
|
un = numpy.empty(nx)
|
|||
|
|
t = 0
|
|||
|
|
|
|||
|
|
u = numpy.asarray([ufunc(t, x0, nu) for x0 in x])
|
|||
|
|
|
|||
|
|
for n in range(nt):
|
|||
|
|
un = u.copy()
|
|||
|
|
for i in range(1, nx-1):
|
|||
|
|
u[i] = un[i] - un[i] * dt / dx *(un[i] - un[i-1]) + nu * dt / dx**2 *\
|
|||
|
|
(un[i+1] - 2 * un[i] + un[i-1])
|
|||
|
|
u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu * dt / dx**2 *\
|
|||
|
|
(un[1] - 2 * un[0] + un[-2])
|
|||
|
|
u[-1] = u[0]
|
|||
|
|
|
|||
|
|
u_analytical = numpy.asarray([ufunc(nt * dt, xi, nu) for xi in x])
|
|||
|
|
|
|||
|
|
pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
pyplot.plot(x,u, marker='o', lw=2, label='Computational')
|
|||
|
|
pyplot.plot(x, u_analytical, label='Analytical')
|
|||
|
|
pyplot.xlim([0, 2 * numpy.pi])
|
|||
|
|
pyplot.ylim([0, 10])
|
|||
|
|
pyplot.legend();
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/cE7P8B/btq9dF2BVpt/J6GbNhRT4dX1nfB2GPkurK/img.png" >}}
|
|||
|
|
|
|||
|
|
## 2-D Linear Convection
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x} + c\frac{\partial u}{\partial y} = 0 $$
|
|||
|
|
|
|||
|
|
$$ u_{i,j}^{n+1} = u_{i,j}^n-c \frac{\Delta t}{\Delta x}(u_{i,j}^n-u_{i-1,j}^n)-c \frac{\Delta t}{\Delta y}(u_{i,j}^n-u_{i,j-1}^n) $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D ##New Library required for projected 3d plots
|
|||
|
|
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
%matplotlib inline
|
|||
|
|
|
|||
|
|
###variable declarations
|
|||
|
|
nx = 81
|
|||
|
|
ny = 81
|
|||
|
|
nt = 100
|
|||
|
|
c = 1
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
dy = 2 / (ny - 1)
|
|||
|
|
sigma = .2
|
|||
|
|
dt = sigma * dx
|
|||
|
|
|
|||
|
|
x = numpy.linspace(0, 2, nx)
|
|||
|
|
y = numpy.linspace(0, 2, ny)
|
|||
|
|
|
|||
|
|
u = numpy.ones((ny, nx)) ##create a 1xn vector of 1's
|
|||
|
|
un = numpy.ones((ny, nx)) ##
|
|||
|
|
|
|||
|
|
###Assign initial conditions
|
|||
|
|
|
|||
|
|
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
|
|||
|
|
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
|
|||
|
|
|
|||
|
|
###Plot Initial Condition
|
|||
|
|
##the figsize parameter can be used to produce different sized images
|
|||
|
|
fig = pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
ax = fig.gca(projection='3d')
|
|||
|
|
X, Y = numpy.meshgrid(x, y)
|
|||
|
|
surf = ax.plot_surface(X, Y, u[:], cmap=cm.viridis)
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/MH3sO/btq9fmIfvXs/Tbral2sgJxUQHfJgs1hG61/img.png" >}}
|
|||
|
|
|
|||
|
|
## 2-D Convection
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = 0 $$$$ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = 0 $$$$ u_{i,j}^{n+1} = u_{i,j}^n - u_{i,j} \frac{\Delta t}{\Delta x} (u_{i,j}^n-u_{i-1,j}^n) - v_{i,j}^n \frac{\Delta t}{\Delta y} (u_{i,j}^n-u_{i,j-1}^n) $$$$ v_{i,j}^{n+1} = v_{i,j}^n - u_{i,j} \frac{\Delta t}{\Delta x} (v_{i,j}^n-v_{i-1,j}^n) - v_{i,j}^n \frac{\Delta t}{\Delta y} (v_{i,j}^n-v_{i,j-1}^n) $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
import numpy
|
|||
|
|
%matplotlib inline
|
|||
|
|
|
|||
|
|
###variable declarations
|
|||
|
|
nx = 101
|
|||
|
|
ny = 101
|
|||
|
|
nt = 80
|
|||
|
|
c = 1
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
dy = 2 / (ny - 1)
|
|||
|
|
sigma = .2
|
|||
|
|
dt = sigma * dx
|
|||
|
|
|
|||
|
|
x = numpy.linspace(0, 2, nx)
|
|||
|
|
y = numpy.linspace(0, 2, ny)
|
|||
|
|
|
|||
|
|
u = numpy.ones((ny, nx)) ##create a 1xn vector of 1's
|
|||
|
|
v = numpy.ones((ny, nx))
|
|||
|
|
un = numpy.ones((ny, nx))
|
|||
|
|
vn = numpy.ones((ny, nx))
|
|||
|
|
|
|||
|
|
###Assign initial conditions
|
|||
|
|
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
|
|||
|
|
u[int(.5 / dy):int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2
|
|||
|
|
##set hat function I.C. : v(.5<=x<=1 && .5<=y<=1 ) is 2
|
|||
|
|
v[int(.5 / dy):int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2
|
|||
|
|
fig = pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
ax = fig.gca(projection='3d')
|
|||
|
|
X, Y = numpy.meshgrid(x, y)
|
|||
|
|
|
|||
|
|
ax.plot_surface(X, Y, u, cmap=cm.viridis, rstride=2, cstride=2)
|
|||
|
|
ax.set_xlabel('$x$')
|
|||
|
|
ax.set_ylabel('$y$');
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/bbiKmO/btq9hyAHR99/3KeBvxXPvCzXXYqrTlSj9k/img.png" >}}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
## 2D Diffusion
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t} = \nu \frac{\partial ^2 u}{\partial x^2} + \nu \frac{\partial ^2 u}{\partial y^2} $$$$ \begin{split}u_{i,j}^{n+1} = u_{i,j}^n &+ \frac{\nu \Delta t}{\Delta x^2}(u_{i+1,j}^n - 2 u_{i,j}^n + u_{i-1,j}^n) \\&+ \frac{\nu \Delta t}{\Delta y^2}(u_{i,j+1}^n-2 u_{i,j}^n + u_{i,j-1}^n)\end{split} $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots
|
|||
|
|
%matplotlib inline
|
|||
|
|
###variable declarations
|
|||
|
|
nx = 31
|
|||
|
|
ny = 31
|
|||
|
|
nt = 17
|
|||
|
|
nu = .05
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
dy = 2 / (ny - 1)
|
|||
|
|
sigma = .25
|
|||
|
|
dt = sigma * dx * dy / nu
|
|||
|
|
|
|||
|
|
x = numpy.linspace(0, 2, nx)
|
|||
|
|
y = numpy.linspace(0, 2, ny)
|
|||
|
|
|
|||
|
|
u = numpy.ones((ny, nx)) # create a 1xn vector of 1's
|
|||
|
|
un = numpy.ones((ny, nx))
|
|||
|
|
|
|||
|
|
###Assign initial conditions
|
|||
|
|
# set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
|
|||
|
|
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
|
|||
|
|
|
|||
|
|
###Run through nt timesteps
|
|||
|
|
def diffuse(nt):
|
|||
|
|
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
|
|||
|
|
|
|||
|
|
for n in range(nt + 1):
|
|||
|
|
un = u.copy()
|
|||
|
|
u[1:-1, 1:-1] = (un[1:-1,1:-1] +
|
|||
|
|
nu * dt / dx**2 *
|
|||
|
|
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
|
|||
|
|
nu * dt / dy**2 *
|
|||
|
|
(un[2:,1: -1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))
|
|||
|
|
u[0, :] = 1
|
|||
|
|
u[-1, :] = 1
|
|||
|
|
u[:, 0] = 1
|
|||
|
|
u[:, -1] = 1
|
|||
|
|
|
|||
|
|
|
|||
|
|
fig = pyplot.figure()
|
|||
|
|
ax = fig.gca(projection='3d')
|
|||
|
|
surf = ax.plot_surface(X, Y, u[:], rstride=1, cstride=1, cmap=cm.viridis,
|
|||
|
|
linewidth=0, antialiased=True)
|
|||
|
|
ax.set_zlim(1, 2.5)
|
|||
|
|
ax.set_xlabel('$x$')
|
|||
|
|
ax.set_ylabel('$y$');
|
|||
|
|
|
|||
|
|
diffuse(14)
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/eLwQEW/btq9e0ysgnx/YVHruuNPlMa6pJODaIGJdK/img.png" >}}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
## Burgers' Equation in 2D
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} = \nu \; \left(\frac{\partial ^2 u}{\partial x^2} + \frac{\partial ^2 u}{\partial y^2}\right) $$$$ \frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = \nu \; \left(\frac{\partial ^2 v}{\partial x^2} + \frac{\partial ^2 v}{\partial y^2}\right) $$$$ \begin{split}v_{i,j}^{n+1} = & v_{i,j}^n - \frac{\Delta t}{\Delta x} u_{i,j}^n (v_{i,j}^n - v_{i-1,j}^n) - \frac{\Delta t}{\Delta y} v_{i,j}^n (v_{i,j}^n - v_{i,j-1}^n) \\&+ \frac{\nu \Delta t}{\Delta x^2}(v_{i+1,j}^n-2v_{i,j}^n+v_{i-1,j}^n) + \frac{\nu \Delta t}{\Delta y^2} (v_{i,j+1}^n - 2v_{i,j}^n + v_{i,j-1}^n)\end{split} $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|||
|
|
%matplotlib inline
|
|||
|
|
###variable declarations
|
|||
|
|
nx = 41
|
|||
|
|
ny = 41
|
|||
|
|
nt = 120
|
|||
|
|
c = 1
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
dy = 2 / (ny - 1)
|
|||
|
|
sigma = .0009
|
|||
|
|
nu = 0.01
|
|||
|
|
dt = sigma * dx * dy / nu
|
|||
|
|
|
|||
|
|
|
|||
|
|
x = numpy.linspace(0, 2, nx)
|
|||
|
|
y = numpy.linspace(0, 2, ny)
|
|||
|
|
|
|||
|
|
u = numpy.ones((ny, nx)) # create a 1xn vector of 1's
|
|||
|
|
v = numpy.ones((ny, nx))
|
|||
|
|
un = numpy.ones((ny, nx))
|
|||
|
|
vn = numpy.ones((ny, nx))
|
|||
|
|
comb = numpy.ones((ny, nx))
|
|||
|
|
|
|||
|
|
###Assign initial conditions
|
|||
|
|
|
|||
|
|
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
|
|||
|
|
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
|
|||
|
|
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
|
|||
|
|
v[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
|
|||
|
|
###(plot ICs)
|
|||
|
|
for n in range(nt + 1): ##loop across number of time steps
|
|||
|
|
un = u.copy()
|
|||
|
|
vn = v.copy()
|
|||
|
|
|
|||
|
|
u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
|
|||
|
|
dt / dx * un[1:-1, 1:-1] *
|
|||
|
|
(un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
|
|||
|
|
dt / dy * vn[1:-1, 1:-1] *
|
|||
|
|
(un[1:-1, 1:-1] - un[0:-2, 1:-1]) +
|
|||
|
|
nu * dt / dx**2 *
|
|||
|
|
(un[1:-1,2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
|
|||
|
|
nu * dt / dy**2 *
|
|||
|
|
(un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1]))
|
|||
|
|
|
|||
|
|
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
|
|||
|
|
dt / dx * un[1:-1, 1:-1] *
|
|||
|
|
(vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
|
|||
|
|
dt / dy * vn[1:-1, 1:-1] *
|
|||
|
|
(vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) +
|
|||
|
|
nu * dt / dx**2 *
|
|||
|
|
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
|
|||
|
|
nu * dt / dy**2 *
|
|||
|
|
(vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1]))
|
|||
|
|
|
|||
|
|
u[0, :] = 1
|
|||
|
|
u[-1, :] = 1
|
|||
|
|
u[:, 0] = 1
|
|||
|
|
u[:, -1] = 1
|
|||
|
|
|
|||
|
|
v[0, :] = 1
|
|||
|
|
v[-1, :] = 1
|
|||
|
|
v[:, 0] = 1
|
|||
|
|
v[:, -1] = 1
|
|||
|
|
fig = pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
ax = fig.gca(projection='3d')
|
|||
|
|
X, Y = numpy.meshgrid(x, y)
|
|||
|
|
ax.plot_surface(X, Y, u, cmap=cm.viridis, rstride=1, cstride=1)
|
|||
|
|
ax.plot_surface(X, Y, v, cmap=cm.viridis, rstride=1, cstride=1)
|
|||
|
|
ax.set_xlabel('$x$')
|
|||
|
|
ax.set_ylabel('$y$');
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/PL2CD/btq9fcyC1VV/MC1B8I2YedaaCFr5Lr06KK/img.png" >}}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
## 2D Laplace Equation
|
|||
|
|
|
|||
|
|
$$ \frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = 0 $$$$ p_{i,j}^n = \frac{\Delta y^2(p_{i+1,j}^n+p_{i-1,j}^n)+\Delta x^2(p_{i,j+1}^n + p_{i,j-1}^n)}{2(\Delta x^2 + \Delta y^2)} $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|||
|
|
%matplotlib inline
|
|||
|
|
|
|||
|
|
def plot2D(x, y, p):
|
|||
|
|
fig = pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
ax = fig.gca(projection='3d')
|
|||
|
|
X, Y = numpy.meshgrid(x, y)
|
|||
|
|
surf = ax.plot_surface(X, Y, p[:], rstride=1, cstride=1, cmap=cm.viridis,
|
|||
|
|
linewidth=0, antialiased=False)
|
|||
|
|
ax.set_xlim(0, 2)
|
|||
|
|
ax.set_ylim(0, 1)
|
|||
|
|
ax.view_init(30, 225)
|
|||
|
|
ax.set_xlabel('$x$')
|
|||
|
|
ax.set_ylabel('$y$')
|
|||
|
|
|
|||
|
|
def laplace2d(p, y, dx, dy, l1norm_target):
|
|||
|
|
l1norm = 1
|
|||
|
|
pn = numpy.empty_like(p)
|
|||
|
|
|
|||
|
|
while l1norm > l1norm_target:
|
|||
|
|
pn = p.copy()
|
|||
|
|
p[1:-1, 1:-1] = ((dy**2 * (pn[1:-1, 2:] + pn[1:-1, 0:-2]) +
|
|||
|
|
dx**2 * (pn[2:, 1:-1] + pn[0:-2, 1:-1])) /
|
|||
|
|
(2 * (dx**2 + dy**2)))
|
|||
|
|
|
|||
|
|
p[:, 0] = 0 # p = 0 @ x = 0
|
|||
|
|
p[:, -1] = y # p = y @ x = 2
|
|||
|
|
p[0, :] = p[1, :] # dp/dy = 0 @ y = 0
|
|||
|
|
p[-1, :] = p[-2, :] # dp/dy = 0 @ y = 1
|
|||
|
|
l1norm = (numpy.sum(numpy.abs(p[:]) - numpy.abs(pn[:])) /
|
|||
|
|
numpy.sum(numpy.abs(pn[:])))
|
|||
|
|
|
|||
|
|
return p
|
|||
|
|
|
|||
|
|
nx = 31
|
|||
|
|
ny = 31
|
|||
|
|
c = 1
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
dy = 2 / (ny - 1)
|
|||
|
|
|
|||
|
|
|
|||
|
|
##initial conditions
|
|||
|
|
p = numpy.zeros((ny, nx)) # create a XxY vector of 0's
|
|||
|
|
|
|||
|
|
|
|||
|
|
##plotting aids
|
|||
|
|
x = numpy.linspace(0, 2, nx)
|
|||
|
|
y = numpy.linspace(0, 1, ny)
|
|||
|
|
|
|||
|
|
##boundary conditions
|
|||
|
|
p[:, 0] = 0 # p = 0 @ x = 0
|
|||
|
|
p[:, -1] = y # p = y @ x = 2
|
|||
|
|
p[0, :] = p[1, :] # dp/dy = 0 @ y = 0
|
|||
|
|
p[-1, :] = p[-2, :] # dp/dy = 0 @ y = 1
|
|||
|
|
|
|||
|
|
p = laplace2d(p, y, dx, dy, 1e-4)
|
|||
|
|
|
|||
|
|
plot2D(x, y, p)
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/bxvdGX/btq9goyjEC1/YkjgRDKkIZuAe2isKNsv60/img.png" >}}
|
|||
|
|
|
|||
|
|
|
|||
|
|
|
|||
|
|
## 2D Poisson Equation
|
|||
|
|
|
|||
|
|
$$ \frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b $$$$ p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)} $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|||
|
|
%matplotlib inline
|
|||
|
|
# Parameters
|
|||
|
|
nx = 50
|
|||
|
|
ny = 50
|
|||
|
|
nt = 100
|
|||
|
|
xmin = 0
|
|||
|
|
xmax = 2
|
|||
|
|
ymin = 0
|
|||
|
|
ymax = 1
|
|||
|
|
|
|||
|
|
dx = (xmax - xmin) / (nx - 1)
|
|||
|
|
dy = (ymax - ymin) / (ny - 1)
|
|||
|
|
|
|||
|
|
# Initialization
|
|||
|
|
p = numpy.zeros((ny, nx))
|
|||
|
|
pd = numpy.zeros((ny, nx))
|
|||
|
|
b = numpy.zeros((ny, nx))
|
|||
|
|
x = numpy.linspace(xmin, xmax, nx)
|
|||
|
|
y = numpy.linspace(xmin, xmax, ny)
|
|||
|
|
|
|||
|
|
# Source
|
|||
|
|
b[int(ny / 4), int(nx / 4)] = 100
|
|||
|
|
b[int(3 * ny / 4), int(3 * nx / 4)] = -100
|
|||
|
|
|
|||
|
|
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
|
|||
|
|
|
|||
|
|
def plot2D(x, y, p):
|
|||
|
|
fig = pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
ax = fig.gca(projection='3d')
|
|||
|
|
X, Y = numpy.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$')
|
|||
|
|
|
|||
|
|
plot2D(x, y, p)
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
##
|
|||
|
|
|
|||
|
|
## Cavity Flow with Navier–Stokes
|
|||
|
|
|
|||
|
|
$$ \frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v}=-\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v} $$$$ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right) $$$$ \frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2} = -\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y} \right) $$$$ \begin{split}p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2\left(\Delta x^2+\Delta y^2\right)} \\& -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \\& \times \left[\frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right)-\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} -2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x}-\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]\end{split} $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|||
|
|
%matplotlib inline
|
|||
|
|
nx = 41
|
|||
|
|
ny = 41
|
|||
|
|
nt = 500
|
|||
|
|
nit = 50
|
|||
|
|
c = 1
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
dy = 2 / (ny - 1)
|
|||
|
|
x = numpy.linspace(0, 2, nx)
|
|||
|
|
y = numpy.linspace(0, 2, ny)
|
|||
|
|
X, Y = numpy.meshgrid(x, y)
|
|||
|
|
|
|||
|
|
rho = 1
|
|||
|
|
nu = .1
|
|||
|
|
dt = .001
|
|||
|
|
|
|||
|
|
u = numpy.zeros((ny, nx))
|
|||
|
|
v = numpy.zeros((ny, nx))
|
|||
|
|
p = numpy.zeros((ny, nx))
|
|||
|
|
b = numpy.zeros((ny, nx))
|
|||
|
|
|
|||
|
|
def build_up_b(b, rho, dt, u, v, dx, dy):
|
|||
|
|
|
|||
|
|
b[1:-1, 1:-1] = (rho * (1 / dt *
|
|||
|
|
((u[1:-1, 2:] - u[1:-1, 0:-2]) /
|
|||
|
|
(2 * dx) + (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
|
|||
|
|
((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
|
|||
|
|
2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
|
|||
|
|
(v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
|
|||
|
|
((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2))
|
|||
|
|
|
|||
|
|
return b
|
|||
|
|
|
|||
|
|
def pressure_poisson(p, dx, dy, b):
|
|||
|
|
pn = numpy.empty_like(p)
|
|||
|
|
pn = p.copy()
|
|||
|
|
|
|||
|
|
for q in range(nit):
|
|||
|
|
pn = p.copy()
|
|||
|
|
p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
|
|||
|
|
(pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
|
|||
|
|
(2 * (dx**2 + dy**2)) -
|
|||
|
|
dx**2 * dy**2 / (2 * (dx**2 + dy**2)) *
|
|||
|
|
b[1:-1,1:-1])
|
|||
|
|
|
|||
|
|
p[:, -1] = p[:, -2] # dp/dx = 0 at x = 2
|
|||
|
|
p[0, :] = p[1, :] # dp/dy = 0 at y = 0
|
|||
|
|
p[:, 0] = p[:, 1] # dp/dx = 0 at x = 0
|
|||
|
|
p[-1, :] = 0 # p = 0 at y = 2
|
|||
|
|
|
|||
|
|
return p
|
|||
|
|
|
|||
|
|
def cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu):
|
|||
|
|
un = numpy.empty_like(u)
|
|||
|
|
vn = numpy.empty_like(v)
|
|||
|
|
b = numpy.zeros((ny, nx))
|
|||
|
|
|
|||
|
|
for n in range(nt):
|
|||
|
|
un = u.copy()
|
|||
|
|
vn = v.copy()
|
|||
|
|
|
|||
|
|
b = build_up_b(b, rho, dt, u, v, dx, dy)
|
|||
|
|
p = pressure_poisson(p, dx, dy, b)
|
|||
|
|
|
|||
|
|
u[1:-1, 1:-1] = (un[1:-1, 1:-1]-
|
|||
|
|
un[1:-1, 1:-1] * dt / dx *
|
|||
|
|
(un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
|
|||
|
|
vn[1:-1, 1:-1] * dt / dy *
|
|||
|
|
(un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
|
|||
|
|
dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])))
|
|||
|
|
|
|||
|
|
v[1:-1,1:-1] = (vn[1:-1, 1:-1] -
|
|||
|
|
un[1:-1, 1:-1] * dt / dx *
|
|||
|
|
(vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
|
|||
|
|
vn[1:-1, 1:-1] * dt / dy *
|
|||
|
|
(vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
|
|||
|
|
dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))
|
|||
|
|
|
|||
|
|
u[0, :] = 0
|
|||
|
|
u[:, 0] = 0
|
|||
|
|
u[:, -1] = 0
|
|||
|
|
u[-1, :] = 1 # set velocity on cavity lid equal to 1
|
|||
|
|
v[0, :] = 0
|
|||
|
|
v[-1, :] = 0
|
|||
|
|
v[:, 0] = 0
|
|||
|
|
v[:, -1] = 0
|
|||
|
|
|
|||
|
|
|
|||
|
|
return u, v, p
|
|||
|
|
|
|||
|
|
u = numpy.zeros((ny, nx))
|
|||
|
|
v = numpy.zeros((ny, nx))
|
|||
|
|
p = numpy.zeros((ny, nx))
|
|||
|
|
b = numpy.zeros((ny, nx))
|
|||
|
|
nt = 100
|
|||
|
|
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)
|
|||
|
|
fig = pyplot.figure(figsize=(11,7), dpi=100)
|
|||
|
|
# plotting the pressure field as a contour
|
|||
|
|
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
|
|||
|
|
pyplot.colorbar()
|
|||
|
|
# plotting the pressure field outlines
|
|||
|
|
pyplot.contour(X, Y, p, cmap=cm.viridis)
|
|||
|
|
# plotting velocity field
|
|||
|
|
pyplot.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2])
|
|||
|
|
pyplot.xlabel('X')
|
|||
|
|
pyplot.ylabel('Y');
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/VgtZK/btq9flvNFsA/4d03urU7VcLPRqzS5g40m1/img.png" >}}
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
u = numpy.zeros((ny, nx))
|
|||
|
|
v = numpy.zeros((ny, nx))
|
|||
|
|
p = numpy.zeros((ny, nx))
|
|||
|
|
b = numpy.zeros((ny, nx))
|
|||
|
|
nt = 700
|
|||
|
|
u, v, p = cavity_flow(nt, u, v, dt, dx, dy, p, rho, nu)
|
|||
|
|
fig = pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
|
|||
|
|
pyplot.colorbar()
|
|||
|
|
pyplot.contour(X, Y, p, cmap=cm.viridis)
|
|||
|
|
pyplot.quiver(X[::2, ::2], Y[::2, ::2], u[::2, ::2], v[::2, ::2])
|
|||
|
|
pyplot.xlabel('X')
|
|||
|
|
pyplot.ylabel('Y');
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/boYcRB/btq9fmasDP2/6w1UPPVU4mG7frDcjQtBIK/img.png" >}}
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
fig = pyplot.figure(figsize=(11, 7), dpi=100)
|
|||
|
|
pyplot.contourf(X, Y, p, alpha=0.5, cmap=cm.viridis)
|
|||
|
|
pyplot.colorbar()
|
|||
|
|
pyplot.contour(X, Y, p, cmap=cm.viridis)
|
|||
|
|
pyplot.streamplot(X, Y, u, v)
|
|||
|
|
pyplot.xlabel('X')
|
|||
|
|
pyplot.ylabel('Y');
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/0J2aJ/btq9e57PmFB/UqArsnX9hzJ84H4rW5AtB1/img.png" >}}
|
|||
|
|
|
|||
|
|
## Channel Flow with Navier–Stokes
|
|||
|
|
|
|||
|
|
$$ \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho}\frac{\partial p}{\partial x}+\nu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}\right)+F $$$$ \frac{\partial^2 p}{\partial x^2}+\frac{\partial^2 p}{\partial y^2}=-\rho\left(\frac{\partial u}{\partial x}\frac{\partial u}{\partial x}+2\frac{\partial u}{\partial y}\frac{\partial v}{\partial x}+\frac{\partial v}{\partial y}\frac{\partial v}{\partial y}\right) $$$$ \begin{split}p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2(\Delta x^2+\Delta y^2)} \\& -\frac{\rho\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \\& \times \left[\frac{1}{\Delta t} \left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} + \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right) - \frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} - 2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x} - \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]\end{split} $$
|
|||
|
|
|
|||
|
|
```python
|
|||
|
|
import numpy
|
|||
|
|
from matplotlib import pyplot, cm
|
|||
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|||
|
|
%matplotlib inline
|
|||
|
|
|
|||
|
|
def build_up_b(rho, dt, dx, dy, u, v):
|
|||
|
|
b = numpy.zeros_like(u)
|
|||
|
|
b[1:-1, 1:-1] = (rho * (1 / dt * ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx) +
|
|||
|
|
(v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) -
|
|||
|
|
((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx))**2 -
|
|||
|
|
2 * ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dy) *
|
|||
|
|
(v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dx))-
|
|||
|
|
((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy))**2))
|
|||
|
|
|
|||
|
|
# Periodic BC Pressure @ x = 2
|
|||
|
|
b[1:-1, -1] = (rho * (1 / dt * ((u[1:-1, 0] - u[1:-1,-2]) / (2 * dx) +
|
|||
|
|
(v[2:, -1] - v[0:-2, -1]) / (2 * dy)) -
|
|||
|
|
((u[1:-1, 0] - u[1:-1, -2]) / (2 * dx))**2 -
|
|||
|
|
2 * ((u[2:, -1] - u[0:-2, -1]) / (2 * dy) *
|
|||
|
|
(v[1:-1, 0] - v[1:-1, -2]) / (2 * dx)) -
|
|||
|
|
((v[2:, -1] - v[0:-2, -1]) / (2 * dy))**2))
|
|||
|
|
|
|||
|
|
# Periodic BC Pressure @ x = 0
|
|||
|
|
b[1:-1, 0] = (rho * (1 / dt * ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx) +
|
|||
|
|
(v[2:, 0] - v[0:-2, 0]) / (2 * dy)) -
|
|||
|
|
((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx))**2 -
|
|||
|
|
2 * ((u[2:, 0] - u[0:-2, 0]) / (2 * dy) *
|
|||
|
|
(v[1:-1, 1] - v[1:-1, -1]) / (2 * dx))-
|
|||
|
|
((v[2:, 0] - v[0:-2, 0]) / (2 * dy))**2))
|
|||
|
|
|
|||
|
|
return b
|
|||
|
|
|
|||
|
|
def pressure_poisson_periodic(p, dx, dy):
|
|||
|
|
pn = numpy.empty_like(p)
|
|||
|
|
|
|||
|
|
for q in range(nit):
|
|||
|
|
pn = p.copy()
|
|||
|
|
p[1:-1, 1:-1] = (((pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2 +
|
|||
|
|
(pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2) /
|
|||
|
|
(2 * (dx**2 + dy**2)) -
|
|||
|
|
dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, 1:-1])
|
|||
|
|
|
|||
|
|
# Periodic BC Pressure @ x = 2
|
|||
|
|
p[1:-1, -1] = (((pn[1:-1, 0] + pn[1:-1, -2])* dy**2 +
|
|||
|
|
(pn[2:, -1] + pn[0:-2, -1]) * dx**2) /
|
|||
|
|
(2 * (dx**2 + dy**2)) -
|
|||
|
|
dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, -1])
|
|||
|
|
|
|||
|
|
# Periodic BC Pressure @ x = 0
|
|||
|
|
p[1:-1, 0] = (((pn[1:-1, 1] + pn[1:-1, -1])* dy**2 +
|
|||
|
|
(pn[2:, 0] + pn[0:-2, 0]) * dx**2) /
|
|||
|
|
(2 * (dx**2 + dy**2)) -
|
|||
|
|
dx**2 * dy**2 / (2 * (dx**2 + dy**2)) * b[1:-1, 0])
|
|||
|
|
|
|||
|
|
# Wall boundary conditions, pressure
|
|||
|
|
p[-1, :] =p[-2, :] # dp/dy = 0 at y = 2
|
|||
|
|
p[0, :] = p[1, :] # dp/dy = 0 at y = 0
|
|||
|
|
|
|||
|
|
return p
|
|||
|
|
|
|||
|
|
##variable declarations
|
|||
|
|
nx = 41
|
|||
|
|
ny = 41
|
|||
|
|
nt = 10
|
|||
|
|
nit = 50
|
|||
|
|
c = 1
|
|||
|
|
dx = 2 / (nx - 1)
|
|||
|
|
dy = 2 / (ny - 1)
|
|||
|
|
x = numpy.linspace(0, 2, nx)
|
|||
|
|
y = numpy.linspace(0, 2, ny)
|
|||
|
|
X, Y = numpy.meshgrid(x, y)
|
|||
|
|
|
|||
|
|
|
|||
|
|
##physical variables
|
|||
|
|
rho = 1
|
|||
|
|
nu = .1
|
|||
|
|
F = 1
|
|||
|
|
dt = .01
|
|||
|
|
|
|||
|
|
#initial conditions
|
|||
|
|
u = numpy.zeros((ny, nx))
|
|||
|
|
un = numpy.zeros((ny, nx))
|
|||
|
|
|
|||
|
|
v = numpy.zeros((ny, nx))
|
|||
|
|
vn = numpy.zeros((ny, nx))
|
|||
|
|
|
|||
|
|
p = numpy.ones((ny, nx))
|
|||
|
|
pn = numpy.ones((ny, nx))
|
|||
|
|
|
|||
|
|
b = numpy.zeros((ny, nx))
|
|||
|
|
|
|||
|
|
udiff = 1
|
|||
|
|
stepcount = 0
|
|||
|
|
|
|||
|
|
while udiff > .001:
|
|||
|
|
un = u.copy()
|
|||
|
|
vn = v.copy()
|
|||
|
|
|
|||
|
|
b = build_up_b(rho, dt, dx, dy, u, v)
|
|||
|
|
p = pressure_poisson_periodic(p, dx, dy)
|
|||
|
|
|
|||
|
|
u[1:-1, 1:-1] = (un[1:-1, 1:-1] -
|
|||
|
|
un[1:-1, 1:-1] * dt / dx *
|
|||
|
|
(un[1:-1, 1:-1] - un[1:-1, 0:-2]) -
|
|||
|
|
vn[1:-1, 1:-1] * dt / dy *
|
|||
|
|
(un[1:-1, 1:-1] - un[0:-2, 1:-1]) -
|
|||
|
|
dt / (2 * rho * dx) *
|
|||
|
|
(p[1:-1, 2:] - p[1:-1, 0:-2]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])) +
|
|||
|
|
F * dt)
|
|||
|
|
|
|||
|
|
v[1:-1, 1:-1] = (vn[1:-1, 1:-1] -
|
|||
|
|
un[1:-1, 1:-1] * dt / dx *
|
|||
|
|
(vn[1:-1, 1:-1] - vn[1:-1, 0:-2]) -
|
|||
|
|
vn[1:-1, 1:-1] * dt / dy *
|
|||
|
|
(vn[1:-1, 1:-1] - vn[0:-2, 1:-1]) -
|
|||
|
|
dt / (2 * rho * dy) *
|
|||
|
|
(p[2:, 1:-1] - p[0:-2, 1:-1]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])))
|
|||
|
|
|
|||
|
|
# Periodic BC u @ x = 2
|
|||
|
|
u[1:-1, -1] = (un[1:-1, -1] - un[1:-1, -1] * dt / dx *
|
|||
|
|
(un[1:-1, -1] - un[1:-1, -2]) -
|
|||
|
|
vn[1:-1, -1] * dt / dy *
|
|||
|
|
(un[1:-1, -1] - un[0:-2, -1]) -
|
|||
|
|
dt / (2 * rho * dx) *
|
|||
|
|
(p[1:-1, 0] - p[1:-1, -2]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(un[1:-1, 0] - 2 * un[1:-1,-1] + un[1:-1, -2]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(un[2:, -1] - 2 * un[1:-1, -1] + un[0:-2, -1])) + F * dt)
|
|||
|
|
|
|||
|
|
# Periodic BC u @ x = 0
|
|||
|
|
u[1:-1, 0] = (un[1:-1, 0] - un[1:-1, 0] * dt / dx *
|
|||
|
|
(un[1:-1, 0] - un[1:-1, -1]) -
|
|||
|
|
vn[1:-1, 0] * dt / dy *
|
|||
|
|
(un[1:-1, 0] - un[0:-2, 0]) -
|
|||
|
|
dt / (2 * rho * dx) *
|
|||
|
|
(p[1:-1, 1] - p[1:-1, -1]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(un[1:-1, 1] - 2 * un[1:-1, 0] + un[1:-1, -1]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(un[2:, 0] - 2 * un[1:-1, 0] + un[0:-2, 0])) + F * dt)
|
|||
|
|
|
|||
|
|
# Periodic BC v @ x = 2
|
|||
|
|
v[1:-1, -1] = (vn[1:-1, -1] - un[1:-1, -1] * dt / dx *
|
|||
|
|
(vn[1:-1, -1] - vn[1:-1, -2]) -
|
|||
|
|
vn[1:-1, -1] * dt / dy *
|
|||
|
|
(vn[1:-1, -1] - vn[0:-2, -1]) -
|
|||
|
|
dt / (2 * rho * dy) *
|
|||
|
|
(p[2:, -1] - p[0:-2, -1]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(vn[1:-1, 0] - 2 * vn[1:-1, -1] + vn[1:-1, -2]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(vn[2:, -1] - 2 * vn[1:-1, -1] + vn[0:-2, -1])))
|
|||
|
|
|
|||
|
|
# Periodic BC v @ x = 0
|
|||
|
|
v[1:-1, 0] = (vn[1:-1, 0] - un[1:-1, 0] * dt / dx *
|
|||
|
|
(vn[1:-1, 0] - vn[1:-1, -1]) -
|
|||
|
|
vn[1:-1, 0] * dt / dy *
|
|||
|
|
(vn[1:-1, 0] - vn[0:-2, 0]) -
|
|||
|
|
dt / (2 * rho * dy) *
|
|||
|
|
(p[2:, 0] - p[0:-2, 0]) +
|
|||
|
|
nu * (dt / dx**2 *
|
|||
|
|
(vn[1:-1, 1] - 2 * vn[1:-1, 0] + vn[1:-1, -1]) +
|
|||
|
|
dt / dy**2 *
|
|||
|
|
(vn[2:, 0] - 2 * vn[1:-1, 0] + vn[0:-2, 0])))
|
|||
|
|
|
|||
|
|
|
|||
|
|
# Wall BC: u,v = 0 @ y = 0,2
|
|||
|
|
u[0, :] = 0
|
|||
|
|
u[-1, :] = 0
|
|||
|
|
v[0, :] = 0
|
|||
|
|
v[-1, :]=0
|
|||
|
|
|
|||
|
|
udiff = (numpy.sum(u) - numpy.sum(un)) / numpy.sum(u)
|
|||
|
|
stepcount += 1
|
|||
|
|
|
|||
|
|
fig = pyplot.figure(figsize = (11,7), dpi=100)
|
|||
|
|
pyplot.quiver(X[::3, ::3], Y[::3, ::3], u[::3, ::3], v[::3, ::3]);
|
|||
|
|
```
|
|||
|
|
|
|||
|
|
{{< figure src="https://blog.kakaocdn.net/dn/du6hla/btq9fdKZP6o/ifKi67Tsr8khMmReNSHn5K/img.png" >}}
|
|||
|
|
|
|||
|
|
출처> [CFD Python: 12 steps to Navier-Stokes :: Lorena A. Barba Group (lorenabarba.com)](https://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/)
|
|||
|
|
|