こだっく No title
No License Python
2021年06月03日
Copy Clone
import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
pi=math.pi
Re=float(10.0)
N=20
M=50
r=1.72
T=int(1000)
dx=1.0/N
dy=1.0/N
dt=0.001
nu=np.zeros((N+2,N+2))
ou=np.zeros((N+2,N+2))
us=np.zeros((N+2,N+2))
nv=np.zeros((N+2,N+2))
ov=np.zeros((N+2,N+2))
vs=np.zeros((N+2,N+2))
pre=np.zeros((N+2,N+2))


for t in range(T):
    
    for i in range(N+2):
        ou[i][0]=-ou[i][1]
        ou[i][N+1]=2*math.sin(i*dx*pi)*math.sin(i*dx*pi)-ou[i][N-1]
        ov[0][i]=-ov[1][i]
        ov[N+1][i]=-ov[N][i]

    #u^*,v^*の計算
    for i in range(0,N):
        for j in range(1,N):
            vs[i+1][j]=ov[i+1][j]+dt*(-0.25*(ou[i][j+1]+ou[i][j]+ou[i+1][j+1]+ou[i+1][j])*(ov[i+2][j]-ov[i][j])/(2*dx)-ov[i][j]*(ov[i+1][j+1]-ov[i+1][j-1])/(2*dx))+(dt/(Re*dx*dx))*(-4*ov[i+1][j]+ov[i+2][j]+ov[i][j]+ov[i+1][j+1]+ov[i+1][j-1])
    
    for i in range(1,N):
        for j in range(0,N):
            us[i][j+1]=ou[i][j+1]+dt*(-ou[i][j+1]*(ou[i+1][j+1]-ou[i-1][j+1])/(2*dx)-0.25*(ov[i+1][j]+ov[i][j]+ov[i+1][j+1]+ov[i][j+1])*(ou[i][j+2]-ou[i][j])/(2*dx))+(dt/(Re*dx*dx))*(-4*ou[i][j+1]+ou[i+1][j+1]+ou[i-1][j+1]+ou[i][j+2]+ou[i][j])
    

 #圧力の境界条件   
    for i in range(N):
        pre[0][i+1]=pre[1][i+1]-(dx/Re)*(-ou[3][i+1]+4*ou[2][i+1]-5*ou[1][i+1]+2*ou[0][i+1])
        pre[N+1][i+1]=pre[N][i+1]-(dx/Re)*(-ou[N-2][i+1]+4*ou[N-1][i+1]-5*ou[N][i+1]+2*ou[N+1][i+1])
        pre[i+1][0]=pre[i+1][1]-(dx/Re)*(-ov[i+1][3]+4*ov[i+1][2]-5*ov[i+1][1]+2*ov[i+1][0])
        pre[i+1][N+1]=pre[i+1][N]-(dx/Re)*(-ov[i+1][N-2]+4*ov[i+1][N-1]-5*ov[i+1][N]+2*ov[i+1][N+1])
    pre[0][0]=pre[1][0]
    pre[0][N+1]=pre[0][N]
    pre[N+1][0]=pre[N][0]
    pre[N+1][N+1]=pre[N][N+1]
               #緩和法
    for k in range(M):
        for i in range(0,N):
            for j in range(0,N):
                pre[i+1][j+1]=(-Re*dx/4*dt)*(us[i+1][j+1]-us[i][j+1]+vs[i+1][j+1]-vs[i+1][j])+0.25*(pre[i+2][j+1]+pre[i][j+1]+pre[i+1][j+2]+pre[i+1][j])
            #n+1における速度の計算
    for i in range(1,N):
        for j in range(0,N):
            nv[i+1][j]=vs[i+1][j]-dt*(pre[i+1][j+1]-pre[i+1][j])/(dx)

    for i in range(0,N):
        for j in range(1,N):
            nv[i+1][j]=vs[i+1][j]-dt*(pre[i+1][j+1]-pre[i+1][j])/(dx)

    
        #古い速度と新しい速度を更新
    for i in range(0,N+2):
        for j in range(0,N+2):
             ou[i][j]=nu[i][j]
             ov[i][j]=nv[i][j]
    #print(pre[15][15])
    #print(pre[14][14])


x=np.arange(0,dx+1.0+dx,dx)
y=np.arange(0,dx+1.0+dx,dx)
X,Y=np.meshgrid(x,y,indexing='ij')

fig=plt.figure()
ax=Axes3D(fig)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("v")


