anonymous No title
No License Python
2021年06月01日
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+1,N+1))
ou=np.zeros((N+1,N+1))
us=np.zeros((N+1,N+1))
nv=np.zeros((N+1,N+1))
ov=np.zeros((N+1,N+1))
vs=np.zeros((N+1,N+1))
pre=np.zeros((N+1,N+1))

#初期条件セッツ
for i in range(N):
    ou[i][N-1]=math.sin(pi*i*dx)*math.sin(pi*i*dx)
    nu[i][N-1]=math.sin(pi*i*dx)*math.sin(pi*i*dx)
    us[i][N-1]=math.sin(pi*i*dx)*math.sin(pi*i*dx)
    

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


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


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

x=np.arange(0,1.0+dx,dx)
y=np.arange(0,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, nv)
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+1,N+1))
ou=np.zeros((N+1,N+1))
us=np.zeros((N+1,N+1))
nv=np.zeros((N+1,N+1))
ov=np.zeros((N+1,N+1))
vs=np.zeros((N+1,N+1))
pre=np.zeros((N+1,N+1))

#初期条件セッツ
for i in range(N):
    ou[i][N-1]=math.sin(pi*i*dx)*math.sin(pi*i*dx)
    nu[i][N-1]=math.sin(pi*i*dx)*math.sin(pi*i*dx)
    us[i][N-1]=math.sin(pi*i*dx)*math.sin(pi*i*dx)
    

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


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


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

x=np.arange(0,1.0+dx,dx)
y=np.arange(0,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, nv)
plt.show()
No one still commented. Please first comment.
Output