LBM BGK HELP !!!

[jeff_yoan90]

Hi all!rnrnrnI just started with the LBM. I correctly simulated 2D Poiseuille flow with a pressure gradient between the inlet and outlet (rho_in rho_out = 1.1 and = 1.0). For this, I used two distribution functions f [lx] [ly] [q] and f_old [lx] [ly] [q].rnrnIn the public void initialize (), I initialized. And the function void collision_streaming_bgk (), in the collision, I used f [x] [y] [k] éffectuer calculation, and the streaming part and the boundary conditions (pressure, bounce-back) I used f_old [x] [y] [k]. And prior to recalculate the macroscopic variables (density rho, and velocity_x velocity_y), I made a swapping between f_old [x] [y] [k] and f [x] [y] [k].rnrnI would like to now work with only one distribution function f [x] [y] [k], without using f_old [x] [y] [k], do you have any suggestions for me as that should proceed.rnrnThank you for help!!!rnrntypedef double T;rnrnT *** f, ***f_old; // LBM populations (old and new)rnT *feq; // Equilibrium pouplationrnT **rho; // Fluid densityrnT **velocity_x; // Fluid velocity x-componentrnT **velocity_y; // Fluid velocity y-componentrnrnrn// This function computes the equilibrium function from the fluid density and velocity.rnrnvoid equilibrium(T density, T vel_x, T vel_y) {rnrnT f1 = (T)3;rnT f2 = (T)9/(T)2;rnT f3 = (T)3/(T)2;rnrnfeq[0] = weights[0] * density * ((T)1 – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[1] = weights[1] * density * ((T)1 + f1 * vel_x + f2 * SQ(vel_x) – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[2] = weights[2] * density * ((T)1 + f1 * vel_y + f2 * SQ(vel_y) – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[3] = weights[3] * density * ((T)1 – f1 * vel_x + f2 * SQ(-vel_x) – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[4] = weights[4] * density * ((T)1 – f1 * vel_y + f2 * SQ(-vel_y) – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[5] = weights[5] * density * ((T)1 + f1 * ( vel_x + vel_y) + f2 * SQ( vel_x + vel_y) – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[6] = weights[6] * density * ((T)1 + f1 * (-vel_x + vel_y) + f2 * SQ(-vel_x + vel_y) – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[7] = weights[7] * density * ((T)1 – f1 * ( vel_x + vel_y) + f2 * SQ(-vel_x – vel_y) – f3 * (SQ(vel_x) + SQ(vel_y)));rnfeq[8] = weights[8] * density * ((T)1 + f1 * ( vel_x – vel_y) + f2 * SQ( vel_x – vel_y) – f3 * (SQ(vel_x) + SQ(vel_y)));rn}rnrnrnrn// Allocate memory for the fluid density, velocity and Initialize the fluid density and velocity. Start with unit density and zero velocity.rnrnvoid initialize() {rnrnrho = new T*[lx];rnvelocity_x = new T*[lx];rnvelocity_y = new T*[lx];rnrnfor (int X = 0; X < lx; ++X) {rnrnrho[X] = new T[ly];rnrnvelocity_x[X] = new T[ly];rnvelocity_y[X] = new T[ly];rn}rnrnrnfor (int X = 0; X < lx; ++X) {rnrnfor (int Y = 0; Y < ly; ++Y) {rnrnrho[X][Y] = 1.;rnrnvelocity_x[X][Y] = 0.;rnvelocity_y[X][Y] = 0.;rn}rn}rnrnrn// Allocate memory for the populations and Initialise density distribution function with equilibrium fo zero density.rnrnfeq = new T[Q];rnf = new T**[lx];rnf_old = new T**[lx];rnrnfor (int X = 0; X < lx; ++X) {rnrnf[X] = new T*[ly];rnf_old[X] = new T*[ly];rnrnfor (int Y = 0; Y < ly; ++Y) {rnrnf[X][Y] = new T[Q];rnf_old[X][Y] = new T[Q];rnrnfor (int iD = 0; iD < Q; ++iD) {rnrnfeq[iD] = 0.;rnrnf[X][Y][iD] = 0.;rnf_old[X][Y][iD] = 0.;rn}rn}rn}rnrnrnfor (int X = 0; X < lx; ++X) {rnrnfor (int Y = 0; Y < ly; ++Y) {rnrnequilibrium(rho[X][Y], velocity_x[X][Y], velocity_y[X][Y]);rnrnfor (int iD = 0; iD < Q; ++iD) {rnrnf[X][Y][iD] = feq[iD];rnf_old[X][Y][iD] = feq[iD];rn}rn}rn}rn}rnrnrnrnrnrn/// Princal Lattice Boltzmann Method LBM BGK: Computed Collision, streaming and boundaries conditions steps.rnrnvoid collision_streaming_bgk() {rnrnrn// Computation Collision, streaming and Boundaries conditions steps.rnrnfor (int X = 0; X < lx; ++X) {rnrnint Xn = (X > 0)?(X-1):(lx-1);rnint Xp = (X < lx-1)?(X+1):(0);rnrnfor (int Y = 0; Y < ly; ++Y) {rnrnint Yn = (Y > 0)?(Y-1):(ly-1);rnint Yp = (Y < ly-1)?(Y+1):(0);rnrn// Compute equilibriumrn// The pequilibrium population are computed.rnrnequilibrium(rho[X][Y], velocity_x[X][Y], velocity_y[X][Y]);rnrnrn// Compute new populationrn// This is the lattice boltzmann equation (combined collision and streaming) including external forcing.rnrnf[X][Y][0] = f[X][Y][0] * ((T)1 – omega) + feq[0] * omega;rnf[X][Y][1] = f[X][Y][1] * ((T)1 – omega) + feq[1] * omega;rnf[X][Y][2] = f[X][Y][2] * ((T)1 – omega) + feq[2] * omega;rnf[X][Y][3] = f[X][Y][3] * ((T)1 – omega) + feq[3] * omega;rnf[X][Y][4] = f[X][Y][4] * ((T)1 – omega) + feq[4] * omega;rnf[X][Y][5] = f[X][Y][5] * ((T)1 – omega) + feq[5] * omega;rnf[X][Y][6] = f[X][Y][6] * ((T)1 – omega) + feq[6] * omega;rnf[X][Y][7] = f[X][Y][7] * ((T)1 – omega) + feq[7] * omega;rnf[X][Y][8] = f[X][Y][8] * ((T)1 – omega) + feq[8] * omega;rnrnrn// Streaming stepsrnrnf_old[X][Y][0] = f[X][Y][0];rnf_old[Xp][Y][1] = f[X][Y][1];rnf_old[X][Yp][2] = f[X][Y][2];rnf_old[Xn][Y][3] = f[X][Y][3];rnf_old[X][Yn][4] = f[X][Y][4];rnf_old[Xp][Yp][5] = f[X][Y][5];rnf_old[Xn][Yp][6] = f[X][Y][6];rnf_old[Xn][Yn][7] = f[X][Y][7];rnf_old[Xp][Yn][8] = f[X][Y][8];rnrn}rn}rnrnrnrnrn// Dirichlet Pressure Boundaries conditionsrn// Zou and He computed Boundaries conditionsrnrn// Inlet (South side)rnrnfor (int X = 0; X < lx; ++X) {rnrnint Y = 0;rnrnrho[X][Y] = rho_in;rnrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = rho_in – (f_old[X][Y][0] + f_old[X][Y][1] + f_old[X][Y][3] + (T)2 * (f_old[X][Y][4] + f_old[X][Y][7] + f_old[X][Y][8]));rnrnf_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];rnf_old[X][Y][5] = f_old[X][Y][7] + (T)1/(T)6 * velocity_y[X][Y] – (T)1/(T)2 * (f_old[X][Y][1]-f_old[X][Y][3]);rnf_old[X][Y][6] = f_old[X][Y][8] + (T)1/(T)6 * velocity_y[X][Y] – (T)1/(T)2 * (f_old[X][Y][3]-f_old[X][Y][1]);rn}rnrnrn// Outlet (East side)rnrnfor (int X = 0; X < lx; ++X) {rnrnint Y = ly-1;rnrnrho[X][Y] = rho_out;rnrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = -rho_out + (f_old[X][Y][0] + f_old[X][Y][1] + f_old[X][Y][3] + (T)2 * (f_old[X][Y][2] + f_old[X][Y][5] + f_old[X][Y][6]));rnrnf_old[X][Y][4] = f_old[X][Y][2] – (T)2/(T)3 * velocity_y[X][Y];rnf_old[X][Y][7] = f_old[X][Y][5] – (T)1/(T)6 * velocity_y[X][Y] + (T)1/(T)2 * (f_old[X][Y][1]-f_old[X][Y][3]);rnf_old[X][Y][8] = f_old[X][Y][6] – (T)1/(T)6 * velocity_y[X][Y] + (T)1/(T)2 * (f_old[X][Y][3]-f_old[X][Y][1]);rn}rnrnrn// Bounce-backrn// Due to the presence of the rigid walls at x = 0 and x = lx-1, the population have to be bounced back.rnrn// left wall (West side)rnrnfor (int Y = 1; Y < ly-1; ++Y) {rnrnint X = 0;rnrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = (T)0;rnrnrho[X][Y] = velocity_x[X][Y]+f_old[X][Y][0]+f_old[X][Y][2]+f_old[X][Y][4]+(T)2*(f_old[X][Y][3]+f_old[X][Y][6]+f_old[X][Y][7]);rnf_old[X][Y][1] = f_old[X][Y][3]+(T)2/(T)3*velocity_x[X][Y];rnf_old[X][Y][5] = f_old[X][Y][7]+(T)1/(T)6*velocity_x[X][Y]+(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])+(T)1/(T)2*velocity_y[X][Y];rnf_old[X][Y][8] = f_old[X][Y][6]+(T)1/(T)6*velocity_x[X][Y]-(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])-(T)1/(T)2*velocity_y[X][Y];rn}rnrnrnrn// right wall (East side)rnrnfor (int Y = 1; Y < ly-1; ++Y) {rnrnint X = lx-1;rnrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = (T)0;rnrnrho[X][Y] = -velocity_x[X][Y]+f_old[X][Y][0]+f_old[X][Y][2]+f_old[X][Y][4]+(T)2*(f_old[X][Y][1]+f_old[X][Y][5]+f_old[X][Y][8]);rnf_old[X][Y][3] = f_old[X][Y][1]-(T)2/(T)3*velocity_x[X][Y];rnf_old[X][Y][7] = f_old[X][Y][5]-(T)1/(T)6*velocity_x[X][Y]-(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])-(T)1/(T)2*velocity_y[X][Y];rnf_old[X][Y][6] = f_old[X][Y][8]-(T)1/(T)6*velocity_x[X][Y]+(T)1/(T)2*(f_old[X][Y][4]-f_old[X][Y][2])+(T)1/(T)2*velocity_y[X][Y];rn}rnrnrnrnrn// Corners nodes treatmentrnrn// Bottom left corner nodern{rnint X = 0;rnint Y = 0;rnrnrho[X][Y] = rho[X][Y+1]; // Extrapolation 1st orderrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = (T)0;rnrnf_old[X][Y][1] = f_old[X][Y][3] + (T)2/(T)3 * velocity_x[X][Y];rnf_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];rnf_old[X][Y][5] = f_old[X][Y][7] + (T)1/(T)6 * (velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][6] = (T)1/(T)12 * (-velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][8] = (T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);rnf_old[X][Y][0] = (T)0;rnrnf_old[X][Y][0] = rho[X][Y] – (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);rn}rnrnrn// Top left corner nodern{rnint X = 0;rnint Y = ly-1;rnrnrho[X][Y] = rho[X][Y-1]; // Extrapolation 1st orderrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = (T)0;rnrnf_old[X][Y][1] = f_old[X][Y][3] + (T)2/(T)3 * velocity_x[X][Y];rnf_old[X][Y][4] = f_old[X][Y][2] – (T)2/(T)3 * velocity_y[X][Y];rnf_old[X][Y][8] = f_old[X][Y][6] + (T)1/(T)6 * (velocity_x[X][Y]-velocity_y[X][Y]);rnf_old[X][Y][5] = (T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][7] = -(T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][0] = (T)0;rnrnf_old[X][Y][0] = rho[X][Y] – (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);rn}rnrnrnrnrn// Bottom right corner nodern{rnint X = lx-1;rnint Y = 0;rnrnrho[X][Y] = rho[X][Y+1]; // Extrapolation 1st orderrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = (T)0;rnrnf_old[X][Y][3] = f_old[X][Y][1] – (T)2/(T)3 * velocity_x[X][Y];rnf_old[X][Y][2] = f_old[X][Y][4] + (T)2/(T)3 * velocity_y[X][Y];rnf_old[X][Y][6] = f_old[X][Y][8] + (T)1/(T)6 * (-velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][5] = (T)1/(T)12 * (velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][7] = -(T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);rnf_old[X][Y][0] = (T)0;rnrnf_old[X][Y][0] = rho[X][Y] – (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);rnrn}rnrnrn// Top right corner nodern{rnint X = lx-1;rnint Y = ly-1;rnrnrho[X][Y] = rho[X][Y-1]; // Extrapolation 1st orderrnvelocity_x[X][Y] = (T)0;rnvelocity_y[X][Y] = (T)0;rnrnf_old[X][Y][3] = f_old[X][Y][1] – (T)2/(T)3 * velocity_x[X][Y];rnf_old[X][Y][4] = f_old[X][Y][2] – (T)2/(T)3 * velocity_y[X][Y];rnf_old[X][Y][7] = f_old[X][Y][5] – (T)1/(T)6 * (velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][6] = (T)1/(T)12 * (-velocity_x[X][Y]+velocity_y[X][Y]);rnf_old[X][Y][8] = (T)1/(T)12 * (velocity_x[X][Y]-velocity_y[X][Y]);rnf_old[X][Y][0] = (T)0;rnrnf_old[X][Y][0] = rho[X][Y] – (f_old[X][Y][1]+f_old[X][Y][2]+f_old[X][Y][3]+f_old[X][Y][4]+f_old[X][Y][5]+f_old[X][Y][6]+f_old[X][Y][7]+f_old[X][Y][8]);rnrn}rnrnrn// Swap popualationsrn// The present code used old and new populations which are swapped at the beginning of each time step.rn// This way, the old populations are not overwritten during propagation.rn// The resulting code is easier to write and debug.rnrnT ***swap_tmp;rnrnswap_tmp = f_old;rnf_old = f;rnf = swap_tmp;rnrnrn// Compute fluid density and velocityrn// This function computes the fluid density and velocity from the populations.rnrnmomenta();rn}rnrnrnrnrn// Compute fluid density and velocityrnrnvoid momenta() {rnrnfor (int X = 0; X < lx; ++X) {rnrnfor (int Y = 0; Y < ly; ++Y) {rnrnrho[X][Y] = f[X][Y][0] + f[X][Y][1] + f[X][Y][2] + f[X][Y][3] + f[X][Y][4] + f[X][Y][5] + f[X][Y][6] + f[X][Y][7] + f[X][Y][8];rnvelocity_x[X][Y] = (f[X][Y][1] + f[X][Y][5] + f[X][Y][8] – (f[X][Y][3] + f[X][Y][6] + f[X][Y][7])) /rho[X][Y];rnvelocity_y[X][Y] = (f[X][Y][2] + f[X][Y][5] + f[X][Y][6] – (f[X][Y][4] + f[X][Y][7] + f[X][Y][8])) /rho[X][Y];rn}rn}rn}

[jo]

Hi,rnrnhave a look at http://www.sciencedirect.com/science/article/pii/S0898122112003835.rnrnJo

[Author]

Posts