Skip to content

Re: Simulation of Hydrophobic surfaces

OpenLB – Open Source Lattice Boltzmann Code › Forums › on Lattice Boltzmann Methods › General Topics › Simulation of Hydrophobic surfaces › Re: Simulation of Hydrophobic surfaces

March 3, 2017 at 5:31 pm #2540
shaunellis22
Member

Dear Mathias,

As a starting point I have decided to edit the multiComponent3d such that I can simulate a fluid in pockets on the upper and lower surface of a channel, with a different fluid in the center of a channel (This fluid is the fluid which will experience the force). I have simulated this by renaming the fluid in pockets with a number 3 and leaving the main flow as 1. However I am getting an averageRho for both fluids as ‘nan’, and the output is not as expected. Please could you point me in the right direction at where I am going wrong. Here is the code I am using:

#include “olb3D.h”
#include “olb3D.hh” // use only generic version!
#include <cstdlib>
#include <iostream>

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

typedef double T;
#define DESCRIPTOR ShanChenDynOmegaForcedD3Q19Descriptor

// Parameters for the simulation setup
const int maxIter = 1000;
const T lx1 = 1.0; // length of step
const T ly1 = 0.25; // height of step
const T lx0 = 18.0; // length of channel
const T ly0 = 1.5; // height of channel
const T lz0 = 1.5; // width of channel
const T lx2 = 2.0; // length of gap 1
const T lx3 = 4.0;
const T lx4 = 6.0;
const T lx5 = 8.0;
const T lx6 = 10.0;
const T lx7 = 12.0;
const T lx8 = 14.0;
const T lx9 = 16.0;
const T ly2 = 1.25;
const T lx18 = 1.0;
const T lx10 = 3.0;
const T lx12 = 7.0;
const T lx13 = 9.0;
const T lx14 = 11.0;
const T lx15 = 13.0;
const T lx16 = 15.0;
const T lx17 = 17.0;

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

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

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

///bottom row

Vector<T,3> extend(lx1, ly1, lz0);
Vector<T,3> origin;
IndicatorCuboid3D<T> cuboid2(extend, origin);

Vector<T,3> origin1(lx2, 0, 0);
IndicatorCuboid3D<T> cuboid3(extend, origin1);

Vector<T,3> origin2(lx3, 0, 0);
IndicatorCuboid3D<T> cuboid4(extend, origin2);

Vector<T,3> origin3(lx4, 0, 0);
IndicatorCuboid3D<T> cuboid5(extend, origin3);

Vector<T,3> origin4(lx5, 0, 0);
IndicatorCuboid3D<T> cuboid6(extend, origin4);

Vector<T,3> origin5(lx6, 0, 0);
IndicatorCuboid3D<T> cuboid7(extend, origin5);

Vector<T,3> origin6(lx7, 0, 0);
IndicatorCuboid3D<T> cuboid8(extend, origin6);

Vector<T,3> origin7(lx8, 0, 0);
IndicatorCuboid3D<T> cuboid9(extend, origin7);

Vector<T,3> origin8(lx9, 0, 0);
IndicatorCuboid3D<T> cuboid10(extend, origin8);

///top row

Vector<T,3> origin9(0, ly2, 0) ;
IndicatorCuboid3D<T> cuboid11(extend, origin9);

Vector<T,3> origin10(lx2, ly2, 0);
IndicatorCuboid3D<T> cuboid12(extend, origin10);

Vector<T,3> origin11(lx3, ly2, 0);
IndicatorCuboid3D<T> cuboid13(extend, origin11);

Vector<T,3> origin12(lx4, ly2, 0);
IndicatorCuboid3D<T> cuboid14(extend, origin12);

Vector<T,3> origin13(lx5, ly2, 0);
IndicatorCuboid3D<T> cuboid15(extend, origin13);

Vector<T,3> origin14(lx6, ly2, 0);
IndicatorCuboid3D<T> cuboid16(extend, origin14);

Vector<T,3> origin15(lx7, ly2, 0);
IndicatorCuboid3D<T> cuboid17(extend, origin15);

Vector<T,3> origin16(lx8, ly2, 0);
IndicatorCuboid3D<T> cuboid18(extend, origin16);

Vector<T,3> origin17(lx9, ly2, 0);
IndicatorCuboid3D<T> cuboid19(extend, origin17);

//bottom row gas

