ME 602 Computational Fluid Dynamics Assignment 02
Submitted By:
#include #include #include #include using namespace std;
enum {imax=129,jmax=129,Re=100,itr_si_max=10,itr_w_max =2,U=1}; //imax = max grid pt. starting frm 1 in x dir //jmax=max grid pt. starting frm 1 in y dir //U =non dimentional vel of lid (at the top wall) //Re=reynolds no fr the expt. //itr_si_ma=no of itr for si //itr_w_max=no of itr fr omega long double si[150][150], u[150][150], v[150][150], omeg[150][150], Er,dx, dy, omeg_prev, si_prev; int i, j,itr_si, itr_w, counter=0;
//calc the grid spacing void calc_grid_spacing() { dx = 1.0/(imax-1); dy = 1.0/(jmax-1); }
//setting the initial boundry conditions void initialize_BC() { for(i=1;i<=imax;i++) { si[i][1] = 0; si[i][jmax] = 0; u[i][1] = 0; u[i][jmax] = U; v[i][1] = 0; v[i][jmax] = 0; omeg[i][1] = 0; omeg[i][jmax] = -2*U/dy; }
for(j=1;j<=jmax;j++) { si[1][j] = 0; si[imax][j] = 0; u[1][j] = 0; v[1][j] = 0; u[imax][j] = 0; v[imax][j] = 0; omeg[1][j] = 0; omeg[imax][j] = 0; } }
//initializing parameters in the soln domain void init_sol_domain() { for(j=2;j1;j--) { for(i=imax-1;i>1;i--) { si_prev = si[i][j]; si[i][j] = (si[i+1][j] + si[i-1][j] + si[i][j+1] +si[i][j-1])/4 + (omeg[i][j]*dx*dy)/4;
Er += fabs(si[i][j]) - fabs(si_prev); } } } //iterating for the value of u & v for(j=jmax-1;j>1;j--) { for(i=imax-1;i>1;i--) { u[i][j] = (si[i][j+1] - si[i][j-1])/(2*dy); v[i][j] = (si[i-1][j] - si[i+1][j])/(2*dx); } }
//iteration for the value of omega over the domain for(itr_w=0;itr_w1;j--) { for(i=imax-1;i>1;i--) { omeg_prev = omeg[i][j]; omeg[i][j] = (omeg[i+1][j]*(1-(u[i][j]*dx*Re)/2) + omeg[i-1][j]*(1+(u[i][j]*dx*Re)/2) +\ omeg[i][j+1]*(1-(v[i][j]*dx*Re)/2) + omeg[i][j-1]*(1+(v[i][j]*dx*Re/2)))/4; } } }
//updating the boundry cond. over the walls for (i=2;i=Ermax || counter<200);
}
void output(void) { FILE *fp_u, *fp_tab, *fp_v, *fp_si;
fp_u = fopen("u.txt", "w"); fp_v = fopen("v.txt", "w"); fp_si = fopen("si.txt", "w"); fp_tab = fopen("table.txt", "w"); printf("no of itr=%d\n", counter); //printing the required values to file //for printing u along vertical line through geometric centre of the cavity fprintf(fp_u,"y\t\tu\n\n"); for(j=1;j<=jmax;j++) { fprintf(fp_u,"%lf\t%lf\n",(j-1)*dy,u[65][j]); }
//for v along the horizontal line through geometric centre of the cavity fprintf(fp_v,"x\t\tv\n\n"); for(j=1;j<=jmax;j++) { fprintf(fp_v,"%lf\t%lf\n",(j-1)*dy,v[j][65]); } //for si contour fprintf(fp_si,"i\t\tj\t\tpsi\n\n"); for(j=1;j<=jmax;j++) { for(i=1;i<=imax;i++) { fprintf(fp_si,"%lf\t%lf\t%lf\n",(i-1)*dy,(j-1)*dy,si[i][j]); } }
fclose(fp_u); fclose(fp_v); fclose(fp_si);
//printing the output for comparison with ghia et al. fprintf(fp_tab," u-velocity\n\n"); fprintf(fp_tab,"129\t%lf\t\n",u[65][129]); fprintf(fp_tab,"126\t%lf\t\n",u[65][126]); fprintf(fp_tab,"125\t%lf\t\n",u[65][125]); fprintf(fp_tab,"124\t%lf\t\n",u[65][124]); fprintf(fp_tab,"123\t%lf\t\n",u[65][123]); fprintf(fp_tab,"110\t%lf\t\n",u[65][110]); fprintf(fp_tab,"95\t%lf\t\n",u[65][95]); fprintf(fp_tab,"80\t%lf\t\n",u[65][80]); fprintf(fp_tab,"65\t%lf\t\n",u[65][65]); fprintf(fp_tab,"59\t%lf\t\n",u[65][59]); fprintf(fp_tab,"37\t%lf\t\n",u[65][37]); fprintf(fp_tab,"23\t%lf\t\n",u[65][23]); fprintf(fp_tab,"14\t%lf\t\n",u[65][14]); fprintf(fp_tab,"10\t%lf\t\n",u[65][10]); fprintf(fp_tab,"9\t%lf\t\n",u[65][9]); fprintf(fp_tab,"8\t%lf\t\n",u[65][8]); fclose(fp_tab); } void main() { calc_grid_spacing(); initialize_BC(); init_sol_domain(); stream(); output(); system("PAUSE"); }
