CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   Main CFD Forum (http://www.cfd-online.com/Forums/main/)
-   -   second and third order advection term discritisation in C (http://www.cfd-online.com/Forums/main/91203-second-third-order-advection-term-discritisation-c.html)

yosuke1984 August 3, 2011 07:51

second and third order advection term discritisation in C
 
Hi all

I'm testing discritisation schemes for advection term, but struggling. Some how, second upwind and third upwind method diverges. I am using identical grid arrangement and courant conditions for all, and others ( quick, first upwind, kk-scheme) worked.

My reference says, higher order scheme gives smaller numerical diffusion but from what I can see from the result, it's opposite. The discritised equations were worked out from the expression on wikipedia and confirmed the reference as well.

Can any body point me to the right direction.

Here is my code and the result.

http://www.geocities.jp/yosuke_31415...rst_upwind.gif
http://www.geocities.jp/yosuke_31415...rst_upwind.gif
fig 1: first upwind

http://www.geocities.jp/yosuke_31415...ond_upwind.gif
http://www.geocities.jp/yosuke_31415...ond_upwind.gifhttp://www.geocities.jp/yosuke_3141592/pict/second_upwind.gif
fig 2: second upwind

Best Regards

y
---------------------------------------------------------------------------------------------------------------
#include<stdio.h>
#include<math.h>
#define NODES 128
#define V 2.0;

void init( const double *dt, const double *dx, double *t, double *u, double *f );

void advec_diff( FILE *wfp, const double *dt, const double *dx, const double *nu, double *u, double *f );
double second_upwind( int *i, double *f, double *u, const double *dx );

int main(){

char filename[128];
double t=0;
const double t_end = 0.5;
const double dt=0.01;
const double dx=0.1;
const double nu=0.0;
double u[NODES] ;
double f[NODES] ;

init( &dt, &dx, &t, u, f);


printf("C=%lf\n", u[0]*dt/dx);


FILE *wfp;

for( t=0; t<t_end; t+=dt){
sprintf( filename, "result-%lf.data", t);
if( (wfp=fopen( filename, "w")) == NULL ){
printf("File write ERROR\n");
}
advec_diff( wfp, &dt, &dx, &nu, u, f );
}
fclose(wfp);
return 0;}

void init( const double *dt, const double *dx, double *t, double *u, double *f ){

int i;
*t = 0.0;
for( i=0; i<NODES; i++){
u[i] = V;
f[i] = 0.0;
}
f[NODES/2-1] = 1.0;
f[NODES/2-2] = 1.0;
f[NODES/2-3] = 1.0;
f[NODES/2-4] = 1.0;
f[NODES/2-5] = 1.0;
}

void advec_diff( FILE *wfp, const double *dt, const double *dx, const double *nu, double *u, double *f ){

int i;
double dmy_f[NODES];
for( i=2; i<NODES-2; i++){
dmy_f[i] = f[i] + *dt * ( -second_upwind( &i, f, u, &*dx ) + *nu*( f[i-1] - 2.0*f[i] + f[i+1]) /( *dx**dx ) );
fprintf( wfp, " %lf\n", dmy_f[i]);
}
for( i=2; i<NODES-2; i++){
f[i] = dmy_f[i];
}
}

double second_upwind( int *i, double *f, double *u, const double *dx ){

double tmp;
tmp = u[*i]*( f[*i-2]-4.0*f[*i-1]+4.0*f[*i+1]-f[*i+2] )/(4.0**dx);
- fabs( u[*i] )*( -f[*i-2]+4.0*f[*i-1]-6.0*f[*i]+4.0*f[*i+1]-f[*i+2] )/ (4.0**dx);
return tmp;
}


All times are GMT -4. The time now is 16:31.