论坛 › CFD实战技术 › [原创]Ausm格式的简单源代码

iceArcher 发表于 2003-9-30 15:44:03

[原创]Ausm格式的简单源代码

    自己编的很简单的CFD二维求解程序，原先是学习Ausm格式用的。说它简单，因为除了最必须的东东，其余一概没有，刚好拿来玩。况且这个结果的精度不高。至于计算结果的正确性如何，偶也8晓得啦。*^_^*
    大家有空看看，最好对程序提些意见。可别用鸡蛋砸偶呵。
//---------------------------------------------------------------------------
//filename: TEuler2d-mod.cpp
//written by : Dr.Tinghe Lee
//this is additional cpp file to assume the right of my code
//a very neat code for ausm scheme test usage!
//---------------------------------------------------------------------------
#include <fstream>
#include <iostream>
#include <cmath>
#define gama 1.4
using namespace std;
//---------------------------------------------------------------------------
// predeclare variables
int im,jm,imax,jmax,it,itmax;
double rm,cfl;
double *x,*y,*rr,*ru,*rv,*re,*pp,*flxrj,*flxuj,*flxvj,*flxej,*dtime;
double *vol,*rrn,*run,*rvn,*ren,*flxri,*flxui,*flxvi,*flxei;
double *rrt,*rut,*rvt,*ret,p0,r0,u0,v0,e0,c0,q0;
//---------------------------------------------------------------------------
// predeclare functions
void input(),initflow(),bound(),rnk(),resid(),output(),farbnd(),walbnd(),
   timestep(),tecgrid();
