Skip to content

Thermal flow with higher temperature differences without gravity

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Thermal flow with higher temperature differences without gravity

Viewing 2 posts - 1 through 2 (of 2 total)
  • Author
    Posts
  • #7118
    rompak
    Participant

    I am expanding the code available in file RayleighBenard2D to model thermoacoustic engine. Everything works fine, but I have encountered 2 issues:
    – I have to apply big temperature differences e. g. of order 400K, so the Boussinesq approximation probably is not valid in my case,
    – I have to do the simulation without gravity i. e. I have to put this simulation in XY plane when gravity is in Z direction, so that the gravity doesnt interfere with XY behaviour.

    I attach the part of the code responsible for these parts. Do you have any ideas how to solve those problems? I think, that problem lies somewhere in the following lines:
    std::vector<T> dir{0.0, 1.};

    T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
    converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();

    NavierStokesAdvectionDiffusionCouplingGenerator2D<T,NSDESCRIPTOR> coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(ly), boussinesqForcePrefactor, converter.getLatticeTemperature(Tcold), 1., dir);

    Thank you in advance.

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

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

    ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> converter(
    (T) 0.1/N, // physDeltaX
    (T) 0.1 / (1e-5 / 0.1 * util::sqrt( Ra / Pr)) * 0.1 / N, // physDeltaT = charLatticeVelocity / charPhysVelocity * physDeltaX
    (T) 0.1, // charPhysLength
    (T) 1e-5 / 0.1 * util::sqrt( Ra / Pr ), // charPhysVelocity
    (T) 1e-5, // physViscosity
    (T) 1.0, // physDensity
    (T) 0.03, // physThermalConductivity
    (T) Pr * 0.03 / 1e-5 / 1.0, // physSpecificHeatCapacity
    (T) Ra * 1e-5 * 1e-5 / Pr / 9.81 / (Thot – Tcold) / util::pow(0.1, 3), // physThermalExpansionCoefficient
    (T) Tcold, // charPhysLowTemperature
    (T) Thot // charPhysHighTemperature
    );
    converter.print();

    /// === 2nd Step: Prepare Geometry ===
    std::vector<T> extend(2,T());
    extend[0] = lx;
    extend[1] = ly;
    std::vector<T> origin(2,T());
    IndicatorCuboid2D<T> cuboid(extend, origin);

    /// Instantiation of a cuboidGeometry with weights
    #ifdef PARALLEL_MODE_MPI
    const int noOfCuboids = singleton::mpi().getSize();
    #else
    const int noOfCuboids = 7;
    #endif
    CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), noOfCuboids);

    cuboidGeometry.setPeriodicity(false, false);

    HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);

    SuperGeometry<T,2> superGeometry(cuboidGeometry, loadBalancer);

    prepareGeometry(superGeometry, converter);

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

    SuperLattice<T, TDESCRIPTOR> ADlattice(superGeometry);
    SuperLattice<T, NSDESCRIPTOR> NSlattice(superGeometry);

    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
    // This coupling must be necessarily be put on the Navier-Stokes lattice!!
    // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//

    std::vector<T> dir{0.0, 1.};

    T boussinesqForcePrefactor = 9.81 / converter.getConversionFactorVelocity() * converter.getConversionFactorTime() *
    converter.getCharPhysTemperatureDifference() * converter.getPhysThermalExpansionCoefficient();

    NavierStokesAdvectionDiffusionCouplingGenerator2D<T,NSDESCRIPTOR> coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(ly), boussinesqForcePrefactor, converter.getLatticeTemperature(Tcold), 1., dir);

    NSlattice.addLatticeCoupling(coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 2, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 3, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 4, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 5, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 6, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 7, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 8, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 9, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 10, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 11, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 12, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 13, coupling, ADlattice);
    NSlattice.addLatticeCoupling(superGeometry, 14, coupling, ADlattice);

    prepareLattice(converter, NSlattice, ADlattice, superGeometry);

    /// === 4th Step: Main Loop with Timer ===
    util::Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() );
    timer.start();

    util::ValueTracer<T> converge(converter.getLatticeTime(50.),epsilon);
    for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT); ++iT) {

    if (converge.hasConverged()) {
    clout << “Simulation converged.” << std::endl;
    getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());

    clout << “Time ” << iT << “.” << std::endl;

    break;
    }

    /// === 5th Step: Definition of Initial and Boundary Conditions ===
    setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry);

    /// === 6th Step: Collide and Stream Execution ===
    ADlattice.collideAndStream();
    NSlattice.collideAndStream();

    NSlattice.executeCoupling();

    /// === 7th Step: Computation and Output of the Results ===
    getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, converge.hasConverged());
    converge.takeValue(ADlattice.getStatistics().getAverageEnergy(),true);
    }

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

    #7122
    mathias
    Keymaster

    I agree with your statement that Boussinesq approximation is not valid for your case. Other compressible modells we do not have implemented in OpenLB — so that would have do be done first.

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