顶盖驱动腔流模拟故障排查：首行更新/涡量内部更新异常

我正在模拟顶盖驱动腔流，模拟流程如下：

• 初始化边界条件
• 采用离散化形式结合SOR算法更新流函数（满足收敛条件时停止）
• 更新涡量边界条件
• 更新涡量内部场
• 循环执行上述更新步骤

后续需要获取速度场与压力场，但当前模拟存在问题：此前网格仅首行更新，调试后流函数更新恢复正常，但涡量内部仍未正确更新。以下是我的Python实现代码：

import numpy as np
import matplotlib.pyplot as plt
import copy

def Lid_Driven_Cavity_Flow(N, u, w, omega, dt, V_lid, Re):
    """ 
    Function to run simulation of lid driven cavity flow
    """
    h=1/N

    def BC_streamfunction():
        """
        Impose boundary conditions on streamfunction
        """
        u[0,:] = V_lid   #lid
        u[1:,N-1] = 0    #right wall
        u[1:,0] = 0      #left wall
        u[N-1,1:N] = 0   #bottom
    
    def BC_vorticity():
        """
        Impose boundary conditions on vorticity
        """
        w[0,:] = (2/h**2) * (u[0,:] - u[1,:]) + (2/h) * V_lid    #lid
        w[1:,N-1] = (2/h**2) * (u[1:N,N-1] - u[1:N,N-2])         #right wall
        w[1:,0] = (2/h**2) * (u[1:N,0] - u[1:N,1])               #left wall
        w[N-1,1:N] = (2/h**2) * (u[N-1,1:N] - u[N-2,1:N])        #bottom
    
    def Update_Streamfunction():   
        """
        This function updates the streamfunction on the interior
        """     
        while True:
            u_prev = copy.deepcopy(u)
            for i in range(1,N-1):     #columns
                for j in range(1,N-1): #rows
                    #print("coordinates are", j, i)
                    r1 = (1/4 * (w[j,i]*h**2 + u[j,i+1] + u[j,i-1] +           \
                          u[j+1,i] + u[j-1,i])) - u_prev[j,i]
                    u[j,i] += omega * r1
                    #print(r1)
                    #print(u)
        
            if np.max(((1/N**2)*np.abs(u_prev - u))) < 10e-6: #Convergence check
                break

    def Update_Vorticity_Interior():
        """
        This function updates the vorticity field on the interior
        """
        while True: 
            w_prev = copy.deepcopy(w)
            for i in range(1,N-1): #N-2 is on the wall, but range doesn't 
                                   #include N-1
                for j in range(1,N-1):
                    a = (-1)/(4*h**2) * (u[j+1,i] - u[j-1,i]) * (w[j,i+1] -    \
                         w[j,i-1])
                    b = (1/4*h**2) * (u[j,i+1] - u[j,i-1]) * (w[j+1,i] -       \
                         w[j-1,i])
                    c = (1/Re*h**2) * (u[j,i+1] - 4*u[j,i] + u[j,i-1] +        \
                         u[j+1,i] + u[j-1,i])

                    r2 = omega * dt * (a + b + c) - w_prev[j,i]

                    w[j,i] += r2

            if np.max(((1/N**2)*np.abs(w_prev - w))) < 10e-6: #Convergence check
                break
    
    """
    Below, there is a while loop that iterates over a number of timesteps.
    First the boundary conditions on the streamfunction and vorticity are 
    intialized. 
    Then for each timestep the streamfunction is updated, then the boundary
    conditions on the interior are updated and lastly the vorticity interior 
    is updated. Then the loop starts over again. 
    """
    t = 0
    tmax = 10
    steps = int(tmax/dt)
    BC_streamfunction()
    #print("After running BC_streamfunction, the matrices are \n", u,"\n", w)
    BC_vorticity()
    #print("After running BC_vorticity, the matrices are \n", u,"\n", w)
    while t < steps:
       t += 1
       Update_Streamfunction()
       #print("After running Update_Streamfunction, the matrices are \n", u, "\n", w)
       BC_vorticity()
       #print("After running BC_vorticity, the matrices are \n", u,"\n", w)
       Update_Vorticity_Interior()
       #print("After running Update_Vorticity_Interior, the matrices are \n", u,"\n", w)

    def Velocity(): 
        """
        Calculate velocity in x- and y-direction, and use the length of the 
        tangent vector as value for the velocity on a grid point.
        """
        v = np.zeros([N,N])
        for i in range(0,N-1):
            for j in range(0,N-1):
                v_y = (-1) * (u[j,i+1] + u[j,i-1]) / (2*h)
                v_x = (u[j+1,i] - u[j-1,i]) / (2*h)
                v[j,i] = np.sqrt(v_y**2+v_x**2)
        return v
    
    def ContourPlot(matrix, title, xtitle, ytitle, number_levels):
        """
        This funciton produces a contourplot. 
        """
        x = range(0,N)
        y = range(0,N)
        X, Y = np.meshgrid(x,y)
        Z = matrix[X,Y]

        fig, ax = plt.subplots()
        CS = ax.contour(X, Y, Z, levels=number_levels)
        ax.set_title(title)
        ax.set_xlabel(xtitle)
        ax.set_ylabel(ytitle)

    ContourPlot(u, 'Streamfunction', 'j-direction', 'i-direction', 100)
    ContourPlot(w, 'Vorticity', 'j-direction', 'i-direction', 100)
    ContourPlot(Velocity(), 'Velocityfield', 'j-direction', 'i-direction', 100)

    def Pressure():
        """
        Here, the pressure will be calculated.
        Still working on this function.
        """
        return 0

