涡度函数产生特殊结果

Question

我正在研究盖子驱动的空腔流动问题，使用涡度纳维-斯托克斯方程的离散版本：

出于某种原因，我的流函数似乎正在工作（至少我相信是这样），但我认为我的涡度图有点奇怪。我根据建议改变了一些事情，现在我得到了......有趣的结果。

完整代码：

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)")

现在我的最终绘图（流和涡度）如下所示：

老实说，我不确定涡度图是怎么回事。它看起来仍然没有达到预期的涡度？难道我错了？

Answer 1

您的主要问题是您没有正确地对涡度方程进行时间步进，并且您的 x,y 坐标与传统的 [i,j] 矩阵表示法不太相符。

您会发现很难计算平流项或最终速度场、将正确的墙设置为移动边界或绘制最终图，除非您按照如下方式对齐坐标和索引：

y,j
^
|
|
|
+--------> x, i

为了加快求解速度，我将 N 减少到 20，并将 dt 增加到 0.001。如果你想增加 N，那么这将减少空间步长，并且你将不得不减少 dt 来补偿。您可能还需要降低欧米伽的值； （请注意，我是放松不足，而不是过度放松）。

我简化了您的绘图例程。

import numpy as np
import matplotlib.pyplot as plt

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


#Initialising Arrays:
u  = np.zeros((N,N))
w  = np.zeros((N,N))
Vx = np.zeros((N,N));    Vx[:,-1] = Vwall        # Top wall
Vy = np.zeros((N,N))
p  = np.zeros((N,N))


# Update stream function  ( Laplacian(psi) = -w )
def updateStreamfn( u, w, h, tolerance=1e-6, maxIter=1000, omega=0.8 ):
    maxDiff = tolerance + 1
    iter = 0

    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):
                u[i,j] = (1-omega) * uOld[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] )
      
        maxDiff = np.max( np.abs( u - uOld ) )
        iter += 1
    
        if iter >= maxIter:
            print("Reached max iterations without convergence.")
            break

    return u


# Update stream function  ( dw/dt + Vx.dw/dx + Vy.dw/dy = (1/Re)Laplacian(w)
def updateWallVorticity( w, u, Vwall, h ):
    w[: ,-1] = (2/h**2) * ( u[: ,-1] - u[: ,-2] ) - (2/h) * Vwall
    w[: ,0 ] = (2/h**2) * ( u[: ,0 ] - u[: ,1 ] )
    w[0 ,: ] = (2/h**2) * ( u[0 ,: ] - u[1 ,: ] )
    w[-1,: ] = (2/h**2) * ( u[-1,: ] - u[-2,: ] )


def updateInteriorVorticity( w, u, dt, Re, h, tolerance=1e-6, maxIter=1000, omega=0.8 ):
    wPrev = np.copy(w)                 # value at end of last TIMESTEP
    maxDiff = tolerance + 1
    iter = 0

    while maxDiff > tolerance and iter < maxIter:
       wOld = np.copy(w)               # value at end of last ITERATION
    
       for i in range( 1, w.shape[0] - 1 ):
           for j in range( 1, w.shape[1] - 1 ):
               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] + w[i-1,j] + w[i,j+1] + w[i,j-1] - 4*w[i,j] ) / h**2
               w[i,j] = (1-omega) * w[i,j]   +   omega * ( wPrev[i,j] + dt * ( -advectionx + advectiony + diffusion) )
        
       maxDiff = np.max(np.abs(w - wOld))
       iter += 1

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



def plotStreamfn( u, title="Stream Function (u)", time_step=None ):
    plt.figure( figsize=(8, 6) )
    x = np.linspace( 0, L, u.shape[0] )
    y = np.linspace( 0, L, u.shape[1] )
    Y, X = np.meshgrid( x, y )
    plt.contourf( X, Y, u, levels=20, cmap='viridis' )
    plt.colorbar( label='Stream Function Value' )
    plt.contour( X, Y, u, levels=10, colors='white', linewidths=1.0, linestyles='solid' )

    if time_step is not None:
        plt.title(f"{title} at Time Step {time_step}")
    else:
        plt.title(title)
        Vx[1:-1, 1:-1] =  ( u[ 1 : -1 , 2 :    ] - u[ 1: -1 ,   : -2 ] ) / (2 * h)    #  d(psi)/dy
        Vy[1:-1, 1:-1] = -( u[ 2 :    , 1 : -1 ] - u[  : -2 , 1 : -1 ] ) / (2 * h)    # -d(psi)/dx
        plt.quiver( X, Y, Vx, Vy )
    plt.show()


#Vorticity plot:
def plotVorticity( w, title="Vorticity (w)", time_step=None ):
    plt.figure( figsize=(8, 6) )
    x = np.linspace( 0, L, u.shape[0] )
    y = np.linspace( 0, L, u.shape[1] )
    Y, X = np.meshgrid( x, y )
    contour = plt.contourf( X, Y, w, levels=50, cmap='viridis' )
    plt.colorbar(contour, label='Vorticity Value')
    plt.contour( X, Y, w, levels=10, colors='white', linewidths=1.0, linestyles='solid' )
    
    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 = 25

for n in range(numTimesteps):
    u = updateStreamfn(u, w, h)
    updateWallVorticity(w, u, Vwall, h)
    w = updateInteriorVorticity(w, u, dt, Re, h)
    
    if n % visualInter == 0:
       print(f"Time Step: {n}")
       plotStreamfn(u, title=f"Stream Function (u)", time_step=n)
       plotVorticity(w, title=f"Vorticity (w)", time_step=n)

plotStreamfn( u, title="Final Stream Function (u)" )
plotVorticity( w, title="Final Vorticity (w)" )