diff options
| -rw-r--r-- | CMakeLists.txt | 12 | ||||
| -rw-r--r-- | channel.cc | 115 | ||||
| -rw-r--r-- | lid_driven_cavity.cc | 3 | ||||
| -rw-r--r-- | lid_driven_cavity_with_obstacles.cc | 1 |
4 files changed, 126 insertions, 5 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 4248a32..bf8dc7b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,6 +6,7 @@ add_library( src/lbm.cc src/fluid_buffer.cc src/boundary_conditions.cc + src/box_obstacle.cc ) target_compile_features( @@ -39,10 +40,19 @@ target_link_libraries( add_executable( lid_driven_cavity_with_obstacles lid_driven_cavity_with_obstacles.cc - src/box_obstacle.cc ) target_link_libraries( lid_driven_cavity_with_obstacles boltzbub ) + +add_executable( + channel + channel.cc +) + +target_link_libraries( + channel + boltzbub +) diff --git a/channel.cc b/channel.cc new file mode 100644 index 0000000..837b355 --- /dev/null +++ b/channel.cc @@ -0,0 +1,115 @@ +#include <iostream> +#include <vector> +#include <algorithm> + +#include "lbm.h" +#include "boundary_conditions.h" +#include "box_obstacle.h" + +constexpr std::size_t dimX = 500; +constexpr std::size_t dimY = 40; + +constexpr double uWall = 0.2; +constexpr double reynolds = 500; + +constexpr double tau = 3. * uWall * (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}, +}; + +void init() { + for ( std::size_t x = 0; x < dimX; ++x ) { + for ( std::size_t y = 0; y < dimY; ++y ) { + fluid.density(x,y) = 1.0; + fluid.velocity(x,y) = { 0.0, 0.0 }; + + static_cast<Cell&>(pop.curr(x,y)).equilibrize( + fluid.density(x,y), fluid.velocity(x,y)); + static_cast<Cell&>(pop.prev(x,y)).equilibrize( + fluid.density(x,y), fluid.velocity(x,y)); + } + } +} + +void computeLbmStep() { + pop.swap(); + + for ( std::size_t x = 0; x < dimX; ++x ) { + for ( std::size_t y = 0; y < dimY; ++y ) { + if ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) { + return !o.isInside(x, y); + }) ) { + streamFluidCell(pop, x, y); + } + } + } + + // 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 ) { + computeWallCell(pop, {x, 0 }, { 0, 1}); + computeWallCell(pop, {x, dimY-1}, { 0,-1}); + } + + // obstacles + for ( const auto& box : obstacles ) { + box.applyBoundary(pop); + } + + 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 ( std::all_of(obstacles.cbegin(), obstacles.cend(), [x, y](const auto& o) { + return !o.isInside(x, y); + }) ) { + collideFluidCell(omega, pop, fluid, x, y); + } + } + } +} + +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; + + for ( std::size_t t = 0; t <= 10000; ++t ) { + computeLbmStep(); + + if ( t % 100 == 0 ) { + std::cout << "."; + std::cout.flush(); + fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk"); + } + } + + std::cout << std::endl; +} diff --git a/lid_driven_cavity.cc b/lid_driven_cavity.cc index ca4b030..79f330b 100644 --- a/lid_driven_cavity.cc +++ b/lid_driven_cavity.cc @@ -1,7 +1,4 @@ #include <iostream> -#include <vector> -#include <string> -#include <algorithm> #include "lbm.h" #include "boundary_conditions.h" diff --git a/lid_driven_cavity_with_obstacles.cc b/lid_driven_cavity_with_obstacles.cc index 092814e..55f519e 100644 --- a/lid_driven_cavity_with_obstacles.cc +++ b/lid_driven_cavity_with_obstacles.cc @@ -1,6 +1,5 @@ #include <iostream> #include <vector> -#include <string> #include <algorithm> #include "lbm.h" |