Liquid-Solid Conjugate Heat Transfer

[wf_guoyl]

Hello, everyone! I have some trouble with the liquid-solid conjugate heat transfer. I have defined two descriptors: NSDESCRIPTOR and TDESCRIPTOR. During thermal simulations, I define two different “omega” to simulate transfer of heat in fluids and solids. The solid part is called “obstacle”. Here is part of my code:

/***************************************************/
prepare Lattice:

Vector<SimType,2> center( lx0/2, ly0/2 );
IndicatorCircle2D<SimType> obstacle( center, 25*conversionLength );

auto obstacleIndicator = superGeometry.getMaterialIndicator({materialObstacle});

// define NS Dynamics
NSlattice.defineDynamics( superGeometry, materialObstacle, &olb::instances::getNoDynamics<SimType, NSDESCRIPTOR>() );
offBC.addZeroVelocityBoundary( superGeometry, materialObstacle, obstacle );

// define AD Dynamics
ADlattice.defineDynamics( fluidIndicator, &advectionDiffusionBulkDynamics );
ADlattice.defineDynamics( obstacleIndicator, &advectionDiffusionObstacleDynamics );

// fluid init
There’s no need to change it.

// thermal init
ADlattice.defineRho( obstacleIndicator, T_init );
ADlattice.iniEquilibrium( obstacleIndicator, T_init, zeroVelocity );

/*************************************************/

But the simulation didn’t work. The error message is “density is NAN”. Specifically, four white blocks appear on the top, bottom, left and right endpoints of the circle obstacle. I don’t know what caused it. Does anyone want to help me? If someone can provide some information on fluid-solid conjugate heat transfer, I would appreciate it very much. My email is wfguoyl@163.com .

[Adrian]

The reason for this is likely an issue in the boundary setup, i.e. not all links between undefined no-dynamics nodes and fluid nodes are reconstructed correctly. This is a relatively common problem that can occur when using Bouzidi boundaries. One trick that helps surprisingly often is to change the discretization resolution by e.g. a single cell. Do you get any warnings printed during simulation setup?

If you can post your full code I’d be happy to take a closer look.

[wf_guoyl]

Hello, Adrian:

Thank you very much for your reply!
Because there are a lot of documents involved, I have sent the documents to your email.
My model is modified based on Julius Weinmiller’s work, which can be found on https://www.openlb.net/forum/users/julius-weinmiller/.

In the code below, I adjust the “Bouzidi” boundary condition to “BounceBack” boundary condition. But the solid part “obstacle” still does not conduct heat.

The following is the main code:

/*****************************************************/

/* OpenLB Extension:
* Phase-change multiphase flow
* Written by Julius Weinmiller for his thesis
* https://www.openlb.net/forum/users/julius-weinmiller/
*
* The most recent release of OpenLB can be downloaded at
* <http://www.openlb.net/&gt;
*/

/* D2Law.cpp:
* This simulation will combine both AD and NS lattice for a multiphase problem
* A droplet is spended without body forces
* periodic boundaries are applied in all three axis
* The termpature is increased and evaporation rate recorded
* The results of the simulations show that the simulations follow the
* well known D2 evaporative law.
*/

#include “olb2D.h”
#include “olb2D.hh” // include full template code
#include “phaseChangeExtension.h”
#include “phaseChangeExtension.hh”

using namespace olb;
#define ThermalMRT

typedef double SimType;
#define NSDESCRIPTOR olb::descriptors::D2Q9<olb::descriptors::VELOCITY, \
olb::descriptors::FORCE, \
olb::descriptors::EXTERNAL_FORCE, \
olb::descriptors::PSI_PSEUDO_RHO>

#ifdef ThermalMRT
#define TDESCRIPTOR olb::descriptors::D2Q5<olb::descriptors::tag::MRT, \
olb::descriptors::VELOCITY, \
olb::descriptors::PHI_THERMAL, \
olb::descriptors::PREV_T_V, \
olb::descriptors::PREV_PHI>
#else
#define TDESCRIPTOR olb::descriptors::D2Q5<olb::descriptors::VELOCITY, \
olb::descriptors::PHI_THERMAL, \
olb::descriptors::PREV_T_V, \
olb::descriptors::PREV_PHI>
#endif

// Parameters for the simulation setup

// Resolutions
const int Nx = 100; // spacial resolution of the model
const SimType Nt = 7; // temporal resolution of the model

