aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--boltz.cc55
1 files changed, 41 insertions, 14 deletions
diff --git a/boltz.cc b/boltz.cc
index 1acebc2..21f5734 100644
--- a/boltz.cc
+++ b/boltz.cc
@@ -81,7 +81,12 @@ class CellBuffer {
constexpr std::size_t dimX = 128;
constexpr std::size_t dimY = 128;
-constexpr double tau = 0.8;
+constexpr double lidVelocityX = 0.1;
+constexpr double reynolds = 1000;
+constexpr double nu = lidVelocityX * dimX / reynolds;
+constexpr double tau = 3*nu + 0.5;
+
+constexpr double omega = 1. / tau;
CellBuffer pop(dimX, dimY);
@@ -90,8 +95,6 @@ Velocity velocity[dimX][dimY];
Force force [dimX][dimY];
void computeFluidCell(std::size_t x, std::size_t y) {
- //velocity[x][y] += tau / density[x][y] * force[x][y];
-
// compute equilibrium
Cell eq;
eq.equilibrize(density[x][y], velocity[x][y]);
@@ -99,7 +102,7 @@ void computeFluidCell(std::size_t x, std::size_t y) {
// collide (BGK, relax towards equilibrium) & stream
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) + 1./tau * (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) + omega * (eq.get(i,j) - pop.prev(x,y).get(i,j));
}
}
}
@@ -131,19 +134,39 @@ void computeWallCell(std::size_t x, std::size_t y, Vector<int> normal) {
pop.curr(x,y).get(neighborB[0], neighborB[1]) = pop.curr(x-neighborB[0], y-neighborB[1]).get(-neighborB[0], -neighborB[1]);
}
+void computeHorizontalMovingWallCell(std::size_t x, std::size_t y, Vector<int> normal, double vX) {
+ const auto [neighborA, neighborB] = neighbors(normal);
+
+ const double rho = pop.curr(x-neighborA[0], y-neighborA[1]).get(-1,0) + pop.curr(x-normal[0],y-normal[1]).get(0,0) + pop.curr(x-neighborB[0], y-neighborB[1]).get(1,0)
+ + 2.*(
+ pop.curr(x-neighborA[0], y-neighborA[1]).get(-neighborA[0], -neighborA[1]) +
+ pop.curr(x-normal[0] , y-normal[1] ).get(-normal[0] , -normal[1] ) +
+ pop.curr(x-neighborB[0], y-neighborB[1]).get(-neighborB[0], -neighborB[1])
+ );
+ const double hor = pop.curr(x-neighborA[0], y-neighborA[1]).get(-1,0) - pop.curr(x-neighborB[0], y-neighborB[1]).get(1,0);
+
+ pop.curr(x,y).get(neighborA[0], neighborA[1]) = pop.curr(x-neighborA[0], y-neighborA[1]).get(-neighborA[0], -neighborA[1]) + 0.5*hor - 0.5*vX*rho;
+ pop.curr(x,y).get(normal[0] , normal[1] ) = pop.curr(x-normal[0] , y-normal[1] ).get(-normal[0] , -normal[1] );
+ pop.curr(x,y).get(neighborB[0], neighborB[1]) = pop.curr(x-neighborB[0], y-neighborB[1]).get(-neighborB[0], -neighborB[1]) - 0.5*hor + 0.5*vX*rho;
+
+ //pop.curr(x,y).get(neighborA[0], neighborA[1]) = pop.curr(x-neighborA[0], y-neighborA[1]).get(-neighborA[0], -neighborA[1]) - 6 * weight.get( 1, 1) * density[x][y] * vX;
+ //pop.curr(x,y).get(normal[0] , normal[1] ) = pop.curr(x-normal[0] , y-normal[1] ).get(-normal[0] , -normal[1] );
+ //pop.curr(x,y).get(neighborB[0], neighborB[1]) = pop.curr(x-neighborB[0], y-neighborB[1]).get(-neighborB[0], -neighborB[1]) + 6 * weight.get(-1, 1) * density[x][y] * vX;
+}
+
void computeLbmStep(std::size_t t) {
pop.swap();
for ( std::size_t x = 1; x < dimX - 1; ++x ) {
for ( std::size_t y = 1; y < dimY - 1; ++y ) {
- if ( x <= 20 || x >= 40 || y <= 20 || y >= 80 ) {
+ //if ( x <= 20 || x >= 40 || y <= 20 || y >= 80 ) {
computeFluidCell(x, y);
- }
+ //}
}
}
// obstacle
- {
+ /*{
for ( std::size_t x = 21; x < 40; ++x ) {
computeWallCell(x, 80, { 0, 1});
computeWallCell(x, 20, { 0,-1});
@@ -157,12 +180,12 @@ void computeLbmStep(std::size_t t) {
computeWallCell(40, 20, { 1,-1});
computeWallCell(20, 80, {-1, 1});
computeWallCell(20, 20, {-1,-1});
- }
+ }*/
// straight wall cell bounce back
for ( std::size_t x = 2; x < dimX - 2; ++x ) {
- computeWallCell(x, 1, { 0, 1});
- computeWallCell(x, dimY-2,{ 0,-1});
+ computeWallCell(x, 1, { 0, 1});
+ computeHorizontalMovingWallCell(x, dimY-2, { 0,-1}, lidVelocityX);
}
for ( std::size_t y = 2; y < dimY - 2; ++y ) {
computeWallCell(1, y, { 1, 0});
@@ -196,11 +219,11 @@ void init() {
}
}
- for ( std::size_t y = 50; y < dimY-50; ++y ) {
+ /*for ( std::size_t y = 50; y < dimY-50; ++y ) {
for ( std::size_t x = 50; x < dimX-50; ++x ) {
density[x][y] = 2.0;
}
- }
+ }*/
}
void writeCurrentStateAsVTK(int time) {
@@ -248,10 +271,14 @@ void writeCurrentStateAsVTK(int time) {
int main() {
init();
- for ( std::size_t t = 0; t < 1000; ++t ) {
+ std::cout << "u: " << lidVelocityX << std::endl;
+ std::cout << "tau: " << tau << std::endl;
+
+ for ( std::size_t t = 0; t < 20000; ++t ) {
computeLbmStep(t);
- if ( t % 2 == 0 ) {
+ if ( t % 1000 == 0 ) {
+ std::cout << "Writing " << t << "." << std::endl;
writeCurrentStateAsVTK(t);
}
}