b一维可压ns方程b0814093738.docx
《b一维可压ns方程b0814093738.docx》由会员分享，可在线阅读，更多相关《b一维可压ns方程b0814093738.docx（24页珍藏版）》请在冰豆网上搜索。
b一维可压ns方程b0814093738
附录B一维可压缩黏性流动问题的数值解法与计算程序
一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风型TVD差分算法进行数值求解。
同时，为了初学者入门和练习方便,这里给出了由C语言和Fortran7语言编写的、计算一维可压缩黏性流动问题的计算程序，供大家学习参考。
B-1利用二阶迎风型TVD差分格式求解一维可压缩黏性流动问题1.一维可压缩黏性流动问题
在两端开口管道中充满了可压缩黏性流体。
当黏性流体以超声速
从左向右运动时，一定会在管道中形成一道正激波，如图B.1所示。
1、Ui、Pi和2、U2、P2分别为激波波前和波后的参数。
该问题可简化为一维可压缩黏性流动问题。
当数值解达到稳定时，在管道中可求解得到一道稳定的激波。
r
1
|-
k
1
p2—1
1
1U{-
1
►
<*2—
1
1
II
1耳-
1
—3P*-
P1—*
1
1
1
1
1
1
X
►
0I
2.基本方程组、初始条件和边界条件
设流体是黏性流体。
一维可压缩黏性流动问题，在数学上可以用
T
E
CT
12
II
Pr
Cp
p2,
M2
E
Cv1
u
2
Pr
k
UL
1
Cp
T
Re,
cv
9j
a
1
M2
Cv
M
为:
f—（B.1）
xx
其中
（B.1b）
其中、up和E分别是量纲为一的密度、速度、压力和单位体积总能，s为流体的黏性项。
Pr为普朗特数（此处公式中Cp、k、是有量纲量），Re为雷诺数，c为比定容热容，Cp为比定压热容，Cp,cv是量纲为一的量，称为气体绝热指数，a为当地声速。
初始条件：
在t0时刻，x,01，其他物理量采用线性插值得
到。
右边界x1处：
x,t
0
x1
x
u1,t
2
M2
1
1..2—M
1
T1,t2M2
1
11
2
11
1M2
3.二阶精度迎风型TVD差分格式
一维N-S方程组中的对流项采用
Harten的二阶精度迎风型
分格式：
n1n.n
ujujrh1
JJj2
hn
J
s
1
2X
t（B.4）
1
'ji2fjfj1Rj2°j2（B.4a）
其中r。
向量①在第l个特征方向上分量为:
x
iii
.1gjgj1Q
j2
l
a1
j2
l
.1.1
jTjT
22
（B.4b）
gjminmod11
J2
l
1,j2
ll
11
j-j-
22
（B.4c）
（B.3）
TVD差
ll
gj1gj
l
.1j
2
0
（B.4d）
0,
流通量矢量f的非线性Jiacobia系数矩阵A=Rar-1为:
非线性Jiacobia系数矩阵A的特征值为:
（B.6）
非线性Jiacobia系数矩阵A的右特征矢量R、R-1为:
一维N-S方程组中的黏性项
1
1
1
Ru
uau
a
（B.7）
12u
2
huah
ua
1
1u2
1
u
1
2a2
2a
2a
R1
u
1u2
1u
1
1
2a2（B.8）
2a
4a2
2a2
2a
u
1u2
1u
1
1
2a
4a2
2a2
2a
2a2
t采用二阶精度中心差分格式
4.计算结果分析
采用用C语言和Fortran77语言对一维可压缩黏性流动问题编制了计算程序,并对雷诺数Re50的流动进行了计算，计算结果如图B.2和图B.3所示。
图B.2Fortran7语言程序得到图B.3C语言程序得到
图B.2和图B.3是计算达到稳定后激波间断位置和密度（）、速度（u）、压
的结果的结果
力（P）和单位质量内能（e）的分布。
由上述计算结果中可以看出，采用二阶精度迎风型TVD差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的结果是完全一致的。
计算结果表明，迎风型TVD差分格式能够精确地捕捉激波间断，计算效果较好。
由于本问题中黏性较大（Re50），所以计算得到的激波比较
光滑，有一定的宽度。
一维可压缩黏性流动问题的解是连续、光滑的
B-2一维可压缩黏性流动问题的数值计算源程序
1.C语言源程序
//UpwindTVD_1D.c
//
//二阶迎风型TVD差分格式求解一维可压缩黏性流动问题（C语言版本）
//
#include"stdio.h"
#include"math.h"
#definegama1.4
#defineTt5.0
#defineim201//网格数
//全局变量：
doubleQ[3][im],Qold[3][im];//Q:
[rou,rou*u,E]doublerou[im],u[im],p[im],T[im],E[im],a[im];doublePr,Re,cv,cp,Ma,dx,dt;
//voidinitial（）
{
doublexl,xr,x;
doubleul,Tl,ur,Tr;//进出口的u，T值inti;
dx=1.0/（im-1）;
dt=1.0e-6;
Pr=0.72;
Re=50.0;
Ma=2.0;
cv=1.0/（gama*（gama-1.0）*Ma*Ma）;
cp=gama*cv;
xl=0.0;
xr=1.0;
ul=1.0;
Tl=1.0;
ur=（2.0/（gama-1.0）+Ma*Ma）/（（gama+1.0）/（gama-1.0）*Ma*
Ma）;
Tr=2.0*gama/（gama+1.0）*Ma*Ma-（gama-1.0）/（gama+1.0）;
Tr=Tr*（（gama-1.0）/（gama+1.0）+2.0/（gama+1.0）/Ma/Ma）;
for（i=0;i<=im-1;i++）
{
x=i*dx;
rou[i]=1.0;
u[i]=（x-xl）/（xr-xl）*（ur-ul）+ul;
T[i]=（x-xl）/（xr-xl）*（Tr-Tl）+Tl;
p[i]=rou[i]*T[i]/（gama*Ma*Ma）;
E[i]=rou[i]*（cv*T[i]+0.5*u[i]*u[i]）;
}
for（i=0;i{
Q[0][i]=rou[i];
Q[1][i]=rou[i]*u[i];
Q[2][i]=E[i];
}
}
//doubleqz（doublex）
{
doubleresult;
if（fabs（x）>0.1）
{
result=fabs（x）;
}
else
{
result=（x*x/0.1+0.1）/2.0;
}
returnresult;
}
//
doubleminmod（doublex,doubley）
{
doubleresult;
if（x*y<=0.0）
{
result=0.0;
}
elseif（fabs（x）>fabs（y））
{
result=y;
}
elseif（fabs（x）{
result=x;
}
returnresult;
}
voidQ2U（）
{
inti;
for（i=0;i{
rou[i]=Q[0][i];
u[i]=Q[1][i]/rou[i];
E[i]=Q[2][i];
T[i]=（E[i]/rou[i]-0.5*u[i]*u[i]）/cv;
p[i]=rou[i]*T[i]/（gama*Ma*Ma）;a[i]=sqrt（T[i]）/Ma;
}
}
//voidUpwindTVD_1D_Solver（）
{
inti,l,k;
doubledq[3][im],lamda[3][im],sigma[3][im],f[3][im];doublefl[3],fwave[3][im],fr[3],g[3][im],Gv[3][im];doublealfa[3][im],beta[3][im],Rn[3][3][im],R[3][3][im];
doubleut,at;//临时变量
Q2U（）;
for（i=0;if[0][i]=rou[i]*u[i];
f[1][i]=rou[i]*u[i]*u[i]+p[i];
f[2][i]=u[i]*（E[i]+p[i]）;
}
for（i=0;i<=im-2;i++）
{
lamda[0][i]=（u[i]+u[i+1]）/2.0;
lamda[1][i]=（u[i]-a[i]+u[i+1]-a[i+1]）/2.0;
lamda[2][i]=（u[i]+a[i]+u[i+1]+a[i+1]）/2.0;
}
for（i=0;i<=im-2;i++）
for（l=0;l<=2;l++）
{
dq[l][i]=Q[l][i+1]-Q[l][i];
sigma[l][i]=qz（lamda[l][i]）/2.0;
}
//计算左右特征向量
for（i=0;i<=im-2;i++）
{
ut=（u[i]+u[i+1]）/2.0;
at=（a[i]+a[i+1]）/2.0;
R[0][0][i]=1.0;
R[0][1][i]=1.0;
R[0][2][i]=1.0;
R[1][0][i]=ut;
R[1][1][i]=ut-at;
R[1][2][i]=ut+at;
R[2][0][i]=ut*ut/2.0;
R[2][1][i]=at*at/（gama-1.0）+ut*ut/2.0-ut*at;
R[2][2][i]=at*at/（gama-1.0）+ut*ut/2.0+ut*at;
Rn[0][0][i]=1.0-（gama-1.0）*ut*ut/2.0/at/at;
Rn[0][1][i]=（gama-1.0）*ut/at/at;
Rn[0][2][i]=-（gama-1.0）/at/at;
Rn[1][0][i]=ut/2.0/at+（gama-1.0）*ut*ut/4.0/at/at;
Rn[1][1][i]=-（gama-1.0）*ut/2./at/at-1.0/2.0/at;
Rn[1][2][i]=（gama-1.0）/2.0/at/at;
Rn[2][0][i]=-ut/2.0/at+（gama-1.0）*ut*ut/4.0/at/at;
Rn[2][1][i]=（1.0-gama）*ut/2./at/at+1.0/2.0/at;
Rn[2][2][i]=（gama-1.0）/2./at/at;
}
for（i=0;i<=im-1;i++）
for（l=0;l<=2;l++）
{
alfa[l][i]=0.0;
//计算alfa
for（i=0;i<=im-2;i++）
for（l=0;l<=2;l++）
for（k=0;k<=2;k++）
{alfa[l][i]=alfa[l][i]+Rn[l][k][i]*dq[k][i];
}
//计算g
for（i=1;i<=im-2;i++）
for（l=0;l<=2;l++）
{
g[l][i]=minmod（sigma[l][i]*alfa[l][i],sigma[l][i-1]*
alfa[l][i-1]）;
}
for（i=0;i<=im-2;i++）
for（l=0;l<=2;l++）
{
if（alfa[l][i]!
=0.0）
}
else
{
beta[l][i]=0.0;
}
}
for（i=0;i<=im-2;i++）
for（l=0;l<=2;l++）
for（k=0;k<=2;k++）
{
fwave[l][i]=fwave[l][i]+0.5*R[l][k][i]*（g[k][i]+g[k]
[i+1]
-qz（lamda[k][i]+beta[k][i]）*alfa[k][i]）;
}
for（i=1;i<=im-2;i++）
for（l=0;l<=2;l++）
{fl[l]=0.5*（f[l][i-1]+f[l][i]）+fwave[l][i-1];fr[l]=0.5*（f[l][i+1]+f[l][i]）+fwave[l][i];
}
//黏性项的处理
for（i=1;i<=im-2;i++）
{
Gv[0][i]=0.0;
Gv[1][i]=（T[i]*（u[i+1]-2.*u[i]+u[i-1]）+（T[i+1]-T[i-1]）
*（u[i+1]-u[i-1]）/4.）*4./3./Re/（dx*dx）;
Gv[2][i]=（（T[i]*（u[i+1]-2.*u[i]+u[i-1]）*u[i]+（T[i+1]-T[i-1]）*
（u[i+1]-u[i-1]）*u[i]/4.+T[i]*（u[i+1]-u[i-1]）
*（u[i+1]-u[i-1]）/4.）*4./3.+（T[i]*（T[i+1]-2.*T[i]+T[i-1]）
+（T[i+1]-T[i-1]）*（T[i+1]-T[i-1]）/4.）*cp/Pr）/Re/（dx*dx）;
}for（i=1;i<=im-2;i++）for（l=0;l<=2;l++）
Q[l][i]=Q[l][i]+Gv[l][i]*dt;
}
//边界条件
Q[0][im-1]=Q[0][im-2];
Q[1][im-1]=Q[0][im-1]*u[im-1];
Q[2][im-1]=Q[0][im-1]*（cv*T[im-1]+0.5*u[im-1]*u[im-1]）;
}
//doubleerror（doubleQ1[3][im],doubleQ2[3][im]）
{
inti,l;
doubleerr,ddu;
err=1.0e-10;
for（i=0;ifor（l=0;l<3;l++）
{
if（Q2[l][i]!
=0.0）
{
ddu=fabs（（Q2[l][i]-Q1[l][i]）/Q2[l][i]）;
}
else
{
ddu=fabs（Q2[l][i]-Q1[l][i]）;
}
if（err}
returnerr;
}
//voidoutput（）
{
doublex,rou1,u1,E1,p1,T1;//临时变量
inti;
FILE*fp;
fp=fopen（"result.dat","w"）;
fprintf（fp,"variables=x,rou,u,p,e\n"）;
for（i=0;i{
x=i*dx;
rou1=Q[0][i];
u1=Q[1][i]/rou1;
E1=Q[2][i];
T1=（E1/rou1-0.5*u1*u1）/cv;
p1=rou1*T1/（gama*Ma*Ma）;
fprintf（fp,"%15.8e%15.8e%15.8e%15.8e%15.8e\n",
x,rou1,u1,p1,cv*T1）;
}
fclose（fp）;
}
//main（）
{
doublet,err;
intn,i,l;
initial（）;
t=0.0;
n=0;
err=1.0e-7;
while（t1.0e-8）
t=t+dt;
n=n+1;
for（i=0;ifor（l=0;l<3;l++）
{
Qold[l][i]=Q[l][i];
}
UpwindTVD_1D_Solver（）;
err=error（Q,Qold）;
if（n/10000*10000==n）
{
printf（"%10dtime:
%16.9eerror:
%16.9e\n",n,t,err）;
output（）;
}
}
output（）;
printf（"Programend!
\n"）;
}
2.Fortran77语言源程序
!
UpwindTVD_1D.f
!
!
二阶迎风型TVD差分格式求解一维可压缩黏性流动问题!
（Fortran77语言版本）
!
programUPWIND_TVD_1D
implicitreal*8（a-h,o-z）parameter（mx=201,Tt=5.0）
dimensionQ（3,mx）,Qold（3,mx）!
Q:
[rou,rou*u,E]dimensionrou（mx）,u（mx）,p（mx）,T（mx）,E（mx）
real*8Ma
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imim=mx
callInitialize（Q）
time=0.0
n=0
1n=n+1
time=time+dt
doi=1,im
dol=1,3
Qold（l,i）=Q（l,i）
enddo
enddo
callUpwindTVD_1D_Solver（Q）write（*,20）n,time,error（Q,Qold）
callOutPut（Q）
endif
20format（1x,i10,'time=',e16.9,'error=',e16.9,1x）.and.
goto1
endif
callOutPut（Q）
write（*,*）'Programend!
'
end
!
subroutineInitialize（Q）
implicitreal*8（a-h,o-z）
real*8Ma
common/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimensionrou（im）,u（im）,p（im）,T（im）,E（im）,Q（3,im）
Re=50.0!
雷诺数
Ma=2.0!
马赫数
pr=0.72!
普朗特数
gama=1.4!
气体常数
cv=1.0/（gama*（gama-1.0）*Ma**2）!
比定容热容
cp=gama*cv!
比定压热容
dt=1.0e-6
dx=1.0/（1.0*（im-1））!
空间步长
xl=0.0
xr=1.0
rou0=1.0
u0=1.0
T0=1.0
temp=（gama+1.）/（gama-1.）
u1=（2./（gama-1.）+Ma**2）/（temp*Ma**2）
T1=（2.*gama/（gama+1）*Ma**2-1./temp）*（1./temp+2./（gam
a+1.）/Ma**2）
doi=1,im
rou（i）=rou0
xi=（i-1）*dx
u（i）=（xi-xl）/（xr-xl）*（u1-u0）+u0
T（i）=（xi-xl）/（xr-xl）*（T1-T0）+T0
P（i）=rou（i）*T（i）/（gama*Ma**2）
E（i）=rou（i）*（cv*T（i）+0.5*u（i）**2）
enddo
doi=1,im
Q（1,i）=rou（i）
Q（2,i）=rou（i）*u（i）
Q（3,i）=E（i）
enddo
end
!
subroutineQ2U（Q,rou,u,p,T,E,a）implicitreal*8（a-h,o-z）real*8Macommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimensionQ（3,im）,rou（im）,u（im）dimensionp（im）,T（im）,E（im）,a（im）doi=1,imrou（i）=Q（1,i）u（i）=Q（2,i）/Q（1,i）
E（i）=Q（3,i）T（i）=（E（i）/rou（i）-0.5*u（i）**2）/cvp（i）=rou（i）*T（i）/（gama*Ma**2）a（i）=dsqrt（T（i））/Maenddo
end
!
subroutineUpwindTVD_1D_Solver（Q）
implicitreal*8（a-h,o-z）
real*8Ma,lamda,minmodcommon/para_def/Sf,Re,Ma,pr,gama,cp,cv,dx,dt,imdimensionQ（3,im）,Qold（3,im）,dq（3,im）,rou（im）,u（im）,p（im）,T（im）
dimensionE（im）,a（im）,f（3,im）,fwave（3,im）,fr（3）,fl（3）dimensionGv（3,im）,g（3,im）,lamda（3,im）,sigma（3,im）dimensionalfa（3,im）,Rn（3,3,im）,R（3,3,im）,beta（3,im）callQ2U（Q,rou,u,p,T,E,a）
doi=1,im
f（1,i）=rou（i）*u（i）
f（2,i）=rou（i）*u（i）**2+p（i）
f（3,i）=u（i）*（E（i）+p（i））
enddo
doi=1,im-1
lamda（1,i）=（u（i）+u（i+1））/2.0
lamda（2,i）=（u（i）-a（i）+u（i+1）-a（i+1））/2.0
lamda（3,i）=（u（i）+a（i）+u（i+1）+a（i+1））/2.0
enddo
doi=1,im-1
dol=1,3
dq（l,i）=（Q（l,i+1）-Q（l,i））sigma（l,i）=qz（lamda（l,i））/2.0!
定常情况enddo
enddo
doi=1,im-1
ut=（u（i）+u（i+1））/2.0!
ut临时变量
at=（a（i）+a（i+1））/2.0!
at临时变量
R（1,1,i）=1.0
R（1,2,i）=1.0
R（1,3,i）=1.0
R（2,1,i）=ut
R（2,2,i）=ut-at
R（2,3,i）=ut+at
R（3,1,i）=ut*ut/2.0
R（3,2,i）=at*at/（gama-1.0）+ut*ut/2.0-ut*atR（3,3,i）=at*at/（gama-1.0）+ut*ut/2.0+ut*atRn（1,1,i）=1.0-（gama-1.0）*ut*ut/2.0/at/atRn（1,2,i）=（gama-1.0）*ut/at/at
Rn（1,3,i）=-（gama-1.0）/at/at
Rn（2,1,i）=ut/2.0/at+（gama-1.0）*ut*ut/4.0/at/at
Rn（2,2,i）=-（gama-1.0）*ut/2./at/at-1.0/2.0/at
Rn（2,3,i）=（gam
升级会员 