顶盖驱动流涡量计算结果异常，请求技术分析

我正在开展顶盖驱动腔流问题的数值模拟工作，采用离散形式的涡量Navier-Stokes方程进行计算：

目前流函数的计算结果看似正常，但涡量图的表现不符合预期。根据建议调整代码后，仍得到异常的涡量结果。

完整代码 #

import numpy as np
import matplotlib.pyplot as plt

#Initialising Parameters:
N = 50
L = 1.0
Re = 5
dt = 0.0001
Vwall = 1.0
h = L/(N-1)


#Initialising Arrays:
u = np.zeros((N, N))
w = np.zeros((N,N))
Vx = np.zeros((N, N))
Vy = np.zeros((N, N))
P = np.zeros((N, N))

#Setting Boundary Conditions (Velocity):
Vx[:, 0] = 0  #Left wall
Vx[:, -1] = 0  #Right wall
Vx[0, :] = 0  #Bottom wall
Vx[-1, :] = 1  #Top wall

Vy[:, 0] = 0  #Left wall
Vy[:, -1] = 0  #Right wall
Vy[0, :] = 0  #Bottom wall
Vy[-1, :] = 0  #Top wall


#Update stream function using 1a (Using SOR):
def updateStreamfn(u, w, h, tolerance=1e-6, maxIter=1000, omega=1.3):
  #SOR for solving Poisson Equation (iterating using while & for loops):
  maxDiff = tolerance + 1
  iter = 0

  #New while loop replacing for loop:
  while maxDiff > tolerance and iter < maxIter:
    uOld = np.copy(u)

    for i in range(1, u.shape[0] - 1):
      for j in range(1, u.shape[1] - 1):

        #Update Equation (using discretised 1a):
        u[i, j] = (1-omega) * u[i, j] + (omega * 0.25) *(
            u[i+1, j] + u[i-1, j] + u[i, j+1] + u[i, j-1] + (h**2 * w[i, j])
            )
      
      #Break at convergence after maximum difference value reached:
    maxDiff = np.max(np.abs(u - uOld))
    iter += 1
    
    if iter >= maxIter:
      print("Reached max iterations without convergence.")
      break

  return u


#Update vorticity on wall with eqn 10 for all boundaries:
def updateWallVorticity(w, u, Vwall, h):
  N = w.shape[0]
  
  #Top wall (y=max at moving lid)
  for i in range(N):
    w[i, 0] = -((2 / h**2) * (u[i, 1] - u[i, 0]) + (2 * Vwall / h))

  #Bottom wall (y=0)
  for i in range(N):
    w[i, N-1] = -((2 / h**2) * (u[i, N-2] - u[i, N-1]))

  #Left wall (x=0)
  for j in range(N):
    w[0, j] = -((2 / h**2) * (u[1, j] - u[0, j]))

  #Right wall (x=max)
  for j in range(N):
    w[N-1, j] = -((2 / h**2) * (u[N-2, j] - u[N-1, j]))


#Update vorticity in interior with 1b:
def updateInteriorVorticity(w, u, dt, Re, h, tolerance=1e-6, maxIter=1000, omega=1.3):
  wNew = np.copy(w)
  maxDiff = tolerance + 1
  iter_count = 0

  while maxDiff > tolerance and iter_count < maxIter:
    wOld = np.copy(wNew)
    
    #For loop to update through iterations:
    for i in range(1, w.shape[0] - 1):
      for j in range(1, w.shape[1] - 1):
        
        #Advection and diffusion terms:
        advectionx = ((u[i, j+1] - u[i, j-1]) / (2 * h) * ((w[i+1, j] - w[i-1, j]) / (2 * h)))
        advectiony = ((u[i+1, j] - u[i-1, j]) / (2 * h)) * ((w[i, j+1] - w[i, j-1]) / (2 * h))
        
        diffusion = (1 / Re) * ((w[i+1, j] - (2 * w[i, j]) + w[i-1, j]) / h**2 + (
            w[i, j+1] - (2 * w[i, j]) + w[i, j-1]) / h**2)
        
        #Update formula:
        wNew[i, j] = (1 - omega) * w[i, j] + omega * (w[i, j] + dt * (-advectionx + advectiony + diffusion))
        
    #Convergence check
    maxDiff = np.max(np.abs(wNew - wOld))
    iter_count += 1

    if iter_count >= maxIter:
      print("Vorticity update reached max iterations without full convergence.")
      break
  
    return wNew


#Task 4 plots:
def plotStreamfnFin(u, Vx, Vy, title="Stream Function (u)", time_step=None):
    plt.figure(figsize=(8, 6))
    plt.contourf(u, levels=50, cmap='viridis')
    plt.colorbar(label='Stream Function Value')

    #Plotting time steps:
    if time_step is not None:
        plt.title(f"{title} at Time Step {time_step}")
    else:
        plt.title(title)
    plt.xlabel("Grid X")
    plt.ylabel("Grid Y")
    
    #Streamlines to indicate flow direction:
    Y, X = np.mgrid[0:u.shape[0], 0:u.shape[1]]
    plt.streamplot(X, Y, Vx, Vy, color='white', linewidth=0.5)
    plt.show()

def plotStreamfnInt(u, title="Stream Function (u)", time_step=None):
    plt.figure(figsize=(8, 6))
    plt.contourf(u, levels=50, cmap='viridis')
    plt.colorbar(label='Stream Function Value')

    #Plotting time steps:
    if time_step is not None:
        plt.title(f"{title} at Time Step {time_step}")
    else:
        plt.title(title)
    plt.xlabel("Grid X")
    plt.ylabel("Grid Y")


#Vorticity plot:
def plotVorticity(w, title="Vorticity (w)", time_step=None):
    plt.figure(figsize=(8, 6))
    contour = plt.contourf(w, levels=50, cmap='viridis')
    plt.colorbar(contour, label='Vorticity Value')
    
    #Overlay contour lines:
    plt.contour(w, levels=10, colors='white', linewidths=0.5, linestyles='solid')
    
    #Title and labels:
    if time_step is not None:
        plt.title(f"{title} at Time Step {time_step}")
    else:
        plt.title(title)
    plt.xlabel("Grid X")
    plt.ylabel("Grid Y")
    plt.show()


#Main time-stepping loop:
numTimesteps = 100
visualInter = 10

for n in range(numTimesteps):
    #Update stream function:
    u = updateStreamfn(u, w, h)
    
    #Update vorticity on boundaries:
    updateWallVorticity(w, u, Vwall, h)

    #Update vorticity in interior:
    w = updateInteriorVorticity(w, u, dt, Re, h)
    
    if n % visualInter == 0:
      print(f"Time Step: {n}")

      #Plotting stream function with streamlines:
      plotStreamfnInt(u, title=f"Stream Function (u)", time_step=n)

      #Plotting the vorticity:
      plotVorticity(w, title=f"Vorticity (w)", time_step=n)
    

#Calculating Vx and Vy:
      Vx[1:-1, 1:-1] = (u[1:-1, 2:] - u[1:-1, :-2]) / (2 * h)
      Vy[1:-1, 1:-1] = -(u[2:, 1:-1] - u[:-2, 1:-1]) / (2 * h)  
  
#Plotting the results at the final timestep
plotStreamfnFin(u, Vx, Vy, title="Final Stream Function (u)")
plotVorticity(w, title="Final Vorticity (w)")

计算结果 #

流函数图（含流线）：

涡量图：

我无法判断涡量图异常的原因，恳请协助排查问题并给出解决建议。

内容的提问来源于stack exchange，提问作者Lucifer Addams