一维黎曼问题数值解与计算程序.docx
- 文档编号:715067
- 上传时间:2022-10-12
- 格式:DOCX
- 页数:16
- 大小:87.22KB
一维黎曼问题数值解与计算程序.docx
《一维黎曼问题数值解与计算程序.docx》由会员分享,可在线阅读,更多相关《一维黎曼问题数值解与计算程序.docx(16页珍藏版)》请在冰豆网上搜索。
一维黎曼问题数值解与计算程序
一维问题数值解与计算程序
一维问题,即激波管问题,是一个典型的一维可压缩无黏气体动力学问题,并有解析解。
对它采用二阶精度两步差分格式进行数值求解。
同时,为了初学者入门和练习方便,这里给出了用语言和编写的计算一维问题的计算程序,供大家学习参考。
A-1利用两步差分格式求解一维问题
1.一维问题
图A.1激波管问题示意图
一维问题实际上就是激波管问题。
激波管是一根两端封闭、内部充满气体的直管,如图A.1所示。
在直管中由一薄膜将激波管隔开,在薄膜两侧充有均匀理想气体(可以是同一种气体,也可以是不同种气体),薄膜两侧气体的压力、密度不同。
当时,气体处于静止状态。
当时,薄膜瞬时突然破裂,气体从高压端冲向低压端,同时在管内形成激波、稀疏波和接触间断等复杂波系。
2.基本方程组、初始条件和边界条件
设气体是理想气体。
一维问题在数学上可以用一维可压缩无黏气体方程组来描述。
在直角坐标系下量纲为一的一维方程组为:
(A.1)
其中(A.2)
这里、、、分别是流体的密度、速度、压力和单位体积总能。
理想气体状态方程:
(A.3)
初始条件:
;。
边界条件:
和处为自由输出条件,,。
3.二阶精度差分格式
两步差分格式:
(A.4)
其中。
计算实践表明,两步差分格式不能抑制激波附近非物
理振荡。
因此在计算激波时,必须采用人工黏性滤波方法:
(A.5)
为了在激波附近人工黏性起作用,而在光滑区域人工黏性为零,需要引入一个与密度(或者压力)相关的开关函数:
(A.6)
由式(A.6)可以看出,在光滑区域,密度变化很缓,因此值也很小;而在激波附近密度变化很陡,值就很大。
带有开关函数的前置人工黏性滤波方法为:
(A.7)
其中参数往往需要通过实际试算来确定,也可采用线性近似方法得到:
(A.8)
由于声速不会超过3,所以取,在本计算中取。
4.计算结果分析
计算分别采用标准的语言和语言编写程序。
计算中网格数取,计算总时间为。
计算得到在时刻的密度、速度和压力分布如图A.2(语言计算结果)和图A.3(语言计算结果)所示。
采用两种不同语言编写程序所得到的计算结果完全吻合。
从图A.2和图A.3中可以发现,两步差分格式能很好地捕捉激波,计算得到的激波面很陡、很窄,计算激波精度是很高的。
采用带开关函数的前置人工滤波法能消除激波附近的非物理振荡,计算效果很好。
从图A.2和图A.3中可以看出通过激波后气体的密度、压力和速度都是增加的;在压力分布中存在第二个台阶,表明在这里存在一个接触间断,在接触间断两侧压力是有间断的,而密度和速度是相等的。
这个计算结果正确地反映了一维问题的物理特性,并被激波管实验所验证。
图A.2采用语言程序得到的一维问题密度、速度和压力分布
图A.3采用语言程序得到的一维问题密度、速度和压力分布
A-2一维问题数值计算源程序
1.语言源程序
//MacCormack1D.cpp:
定义控制台应用程序的入口点。
//
/*-----------------------------------------------------------------------------------------------------
*利用差分格式求解一维激波管问题(语言版本)
*-------------------------------------------------------------------------------------------------------*/
//#include"stdafx.h"
#include
#include
#include
#defineGAMA1.4//气体常数
#definePI3.141592654
#defineL2.0//计算区域
#defineTT0.4//总时间
#defineSf0.8//时间步长因子
#defineJ1000//网格数
//全局变量
doubleU[J+2][3],Uf[J+2][3],Ef[J+2][3];
/*-------------------------------------------------------
计算时间步长
入口:
U,当前物理量,dx,网格宽度;
返回:
时间步长。
---------------------------------------------------------*/
doubleCFL(doubleU[J+2][3],doubledx)
{
inti;
doublemaxvel,p,u,vel;
maxvel=1e-100;
for(i=1;i<=J;i++)
{
u=U[i][1]/U[i][0];
p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u);
vel=sqrt(GAMA*p/U[i][0])+fabs(u);
if(vel>maxvel)maxvel=vel;
}
returnSf*dx/maxvel;
}
/*-------------------------------------------------------
初始化
入口:
无;
出口:
U,已经给定的初始值,
dx,网格宽度。
---------------------------------------------------------*/
voidInit(doubleU[J+2][3],double&dx)
{
inti;
doublerou1=1.0,u1=0.0,p1=1.0;//初始条件
doublerou2=0.125,u2=0.0,p2=0.1;
dx=L/J;
for(i=0;i<=J/2;i++)
{
U[i][0]=rou1;
U[i][1]=rou1*u1;
U[i][2]=p1/(GAMA-1)+rou1*u1*u1/2;
}
for(i=J/2+1;i<=J+1;i++)
{
U[i][0]=rou2;
U[i][1]=rou2*u2;
U[i][2]=p2/(GAMA-1)+rou2*u2*u2/2;
}
}
/*-------------------------------------------------------
边界条件
入口:
dx,网格宽度;
出口:
U,已经给定的边界。
---------------------------------------------------------*/
voidbound(doubleU[J+2][3],doubledx)
{
intk;
//左边界
for(k=0;k<3;k++)U[0][k]=U[1][k];
//右边界
for(k=0;k<3;k++)U[J+1][k]=U[J][k];
}
/*-------------------------------------------------------
根据U计算E
入口:
U,当前U矢量;
出口:
E,计算得到的E矢量,
U、E的定义见Euler方程组。
---------------------------------------------------------*/
voidU2E(doubleU[3],doubleE[3])
{
doubleu,p;
u=U[1]/U[0];
p=(GAMA-1)*(U[2]-0.5*U[1]*U[1]/U[0]);
E[0]=U[1];
E[1]=U[0]*u*u+p;
E[2]=(U[2]+p)*u;
}
/*-------------------------------------------------------
一维差分格式求解器
入口:
U,上一时刻的U矢量,Uf、Ef,临时变量,
dx,网格宽度,dt,时间步长;
出口:
U,计算得到的当前时刻U矢量。
---------------------------------------------------------*/
voidMacCormack_1DSolver(doubleU[J+2][3],doubleUf[J+2][3],doubleEf[J+2][3],doubledx,doubledt)
{
inti,k;
doubler,nu,q;
r=dt/dx;
nu=0.25;
for(i=1;i<=J;i++)
{
q=fabs(fabs(U[i+1][0]-U[i][0])-fabs(U[i][0]-U[i-1][0]))
/(fabs(U[i+1][0]-U[i][0])+fabs(U[i][0]-U[i-1][0])+1e-100);//开关函数
for(k=0;k<3;k++)
Ef[i][k]=U[i][k]+0.5*nu*q*(U[i+1][k]-2*U[i][k]+U[i-1][k]);//人工黏性项
}
for(k=0;k<3;k++)
for(i=1;i<=J;i++)U[i][k]=Ef[i][k];
for(i=0;i<=J+1;i++)U2E(U[i],Ef[i]);
for(i=0;i<=J;i++)
for(k=0;k<3;k++)
Uf[i][k]=U[i][k]-r*(Ef[i+1][k]-Ef[i][k]);//U(n+1/2)(i+1/2)
for(i=0;i<=J;i++)U2E(Uf[i],Ef[i]);//E(n+1/2)(i+1/2)
for(i=1;i<=J;i++)
for(k=0;k<3;k++)
U[i][k]=0.5*(U[i][k]+Uf[i][k])-0.5*r*(Ef[i][k]-Ef[i-1][k]);//U(n+1)(i)
}
/*-------------------------------------------------------
输出结果,用数据格式画图
入口:
U,当前时刻U矢量,dx,网格宽度;
出口:
无。
---------------------------------------------------------*/
voidOutput(doubleU[J+2][3],doubledx)
{
inti;
FILE*fp;
doublerou,u,p;
fp=fopen("result.txt","w");
for(i=0;i<=J+1;i++)
{
rou=U[i][0];
u=U[i][1]/rou;
p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u);
fprintf(fp,"%20f%20.10e%20.10e%20.10e%20.10e\n",i*dx,rou,u,p,U[i][2]);
}
fclose(fp);
}
/*-------------------------------------------------------
主函数
入口:
无;
出口:
无。
-----------------------
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 黎曼 问题 数值 计算 程序