Gallium Melting from 2D to 3D

[Jijo]

Hello everyone,

I want to change the galliummtling2d example to 3d and I want to ask if the class (TotalEnthalpyPhaseChangeCouplingGenerator2D) is available in 3D? Also, can I have some steps regarding this conversion.

Your reply is very appreciated.

[ramirofreile]

Hello Anas,

There is a coupling for the total enthalpy in 3D. I suggest you always look for classes in the Doxygen search engine:

https://www.openlb.net/DoxyGen/html/index.html

The conversion should be similar than how you would convert some of the examples in the laminar folder, e.g. cavity2D->3D, poiseuille 2D->3D.

Ramiro Freile

[Jijo]

Hey Ramiro Freile,

Thank you for your help. I have done the conversion but the output text only shows 1 material number rather than 4. Can you please check out my code? I would appreciate that.

The code:
#include “olb3D.h”
#include “olb3D.hh”

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

using T = double;

using NSDESCRIPTOR = D3Q19<POROSITY,VELOCITY_SOLID,FORCE>;
using TDESCRIPTOR = D3Q7<VELOCITY,TEMPERATURE>;

using TotalEnthalpyAdvectionDiffusionDynamics = TotalEnthalpyAdvectionDiffusionTRTdynamics<T,TDESCRIPTOR>;

// Parameters for the simulation setup
const T lx = 88.9e-3; // length of the channel
const T ly = 63.5e-3; // height of the channel
const T lz = 38.1e-3;
int N = 128; // resolution of the model
T tau = 0.51; // relaxation time
const T Ra = 1.656e6; // Rayleigh number
const T Pr = 0.0216; // Prandtl number
const T Ste = 0.039; // Stephan number
const T maxPhysT = 1140.; // simulation time

const T Tcold = 0;
const T Tmelt = (302.8 – 301.3)/(311.0 – 301.3);
const T Thot = 1;

const T lambda_s = 40.6; // W / m K
const T lambda_l = 25.3; // W / m K
const T R_lambda = lambda_s/lambda_l;

const T cp_s = 1.0; // J / kg K
const T cp_l = 0.96; // J / kg K
const T R_cp = cp_s/cp_l;

//for this case, the harmonic mean (cp_ref) is applicable
const T cp_ref = 2.0 * cp_s * cp_l / (cp_s + cp_l); // J / kg K

const T R_alpha = lambda_s / lambda_l * cp_l / cp_s;

const T L = cp_l * (Thot – Tmelt) / Ste; // J / kg

T lattice_Hcold, lattice_Hhot;
T physDeltaX, physDeltaT;

/// Stores geometry information in form of material numbers
void prepareGeometry(SuperGeometry<T,3>& superGeometry,
ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter)
{

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

superGeometry.rename(0,4);

std::vector<T> extend(3,T());
extend[0] = lx;
extend[1] = ly;
extend[2] = lz;

std::vector<T> origin(3,T());
origin[0] = converter.getPhysLength(1);
origin[1] = 0.5*converter.getPhysLength(1);
origin[2] = 0.0;
IndicatorCuboid2D<T> cuboid2(extend, origin);

std::vector<T> extendwallleft(3,T(0));
extendwallleft[0] = converter.getPhysLength(1);
extendwallleft[1] = ly;
extendwallleft[2] = 0.1;

std::vector<T> originwallleft(3,T(0));
originwallleft[0] = 0.0;
originwallleft[1] = 0.0;
originwallleft[2] = 0.0;
IndicatorCuboid3D<T> wallleft(extendwallleft, originwallleft);

std::vector<T> extendwallright(3,T(0));
extendwallright[0] = converter.getPhysLength(1);
extendwallright[1] = ly;
extendwallright[2] = 0.1;

std::vector<T> originwallright(3,T(0));
originwallright[0] = lx+converter.getPhysLength(1);
originwallright[1] = 0.0;
originwallright[2] = 0.0;

IndicatorCuboid3D<T> wallright(extendwallright, originwallright);

superGeometry.rename(4,2,1,wallleft);
superGeometry.rename(4,3,1,wallright);

/// 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( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
SuperGeometry<T,3>& superGeometry )
{

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

T omega = converter.getLatticeRelaxationFrequency();
T Tomega = converter.getLatticeThermalRelaxationFrequency();

NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);
NSlattice.defineDynamics<ForcedPSMBGKdynamics>(superGeometry.getMaterialIndicator({1, 2, 3, 4}));

