Skip to content

Simulation of Hydrophobic surfaces

Viewing 6 posts - 1 through 6 (of 6 total)
  • Author
    Posts
  • #1895
    shaunellis22
    Member

    I am new with openlb and I am currently tasked with simulating flow over hydrophobic surfaces. It is easy enough to simulate flow through a cuboid, but how would I simulate gas pockets on the upper and lower surfaces of the cuboid?

    Many thanks,

    #2538
    mathias
    Keymaster

    Dear shaunellis-22,

    Have a look at our 2-phase and 2-component examples in the example section.

    Best
    Mathias

    #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();
    }

    #2544
    mathias
    Keymaster

    Dear shaunellis-22,

    your simulation diverged. Did you get nan straight after a few time steps? It might be that you set some unqhysical boundary conditions or the chosen force is too high. Also too high velocities yield to unstabliities. I would start with the provided example and change one thing after another checking all the time that the simulation is stable.

    Best
    Mathias

    #2549
    shaunellis22
    Member

    Hi Mathias,

    Thank you for your help, I have now simulated the geometry of the sytsem under a force. I have subsequently been trying to simulate poiseuille flow through the geometry rather than the fluid undergoing a force, however have been struggling to do so. Do you have any tips on this?

    Many thanks,

    Shaun

    #2547
    mathias
    Keymaster

    I would start with changing one of the multi-component examples by changing boundary conditions as well as putting new geometry features in.

    Best
    Mathias

Viewing 6 posts - 1 through 6 (of 6 total)
  • You must be logged in to reply to this topic.