Vector<T,3> origin18(lx10, 0, 0);
IndicatorCuboid3D<T> gas(extend, origin18);

Vector<T,3> origin19(lx11, 0, 0);
IndicatorCuboid3D<T> gas1(extend, origin19);

Vector<T,3> origin20(lx12, 0, 0);
IndicatorCuboid3D<T> gas2(extend, origin20);

Vector<T,3> origin21(lx13, 0, 0);
IndicatorCuboid3D<T> gas3(extend, origin21);

Vector<T,3> origin22(lx14, 0, 0);
IndicatorCuboid3D<T> gas4(extend, origin22);

Vector<T,3> origin23(lx15, 0, 0);
IndicatorCuboid3D<T> gas5(extend, origin23);

Vector<T,3> origin24(lx16, 0, 0);
IndicatorCuboid3D<T> gas6(extend, origin24);

Vector<T,3> origin25(lx17, 0, 0);
IndicatorCuboid3D<T> gas7(extend, origin25);

Vector<T,3> origin26(lx18, 0, 0);
IndicatorCuboid3D<T> gas17(extend, origin26);

// Top row gas

Vector<T,3> origin27(lx10, ly2, 0) ;
IndicatorCuboid3D<T> gas8(extend, origin27);

Vector<T,3> origin28(lx11, ly2, 0);
IndicatorCuboid3D<T> gas9(extend, origin28);

Vector<T,3> origin29(lx12, ly2, 0);
IndicatorCuboid3D<T> gas10(extend, origin29);

Vector<T,3> origin30(lx13, ly2, 0);
IndicatorCuboid3D<T> gas11(extend, origin30);

Vector<T,3> origin31(lx14, ly2, 0);
IndicatorCuboid3D<T> gas12(extend, origin31);

Vector<T,3> origin32(lx15, ly2, 0);
IndicatorCuboid3D<T> gas13(extend, origin32);

Vector<T,3> origin33(lx16, ly2, 0);
IndicatorCuboid3D<T> gas14(extend, origin33);

Vector<T,3> origin34(lx17, ly2, 0);
IndicatorCuboid3D<T> gas15(extend, origin34);

Vector<T,3> origin35(lx18, ly2, 0);
IndicatorCuboid3D<T> gas16(extend, origin35);

superGeometry.rename(1,2,cuboid2);
superGeometry.rename(1,2,cuboid3);
superGeometry.rename(1,2,cuboid4);
superGeometry.rename(1,2,cuboid5);
superGeometry.rename(1,2,cuboid6);
superGeometry.rename(1,2,cuboid7);
superGeometry.rename(1,2,cuboid8);
superGeometry.rename(1,2,cuboid9);
superGeometry.rename(1,2,cuboid10);
superGeometry.rename(1,2,cuboid11);
superGeometry.rename(1,2,cuboid12);
superGeometry.rename(1,2,cuboid13);
superGeometry.rename(1,2,cuboid14);
superGeometry.rename(1,2,cuboid15);
superGeometry.rename(1,2,cuboid16);
superGeometry.rename(1,2,cuboid17);
superGeometry.rename(1,2,cuboid18);
superGeometry.rename(1,2,cuboid19);

superGeometry.rename(1,3,gas);
superGeometry.rename(1,3,gas1);
superGeometry.rename(1,3,gas2);
superGeometry.rename(1,3,gas3);
superGeometry.rename(1,3,gas4);
superGeometry.rename(1,3,gas5);
superGeometry.rename(1,3,gas6);
superGeometry.rename(1,3,gas7);
superGeometry.rename(1,3,gas17);
superGeometry.rename(1,3,gas8);
superGeometry.rename(1,3,gas9);
superGeometry.rename(1,3,gas10);
superGeometry.rename(1,3,gas11);
superGeometry.rename(1,3,gas12);
superGeometry.rename(1,3,gas13);
superGeometry.rename(1,3,gas14);
superGeometry.rename(1,3,gas15);
superGeometry.rename(1,3,gas16);

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

superGeometry.print();

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

