From 369159ceb811cec067f194d75ffef7950b20e0ac Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Tue, 19 Jun 2018 22:03:25 +0200 Subject: Split collide and stream steps Simplifies correct implementation of Zou He boundary condition (used for the lid in lid driven cavity). Fix tau using lid velocity and Reynolds number. Simulation now seems to converge for Reynolds numbers greater than 100 if lattice resolution is sufficient. --- boltz.cc | 118 ++++++++++++++++++++++++++++++++++++--------------------------- 1 file changed, 67 insertions(+), 51 deletions(-) diff --git a/boltz.cc b/boltz.cc index 3949efa..a361a58 100644 --- a/boltz.cc +++ b/boltz.cc @@ -92,13 +92,32 @@ class CellBuffer { } }; +std::pair, Vector> neighbors(Vector v) { + if ( v[0] == 0 ) { + return { + { -1, v[1] }, + { 1, v[1] } + }; + } else if ( v[1] == 0 ) { + return { + { v[0], -1 }, + { v[0], 1 } + }; + } else { + return { + { 0, v[1] }, + { v[0], 0 } + }; + } +} + constexpr std::size_t dimX = 256; -constexpr std::size_t dimY = 256; +constexpr std::size_t dimY = dimX; -constexpr double tau = 2./3.; -constexpr int reynolds = 1000; -constexpr double lidVelocityX = reynolds * (2.*tau - 1.) / (6.*(dimX - 1.)); +constexpr double uLid = 0.1; +constexpr double reynolds = 1000; +constexpr double tau = 3 * uLid * (dimX-1) / reynolds + 0.5; constexpr double omega = 1. / tau; CellBuffer pop(dimX, dimY); @@ -106,47 +125,48 @@ CellBuffer pop(dimX, dimY); Density density [dimX][dimY]; Velocity velocity[dimX][dimY]; -void computeFluidCell(std::size_t x, std::size_t y) { - // compute equilibrium - Cell eq; - eq.equilibrize(density[x][y], velocity[x][y]); - - // collide (BGK, relax towards equilibrium) & stream +void streamFluidCell(std::size_t x, std::size_t y) { if ( x != 0 && x != dimX - 1 && y != 0 && y != dimY -1 ) { + // stream internal cells 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)); + pop.curr(x+i,y+j).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 + // stream boundary cells, + // missing 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)); + pop.curr(x+i,y+j).get(i,j) = pop.prev(x,y).get(i,j); } } } } } -std::pair, Vector> neighbors(Vector v) { - if ( v[0] == 0 ) { - return { - { -1, v[1] }, - { 1, v[1] } - }; - } else if ( v[1] == 0 ) { - return { - { v[0], -1 }, - { v[0], 1 } - }; +void collideFluidCell(std::size_t x, std::size_t y) { + // compute equilibrium + Cell eq; + eq.equilibrize(density[x][y], velocity[x][y]); + + // collide (BGK, relax towards equilibrium) + 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,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j)); + } + } } else { - return { - { 0, v[1] }, - { v[0], 0 } - }; + // partial collide for boundary cells + 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,y).get(i,j) = pop.curr(x,y).get(i,j) + omega * (eq.get(i,j) - pop.curr(x,y).get(i,j)); + } + } + } } } @@ -158,14 +178,6 @@ void computeWallCell(Vector cell, Vector 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).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); @@ -184,17 +196,18 @@ void computeZouHeVelocityWallCell(Vector cell, Vector normal, void computeLbmStep(std::size_t t) { pop.swap(); - for ( std::size_t x = 1; x < dimX-1; ++x ) { + for ( std::size_t x = 0; x < dimX; ++x ) { for ( std::size_t y = 0; y < dimY; ++y ) { - computeFluidCell(x, y); + streamFluidCell(x, y); } } // straight wall cell bounce back for ( std::size_t x = 0; x < dimX; ++x ) { + computeZouHeVelocityWallCell({x, dimY-1}, { 0,-1}, uLid); + } + for ( std::size_t x = 1; x < dimX-1; ++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 = 1; y < dimY-1; ++y ) { computeWallCell({0, y}, { 1, 0}); @@ -202,16 +215,15 @@ void computeLbmStep(std::size_t t) { } // edge wall cell bounce back - //computeWallCell({0, dimY-1}, { 1,-1}); - //computeWallCell({dimX-1, dimY-1}, {-1,-1}); - //computeWallCell({0, 0 }, { 1, 1}); - //computeWallCell({dimX-1, 0 }, {-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 ) { for ( std::size_t y = 0; y < dimY; ++y ) { density[x][y] = pop.curr(x,y).sum(); velocity[x][y] = pop.curr(x,y).velocity(density[x][y]); + + collideFluidCell(x, y); } } } @@ -273,15 +285,19 @@ void writeCurrentStateAsVTK(int time) { int main() { init(); - std::cout << "tau: " << tau << std::endl; - std::cout << "u: " << lidVelocityX << std::endl; + std::cout << "Re: " << reynolds << std::endl; + std::cout << "uLid: " << uLid << std::endl; + std::cout << "tau: " << tau << std::endl; - for ( std::size_t t = 0; t <= 40000; ++t ) { + for ( std::size_t t = 0; t <= 100000; ++t ) { computeLbmStep(t); - if ( t < 100 || t % 1000 == 0 ) { - std::cout << "Writing " << t << "." << std::endl; + if ( (t < 1000 && t % 10 == 0) || t % 1000 == 0 ) { + std::cout << "."; + std::cout.flush(); writeCurrentStateAsVTK(t); } } + + std::cout << std::endl; } -- cgit v1.2.3