ADlattice.defineDynamics<NoDynamics>(superGeometry, 0);
ADlattice.defineDynamics<TotalEnthalpyAdvectionDiffusionDynamics>(superGeometry.getMaterialIndicator({1, 2, 3}));
ADlattice.defineDynamics<BounceBack>(superGeometry, 4);

/// sets boundary
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry.getMaterialIndicator({2, 3}));
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({2, 3}));

NSlattice.setParameter<descriptors::OMEGA>(omega);

ADlattice.setParameter<descriptors::OMEGA>(Tomega);
ADlattice.setParameter<collision::TRT::MAGIC>(0.25);
ADlattice.setParameter<TotalEnthalpy::T_S>(Tmelt);
ADlattice.setParameter<TotalEnthalpy::T_L>(Tmelt);
ADlattice.setParameter<TotalEnthalpy::CP_S>(cp_s);
ADlattice.setParameter<TotalEnthalpy::CP_L>(cp_l);
ADlattice.setParameter<TotalEnthalpy::LAMBDA_S>(cp_ref / descriptors::invCs2<T,TDESCRIPTOR>() * (converter.getLatticeThermalRelaxationTime() – 0.5) * R_lambda);
ADlattice.setParameter<TotalEnthalpy::LAMBDA_L>(cp_ref / descriptors::invCs2<T,TDESCRIPTOR>() * (converter.getLatticeThermalRelaxationTime() – 0.5));
ADlattice.setParameter<TotalEnthalpy::L>(L);

/// define initial conditions
AnalyticalConst3D<T,T> rho(1.);
AnalyticalConst3D<T,T> u0(0.0, 0.0, 0.0);
AnalyticalConst3D<T,T> T_cold(lattice_Hcold);
AnalyticalConst3D<T,T> T_hot(lattice_Hhot);

/// for each material set Rho, U and the Equilibrium
NSlattice.defineRhoU(superGeometry.getMaterialIndicator({1, 2, 3, 4}), rho, u0);
NSlattice.iniEquilibrium(superGeometry.getMaterialIndicator({1, 2, 3, 4}), rho, u0);

ADlattice.defineField<descriptors::VELOCITY>(superGeometry.getMaterialIndicator({1, 2, 3}), u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({1, 3}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({1, 3}), T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);

/// Make the lattice ready for simulation
NSlattice.initialize();
ADlattice.initialize();

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

void setBoundaryValues( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
int iT, SuperGeometry<T,3>& superGeometry)
{

// nothing to do here

}

void getResults( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice, int iT,
SuperGeometry<T,3>& superGeometry,
util::Timer<T>& timer,
bool converged)
{

OstreamManager clout(std::cout,”getResults”);

SuperVTMwriter3D<T> vtkWriter(“thermalNaturalConvection3D”);
SuperLatticeGeometry3D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeField3D<T, TDESCRIPTOR, VELOCITY> velocity(ADlattice);
SuperLatticePhysPressure3D<T, NSDESCRIPTOR> pressure(NSlattice, converter);

SuperLatticeDensity3D<T, TDESCRIPTOR> enthalpy(ADlattice);
enthalpy.getName() = “enthalpy”;
SuperLatticeField3D<T, NSDESCRIPTOR, POROSITY> liquid_frac(NSlattice);
liquid_frac.getName() = “liquid fraction”;
SuperLatticeField3D<T, TDESCRIPTOR, TEMPERATURE> temperature(ADlattice);
temperature.getName() = “temperature”;
SuperLatticeField3D<T, NSDESCRIPTOR, FORCE> force(NSlattice);
force.getName() = “force”;
vtkWriter.addFunctor( geometry );
vtkWriter.addFunctor( pressure );
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( enthalpy );
vtkWriter.addFunctor( liquid_frac );
vtkWriter.addFunctor( temperature );
vtkWriter.addFunctor( force );

const int vtkIter = converter.getLatticeTime(0.5);

if (iT == 0) {
/// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry3D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeCuboid3D<T, NSDESCRIPTOR> cuboid(NSlattice);
SuperLatticeRank3D<T, NSDESCRIPTOR> rank(NSlattice);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);

vtkWriter.createMasterFile();
}

/// Writes the VTK files
if (iT % vtkIter == 0 || converged) {

timer.update(iT);
timer.printStep();

/// NSLattice statistics console output
NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));

