Problem on temperature diffusion on galliumMelting3D

[qiong]

Hi,

I changed the galliumMelting2D modeling to 3D but it seemed that something is wrong. The modeling can work and I can only see slight temperature change from right and left faces towards inside. It seems that the temperature does not diffuse from the hot wall to the cold wall. Can anyone help me figure it out?
Thanks

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

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

using T = FLOATING_POINT_TYPE;

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

using TotalEnthalpyAdvectionDiffusionDynamics = TotalEnthalpyAdvectionDiffusionTRTdynamics<T,TDESCRIPTOR>;

const T lx = 90e-3; // length of the channel
const T ly = 90e-3; // height of the channel
const T lz = 90e-3;
int N = 100; // resolution of the model
T tau = 0.53; // relaxation time

const T Ra = 2e8; // 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;

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 density = 1.; // kg / m^3 //?

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

T lattice_Hcold, lattice_Hhot;
T physDeltaX, physDeltaT;

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.5*converter.getPhysLength(1);
IndicatorCuboid3D<T> cuboid3(extend, origin);

superGeometry.rename(4,1,cuboid3);

std::vector<T> extendwallleft(3,T(0));
extendwallleft[0] = converter.getPhysLength(1);
extendwallleft[1] = ly;
extendwallleft[2] = lz;
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] = lz;
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);

superGeometry.clean();

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<ForcedPSMBGKdynamics>(superGeometry.getMaterialIndicator({1, 2, 3, 4}));
ADlattice.defineDynamics<TotalEnthalpyAdvectionDiffusionDynamics>(superGeometry.getMaterialIndicator({1, 2, 3}));

setBounceBackBoundary(ADlattice, superGeometry, 4);

setRegularizedTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry.getMaterialIndicator({2, 3}));
setInterpolatedVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({2, 3, 4}));

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

ADlattice.setParameter<descriptors::OMEGA>(Tomega);
ADlattice.setParameter<collision::TRT::MAGIC>(T(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);

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

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

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)
{

}

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(“galliumMelting3d”);
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) {

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

if (iT % vtkIter == 0 || converged) {

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

NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));

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[])
{

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;

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

CuboidGeometry3D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), singleton::mpi().getSize());

HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);

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

prepareGeometry(superGeometry, converter);

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

std::vector<T> dir{0.0, 0.0, 1.0};
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(converter, NSlattice, ADlattice, superGeometry);

util::Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() );
timer.start();

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

setBoundaryValues(converter, NSlattice, ADlattice, iT, superGeometry);

NSlattice.executeCoupling();
NSlattice.collideAndStream();
ADlattice.collideAndStream();

getResults(converter, NSlattice, ADlattice, iT, superGeometry, timer, false);
}

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

}

[stephan]

Dear qiong,

thanks for posting. Unfortunately, I cannot see anything directly from looking at your code and we cannot offer debugging of your complete application.
However, you could double-check that the boundary conditions are all set correctly and geometry indicators are set properly (maybe compare to another 3D example that uses an easier model).
Also, please consider to have a look here https://www.sciencedirect.com/science/article/pii/S0017931019361927.

We have a spring school coming up next year, where we could support you with implementing own setups and models (see https://www.openlb.net/spring-school-2025/ for more information).

BR
Stephan

[Author]

Posts