CFD圆柱绕流压力场异常解决：交错网格与泊松方程详解

2025-03-17 07:53:00

CFD：基于交错网格和 Navier-Stokes 方程模拟圆盘绕流

问题是这样的：我尝试用 Python 编写一段程序，模拟在一个交错网格上的圆柱绕流。这段代码在模拟顶盖驱动流和流道流动 (没有圆柱) 的情况下都能正常运行。但是当我加入圆柱体后，压力场突然变得很奇怪, 特别是圆柱体所在的那些行，压力看起来直接消失了, 和周围的压力场不连续。即使网格很大、圆柱体很小，而且位置离边界很远，从第一次迭代开始，压力场就已经不对劲了。我已经尝试了各种方法, 但依然没解决, 非常苦恼, 求助!

问题原因分析

这个问题很可能由以下几个方面的原因共同导致：

圆柱边界条件处理不当: 在交错网格上，速度和压力定义在不同的位置。在处理圆柱这种固体边界时，需要特别小心地施加无滑移和无穿透边界条件。代码中的 Apply_Velocity_Boundry_Conditions 函数似乎尝试处理了圆柱边界，但可能存在逻辑错误或者遗漏了某些情况, 尤其是拐角处的处理。
压力泊松方程求解问题： 压力泊松方程的求解是整个流程中的关键一步。代码中使用 LU 分解来求解线性方程组，但在圆柱存在的情况下，系数矩阵 A 的构建可能存在问题, 这可能是导致圆柱附近的压力出现异常的原因。特别是与flag标记相关的部分。
数值不稳定性/数值耗散： 如果时间步长 dt 过大，或者对流项的离散格式（代码中没明确给出对流项的差分格式，用了自定义的函数例如du2_dx, duv_dy等等）不合适, 则模拟可能出现数值震荡或不稳定，尤其是在高雷诺数情况下。即使开始几步看上去正常，随着时间推移，误差也会逐渐积累，最终导致压力场的错误。
flag矩阵的更新 。在圆柱附近的flag矩阵用于识别流体区域与固体区域, 不正确的flag矩阵也会造成压力和速度数值的不正确。

解决方案

下面分条列出可能的解决方案，并对每种方案进行详细说明：

1. 仔细检查和修正圆柱边界条件

原理和作用:

圆柱表面的边界条件必须同时满足：

无滑移 (No-slip): 流体在固体壁面上的速度为零 (相对于壁面)。
无穿透 (No-penetration): 流体不能流入或流出固体壁面。

在交错网格中，这意味着：

紧邻圆柱表面的、平行于壁面的速度分量必须为零。
紧邻圆柱表面的、垂直于壁面的速度分量必须通过插值或者特殊处理，确保流体不会穿透壁面。

代码示例/修正 (在 Apply_Velocity_Boundry_Conditions 函数内):

def Apply_Velocity_Boundry_Conditions(U,V,U_avg=1.5):
    # ... (之前的边界条件代码保持不变) ...
     #Left Boundry:
    V[:,0] = -V[:,1]                # No - Slip
                       # No - Penetration
    U_max= 1.5 *U_avg
    for j in range(1,j_max+1):
        y = j*dy
        U[j,0] = 4*U_max * (y*(Ly-y)) / Ly**2

    #Right Boundry:
    V[:,-1] = V[:,-2]               # No - Slip
    U[:,-2] = U[:,-3]               # No - Penetration
    

    #Upper Boundry:
    U[-1,:] = -U[-2,:]              # No - Slip
    V[-2,:] = 0                     # No - Penetration

    #Buttom Boundry:
    U[0,:] = -U[1,:]                # No - Slip
    V[0,:] = 0                      # No - Penetration

    #Cylinder boundry condition!!!!!!!
    for j in range(1, j_max):
        for i in range(1, i_max):
            if flag[j, i] == 0:  # Cylinder cell， 现在我们直接找障碍物cell
                # Enforce no penetration by setting normal velocities properly
                # West
                if flag[j, i-1] == 1:  # 西侧是流体
                    U[j,i-1] = -U[j,i]   #保证在(j,i-0.5)处u=0 
                # East
                if flag[j, i+1] == 1:  # 东侧是流体
                   U[j,i] = - U[j,i+1]
                #South
                if flag[j-1, i] == 1:  # 南侧是流体
                    V[j-1,i] = -V[j,i]  #保证在(j-0.5,i)处v=0
                # North
                if flag[j+1, i] == 1:
                   V[j,i] = -V[j+1,i]

    return U, V

主要改动与解释

原代码是寻找流体cell, 并假设这个流体cell临近障碍物cell，现在改为直接寻找flag=0的固体cell,然后根据上下左右的flag判断哪一边是流体。这样更直接和准确。
U[j, i-1] = -U[j,i] 以及类似的处理，保证了在障碍物边界的u,v值为0。

进阶/补充:
上面代码是最基础的处理,还有很多进阶方法：

Immersed Boundary Method (IBM): 一种更高级的处理复杂边界的方法。
Ghost Cell Method: 通过“镜像”固体单元内部的速度来设置边界条件。

2. 重新审查压力泊松方程的离散化

原理和作用：

压力泊松方程：

∇²P = ρ(∂(∇⋅u )/∂t - ∇⋅(u ⊗u ) + ∇⋅ν∇u )

