返回

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

python

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

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

问题原因分析

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

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

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

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

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

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

  1. 系数矩阵 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
  1. RHS 的计算: 检查 rhs (即方程的右手边) 的计算是否正确, 尤其要注意在圆柱边界附近, F_nG_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) )
  1. 边界附近的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相关的错误检测.

总结

按顺序尝试以下方法:

  1. 首先, 重点检查圆柱边界条件, 应用本文中修改后的 Apply_Velocity_Boundry_Conditions.
  2. 审查压力泊松方程中的A矩阵与flag矩阵部分, 和rhs 的计算.
  3. 修改flag矩阵的更新方法 (可以先尝试删除相关代码), 或者在后续计算时考虑 flag 的改变.
  4. 减小时间步长, 修改对流项差分格式。

祝你好运, 希望这些建议能帮你解决圆柱绕流模拟中的压力问题!