#include<stdio.h> #include<math.h> /*------ Flow in a Square Channel with a particle in between ---------*/ #define NF 1 #define Re 40.0 #define u_inlet 0.1036333333 #define NX 801 //no. of nodes in the x-direction #define NY 801 //no. of nodes in the y-direction #define a 40.0 //particle diameter #define radius 20.0 //particle radius #define g 400.0 // x-coordinate of the center of the spherical particle #define f 400.0 // y-coordinate of the center of the spherical particle #define lmtheta 2.25 // angular separation of lagrangian marker #define Q 9 // no. of direction of velocity vectors (D2Q9 model implemented) #define T 100001 // no. of time steps #define tau 0.8109 // dimensionless relaxation time due to fluid–particle collisions #define visc (tau-0.5)/3.0 //viscosity #define rho0 1.0 // density initialization //#define Ht NY-1.0 //height in vertical direction #define pie 3.141592653589 //value of pie #define w0 4.0/9.0 // weighting coefficient #define w1 1.0/9.0 // weighting coefficients for i = 1,2,3,4 #define w2 1.0/36.0 // weighting coefficients for i = 5,6,7,8 //#define Ub 0.0 #define h 1 //double CFX; //int k; int cx [Q]={0,1,0,-1,0,1,-1,-1,1}; //discrete velocity vectors int cy [Q]={0,0,1,0,-1,1,1,-1,-1}; // discrete velocity vectors double t [Q]={w0,w1,w1,w1,w1,w2,w2,w2,w2}; int next_x [NX][Q]; int next_y [NY][Q]; double SumFX; double r; double FX [NX][NY]; double FY [NX][NY]; double f1 [NX][NY][Q]; //main probability distribution function double f2 [NX][NY][Q]; //temporary PDF double rho [NX][NY]; double ux [NX][NY]; double uy [NX][NY]; double U [NX][NY]; double V [NX][NY]; //double ux_eq [NX][NY]; //double uy_eq [NX][NY]; double bnode [NX][NY]; //double ux_exact [NY]; void tables(); void initialise(); void cal_obs(); void IBM(); void collision(); void streaming(); void data(int, int); void cal_Cd(int, int); //void exact(); int main() { int count = 0, k; tables(); initialise(); for(k=0; k<T; k++) { printf("k=%d\r\t",k); cal_obs(); IBM(); collision(); streaming(); /*if(k%50000==0) { data(count, k); }*/ if(k==100000) { cal_Cd(count, k); } } return 0; } void tables() { int x,y,i; for(x=0; x<NX; x++) { for(y=0; y<NY; y++) { for(i=0; i<Q; i++) { next_x[x][i] = x + cx[i]; if(next_x[x][i]<0) next_x [x][i] = NX-1; if(next_x[x][i]>(NX-1)) next_x[x][i] = 0; next_y[y][i] = y + cy[i]; if(next_y[y][i]<0) next_y [y][i] = NY-1; if(next_y[y][i]>(NY-1)) next_y[y][i] = 0; } } } return ; } void initialise() { int x,y,i; for(x=0; x<NX; x++) { for(y=0; y<NY; y++) { rho[x][y] = rho0; ux[x][y] = 0.0; uy[x][y] = 0.0; for(i=0; i<Q; i++) { f1[x][y][i] = rho0*t[i]; } } } for(x=0; x<NX; x++) { for(y=0; y<NY; y++) { bnode[x][y] = 0.0; if(y==0 || y==(NY-1)) { bnode[x][y] = 1.0; } } } return ; } double cal_dela(int x, int y, double bX, double bY) { double rX = x - bX; //x-coordinate difference between the current marker and the current node double rY = y - bY; //y-coordinate difference between the current marker and the current node double arX = fabsf(rX); //absolute value of 'r' double arY = fabsf(rY); // ,, double delX = 0, delY = 0; //calculating delta function if(arX >= 0 && arX < 1) delX = 1/8.0*(3-(2*arX)+sqrt(1+(4*arX)-(4*arX*arX))); if(arX >= 1 && arX < 2) delX = 1/8.0*(5-(2*arX)-sqrt(-7+(12*arX)-(4*arX*arX))); else // |r| >= 2 delX = 0; if(arY >= 0 && arY < 1) //calculating delta function (y); delY = 1/8.0*(3-(2*arY)+sqrt(1+(4*arY)-(4*arY*arY))); if(arY >= 1 && arY < 2) delY = 1/8.0*(5-(2*arY)-sqrt(-7+(12*arY)-(4*arY*arY))); else // |r| >= 2 delY = 0; double D = delX * delY; //calculation of Eq. No. 26 return D; } double delta2(int x, int y, double bX, double bY) { double D; double rX = x - bX; //x-coordinate difference between the current marker and the current node double rY = y - bY; //y-coordinate difference between the current marker and the current node double arX = fabs(rX); //absolute value of 'r' double arY = fabs(rY); double delX = 0.0; double delY = 0.0; if(arX <= 1) delX = 1.0 - arX; else delX = 0.0; if(arY <= 1) delY = 1.0 - arY; else delY = 0.0; D = delX * delY; return D; } void cal_obs() { int x,y,i; //double dense; for(x=0; x<NX; x++) { for (y=0; y<NY; y++) { rho[x][y]=0.0; ux[x][y]=0.0; uy[x][y]=0.0; //if(bnode[x][y] !=1.0) { for (i=0; i<Q; i++) { rho[x][y] += f1[x][y][i]; ux[x][y] += f1[x][y][i]*cx[i]; uy[x][y] += f1[x][y][i]*cy[i]; } //rho[x][y] = dense; ux[x][y] = ux[x][y]/rho[x][y]; uy[x][y] = uy[x][y]/rho[x][y]; } } } return; } void IBM() { double dels = (a * pie * lmtheta)/360.0; //arc length calculation int b = 0; int markerCount = 360/lmtheta; //calculating no. of markers double UbX[markerCount], UbY[markerCount]; //velocities at markers double FbX[markerCount], FbY[markerCount]; //forces at markers double Fx[NX][NY]; float Fy[NX][NY]; //forces at grids for(int i=0; i<markerCount; ++i) UbX[i] = 0.0, //initialising to zero UbY[i] = 0.0, FbX[i] = 0.0, FbY[i] = 0.0; for(int i=0; i<NX; i++) for(int j=0; j<NY; j++) { Fx[i][j] = 0.0; Fy[i][j] = 0.0; } //double radius = a/2.0; //iterate all the markers one by one for(double theta=0; theta<=360; theta+=lmtheta, b+=1) { double radiantheta = theta * (pie / 180.0); //converting angle to radian double bX = (g + (a/2.0)*cos(radiantheta)); //calculating x-cordinate of lagrangian marker double bY = (f + (a/2.0)*sin(radiantheta)); //calculating y-cordinate of lagrangian marker int rbx = round(bX), rby = round(bY); //rounding it to nearest integer // for(int x=0; x<NX; x++) // { // for(int y=0; y<NY; y++) // { int vicinity_map[4][2] = { // x, y {-1, 0}, //left {+1, 0}, //right {0, +1}, //down {0, -1}, //up }; int m=0; while(m != NF) { //iterate over each direction int rx, ry; //temp variables to store the coords of nodes in the vicinity of current marker for(int i=0; i<4; ++i) { //add in either(x or y) coord to move in one of the 4 desired directions //Note: we are not skipping any nearby node in any of the 4 directions rx = round(bX + vicinity_map[i][0]); ry = round(bY + vicinity_map[i][1]); double D = delta2(rx, ry, bX, bY); UbX[b] += (ux[rx][ry] * D * pow(h, 2)); //unforced velocity interpolation on markers UbY[b] += (uy[rx][ry] * D * pow(h, 2)); } //OLD CODE // //for (int i = rbx - 1; i <= rbx + 1; i++) //{ // for (int j = rby - 1; j <= rby + 1; j++) // { // double distance = sqrt(pow(i-g, 2) + pow(j-f, 2)); // // if(distance > radius) //checking if the point is outside the circle // { // UbX[b] += (ux[i][j] * D * pow(h, 2)); // UbY[b] += (uy[i][j] * D * pow(h, 2)); // } // } //} int Ub = 0; //calculating boundary force on markers // printf("f: %d, %d: %f\r\n", rbx, rby, rho[rbx][rby]); FbX[b] = (2*rho[rbx][rby] * Ub - UbX[b]) ; //x-component of force on the marker from it's vicinity FbY[b] = (2*rho[rbx][rby] * Ub - UbY[b]) ; //y-component for(int i=0; i<4; ++i) { //add in either(x or y) coord to move in one of the 4 desired directions //Note: we are not skipping any nearby node in any of the 4 directions rx = round(bX + vicinity_map[i][0/*delta in x-corrd*/]); ry = round(bY + vicinity_map[i][1/*delta in y-coord*/]); double D = delta2(rx, ry, bX, bY); Fx[rx][ry] += FbX[b] * D * dels; //updating the force on neraby nodes of current markers Fy[rx][ry] += FbY[b] * D * dels; } ///OLD CODE //for (int i = rbx - 1; i <= rbx + 1; i++) //{ // for (int j = rby - 1; j <= rby + 1; j++) // { // double distance = sqrt(pow(i-g, 2) + pow(j-f, 2)); // // if(distance > radius) //checking if the point is outside the circle // { // Fx[i][j] += FbX[b] * D * dels; //updating the force on nearby nodes of the current marker // Fy[i][j] += FbY[b] * D * dels; // } // } //} // } // } m = m + 1; } } int dt = 1; //CFX = 0; //SumFX = 0.0; //calculating updated velocity on grids for(int i=0; i<NX; ++i) { for(int j=0; j<NY; ++j) { // printf("s: %d, %d: %f\r\n", i, j, rho[i][j]); ux[i][j] = ux[i][j] + (dt/(2*rho[i][j]))*Fx[i][j]; uy[i][j] = uy[i][j] + (dt/(2*rho[i][j]))*Fy[i][j]; FX[i][j] = Fx[i][j]; FY[i][j] = Fy[i][j]; //SumFX += FX[i][j]; // printf("ff: %f, %f", Fx[i][j], Fy[i][j]); } } return; } void collision() { int x,y,i; double udotc, u2, feq, force_term; for(x=0; x<NX; x++) { for(y=0; y<NY; y++) { U[x][y]=ux[x][y]; V[x][y]=uy[x][y]; if(bnode[x][y] != 1.0) { for(i=0; i<Q; i++) { //U[x][y] = ux[x][y]+(fx*0.5)/rho[x][y]; //V[x][y] = uy[x][y]+(fy*0.5)/rho[x][y]; udotc = U[x][y]*cx[i]+V[x][y]*cy[i]; u2 = U[x][y]*U[x][y]+V[x][y]*V[x][y]; force_term = (1.0-0.5/tau)*t[i]*((3.0*(cx[i]-U[x][y])+9.0*cx[i]*(cx[i]*U[x][y]+cy[i]*V[x][y]))*FX[x][y]) + ((3.0*(cy[i]-V[x][y])+9.0*cy[i]*(cx[i]*U[x][y]+cy[i]*V[x][y]))*FY[x][y]) ; feq = rho[x][y]*t[i]*(1+3.0*udotc+4.5*udotc*udotc-1.5*u2); //f1[x][y][i] = f1[x][y][i]-(1.0/tau)*(f1[x][y][i]-feq); f1[x][y][i] = f1[x][y][i]-(1.0/tau)*(f1[x][y][i]-feq) +(force_term); // Compute force for each direction //for (int i = 0; i < Q; i++) //{ // Compute deviation of velocity from discrete velocity //double dev = (c[i].x - u.x) / c*c; // Compute dot product of e_i with u //double dot = e[i].x * u.x + e[i].y * u.y + e[i].z * u.z; //dot /= c * c*c*c; // Compute dot product of previous result with e_i //dot *= e[i]; // Compute weight factor for direction i //double weight = w[i]; // Compute force for direction i //double force = (1.0 - 0.5 * tau) * weight * (3.0 * dev + 9.0 * dot / (c * c * c * c)); // Add force to total force //F.x += force * e[i].x; //F.y += force * e[i].y; //F.z += force * e[i].z; } } } } return ; } void streaming() { int x,y,i,newx,newy; double rho_inlet, u_outlet; /*for(x=0; x<NX; x++) { for(y=0; y<NY; y++) { //if(bnode[x][y]!=1.0) { for(i=0; i<Q; i++) { newx = next_x[x][i]; newy = next_y[y][i]; f1[newx][newy][i] = f2[x][y][i]; } } } }*/ for(x=0; x<NX; x++) { for (y=0; y<NY; y++) { for (i=0; i<Q; i++) { newx=next_x[x][i]; newy=next_y[y][i]; f2[newx][newy][i] = f1[x][y][i]; } } } for(x=0; x<NX; x++) { for (y=0; y<NY; y++) { for (i=0; i<Q; i++) { f1[x][y][i]=f2[x][y][i]; } } } /*------ Standard bounce-back scheme -----------*/ for(x=0; x<NX; x++) { for(y=0; y<NY; y++) { if(y==0) { f1[x][y][2] = f1[x][y][4]; f1[x][y][5] = f1[x][y][7]; f1[x][y][6] = f1[x][y][8]; } if(y==NY-1) { f1[x][y][4] = f1[x][y][2]; f1[x][y][7] = f1[x][y][5]; f1[x][y][8] = f1[x][y][6]; } } } for (y=0; y<NY; y++) { rho_inlet = (f1[0][y][0]+f1[0][y][2]+f1[0][y][4]+2.0*(f1[0][y][3]+f1[0][y][6]+f1[0][y][7]))/(1.0-u_inlet); f1[0][y][1] = f1[0][y][3]+(2.0/3.0)*rho_inlet*u_inlet; f1[0][y][5] = f1[0][y][7]-(1.0/2.0)*(f1[0][y][2]-f1[0][y][4])+(1.0/6.0)*rho_inlet*u_inlet; f1[0][y][8] = f1[0][y][6]+(1.0/2.0)*(f1[0][y][2]-f1[0][y][4])+(1.0/6.0)*rho_inlet*u_inlet; } for (y=0; y<NY; y++) { for(x=NX-1; x<NX; x++) { u_outlet = -1.0 + (f1[x][y][0]+f1[x][y][2]+f1[x][y][4]+2.0*(f1[x][y][1]+f1[x][y][5]+f1[x][y][8]))/rho0; f1[x][y][3] = f1[x][y][1]-(2.0/3.0)*rho0*u_outlet; f1[x][y][7] = f1[x][y][5]-(1.0/6.0)*rho0*u_outlet+(1.0/2.0)*(f1[x][y][2]-f1[x][y][4]); f1[x][y][6] = f1[x][y][8]-(1.0/6.0)*rho0*u_outlet-(1.0/2.0)*(f1[x][y][4]-f1[x][y][2]); } } return; } void data(int data_counter, int id) { int x, y; //double Cd; //double SumFX; FILE *f3; char phase [100]; /* Data for Tecplot*/ //SumFX = 0.0; sprintf ( phase, "IBM_phase%d.dat",id); f3 = fopen(phase, "w"); fprintf( f3, "Variables = X, Y, RHO, U, V, Fx, Fy \n"); fprintf(f3,"Zone I= %d,J= %d,F=POINT\n", NX, NY); for (y = NY - 1; y > -1; y--) { for (x = 0; x < NX; x++) { //SumFX += FX[x][y]; fprintf(f3,"%d\t%d\t%lf\t%lf\t%lf\t%0.10lf\t%0.10lf\n",x,y,rho[x][y],U[x][y],V[x][y],FX[x][y],FY[x][y]); fprintf(f3, "\n"); } } fflush(f3); fclose(f3); return; } void cal_Cd(int data_counter, int id) { int k,x,y; double SumFX; double Cd; FILE *f3; char phase [100]; /* Data for Tecplot*/ SumFX = 0.0; sprintf ( phase, "Cd%d.dat",id); f3 = fopen(phase, "w"); fprintf( f3, "Variables = k, Cd \n"); fprintf(f3,"Zone I= %d,J= %d,F=POINT\n", NX, NY); for (y = NY - 1; y > -1; y--) { for (x = 0; x < NX; x++) { SumFX += FX[x][y]; } } for(k=99999; k<T; k++) { Cd = -(SumFX)/(pow(u_inlet,2) * radius); fprintf(f3,"%d\t%lf\t%lf\t\n", k,Cd, SumFX); } fflush(f3); fclose(f3); return; } /* void data(int, int) { { //store output in a file FILE *f3; f3 = fopen("IB-LBM_Output.txt", "w"); fprintf( f3, "Variables = X, Y, RHO, U, V \n"); //double radius = a/2.0; for (int x = 0; x < NX; x++) { for (int y=0; y<NY; ++y) { //double distance = sqrt(pow(x-g, 2) + pow(x-f, 2)); //only if the node is outside the circle, we are writing to the file //if(distance > radius) { fprintf(f3,"%d\t%d\t%lf\t%lf\t%lf\n",x,y,rho[x][y],ux[x][y],uy[x][y]); fprintf(f3, "\n"); //} } } fflush(f3); fclose(f3); return; } }*/
Preview:
downloadDownload PNG
downloadDownload JPEG
downloadDownload SVG
Tip: You can change the style, width & colours of the snippet with the inspect tool before clicking Download!
Click to optimize width for Twitter