diff options
| -rw-r--r-- | channel.cc | 55 |
1 files changed, 22 insertions, 33 deletions
@@ -9,18 +9,19 @@ constexpr std::size_t dimX = 500; constexpr std::size_t dimY = 40; -constexpr double uWall = 0.2; -constexpr double reynolds = 500; +constexpr double uInflow = 0.1; +constexpr double reynolds = 1000; -constexpr double tau = 3. * uWall * (dimX-1) / reynolds + 0.5; +constexpr double tau = 3. * uInflow * (dimX-1) / reynolds + 0.5; constexpr double omega = 1. / tau; DataCellBuffer pop(dimX, dimY); FluidBuffer fluid(dimX, dimY); std::vector<BoxObstacle> obstacles{ - {300, 0, 320, 25}, - {340, 15, 360, 39}, + {100, 0, 120, 25}, + {140, 15, 160, 39}, + {300, 15, 340, 25}, }; void init() { @@ -50,25 +51,7 @@ void computeLbmStep() { } } - // periodic boundary - for ( std::size_t y = 1; y < dimY-1; ++y ) { - pop.curr(0,y+1).get(1, 1) = pop.prev(dimX-1,y).get( 1, 1); - pop.curr(0,y ).get(1, 0) = pop.prev(dimX-1,y).get( 1, 0); - pop.curr(0,y-1).get(1,-1) = pop.prev(dimX-1,y).get( 1,-1); - } - for ( std::size_t y = 1; y < dimY-1; ++y ) { - pop.curr(dimX-1,y+1).get(-1, 1) = pop.prev(0,y).get(-1, 1); - pop.curr(dimX-1,y ).get(-1, 0) = pop.prev(0,y).get(-1, 0); - pop.curr(dimX-1,y-1).get(-1,-1) = pop.prev(0,y).get(-1,-1); - } - - // straight wall cell bounce back - for ( std::size_t x = 0; x < 100; ++x ) { - computeZouHeVelocityWallCell(pop, {x, 0 }, { 0, 1}, uWall); - computeZouHeVelocityWallCell(pop, {x, dimY-1}, { 0,-1}, uWall); - } - - for ( std::size_t x = 100; x < dimX-1; ++x ) { + for ( std::size_t x = 0; x < dimX; ++x ) { computeWallCell(pop, {x, 0 }, { 0, 1}); computeWallCell(pop, {x, dimY-1}, { 0,-1}); } @@ -80,9 +63,15 @@ void computeLbmStep() { for ( std::size_t x = 0; x < dimX; ++x ) { for ( std::size_t y = 0; y < dimY; ++y ) { - Cell& cell = static_cast<Cell&>(pop.curr(x,y)); - fluid.density(x,y) = cell.sum(); - fluid.velocity(x,y) = cell.velocity(fluid.density(x,y)); + if ( x == 0 ) { + // inflow + fluid.density(x,y) = 1.0; + fluid.velocity(x,y) = { uInflow, 0 }; + } else { + Cell& cell = static_cast<Cell&>(pop.curr(x,y)); + fluid.density(x,y) = cell.sum(); + fluid.velocity(x,y) = cell.velocity(fluid.density(x,y)); + } if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) { return !o.isInside(x, y); @@ -96,15 +85,15 @@ void computeLbmStep() { int main() { init(); - std::cout << "Re: " << reynolds << std::endl; - std::cout << "uWall: " << uWall << std::endl; - std::cout << "tau: " << tau << std::endl; - std::cout << "omega: " << omega << std::endl; + std::cout << "Re: " << reynolds << std::endl; + std::cout << "uInflow: " << uInflow << std::endl; + std::cout << "tau: " << tau << std::endl; + std::cout << "omega: " << omega << std::endl; - for ( std::size_t t = 0; t <= 10000; ++t ) { + for ( std::size_t t = 0; t <= 5000; ++t ) { computeLbmStep(); - if ( t % 100 == 0 ) { + if ( t % 50 == 0 ) { std::cout << "."; std::cout.flush(); fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk"); |