//---------------------------------------------------------------------------
template<class Type> Type mmin(Type a,Type b)
{
    return (a<=b)?a:b;
}
//---------------------------------------------------------------------------
template<class Type> Type mmax(Type a,Type b)
{
    return (a>=b)?a:b;
}
//---------------------------------------------------------------------------
int main()
{
    rm= 3.0;
    cfl = 1.1;
    itmax = 2000;
   
    input();
    tecgrid();
    initflow();
    bound();
    for(it=0;it<itmax;it++)
    {
      rnk();
      resid();
    }
    output();
    delete [] vol; delete [] pp; delete [] rr; delete [] ru; delete [] rv; delete [] re;
    delete [] rrn; delete [] run; delete [] rvn; delete [] ren; delete [] dtime;
    delete [] rrt; delete [] rut; delete [] rvt; delete [] ret;
    delete [] flxri; delete [] flxui; delete [] flxvi; delete [] flxei;
    delete [] flxrj; delete [] flxuj; delete [] flxvj; delete [] flxej;
    delete [] x; delete [] y;
    return 0;
}
//---------------------------------------------------------------------------
void input()
{
    int i,j,ii,i1,i2,i3,i4;
    double rx1,rx2,ry1,ry2,vot;
    fstream mfile;
    mfile.open("gri2d.dat",ios::in);
    mfile>> im >> jm;
    imax = im+2;
    jmax = jm+2;
    x = new double ;
    y = new double ;
    for(i=1;i<=im;i++)
    for(j=1;j<=jm;j++)
    {
      ii = i*jmax+j;
      mfile>> x>>y;
    }
    mfile.close();
    vol = new double ;
    rr= new double ;
    ru= new double ;
    rv= new double ;
    re= new double ;
    pp= new double ;
    flxri = new double ;
    flxui = new double ;
    flxvi = new double ;
    flxei = new double ;
    flxrj = new double ;
    flxuj = new double ;
    flxvj = new double ;
    flxej = new double ;
    rrn= new double ;
    run= new double ;
    rvn= new double ;
    ren= new double ;
    rrt= new double ;
    rut= new double ;
    rvt= new double ;
    ret= new double ;
    dtime = new double ;
   
    for(i=1;i<im+1;i++)
    for(j=1;j<jm+1;j++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      i3 = (i+1)*jmax+j+1;
      i4 = i*jmax+j+1;
      rx1 = x-x;
      ry1 = y-y;
      rx2 = x-x;
      ry2 = y-y;
      vot = rx1*ry2-rx2*ry1;
      rx2 = rx1;
      ry2 = ry1;
      rx1 = x-x;
      ry1 = y-y;
      vol = 0.5*(rx1*ry2-rx2*ry1+vot);
    }
}
//---------------------------------------------------------------------------
void initflow()
{
    int i,j,ii;
    p0 = 1;
    r0 = 1;
    c0 = sqrt(gama*p0/r0);
    q0 = rm*c0;
    u0 = q0;
    v0 = 0.0;
    e0 = p0/((gama-1.0)*r0)+q0*q0*0.5;
    for(i=0;i<imax;i++)
    for(j=0;j<jmax;j++)
    {
      ii = i*jmax+j;
      rr = r0;
      ru = r0*u0;
      rv = r0*v0;
      re = r0*e0;
      pp = p0;
    }
}
//---------------------------------------------------------------------------
void farbnd()
{
    int i,j,ii,i1;
    for(i=0;i<=0;i++)
    for(j=1;j<jm;j++)
    {
      ii = i*jmax+j;
      rr = r0;
      ru = r0*u0;
      rv = r0*v0;
      re = r0*e0;
      pp = p0;
    }
    for(i=im;i<=im;i++)
    for(j=1;j<jm;j++)
    {
      ii = i*jmax+j;
      i1 = (i-1)*jmax+j;
      rr = rr;
      ru = ru;
      rv = rv;
      re = re;
      pp = pp;
    }
}
//---------------------------------------------------------------------------
void walbnd()
{
    int i,j,ii,i1;
    for(i=1;i<im;i++)
    for(j=0;j<=0;j++)
    {
      ii = i*jmax+j;
      i1 = i*jmax+j+1;
      rr = 0;
      ru = -ru;
      rv = -rv;
      re = 0;
      pp = pp;
    }
    for(i=1;i<im;i++)
    for(j=jm;j<=jm;j++)
    {
      ii = i*jmax+j;
      i1 = i*jmax+j-1;
      rr = 0;
      ru = -ru;
      rv = -rv;
      re = 0;
      pp = pp;
    }
}
//---------------------------------------------------------------------------
void flux()
{
    int i,j,i1,i2,i3,i4,il,ir;
    double rrl,rul,rvl,rel,rhl,rrr,rur,rvr,rer,rhr,ppl,ppr;
    double xn,yn,al,ar,area,ql,qr,ml,mr,ms,pl,pr,alf,arf,ps,as;
   
    for(i=1;i<=im;i++)
    for(j=1;j<jm;j++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      i3 = (i+1)*jmax+j+1;
      i4 = i*jmax+j+1;
      il = (i-1)*jmax+j;
      ir = i1;
      xn = y-y;
      yn = x-x;
      area = sqrt(xn*xn+yn*yn);
      rrl = rr;
      rul = ru;
      rvl = rv;
      rel = re;
      ppl = pp;
      rhl = rel+ppl;
      rrr = rr;
      rur = ru;
      rvr = rv;
      rer = re;
      ppr = pp;
      rhr = rer+ppr;
      al = sqrt(2.0*(gama-1.0)*rhl/rrl/(gama+1.0));
      ar = sqrt(2.0*(gama-1.0)*rhr/rrr/(gama+1.0));
      ql = (rul*xn+rvl*yn)/rrl/area;
      qr = (rur*xn+rvr*yn)/rrr/area;
      al = al*al/mmax(fabs(ql),al);
      ar = ar*ar/mmax(fabs(qr),ar);
      as = mmin(al,ar);
      ml = ql/as;
      mr = qr/as;
      if(fabs(ml)<1.0) ml = 0.25*(ml+1.0)*(ml+1.0);
      if(ml<=-1.0) ml = 0.0;
      if(fabs(mr)<1.0) mr = -0.25*(mr-1.0)*(mr-1.0);
      if(mr>=1.0) mr = 0.0;
      pl = ppl;
      pr = ppr;
      if(fabs(ml)<1.0) pl = 0.25*ppl*(ml+1.0)*(ml+1.0)*(2.0-ml);
      if(ml<=-1.0) pl = 0.0;
      if(fabs(mr)<1.0) pr = 0.25*ppr*(mr-1.0)*(mr-1.0)*(2.0+mr);
      if(mr>=1.0) pr = 0.0;
      ps= pl+pr;
      ms= ml+mr;
      alf = 0.5*ms*area;
      arf = 0.5*(fabs(ms)+0.5*(mr-1.0)*(mr-1.0))*area;
      flxri = alf*(rrl*al+rrr*ar)-arf*(rrr*ar-rrl*al);
      flxui = alf*(rul*al+rur*ar)-arf*(rur*ar-rul*al)+ps*xn;
      flxvi = alf*(rvl*al+rvr*ar)-arf*(rvr*ar-rvl*al)+ps*yn;
      flxei = alf*(rhl*al+rhr*ar)-arf*(rhr*ar-rhl*al);
    }
    for(j=2;j<jm;j++)
    for(i=1;i<im;i++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      i3 = (i+1)*jmax+j+1;
      i4 = i*jmax+j+1;
      il = i*jmax+j-1;
      ir = i1;
      xn = y-y;
      yn = x-x;
      area = sqrt(xn*xn+yn*yn);
      rrl = rr;
      rul = ru;
      rvl = rv;
      rel = re;
      ppl = pp;
      rhl = rel+ppl;
      rrr = rr;
      rur = ru;
      rvr = rv;
      rer = re;
      ppr = pp;
      rhr = rer+ppr;
      al = sqrt(2.0*(gama-1.0)*rhl/rrl/(gama+1.0));
      ar = sqrt(2.0*(gama-1.0)*rhr/rrr/(gama+1.0));
      ql = (rul*xn+rvl*yn)/rrl/area;
      qr = (rur*xn+rvr*yn)/rrr/area;
      al = al*al/mmax(fabs(ql),al);
      ar = ar*ar/mmax(fabs(qr),ar);
      as = mmin(al,ar);
      ml = ql/as;
      mr = qr/as;
      if(fabs(ml)<1.0) ml = 0.25*(ml+1.0)*(ml+1.0);
      if(ml<=-1.0) ml = 0.0;
      if(fabs(mr)<1.0) mr = -0.25*(mr-1.0)*(mr-1.0);
      if(mr>=1.0) mr = 0.0;
      pl = ppl;
      pr = ppr;
      if(fabs(ml)<1.0) pl = 0.25*ppl*(ml+1.0)*(ml+1.0)*(2.0-ml);
      if(ml<=-1.0) pl = 0.0;
      if(fabs(mr)<1.0) pr = 0.25*ppr*(mr-1.0)*(mr-1.0)*(2.0+mr);
      if(mr>=1.0) pr = 0.0;
      ps= pl+pr;
      ms= ml+mr;
      alf = 0.5*ms*area;
      arf = 0.5*(fabs(ms)+0.5*(mr-1.0)*(mr-1.0))*area;
   
      flxrj = alf*(rrl*al+rrr*ar)-arf*(rrr*ar-rrl*al);
      flxuj = alf*(rul*al+rur*ar)-arf*(rur*ar-rul*al)+ps*xn;
      flxvj = alf*(rvl*al+rvr*ar)-arf*(rvr*ar-rvl*al)+ps*yn;
      flxej = alf*(rhl*al+rhr*ar)-arf*(rhr*ar-rhl*al);
    }
    for(j=1;j<=1;j++)
    for(i=1;i<im;i++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      xn = y-y;
      yn = x-x;
      ps= pp;
      flxrj =0;
      flxuj =ps*xn;
      flxvj =ps*yn;
      flxej =0;
    }
    for(j=jm;j<=jm;j++)
    for(i=1;i<im;i++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      i3 = i*jmax+j-1;
      xn = y-y;
      yn = x-x;
      ps= pp;
      flxrj =0;
      flxuj =ps*xn;
      flxvj =ps*yn;
      flxej =0;
    }
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      i4 = i*jmax+j+1;
      rrn = -flxri+flxri-flxrj+flxrj;
      run = -flxui+flxui-flxuj+flxuj;
      rvn = -flxvi+flxvi-flxvj+flxvj;
      ren = -flxei+flxei-flxej+flxej;
    }
}
//---------------------------------------------------------------------------
void rnk()
{
    int i,j,n,ii;
    double rq;
    double alf={0.25,1.0/6.0,3.0/8.0,0.5,1.0};
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
      ii = i*jmax+j;
      rrt = rr;
      rut = ru;
      rvt = rv;
      ret = re;
    }
    for(n=0;n<5;n++)
    {
      timestep();
      flux();
      for(i=1;i<im;i++)
      for(j=1;j<jm;j++)
      {
            ii = i*jmax+j;
            rr = rrt - alf*dtime*rrn/vol;
            ru = rut - alf*dtime*run/vol;
            rv = rvt - alf*dtime*rvn/vol;
            re = ret - alf*dtime*ren/vol;
            rq = ru*ru+rv*rv;
            pp = (gama-1.0)*(re-rq/rr*0.5);
            if(pp<0||rr<0)
            {
                rr = r0;
                ru = r0*u0;
                rv = r0*v0;
                re = r0*e0;
                pp = p0;
            }
      }
      bound();
    }
}
//---------------------------------------------------------------------------
void bound()
{
    farbnd();
    walbnd();
}
//---------------------------------------------------------------------------
void timestep()
{
    int i,j,i1,i2,i3,i4;
    double xt,yt,xn,yn,a1,a2,c,uu,vv;
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      i3 = (i+1)*jmax+j+1;
      i4 = i*jmax+j+1;
      c= sqrt(gama*pp/rr);
      xt = 0.5*(x-x+x-x);
      yt = 0.5*(y-y+y-y);
      xn = yt;
      yn = -xt;
      a1 = sqrt(xn*xn+yn*yn);
      uu = (ru*xn+rv*yn)/rr/a1;
      xt = 0.5*(x-x+x-x);
      yt = 0.5*(y-y+y-y);
      xn = -yt;
      yn =xt;
      a2 = sqrt(xn*xn+yn*yn);
      vv = (ru*xn+rv*yn)/rr/a2;
      dtime = cfl*vol/(fabs(uu)+fabs(vv)+c*(a1+a2));
    }
}
//---------------------------------------------------------------------------
void resid()
{
    int i,j,ii,ierrm,jerrm;
    double err,errmax,dskip;
    err = 0;
    errmax = 0;
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
      ii = i*jmax+j;
      dskip = rr-rrt;
      err = dskip*dskip+err;
      if(fabs(dskip)>errmax)
      {
            errmax = fabs(dskip);
            ierrm = i;
            jerrm = j;
      }
    }
    err = sqrt(err/(im-1)/(jm-1));
    cout<<"it="<<it<<" err="<<err<<" irm="<<ierrm<<" jrm="<<jerrm
      <<" errm="<<errmax<<endl;
}
//---------------------------------------------------------------------------
void output()
{
    int i,j,i1,i2,i3,i4;
    double uu,vv,qq,cc,ma;
    double *xc,*yc;
    xc = new double;
    yc = new double;
    for(i=1;i<im;i++)
    for(j=1;j<jm;j++)
    {
      i1 = i*jmax+j;
      i2 = (i+1)*jmax+j;
      i3 = (i+1)*jmax+j+1;
      i4 = i*jmax+j+1;
      xc = 0.25*(x+x+x+x);
      yc = 0.25*(y+y+y+y);
    }
    fstream mfile;
    mfile.open("tecress.dat",ios::out);
    mfile<<"variables=\"x\",\"y\",\"ma\",\"r\",\"p\",\"u\",\"v\"";
    mfile<<endl;
    mfile<<"zone i="<<im-1<<" j="<<jm-1<<" f=point"<<endl;
    for(j=1;j<jm;j++)
    for(i=1;i<im;i++)
    {
      i1 = i*jmax+j;
      uu = ru/rr;
      vv = rv/rr;
      qq = sqrt(uu*uu+vv*vv);
      cc = sqrt(gama*pp/rr);
      ma = qq/cc;
      mfile<<xc<<" "<<yc<<" "<<ma<<" "<<rr<<" ";
      mfile<<pp<<" "<<ru/rr<<" "<<rv/rr<<endl;
    }
    mfile.close();
    delete [] xc;
    delete [] yc;
}
//---------------------------------------------------------------------------
void tecgrid()
{
    int i,j,ii;
    fstream mfile;
    mfile.open("tecMe.dat",ios::out);
    mfile<<"variables=\"x\",\"y\""<<endl;
    mfile<<"zone i="<<im<<" j="<<jm<<" f=point"<<endl;
    for(j=1;j<=jm;j++)
    for(i=1;i<=im;i++)
    {
      ii = i*jmax+j;
      mfile<<x<<" "<<y<<endl;
    }
    mfile.close();
}
//---------------------------------------------------------------------------

