//
// main.c
// Re40Dia20NF1D2
//
// Created by Siddhartha Shandilya on 15/04/23.
//
#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.2093333333
#define NX 801 //no. of nodes in the x-direction
#define NY 801 //no. of nodes in the y-direction
#define a 20.0 //particle diameter
#define radius 10.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 1001 // no. of time steps
#define tau 0.814 // 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 Cd;
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 data_Cd(int, int);
//void exact();
int main()
{
int count = 0;
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==1000)
{
data_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];
double SumFX; //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];
Cd = -(SumFX)/(pow(u_inlet,2) * radius);
printf("%d\t%lf\t\n", k,Cd);
// 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 data_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 = Time Steps, Cd \n");
//fprintf(f3,"Zone I= %d,J= %d,F=POINT\n", NX, NY);
for(k=0; k<T; k++)
{
fprintf(f3,"%d\t%lf\t\n", k,Cd);
}
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;
}
}*/