From 2bcbae122b141147e472bdbbad219c3571f71d0b Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Tue, 16 Oct 2018 12:44:24 +0200 Subject: Implement moving wall / velocity Dirichlet boundary condition Usable as both the inflow condition of the channel example and the top wall of a lid driven cavity. --- README.md | 2 +- channel.cc | 30 ++++++++++++++++++------------ lid_driven_cavity.cc | 8 +++++--- src/boundary_conditions.cc | 17 +++++++++++++++++ src/boundary_conditions.h | 2 ++ src/vector.h | 5 +++++ 6 files changed, 48 insertions(+), 16 deletions(-) diff --git a/README.md b/README.md index 385627b..5ab879a 100644 --- a/README.md +++ b/README.md @@ -10,7 +10,7 @@ | - | - | | `lid_driven_cavity.cc` | Lid driven cavity using Zou/He boundary conditions for the top wall and simple bounce back for all other walls | | `lid_driven_cavity_with_obstacles.cc` | Same as `lid_driven_cavity.cc` but includes a grid of boxes to make things more interesting | -| `channel.cc` | Channel flow some obstacles and Dirichlet inflow | +| `channel.cc` | Channel flow with some obstacles and Dirichlet inflow condition | ## Build diff --git a/channel.cc b/channel.cc index 3e7641b..b1a6166 100644 --- a/channel.cc +++ b/channel.cc @@ -9,8 +9,8 @@ constexpr std::size_t dimX = 500; constexpr std::size_t dimY = 40; -constexpr double uInflow = 0.1; -constexpr double reynolds = 1000; +constexpr double uInflow = 0.01; +constexpr double reynolds = 100; constexpr double tau = 3. * uInflow * (dimX-1) / reynolds + 0.5; constexpr double omega = 1. / tau; @@ -56,6 +56,10 @@ void computeLbmStep() { computeWallCell(pop, {x, dimY-1}, { 0,-1}); } + for ( std::size_t y = 1; y < dimY-1; ++y ) { + computeMovingWallCell(pop, {0,y}, {1,0}, {uInflow,0}); + } + // obstacles for ( const auto& box : obstacles ) { box.applyBoundary(pop); @@ -63,16 +67,18 @@ void computeLbmStep() { for ( std::size_t x = 0; x < dimX; ++x ) { for ( std::size_t y = 0; y < dimY; ++y ) { - if ( x == 0 && y > 0 && y < dimY-1 ) { - // inflow - fluid.density(x,y) = 1.0; - fluid.velocity(x,y) = { uInflow, 0 }; - } else { - Cell& cell = static_cast(pop.curr(x,y)); - fluid.density(x,y) = cell.sum(); - fluid.velocity(x,y) = cell.velocity(fluid.density(x,y)); + Cell& cell = static_cast(pop.curr(x,y)); + + // bulk density + fluid.density(x,y) = cell.sum(); + + // outflow + if ( x == dimX-1 && y > 0 && y < dimY-1 ) { + fluid.density(x,y) = 1.0; } + 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); }) ) { @@ -90,10 +96,10 @@ int main() { std::cout << "tau: " << tau << std::endl; std::cout << "omega: " << omega << std::endl; - for ( std::size_t t = 0; t <= 5000; ++t ) { + for ( std::size_t t = 0; t <= 10000; ++t ) { computeLbmStep(); - if ( t % 100 == 0 ) { + if ( t % 1000 == 0 ) { std::cout << "."; std::cout.flush(); fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk"); diff --git a/lid_driven_cavity.cc b/lid_driven_cavity.cc index 79f330b..7a978ad 100644 --- a/lid_driven_cavity.cc +++ b/lid_driven_cavity.cc @@ -38,10 +38,12 @@ void computeLbmStep() { } } - // straight wall cell bounce back + // moving top wall for ( std::size_t x = 0; x < dimX; ++x ) { - computeZouHeVelocityWallCell(pop, {x, dimY-1}, { 0,-1}, uLid); + computeMovingWallCell(pop, {x, dimY-1}, {0, -1}, {uLid, 0}); } + + // straight wall cell bounce back for ( std::size_t x = 1; x < dimX-1; ++x ) { computeWallCell(pop, {x, 0}, { 0, 1}); } @@ -76,7 +78,7 @@ int main() { for ( std::size_t t = 0; t <= 10000; ++t ) { computeLbmStep(); - if ( t % 100 == 0 ) { + if ( t % 1000 == 0 ) { std::cout << "."; std::cout.flush(); fluid.writeAsVTK("result/data_t" + std::to_string(t) + ".vtk"); diff --git a/src/boundary_conditions.cc b/src/boundary_conditions.cc index 3488463..08b157c 100644 --- a/src/boundary_conditions.cc +++ b/src/boundary_conditions.cc @@ -1,5 +1,7 @@ #include "boundary_conditions.h" +#include "lbm.h" + std::pair, Vector> neighbors(Vector v) { if ( v[0] == 0 ) { return { @@ -27,6 +29,21 @@ void computeWallCell(DataCellBuffer& pop, Vector cell, Vector pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB); } +void computeMovingWallCell(DataCellBuffer& pop, Vector cell, Vector normal, Vector u) { + const auto [neighborA, neighborB] = neighbors(normal); + + 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).get(-neighborA) + + pop.curr(cell).get(-normal ) + + pop.curr(cell).get(-neighborB) + ); + + pop.curr(cell).get(neighborA) = pop.curr(cell).get(-neighborA) - (6. * weight.get(-neighborA) * rho * (-neighborA * u)); + pop.curr(cell).get(normal ) = pop.curr(cell).get(-normal ) - (6. * weight.get(-normal) * rho * (-normal * u)); + pop.curr(cell).get(neighborB) = pop.curr(cell).get(-neighborB) - (6. * weight.get(-neighborB) * rho * (-neighborB * u)); +} + void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector cell, Vector normal, double vX) { const auto [neighborA, neighborB] = neighbors(normal); diff --git a/src/boundary_conditions.h b/src/boundary_conditions.h index 6b7de19..5f528e1 100644 --- a/src/boundary_conditions.h +++ b/src/boundary_conditions.h @@ -5,4 +5,6 @@ void computeWallCell(DataCellBuffer& pop, Vector cell, Vector normal); +void computeMovingWallCell(DataCellBuffer& pop, Vector cell, Vector normal, Vector u); + void computeZouHeVelocityWallCell(DataCellBuffer& pop, Vector cell, Vector normal, double vX); diff --git a/src/vector.h b/src/vector.h index 5e105c4..839707f 100644 --- a/src/vector.h +++ b/src/vector.h @@ -52,6 +52,11 @@ Vector operator*(T scalar, const Vector& v) { }; } +template +decltype(T{}*W{}) operator*(const Vector& a, const Vector& b) { + return a[0]*b[0] + a[1]*b[1]; +} + template Vector operator-(const Vector& a, const Vector& b) { return Vector{ -- cgit v1.2.3