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
相关的错误检测.
总结
按顺序尝试以下方法:
- 首先, 重点检查圆柱边界条件, 应用本文中修改后的
Apply_Velocity_Boundry_Conditions
. - 审查压力泊松方程中的A矩阵与flag矩阵部分, 和rhs 的计算.
- 修改flag矩阵的更新方法 (可以先尝试删除相关代码), 或者在后续计算时考虑
flag
的改变. - 减小时间步长, 修改对流项差分格式。
祝你好运, 希望这些建议能帮你解决圆柱绕流模拟中的压力问题!