b一维可压ns方程b文档格式.docx
《b一维可压ns方程b文档格式.docx》由会员分享，可在线阅读，更多相关《b一维可压ns方程b文档格式.docx（14页珍藏版）》请在冰豆网上搜索。
为普朗特数（此处公式中是有量纲量），为雷诺数，为比定容热容，为比定压热容，是量纲为一的量，称为气体绝热指数，为当地声速。
求解区域为。
取。
初始条件：
在时刻，，其他物理量采用线性插值得到。
边界条件：
左边界处：
右边界处：
3.二阶精度迎风型差分格式
一维方程组中的对流项采用的二阶精度迎风型差分格式：
（B.4a）
其中。
向量在第个特征方向上分量为：
（B.4c）
（B.4f）
（B.4g）
流通量矢量的非线性系数矩阵为：
非线性系数矩阵的特征值为：
非线性系数矩阵的右特征矢量为：
一维方程组中的黏性项采用二阶精度中心差分格式。
4.计算结果分析
采用用语言和语言对一维可压缩黏性流动问题编制了计算程序，并对雷诺数的流动进行了计算，计算结果如图和图所示。
图和图是计算达到稳定后激波间断位置和密度、速度、压力和单位质量内能的分布。
由上述计算结果中可以看出，采用二阶精度迎风型差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的结果是完全一致的。
计算结果表明，迎风型差分格式能够精确地捕捉激波间断，计算效果较好。
由于本问题中黏性较大，所以计算得到的激波比较光滑，有一定的宽度。
一维可压缩黏性流动问题的解是连续、光滑的。
B-2一维可压缩黏性流动问题的数值计算源程序
1.语言源程序
urat/at;
Rn[1][2][i]=/at/at;
Rn[2][0][i]=-ut/at+*ut*ut/at/at;
Rn[2][1][i]=*ut/2./at/at+at;
Rn[2][2][i]=/2./at/at;
}
for（i=0;
i<
=im-1;
i++）
for（l=0;
l<
=2;
l++）
{
alfa[l][i]=;
fwave[l][i]=;
g[l][i]=;
}
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;
=im-2;
for（l=0;
{
Q[l][i]=Q[l][i]+Gv[l][i]*dt;
言源程序
!
-----------------------------------------------------------------------
二阶迎风型差分格式求解一维可压缩黏性流动问题
（语言版本）
programUPWIND_TVD_1D
implicitreal*8（a-h,o-z）
parameter（mx=201,Tt=
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,im
im=mx
callInitialize（Q）
time=
n=0
1n=n+1
time=time+dt
doi=1,im
dol=1,3
Qold（l,i）=Q（l,i）
enddo
callUpwindTVD_1D_Solver（Q）
write（*,20）n,time,error（Q,Qold）
callOutPut（Q）
endif
20format（1x,i10,'
time='
,'
error='
,1x）
.and.
goto1
write（*,*）'
Programend!
'
end
--------------------------------------------------------------------
subroutineInitialize（Q）
dimensionrou（im）,u（im）,p（im）,T（im）,E（im）,Q（3,im）
Re=!
雷诺数
Ma=!
马赫数
pr=!
普朗特数
gama=!
气体常数
cv=（gama**Ma**2）!
比定容热容
cp=gama*cv!
比定压热容
dt=
dx=*（im-1））!
空间步长
xl=
xr=
rou0=
u0=
T0=
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./（gama+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）+*u（i）**2）
Q（1,i）=rou（i）
Q（2,i）=rou（i）*u（i）
Q（3,i）=E（i）
subroutineQ2U（Q,rou,u,p,T,E,a）
dimensionQ（3,im）,rou（im）,u（im）
dimensionp（im）,T（im）,E（im）,a（im）
rou（i）=Q（1,i）
u（i）=Q（2,i）/Q（1,i）
E（i）=Q（3,i）
T（i）=（E（i）/rou（i）*u（i）**2）/cv
p（i）=rou（i）*T（i）/（gama*Ma**2）
a（i）=dsqrt（T（i））/Ma
subroutineUpwindTVD_1D_Solver（Q）
real*8Ma,lamda,minmod
dimensionQ（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）
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））
doi=1,im-1
lamda（1,i）=（u（i）+u（i+1））/
lamda（2,i）=（u（i）-a（i）+u（i+1）-a（i+1））/
lamda（3,i）=（u（i）+a（i）+u（i+1）+a（i+1））/
dq（l,i）=（Q（l,i+1）-Q（l,i））
sigma（l,i）=qz（lamda（l,i））/!
定常情况
doi=1,im-1
ut=（u（i）+u（i+1））/!
ut临时变量
at=（a（i）+a（i+1））/!
at临时变量
R（1,1,i）=
R（1,2,i）=
R（1,3,i）=
R（2,1,i）=ut
R（2,2,i）=ut-at
R（2,3,i）=ut+at
R（3,1,i）=ut*ut/
R（3,2,i）=at*at/+ut*ut/*at
R（3,3,i）=at*at/+ut*ut/+ut*at
Rn（1,1,i）=*ut*ut/at/at
Rn（1,2,i）=*ut/at/at
Rn（1,3,i）=-/at/at
Rn（2,1,i）=ut/at+*ut*ut/at/at
Rn（2,2,i）=-*ut/2./at/at
Rn（2,3,i）=/at/at
Rn（3,1,i）=-ut/at+*ut*ut/at/at
Rn（3,2,i）=*ut/2./at/at+at
Rn（3,3,i）=/2./at/at
enddo
alfa（l,i）=
fwave（l,i）=
g（l,i）=
dok=1,3
alfa（l,i）=alfa（l,i）+Rn（l,k,i）*dq（k,i）
升级会员 