Thermal flow with higher temperature differences without gravity

[rompak]

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

[mathias]

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.

[Author]

Posts