const SimType physLength = 100e-6; // m
const SimType lx0 = 2*physLength; // length of channel
const SimType ly0 = 2*physLength; // height of channel
const SimType conversionLength = physLength/Nx ;

const SimType sphereRadius = physLength/5.;
const SimType sphereInterfaceWidth = sphereRadius/4;
const SimType sphereOrigin2[2] = {lx0/2., ly0/2.};
const SimType sphereOrigin[2] = {lx0/2., sphereRadius+sphereInterfaceWidth};
const SimType sphereRadius2 = 25*conversionLength;
const SimType sphereInterfaceWidth2 = sphereRadius2/4;

const SimType rhoVapourMultiplier = 1.0;
const SimType rhoLiquidMultiplier = 0.99;

const SimType physViscosity = 1e-5; // m2 / s

const SimType maxVelocity = 20;
const SimType charTime = physLength/maxVelocity;

const SimType endTimeInit = 10*charTime; // Time to allow droplet to come to rest
const SimType endTimeEvap = 2010.0*charTime; // max. simulation time in s, SI unit

const SimType dx = physLength/ Nx;
const SimType dt = charTime/Nx/Nt;

////
// Thermal constants
const SimType TempCrit = 647; // [K] critical temperature of water in Kelvin

const SimType specificGasConstantVapour = 461.52;
const SimType propertyTemperatureRatio = 0.6;
const SimType physWaterTemp = TempCrit * propertyTemperatureRatio – 273.15;
// http://www.thermopedia.com/content/1150/
// Table 8 for thermal conductivity
// Table 4 for specific heat at constant pressure
// Table 1 for density
const SimType physDensityLiquid = 947.13;
const SimType specificHeatCapacityVolumeLiquid = 3.8e3; // J / kg K
const SimType thermalConductivityLiquid = 682e-3; // W / m K
const SimType physDensityVapour = 0.9643;
const SimType specificHeatCapacityVolumeVapour = 2.3e3 – specificGasConstantVapour; // Cv = Cp – Rsp REAL
const SimType thermalConductivityVapour = 31e-3; // W / m K

const SimType liquidDensityRatio = physDensityLiquid/8.72518;

const SimType specificHeatCapacityVolume = 4e3; // J / kg K
const SimType thermalConductivity = 4; // W / m K

const SimType wallRho = 5.0;
const SimType initLiquidHeightFraction = 0.4;

const int iTsmoothUpdate = 10;
const int rhoRampIter = 5000;

const int materialBase = 0;
const int materialFluid = 1;
const int materialLowerWall = 2;
const int materialUpperWall = 3;
const int materialInflow = 4;
const int materialOutflow = 5;
const int materialObstacle = 6;

void prepareGeometry( SuperGeometry2D<SimType>& superGeometry )
{

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

superGeometry.rename( materialBase, materialLowerWall );

Vector<SimType,2> extend( lx0, ly0 );
Vector<SimType,2> origin;

// create fluid and lower wall
superGeometry.rename( materialLowerWall, materialFluid, 0, 1 );

// create upper wall
origin = { 0, ly0-conversionLength };
extend = { lx0, conversionLength };
olb::IndicatorCuboid2D<SimType> upperBC( extend, origin );
superGeometry.rename( materialLowerWall, materialUpperWall, materialFluid, upperBC );

// create inflow
origin = { 0, 0 };
extend = { conversionLength, ly0-conversionLength };
olb::IndicatorCuboid2D<SimType> inflow( extend, origin );
superGeometry.rename( materialFluid, materialInflow, inflow );

// create outflow
origin = { lx0-conversionLength, 0 };
extend = { conversionLength, ly0-conversionLength };
olb::IndicatorCuboid2D<SimType> outflow( extend, origin);
superGeometry.rename( materialFluid, materialOutflow, outflow);

// creat obstacle
origin = { lx0/2-50*conversionLength, ly0/2 };
extend = { 100*conversionLength, 25*conversionLength };
olb:: IndicatorCuboid2D<SimType> obstacle( extend, origin);
superGeometry.rename( materialFluid, materialObstacle, obstacle);

// 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;
}

