Welcome! Log In Create A New Profile


[SOLVED] Can't get rid of a stubborn fixed temperature column of nodes!

Posted by eduardowojeduardowoj  
[SOLVED] Can't get rid of a stubborn fixed temperature column of nodes!
May 18, 2017 01:02PM
EDIT I posted the solution below my original text, and removed the original, bugged code.

Hi guys,

I really don't know how to express my problem better than how's written on the subject. I got sdea's sample advection-diffusion code for the advective transport of a component (http://www.palabos.org/forum/read.php?4,8390) and modified it to consider two transported components - my ultimate goal with this is to replicate Duncan & Toor's experiment with LBM.

Since heat and mass transfer are analogous, where it says "temperature" I'm considering "concentration". My full code is below the text. Sorry for the huge amount of comments, I've been testing a lot of stuff trying to figure out how to solve this problem.

I was able to pinpoint the source of the problem. On the function concentrationSetup (again, concentration is modelled as temperature, due to the mathematical analogy between heat and mass transfer), there's the definition onf the boundary conditions:

Language: C++
Box2D leftwall(0, 0 , 0,ny-1); Box2D rightwall(nx ,nx, 0,ny-1);     // Set the boundary-conditions for lattice1 BC.addTemperatureBoundary1N(leftwall, lattice1); BC.addTemperatureBoundary1P(rightwall, lattice1); // Set the boundary-conditions for lattice2 BC.addTemperatureBoundary1N(leftwall, lattice2); BC.addTemperatureBoundary1P(rightwall, lattice2);

The problem lies within the leftwall boundary condition. With this definition, and initializing the concentration (temperature) of 0.5 (green) for each components at the first and last sixth of the channel and 0 (blue) for the rest (reffer to the full code below), this leftwall column of nodes initializes to 1.0 (dark red):

A bit hard to see, but it's the very first column of nodes, on the left.

When I change leftwall a bit to

Language: C++
Box2D leftwall(10, 10 , 0,ny-1);

the red column shifts its position to the tenth column:

And when I increase its width, the red column gets a bit fatter too:

Language: C++
Box2D leftwall(10, 15, 0, ny-1)

And this column is present also on the snapshots for the concentration of the second component, represented by lattice2:

This is obviously a problem, because this spurious concentration/temperature gets advected with the rest of the fluid, and I don't want this behaviour! For example, this is the following iteration of the program:

This is driving me nuts!!



Here's what I did. I starter from "scratch", from sdea's code. Configured two components as I wanted and I still had the fixed-values column problem. Then I decided to try something a little bit different: I commented the boundary condition definition lines, and replaced these definitions with explicit definitions over the region I wanted. And, voilĂ , problem solved!

However, I am not very comfortable with simply commenting boundary conditions definitions. So this is something I must correct in the future. I set up the advective field with a null velocity, i.e., it's a still fluid, and configured three diffusive components. I then initialized random velocities on each lattice site for these components, to try to mimic the thermodynamic definition of temperature as thermal kinetic energy of the molecules and put this to run, writing the concentration profile of each component (at constant y = ny/2). Below is an animation of the evolution of the concentration profiles, showing that the components do diffuse on the system.

I just need to find how to couple them now.

The code is:

