aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--channel.cc55
1 files changed, 22 insertions, 33 deletions
diff --git a/channel.cc b/channel.cc
index 837b355..4307dfd 100644
--- a/channel.cc
+++ b/channel.cc
@@ -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");