// Set up the geometry of the simulation
void prepareLattice( olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> const& converter,
olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
olb::SuperLattice2D<SimType, TDESCRIPTOR>& ADlattice,
olb::Dynamics<SimType, NSDESCRIPTOR>& bulkDynamics,
olb::Dynamics<SimType, TDESCRIPTOR>& advectionDiffusionBulkDynamics,
olb::Dynamics<SimType, TDESCRIPTOR>& advectionDiffusionObstacleDynamics,
olb::BounceBack<SimType, NSDESCRIPTOR>& bounceBackUpperWall,
olb::BounceBack<SimType, NSDESCRIPTOR>& bounceBackLowerWall,
olb::sOnLatticeBoundaryCondition2D<SimType,TDESCRIPTOR>& TboundaryCondition,
olb::sOnLatticeBoundaryCondition2D<SimType,NSDESCRIPTOR>& NSboundaryCondition,
SuperGeometry2D<SimType>& superGeometry,
SimType TempInitial,
SimType rhoVapor, SimType rhoLiquid )
{

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

const SimType omega = converter.getLatticeRelaxationFrequency();
const SimType Tomega = converter.getLatticeThermalRelaxationFrequency();

clout << “NS omega = ” << omega << endl;
clout << “T omega = ” << Tomega << endl;

auto bulkIndicator = superGeometry.getMaterialIndicator({materialFluid, materialLowerWall, materialUpperWall, materialInflow, materialOutflow, materialObstacle});
auto fluidIndicator = superGeometry.getMaterialIndicator({materialFluid, materialInflow, materialOutflow});
auto wallIndicator = superGeometry.getMaterialIndicator({materialUpperWall, materialLowerWall});
auto obstacleIndicator = superGeometry.getMaterialIndicator({materialObstacle});

clout << “Define NS Dynamics …” << endl;
NSlattice.defineDynamics( superGeometry, materialBase, &olb::instances::getNoDynamics<SimType, NSDESCRIPTOR>() );
NSlattice.defineDynamics( fluidIndicator, &bulkDynamics );
NSlattice.defineDynamics( superGeometry, materialUpperWall, &bounceBackUpperWall );
NSlattice.defineDynamics( superGeometry, materialLowerWall, &bounceBackLowerWall );
NSlattice.defineDynamics( superGeometry, materialObstacle, &olb::instances::getBounceBack<SimType, NSDESCRIPTOR>() );

clout << “Define AD Dynamics …” << endl;
ADlattice.defineDynamics( superGeometry, materialBase, &olb::instances::getNoDynamics<SimType, TDESCRIPTOR>() );
ADlattice.defineDynamics( fluidIndicator, &advectionDiffusionBulkDynamics );
ADlattice.defineDynamics( wallIndicator, &advectionDiffusionBulkDynamics );
ADlattice.defineDynamics( obstacleIndicator, &advectionDiffusionObstacleDynamics );

TboundaryCondition.addTemperatureBoundary(superGeometry, materialUpperWall, Tomega);
TboundaryCondition.addTemperatureBoundary(superGeometry, materialLowerWall, Tomega);
TboundaryCondition.addTemperatureBoundary(superGeometry, materialInflow, Tomega);
TboundaryCondition.addTemperatureBoundary(superGeometry, materialOutflow, Tomega);

clout << “Define Initial Conditions …” << endl;
// Fluid Initial conditions
olb::AnalyticalConst2D<SimType,SimType> zeroVelocity( 0,0 );

AnalyticalConst2D<SimType,SimType> noiseVap( rhoVapor*0.1 );
AnalyticalConst2D<SimType,SimType> noiseLiq( rhoLiquid*0.1 );
AnalyticalRandom2D<SimType,SimType> random;

SimType estimatedAdditionalAverageVapour = 2*(wallRho – rhoVapor)*conversionLength/ly0;
clout << “estimatedAdditionalAverageVapour: ” << estimatedAdditionalAverageVapour << std::endl;
// this estimate assumes that 2 full fluid nodes equal to the wall density are needed to compensate for the wall density gradient
// This additional density is then spread over the entire y domain

olb::AnalyticalConst2D<SimType,SimType> constRhoVapour( rhoVapor+estimatedAdditionalAverageVapour );
olb::AnalyticalConst2D<SimType,SimType> constRhoLiquid( rhoLiquid * rhoLiquidMultiplier );
olb::AnalyticalConst2D<SimType,SimType> constRhoUpperBC( rhoVapor );
olb::AnalyticalConst2D<SimType,SimType> constRhoWall( rhoVapor );

// Temperature Inital Conditions
olb::AnalyticalConst2D<SimType,SimType> T_init( converter.getLatticeTemperature(TempInitial));

// Random density distribution across inlet
olb::AnalyticalIdentity2D<SimType,SimType> noiseIndicatorVap( random*noiseVap );
olb::AnalyticalIdentity2D<SimType,SimType> noiseIndicatorLiq( random*noiseLiq );
olb::AnalyticalIdentity2D<SimType,SimType> rhoBase( constRhoVapour );

olb::SmoothIndicatorCircle2D<SimType,SimType> sphere( sphereOrigin, sphereRadius, sphereInterfaceWidth);

olb::AnalyticalIdentity2D<SimType,SimType> sphereIndicator( sphere );

olb::AnalyticalIdentity2D<SimType,SimType> rho3( rhoBase + noiseIndicatorVap + sphereIndicator*( constRhoLiquid – constRhoVapour + noiseIndicatorLiq) );

//Initialize all values of distribution functions to their local equilibrium

// Setting initial psi values
olb::AnalyticalConst2D<SimType,SimType> initialPsi( 1.0 );
NSlattice.defineField<olb::descriptors::PSI_PSEUDO_RHO>( fluidIndicator, initialPsi ) ;

// Fluid init
NSlattice.defineRhoU( superGeometry, materialFluid, rho3, zeroVelocity );
NSlattice.iniEquilibrium( superGeometry, materialFluid, rho3, zeroVelocity );

NSlattice.defineRhoU( superGeometry, materialUpperWall, constRhoUpperBC, zeroVelocity );
NSlattice.iniEquilibrium( superGeometry, materialUpperWall, constRhoUpperBC, zeroVelocity );

NSlattice.defineRhoU( superGeometry, materialLowerWall, constRhoWall, zeroVelocity );
NSlattice.iniEquilibrium( superGeometry, materialLowerWall, constRhoWall, zeroVelocity );

NSlattice.defineRhoU( superGeometry, materialInflow, rho3, zeroVelocity );
NSlattice.iniEquilibrium( superGeometry, materialInflow, rho3, zeroVelocity );

NSlattice.defineRhoU( superGeometry, materialOutflow, rho3, zeroVelocity );
NSlattice.iniEquilibrium( superGeometry, materialOutflow, rho3, zeroVelocity );

// Thermal init
ADlattice.defineRho( superGeometry, materialFluid, T_init );
ADlattice.iniEquilibrium( superGeometry, materialFluid, T_init, zeroVelocity );

ADlattice.defineRho( superGeometry, materialUpperWall, T_init );
ADlattice.iniEquilibrium( superGeometry, materialUpperWall, T_init, zeroVelocity );

ADlattice.defineRho( superGeometry, materialLowerWall, T_init );
ADlattice.iniEquilibrium( superGeometry, materialLowerWall, T_init, zeroVelocity );

ADlattice.defineRho( superGeometry, materialInflow, T_init );
ADlattice.iniEquilibrium( superGeometry, materialInflow, T_init, zeroVelocity );

ADlattice.defineRho( superGeometry, materialOutflow, T_init );
ADlattice.iniEquilibrium( superGeometry, materialOutflow, T_init, zeroVelocity );

ADlattice.defineRho( obstacleIndicator, T_init );
ADlattice.iniEquilibrium( obstacleIndicator, T_init, zeroVelocity );

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

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

// Slowly guides wall density to wanted density
void setBoundaryValues( olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> &converter,
olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
olb::SuperLattice2D<SimType, TDESCRIPTOR>& ADlattice,
int iT,
SuperGeometry2D<SimType>& superGeometry,
olb::BounceBackVariableRho<SimType, NSDESCRIPTOR> &bounceBackWall,
SimType wallRhoMax, SimType rhoVapour,
SimType newTempRatio )
{
olb::OstreamManager clout( std::cout,”setBoundaryValues” );

if ( iT%iTsmoothUpdate==0 ) {

SimType frac[1]= {};

// if below iTmaxStart, get percentage ramp up, else 1
if (iT<= rhoRampIter){
// Smooth start curve, polynomial
olb::PolynomialStartScale<SimType,int> startScale( rhoRampIter, SimType( 1 ) );
int iTvec[1]= {iT};
startScale( frac,iTvec );
}
else {
frac[0] = 1;
}
bounceBackWall.setDensity( frac[0]*(wallRhoMax-rhoVapour) + rhoVapour );
}

}

// Output to console and files
void getResults( olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
olb::SuperLattice2D<SimType, TDESCRIPTOR>& ADlattice,
olb::SuperGeometry2D<SimType>& superGeometry,
olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> &converter,
int iT,
Timer<SimType>& timer,
bool forcedSave=false)
{
olb::OstreamManager clout( std::cout,”getResults” );

olb::SuperVTMwriter2D<SimType> vtmWriterNS( “fluidLatticeResults” );
olb::SuperLatticePhysVelocity2D<SimType, NSDESCRIPTOR> velocity( NSlattice, converter );
olb::SuperLatticeDensity2D<SimType, NSDESCRIPTOR> density( NSlattice );
olb::SuperLatticeVelocity2D<SimType, NSDESCRIPTOR> latticeVelocity( NSlattice );
olb::SuperLatticePhysPressure2D<SimType, NSDESCRIPTOR> pressure( NSlattice, converter );

olb::SuperVTMwriter2D<SimType> vtmWriterThermal( “thermalLatticeResults” );
olb::SuperLatticeDensity2D<SimType, TDESCRIPTOR> latticeTemperature( ADlattice );
olb::SuperLatticePhysTemperature2D<SimType, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);

vtmWriterNS.addFunctor( velocity );
vtmWriterNS.addFunctor( latticeVelocity );
vtmWriterNS.addFunctor( pressure );
vtmWriterNS.addFunctor( density );
vtmWriterThermal.addFunctor( latticeTemperature );
vtmWriterThermal.addFunctor( temperature );

const int vtkIter = converter.getLatticeTime( 2500*charTime/Nx/Nt );
const int statIter = converter.getLatticeTime( 5000*charTime/Nx/Nt );
const int ppmIter = converter.getLatticeTime( 2500*charTime/Nx/Nt );

if ( iT==0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
olb::SuperLatticeGeometry2D<SimType, NSDESCRIPTOR> geometry( NSlattice, superGeometry );
olb::SuperLatticeCuboid2D<SimType, NSDESCRIPTOR> cuboid( NSlattice );
olb::SuperLatticeRank2D<SimType, NSDESCRIPTOR> rank( NSlattice );
vtmWriterNS.write( geometry );
vtmWriterNS.write( cuboid );
vtmWriterNS.write( rank );
vtmWriterNS.createMasterFile();
vtmWriterThermal.createMasterFile();
}

// Writes the vtk files
if ( iT%vtkIter==0 || forcedSave) {
vtmWriterNS.write( iT );
vtmWriterThermal.write( iT );
}

// Writes the ppm files
if ( iT%ppmIter==0 || (iT<=50 && (iT%10==0) ) || forcedSave) {
olb::BlockReduction2D2D<SimType> planeReductionRho( density );
olb::BlockReduction2D2D<SimType> planeReductionT( temperature );
olb::heatmap::plotParam<SimType> thermalPlotParam;
thermalPlotParam.name = “physTemperature”;
thermalPlotParam.minValue = converter.getCharPhysLowTemperature()-10;
thermalPlotParam.maxValue = converter.getCharPhysHighTemperature()+10;
// write output as JPEG
olb::heatmap::write(planeReductionRho, iT);
olb::heatmap::write(planeReductionT, iT, thermalPlotParam);
}

// Writes output on the console
if ( iT%statIter==0 || forcedSave ) {
// Timer console output
timer.update( iT );
timer.printStep();

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

void writeGeometry( olb::SuperLattice2D<SimType, NSDESCRIPTOR>& NSlattice,
olb::SuperGeometry2D<SimType>& superGeometry)
{
olb::SuperVTMwriter2D<SimType> vtmWriter( “d2law2D” );
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
olb::SuperLatticeGeometry2D<SimType, NSDESCRIPTOR> geometry( NSlattice, superGeometry );
olb::SuperLatticeCuboid2D<SimType, NSDESCRIPTOR> cuboid( NSlattice );
olb::SuperLatticeRank2D<SimType, NSDESCRIPTOR> rank( NSlattice );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );
//vtmWriter.createMasterFile();
}

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

// === 1st Step: Initialization ===
bool verbosity = false;
olb::olbInit( &argc, &argv, verbosity );
olb::OstreamManager clout( std::cout,”main” );

// == Define simulations constants
const SimType sigma = 0.35;
const SimType eosA = 0.5/49.0;
std::vector<SimType> rho0 ={1.0};
const SimType G = -1;

SimType molarMass = 18.0152;

const SimType temperatureRatio = 0.7;
const SimType newTempRatio = 0.8;

SimType physThermalConductivity;
for( physThermalConductivity = 100.0e-3; physThermalConductivity<=500.01e-3; physThermalConductivity+=200.0e-3 ){

std::ostringstream outputDirStream;
outputDirStream << std::setprecision(3);
outputDirStream << “./results2D/”;
mkdir( outputDirStream.str().c_str(), S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH );
outputDirStream << “lambda_” << physThermalConductivity << “/”;
olb::singleton::directories().setOutputDir( outputDirStream.str() );

// === 2nd Step: Prepare Geometry ===
olb::Vector<SimType,2> extend( lx0, ly0 );
olb::Vector<SimType,2> origin;
olb::IndicatorCuboid2D<SimType> cuboid( extend, origin );

// Instantiation of a cuboidGeometry with weights
#ifdef PARALLEL_MODE_MPI
const int noOfCuboids = olb::singleton::mpi().getSize();
#else
const int noOfCuboids = 1;
#endif
olb::CuboidGeometry2D<SimType> cuboidGeometry( cuboid, conversionLength, noOfCuboids );
cuboidGeometry.setPeriodicity( false, false );

// Instantiation of a loadBalancer
olb::HeuristicLoadBalancer<SimType> loadBalancer( cuboidGeometry );

// Instantiation of a superGeometry
olb::SuperGeometry2D<SimType> superGeometry( cuboidGeometry, loadBalancer, 2 );

prepareGeometry( superGeometry );

// === 3rd Step: Prepare Lattice ===
olb::SuperLattice2D<SimType, TDESCRIPTOR> ADlattice(superGeometry);
olb::SuperLattice2D<SimType, NSDESCRIPTOR> NSlattice(superGeometry);
writeGeometry(NSlattice, superGeometry);

// Pseudopotential multiphase single component dynamics
// A dimensionless R is needed to recover correct latent heat
SimType Rdimless = 8314.0/molarMass /(dx*dx/dt/dt) * TempCrit;

olb::PengRobinsonExt<SimType,SimType> interactionPotential( G, 0.3443, eosA, 2./21., temperatureRatio, Rdimless );
interactionPotential.write();

SimType rho[2];
olb::MaxwellConstruction<SimType, decltype(interactionPotential)> maxwellConstruction( interactionPotential );
maxwellConstruction(rho, temperatureRatio); //temperatureRatio
SimType densityLatticeVapour = rho[1];
SimType densityLatticeLiquid = rho[0];
clout << “Density liquid: ” << densityLatticeLiquid << std::endl;
clout << “Density vapour: ” << densityLatticeVapour << std::endl;

// == Define converter
olb::ThermalUnitConverter<SimType, NSDESCRIPTOR, TDESCRIPTOR> converter(
(SimType) dx, // physDeltaX: spacing between two lattice cells in __m__
(SimType) dt, // physDeltaT = charLatticeVelocity / charPhysVelocity * dx
(SimType) physLength, // charPhysLength: reference length of simulation geometry
(SimType) maxVelocity, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(SimType) physViscosity, // physViscosity: physical kinematic viscosity in __m^2 / s__
(SimType) liquidDensityRatio, // physDensity: physical density in __kg / m^3__
(SimType) thermalConductivity, // physThermalConductivity
(SimType) specificHeatCapacityVolume, // specificHeatCapacity
(SimType) 0.0, // physThermalExpansionCoefficient
(SimType) temperatureRatio*TempCrit, // charPhysLowTemperature
(SimType) newTempRatio*TempCrit // charPhysHighTemperature
);
// Prints the converter log as console output
converter.print();
converter.write();

olb::FluidProperties<SimType> fluidProperties(
densityLatticeVapour, densityLatticeLiquid,
converter.getLatticeSpecificHeatCapacity( specificHeatCapacityVolume ),
converter.getLatticeSpecificHeatCapacity( specificHeatCapacityVolume ),
converter.getLatticeThermalConductivity( physThermalConductivity ),
converter.getLatticeThermalConductivity( physThermalConductivity ),
converter.getCharPhysLowTemperature(), converter.getCharPhysHighTemperature(), TempCrit );

fluidProperties.write();

// Modified Guo forcing dynamics
olb::ForcedModifiedGuoBGKdynamics<SimType, NSDESCRIPTOR> fluidSimulationDynamics (
converter.getLatticeRelaxationFrequency(),
olb::instances::getBulkMomenta<SimType,NSDESCRIPTOR>(), sigma
);

// Modified Advection Diffusion dynamics
olb::Dynamics<SimType, TDESCRIPTOR>* thermalSimulationDynamics;
olb::Dynamics<SimType, TDESCRIPTOR>* obstacleThermalSimulationDynamics;
#ifdef ThermalMRT
thermalSimulationDynamics = new olb::PseudopotentialAdvectionDiffusionMRTdynamics<SimType, TDESCRIPTOR>(
converter.getLatticeThermalRelaxationFrequency(),
olb::instances::getBulkMomenta<SimType,TDESCRIPTOR>()
);

obstacleThermalSimulationDynamics = new olb::PseudopotentialAdvectionDiffusionMRTdynamics<SimType, TDESCRIPTOR>(
converter.getLatticeThermalRelaxationFrequency()-0.2,
olb::instances::getBulkMomenta<SimType,TDESCRIPTOR>()
);

#else
thermalSimulationDynamics = new olb::PseudopotentialAdvectionDiffusionBGKdynamics<SimType, TDESCRIPTOR>(
converter.getLatticeThermalRelaxationFrequency(),
olb::instances::getBulkMomenta<SimType,TDESCRIPTOR>()
);
#endif

// boundary condition
olb::sOnLatticeBoundaryCondition2D<SimType,TDESCRIPTOR> TboundaryCondition(ADlattice);
olb::createAdvectionDiffusionBoundaryCondition2D<SimType,TDESCRIPTOR>(TboundaryCondition);

// A bounce-back node with fictitious density
// which is experienced by the partner fluid
olb::BounceBackVariableRho<SimType, NSDESCRIPTOR> bounceBackUpperWall(densityLatticeVapour);

// A bounce-back node with fictitious density
// which is experienced by the partner fluid
olb::BounceBackVariableRho<SimType, NSDESCRIPTOR> bounceBackLowerWall(densityLatticeVapour);

olb::sOnLatticeBoundaryCondition2D<SimType,NSDESCRIPTOR> NSboundaryCondition(NSlattice);

prepareLattice( converter,
NSlattice, ADlattice,
fluidSimulationDynamics, *thermalSimulationDynamics,
*obstacleThermalSimulationDynamics,
bounceBackUpperWall, bounceBackLowerWall,
TboundaryCondition, NSboundaryCondition,
superGeometry,
TempCrit * temperatureRatio,
densityLatticeVapour, densityLatticeLiquid );

// === 4th Step: Prepare Lattice Coupling
olb::LatticeCouplingGenerator2D<SimType, NSDESCRIPTOR> *thermalMultiphaseCoupling;
#ifdef ThermalMRT
thermalMultiphaseCoupling = new olb::CombinedShanChenThermalMRTCouplingGenerator2D<SimType, NSDESCRIPTOR>(
0, converter.getLatticeLength(lx0),
0, converter.getLatticeLength(ly0),
G, converter.getLatticeThermalRelaxationFrequency(),
0, 0, {0,1},
fluidProperties, interactionPotential);
#else
thermalMultiphaseCoupling = new olb::CombinedShanChenThermalCouplingGenerator2D<SimType, NSDESCRIPTOR>(
0, converter.getLatticeLength(lx0),
0, converter.getLatticeLength(ly0),
G, converter.getLatticeThermalRelaxationFrequency(),
0, 0, {0,1},
fluidProperties, interactionPotential);
#endif

NSlattice.addLatticeCoupling(*thermalMultiphaseCoupling, ADlattice);

// === 5th Step: Main Loop with Timer ===
clout << “starting simulation…” << endl;
olb::util::Timer<SimType> timer( converter.getLatticeTime( endTimeEvap ), superGeometry.getStatistics().getNvoxel() );
timer.start();
int latticeInitTime = converter.getLatticeTime( endTimeInit );

//setBoundaryValues( converter, NSlattice, ADlattice, 0, superGeometry, temperatureRatio);
getResults(NSlattice, ADlattice, superGeometry, converter, 0, timer );

olb::AnalyticalConst2D<SimType,SimType> zeroVelocity( 0,0 );
olb::AnalyticalConst2D<SimType,SimType> T_const( converter.getLatticeTemperature(TempCrit * temperatureRatio));

// ******************the first loop******************
for ( int iT = 0; iT < latticeInitTime; ++iT ){

// === 6th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( converter, NSlattice, ADlattice, iT, superGeometry, bounceBackLowerWall, wallRho, densityLatticeVapour, temperatureRatio);
// === 7th Step: Collide and Stream Execution ===

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

// === Reset temperature field
ADlattice.defineRho(superGeometry, materialFluid, T_const);
ADlattice.iniEquilibrium(superGeometry, materialFluid, T_const, zeroVelocity);

ADlattice.defineRho(superGeometry, materialInflow, T_const);
ADlattice.iniEquilibrium(superGeometry, materialInflow, T_const, zeroVelocity);

ADlattice.defineRho(superGeometry, materialOutflow, T_const);
ADlattice.iniEquilibrium(superGeometry, materialOutflow, T_const, zeroVelocity);

// === 8th Step: Computation and Output of the Results ===
getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer );
if ( std::isnan( NSlattice.getStatistics().getAverageRho() ) ) {
clout << “Terminating this run at step ” << iT+1 << ” => density is NAN” << std::endl;
getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer, true );
break;
}
}

