From c56162c19cca84301813aff9b1326cb2c686bf48 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Thu, 5 Apr 2018 22:45:57 +0200 Subject: Lid driven cavity with bounce-back, Zou He boundary conditions --- boltz.cc | 133 +++++++++++++++++++++++++++++---------------------------------- 1 file changed, 61 insertions(+), 72 deletions(-) diff --git a/boltz.cc b/boltz.cc index 29f8518..3796b87 100644 --- a/boltz.cc +++ b/boltz.cc @@ -27,7 +27,6 @@ struct DataCell { }; using Velocity = Vector; -using Force = Vector; using Density = double; constexpr DataCell weight{ @@ -93,13 +92,12 @@ class CellBuffer { } }; -constexpr std::size_t dimX = 128; -constexpr std::size_t dimY = 128; +constexpr std::size_t dimX = 256; +constexpr std::size_t dimY = 256; -constexpr double lidVelocityX = 0.05; -constexpr double reynolds = 1000; -constexpr double nu = lidVelocityX * dimX / reynolds; -constexpr double tau = 3*nu + 0.5; +constexpr double tau = 2./3.; +constexpr int reynolds = 1000; +constexpr double lidVelocityX = reynolds * (2.*tau - 1.) / (6.*(dimX - 1.)); constexpr double omega = 1. / tau; @@ -107,7 +105,6 @@ CellBuffer pop(dimX, dimY); Density density [dimX][dimY]; Velocity velocity[dimX][dimY]; -Force force [dimX][dimY]; void computeFluidCell(std::size_t x, std::size_t y) { // compute equilibrium @@ -115,9 +112,21 @@ void computeFluidCell(std::size_t x, std::size_t y) { eq.equilibrize(density[x][y], velocity[x][y]); // collide (BGK, relax towards equilibrium) & stream - for ( int i = -1; i <= 1; ++i ) { - for ( int j = -1; j <= 1; ++j ) { - pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j)); + if ( x != 0 && x != dimX - 1 && y != 0 && y != dimY -1 ) { + for ( int i = -1; i <= 1; ++i ) { + for ( int j = -1; j <= 1; ++j ) { + pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j)); + } + } + } else { + // partial collide and stream for boundary cells, + // remaining populations to be determined by boundary conditions + for ( int i = -1; i <= 1; ++i ) { + for ( int j = -1; j <= 1; ++j ) { + if ( x+i >= 0 && x+i <= dimX - 1 && y+j >= 0 && y+j <= dimY - 1 ) { + pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j)); + } + } } } } @@ -144,78 +153,59 @@ std::pair, Vector> neighbors(Vector v) { void computeWallCell(Vector cell, Vector normal) { const auto [neighborA, neighborB] = neighbors(normal); - pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA); - pop.curr(cell).get(normal ) = pop.curr(cell-normal ).get(-normal ); - pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB); + pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA); + pop.curr(cell).get(normal ) = pop.curr(cell).get(-normal ); + pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB); } void computeSimpleVelocityWallCell(Vector cell, Vector normal, double vX) { const auto [neighborA, neighborB] = neighbors(normal); - pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA) - 6 * weight.get( 1, 1) * density[cell[0]][cell[1]] * vX; - pop.curr(cell).get(normal ) = pop.curr(cell-normal ).get(-normal ); - pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB) + 6 * weight.get(-1, 1) * density[cell[0]][cell[0]] * vX; + pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA) - 6 * weight.get( 1, 1) * density[cell[0]][cell[1]] * vX; + pop.curr(cell).get(normal ) = pop.curr(cell).get(-normal ); + pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) + 6 * weight.get(-1, 1) * density[cell[0]][cell[0]] * vX; } void computeZouHeVelocityWallCell(Vector cell, Vector normal, double vX) { const auto [neighborA, neighborB] = neighbors(normal); - const double rho = pop.curr(cell-neighborA).get(-1,0) + pop.curr(cell-normal).get(0,0) + pop.curr(cell-neighborB).get(1,0) + const double rho = pop.curr(cell).get(-1,0) + pop.curr(cell).get(0,0) + pop.curr(cell).get(1,0) + 2.*( - pop.curr(cell-neighborA).get(-neighborA) + - pop.curr(cell-normal ).get(-normal ) + - pop.curr(cell-neighborB).get(-neighborB) + pop.curr(cell).get(-neighborA) + + pop.curr(cell).get(-normal ) + + pop.curr(cell).get(-neighborB) ); - const double hor = pop.curr(cell-neighborA).get(-1,0) - pop.curr(cell-neighborB).get(1,0); - pop.curr(cell).get(neighborA) = pop.curr(cell-neighborA).get(-neighborA) + 0.5*hor - 0.5*vX*rho; - pop.curr(cell).get(normal ) = pop.curr(cell-normal ).get(-normal ); - pop.curr(cell).get(neighborB) = pop.curr(cell-neighborB).get(-neighborB) - 0.5*hor + 0.5*vX*rho; + pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA) + 0.5*(pop.curr(cell).get( 1,0) - pop.curr(cell).get(-1,0) - vX*rho); + pop.curr(cell).get(normal ) = pop.curr(cell).get(-normal ); + pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) + 0.5*(pop.curr(cell).get(-1,0) - pop.curr(cell).get( 1,0) + vX*rho); } void computeLbmStep(std::size_t t) { pop.swap(); - for ( std::size_t x = 1; x < dimX - 1; ++x ) { - for ( std::size_t y = 1; y < dimY - 1; ++y ) { - //if ( x <= 20 || x >= 40 || y <= 20 || y >= 80 ) { - computeFluidCell(x, y); - //} + for ( std::size_t x = 1; x < dimX-1; ++x ) { + for ( std::size_t y = 0; y < dimY; ++y ) { + computeFluidCell(x, y); } } - // obstacle - /*{ - for ( std::size_t x = 21; x < 40; ++x ) { - computeWallCell(x, 80, { 0, 1}); - computeWallCell(x, 20, { 0,-1}); - } - for ( std::size_t y = 21; y < 80; ++y ) { - computeWallCell(40, y, { 1, 0}); - computeWallCell(20, y, {-1, 0}); - } - - computeWallCell(40, 80, { 1, 1}); - computeWallCell(40, 20, { 1,-1}); - computeWallCell(20, 80, {-1, 1}); - computeWallCell(20, 20, {-1,-1}); - }*/ - // straight wall cell bounce back - for ( std::size_t x = 2; x < dimX - 2; ++x ) { - computeWallCell({x, 1}, { 0, 1}); - computeZouHeVelocityWallCell({x, dimY-2}, { 0,-1}, lidVelocityX); + for ( std::size_t x = 0; x < dimX; ++x ) { + computeWallCell({x, 0}, { 0, 1}); + computeZouHeVelocityWallCell({x, dimY-1}, { 0,-1}, lidVelocityX); + //computeSimpleVelocityWallCell({x, dimY-1}, { 0,-1}, lidVelocityX); } - for ( std::size_t y = 2; y < dimY - 2; ++y ) { - computeWallCell({1, y}, { 1, 0}); - computeWallCell({dimX-2, y}, {-1, 0}); + for ( std::size_t y = 1; y < dimY-1; ++y ) { + computeWallCell({0, y}, { 1, 0}); + computeWallCell({dimX-1, y}, {-1, 0}); } // edge wall cell bounce back - computeWallCell({1, 1 }, { 1, 1}); - computeWallCell({1, dimY-2}, { 1,-1}); - computeWallCell({dimX-2, 1 }, {-1, 1}); - computeWallCell({dimX-2, dimY-2}, {-1,-1}); + //computeWallCell({0, dimY-1}, { 1,-1}); + //computeWallCell({dimX-1, dimY-1}, {-1,-1}); + //computeWallCell({0, 0 }, { 1, 1}); + //computeWallCell({dimX-1, 0 }, {-1, 1}); // update density, velocity field for ( std::size_t x = 0; x < dimX; ++x ) { @@ -229,9 +219,8 @@ void computeLbmStep(std::size_t t) { void init() { for ( std::size_t x = 0; x < dimX; ++x ) { for ( std::size_t y = 0; y < dimY; ++y ) { - density[x][y] = 1.0; velocity[x][y] = { 0.0, 0.0 }; - force[x][y] = { 0.0, 0.0 }; + density[x][y] = 1.0; pop.curr(x,y).equilibrize(density[x][y], velocity[x][y]); pop.prev(x,y).equilibrize(density[x][y], velocity[x][y]); @@ -247,33 +236,33 @@ void writeCurrentStateAsVTK(int time) { fout << "lbm_output\n"; fout << "ASCII\n"; fout << "DATASET RECTILINEAR_GRID\n"; - fout << "DIMENSIONS " << dimX - 2 << " " << dimY - 2 << " 1" << "\n"; + fout << "DIMENSIONS " << dimX << " " << dimY << " 1" << "\n"; - fout << "X_COORDINATES " << dimX - 2 << " float\n"; - for( std::size_t x = 1; x < dimX - 1; ++x ) { + fout << "X_COORDINATES " << dimX << " float\n"; + for( std::size_t x = 0; x < dimX; ++x ) { fout << x << " "; } - fout << "\nY_COORDINATES " << dimY - 2 << " float\n"; - for( std::size_t y = 1; y < dimY - 1; ++y ) { + fout << "\nY_COORDINATES " << dimY << " float\n"; + for( std::size_t y = 0; y < dimY; ++y ) { fout << y << " "; } fout << "\nZ_COORDINATES " << 1 << " float\n"; fout << 0 << "\n"; - fout << "POINT_DATA " << (dimX - 2) * (dimY - 2) << "\n"; + fout << "POINT_DATA " << dimX * dimY << "\n"; fout << "VECTORS velocity float\n"; - for ( std::size_t y = 1; y < dimY - 1; ++y ) { - for ( std::size_t x = 1; x < dimX - 1; ++x ) { + for ( std::size_t y = 0; y < dimY; ++y ) { + for ( std::size_t x = 0; x < dimX; ++x ) { fout << velocity[x][y][0] << " " << velocity[x][y][1] << " 0\n"; } } fout << "SCALARS density float 1\n"; fout << "LOOKUP_TABLE default\n"; - for ( std::size_t y = 1; y < dimY - 1; ++y ) { - for ( std::size_t x = 1; x < dimX - 1; ++x ) { + for ( std::size_t y = 0; y < dimY; ++y ) { + for ( std::size_t x = 0; x < dimX; ++x ) { fout << density[x][y] << "\n"; } } @@ -284,13 +273,13 @@ void writeCurrentStateAsVTK(int time) { int main() { init(); - std::cout << "u: " << lidVelocityX << std::endl; std::cout << "tau: " << tau << std::endl; + std::cout << "u: " << lidVelocityX << std::endl; - for ( std::size_t t = 0; t < 20000; ++t ) { + for ( std::size_t t = 0; t <= 40000; ++t ) { computeLbmStep(t); - if ( t % 1000 == 0 ) { + if ( t < 100 || t % 1000 == 0 ) { std::cout << "Writing " << t << "." << std::endl; writeCurrentStateAsVTK(t); } -- cgit v1.2.3