/// ADLattice statistics console output
ADlattice.getStatistics().print(iT,converter.getPhysTime(iT));

if ( NSlattice.getStatistics().getAverageRho() != NSlattice.getStatistics().getAverageRho() or ADlattice.getStatistics().getAverageRho() != ADlattice.getStatistics().getAverageRho() ) {
clout << “simulation diverged! stopping now.” << std::endl;
exit(1);
}
vtkWriter.write(iT);
}
}

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

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

T char_lattice_u = 0.2;

if (argc >= 2) {
N = atoi(argv[1]);
}
if (argc >= 3) {
tau = atof(argv[2]);
}
if (argc >= 4) {
char_lattice_u = atof(argv[3]);
}

const T char_u = util::sqrt( 9.81 * 1.2e-4 * (311. – 302.8) * 6093. );
const T conversion_u = char_u / char_lattice_u;

physDeltaX = lx / N;
physDeltaT = physDeltaX / conversion_u;
physDeltaT = 6093. / 1.81e-3 / descriptors::invCs2<T,NSDESCRIPTOR>() * (tau – 0.5) * physDeltaX * physDeltaX;

lattice_Hcold = cp_s * Tcold;
lattice_Hhot = cp_l * Thot;

clout << “H_cold ” << lattice_Hcold << ” H_hot ” << lattice_Hhot << std::endl;

ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const converter(
(T) physDeltaX, // physDeltaX
(T) physDeltaT, // physDeltaT
(T) lx, // charPhysLength
(T) char_u, // charPhysVelocity
(T) 1.81e-3 / 6093., // physViscosity
(T) 6093., // physDensity
(T) 32., // physThermalConductivity
(T) 381., // physSpecificHeatCapacity
(T) 1.2e-4, // physThermalExpansionCoefficient
(T) 302.8, // charPhysLowTemperature
(T) 311 // charPhysHighTemperature
);
converter.print();
clout << “lattice cp ” << converter.getLatticeSpecificHeatCapacity(cp_l) << std::endl;

/// === 2nd Step: Prepare Geometry ===

std::vector<T> extend(3,T());
extend[0] = lx + 2*converter.getPhysLength(1);
extend[1] = ly + converter.getPhysLength(1);
extend[2] = lz + 2*converter.getPhysLength(1);
std::vector<T> origin(3,T());
IndicatorCuboid3D<T> cuboid(extend, origin);

/// Instantiation of a cuboidGeometry with weights
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize());
cuboidGeometry.setPeriodicity(false, false, true); // x, y, z

/// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);

/// Instantiation of a superGeometry
SuperGeometry<T,3> superGeometry(cuboidGeometry, loadBalancer);

prepareGeometry(superGeometry, converter);

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

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

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

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

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

NSlattice.addLatticeCoupling(superGeometry, 1, coupling, ADlattice);

//prepareLattice and setBoundaryConditions
prepareLattice(converter, NSlattice, ADlattice, superGeometry);

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

for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT)+1; ++iT) {

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

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

/// === 7th Step: Computation and Output of the Results ===
getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, false);
}

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

}

[mathias]

It seems that you did a mistake with the indicators used in the rename functions. You must take care that the domain of interest at the boundary is within the indicator damain. We usually extent the indicator dammain by half a voxel size — see the other examples in OpenLB.

[Jijo]

Dear mathias,

I have modified the code and the geometry seems fine but the results are still inadequate.

Image Link: https://ibb.co/WkVz0jw

Code snippet:

/* Lattice Boltzmann sample, written in C++, using the OpenLB
* library
*
* Copyright (C) 2020 Maximilian Gaedtke, Larissa Dietz
* E-mail contact: info@openlb.net
* The most recent release of OpenLB can be downloaded at
* <http://www.openlb.net/&gt;
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public
* License along with this program; if not, write to the Free
* Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/

/* The solution for the melting problem (solid-liquid phase change)
coupled with natural convection is found using the lattice Boltzmann
method after Rongzong Huang and Huiying Wu (2015)[1]. The equilibrium
distribution function for the temperature is modified in order to deal
with the latent-heat source term. That way, iteration steps or solving
a group of linear equations is avoided, which results in enhanced efficiency.
The phase interface is located by the current total enthalpy, and
its movement is considered by the immersed moving boundary scheme after
Noble and Torczynski (1998)[2]. This method was validated by comparison
with experimental values (e.g. Gau and Viskanta (1986) [3]).

[1] Rongzong Huang, Huiying Wu, Phase interface effects in the total enthalpy-based lattice
Boltzmann model for solid–liquid phase change, Journal of Computational Physics 294 (2015) 346–362.

[2] D. Noble, J. Torczynski, A lattice-Boltzmann method for partially saturated
computational cells, Int. J. Modern Phys. C 9 (8) (1998) 1189–1202.

[3] C. Gau, R. Viskanta, Melting and Solidification of a Pure Metal on a
Vertikal Wall, Journal of Heat Transfer 108(1) (1986): 174–181.
*/

