LBM Poiseuille 2D help

[jeff_yoan90]

Hello allrnrnI’m a beginner with the Lattice Boltzmann method. And I like to have your help, because I have any problems with the calculation of the pressure field in my rectangular duct for boundary conditions of Von Neumann (Velocity inlet/Outlet flow in side East and West):rnrnCBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBCrnI OrnI OrnI OrnCBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBCrnrnrnI: velocity Inlet flowrnO: velocity Outlet flowrnC: corner nodesrnB: Bounce backrnrnI would also like to know how to deal specifically with the corner nodes. Then I’d like to change my boundary conditions. At first, I’d imposed a velocity in inlet and pressure in outlet. Then I’d like an imposed pressure gradient between the inlet and outlet.rnrnHere is my code, I know if I respect much LBM algorithm, but as I have said I am a novice with this method, I come through the forum, solicit your support, your comments, corrections and your explanation.rnrnMy code:rnrn#include <stdlib.h>rn#include <stdio.h>rn#include <math.h>rn#include <string.h>rn#include <unistd.h>rnrnrn// Grid size in x and y dimensionsrn#define lx 301rn#define ly 81rnrn// Lattice Boltzmann D2Q9rn#define Q 9rnrn// D2Q9 weighting factorsrndouble t0 = 4./9.;rndouble t1 = 1./9.;rndouble t2 = 1./36.;rnrn// Fluid density per linkrn#define density 0.1rnrn// BGK collision timerndouble tau;rnrn// Collision parameterrndouble omega;rnrn// Maximum number of iterationrn#define MaxStep 100000001rnrn// Output frequencyrn#define NUMAX 1000rnrn// Length of rectangular cylinderrn#define a 21rn#define b 21rnrn// Distance from inlet boundaryrn#define x1 80rnrn// Coordonate of rectangular pointrnint x2 = x1 + a;rnrn// Computing controlrnint NUMB;rnrn// Iteration counterrnint istep;rnrn// Velocity inletrndouble Uw = 0.02;rnrn// Reynolds number Rerndouble Re = 100.;rnrn// Reference length for Reynolds numberrndouble XYMAX;rnrn// Density function of distributionrndouble f[lx][ly][Q];rnrn// ftemp…. array of temporare storagerndouble ftemp[lx][ly][Q];rnrn// Solid wall nodes arrayrnint wall[lx][ly];rnrnrnrnrnrnrn/**************************************/rn/* Define solid nodes wall boundary */rn/**************************************/rnrnvoid Wall_coordinate(void)rnrn{rnint x,y,obst_x,obst_y,obst_r;rnrnobst_x = (lx-1)/5;rnobst_y = (ly-1)/2;rnobst_r = (ly-1)/10;rnrnrnfor (y = 0; y < ly; y++) {rnrnfor (x = 0; x < lx; x++) {rnrnwall[x][y] = 0;rn}rn}rnrnrnfor (x = 0; x < lx; x++) {rnrnwall[x][0] = 1;rnwall[x][ly-1] = 1;rn}rnrnrnrnrnfor (x = 0; x < lx; x++) {rnrnfor (y = 0; y < ly; y++) {rnrnif (((x – obst_x) * (x – obst_x) + (y – obst_y) * (y – obst_y)) <= (obst_r * obst_r)) wall[x][y] = 1;rn}rn}rnrnrnrn/*for (x = x1; x < x2; x++) {rnrnwall[x][y1] = 1;rnwall[x][y2] = 1;rn}*/rnrnrn/*for (y = y1; y < y2; y++) {rnrnwall[x1][y] = 1;rnwall[x2][y] = 1;rn}*/rnrn}rnrnrnrnrn/*********************************/rn/* Initialize density function */rn/*********************************/rnrnvoid Init_density(void)rnrn{rnint x,y;rndouble w0,w1,w2;rnrn// Compute weighting factorsrnrnw0 = t0 * density;rnw1 = t1 * density;rnw2 = t2 * density;rnrn// Loop over computational domainrnrnfor (x = 0; x < lx; x++) {rnrnfor (y = 0; y < ly; y++) {rnrnf[x][y][0] = w0;rnf[x][y][1] = w1;rnf[x][y][2] = w1;rnf[x][y][3] = w1;rnf[x][y][4] = w1;rnf[x][y][5] = w2;rnf[x][y][6] = w2;rnf[x][y][7] = w2;rnf[x][y][8] = w2;rn}rn}rnrn}rnrnrnrnrnrn/***********************************************/rn/* Inlet & outlet velocity bounday condition */rn/***********************************************/rnrnvoid Inlet_outlet_BC(void)rnrn{rnrn//local variablesrnint i,x,y;rndouble C1,C2,C3,U,V,uv[Q],Dxy;rnrn//square speed of soundrnC1 = 1./3.;rnC2 = 2. * C1 * C1;rnC3 = 2. * C1;rnrnrn// Inlet velocity Bondary Conditionrnrnfor (y = 0; y < ly; y++) {rnrnrnx = 0;rnrnDxy = density;rnrn// give velocity on the boundaryrnU = Uw;rnV = 0.0;rnrn// square velocityrnuv[0] = U * U + V * V;rnuv[1] = U;rnuv[2] = V;rnuv[3] = – U;rnuv[4] = – V;rnuv[5] = U + V;rnuv[6] = – U + V;rnuv[7] = – U – V;rnuv[8] = U – V;rnrnf[x][y][0] = t0 * Dxy * (1. – uv[0]/C3);rnf[x][y][1] = t1 * Dxy * (1. + uv[1]/C1 + uv[1] * uv[1]/C2 – uv[0]/C3);rnf[x][y][2] = t1 * Dxy * (1. + uv[2]/C1 + uv[2] * uv[2]/C2 – uv[0]/C3);rnf[x][y][3] = t1 * Dxy * (1. + uv[3]/C1 + uv[3] * uv[3]/C2 – uv[0]/C3);rnf[x][y][4] = t1 * Dxy * (1. + uv[4]/C1 + uv[4] * uv[4]/C2 – uv[0]/C3);rnf[x][y][5] = t2 * Dxy * (1. + uv[5]/C1 + uv[5] * uv[5]/C2 – uv[0]/C3);rnf[x][y][6] = t2 * Dxy * (1. + uv[6]/C1 + uv[6] * uv[6]/C2 – uv[0]/C3);rnf[x][y][7] = t2 * Dxy * (1. + uv[7]/C1 + uv[7] * uv[7]/C2 – uv[0]/C3);rnf[x][y][8] = t2 * Dxy * (1. + uv[8]/C1 + uv[8] * uv[8]/C2 – uv[0]/C3);rn}rnrnrnrn// Outlet velocity Bondary Conditionrnrnfor (y = 0; y < ly; y++) {rnrnrnx = lx-1;rnrnDxy = density;rnrn// give velocity on the boundaryrnU = Uw;rnV = 0.0;rnrn// square velocityrnuv[0] = U * U + V * V;rnuv[1] = U;rnuv[2] = V;rnuv[3] = – U;rnuv[4] = – V;rnuv[5] = U + V;rnuv[6] = – U + V;rnuv[7] = – U – V;rnuv[8] = U – V;rnrnf[x][y][0] = t0 * Dxy * (1. – uv[0]/C3);rnf[x][y][1] = t1 * Dxy * (1. + uv[1]/C1 + uv[1] * uv[1]/C2 – uv[0]/C3);rnf[x][y][2] = t1 * Dxy * (1. + uv[2]/C1 + uv[2] * uv[2]/C2 – uv[0]/C3);rnf[x][y][3] = t1 * Dxy * (1. + uv[3]/C1 + uv[3] * uv[3]/C2 – uv[0]/C3);rnf[x][y][4] = t1 * Dxy * (1. + uv[4]/C1 + uv[4] * uv[4]/C2 – uv[0]/C3);rnf[x][y][5] = t2 * Dxy * (1. + uv[5]/C1 + uv[5] * uv[5]/C2 – uv[0]/C3);rnf[x][y][6] = t2 * Dxy * (1. + uv[6]/C1 + uv[6] * uv[6]/C2 – uv[0]/C3);rnf[x][y][7] = t2 * Dxy * (1. + uv[7]/C1 + uv[7] * uv[7]/C2 – uv[0]/C3);rnf[x][y][8] = t2 * Dxy * (1. + uv[8]/C1 + uv[8] * uv[8]/C2 – uv[0]/C3);rn}rnrn}rnrnrnrnrn/*************************************************************/rn/* Propagate fluid densities to their next neighbour nodes */rn/*************************************************************/rnrnvoid Stream(void)rnrn{rnrn//local variablesrnint x,y,xn,xp,yn,yp;rnrn//loop over all frnrnfor (x = 0; x < lx; x++) {rnrnxn = (x > 0)?(x-1):(lx-1);rnxp = (x < lx-1)?(x+1):(0);rnrnrnfor (y = 0; y < ly; y++) {rnrnyn = (y > 0)?(y-1):(ly-1);rnyp = (y < ly-1)?(y+1):(0);rnrnftemp[x][y][0] = f[x][y][0];rnftemp[xp][y][1] = f[x][y][1];rnftemp[x][yp][2] = f[x][y][2];rnftemp[xn][y][3] = f[x][y][3];rnftemp[x][yn][4] = f[x][y][4];rnftemp[xp][yp][5] = f[x][y][5];rnftemp[xn][yp][6] = f[x][y][6];rnftemp[xn][yn][7] = f[x][y][7];rnftemp[xp][yn][8] = f[x][y][8];rn}rn}rnrn}rnrnrnrnrnrn/*************************************************************/rn/* Fluid densities are rotated. By the next propagation step */rn/* this results in a bounce back from obstacle nodes. */rn/*************************************************************/rnrnvoid Bounceback(void)rnrn{rnrn//local variablesrnint x,y;rnrn//loop over all frnrnfor (x = 0; x < lx; x++) {rnrnfor (y = 0; y < ly; y++) {rnrnif (wall[x][y] == 1) {rnrnf[x][y][1] = ftemp[x][y][3];rnf[x][y][2] = ftemp[x][y][4];rnf[x][y][3] = ftemp[x][y][1];rnf[x][y][4] = ftemp[x][y][2];rnf[x][y][5] = ftemp[x][y][7];rnf[x][y][6] = ftemp[x][y][8];rnf[x][y][7] = ftemp[x][y][5];rnf[x][y][8] = ftemp[x][y][6];rn}rn}rn}rn}rnrnrnrnrn/***************************************/rn/* One-step density Collision process */rn/***************************************/rnrnvoid Collision(void)rnrn{rn//local variablesrnint i,x,y;rndouble C1,C2,C3,U,V,uv[Q],Dxy,feq[Q];rnrn//square speed of soundrnC1 = 1./3.;rnC2 = 2. * C1 * C1;rnC3 = 2. * C1;rnrnfor (x = 0; x < lx; x++) {rnrnfor (y = 0; y < ly; y++) {rnrnDxy = 0.;rnrn//only fluid nodes are considered herernrnif (wall[x][y] == 0) {rnrn// calculate local density Dxyrnfor (i = 0; i < Q; i++) { Dxy += ftemp[x][y]; }rnrn// calculate velocity componentsrnU = (ftemp[x][y][1] + ftemp[x][y][5] + ftemp[x][y][8] – (ftemp[x][y][3] + ftemp[x][y][6] + ftemp[x][y][7])) / Dxy;rnV = (ftemp[x][y][2] + ftemp[x][y][5] + ftemp[x][y][6] – (ftemp[x][y][4] + ftemp[x][y][7] + ftemp[x][y][8])) / Dxy;rnrn// square velocityrnuv[0] = U * U + V * V;rnuv[1] = U;rnuv[2] = V;rnuv[3] = – U;rnuv[4] = – V;rnuv[5] = U + V;rnuv[6] = – U + V;rnuv[7] = – U – V;rnuv[8] = U – V;rnrnfeq[0] = t0 * Dxy * (1. – uv[0]/C3);rnfeq[1] = t1 * Dxy * (1. + uv[1]/C1 + uv[1] * uv[1]/C2 – uv[0]/C3);rnfeq[2] = t1 * Dxy * (1. + uv[2]/C1 + uv[2] * uv[2]/C2 – uv[0]/C3);rnfeq[3] = t1 * Dxy * (1. + uv[3]/C1 + uv[3] * uv[3]/C2 – uv[0]/C3);rnfeq[4] = t1 * Dxy * (1. + uv[4]/C1 + uv[4] * uv[4]/C2 – uv[0]/C3);rnfeq[5] = t2 * Dxy * (1. + uv[5]/C1 + uv[5] * uv[5]/C2 – uv[0]/C3);rnfeq[6] = t2 * Dxy * (1. + uv[6]/C1 + uv[6] * uv[6]/C2 – uv[0]/C3);rnfeq[7] = t2 * Dxy * (1. + uv[7]/C1 + uv[7] * uv[7]/C2 – uv[0]/C3);rnfeq[8] = t2 * Dxy * (1. + uv[8]/C1 + uv[8] * uv[8]/C2 – uv[0]/C3);rnrnfor (i = 0; i < Q; i++) {f[x][y] = (1. – omega) * ftemp[x][y] + omega * feq;}rnrnrnrnrn}rn}rn}rnrnrn}rnrnrnrnrn/*********************************************************************************************/rn/* Compute total density to verify of the system that no divergence occurs when iterating */rn/*********************************************************************************************/rnrnvoid check_density(void)rnrn{rn//local variablesrnint x,y,iLB;rndouble sum = 0.;rnrnfor (x = 0; x < lx; x++)rnrn{rnrnfor (y = 0; y < ly; y++)rnrn{rnfor (iLB = 0; iLB < Q; iLB++){ sum +=f[x][y][iLB];}rn}rn}rnrnprintf(“”–> Total density in the system %lfn””,sum);rn}rnrnrnrn/********************************************************/rn/* Output of results to file ‘ densities000000.vtyk ‘ */rn/********************************************************/rnrnvoid Write_densities(void)rnrn{rn//local variablesrnint i,x,y;rndouble U,V,P,W,Dxy,Dxy0,cs,pasxyz;rnrn//square speed of soundrncs = 1./3.;rnrnpasxyz = 1./lx;rnrnprintf(“”istep = %dn””,istep);rnrnchar nomfic[25];rnrnFILE * sortie;rnsprintf(nomfic,””Visual/densities%.6i.vtk””,istep);rnsortie = fopen(nomfic, “”w””);rnrnfprintf(sortie,””# vtk DataFile Version 2.0n””);rnfprintf(sortie,””Sortie domaine LB+LINK t %dn””,istep);rnfprintf(sortie,””ASCIIn””);rnfprintf(sortie,””DATASET RECTILINEAR_GRIDn””);rnfprintf(sortie,””DIMENSIONS %d %d 1n””,lx,ly);rnfprintf(sortie,””X_COORDINATES %d floatn””,lx);rnrnfor (i = 0; i < lx; i++) { fprintf(sortie,””%e “”,(float)i * pasxyz);}rnrnfprintf(sortie,””n””);rnfprintf(sortie,””Y_COORDINATES %d floatn””,ly);rnrnfor (i = 0; i < ly; i++) { fprintf(sortie,””%e “”,(float)i * pasxyz);}rnrnfprintf(sortie,””n””);rnfprintf(sortie,””Z_COORDINATES 1 floatn””);rnfprintf(sortie,””0n””);rnrnfprintf(sortie,””POINT_DATA %dn””,lx*ly);rnfprintf(sortie,””SCALARS Pression float 1n””);rnfprintf(sortie,””LOOKUP_TABLE defaultn””);rnrnfor (y = 0; y < ly; y++) {rnrnfor (x = 0; x < lx; x++) {rnrnDxy = 0.;rnrn// pressure = average pressurernif (wall[x][y] == 1) { P = 0.;}rnrnelse {rnrn// calculate local Density Dxyrnfor (i = 0; i < Q; i++) { Dxy += f[x][y];}rnrn// calculate pressurernP = Dxy * cs;rn}rnrnfprintf(sortie,””%.4lfn””,P);rn}rn}rnrnrnfprintf(sortie,””VECTORS VecVelocity floatn””);rnrnfor (y = 0; y < ly; y++)rnrn{rnfor (x = 0; x < lx; x++)rnrn{rnDxy = 0.;rnrnif (wall[x][y] == 1)rnrn{rnU = 0.;rnV = 0.;rn}rnrnelsernrn{rn// calculate local Density Dxyrnfor (i = 0; i < Q; i++) { Dxy += f[x][y];}rnrn// calculate velocity componentsrnrnU = (f[x][y][1] + f[x][y][5] + f[x][y][8] – (f[x][y][3] + f[x][y][6] + f[x][y][7]))/Dxy;rnV = (f[x][y][2] + f[x][y][5] + f[x][y][6] – (f[x][y][4] + f[x][y][7] + f[x][y][8]))/Dxy;rn}rnrnfprintf(sortie,””%.4lf %.4lf 0.n””,U,V);rn}rn}rnrnrnfclose(sortie);rnrn}rnrnrnrnrnrnrn// Principal function main()rnrnint main(int argc, char** argv)rnrn{rnrnsystem(“”mkdir Visual””);rnrn// Begin Initialisationrnrn// Call wall coordinate functionrnWall_coordinate();rnrn// Call init density functionrnInit_density();rnrnXYMAX = x2 – x1;rnrntau = Uw * XYMAX * 3./Re + 1./2.;rnomega = 1./tau;rnRe = Uw * XYMAX * 3./(1./omega – 1./2.);rnrnprintf(“”n””);rnprintf(“”BGK collision timen””);rnprintf(“”tau = %.4lfn””,tau);rnprintf(“”n””);rnprintf(“”Reynolds numbern””);rnprintf(“”Re = %.2lfn””,Re);rnprintf(“”n””);rnprintf(“”Grid size dimensionsn””);rnprintf(“”L = %d , H = %dn””,lx,ly);rnprintf(“”n””);rnrn// Initial call check densitiesrncheck_density();rnrn// End Initialisation !!!rnrnrnrn// Begin IterationsrnrnNUMB = 0;rnrnfor (istep = 0; istep < MaxStep; istep++) {rnrnif (NUMB <= NUMAX) NUMB += 1;rnrn// Call Inlet-outlet boundary conditionrnInlet_outlet_BC();rnrn// Call Stream density propagationrnStream();rnrn// Call bounc back from solid wallrnBounceback();rnrn// Call density CollisionrnCollision();rnrn/* each MaxStep/10 iteration … */rnrnif (NUMB > NUMAX) {rnrn/* Output the data velocity, vorticity and pressure */rnrn// Call write output files vtkrnWrite_densities();rnrn// Call check densitiesrncheck_density();rnrnNUMB = 1;rn}rnrn}rnrn// End Iterations !!!rnrnrnreturn 0;rnrn}rnrnthank you very much !!!

[mathias]

Dear Jeff,rnrnat first I would try simple Bounce Back as boundary condition for all corners. I further recomment to have a look at papers of Michael Junk et al.. I remember some very helpful papers of him about that topic.rnrnMathias

[Author]

Posts