olb::AnalyticalConst2D<SimType,SimType> T_new( converter.getLatticeTemperature(TempCrit * newTempRatio));
olb::AnalyticalConst2D<SimType,SimType> T_newHigher( converter.getLatticeTemperature(TempCrit * 0.9));

ADlattice.defineRho(superGeometry, materialUpperWall, T_new);
ADlattice.defineRho(superGeometry, materialLowerWall, T_new);
ADlattice.defineRho(superGeometry, materialInflow, T_newHigher);
ADlattice.defineRho(superGeometry, materialOutflow, T_new);

// ******************the second loop******************
for ( int iT = latticeInitTime; iT < converter.getLatticeTime( endTimeEvap ); ++iT ) {

// === 6th Step: Definition of Initial and Boundary Conditions ===

// === 7th Step: Collide and Stream Execution ===

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

// === 8th Step: Computation and Output of the Results ===
getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer );
if ( std::isnan( NSlattice.getStatistics().getAverageRho() ) ) {
clout << “Terminating this run at step ” << iT+1 << ” => density is NAN” << std::endl;
getResults(NSlattice, ADlattice, superGeometry, converter, iT+1, timer, true );
break;
}
}
delete thermalSimulationDynamics; // free previous alloc to avoid memory leak
delete thermalMultiphaseCoupling;