#include “olb3D.h”
#include “olb3D.hh”

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

using T = double;

using NSDESCRIPTOR = D3Q19<POROSITY,VELOCITY_SOLID,FORCE>;
using TDESCRIPTOR = D3Q7<VELOCITY,TEMPERATURE>;

using TotalEnthalpyAdvectionDiffusionDynamics = TotalEnthalpyAdvectionDiffusionTRTdynamics<T,TDESCRIPTOR>;

// Parameters for the simulation setup
const T lx = 88.9e-3; // length of the channel
const T ly = 63.5e-3; // height of the channel
const T lz = 38.1e-3;
int N = 128; // resolution of the model
T tau = 0.51; // relaxation time
const T Ra = 1.656e6; // Rayleigh number
const T Pr = 0.0216; // Prandtl number
const T Ste = 0.039; // Stephan number
const T maxPhysT = 1140.; // simulation time

const T Tcold = 0.5;
const T Tmelt = (302.8 – 301.3)/(311.0 – 301.3) + 0.5;
const T Thot = 1.5;

const T lambda_s = 33.5; // W / m K
const T lambda_l = 32.0; // W / m K
const T R_lambda = lambda_s/lambda_l;

const T cp_s = 1.0; // J / kg K
const T cp_l = 1.0; // J / kg K
const T R_cp = cp_s/cp_l;

//for this case, the harmonic mean (cp_ref) is applicable
const T cp_ref = 2.0 * cp_s * cp_l / (cp_s + cp_l); // J / kg K

const T R_alpha = lambda_s / lambda_l * cp_l / cp_s;

const T L = cp_l * (Thot – Tmelt) / Ste; // J / kg

T lattice_Hcold, lattice_Hhot;
T physDeltaX, physDeltaT;

/// Stores geometry information in form of material numbers
void prepareGeometry(SuperGeometry<T,3>& superGeometry,
ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter)
{

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

superGeometry.rename(0,4);

std::vector<T> extend(3,T());
extend[0] = lx;
extend[1] = ly;
extend[2] = lz;

std::vector<T> origin(3,T());
origin[0] = converter.getPhysLength(1);
origin[1] = converter.getPhysLength(1);
origin[2] = converter.getPhysLength(1);

IndicatorCuboid3D<T> cuboid2(extend, origin);

superGeometry.rename(4,1,cuboid2);

std::vector<T> extendwallleft(3,T(0));
extendwallleft[0] = converter.getPhysLength(1);
extendwallleft[1] = ly + converter.getPhysLength(1);
extendwallleft[2] = lz + converter.getPhysLength(1);

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

IndicatorCuboid3D<T> wallleft(extendwallleft, originwallleft);

std::vector<T> extendwallright(3,T(0));
extendwallright[0] = converter.getPhysLength(1);
extendwallright[1] = ly + converter.getPhysLength(1);
extendwallright[2] = lz + converter.getPhysLength(1);

std::vector<T> originwallright(3,T(0));
originwallright[0] = lx+converter.getPhysLength(1);
originwallright[1] = 0.0;
originwallright[2] = 0.0;

IndicatorCuboid3D<T> wallright(extendwallright, originwallright);

superGeometry.rename(4,2,1, wallleft);
superGeometry.rename(4,3,1, wallright);

/// 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( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
SuperGeometry<T,3>& superGeometry )
{

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

T omega = converter.getLatticeRelaxationFrequency();
T Tomega = converter.getLatticeThermalRelaxationFrequency();

ADlattice.defineDynamics<NoDynamics>(superGeometry, 0);
NSlattice.defineDynamics<NoDynamics>(superGeometry, 0);

NSlattice.defineDynamics<ForcedPSMBGKdynamics>(superGeometry.getMaterialIndicator({1, 2, 3, 4}));

ADlattice.defineDynamics<TotalEnthalpyAdvectionDiffusionDynamics>(superGeometry.getMaterialIndicator({1, 2, 3}));
ADlattice.defineDynamics<BounceBack>(superGeometry, 4);

/// sets boundary
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry.getMaterialIndicator({2, 3}));
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({2, 3}));