boling2000 发表于 2003-10-2 00:15:57

[原创]Ausm格式的简单源代码

很不错呀，编程风格也很好，看来功底很深厚！！！
佩服中。。。

tigerwoods 发表于 2003-10-6 18:06:38

[原创]Ausm格式的简单源代码

不错，不错。
虽然激波形状很粗，但网格也稀。

solomonzj 发表于 2003-10-31 10:55:45

[原创]Ausm格式的简单源代码

大侠能不能对 Ausm 格式，进行一点介绍亚！
渴望了解！有什么优点吗？
谢谢了

songke21 发表于 2004-3-15 18:27:40

[原创]Ausm格式的简单源代码

楼主，我也在搞AUSM格式，不少地方弄不明白。请问你的格式里用的是什么限制器？我想用Van Albada限制器，但是其具体形式看不太明白，请指教。

wangdingxi 发表于 2004-4-12 08:28:49

[原创]Ausm格式的简单源代码

最近我也了解了一下AUSM格式，AUSM是Advection Upstream Splitting Method,现在AUSM已经发展成一系列格式，如AUSM+,AUSMD,AUSMV,AUSMDV等等。在原理上AUSM格式和Van Leer flux splitting差不多，都是属于矢通量分裂。在AUSM中，将压力进行了分裂，但是Van Leer flux splitting没有。详细可以看看这个链接http://www.chimeracfd.com/java/gryphon/fluxausm.html

laser2k 发表于 2006-1-6 17:01:55

[原创]Ausm格式的简单源代码

[这个贴子最后由laser2k在 2006/01/06 05:04pm 第 1 次编辑]

   

zhvickie 发表于 2006-1-10 13:46:17

[原创]Ausm格式的简单源代码

有没有例子数据？

oldcat 发表于 2006-1-31 12:29:32

[原创]Ausm格式的简单源代码

楼主能把程序拿出来共享，应该鼓励。
我也觉得例子没有算好，所以看了一下程序，发现这个程序只能用来算正交网格，lz给出的例子不恰当，因为网格不是正交的。

gotowuhan 发表于 2006-4-5 08:22:26

[原创]Ausm格式的简单源代码

鼓励原创

applytina 发表于 2013-7-11 21:06:51

回复 1# iceArcher 的帖子

这个程序是不是没贴完，还是要从别的地方看。正在学习ausm格式。对其中一些点很是不理解。

kaiqun98 发表于 2021-5-1 22:18:11

本帖最后由 kaiqun98 于 2021-5-1 22:19 编辑

感谢楼主分享

查看完整版本: [原创]Ausm格式的简单源代码