算例收集2：二维平板算例
1.问题描述
平板上的超声速流动：完整的Navier-Stokes方程的数值求解。
考虑零攻角尖前缘平板上的超声速流动，平板长度为L，如下图所示，层流边界层在平板前缘产生，并且在低Reynolds数时保持层流。由于粘性边界层的存在，平板好像具有一定曲率一样，因此会在前缘产生弯曲的激波。
1.1控制方程
二维的Navier-Stokes方程形式如下（忽略体积力和体积热）：
$ \frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x}(\rho u)+\frac{\partial }{\partial y}(\rho v)=0$
$ \frac{\partial }{\partial t}(\rho u)+\frac{\partial }{\partial x}(\rho u^{2}+p-\tau _{xx})+\frac{\partial }{\partial y}(\rho uv-\tau _{xy})=0$
$ \frac{\partial }{\partial t}(\rho v)+\frac{\partial }{\partial x}(\rho uv-\tau _{xy})+\frac{\partial }{\partial y}(\rho v^{2}+p-\tau _{yy})=0$
\(\frac{\partial }{\partial t}(E_{t})+\frac{\partial }{\partial x}[(E_{t}+p)u+q_x-u\tau _{xx}-v\tau_{xy}]+\frac{\partial }{\partial y}[(E_{t}+p)v+q_y-u\tau _{yx}-v\tau_{yy}]=0\)
以上的方程中,\(E_{t}\)是单位体积动能和内能的和，定义如下：
\(E_t=\rho(e+\frac{V^2}{2})\)
剪应力和正应力由速度梯度得到：
\(\tau_{xx}=\lambda (\nabla\bullet \mathbf{V})+2\mu \frac{\partial u}{\partial x}\)
\(\tau_{yy}=\lambda (\nabla\bullet \mathbf{V})+2\mu \frac{\partial v}{\partial y}\)
\(\tau_{xy}=\tau_{yx}=\mu (\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})\)
导热矢量由傅里叶定律得到：
\(q_x=-k\frac{\partial T}{\partial x}\)
\(q_y=-k\frac{\partial T}{\partial y}\)
上述方程组共有9个未知数\(\rho ,u,v\left | V\right |,p,T,e,\mu,k\)
为使方程封闭还需要5个方程：
1）完全气体状态方程
\(p=\rho RT\)
2）内能
\(e=c_vT\)
3）合速度
\(\left | V\right |=\sqrt{u^2+v^2}\)
4）由Sutherland公式计算粘性系数
\(\mu =\mu _0(\frac{T}{T_0})^{3/2}\frac{T_0+110}{T+110}\)
其中\(\mu _0\)和\(T_0\)是标准海平面处的参考值。
5）由Prandtl数计算热传导系数
\(Pr=0.71=\frac{\mu c_p}{k}\)
现将9个方程改写为矢量形式：
$ \frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{E}}{\partial x}+\frac{\partial \mathbf{F}}{\partial y}=0$
\(\mathbf{U}=\begin{Bmatrix}
\rho \\\rho u
\\\rho v
\\E_t
\end{Bmatrix}\),\(\mathbf{E}=\begin{Bmatrix}
\rho u \\\rho u^2+p-\tau_{xx}
\\\rho uv-\tau_{xy}
\\(E_t+p)u-u\tau_{xx}-v\tau_{xy}+q_x
\end{Bmatrix}\),\(\mathbf{F}=\begin{Bmatrix}
\rho v \\\rho uv-\tau_{yx}
\\\rho v^2+p-\tau_{yy}
\\(E_t+p)v-u\tau_{yx}-v\tau_{yy}+q_y
\end{Bmatrix}\)
1.2 初始条件和边界条件
计算域如下图所示，平板长度为0.00001m，来流马赫数4，Re约为1000。
类型1：前缘((IMIN,JMIN)或(1,1))处给定无滑移边界 (u_((1,1))=0,0) , 温度 T_((1,1) 和压力p,分别为自由来流值。
类型2：左边界处（不包括前缘）和上边界处的x方向速度为u,温度和压力假设 为自由来流值；y方向速度v设为0。
类型3：平板表面处，速度为无滑移条件(u=v=0.0)。温度（除了前缘处）等于 壁面温度T。壁面的压力（除了前缘处）通过内点(j=2,j=3)的值外插得到。 例如：$p_{(i,1)}=2p_{(i,2)}-p_{(i,3)}$
类型4：最后，右边界（不包括JMIN=1,JMAX=70）的所有参数由j位置相同 的两个内点外插得到。 例如：$u_{(IMAX,j)}=2u_{(IMAX-1,j)}-u_{(IMAX-2,j)}$
1.3时间步长
已知平板长度为LHORI，x方向步长表示为：
\(\Delta x=\frac{LHORI}{IMAX-1}\)
由Blasius尾缘的计算可以预测，计算区域y方向的高度至少为边界层5倍以
上，才能满足计算。
\(LVERT=5\times \delta ,\delta = \frac{5LHORI}{\sqrt{R_{eL}}}\)
\(\Delta y=\frac{LVERT}{JMAX-1}\)
随后计算每个点的网格Re数：
因为使用的方法是显式格式，时间步长由稳定条件决定。为了确定时间步
长，使用下面的Courant-Friedrichs-Lewy(CFL)条件 ：
\(\Delta t=min[K(\Delta t_{CFL})_{i,j}]\),一般取 0.5≤K≤0.8
2.数值方法
计算方法还是采用格点格式的有限差分法。
这里讨论的通量分裂格式（也就是之前说的重构方法）采用MacCormack方法，对于插值格式没用采用任何特殊的格式，主要是一阶的前差、后差和中心差分，最后通过MacCormack方法的整合达到空间二阶精度。
MacCormack方法简单来就4步：
1）由已知的t时刻流场，通过前面给出的矢量形式控制方程右边的空间向前
差分得\((\frac{\partial \textbf{U}}{\partial t})_{i,j}^{t}\)
2) 由第1步的计算值求解t+Δt时刻的预测值
3）由第2步计算的预测值使用后向差分计算得到\(\overline{(\frac{\partial \textbf{U}}{\partial t})}_{i,j}^{t+\Delta t}\)
4）最后
3.程序结构
求解平板流动问题的NS方程的程序结构如下
1）设定初始条件
Ma=4.0; %mach number
lhori=1E-5; %plate length
a_far=340.28; %sonic speed
p_far=101325.0; %pressure
T_far=288.16; %temperature
T_w=T_far;
gama=1.4; %heat ratio
PR=0.71; %prandtl number
mu_ref=1.7849E-5; %viscosity ref
T_ref=288.16; %temperature ref
mu_far=mu_ref*(T_far/T_ref)^(3/2)*((T_ref+110)/(T_far+110)); %viscosity
R_con=287; %gas constant
rho_far=p_far/(R_con*T_far); %density
v_far=Ma*sqrt(gama*R_con*T_far); %velocity
Re_far=rho_far*v_far*lhori/mu_ref; %Re number
c_v=R_con/(gama-1);
c_p=gama*c_v;
e_far=c_v*T_far;
点击查看代码
%grid
nx=70;
ny=70;
dx=lhori/(nx-1);
delta=5*lhori/sqrt(Re_far);
levert=5*delta;
dy=levert/(ny-1);
x=zeros(nx,ny);
y=zeros(nx,ny);
for j=1:ny
for i=1:nx
x(i,j)=dx*(i-1);
y(i,j)=dy*(j-1);
end
end
%calculate
q=zeros(4,nx,ny);
%initialize
for j=1:ny
for i=1:nx
q(1,i,j)=v_far; %u
q(2,i,j)=0.0; %v
q(3,i,j)=p_far; %p
q(4,i,j)=T_far; %T
end
end
%wall
for i=1:nx
q(1,i,1)=0.0;
q(2,i,1)=0.0;
q(3,i,1)=p_far;
q(4,i,1)=T_w;
end
2）TSTEP:确定合理的时间步长
function DT=Timestep(nx,ny,dx,dy,gama,q0)
PR=0.71;
vv=zeros(nx,ny);
for j=1:ny
for i=1:nx
mu=Viscosity(q0(4,i,j));
vv(i,j)=4/3*mu*(gama*mu/PR)/(q0(3,i,j)/(287*q0(4,i,j)));
end
end
v_ij=max(max(vv));
dt=zeros(nx,ny);
for j=1:ny
for i=1:nx
aa=sqrt(gama*287*q0(4,i,j));
dt(i,j)=(abs(q0(1,i,j))/dx+abs(q0(2,i,j))/dy...
+aa*sqrt(1/(dx^2)+1/(dy^2))+2*v_ij*(1/(dx^2)+1/(dy^2)))^(-1);
end
end
K=0.75;
DT=K*min(min(dt));
end
3)MAC(MacCormack):使用预测一修正方法更新(i,j)点的流场参数
function U4=MacCormack(nx,ny,dx,dy,c_v,p_f,T_f,T_w,u_f,dt,q1)
点击查看代码
function U4=MacCormack(nx,ny,dx,dy,c_v,p_f,T_f,T_w,u_f,dt,q1)
%predict step
U1=turnvar1(nx,ny,c_v,q1);
E1=turnvar2(nx,ny,dx,dy,c_v,U1,2.0,2.0);
F1=turnvar3(nx,ny,dx,dy,c_v,U1,4.0,2.0);
U_pre=zeros(4,nx,ny); %x--forwad y--forwad
for j=2:ny-1
for i=2:nx-1
U_pre(:,i,j)=U1(:,i,j)-(dt/dx)*(E1(:,i+1,j)-E1(:,i,j))...
-(dt/dy)*(F1(:,i,j+1)-F1(:,i,j));
end
end
U2=BC(nx,ny,p_f,T_f,T_w,u_f,c_v,U_pre);
%correct step
E2=turnvar2(nx,ny,dx,dy,c_v,U2,1.0,1.0);
F2=turnvar3(nx,ny,dx,dy,c_v,U2,3.0,1.0);
U3=zeros(4,nx,ny); %x--back y--baCK
for j=2:ny-1
for i=2:nx-1
U3(:,i,j)=0.5*(U1(:,i,j)+U2(:,i,j)...
-(dt/dx)*(E2(:,i,j)-E2(:,i-1,j))-(dt/dy)*(F2(:,i,j)-F2(:,i,j-1)));
end
end
U4=BC(nx,ny,p_f,T_f,T_w,u_f,c_v,U3);
end
4）CONVER:检测流场是否收敛
%residual
sum=0.0;
for j=1:ny
for i=1:nx
rho1=q(3,i,j)/(R_con*q(4,i,j));
rho2=q_new(3,i,j)/(R_con*q_new(4,i,j));
sum=sum+(rho2-rho1)^2;
end
end
res=sqrt(sum/((nx-1)*(ny-1)));
if(res<1E-15)
break
end
5)DYNVIS和THERMC是函数子程序,调用它们来求解(i,j)点的动力黏性系数和热传导系数。主程序只在流场初始化的过程中调用这些子程序。而MAC每次调用时都使用这些函数程。
function mu=Viscosity(T)
mu_R=1.7849E-5;
T_R=288.16;
mu=mu_R*power((T/T_R),3/2)*((T_R+110)/(T+110));
end
function K=Heatconduc(T,c_p)
PR=0.71;
mu=Viscosity(T);
K=mu*c_p/PR;
end
6）黏性影响由以下5个函数描述：TAUXX,TAUXY,TAUYY,QX和QY。当需要确定某剪切力或热传导项时，对应的函数就被调用。
function tao_diag=tao(nx,ny,dx,dy,q0,index,dir)
点击查看代码
phiux=zeros(nx,ny);
phivy=zeros(nx,ny);
if(index==1.0) %x-forward y-center
for j=1:ny
for i=1:nx
if(i==nx)
phiux(i,j)=(q0(1,i,j)-q0(1,i-1,j))/dx;
else
phiux(i,j)=(q0(1,i+1,j)-q0(1,i,j))/dx;
end
if(j==1)
phivy(i,j)=(q0(2,i,j+1)-q0(2,i,j))/dy;
elseif(1<j)&&(j<ny)
phivy(i,j)=(q0(2,i,j+1)-q0(2,i,j-1))/(2*dy);
else
phivy(i,j)=(q0(2,i,j)-q0(2,i,j-1))/dy;
end
end
end
elseif(index==2.0) %x-back y-center
for j=1:ny
for i=1:nx
if(i==1)
phiux(i,j)=(q0(1,i+1,j)-q0(1,i,j))/dx;
else
phiux(i,j)=(q0(1,i,j)-q0(1,i-1,j))/dx;
end
if(j==1)
phivy(i,j)=(q0(2,i,j+1)-q0(2,i,j))/dy;
elseif(1<j)&&(j<ny)
phivy(i,j)=(q0(2,i,j+1)-q0(2,i,j-1))/(2*dy);
else
phivy(i,j)=(q0(2,i,j)-q0(2,i,j-1))/dy;
end
end
end
elseif(index==3.0) %x-center y-forward
for j=1:ny
for i=1:nx
if(i==1)
phiux(i,j)=(q0(1,i+1,j)-q0(1,i,j))/dx;
elseif(1<i)&&(i<nx)
phiux(i,j)=(q0(1,i+1,j)-q0(1,i-1,j))/(2*dx);
else
phiux(i,j)=(q0(1,i,j)-q0(1,i-1,j))/dx;
end
if(j==ny)
phivy(i,j)=(q0(2,i,j)-q0(2,i,j-1))/dy;
else
phivy(i,j)=(q0(2,i,j+1)-q0(2,i,j))/dy;
end
end
end
else %x-center y-back
for j=1:ny
for i=1:nx
if(i==1)
phiux(i,j)=(q0(1,i+1,j)-q0(1,i,j))/dx;
elseif(1<i)&&(i<nx)
phiux(i,j)=(q0(1,i+1,j)-q0(1,i-1,j))/(2*dx);
else
phiux(i,j)=(q0(1,i,j)-q0(1,i-1,j))/dx;
end
if(j==1)
phivy(i,j)=(q0(2,i,j+1)-q0(2,i,j))/dy;
else
phivy(i,j)=(q0(2,i,j)-q0(2,i,j-1))/dy;
end
end
end
end
tao_diag=zeros(nx,ny);
if(dir==1.0) %t_xx
for j=1:ny
for i=1:nx
mu=Viscosity(q0(4,i,j));
tao_diag(i,j)=2*mu*phiux(i,j)-2/3*mu*(phiux(i,j)+phivy(i,j));
end
end
end
if(dir==2.0) %t_yy
for j=1:ny
for i=1:nx
mu=Viscosity(q0(4,i,j));
tao_diag(i,j)=2*mu*phivy(i,j)-2/3*mu*(phiux(i,j)+phivy(i,j));
end
end
end
function Q_x=Qx(nx,ny,dx,c_p,q0,index)
点击查看代码
Q_x=zeros(nx,ny);
if(index==1.0) %x--forward;
for j=1:ny
for i=1:nx
kk=Heatconduc(q0(4,i,j),c_p);
if(i==nx)
Q_x(i,j)=-kk*(q0(4,i,j)-q0(4,i-1,j))/dx;
else
Q_x(i,j)=-kk*(q0(4,i+1,j)-q0(4,i,j))/dx;
end
end
end
end
if(index==2.0) %x--back;
for j=1:ny
for i=1:nx
kk=Heatconduc(q0(4,i,j),c_p);
if(i==1)
Q_x(i,j)=-kk*(q0(4,i+1,j)-q0(4,i,j))/dx;
else
Q_x(i,j)=-kk*(q0(4,i,j)-q0(4,i-1,j))/dx;
end
end
end
end
7）当内点的流场参数确定后（不论在预测或是修正步中），边界条件通过调用子程序BC确定。
function UU=BC(nx,ny,p_f,T_f,T_w,u_f,c_v,U0)
点击查看代码
UU=zeros(4,nx,ny);
for j=2:ny-1
for i=2:nx-1
UU(:,i,j)=U0(:,i,j);
end
end
%BC---1
u=0.0;
v=0.0;
p=p_f;
T=T_f;
rho=p/(287*T);
e=c_v*T;
UU(1,1,1)=rho;
UU(2,1,1)=rho*u;
UU(3,1,1)=rho*v;
UU(4,1,1)=rho*(e+(u^2+v^2)/2);
%BC---2
for j=2:ny
u=u_f;
v=0.0;
p=p_f;
T=T_f;
rho=p/(287*T);
e=c_v*T;
UU(1,1,j)=rho;
UU(2,1,j)=rho*u;
UU(3,1,j)=rho*v;
UU(4,1,j)=rho*(e+(u^2+v^2)/2);
end
for i=1:nx
u=u_f;
v=0.0;
p=p_f;
T=T_f;
rho=p/(287*T);
e=c_v*T;
UU(1,i,ny)=rho;
UU(2,i,ny)=rho*u;
UU(3,i,ny)=rho*v;
UU(4,i,ny)=rho*(e+(u^2+v^2)/2);
end
%BC---4
qq=zeros(4,nx,ny);
for j=2:ny-1
for i=2:nx-1
u=U0(2,i,j)/U0(1,i,j);
v=U0(3,i,j)/U0(1,i,j);
ee=U0(4,i,j)/U0(1,i,j)-(u^2+v^2)/2.0;
T=ee/c_v;
p=U0(1,i,j)*287*T;
qq(1,i,j)=u;
qq(2,i,j)=v;
qq(3,i,j)=p;
qq(4,i,j)=T;
end
end
for j=2:ny-1
u=2*qq(1,nx-1,j)-qq(1,nx-2,j);
v=2*qq(2,nx-1,j)-qq(2,nx-2,j);
p=2*qq(3,nx-1,j)-qq(3,nx-2,j);
T=2*qq(4,nx-1,j)-qq(4,nx-2,j);
rho=p/(287*T);
e=c_v*T;
UU(1,nx,j)=rho;
UU(2,nx,j)=rho*u;
UU(3,nx,j)=rho*v;
UU(4,nx,j)=rho*(e+(u^2+v^2)/2);
end
%BC---3
for i=2:nx
u=0.0;
v=0.0;
if(i<nx)
p=2*qq(3,i,2)-qq(3,i,2);
else
e1=UU(4,i,2)/UU(1,i,2)-((UU(2,i,2)/UU(1,i,2))^2+(UU(3,i,2)/UU(1,i,2))^2)/2.0;
T1=e1/c_v;
p1=UU(1,i,2)*287*T1;
e2=UU(4,i,3)/UU(1,i,3)-((UU(2,i,3)/UU(1,i,3))^2+(UU(3,i,3)/UU(1,i,3))^2)/2.0;
T2=e2/c_v;
p2=UU(1,i,3)*287*T2;
p=2*p1-p2;
end
T=T_w;
rho=p/(287*T);
e=c_v*T;
UU(1,i,1)=rho;
UU(2,i,1)=rho*u;
UU(3,i,1)=rho*v;
UU(4,i,1)=rho*(e+(u^2+v^2)/2);
end
4.计算结果
流场的无量纲压力分布云图：
无量纲温度分布云图：
平板表面的无量纲压力分布：
本文计算时壁面的边界条件为等温边界，也就是常温边界条件。本文计算的壁面压力分布整体上与文献基本一致，除了在平板前缘处压力峰值上，本文计算的无量纲压力峰值在2.5，而文献的结果大约在2.8~2.9左右。总体上可以说明算法可行性😝 😝 😝。
参考文献： John D. Anderson.计算流体力学入门[M].北京:清华大学出版社,2010.
浙公网安备 33010602011771号