计算科学 - 泰勒格林涡旋问题的二维泊松求解器 - 吾爱随笔录

我正在尝试使用 RK3 编写 2D Navier Stokes 求解器以在 python 中进行时间推进。为了调试，为了简单起见，我已将 RK3 转换为 Euler 步骤。检查我的预测速度的差异，我发现差异接近于零是可以接受的。然后我去我的压力求解器，解决方案不正确。它甚至不遵循与精确解决方案相同的模式。谁能看到我的错误在哪里？

我的代码不是最有效的，但我并不担心目前只是试图获得物理上正确的答案。u、v 和 p 的解析解由

u = - e^{- 2 t} \cos x \sin y

$u = -e^{-2t}\cos x \sin y$

v = - e^{- 2 t} \sin x \cos y

$v = -e^{-2t}\sin x \cos y$

p = - \frac{e^{- 4 t}}{4} (\cos 2 x + \sin 2 y)

$p = -\frac{e^{-4t}}{4} \left(\cos 2x + \sin 2y \right)$

def Poisson_Solver_Neumann(u, v, Nx, Ny, dt, dx, dy, T, xp, yp):
"""Solves the 2D Poisson equation implicitly on a staggered grid
using Neumann Boundary Conditions

Params:
------
u, v           2D array of float, x and y velocities
Nx, Ny         float, Number of segments in x and y
dt             float, time step size
T              float, current time
X, Y           2D array of float, meshgrid

Returns:
-------
ANeum        2D array of float, A matrix with Neumann conditions
f_RHSn       1D array of float, f(x,y) for Neumann conditions
"""
#Building A
ANeum = np.zeros((Nx*Ny,Nx*Ny),dtype=float)

a = Nx * Ny
b = Nx * Ny
c_int = -4 * (Nx * Ny)
c_edge = -3 * (Nx * Ny)
c_corner = -2 * (Nx * Ny)
d = Nx * Ny
e = Nx * Ny

#Set corner points
ANeum[0,0] = c_corner
ANeum[-1,-1] = c_corner
ANeum[Nx*Ny-Ny,Nx*Ny-Ny] = c_corner
ANeum[Ny-1,Ny-1] = c_corner

#Set edges in first block
for j in range(1,Ny-1):
    ANeum[j,j] = c_edge
    j +=j

#Set edges in last block
for j in range((Nx*Ny)-Ny,(Nx*Ny)-2):
    ANeum[j+1,j+1] = c_edge
    j +=j

#Set edges along main diagonal except for first block
for j in range(Nx+1,Ny*Nx):
    if j % Nx ==0:
        ANeum[j-1,j-1] = c_edge
    j +=j

#Set edges on main diagonal except for last block
for j in range(Nx,(Ny*Nx)-Nx):
    if j % Nx ==0:
        ANeum[j,j] = c_edge
    j +=j

#Second diagonal above and below diagonal
for j in range(Ny,Ny*Nx):
    ANeum[j,j-Ny] = a
    ANeum[j-Nx,j] = e
    j +=j

#first diagonal below main diagonal   
for j in range(1,Ny*Nx):
    if j % Ny ==0:
        ANeum[j,j-1] = 0
    else:
        ANeum[j,j-1] = b
    j +=j

#first diagonal above main diagonal
for j in range(0,Ny*Nx):
    if j % Nx ==0:
        ANeum[j-1,j] = 0
    else:
        ANeum[j-1,j] = d
    j +=j

#Main Diagonal
for j in range(0,Ny*Nx):
    if ANeum[j,j] ==0:
        ANeum[j,j] = c_int

#Setting random point to Dirichlet
ANeum[0,:] = 0
ANeum[0,0] = 1

return ANeum

我的欧拉步骤在这里：

for t in range(0,nt):
un = np.empty_like(u)
vn = np.empty_like(v)
#copy info from last loop, needed for current loop into variable_n
un = u.copy()
vn = v.copy()
#Pn = p.copy()
T = T_0 + t*dt

####STEP 1
### t -----> t + (1/3)*dt
un, vn = BCs(un, vn)
G1 = Build_G1(un, vn, dx, dy, nu)
G2 = Build_G2(un, vn, dx, dy, nu)
u1 = un + (dt*G1) #put 1/3 back in
v1 = vn + (dt*G2)
u1, v1 = BCs(u1, v1)
u_div = (u1[1:-1,2:-1] - u1[1:-1,1:-2])/dx
v_div = (v1[2:-1,1:-1] - v1[1:-2,1:-1])/dy
#Solve for Pressure Field
b1 = np.zeros((Ny+2,Nx+2), dtype=float)
b1[1:-1,1:-1] = (u1[1:-1,2:-1] - u1[1:-1,1:-2])/dx +\
    (v1[2:-1,1:-1] - v1[1:-2,1:-1])/dy
b1 = np.reshape(b1[1:-1,1:-1], Ny*Nx)
A1 = Poisson_Solver_Neumann(u1, v1, Nx, Ny, dt, dx, dy, T, xp, yp)
b1 = b1[:]*(1.0/dt)
temp1 = np.linalg.solve(A1,b1)

p_star = np.reshape(temp1, (Ny,Nx))


F1_Pstar, F2_Pstar = Pressure_Terms(p_star, dx, dy)

#calc predictor velocities
u_star = u1.copy()
v_star = v1.copy()

u_star[1:-1,2:-2] = u1[1:-1,2:-2] - (dt*F1_Pstar)
v_star[2:-2,1:-1] = v1[2:-2,1:-1] - (dt*F2_Pstar)

u = u_star.copy()       
v = v_star.copy()

更新：我的代码现在正在第一步工作，但随后爆炸了。

编辑：一些代码更改