timer.stop();
timer.printSummary();
} // end for loop
}

[wf_guoyl]

Now let me recap the question, I define two different “omega” to simulate transfer of heat in fluids and solids(obstacle). That is “obstacleThermalSimulationDynamics” and “thermalSimulationDynamics”. The program can run. But heat transfer can not occur in the solid part (obstacle). If anyone could help me, I would be very grateful.

[Adrian]

I received your mail and was able to compile and run the cuboid case.

Contrary to what you want to achieve (as I understand it) the advectionDiffusionObstacleDynamics are not actually applied to the entire obstacle cuboid on the AD lattice. The reason for this is that a material indicator for number 6 is used but this is only set on the outer boundary of the obstacle due to the default geometry cleaning. You can see this in the geometry VTK output.

I removed the clean calls and the obstacle now takes part in the AD dynamics. Alternatively one could use the original analytical obstacle indicator instead of a material indicator.

• This reply was modified 2 years, 7 months ago by Adrian.

[wf_guoyl]

Hello, Adrian

Thank you very much for your answer! I delete “superGeometry.clean()” and “superGeometry.innerClean()”, and the program can run very well! But I don’t understand where the “original analytical obstacle indicator” is. Do you mean “SmoothIndicatorCuboid2D”? Are there any examples of this in OpenLB?

Thanks again for your kind help!

[Adrian]

Good to hear!

I mean the indicator you use to set the material number of the obstacle:

olb:: IndicatorCuboid2D<SimType> obstacle( extend, origin);
superGeometry.rename( materialFluid, materialObstacle, obstacle);

You could pass this indicator directly to defineDynamics. But it makes sense to maintain a 1:1 relation between material and dynamics here, so no need to change it if it works.

[wf_guoyl]

Thank you for your quick reply!

OpenLB is an excellent piece of software, thanks again to you and all the developers!

Best wishes!

[Author]

Posts