void prepareLattice(SuperLattice3D<T, DESCRIPTOR>& sLatticeOne,
SuperLattice3D<T, DESCRIPTOR>& sLatticeTwo,
Dynamics<T, DESCRIPTOR>& bulkDynamics1,
Dynamics<T, DESCRIPTOR>& bulkDynamics2,
Dynamics<T, DESCRIPTOR>& bounceBackRho0,
Dynamics<T, DESCRIPTOR>& bounceBackRho1,
SuperGeometry3D<T>& superGeometry) {

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

// The setup is: periodicity along horizontal direction, bounce-back on top
// and bottom. The upper half is initially filled with fluid 1 + random noise,
// and the lower half with fluid 2. Only fluid 1 experiences a forces,
// directed downwards.

/// define lattice Dynamics
sLatticeOne.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());
sLatticeTwo.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, DESCRIPTOR>());

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

sLatticeOne.defineDynamics(superGeometry, 3, &bounceBackRho0);
sLatticeTwo.defineDynamics(superGeometry, 3, &bounceBackRho1);
sLatticeOne.defineDynamics(superGeometry, 4, &bounceBackRho1);
sLatticeTwo.defineDynamics(superGeometry, 4, &bounceBackRho0);

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

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

if(iT==0) {

AnalyticalConst3D<T,T> noise(4.e-2);
std::vector<T> v(3,T());
AnalyticalConst3D<T,T> zeroV(v);
AnalyticalConst3D<T,T> zero(0.);
AnalyticalLinear3D<T,T> one(0.,-force*DESCRIPTOR<T>::invCs2,0.,0.98+force*ly0*DESCRIPTOR<T>::invCs2);
AnalyticalConst3D<T,T> onePlus(0.98+force*ly0/2.*DESCRIPTOR<T>::invCs2);
AnalyticalRandom3D<T,T> random;
AnalyticalIdentity3D<T,T> randomOne(random*noise+one);
AnalyticalIdentity3D<T,T> randomPlus(random*noise+onePlus);
std::vector<T> F(3,T());
F[1] = -force;
AnalyticalConst3D<T,T> f(F);

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

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

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

/*sLatticeOne.defineRhoU(superGeometry, 3, zero, zeroV);
sLatticeOne.iniEquilibrium(superGeometry, 3, zero, zeroV);
sLatticeOne.defineExternalField(superGeometry, 3,
DESCRIPTOR<T>::ExternalField::externalForceBeginsAt,
DESCRIPTOR<T>::ExternalField::sizeOfExternalForce, f );
sLatticeTwo.defineRhoU(superGeometry, 3, one, zeroV);
sLatticeTwo.iniEquilibrium(superGeometry, 3, one, zeroV);

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

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

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

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

const int vtkIter = 500;
const int statIter = 10;

if (iT==0) {
/// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry3D<T, DESCRIPTOR> geometry(sLatticeOne, superGeometry);
SuperLatticeCuboid3D<T, DESCRIPTOR> cuboid(sLatticeOne);
SuperLatticeRank3D<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;
SuperLatticeVelocity3D<T, DESCRIPTOR> velocity(sLatticeOne);
SuperLatticeDensity3D<T, DESCRIPTOR> density(sLatticeOne);
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( density );
vtkWriter.write(iT);

BlockLatticeReduction3D<T, DESCRIPTOR> planeReduction(density, 0, 0, -1 );
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 = 7./(T)ly0/(T)ly0;

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

#ifdef PARALLEL_MODE_MPI
CuboidGeometry3D<T> cGeometry(0, 0, 0, 1, lx0, ly0, lz0, singleton::mpi().getSize());
#else
CuboidGeometry3D<T> cGeometry(0, 0, 0, 1, lx0, ly0, lz0, 1);
#endif

cGeometry.setPeriodicity(true, false, true);

HeuristicLoadBalancer<T> loadBalancer(cGeometry);

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

prepareGeometry(superGeometry);

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

SuperLattice3D<T, DESCRIPTOR> sLatticeOne(superGeometry);
SuperLattice3D<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.001 );

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

sLatticeOne.addLatticeCoupling(superGeometry, 1, coupling, sLatticeTwo);
sLatticeOne.addLatticeCoupling(superGeometry, 3, coupling, sLatticeTwo);
//sLatticeOne.addLatticeCoupling(superGeometry, 3, coupling, sLatticeTwo);
//sLatticeOne.addLatticeCoupling(superGeometry, 4, coupling, sLatticeTwo);

prepareLattice(sLatticeOne, sLatticeTwo, bulkDynamics1, bulkDynamics2,
bounceBackRho0, bounceBackRho1, 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();
//sLatticeTwo.executeCoupling();

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

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