在交错网格上, 对上式进行离散，并整理得到一个线性方程组 Ap * = b ,其中:

p 是一个包含所有网格点压力值的向量。
A 是一个稀疏矩阵，其系数由离散格式 (例如中心差分) 和网格间距决定。
b 是一个向量，包含了速度场的信息 (RHS 右手项)。

正确的离散化是获得准确压力场的关键。

代码示例/检查点 (在 Solving_For_Pressure 函数内)：

系数矩阵 A 的构建: 确保系数矩阵 A 的构建正确反映了圆柱边界附近的情况. 特别关注和 flag 相关的部分. 因为flag 定义了圆柱体的位置. 系数矩阵 A 要根据 flag 来确定哪些点是流体点，哪些点是固体点。固体点上的压力泊松方程不需要求解，可以直接设置为 0 或其他固定值 (例如，设置障碍物内部所有压力都为 0)。

  for i in range(1,i_max+1):
        for j in range(1,j_max+1):
            if flag[j,i]==1:
                #lhs_value = ( P_it1[j, i] - factor * (( P_it1[j, i-1]) / (dx**2) + ( P_it1[j-1, i]) / (dy** 2) ) )
                A[k(j,i),k(j,i)] = 1 + 2*factor/(dx*dx) + 2*factor/(dy*dy)  # 中心点系数
                A[k(j,i),k(j,i-1)] = -factor / dx**2       # West
                A[k(j,i),k(j-1,i)] = -factor / dy**2       # South
                A[k(j,i),k(j,i+1)] = -factor/dx**2         # East, 需保证i+1不会超出索引
                A[k(j,i),k(j+1,i)] = -factor/dy**2        #North, 需保证j+1不会超出索引
                # Right-hand side
                rhs[j,i] = (1/dt) * ( ((F_n[j,i]-F_n[j,i-1])/dx) +((G_n[j,i]-G_n[j-1,i])/dy) )
                b[k(j,i)] =  rhs[j, i] #将常数项放入b向量

            if flag[j,i]==0:
                A[k(j,i), :] = 0  # Zero out the entire row (no influence)
                A[k(j,i), k(j,i)] = 1  # Identity row (pressure fixed)
                b[k(j,i)] = 0  # Ensure no unwanted pressure values

RHS 的计算: 检查 rhs (即方程的右手边) 的计算是否正确, 尤其要注意在圆柱边界附近, F_n 和 G_n 的处理是否和速度边界条件一致.

rhs[j,i] = (1/dt) * ( ((F_n[j,i]-F_n[j,i-1])/dx) +((G_n[j,i]-G_n[j-1,i])/dy) )

边界附近的A矩阵特殊处理 对于固体格点附近, 由于在固体内部压力不需要更新，在固体cell临近流体的地方，A矩阵的构造要做对应的改变, 例如单边差分。这也是上面代码已经尝试解决的问题.

3. 减小时间步长或采用更稳定的数值格式

原理和作用:

时间步长 (dt): 根据 CFL (Courant-Friedrichs-Lewy) 条件，时间步长必须足够小，才能保证数值稳定性。CFL 条件通常要求：
```
dt ≤ C * (min(dx, dy) / max(|u|, |v|))
```
其中 C 是一个小于 1 的常数 (通常取 0.1-0.5)。 calculate_time_step 函数实现了CFL 条件, 但是其中系数 tau=0.2, 已经较小。可以再适当调小tau尝试.
对流项的离散格式: 代码里自定义函数部分,如 du2_dx, duv_dy, 没有给出具体实现，不知道使用的是几阶精度的差分。如果是一阶迎风格式，虽然稳定性好，但数值耗散大, 可能导致精度下降。可以改用中心差分 (二阶精度，但稳定性差) 或其他高阶格式 (例如 QUICK, WENO)。

对于 du2_dx, 一个二阶中心差分的例子:

    def du2_dx(u, j, i, dx):
      d_u2_dx = (u[j, i+1]**2 - u[j,i-1]** 2)/(4*dx) + (u[j,i]*(u[j,i+1] - 2*u[j,i] + u[j,i-1]))/(2 * dx)
      return  d_u2_dx

4. flag 矩阵后处理

原理及作用
原代码中，有对flag矩阵进行后处理，以删除与3个以上流体网格接壤的障碍物网格。但实际上可能会过多的改变了flag的值。而且，后续计算中使用的 num_of_tot_cells_x, num_of_tot_cells_y没有考虑到flag矩阵的修改, 可能会导致计算上的偏差.

代码修改/建议

可以考虑删除flag矩阵后处理相关的代码:

    # Post-processing: Remove inadmissible boundary cells (delete any cylinder cells which border 3 or more fluid cells)
    #for i in range(1, num_of_tot_cells_x):  # Avoid ghost cells
       # for j in range(1, num_of_tot_cells_y):
          #  if flag[j, i] == 0:  # Only check obstacle cells
               # # Count adjacent fluid cells
               # adjacent_fluid = sum([flag[j, i-1], flag[j, i+1], flag[j-1, i], flag[j+1, i]])

               # if adjacent_fluid >= 3:  # More than 2 adjacent fluid cells (inadmissible boundary)
                #    flag[j, i] = 1  # Convert back to fluid to avoid issues

及其后续的 invalid_boundary_found相关的错误检测.