Help us understand the problem. What is going on with this article?

# 非定常クエット流の計算

More than 1 year has passed since last update.

## はじめに

$$\frac{\partial u}{\partial t} = \nu \frac{\partial^2 u}{\partial y^2} \tag{1}$$

ここで$\nu$は流体の動粘性係数、$u$は$x$方向の速度成分である。

$${\displaystyle u(y,0)=0,\quad 0<y<h,} \tag{2}$$

$${\displaystyle u(0,t)=0,\quad u(h,t)=U,\quad t>0.} \tag{3}$$

この時、運動方程式(1)は解析的に解くことができる。解析解は以下の様に表される。

{\displaystyle u(y,t)=U{\frac {y}{h}}-{\frac {2U}{\pi }}\sum _{n=1}^{\infty }{\frac {1}{n}}e^{-n^{2}\pi ^{2}\nu t/h^{2}}\sin \left[n\pi \left(1-{\frac {y}{h}}\right)\right]} \tag{4}


# Pythonによるシミュレーション

import numpy as np
import matplotlib.pyplot as plt

''' Analytical solution of Couette flow '''
def u_exact(U,h,nu,y,t):
u = U * (1-y / h)
for n in range(1,100000):
u -= 2 * U / np.pi * (1/n) *  np.e **(-n**2*np.pi**2*nu*t/h**2)*np.sin(n*np.pi*y/h)
return u

''' Parameters '''
h = 1   # gap between plate
U = 1   # velocity of under wall
nu = 0.0005 # viscosity coefficient

nx = 51 # number of mesh

y = h * np.linspace(0,1,nx) # y coordinates
dy = h / (nx-1) # mesh size
dt = 0.02   # timestep
nend = 20000    # number of iteration

''' Initial conditions '''
u = np.zeros(nx)
un = np.zeros(nx)
u1 = np.zeros(nx)

u[0] = U
u1[0] = U
un[0] = U

''' 2nd order Runge-Kutta '''
for nt in range(nend):
u1[1:nx-1] = u[1:nx-1] + 0.5 * dt * nu * (u[0:nx-2] - 2 * u[1:nx-1] + u[2:nx]) / dy**2

#    for j in range(1,nx-1):
#        u1[j] = u[j] + 0.5* dt * nu * (u[j-1] - 2 * u[j] + u[j+1]) / dy**2

un[1:nx-1] = u[1:nx-1] + dt * nu * (u1[0:nx-2] - 2 * u1[1:nx-1] + u1[2:nx]) / dy**2
#    for j in range(1,nx-1):
#        un[j] = u[j] + dt * nu * (u1[j-1] - 2 * u1[j] + u1[j+1]) / dy**2

u = un.copy()

''' Plot '''
if np.mod(nt + 3000, 4000) == 0 and nt < 5000 and nt > 0:
ue = u_exact(U,h,nu,y,dt*nt)

plt.plot(ue,y, color = 'black', label = "Analytical solution")
plt.plot(un,y,"o", color = 'None', markeredgecolor='black', label = "Numerical solution")

elif np.mod(nt + 3000, 4000) == 0 and nt > 0:
ue = u_exact(U,h,nu,y,dt*nt)

plt.plot(ue,y, color = 'black')
plt.plot(un,y,"o", color = 'None', markeredgecolor='black')

plt.xlabel("u")
plt.ylabel("y")
plt.legend()
plt.show()


Why not register and get more from Qiita?
1. We will deliver articles that match you
By following users and tags, you can catch up information on technical fields that you are interested in as a whole
2. you can read useful information later efficiently
By "stocking" the articles you like, you can search right away