Fortran一维Euler方程CFD代码NaN问题技术求助

问题本质 #

计算croe时平方根参数为负，根源是流动变量（密度、压力、内能）出现非物理解，导致声速计算崩溃。

具体修复措施 #

• 物理约束变量范围：在声速计算前强制限制变量的物理下限，避免负数值：

  ! 插入到croe计算前的代码段
  rho = max(rho, 1.0e-6)  ! 防止密度为负或趋近于0
  p = max(p, 1.0e-6)      ! 保证压力为正
  ! 重新计算内能并约束下限
  e_int = p / ((gamma - 1.0d0) * rho)
  e_int = max(e_int, 1.0e-6)
  

• 动态调整时间步长：固定时间步长可能在流动突变区（如激波、喉部）超出数值稳定域，改用动态CFL数：

  ! 替换原固定dt的计算逻辑
  real :: cfl, max_speed, dt
  cfl = 0.4  ! 初始CFL可设为0.4，若仍出问题再降至0.2
  max_speed = maxval(u + sqrt(gamma * p / rho))  ! 域内最大流速+声速
  dt = cfl * dx / max_speed
  

• 修正Roe平均的声速计算：若使用Roe格式，添加熵修正避免跨激波时的数值异常：

  ! 给平方根内的参数添加小量偏移
  real :: delta
  delta = 1.0e-6 * max(abs(p_left - p_right), abs(rho_left - rho_right))
  croe = sqrt( max( (gamma - 1.0d0)*(e_left + e_right - 0.5d0*(u_left**2 + u_right**2)), delta ) )
  

• 检查边界条件：确认喷管进出口边界（如出口压力边界、喉部对称边界）的实现是否正确，避免非物理的边界反馈干扰内部流动。

快速调试手段 #

在迭代到异常步骤时输出变量值，精准定位问题源头：

if (ntime >= 53 .and. k == 2) then
  write(*, '(A,I3,A,I2,A,4ES12.4)') 'ntime=', ntime, ', k=', k, ', rho/u/p/e=', rho, u, p, e
endif

内容的提问来源于stack exchange，提问作者Γεώργιος Γιουρούκος