Lid_Driven_Cavity_Flow(10, np.zeros([10,10]), np.zeros([10,10]), 0.25, 0.01, 1, 100)

问题定位与修复方案 #

1. 涡量更新公式核心错误 #

Update_Vorticity_Interior中的离散逻辑完全偏离了涡量输运方程的正确形式，且混淆了SOR迭代与时间步更新的逻辑：

• 原代码错误地给涡量更新加了while循环迭代，显式时间推进格式不需要内循环迭代
• 对流项、扩散项的离散公式推导错误，更新增量r2的计算逻辑混乱

正确的涡量输运方程离散实现如下：

def Update_Vorticity_Interior():
    """
    显式时间推进更新内部涡量场
    """
    w_prev = copy.deepcopy(w)
    # 从流函数计算速度场u、v
    u_vel = np.zeros_like(u)
    v_vel = np.zeros_like(u)
    for j in range(1, N-1):
        for i in range(1, N-1):
            u_vel[j,i] = (u[j+1,i] - u[j-1,i])/(2*h)  # u = ∂ψ/∂y
            v_vel[j,i] = -(u[j,i+1] - u[j,i-1])/(2*h) # v = -∂ψ/∂x
    
    # 涡量输运方程离散计算
    for j in range(1, N-1):
        for i in range(1, N-1):
            # 对流项：u*∂ω/∂x + v*∂ω/∂y
            conv_x = u_vel[j,i] * (w_prev[j,i+1] - w_prev[j,i-1])/(2*h)
            conv_y = v_vel[j,i] * (w_prev[j+1,i] - w_prev[j-1,i])/(2*h)
            # 扩散项：(1/Re)*∇²ω
            diff = (1/Re) * (w_prev[j,i+1] + w_prev[j,i-1] + w_prev[j+1,i] + w_prev[j-1,i] - 4*w_prev[j,i])/(h**2)
            # 显式时间推进更新
            w[j,i] = w_prev[j,i] + dt*(diff - conv_x - conv_y)

2. 顶盖涡量边界条件符号错误 #

顶盖的涡量边界条件公式应为：
$$\omega_{0,i} = -\frac{2}{h^2}(\psi_{0,i} - \psi_{1,i}) - \frac{2}{h}U_{lid}$$
原代码中+ (2/h)*V_lid符号错误，需改为- (2/h)*V_lid，否则顶盖涡量符号完全反向，影响整个流场。

3. 流函数边界条件设置错误 #

流函数的物理边界条件是：固壁上速度为0，因此流函数的法向导数为0，通常取所有固壁上的流函数为常数0，不需要将顶盖流函数设为V_lid（速度由流函数的导数推导）。修正后的流函数边界条件：

def BC_streamfunction():
    """
    设置流函数边界条件
    """
    u[:,0] = 0    # 左壁
    u[:,N-1] = 0  # 右壁
    u[N-1,:] = 0  # 底壁
    u[0,:] = 0    # 顶盖（流函数取常数0）

4. 速度场计算越界问题 #

Velocity函数中，当i=0或j=0时，u[j,i-1]会触发数组越界，需限制遍历范围为内部点，边界点直接使用已知边界条件：

def Velocity(): 
    """
    计算速度场的大小
    """
    v_mag = np.zeros([N,N])
    v_mag[0,:] = V_lid  # 顶盖速度直接设为边界条件
    # 仅计算内部点速度
    for j in range(1, N-1):
        for i in range(1, N-1):
            u_vel = (u[j+1,i] - u[j-1,i])/(2*h)
            v_vel = -(u[j,i+1] - u[j,i-1])/(2*h)
            v_mag[j,i] = np.sqrt(u_vel**2 + v_vel**2)
    return v_mag