NSlattice.setParameter<descriptors::OMEGA>(omega);

ADlattice.setParameter<descriptors::OMEGA>(Tomega);
ADlattice.setParameter<collision::TRT::MAGIC>(0.25);
ADlattice.setParameter<TotalEnthalpy::T_S>(Tmelt);
ADlattice.setParameter<TotalEnthalpy::T_L>(Tmelt);
ADlattice.setParameter<TotalEnthalpy::CP_S>(cp_s);
ADlattice.setParameter<TotalEnthalpy::CP_L>(cp_l);
ADlattice.setParameter<TotalEnthalpy::LAMBDA_S>(cp_ref / descriptors::invCs2<T,TDESCRIPTOR>() * (converter.getLatticeThermalRelaxationTime() – 0.5) * R_lambda);
ADlattice.setParameter<TotalEnthalpy::LAMBDA_L>(cp_ref / descriptors::invCs2<T,TDESCRIPTOR>() * (converter.getLatticeThermalRelaxationTime() – 0.5));
ADlattice.setParameter<TotalEnthalpy::L>(L);

/// define initial conditions
AnalyticalConst3D<T,T> rho(1.);
AnalyticalConst3D<T,T> u0(0.0, 0.0, 0.0);
AnalyticalConst3D<T,T> T_cold(lattice_Hcold);
AnalyticalConst3D<T,T> T_hot(lattice_Hhot);

/// for each material set Rho, U and the Equilibrium
NSlattice.defineRhoU(superGeometry.getMaterialIndicator({0, 1, 2, 3, 4}), rho, u0);
NSlattice.iniEquilibrium(superGeometry.getMaterialIndicator({0, 1, 2, 3, 4}), rho, u0);

ADlattice.defineField<descriptors::VELOCITY>(superGeometry.getMaterialIndicator({1, 2, 3}), u0);
ADlattice.defineRho(superGeometry.getMaterialIndicator({1, 3}), T_cold);
ADlattice.iniEquilibrium(superGeometry.getMaterialIndicator({1, 3}), T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_hot);
ADlattice.iniEquilibrium(superGeometry, 2, T_hot, u0);

/// Make the lattice ready for simulation
NSlattice.initialize();
ADlattice.initialize();

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

void setBoundaryValues( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice,
int iT, SuperGeometry<T,3>& superGeometry)
{

// nothing to do here

}

void getResults( ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
SuperLattice<T, NSDESCRIPTOR>& NSlattice,
SuperLattice<T, TDESCRIPTOR>& ADlattice, int iT,
SuperGeometry<T,3>& superGeometry,
util::Timer<T>& timer,
bool converged)
{

OstreamManager clout(std::cout,”getResults”);

SuperVTMwriter3D<T> vtkWriter(“thermalNaturalConvection3D”);
SuperLatticeGeometry3D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeField3D<T, TDESCRIPTOR, VELOCITY> velocity(ADlattice);
SuperLatticePhysPressure3D<T, NSDESCRIPTOR> pressure(NSlattice, converter);

SuperLatticeDensity3D<T, TDESCRIPTOR> enthalpy(ADlattice);
enthalpy.getName() = “enthalpy”;
SuperLatticeField3D<T, NSDESCRIPTOR, POROSITY> liquid_frac(NSlattice);
liquid_frac.getName() = “liquid fraction”;
SuperLatticeField3D<T, TDESCRIPTOR, TEMPERATURE> temperature(ADlattice);
temperature.getName() = “temperature”;
SuperLatticeField3D<T, NSDESCRIPTOR, FORCE> force(NSlattice);
force.getName() = “force”;
vtkWriter.addFunctor( geometry );
vtkWriter.addFunctor( pressure );
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( enthalpy );
vtkWriter.addFunctor( liquid_frac );
vtkWriter.addFunctor( temperature );
vtkWriter.addFunctor( force );

const int vtkIter = converter.getLatticeTime(0.5);

if (iT == 0) {
/// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry3D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeCuboid3D<T, NSDESCRIPTOR> cuboid(NSlattice);
SuperLatticeRank3D<T, NSDESCRIPTOR> rank(NSlattice);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);

vtkWriter.createMasterFile();
}

