算例收集7：Rayleigh-Taylor不稳定问题

1.问题描述
Rayleigh-Taylor不稳定性（RTI）算例是计算流体力学（CFD）中验证多相流、界面捕捉和混合过程数值方法的核心基准测试。它模拟了重力场（或加速度场）中密度梯度与压力梯度方向相反时发生的流体失稳现象。该算例包含了流动不连续性和多尺度小结构，模拟了重力场（或加速度场）中密度梯度与压力梯度方向相反时发生的流体失稳现象。
通过密度分层+反向加速度的简单设置，揭示了流体失稳→涡演化→湍流混合的全过程。

1.1控制方程
控制方程采用二维Euler方程，由于算例中包含了上下两层密度不同的气体介质，考虑到重力场对流体运动的影响，在方程右端添加了源项：
\(\frac{ \partial q}{\partial t} + \frac{\partial F}{\partial x} + \frac{\partial G}{\partial y}= S\)
其中：
\(q=\begin{pmatrix} \rho \\\rho u \\\rho v \\\rho e \end{pmatrix},F=\begin{pmatrix} \rho u \\\rho u^2+p \\\rho uv \\u(\rho e+p) \end{pmatrix},G=\begin{pmatrix} \rho v \\\rho uv \\\rho v^2+p \\v(\rho e+p) \end{pmatrix},S=\begin{pmatrix} 0 \\0 \\\rho \\\rho v \end{pmatrix}\)

1.2 初始条件及边界条件
计算区域的大小为[0,0.25]x[0,1]，以y=1/2为分解，下层为重介质，上层为轻介质。上下两层介质的\((\rho , u,v,p)\)分别为：

其中\(c=\sqrt{\gamma p/\rho}\),\(\gamma=5/3\)。
左、右边界设置为反射边界条件，上、下边界的值分别固定为（1, 0, 0, 2.5）和（2, 0, 0, 1）
计算至\(t_{end}==1.95s\)

gama=5.0/3.0
nx=128
ny=512

#CONTIUE
calculated_times = 60000   #10000 + 20000 + 30000

tf=1.95
dt=2.5e-5
num_iter=round(Int,tf/dt) - calculated_times

xb=0
xe=0.25
yb=0
ye=1.0
dx=(xe-xb)/(nx-1)
dy=(ye-yb)/(ny-1)

x=Array{Float64}(undef,nx)
y=Array{Float64}(undef,ny)
qn=Array{Float64}(undef,4,ny,nx)
for i=1:nx
  x[i]=(i-1)*dx
end
for j=1:ny
  y[j]=(j-1)*dy
end

if(calculated_times == 0)
  for i=1:nx
      global pr
      global rho
      global u
      global v
      for j=1:ny
          if(y[j] < 1/2.0)
              pr = 2.0*y[j]+1.0
              rho = 2.0
              a=sqrt(gama*pr/rho)
              u = 0.0
              v = - 0.025*a*cos(8.0*pi*x[i])
              q1=turn_var(gama,rho,u,v,pr)
              qn[:,j,i] = q1
          else
              pr = y[j]+ 3/2.0
              rho = 1.0
              a=sqrt(gama*pr/rho)
              u = 0.0
              v = - 0.025*a*cos(8.0*pi*x[i])
              q2=turn_var(gama,rho,u,v,pr)
              qn[:,j,i] = q2
          end
      end
  end
else
  qn=read_txt(nx,ny,calculated_times)
end
#bottom bc
pr = 1.0
rho = 2.0
u = 0.0
v = 0
q_bc1=turn_var(gama,rho,u,v,pr)
#top bc
pr = 2.5
rho = 1.0
u = 0.0
v = 0
q_bc2=turn_var(gama,rho,u,v,pr)

2.数值方法
采用格点格式有限差分法。
2.1 插值格式
采用5阶精度的WENO-JS格式进行插值，计算单元的左右界面的守恒变量。
2.2 重构格式
采用AUSMPW+格式。
2.3 时间推进
采用3步Runge-Kutta格式。

qt=Array{Float64}(undef,4,ny,nx)
for iter=1:num_iter
  rhs1=rhs(nx,ny,dx,dy,qn,q_bc1,q_bc2)
  for k=1:4
      for j=1:ny
          for i=1:nx
              qt[k,j,i]=qn[k,j,i]+dt*rhs1[k,j,i]
          end           
      end
  end

  rhs2=rhs(nx,ny,dx,dy,qn,q_bc1,q_bc2)
  for k=1:4
      for j=1:ny
          for i=1:nx
              qt[k,j,i]=0.75*qn[k,j,i]+0.25*qt[k,j,i]+0.25*dt*rhs2[k,j,i]
          end           
      end
  end  

  rhs3=rhs(nx,ny,dx,dy,qn,q_bc1,q_bc2)
  for k=1:4
      for j=1:ny
          for i=1:nx
              qt[k,j,i]=(1.0/3.0)*qn[k,j,i]+(2.0/3.0)*qt[k,j,i]+(2.0/3.0)*dt*rhs3[k,j,i]
          end           
      end
  end  
  
  if(mod(iter,5)==0)
      norm = compute_l2norm(nx,ny,qn,qt)
      println("iter==",iter,"  ","norm==",norm)
  end 

  for k=1:4
      for j=1:ny
          for i=1:nx
              qn[k,j,i]=qt[k,j,i]
          end           
      end
  end 
  
end

3.计算结果

计算时间等于1.95s时本文计算方法得到的流场密度分布云图如下：

密度分布云图，计算网格点数为128x512

下面给出参考文献中的计算结果：

在相同网格条件下，文献中的计算结果更为精细，捕捉到很多的细小涡结构，因此本算例采用的计算方法耗散较大，流场中的很多细节还不能表现出来，对于这类流体失稳问题的计算不是很适用。

参考文献
T. Yang, G. Zhao, Q. Zhao.Novel TENO schemes with improved accuracy order based on perturbed polynomial reconstruction[J]. Journal of Computational Physics 488 (2023) 112219.