算例收集1:一维Sod激波管问题
1. 问题描述
一维Sod激波管是一个经典的Riemann问题，其控制方程由Euler方程主导。该问题包含了膨胀波、激波和接触间断，求解该问题的主要目就是验证数值计算方法捕捉和分辨流场中间断的能力和计算格式本身的稳定性，几乎所有的差分格式都是建立在求解Riemann问题基础之上的，尔后向多维拓展。
1.1 控制方程
一维Euler方程为：
\(\frac{ \partial q}{\partial t} + \frac{\partial F}{\partial x} = 0\)
其中：
\(q=\begin{pmatrix}
\rho \\\rho u
\\\rho e
\end{pmatrix},F=\begin{pmatrix}
\rho u \\\rho u^2+p
\\u(\rho e+p)
\end{pmatrix}\)
1.2 初始条件
计算域为\(x\in [1,-1]\)，初始条件由下式定义：
\(\left\{\begin{matrix}
\rho =1,u=0,p=1(x<0.5)
\\\rho =0.125,u=0,p=0.1(x\geqslant 0.5)
\end{matrix}\right.\)
计算至t=0.2s
%initalize
x=zeros(nx,1);
for i=1:nx
x(i)= -1.0 + (i-1.0)*dx;
end
%Sod's Riemann problem
%Left side
rhoL = 1.0;
uL = 0.0;
pL = 1.0;
eL = pL/(rhoL*(gama-1.0)) + 0.5*uL*uL;
%Right side
rhoR = 0.125;
uR = 0.0;
pR = 0.1;
eR = pR/(rhoR*(gama-1.0)) + 0.5*uR*uR;
qn=zeros(nx,3);
for i=1:nx
if(x(i)<0.0)
qn(i,1)=rhoL;
qn(i,2)=rhoL*uL;
qn(i,3)=rhoL*eL;
else
qn(i,1)=rhoR;
qn(i,2)=rhoR*uR;
qn(i,3)=rhoR*eR;
end
end
1.3 间断分布过程
间断分解如下图所示，t=0时刻管内速度为零，隔膜左边为高压P1，右边为低压 P2，突然撤去隔膜，初始的压力间断将以非定常正激波的形式向右传播，同时一个等熵非定常中心稀疏波向左传播。在某t>0时刻，管内流体被分成4个区域，A区和D区是波还没有传播到的区域，其未受扰动保持初始间断的左右状态；B区为中心稀疏波传播过后的区域，B区压力为PB，C区为正激波传播后的区域，C区压力后正激波后压力PC，且 PB=PC，但是 B区和C区的温度和密度是不同的，所以两个区是以一个接触间断分开。
2. 数值方法
数值方法就是把控制方程离散的一个过程，简单来说就是把控制方程中的守恒变量q通过计算网格进行离散，网格点的位置对应于其在流场中的位置，网格点上的变量q存储了流场中该点的流动参数，最后通过时间迭代和边界条件求解，得到对应时刻下的流场。
控制方程的离散一般来说有两步，“插值”和“重构”，“插值”是通过相邻网格点的值插值得到网格界面处i ± 1/2的值，“重构”是指利用插值步得的值，采用各种格式（无论是上风型的或是中心型的）构造数值通量，进行离散求解。
本文使用的离散方法是格点格式的有限差分法。
插值方法采用了MUSCL（Monotone Upstream-centred Scheme for Conservation Laws）格式，根据使用的限制器不同，该插值格式能达到的空间精度不同。
%以左边插值为例
function qL = mus_L(nx,q0)
qL=zeros(nx+1,3);
i=1;
q1=q0(i+2,:);
q2=q0(i+1,:);
q3=q0(i,:);
qL(i,:)=ML(q1,q2,q3);
i=2;
q1=q0(i,:);
q2=q0(i-1,:);
q3=q0(i,:);
qL(i,:)=ML(q1,q2,q3);
for i=3:nx
q1=q0(i-2,:);
q2=q0(i-1,:);
q3=q0(i,:);
qL(i,:)=ML(q1,q2,q3);
end
i=nx+1;
q1=q0(i-2,:);
q2=q0(i-1,:);
q3=q0(i-2,:);
qL(i,:)=ML(q1,q2,q3);
end
function qL=ML(q1,q2,q3)
qL=zeros(1,3);
K=1/3;
dp=delta(q3(1,:),q2(1,:));
dm=delta(q2(1,:),q1(1,:));
rl=zeros(1,3);
rl1=zeros(1,3);
for i=1:3
rl(i)=dp(i)/(dm(i)+1e-5);
rl1(i)=dm(i)/(dp(i)+1e-5);
end
v1=limiter(rl);
v2=limiter(rl1);
for i=1:3
qL(i)=q2(1,i) + 1/4 * ((1-K)*v1(1,i) + (1+K)*rl(1,i)*v2(1,i)) * dm(1,i);
end
end
采用van leer限制器,该格式使空间离散精度为二阶。
function v=limiter(R) %van leer limiter
v=zeros(1,3);
for i=1:3
v(i)=(R(i)+abs(R(i)))/(1+R(i));
end
end
重构方法一共选择了4种进行对比，分别是Vanleer格式、HLLC格式、AUSM+格式、AUSMPW+格式。
一维Van leer格式：
function F=Van Leer(nx,gama,fL,fR,qL,qR)
F=zeros(nx+1,3);
fp=zeros(1,3);
fm=zeros(1,3);
for i=1:nx+1
rhoL=qL(i,1);
uL=qL(i,2)/qL(i,1);
eL=qL(i,3)/qL(i,1);
pL=(gama-1)*rhoL*(eL-0.5*uL^2);
aL=sqrt(abs(gama*pL/rhoL));
MaL=uL/aL;
rhoR=qR(i,1);
uR=qR(i,2)/qR(i,1);
eR=qR(i,3)/qR(i,1);
pR=(gama-1)*rhoR*(eR-0.5*uR^2);
aR=sqrt(abs(gama*pR/rhoR));
MaR=uR/aR;
%split Ma
if(MaL>=1)
MaL_P=MaL;
elseif(abs(MaL)<1)
MaL_P=1/4*(MaL+1)^2;
elseif(MaL<=-1)
MaL_P=0.0;
end
if(MaR>=1)
MaR_m=0.0;
elseif(abs(MaR)<1)
MaR_m=1/4*(MaR-1)^2;
elseif(MaR<=-1)
MaR_m=MaR;
end
Ma_mid=MaL_P+ MaR_m;
if(Ma_mid>=1.0)
fp(1,:)=fL(i,:);
fm(1,:)=[0.0,0.0,0.0];
elseif(abs(Ma_mid)<1.0)
f1=1/4*rhoL*aL*(MaL+1)^2;
fp(1,1)=f1;
fp(1,2)=f1*((gama-1)*uL+2*aL)/gama;
fp(1,3)=f1*((gama-1)*uL+2*aL)^2/(2*(gama^2-1));
f2=-1/4*rhoR*aR*(MaR-1)^2;
fm(1,1)=f2;
fm(1,2)=f2*((gama-1)*uR-2*aR)/gama;
fm(1,3)=f2*((gama-1)*uR-2*aR)^2/(2*(gama^2-1));
elseif(Ma_mid<=-1)
fp(1,:)=[0.0,0.0,0.0];
fm(1,:)=fR(i,:);
end
F(i,:)=fp(1,:)+fm(1,:);
end
end
一维HLLC格式：
function F=hllc(nx,gama,fL,fR,qL,qR)
F=zeros(nx+1,3);
D=zeros(1,3);
D(1)=0;
D(2)=1;
for i=1:nx+1
rhoL=qL(i,1);
uL=qL(i,2)/qL(i,1);
eL=qL(i,3)/qL(i,1);
pL=(gama-1)*rhoL*(eL-0.5*uL^2);
aL=sqrt(abs(gama*pL/rhoL));
rhoR=qR(i,1);
uR=qR(i,2)/qR(i,1);
eR=qR(i,3)/qR(i,1);
pR=(gama-1)*rhoR*(eR-0.5*uR^2);
aR=sqrt(abs(gama*pR/rhoR));
SL=min(uL,uR)-max(aL,aR);
SR=min(uL,uR)+max(aL,aR);
SM=(pR-pL+rhoL*uL*(SL-uL)-rhoR*uR*(SR-uR))/(rhoL*(SL-uL)-rhoR*(SR-uR));
pLR=0.5*(pL+pR+rhoL*(SL-uL)*(SM-uL)+rhoR*(SR-uR)*(SM-uR));
D(3)=SM;
if(SL>=0.0)
F(i,:)=fL(i,:);
elseif(SR<=0.0)
F(i,:)=fR(i,:);
elseif((SL<=0.0)&&(SM>=0.0))
F(i,:)=(SM*(SL*qL(i,:)-fL(i,:))+SL*pLR*D)/(SL-SM);
elseif((SR>=0.0)&&(SM<=0.0))
F(i,:)=(SM*(SR*qR(i,:)-fR(i,:))+SR*pLR*D)/(SR-SM);
end
end
end
一维AUSM+格式：
function F=AUSM_plus(nx,gama,qL,qR)
F=zeros(nx+1,3);
for i=1:nx+1
rhoL=qL(i,1);
uL=qL(i,2)/qL(i,1);
eL=qL(i,3)/qL(i,1);
pL=(gama-1)*rhoL*(eL-0.5*uL^2);
aL=sqrt(2*(gama-1)/(gama+1)*(eL+pL/rhoL));
rhoR=qR(i,1);
uR=qR(i,2)/qR(i,1);
eR=qR(i,3)/qR(i,1);
pR=(gama-1)*rhoR*(eR-0.5*uR^2);
aR=sqrt(2*(gama-1)/(gama+1)*(eR+pR/rhoR));
aL_ave=(aL^2)/max(abs(uL),aL);
aR_ave=(aR^2)/max(abs(uR),aR);
a_mid=min(aL_ave,aR_ave);
MaL=uL/a_mid;
MaR=uR/a_mid;
[phiL,PL_V]=AUSM_turnvar(gama,qL(i,:));
[phiR,PR_V]=AUSM_turnvar(gama,qR(i,:));
% Mach number and pressure splitting functions
M0L=ausmsplit_Ma(MaL);
M0R=ausmsplit_Ma(MaR);
P0L=ausmsplit_P(MaL);
P0R=ausmsplit_P(MaR);
for j=1:3
F(i,j)=a_mid*(M0L(1)*phiL(j)+M0R(2)*phiR(j))+...
P0L(1)*PL_V(j)+P0R(2)*PR_V(j);
end
end
end
一维AUSMPW+格式
function F=AUSMPW_plus(nx,gama,qL,qR)
F=zeros(nx+1,3);
alpha=0;
for i=1:nx+1
rhoL=qL(i,1);
uL=qL(i,2)/qL(i,1);
eL=qL(i,3)/qL(i,1);
pL=(gama-1)*rhoL*(eL-0.5*uL^2);
aL=sqrt(2*(gama-1)/(gama+1)*(eL+pL/rhoL));
rhoR=qR(i,1);
uR=qR(i,2)/qR(i,1);
eR=qR(i,3)/qR(i,1);
pR=(gama-1)*rhoR*(eR-0.5*uR^2);
aR=sqrt(2*(gama-1)/(gama+1)*(eR+pR/rhoR));
aL_ave=aL^2/max(abs(uL),aL);
aR_ave=aR^2/max(abs(uR),aR);
a_mid=min(aL_ave,aR_ave);
MaL=uL/a_mid;
MaR=uR/a_mid;
[phiL,PL_V]=AUSM_turnvar(gama,qL(i,:));
[phiR,PR_V]=AUSM_turnvar(gama,qR(i,:));
% Mach number and pressure splitting functions
M0L=split_Ma(MaL);
M0R=split_Ma(MaR);
P0L=split_P(alpha,MaL);
P0R=split_P(alpha,MaR);
omega=1-power(min(pL/pR,pR/pL),3);
ps=P0L(1)*pL+P0R(2)*pR;
if(abs(MaL)<1)
f0l=pL/ps-1.0;
else
f0l=0.0;
end
if(abs(MaR)<1)
f0r=pR/ps-1.0;
else
f0r=0.0;
end
m_mid=M0L(1)+M0R(2);
if(m_mid>=0)
M0LP_ave=M0L(1)+M0R(2)*((1-omega)*(1+f0r)-f0l);
M0Rm_ave=M0R(2)*omega*(1+f0r);
else
M0LP_ave=M0L(1)*omega*(1+f0l);
M0Rm_ave=M0R(2)+M0L(1)*((1-omega)*(1+f0l)+f0l-f0r);
end
for j=1:3
F(i,j)=a_mid*(M0LP_ave*phiL(j)+M0Rm_ave*phiR(j))+...
P0L(1)*PL_V(j)+P0R(2)*PR_V(j);
end
end
end
时间推进采用显式三阶Runge-Kutta格式。
%三阶Runge-Kutta格式;
for i=1:nt
rhs1=calculate.rhs(nx,dx,qn);
for j=1:nx
qt(j,:)=qn(j,:)+dt*rhs1(j,:);
end
rhs2=calculate.rhs(nx,dx,qt);
for j=1:nx
qt(j,:)=0.75*qn(j,:)+0.25*qt(j,:)+0.25*dt*rhs2(j,:);
end
rhs3=calculate.rhs(nx,dx,qt);
for j=1:nx
qt(j,:)=(1.0/3.0)*qn(j,:)+(2.0/3.0)*qt(j,:)+(2.0/3.0)*dt*rhs3(j,:);
end
qn=qt;
end
function R=rhs(nx,dx,qt)
gama=1.4;
qL=mus_L(nx,qt);
qR=mus_R(nx,qt);
fL=flux(nx,gama,qL);
fR=flux(nx,gama,qR);
%F=hllc(nx,gama,fL,fR,qL,qR);
%F=AUSMPW_plus(nx,gama,qL,qR);
F=AUSM_plus(nx,gama,qL,qR);
%F=VanLeer(nx,gama,fL,fR,qL,qR);
R=zeros(nx,3);
for i=1:nx
R(i,:)=-(F(i+1,:)-F(i,:))/dx;
end
end
3. 计算结果
【左上：密度，右上：速度，左下：能量，右下：压力】
Van Leer通量重构格式的计算结果：
HLLC通量重构格式的计算结果：
AUSM+通量重构格式的计算结果：
AUSMPW+通量重构格式的计算结果：
参考解：
这里可以看出HLLC、AUSM+、AUSMPW+和参考解的计算结果基本吻合，AUSMPW+格式稳定性好，Van Leer格式在A区和B区之间稀疏波传播的空间中出现一定振荡。
参考文献：潘沙.高超声速气动热数值模拟方法及大规模并行计算研究[D].长沙:国防科学技术大学,2010,4.
浙公网安备 33010602011771号