/// Writes the VTK files
if (iT % vtkIter == 0 || converged) {

timer.update(iT);
timer.printStep();

/// NSLattice statistics console output
NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));

/// ADLattice statistics console output
ADlattice.getStatistics().print(iT,converter.getPhysTime(iT));

if ( NSlattice.getStatistics().getAverageRho() != NSlattice.getStatistics().getAverageRho() or ADlattice.getStatistics().getAverageRho() != ADlattice.getStatistics().getAverageRho() ) {
clout << “simulation diverged! stopping now.” << std::endl;
exit(1);
}
vtkWriter.write(iT);
}
}

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

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

T char_lattice_u = 0.2;

if (argc >= 2) {
N = atoi(argv[1]);
}
if (argc >= 3) {
tau = atof(argv[2]);
}
if (argc >= 4) {
char_lattice_u = atof(argv[3]);
}

const T char_u = util::sqrt( 9.81 * 1.2e-4 * (311. – 302.8) * 6093. );
const T conversion_u = char_u / char_lattice_u;

physDeltaX = lx / N;
physDeltaT = physDeltaX / conversion_u;
physDeltaT = 6093. / 1.81e-3 / descriptors::invCs2<T,NSDESCRIPTOR>() * (tau – 0.5) * physDeltaX * physDeltaX;

lattice_Hcold = cp_s * Tcold;
lattice_Hhot = cp_l * Thot;

clout << “H_cold ” << lattice_Hcold << ” H_hot ” << lattice_Hhot << std::endl;

ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> const converter(
(T) physDeltaX, // physDeltaX
(T) physDeltaT, // physDeltaT
(T) lx, // charPhysLength
(T) char_u, // charPhysVelocity
(T) 1.81e-3 / 6093., // physViscosity
(T) 6093., // physDensity
(T) 32., // physThermalConductivity
(T) 381., // physSpecificHeatCapacity
(T) 1.2e-4, // physThermalExpansionCoefficient
(T) Tcold, // charPhysLowTemperature
(T) Thot // charPhysHighTemperature
);
converter.print();
clout << “lattice cp ” << converter.getLatticeSpecificHeatCapacity(cp_l) << std::endl;

/// === 2nd Step: Prepare Geometry ===

std::vector<T> extend(3,T());
extend[0] = lx + 2.0*converter.getPhysLength(1);
extend[1] = ly + 2.0*converter.getPhysLength(1);
extend[2] = lz + 2.0*converter.getPhysLength(1);

std::vector<T> origin(3,T());
IndicatorCuboid3D<T> cuboid(extend, origin);

/// Instantiation of a cuboidGeometry with weights
CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize());

/// Instantiation of a loadBalancer
HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);

/// Instantiation of a superGeometry
SuperGeometry<T,3> superGeometry(cuboidGeometry, loadBalancer);

prepareGeometry(superGeometry, converter);

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

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

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

clout << “wanted value ” << descriptors::invCs2<T,TDESCRIPTOR>() << std::endl;

T boussinesqForcePrefactor = Ra / util::pow(T(N),3) * Pr * util::pow(cp_ref / descriptors::invCs2<T,TDESCRIPTOR>() * (converter.getLatticeThermalRelaxationTime() – 0.5), 2);

clout << “boussinesq ” << Ra / util::pow(T(N), 3) * Pr * lambda_l * lambda_l << std::endl;

TotalEnthalpyPhaseChangeCouplingGenerator3D<T,NSDESCRIPTOR,TotalEnthalpyAdvectionDiffusionDynamics>
coupling(0, converter.getLatticeLength(lx), 0, converter.getLatticeLength(ly), 0, converter.getLatticeLength(lz),
boussinesqForcePrefactor, Tcold, 1., dir);

NSlattice.addLatticeCoupling(superGeometry, 1, coupling, ADlattice);

//prepareLattice and setBoundaryConditions
prepareLattice(converter, NSlattice, ADlattice, superGeometry);

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

for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT)+1; ++iT) {

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

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

/// === 7th Step: Computation and Output of the Results ===
getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, false);
}

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

}

[mathias]

First check if the materials are set correctly using the geometry file and paraview, then check the flow bc and the flow (maybe disable the melting). Finally, check the temerature bc (maybe disable the melting). So, you need to check step by step increasing the complexity..

[Author]

Posts