修正后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt
import copy

def Lid_Driven_Cavity_Flow(N, u, w, omega, dt, V_lid, Re):
    """ 
    顶盖驱动腔流模拟函数
    """
    h = 1/N

    def BC_streamfunction():
        """
        设置流函数边界条件
        """
        u[:,0] = 0    # 左壁
        u[:,N-1] = 0  # 右壁
        u[N-1,:] = 0  # 底壁
        u[0,:] = 0    # 顶盖
    
    def BC_vorticity():
        """
        设置涡量边界条件
        """
        # 顶盖边界条件
        w[0,:] = -2*(u[0,:] - u[1,:])/(h**2) - 2*V_lid/h
        # 右壁边界条件
        w[1:N-1, N-1] = 2*(u[1:N-1, N-1] - u[1:N-1, N-2])/(h**2)
        # 左壁边界条件
        w[1:N-1, 0] = 2*(u[1:N-1, 0] - u[1:N-1, 1])/(h**2)
        # 底壁边界条件
        w[N-1, 1:N-1] = -2*(u[N-1, 1:N-1] - u[N-2, 1:N-1])/(h**2)
    
    def Update_Streamfunction():   
        """
        SOR算法更新内部流函数
        """     
        max_iter = 10000
        tol = 1e-6
        for _ in range(max_iter):
            u_prev = copy.deepcopy(u)
            for i in range(1, N-1):     
                for j in range(1, N-1): 
                    # 流函数与涡量关系：∇²ψ = -ω，离散后更新
                    u[j,i] = omega * ((u[j+1,i] + u[j-1,i] + u[j,i+1] + u[j,i-1] + w[j,i]*h**2)/4) + (1-omega)*u[j,i]
            if np.max(np.abs(u_prev - u)) < tol:
                break

    def Update_Vorticity_Interior():
        """
        显式时间推进更新内部涡量场
        """
        w_prev = copy.deepcopy(w)
        # 计算速度场
        u_vel = np.zeros_like(u)
        v_vel = np.zeros_like(u)
        for j in range(1, N-1):
            for i in range(1, N-1):
                u_vel[j,i] = (u[j+1,i] - u[j-1,i])/(2*h)
                v_vel[j,i] = -(u[j,i+1] - u[j,i-1])/(2*h)
        
        # 更新涡量
        for j in range(1, N-1):
            for i in range(1, N-1):
                conv_x = u_vel[j,i] * (w_prev[j,i+1] - w_prev[j,i-1])/(2*h)
                conv_y = v_vel[j,i] * (w_prev[j+1,i] - w_prev[j-1,i])/(2*h)
                diff = (1/Re) * (w_prev[j,i+1] + w_prev[j,i-1] + w_prev[j+1,i] + w_prev[j-1,i] - 4*w_prev[j,i])/(h**2)
                w[j,i] = w_prev[j,i] + dt*(diff - conv_x - conv_y)
    
    # 初始化边界条件
    BC_streamfunction()
    BC_vorticity()
    
    # 时间步循环
    tmax = 10
    steps = int(tmax/dt)
    for t in range(steps):
        Update_Streamfunction()
        BC_vorticity()
        Update_Vorticity_Interior()
        if t % 100 == 0:
            print(f"完成时间步 {t}/{steps}")

    def Velocity(): 
        """
        计算速度场的大小
        """
        v_mag = np.zeros([N,N])
        v_mag[0,:] = V_lid
        for j in range(1, N-1):
            for i in range(1, N-1):
                u_vel = (u[j+1,i] - u[j-1,i])/(2*h)
                v_vel = -(u[j,i+1] - u[j,i-1])/(2*h)
                v_mag[j,i] = np.sqrt(u_vel**2 + v_vel**2)
        return v_mag
    
    def ContourPlot(matrix, title, xtitle, ytitle, number_levels):
        """
        绘制等高线图
        """
        x = np.linspace(0, 1, N)
        y = np.linspace(0, 1, N)
        X, Y = np.meshgrid(x, y)
        # 反转y轴，匹配物理空间（顶盖在上）
        Y = np.flip(Y)
        Z = np.flip(matrix, axis=0)

        fig, ax = plt.subplots(figsize=(8,6))
        CS = ax.contourf(X, Y, Z, levels=number_levels, cmap='viridis')
        fig.colorbar(CS)
        ax.set_title(title)
        ax.set_xlabel(xt