[ivan]

Any ideas of how to do it, Mathias?

[ivan]

.

[ivan]

Okay, so here is one of the combination I’ve been trying recently with which I simply get a lot of “nan” results. The bold lines are the ones related to the pressure boundary condition I am trying to include and without them the simulation works.

Code:

#include “olb2D.h”
#include “olb2D.hh”
#include <cstdlib>
#include <iostream>

using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace std;

typedef double T;
#define DESCRIPTOR ShanChenDynOmegaForcedD2Q9Descriptor

// Parameters for the simulation setup
const int nx = 100;
const int ny = 250;
const T centerDropletX = (nx-1)/2;
const T centerDropletY = 20.;
const T radiusDroplet = 10.;
const T wallThickness =4.;
const int maxIter = 6900;

/// Stores geometry information in form of material numbers
void prepareGeometry(SuperGeometry2D<T>& superGeometry) {

OstreamManager clout(std::cout,”prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;

// Sets material number for fluid and boundary
superGeometry.rename(0,1);

std::vector<T> origin1(2,T());
origin1[0] = 0;
origin1[1] = 0;

std::vector<T> origin2(2,T());
origin2[0] = 0;
origin2[1] = ny-1;

std::vector<T> origin3(2,T());
origin3[0] = 0;
origin3[1] = wallThickness;

std::vector<T> origin4(2,T());
origin4[0] = nx-wallThickness-1;
origin4[1] = wallThickness;

std::vector<T> extend1(2,T());
extend1[0] = nx;
extend1[1] = wallThickness;

std::vector<T> extend2(2,T());
extend2[0] = wallThickness;
extend2[1] = ny-wallThickness;

std::vector<T>extend3(2,T());
extend3[0] = nx;
extend3[1] = 1;

Vector<T,2> centerDroplet(centerDropletX,centerDropletY);

IndicatorCuboid2D<T> bottom(extend1, origin1);
IndicatorCuboid2D<T> top(extend3, origin2);
IndicatorCuboid2D<T> leftwall(extend2, origin3);
IndicatorCuboid2D<T> rightwall(extend2, origin4);
IndicatorCircle2D<T> droplet(centerDroplet, radiusDroplet);

superGeometry.rename(1,2,droplet);
superGeometry.rename(1,3,bottom);
superGeometry.rename(1,3,leftwall);
superGeometry.rename(1,3,rightwall);
superGeometry.rename(1,4,top);

/// Removes all not needed boundary voxels inside the surface
superGeometry.innerClean();
superGeometry.checkForErrors();

superGeometry.print();

clout << “Prepare Geometry … OK” << std::endl;
}

void prepareLattice(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
Dynamics<T, DESCRIPTOR>& bulkDynamics1,
Dynamics<T, DESCRIPTOR>& bulkDynamics2,
Dynamics<T, DESCRIPTOR>& bounceBackRho0,
Dynamics<T, DESCRIPTOR>& bounceBackRho1,
[b] sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeOne,
sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeOne,
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR>& sBoundaryConditionLatticeTwo,
sOffLatticeBoundaryCondition2D<T,DESCRIPTOR>& offBcLatticeTwo,[/b]
SuperGeometry2D<T>& superGeometry)
{

OstreamManager clout(std::cout,”prepareLattice”);
clout << “Prepare Lattice …” << std::endl;

[b] const T omega1 = 1.0;
const T omega2 = 1.0;[/b]

/// define lattice Dynamics
sLatticeOne.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
sLatticeOne.defineDynamics(superGeometry, 1, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 2, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 3, &bounceBackRho1);
sLatticeOne.defineDynamics(superGeometry, 4, &bulkDynamics1);

[b] sBoundaryConditionLatticeOne.addPressureBoundary(superGeometry, 4, omega1);[/b]

sLatticeTwo.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
sLatticeTwo.defineDynamics(superGeometry, 1, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 2, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 3, &bounceBackRho1);
sLatticeTwo.defineDynamics(superGeometry, 4, &bulkDynamics2);
[b]
sBoundaryConditionLatticeTwo.addPressureBoundary(superGeometry, 4, omega2);[/b]

clout << “Prepare Lattice … OK” << std::endl;
}

void setBoundaryValues(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
T force, int iT, SuperGeometry2D<T>& superGeometry) {

if(iT==0) {

AnalyticalConst2D<T,T> noise(4.e-2);
std::vector<T> v(2,T());
AnalyticalConst2D<T,T> zeroV(v);
AnalyticalConst2D<T,T> zero(0.);
AnalyticalLinear2D<T,T> one(0.,-force*DESCRIPTOR<T>::invCs2,0.98+force*ny*DESCRIPTOR<T>::invCs2);
AnalyticalConst2D<T,T> onePlus(0.98+force*ny/2.*DESCRIPTOR<T>::invCs2);
AnalyticalRandom2D<T,T> random;
AnalyticalIdentity2D<T,T> randomOne(random*noise+one);
AnalyticalIdentity2D<T,T> randomPlus(random*noise+onePlus);
std::vector<T> F(2,T());
F[1] = -force;
AnalyticalConst2D<T,T> f(F);

/// for each material set the defineRhou and the Equilibrium

sLatticeOne.defineRhoU(superGeometry, 1, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 1, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 1,
DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

sLatticeOne.defineRhoU(superGeometry, 2, randomOne, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 2, randomOne, zeroV);
sLatticeOne.defineExternalField(superGeometry, 2,
DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

sLatticeOne.defineRhoU(superGeometry, 4, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 4, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 4,
DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

sLatticeTwo.defineRhoU(superGeometry, 1, randomPlus, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 1, randomPlus, zeroV);

sLatticeTwo.defineRhoU(superGeometry, 2, zero, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 2, zero, zeroV);

sLatticeTwo.defineRhoU(superGeometry, 4, randomPlus, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 4, randomPlus, zeroV);

/// Make the lattice ready for simulation
sLatticeOne.initialize();
sLatticeTwo.initialize();
}
}

void getResults(SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
SuperLattice2D<T, DESCRIPTOR>& sLatticeOne, int iT,
SuperGeometry2D<T>& superGeometry, Timer<T>& timer) {

OstreamManager clout(std::cout,”getResults”);
SuperVTKwriter2D<T> vtkWriter(“rayleighTaylor2dsLatticeOne”);

const int vtkIter = 100;
const int statIter = 10;

if (iT==0) {
/// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry(sLatticeOne, superGeometry);
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid(sLatticeOne);
SuperLatticeRank2D<T, DESCRIPTOR> rank(sLatticeOne);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);
vtkWriter.createMasterFile();
}

// Get statistics
if (iT%statIter==0 && iT > 0) {
// Timer console output
timer.update( iT );
timer.printStep();

clout << “averageRhoFluidOne=” << sLatticeOne.getStatistics().getAverageRho();
clout << “; averageRhoFluidTwo=” << sLatticeTwo.getStatistics().getAverageRho() << std::endl;
}

/// Writes the VTK files
if (iT%vtkIter==0) {
clout << “Writing VTK …” << std::endl;
SuperLatticeVelocity2D<T, DESCRIPTOR> velocity(sLatticeOne);
SuperLatticeDensity2D<T, DESCRIPTOR> density(sLatticeOne);
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( density );
vtkWriter.write(iT);

BlockLatticeReduction2D<T, DESCRIPTOR> planeReduction(density);
BlockGifWriter<T> gifWriter;
gifWriter.write(planeReduction, iT, “density”);

clout << “Writing VTK … OK” << std::endl;
}
}

int main(int argc, char *argv[]) {

/// === 1st Step: Initialization ===

olbInit(&argc, &argv);
singleton::directories().setOutputDir(“./tmp/”);
OstreamManager clout(std::cout,”main”);

const T omega1 = 1.0;
const T omega2 = 1.0;
const T G = 3.;
T force = -30./(T)ny/(T)ny;

/// === 2nd Step: Prepare Geometry ===
/// Instantiation of a cuboidGeometry with weights

#ifdef PARALLEL_MODE_MPI
CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, singleton::mpi().getSize());
#else
CuboidGeometry2D<T> cGeometry(0, 0, 1, nx, ny, 1);
#endif

cGeometry.setPeriodicity(false, false);

HeuristicLoadBalancer<T> loadBalancer(cGeometry);

SuperGeometry2D<T> superGeometry(cGeometry, loadBalancer, 2);

prepareGeometry(superGeometry);

/// === 3rd Step: Prepare Lattice ===

SuperLattice2D<T, DESCRIPTOR> sLatticeOne(superGeometry);
SuperLattice2D<T, DESCRIPTOR> sLatticeTwo(superGeometry);

ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics1 (
omega1, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );
ForcedBGKdynamics<T, DESCRIPTOR> bulkDynamics2 (
omega2, instances::getExternalVelocityMomenta<T,DESCRIPTOR>() );

// A bounce-back node with fictitious density 1,
// which is experienced by the partner fluid
BounceBack<T, DESCRIPTOR> bounceBackRho1( 1. );
// A bounce-back node with fictitious density 0,
// which is experienced by the partner fluid
BounceBack<T, DESCRIPTOR> bounceBackRho0( 0. );

std::vector<T> rho0;
rho0.push_back(1);
rho0.push_back(1);
PsiEqualsRho<T,T> interactionPotential;
ShanChenForcedGenerator2D<T,DESCRIPTOR> coupling(G,rho0,interactionPotential);

[b] // Pressure boundary
sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeOne(sLatticeOne);
createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeOne);
sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeOne(sLatticeOne);

sOnLatticeBoundaryCondition2D<T,DESCRIPTOR> sBoundaryConditionLatticeTwo(sLatticeTwo);
createLocalBoundaryCondition2D<T,DESCRIPTOR>(sBoundaryConditionLatticeTwo);
sOffLatticeBoundaryCondition2D<T, DESCRIPTOR> sOffBoundaryConditionLatticeTwo(sLatticeTwo);[/b]

sLatticeOne.addLatticeCoupling(superGeometry, 1, coupling, sLatticeTwo);
sLatticeOne.addLatticeCoupling(superGeometry, 2, coupling, sLatticeTwo);

prepareLattice(sLatticeOne, sLatticeTwo, bulkDynamics1, bulkDynamics2,
bounceBackRho0, bounceBackRho1,[b] sBoundaryConditionLatticeOne, sOffBoundaryConditionLatticeOne, sBoundaryConditionLatticeTwo, sOffBoundaryConditionLatticeTwo,[/b] superGeometry);

/// === 4th Step: Main Loop with Timer ===
int iT = 0;
clout << “starting simulation…” << endl;
Timer<T> timer( maxIter, superGeometry.getStatistics().getNvoxel() );
timer.start();

for (iT=0; iT<maxIter; ++iT) {

/// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues(sLatticeOne, sLatticeTwo, force, iT, superGeometry);

/// === 6th Step: Collide and Stream Execution ===
sLatticeOne.collideAndStream();
sLatticeTwo.collideAndStream();

sLatticeOne.communicate();
sLatticeTwo.communicate();

sLatticeOne.executeCoupling();

/// === 7th Step: Computation and Output of the Results ===
getResults(sLatticeTwo, sLatticeOne, iT, superGeometry, timer);
}

timer.stop();
timer.printSummary();
}

[ivan]

Okay, I have just noticed that is has to do with this part of the code

Code:

void setBoundaryValues(SuperLattice2D<T, DESCRIPTOR>& sLatticeOne,
SuperLattice2D<T, DESCRIPTOR>& sLatticeTwo,
T force, int iT, SuperGeometry2D<T>& superGeometry) {

if(iT==0) {

AnalyticalConst2D<T,T> noise(4.e-2);
std::vector<T> v(2,T());
AnalyticalConst2D<T,T> zeroV(v);
AnalyticalConst2D<T,T> zero(0.);
AnalyticalLinear2D<T,T> one(0.,-force*DESCRIPTOR<T>::invCs2,0.98+force*ny*DESCRIPTOR<T>::invCs2);
AnalyticalConst2D<T,T> onePlus(0.98+force*ny/2.*DESCRIPTOR<T>::invCs2);
AnalyticalRandom2D<T,T> random;
AnalyticalIdentity2D<T,T> randomOne(random*noise+one);
AnalyticalIdentity2D<T,T> randomPlus(random*noise+onePlus);
std::vector<T> F(2,T());
F[1] = -force;
AnalyticalConst2D<T,T> f(F);

/// for each material set the defineRhou and the Equilibrium

sLatticeOne.defineRhoU(superGeometry, 1, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 1, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 1,DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

sLatticeOne.defineRhoU(superGeometry, 2, randomOne, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 2, randomOne, zeroV);
sLatticeOne.defineExternalField(superGeometry, 2, DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

sLatticeTwo.defineRhoU(superGeometry, 1, randomPlus, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 1, randomPlus, zeroV);

sLatticeTwo.defineRhoU(superGeometry, 2, zero, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 2, zero, zeroV);

Which is atributing the densities and equilibrium. So, what I tried to do is to copy the exact same code for material number 4, once I want it to behave exactly like number 1, so:

Code:

/// for each material set the defineRhou and the Equilibrium

sLatticeOne.defineRhoU(superGeometry, 1, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 1, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 1, DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

[b] sLatticeOne.defineRhoU(superGeometry, 4, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 4, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 4, DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );[/b]

sLatticeOne.defineRhoU(superGeometry, 2, randomOne, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 2, randomOne, zeroV);
sLatticeOne.defineExternalField(superGeometry, 2, DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );

sLatticeTwo.defineRhoU(superGeometry, 1, randomPlus, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 1, randomPlus, zeroV);
[b]
sLatticeTwo.defineRhoU(superGeometry, 4, randomPlus, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 4, randomPlus, zeroV);[/b]

sLatticeTwo.defineRhoU(superGeometry, 2, zero, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 2, zero, zeroV);

And right now I do get the static fluid on material 4 BUT for some reason the droplet of the lighter immiscible fluid (material number 2) is unable to flow through it, as if it was a barrier, or a discontinuity. Any idea of how to make the cuboid with material number 4 a continuous part of the static fluid of the bigger cuboid of material number one?

[ivan]

Thank you, Mathias,

As I am simulating something closer to a recipient, there would be no inlet, just an open air top with a pressure boundary condition. And yes, I have been basing myself in the cylinder example, but the main problem is I do not know which descriptors should be included along the addPressure function in order to make it work properly. Besides that, there are some other obstacles before it.

For example, that’s how I’ve set my Geometry

Code:

// Sets material number for fluid and boundary
superGeometry.rename(0,1);

std::vector<T> origin1(2,T());
origin1[0] = 0;
origin1[1] = 0;

std::vector<T> origin2(2,T());
origin2[0] = 0;
origin2[1] = ny-wallThickness;

std::vector<T> origin3(2,T());
origin3[0] = 0;
origin3[1] = wallThickness;

std::vector<T> origin4(2,T());
origin4[0] = nx-wallThickness-1;
origin4[1] = wallThickness;

std::vector<T> extend1(2,T());
extend1[0] = nx;
extend1[1] = wallThickness;

std::vector<T> extend2(2,T());
extend2[0] = wallThickness;
extend2[1] = ny-wallThickness;

Vector<T,2> centerDroplet(centerDropletX,centerDropletY);

IndicatorCuboid2D<T> bottom(extend1, origin1);
IndicatorCuboid2D<T> top(extend1, origin2);
IndicatorCuboid2D<T> leftwall(extend2, origin3);
IndicatorCuboid2D<T> rightwall(extend2, origin4);
IndicatorCircle2D<T> droplet(centerDroplet, radiusDroplet);

superGeometry.rename(1,2,droplet);
superGeometry.rename(1,3,bottom);
superGeometry.rename(1,3,leftwall);
superGeometry.rename(1,3,rightwall);
superGeometry.rename(1,4,top);

And that’s the dynamics of the material numbers:

Code:

/// define lattice Dynamics
sLatticeOne.defineDynamics(superGeometry, 0, &instances::getBounceBack<T, DESCRIPTOR>());
sLatticeOne.defineDynamics(superGeometry, 1, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 2, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 3, &instances::getBounceBack<T, DESCRIPTOR>());
sLatticeOne.defineDynamics(superGeometry, 4, &instances::getBounceBack<T, DESCRIPTOR>());

sLatticeTwo.defineDynamics(superGeometry, 0, &instances::getBounceBack<T, DESCRIPTOR>());
sLatticeTwo.defineDynamics(superGeometry, 1, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 2, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 3, &instances::getBounceBack<T, DESCRIPTOR>());
sLatticeTwo.defineDynamics(superGeometry, 4, &instances::getBounceBack<T, DESCRIPTOR>());

So basically number 1 is my static fluid, number 2 is the lighter immiscible fluid, number 3 are the walls and number 4 is the open air top, which for now is randomly set with bounce back as well. First thing I should do in order to include the pressure boundary condition is to set material number 4 as the static fluid as well. The problem is for some reason, if I switch the previous code to the following one

Code:

/// define lattice Dynamics
sLatticeOne.defineDynamics(superGeometry, 0, &instances::getBounceBack<T, DESCRIPTOR>());
sLatticeOne.defineDynamics(superGeometry, 1, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 2, &bulkDynamics1);
sLatticeOne.defineDynamics(superGeometry, 3, &instances::getBounceBack<T, DESCRIPTOR>());
[b] sLatticeOne.defineDynamics(superGeometry, 4, &bulkDynamics1);[/b]

sLatticeTwo.defineDynamics(superGeometry, 0, &instances::getBounceBack<T, DESCRIPTOR>());
sLatticeTwo.defineDynamics(superGeometry, 1, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 2, &bulkDynamics2);
sLatticeTwo.defineDynamics(superGeometry, 3, &instances::getBounceBack<T, DESCRIPTOR>());
[b] sLatticeTwo.defineDynamics(superGeometry, 4, &bulkDynamics2);[/b]

The material number 4 is set as the lighter fluid as well. How can I set it as the static fluid without using the same material number?

[ivan]

Hello, everybody,rnrnCurrently I am working with the example Multicomponent2D and by changing some geometry and force aspects I was able to turn it into a simulation of a rising bubble of a lighter fluid in a heavier one, as explained in the topic “”Simulation of a Raising Droplet in a Channel- Help””, and now I would like to convert the whole simulation into physical units. My first thought was to implement an LBConverter class, so I got the archive “”cylinder2d.dat”” to have a look and a doubt came up:rnrn

Code:

LBconverter informationrnrncharacteristical valuesrn———————————————————————-rnDimension(d): dim=2rnCharacteristical length(m): charL=0.1rnCharacteristical speed(m/s): charU=0.2rnCharacteristical time(s): charT=0.5rnDensity factor(kg/m^d): charRho=1rnCharacterestical mass(kg): charMass=0.01rnCharacterestical force(N): charForce=0.004rnCharacterestical pressure(Pa): charPressure=0.04rnPressure level(Pa): pressureLevel=0rnPhys. kinematic viscosity(m^2/s): charNu=0.001rnPhys. dynamic viscosity(N*s/m^2): dynVisco=0.001rn======================================================================rnrnlattice valuesrn———————————————————————-rnDeltaX: deltaX=0.05rnLattice velocity: latticeU=0.02rnDeltaT: deltaT=0.001rnReynolds number: Re=20rnDimlessNu: dNu=0.05rnViscosity for computation: latticeNu=0.02rnRelaxation time: tau=0.56rnRelaxation frequency: omega=1.78571rn======================================================================rnrnconversion factorsrn———————————————————————-rnlatticeL(m): latticeL=0.005rnTime step (s): physTime=0.0005rnVelocity factor(m/s): physVelocity=10rnFlowRate factor(m^d/s): physFlowRate=0.05rnMass factor(kg): physMass=2.5e-05rnForce factor(N): physForce=0.5rnForce factor massless(N/kg): physMasslessForce=20000rnPressure factor(Pa): physPressure=100rnlatticePressure: latticeP=0.01

rnrnCan anybody explain why is latticeNu calculated this way?rnrnlatticeNu=charNu*latticeU/charU/latticeLrnrnBesides that, I was wondering how could I implement the LBconverter once in my example I have two different fluid? Would I have to implement two different ones? If so, how?rnrnrnBest regards,rnrnIvan

[ivan]

Thank you so much, Mathias 🙂

[ivan]

Hello, everybody,rnrnInitially, I just wanted to change its length, for which I thought of editing the number of horizontal nodes nx.rnrn

Code:

[quote]// Parameters for the simulation setuprnconst int nx = 400;rnconst int ny = 200;rnconst int maxIter = 20000;[/quote]

rnrnOnce nx was changed from 400 to 30 I got a simulation in which the fluids don’t move, so I decided to analyse the material numbers to understand the reason of it. Above there is the original code and a few pictures for showing what I understood from it so far (which I assume is wrong).rnrn

Code:

[quote]rn // Sets material number for fluid and boundaryrn superGeometry.rename(0,1);rnrn std::vector<T> origin1(2,T());rn origin1[0] = -2.;rn origin1[1] = -2.;rn std::vector<T> origin2(2,T());rn origin2[0] = -2.;rn origin2[1] = ny/2.;rn std::vector<T> origin3(2,T());rn origin3[0] = -2.;rn origin3[1] = ny-1.;rnrn std::vector<T> extend1(2,T());rn extend1[0] = nx+3.;rn extend1[1] = 2.;rn std::vector<T> extend2(2,T());rn extend2[0] = nx+3.;rn extend2[1] = ny/2.+2.;rnrn IndicatorCuboid2D<T> bottom(extend1, origin1);rn IndicatorCuboid2D<T> upper(extend2, origin2);rn IndicatorCuboid2D<T> top(extend1, origin3);rnrn superGeometry.rename(1,2,upper);rn superGeometry.rename(1,3,bottom);rn superGeometry.rename(2,4,top);[/quote]

rnrnObservation: Instead of 400×200, I used 40×20 cells for the representation just to make it simpler, but I am still talking about the original dimensions. rnrnInitially, we have 40×20 nodes, represented by the yellow cells, which material numbers were renamed from 0 to 1rnrn[img]https://postimg.org/image/hlzuqgiip/[/img]rnrnIt was inserted a Cuboid called Bottom, represented by the orange cellsrnrn[img]https://postimg.org/image/lweim1nlt/[/img]rnrnThe Upper Cuboid was inserted and its represented by the green cellsrnrn[img]https://postimg.org/image/so4xowcld/[/img]rnrnThe Top Cuboid was inserted and its represented by the blue cellsrnrn[img]https://postimg.org/image/60poiqx1d/[/img]rnrnThe material numbers were changed according to the following imagesrnrn[img]https://postimg.org/image/md4bm80qp/[/img]rnrn[img]https://postimg.org/image/sfbycpp6p/[/img]rnrn[img]https://postimg.org/image/5rwp6k9mp/[/img]rnrnrnIf I understood correctly, changing nx from 400 to 30 would not change the interactions between the fluids/the simulation itself, since the horizontal dimension of the the cuboids are given in function of it (nx). Then why doesn’t it work?rnrnExtra questionsrnrnWhen a 400×200 nodes Lattice is set, what does it mean to have something like “”origin1[0]=-2″”? Does it mean this node is two nodes left distant the first one (which I assume is on the left bottom of the lattice)?rnrnIf so, why was used -2 instead of -1 or 0? And why nx+3 instead of nx-origin (in this example, nx+2)?Why do we need those extra nodes before and after the lattice which don not seem to be part of it?rnrnI hope I have made myself clear enough ( I assume there must be a couple mistakes with the vocabulary ) for you guys to understand my point.rnrnBest regards,rnrnIvanrnrnrn

[ivan]

edited

[ivan]

Hello, Alejandro,rnrnSo, basically I’ve changed line 156 fromrnrn

Code:

AnalyticalConst2D<T,T> rhoF(1);

to

Code:

AnalyticalConst2D<T,T> rhoF(1000);

rnrnAnd the converter object fromrnrn

Code:

LBconverter<T> converter(rn (int) 2, // dimrn (T) L, // latticeL_rn (T) 0.02/M, // latticeU_rn (T) 0.2*2.*radiusCylinder/Re, // charNu_rn (T) 2*radiusCylinder, // charL_ = 1rn (T) 0.2 // charU_ = 1rn );rn converter.print();

rnrntornrn

Code:

LBconverter<T> converter(rn (int) 2, // dimrn (T) L, // latticeL_rn (T) 0.02/M, // latticeU_rn (T) 0.000001, // charNu_rn (T) 2*radiusCylinder, // charL_ = 1rn (T) Re*0.000001/(2*radiusCylinder), // charU_ = 1rn (T) 1000. // charRhorn );rn converter.print();

rnrnAnd the images are not being generated.rnrnAtt.rnrn

[ivan]

Hello, Alejandro,rnrnThank you one more time for the reply.rnrnJust let me ask you, did you mean the density of the water (CharRho) to be 1000kg/m³? We’re talking about liquid, right? I this case, would the converter object be this one?rnrn

Code:

LBconverter<T> converter(rn (int) 2, // dimrn (T) L, // latticeL_rn (T) 0.02/M, // latticeU_rn (T) 0.000001, // charNu_rn (T) 2*radiusCylinder, // charL_ = 1rn (T) Re*0.000001/(2*radiusCylinder), // charU_ = 1rn (T) 1000. // charRho : characteristical density of the waterrn );

rnrnrnIf so, how does the cpp recognize the “”(T) 1000.”” line as related to the density? Still, would I have to do the same alterations in the line 156 as explained in the very first post?rnrnBest,rnrnIvan

[ivan]

Hello, Alejandro,rnrnFirst of all, thank you for the fast reply.rnrnBasically, what I meant with “”messed up images”” is this:rnrnrnrnWhich is what happens when I simply switch rho(1) to rho(1000), as previously explained.rnrnAnd, as Im still a beginner, can you make it a little clearer? How do I set the class LBconverter<T>? I mean, the line 303 I mentioned is right here:rnrn

Code:

LBconverter<T> converter(rn (int) 2, // dimrn (T) L, // latticeL_rn (T) 0.02/M, // latticeU_rn (T) 0.2*2.*radiusCylinder/Re, // charNu_rn (T) 2*radiusCylinder, // charL_ = 1rn (T) 0.2 // charU_ = 1rn );rn converter.print();

rnrnIs that what you’re talking about? If so, the doubts remain.rnrnAtt.rnrn

[ivan]

Oh, it did! Thank you very much, Robin!

[ivan]

Hello, Alejandro,rnrnFirst of all, thank you so much for replying and I am sorry for taking long to reply you back, I did not have a look at the forum in the last days. Concerning the problem, in the example “”Cylinder2d””, we have a channel in which there is a flow and a round obstacle/circle/cylinder. The idea was to use two different lattices, as in the example “”Multicomponent2d””. In the first one, I would keep the properties of the original fluid and set the material of the obstacle as a second fluid of density zero; in the second one, the obstacle – which now behaves as a fluid – would receive a different density and the original fluid, zero. I am not sure if that would work but as far as I understood I would have a drop of a second fluid in the main flow ( which is no really the original idea of having a raising drop but also works for my research ). What do you think of that?rnrnBest,rnrnIvan

[ivan]

Yes, I have also been having a look at it for a while, but it is still a little difficult for me to understand all of it. rnrnThanks for the tip, anywayrnrnAtt.

[Author]

Posts