Language: C++
#include "palabos2D.h" #include "palabos2D.hh" #include <iostream> #include <vector> #include <cmath> #include <iostream> #include <fstream> #include <iomanip> #include <string> #include <sstream> #include <cstdio> #include <cstdlib>   using namespace std; using namespace plb;   typedef double T;     #define NSDES descriptors::D2Q9Descriptor #define AVDES descriptors::AdvectionDiffusionD2Q5Descriptor   // ************************************************************************ // Part for the NS part of the code. A classic poiseuille flow is simulated. // Most of the code is taken from Palabos tutorial 1.5 // ************************************************************************   // Poiselle velocity analitical // Velocity on the parabolic Poiseuille profile T poiseuilleVelocity(plint iY, IncomprFlowParam<T> const& parameters) { T y = (T)iY / parameters.getResolution(); return 4.*parameters.getLatticeU() * (y-y*y); }   // A functional, used to initialize the velocity for the boundary conditions template<typename T> class PoiseuilleVelocity { public: PoiseuilleVelocity(IncomprFlowParam<T> parameters_) : parameters(parameters_) { } // This version of the operator returns the velocity only, // to instantiate the boundary condition. void operator()(plint iX, plint iY, Array<T,2>& u) const { u[0] = 0.*poiseuilleVelocity(iY, parameters); u[1] = T(); } // This version of the operator returns also a constant value for // the density, to create the initial condition. void operator()(plint iX, plint iY, T& rho, Array<T,2>& u) const { u[0] = 0.*poiseuilleVelocity(iY, parameters); u[1] = T(); rho = (T)1; } private: IncomprFlowParam<T> parameters; };   // Function to set up the channel void channelSetup ( MultiBlockLattice2D<T,NSDES>& lattice, IncomprFlowParam<T> const& parameters, OnLatticeBoundaryCondition2D<T,NSDES>& boundaryCondition ) { // Create Velocity boundary conditions. boundaryCondition.setVelocityConditionOnBlockBoundaries(lattice);   // Specify the boundary velocity. setBoundaryVelocity ( lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters) );   // Create the initial condition. initializeAtEquilibrium ( lattice, lattice.getBoundingBox(), PoiseuilleVelocity<T>(parameters) );   lattice.initialize(); }   // ************************************************************************ // Advection-diffusion setup ********************************************* // ************************************************************************   void concentrationSetup(MultiBlockLattice2D<T,AVDES> &lattice1, MultiBlockLattice2D<T,AVDES> &lattice2, MultiBlockLattice2D<T,AVDES> &lattice3, OnLatticeAdvectionDiffusionBoundaryCondition2D<T, AVDES> &BC, IncomprFlowParam<T> parameters1, IncomprFlowParam<T> parameters2, IncomprFlowParam<T> parameters3) {   // Get lattice dimensions const plint nx = parameters1.getNx(); const plint ny = parameters1.getNy(); plint xx, yy; T ran1, ran2;   // Set-up the geometry //Box2D topwall(0, nx-1, ny-1, ny-1); //Box2D bottomwall(0, nx-1, 0, 0); Box2D leftwall(0, nx/6, 0,ny-1); Box2D rightwall(nx-nx/6, nx-1, 0,ny-1);   // Set the boundary-conditions for lattice1 //BC.addTemperatureBoundary0P(rightwall, lattice1); //BC.addTemperatureBoundary0N(leftwall, lattice1); // Set the boundary-conditions for lattice2 //BC.addTemperatureBoundary0P(rightwall, lattice2); //BC.addTemperatureBoundary0N(leftwall, lattice2);   // Try to impose a different temperature for lattice1 setBoundaryTemperature(lattice1, rightwall, 0.0); setBoundaryTemperature(lattice1, leftwall, 0.5); // Try to impose a different temperature for lattice2 setBoundaryTemperature(lattice2, rightwall, 0.5); setBoundaryTemperature(lattice2, leftwall, 0.0); // Try to impose a different temperature for lattice3 setBoundaryTemperature(lattice3, rightwall, 0.0); setBoundaryTemperature(lattice3, leftwall, 0.25);       //Box2D domain(0, nx-1, 0, ny-1); // //Array<T,2> u1(1.0,0.); //Array<T,2> u2(-1.0,0.); // //// Init lattice1 //initializeAtEquilibrium(lattice1, domain, 0.0, u1,0.0); //initializeAtEquilibrium(lattice1, leftwall, 0.5,u1,0.5); //// Init lattice2 //initializeAtEquilibrium(lattice2, domain, 0.0, u2,0.0); //initializeAtEquilibrium(lattice2, rightwall, 0.5,u2,0.5);   pcout << "Preparing random velocities on lattices" << endl; for(xx = 0; xx < nx; xx++) { for(yy = 0; yy < ny; yy++) { Box2D domain(xx, xx, yy, yy);   ran1 = 2*((double)random()/RAND_MAX-0.5); ran2 = 2*((double)random()/RAND_MAX-0.5); Array<T,2> u1(ran1,ran2);   // Init lattice1 initializeAtEquilibrium(lattice1, domain, 0.0, u1,0.0); initializeAtEquilibrium(lattice1, leftwall, 0.5,u1,0.5);   // Init lattice2 ran1 = 2*((double)random()/RAND_MAX-0.5); ran2 = 2*((double)random()/RAND_MAX-0.5); Array<T,2> u2(-ran1,-ran2); initializeAtEquilibrium(lattice2, domain, 0.0, u2,0.0); initializeAtEquilibrium(lattice2, rightwall, 0.5,u2,0.5); } }   Box2D domain2(0, nx-1, 0, ny-1); Array<T,2> u3(0.0,0.0); initializeAtEquilibrium(lattice3, domain2, 0.0, u3,0.0); initializeAtEquilibrium(lattice3, leftwall, 0.25,u3,0.25); // Command to init the lattices lattice1.periodicity().toggle(1,true); lattice1.initialize(); lattice2.periodicity().toggle(1,true); lattice2.initialize(); lattice3.periodicity().toggle(1,true); lattice3.initialize(); }   // ************************************************************************ // ************************************************************************ // ************************************************************************       int main(int argc, char *argv[]) { // Init simulation plbInit(&argc, &argv); global::directories().setOutputDir("./tmp/");   // Parameters for the dynamics IncomprFlowParam<T> parameters1( (T) 1e-2, // Reference velocity (the maximum velocity // in the Poiseuille profile) in lattice units. (T) 5000., // Reynolds number 150, // Resolution of the reference length (channel height). 2, // Channel length in dimensionless variables 1. // Channel height in dimensionless variables );   IncomprFlowParam<T> parameters2( (T) 1e-2, (T) 50., 150, 2, 1. );   IncomprFlowParam<T> parameters3( (T) 1e-2, (T) 50., 150, 2, 1. );   stringstream buffer; string filename1, filename2, filename3; plint nx = parameters1.getNx(); const plint maxIter = 30000;     Array<T,2> u2(2.0, 0.0);   // Instantiate the NS lattice MultiBlockLattice2D<T, NSDES> nsLattice(parameters1.getNx(), parameters1.getNy(), new BGKdynamics<T, NSDES>(parameters1.getOmega()));   // Instantiate adv-diff lattice1 MultiBlockLattice2D<T, AVDES> advLattice1(parameters1.getNx(), parameters1.getNy(), new AdvectionDiffusionBGKdynamics<T, AVDES>(1.)); // Instantiate adv-diff lattice2 MultiBlockLattice2D<T, AVDES> advLattice2(parameters1.getNx(), parameters1.getNy(), new AdvectionDiffusionBGKdynamics<T, AVDES>(1.)); // Instantiate adv-diff lattice3 MultiBlockLattice2D<T, AVDES> advLattice3(parameters1.getNx(), parameters1.getNy(), new AdvectionDiffusionBGKdynamics<T, AVDES>(0.7));   // Instantiate the boundary conditions for both OnLatticeBoundaryCondition2D<T,NSDES> *nsBC = createLocalBoundaryCondition2D<T, NSDES>(); OnLatticeAdvectionDiffusionBoundaryCondition2D<T, AVDES> *avBC = createLocalAdvectionDiffusionBoundaryCondition2D<T, AVDES>();   // Set-up channel for NS equations channelSetup(nsLattice, parameters1, *nsBC);   //Set-up concentration field for both lattices concentrationSetup(advLattice1, advLattice2, advLattice3, *avBC, parameters1, parameters2,parameters3);   // Main LBM cycle ImageWriter<T> image("leeloo"); for (plint iT=0; iT <= maxIter; ++iT) {   if(iT%100==0) { pcout << "Navier-Stokes step " << iT << endl; } nsLattice.collideAndStream(); }   // First attempt to couple the two physics latticeToPassiveAdvDiff(nsLattice, advLattice1, advLattice1.getBoundingBox()); latticeToPassiveAdvDiff(nsLattice, advLattice2, advLattice1.getBoundingBox()); latticeToPassiveAdvDiff(nsLattice, advLattice3, advLattice1.getBoundingBox());   for (plint iT=0; iT <= maxIter; ++iT) {   if(iT%100==0) { pcout << "Writing GIF for advection - timestep = " << iT << endl; image.writeScaledGif(createFileName("ADV-conc1_",iT,6), *computeDensity(advLattice1), parameters1.getNx(),parameters1.getNy()); image.writeScaledGif(createFileName("ADV-conc2_",iT,6), *computeDensity(advLattice2), parameters1.getNx(),parameters1.getNy()); image.writeScaledGif(createFileName("ADV-conc3_",iT,6), *computeDensity(advLattice3), parameters1.getNx(),parameters1.getNy());   buffer << "./ascii_output/conc_profile1_"; buffer << setprecision(5) << iT << ".dat"; filename1 = buffer.str(); //pcout << filename << endl; // Debug - Check filename plb_ofstream concProfile1(filename1.c_str()); buffer.str("");   buffer << "./ascii_output/conc_profile2_"; buffer << setprecision(5) << iT << ".dat"; filename2 = buffer.str(); //pcout << filename << endl; // Debug - Check filename plb_ofstream concProfile2(filename2.c_str()); buffer.str("");   buffer << "./ascii_output/conc_profile3_"; buffer << setprecision(5) << iT << ".dat"; filename3 = buffer.str(); //pcout << filename << endl; // Debug - Check filename plb_ofstream concProfile3(filename3.c_str()); buffer.str("");   for(nx = 0; nx < parameters1.getNx(); nx++) {   Box2D concSection(nx, nx, 0.5*parameters1.getNy(), 0.5*parameters1.getNy()); concProfile1 << setprecision(5) << (double)nx/parameters1.getNx() << " " << *computeDensity(advLattice1,concSection) << endl; concProfile2 << setprecision(5) << (double)nx/parameters1.getNx() << " " << *computeDensity(advLattice2,concSection) << endl; concProfile3 << setprecision(5) << (double)nx/parameters1.getNx() << " " << *computeDensity(advLattice3,concSection) << endl; } //image.writeScaledGif(createFileName("vel",iT, 2), *computeVelocityNorm(nsLattice), parameters.getNx(), parameters.getNy()); } advLattice1.collideAndStream(); advLattice2.collideAndStream(); advLattice3.collideAndStream(); }     //pcout << "Hello palabos"; return 0; }

Edited 1 time(s). Last edit at 05/26/2017 01:27PM by eduardowoj.
Sorry, you do not have permission to post/reply in this forum.