#print(psi)
ax.plot_wireframe(X, Y, pre)
plt.show()
import numpy as np
import math
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
pi=math.pi
Re=float(10.0)
N=20
M=50
r=1.72
T=int(1000)
dx=1.0/N
dy=1.0/N
dt=0.001
nu=np.zeros((N+2,N+2))
ou=np.zeros((N+2,N+2))
us=np.zeros((N+2,N+2))
nv=np.zeros((N+2,N+2))
ov=np.zeros((N+2,N+2))
vs=np.zeros((N+2,N+2))
pre=np.zeros((N+2,N+2))


for t in range(T):
    
    for i in range(N+2):
        ou[i][0]=-ou[i][1]
        ou[i][N+1]=2*math.sin(i*dx*pi)*math.sin(i*dx*pi)-ou[i][N-1]
        ov[0][i]=-ov[1][i]
        ov[N+1][i]=-ov[N][i]

    #u^*,v^*の計算
    for i in range(0,N):
        for j in range(1,N):
            vs[i+1][j]=ov[i+1][j]+dt*(-0.25*(ou[i][j+1]+ou[i][j]+ou[i+1][j+1]+ou[i+1][j])*(ov[i+2][j]-ov[i][j])/(2*dx)-ov[i][j]*(ov[i+1][j+1]-ov[i+1][j-1])/(2*dx))+(dt/(Re*dx*dx))*(-4*ov[i+1][j]+ov[i+2][j]+ov[i][j]+ov[i+1][j+1]+ov[i+1][j-1])
    
    for i in range(1,N):
        for j in range(0,N):
            us[i][j+1]=ou[i][j+1]+dt*(-ou[i][j+1]*(ou[i+1][j+1]-ou[i-1][j+1])/(2*dx)-0.25*(ov[i+1][j]+ov[i][j]+ov[i+1][j+1]+ov[i][j+1])*(ou[i][j+2]-ou[i][j])/(2*dx))+(dt/(Re*dx*dx))*(-4*ou[i][j+1]+ou[i+1][j+1]+ou[i-1][j+1]+ou[i][j+2]+ou[i][j])
    

 #圧力の境界条件   
    for i in range(N):
        pre[0][i+1]=pre[1][i+1]-(dx/Re)*(-ou[3][i+1]+4*ou[2][i+1]-5*ou[1][i+1]+2*ou[0][i+1])
        pre[N+1][i+1]=pre[N][i+1]-(dx/Re)*(-ou[N-2][i+1]+4*ou[N-1][i+1]-5*ou[N][i+1]+2*ou[N+1][i+1])
        pre[i+1][0]=pre[i+1][1]-(dx/Re)*(-ov[i+1][3]+4*ov[i+1][2]-5*ov[i+1][1]+2*ov[i+1][0])
        pre[i+1][N+1]=pre[i+1][N]-(dx/Re)*(-ov[i+1][N-2]+4*ov[i+1][N-1]-5*ov[i+1][N]+2*ov[i+1][N+1])
    pre[0][0]=pre[1][0]
    pre[0][N+1]=pre[0][N]
    pre[N+1][0]=pre[N][0]
    pre[N+1][N+1]=pre[N][N+1]
               #緩和法
    for k in range(M):
        for i in range(0,N):
            for j in range(0,N):
                pre[i+1][j+1]=(-Re*dx/4*dt)*(us[i+1][j+1]-us[i][j+1]+vs[i+1][j+1]-vs[i+1][j])+0.25*(pre[i+2][j+1]+pre[i][j+1]+pre[i+1][j+2]+pre[i+1][j])
            #n+1における速度の計算
    for i in range(1,N):
        for j in range(0,N):
            nv[i+1][j]=vs[i+1][j]-dt*(pre[i+1][j+1]-pre[i+1][j])/(dx)

    for i in range(0,N):
        for j in range(1,N):
            nv[i+1][j]=vs[i+1][j]-dt*(pre[i+1][j+1]-pre[i+1][j])/(dx)

    
        #古い速度と新しい速度を更新
    for i in range(0,N+2):
        for j in range(0,N+2):
             ou[i][j]=nu[i][j]
             ov[i][j]=nv[i][j]
    #print(pre[15][15])
    #print(pre[14][14])


x=np.arange(0,dx+1.0+dx,dx)
y=np.arange(0,dx+1.0+dx,dx)
X,Y=np.meshgrid(x,y,indexing='ij')

fig=plt.figure()
ax=Axes3D(fig)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("v")


#print(psi)
ax.plot_wireframe(X, Y, pre)
plt.show()
No one still commented. Please first comment.
Output