CFD Online URL
[Sponsors]
Home > Forums > Main CFD Forum

second and third order advection term discritisation in C

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   August 3, 2011, 08:51
Default second and third order advection term discritisation in C
  #1
New Member
 
yosuke
Join Date: Aug 2011
Posts: 5
Rep Power: 0
yosuke1984 is on a distinguished road
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
fig 1: first upwind


http://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;
}

Last edited by yosuke1984; August 3, 2011 at 09:25.
yosuke1984 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 01:17.