logo资料库

计算流体力学c程序及fortran 程序.doc

发布时间：2022-06-09 发布人：admin 分类：说明书 资料大小：0.32M 资料格式：doc 举报 版权申诉
第1页 / 共17页
第2页 / 共17页
第3页 / 共17页
第4页 / 共17页
第5页 / 共17页
第6页 / 共17页
第7页 / 共17页
第8页 / 共17页
资料共17页，剩余部分请下载后查看
    B-1 利用二阶迎风型差分格式求解一维可压缩黏
    （B.6）
    4.计算结果分析
    B-2 一维可压缩黏性流动问题的数值计算源程序
    Q2U() ;
    ! UpwindTVD_1D.f
    !（语言版本）
    Rn(1,1,i)=1.0-(gama-1.0)*ut*ut/2.0/at/at
    Rn(2,2,i)=-(gama-1.0)*ut/2./at/at-1.0/2.0/at
    Rn(3,1,i)=-ut/2.0/at+(gama-1.0)*ut*ut/4.0/at/at
    附录 B 一维可压缩黏性流动问题的数值解法与计算程序 一维可压缩黏性流动是气体动力学中最经典的黏性流动问题，对它采用迎风 型TVD差分算法进行数值求解。同时，为了初学者入门和练习方便,这里给出了 由 C 语言和 Fortran77语言编写的、计算一维可压缩黏性流动问题的计算程序，供 大家学习参考。 B-1 利用二阶迎风型TVD差分格式求解一维可压缩黏性流动问题 1.一维可压缩黏性流动问题 在两端开口管道中充满了可压缩黏性流体。当黏性流体以超声速从左向右运 动时，一定会在管道中形成一道正激波，如图 B.1 所示。 1 p 2 2 u p、 、 和 2 1 1 、 、 分 u 别为激波波前和波后的参数。该问题可简化为一维可压缩黏性流动问题。当数值 解达到稳定时，在管道中可求解得到一道稳定的激波。 图 B.1 一维可压缩黏性流动问题示意图 2.基本方程组、初始条件和边界条件 设流体是黏性流体。一维可压缩黏性流动问题，在数学上可以用一维可压缩 黏性流动 N-S 方程组来描述。量纲为一的一维 N-S 方程组为： 其中 u  t   f  x   s  x  u     u    E       , f  u  2 u p     E p u            , s       s 1 s 2 s 3                 -B.1- （B.1） （B.1a）             T x 0 4 u   3 Re x  4 u u   3 Re  C K x Pr Rr  p
    p Re  ,  T  2 M   U L      E      C T v  1 2 2 u ,    Pr  , c v  1  1  M 2     ,   c c p v c  p  k  , a  T M  （B.1b） 其中 u p、、 和 E 分别是量纲为一的密度、速度、压力和单位体积总能，s 为流体 的黏性项。 Pr 为普朗特数（此处公式中 pc 、 、 是有量纲量）， Re 为雷诺数， k   vc 为比定容热容， pc 为比定压热容， ,p c c 是量纲为一的量，称为气体绝热指数， v a 为当地声速。 求解区域为 0 x  。取 1 M   2, Re  50, Pr  0.72,   1.4, k    。 T 初始条件：在 0 t  时刻，  x ,0  1  ，其他物理量采用线性插值得到。 边界条件：左边界 x  处：   0 0, t   u  0, t   T  0, t   1 （B.2） 右边界 1x  处：      , x t x  0 x 1  u  1, t   T  1, t      2  1   M 2   M 2   2 1   1   1    1      1     2  （B.3） M     1 1  2  1    M 2      3.二阶精度迎风型TVD差分格式 一维 N-S 方程组中的对流项采用 Harten 的二阶精度迎风型TVD差分格式： 1  u n j  u n j  r    h n j  1 2  h n j  1 2  h j  1 2 1 2    f j  f j 1   R Φ j  1 2 s x          1 2  j  t （B.4） （B.4a） 其中 r  t  x  。向量Φ 在第l 个特征方向上分量 l 为：  g l j  g l j 1   l  j  1 2  Q a   l j  1 2  -B.2-      j l j  1 2  1 2 （B.4b）
    g l j  l min mod     j l j l j l j ,    1 2 1 2     1 2 1 2    g          1 2 0,  l  j  1 2 l  j  0 l  j  1 2  0 l  j  1 2   u l j  1 2 l 1 j  l  j  g l j ,  1 2 l j  R 1 z z  1 2   ,     Q z  z 1 2      ,    2 z       , （B.4c） （B.4d） （B.4e） （B.4g） （B.5） （B.6） （B.7）  0.05    0.5  （B.4f） l  i  1 2   1   2      1   2  Q Q    l  i     l  i  1 2  l r  i      2 1 2 ,    , 1 2    非定常 定常 流通量矢量 f 的非线性Jiacobian系数矩阵 =A R Λ R 为： -1           A u     f u  0   2 2 u 3  1  3    u 0   3 u  2 2 1   2 u a   非线性Jiacobian系数矩阵 A 的特征值为： 0  1 0  2 0 0 u a          2  1 Λ  u , 2 a   1 3   2 u  2 0   0     3  3 ,   u a     1      u  非线性Jiacobian系数矩阵 A 的右特征矢量 -1、R R 为： R         1 u 1 2 2 u 1  u a 1  u a h ua h ua          -B.3-
    2  1 R            1   u 2 a u 2 a  一维 N-S 方程组中的黏性项   4.计算结果分析 t u  2 u 1 2  u   u   1   2 2 a   1   2 4 a   1   2 4 a   a   1   2 2 a  1   2 2 a  1          s 采用二阶精度中心差分格式。  x   2 a 1   2 2 a 1   2 2 a 1 2 a 1 2 a  u     2 u （B.8） 采用用 C 语言和 Fortran77语言对一维可压缩黏性流动问题编制了计算程序， 并对雷诺数 Re  的流动进行了计算，计算结果如图 B.2 和图 B.3 所示。 50 图 B.2 Fortran77 语言程序得到的结果 图 B.3 C 语言程序得到的结果 图 B.2 和图 B.3 是计算达到稳定后激波间断位置和密度 ( ) 、速度 ( )u 、压 力 ( )p 和单位质量内能 ( )e 的分布。由上述计算结果中可以看出，采用二阶精度迎 风型 TVD差分格式计算一维可压缩黏性流动问题得到的数值解和经典文献中的 结果是完全一致的。计算结果表明，迎风型TVD差分格式能够精确地捕捉激波间 断，计算效果较好。由于本问题中黏性较大 ( Re  ，所以计算得到的激波比较 50) 光滑，有一定的宽度。一维可压缩黏性流动问题的解是连续、光滑的。 -B.4-
    B-2 一维可压缩黏性流动问题的数值计算源程序 UpwindTVD_1D.c 1. C 语言源程序 // //------------------------------------------------------------------------ // 二阶迎风型TVD差分格式求解一维可压缩黏性流动问题（ C 语言版本） //------------------------------------------------------------------------- #include "stdio.h" #include "math.h" #define gama 1.4 #define Tt 5.0 #define im 201 //网格数 //全局变量： double Q[3][im],Qold[3][im] ; //Q: [rou, rou*u, E] double rou[im],u[im],p[im],T[im],E[im],a[im] ; double Pr,Re,cv,cp,Ma,dx,dt ; //------------------------------------------------------------------- void initial() { double xl,xr,x; double ul,Tl,ur,Tr;//进出口的 u，T 值 int i; 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) ; -B.5-
    E[i]=rou[i]*(cv*T[i]+0.5*u[i]*u[i]) ; } for(i=0; i0.1) { result=fabs(x) ; } else { result=(x*x/0.1+0.1)/2.0 ; } return result ; } //------------------------------------------------------------------- double minmod(double x,double y) { double result ; if(x*y<=0.0) { result=0.0 ; } else if(fabs(x)>fabs(y)) { result=y ; } else if(fabs(x)
    for(i=0; i
    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 ; fwave[l][i]=0.0 ; g[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) -B.8-
    分享到：
    赞 收藏

    相关推荐

    资料库

    课程资源

    共收录17145份资料，累计13个分类，关注成员有19位，主要包括：PHP,网络管理,网页制作,Java,.Net,数据库,3G/移动开发,C/C++,游戏开发,嵌入式,讲义,软件测试,专业指导

    热门标签

    PHP 网络管理 网页制作 Java .Net 数据库 3G/移动开发 C/C++ 游戏开发 嵌入式 讲义 软件测试 专业指导

    最新资料