Preview:
#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;
        }
}*/
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