Incomplete indicator shape
OpenLB – Open Source Lattice Boltzmann Code › Forums › on OpenLB › General Topics › Incomplete indicator shape
- This topic has 7 replies, 3 voices, and was last updated 3 years, 1 month ago by Adrian.
-
AuthorPosts
-
December 23, 2021 at 12:49 pm #6219MikeParticipant
Dear everyone,
Recently I have tried to simulate flow in finned(sin-shaped fin) channel,but I have difficulty in creating the sin-shaped fin.The post-process image both in openlb and paraview seem to be weird.Here is my code:Indicator:
template <typename S>
IndicatorSin2D<S>::IndicatorSin2D(Vector<S,2> _center, S _radius, S _theta, S _w, S _A)
: center(_center),
radius(_radius),
theta(_theta),
w(_w),
A(_A)
{
this->_myMin = _center – radius;
this->_myMax = _center + radius;
}template <typename S>
bool IndicatorSin2D<S>::operator()(bool output[], const S input[])
{
output[0] = (input[0] <= this->_myMax[0] && input[0] >= this->_myMin[0]) &&
( input[1] <= this->A * sin(this->w * (input[0] – this->theta))
&& input[1] >= 0);
return output[0];
}.cpp
#include “olb2D.h”
#ifndef OLB_PRECOMPILED // Unless precompiled version is used,
#include “olb2D.hh” // include full template code
#endif
#include <vector>
#include <cmath>
#include <iostream>
#include <fstream>using namespace olb;
using namespace olb::descriptors;
using namespace olb::graphics;
using namespace olb::util;
using namespace std;typedef double T;
typedef D2Q9<FORCE> NSDESCRIPTOR;
typedef D2Q5<VELOCITY> TDESCRIPTOR;//#define M_PI 3.14159265358979323846;
// Parameters for the simulation setup
const T lx = 0.005; // length of the channel
const T ly = 0.001; // height of the channel
const int N = 60; // resolution
const T maxPhysT = 0.5; // max. simulation time in s, SI unit
const T epsilon = 1.e-4; // precision of the convergence (residuum)const T Thot = 289; // temperature of the lower wall in Kelvin
const T Tcold = 288; // temperature of the fluid and top wall in Kelvin
//parameters of sinfin
const T amplitude = 0.0002;
const T w = 5e3 * M_PI;
const T theta = 0.0024;/// Stores geometry information in form of material numbers
void prepareGeometry(SuperGeometry2D<T>& superGeometry,
ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter)
{OstreamManager clout(std::cout,”prepareGeometry”);
clout << “Prepare Geometry …” << std::endl;superGeometry.rename(0,2);
//superGeometry.rename(2,1,0,1);
superGeometry.rename(2,1,1,1 );
Vector<T, 2>center(0.5*lx, 0.5*amplitude);
IndicatorSin2D<T>sinfin(center, 0.5*amplitude, theta, w, amplitude);std::vector<T> extend( 2, T(0) );
extend[0] = lx;
extend[1] = converter.getPhysLength(1);
std::vector<T> origin( 2, T(0) );
IndicatorCuboid2D<T> bottom(extend, origin);origin[1] = ly-converter.getPhysLength(1);
IndicatorCuboid2D<T> top(extend, origin);std::vector<T> extendwallleft(2,T(0));
extendwallleft[0] = converter.getPhysLength(1);
extendwallleft[1] = ly;
std::vector<T> originwallleft(2,T(0));
originwallleft[0] = 0.0;
originwallleft[1] = 0.0;
IndicatorCuboid2D<T> inflow(extendwallleft, originwallleft);std::vector<T> extendwallright(2,T(0));
extendwallright[0] = converter.getPhysLength(1);
extendwallright[1] = ly;
std::vector<T> originwallright(2,T(0));
originwallright[0] = lx-converter.getPhysLength(1);
originwallright[1] = 0.0;
IndicatorCuboid2D<T> outflow(extendwallright, originwallright);/// Set material numbers for bottom, top and pertubation
superGeometry.rename(2,6,1,bottom);
// superGeometry.rename(1,5,test);
superGeometry.rename(2,2,1,top);
superGeometry.rename(2,3,1,inflow );
superGeometry.rename(2,4,1,outflow );
superGeometry.rename(1,5,sinfin);/// 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> &converter,
SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
SuperLattice2D<T, TDESCRIPTOR>& ADlattice,
Dynamics<T, NSDESCRIPTOR> &bulkDynamics,
Dynamics<T, TDESCRIPTOR>& advectionDiffusionBulkDynamics,
Dynamics<T, TDESCRIPTOR>& sourcedAdvectionDiffusionBGKdynamics,
SuperGeometry2D<T>& superGeometry )
{OstreamManager clout(std::cout,”prepareLattice”);
T Tomega = converter.getLatticeThermalRelaxationFrequency();
const T omega = converter.getLatticeRelaxationFrequency();
/// define lattice Dynamics
clout << “defining dynamics” << endl;
ADlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, TDESCRIPTOR>());
NSlattice.defineDynamics(superGeometry, 0, &instances::getNoDynamics<T, NSDESCRIPTOR>());
ADlattice.defineDynamics(superGeometry, 1, &advectionDiffusionBulkDynamics);
ADlattice.defineDynamics(superGeometry, 2, &advectionDiffusionBulkDynamics);
ADlattice.defineDynamics(superGeometry, 3, &advectionDiffusionBulkDynamics);
ADlattice.defineDynamics(superGeometry, 4, &advectionDiffusionBulkDynamics);
ADlattice.defineDynamics(superGeometry, 5, &advectionDiffusionBulkDynamics);
ADlattice.defineDynamics(superGeometry, 6, &advectionDiffusionBulkDynamics);
//ADlattice.defineDynamics(superGeometry, 5, &sourcedAdvectionDiffusionBGKdynamics);NSlattice.defineDynamics(superGeometry, 1, &bulkDynamics);
NSlattice.defineDynamics(superGeometry, 2, &instances::getBounceBack<T, NSDESCRIPTOR>());
NSlattice.defineDynamics(superGeometry, 3, &bulkDynamics);//xinjia
NSlattice.defineDynamics(superGeometry, 4, &bulkDynamics);
NSlattice.defineDynamics(superGeometry, 5, &instances::getBounceBack<T, NSDESCRIPTOR>());
NSlattice.defineDynamics(superGeometry, 6, &instances::getBounceBack<T, NSDESCRIPTOR>());
//NSlattice.defineDynamics(superGeometry, 7, &instances::getNoDynamics<T, NSDESCRIPTOR>());
/// sets boundary
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 2);
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 3);
//setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 4);
//xinjia
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 5);
//setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 5);
setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 6);
//setAdvectionDiffusionTemperatureBoundary<T,TDESCRIPTOR>(ADlattice, Tomega, superGeometry, 7);
setLocalPressureBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 4);
//setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry.getMaterialIndicator({2, 3}));
setLocalVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, omega, superGeometry, 3);// Vector<T, 2>center(0.5*lx, 0.5*amplitude);
// IndicatorSin2D<T>sinfin(center, 0.5*amplitude, theta, w, amplitude);
// setBouzidiZeroVelocityBoundary<T,NSDESCRIPTOR>(NSlattice, superGeometry, 5, sinfin);
// setBouzidiVelocityBoundary<T,TDESCRIPTOR>(ADlattice, superGeometry, 5, sinfin);AnalyticalConst2D<T,T> rho(1.225);
AnalyticalConst2D<T,T> u0(0.0, 0.0);
AnalyticalConst2D<T,T> T_cold(converter.getLatticeTemperature(Tcold));
AnalyticalConst2D<T,T> T_hot(converter.getLatticeTemperature(Thot));T maxVelocity = converter.getCharLatticeVelocity();
T distance2Wall = converter.getConversionFactorLength()/2;
Poiseuille2D<T> u( superGeometry, 3, maxVelocity, distance2Wall );/// for each material set Rho, U and the Equilibrium
NSlattice.defineRhoU(superGeometry, 1, rho, u0);
NSlattice.iniEquilibrium(superGeometry, 1, rho, u0);
NSlattice.defineRhoU(superGeometry, 2, rho, u0);
NSlattice.iniEquilibrium(superGeometry, 2, rho, u0);
NSlattice.defineRhoU(superGeometry, 3, rho, u);
NSlattice.iniEquilibrium(superGeometry, 3, rho, u);
NSlattice.defineRhoU( superGeometry, 4, rho, u0 );
NSlattice.iniEquilibrium( superGeometry, 4, rho, u0 );
NSlattice.defineRhoU(superGeometry, 5, rho, u0);
NSlattice.iniEquilibrium(superGeometry, 5, rho, u0);
NSlattice.defineRhoU(superGeometry, 6, rho, u0);
NSlattice.iniEquilibrium(superGeometry, 6, rho, u0);
// NSlattice.defineRhoU(superGeometry, 7, rho, u0);
// NSlattice.iniEquilibrium(superGeometry, 7, rho, u0);ADlattice.defineRho(superGeometry, 1, T_cold);
ADlattice.iniEquilibrium(superGeometry, 1, T_cold, u0);
ADlattice.defineRho(superGeometry, 2, T_cold);
ADlattice.iniEquilibrium(superGeometry, 2, T_cold, u0);
ADlattice.defineRho(superGeometry, 3, T_cold);
ADlattice.iniEquilibrium(superGeometry, 3, T_cold, u);
ADlattice.defineRho(superGeometry, 4, T_cold);
ADlattice.iniEquilibrium(superGeometry, 4, T_cold, u0);
ADlattice.defineRho(superGeometry, 5, T_hot);
ADlattice.iniEquilibrium(superGeometry, 5, T_hot, u0);
ADlattice.defineRho(superGeometry, 6, T_hot);
ADlattice.iniEquilibrium(superGeometry, 6, T_hot, u0);
// ADlattice.defineRho(superGeometry, 7, T_hot);
// ADlattice.iniEquilibrium(superGeometry, 7, T_hot, u0);NSlattice.initialize();
ADlattice.initialize();
clout << “Prepare Lattice … OK” << std::endl;
}void setBoundaryValues(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
SuperLattice2D<T, TDESCRIPTOR>& ADlattice,
int iT, SuperGeometry2D<T>& superGeometry)
{}
void getResults(ThermalUnitConverter<T, NSDESCRIPTOR, TDESCRIPTOR> &converter,
SuperLattice2D<T, NSDESCRIPTOR>& NSlattice,
SuperLattice2D<T, TDESCRIPTOR>& ADlattice, int iT,
SuperGeometry2D<T>& superGeometry,
Timer<T>& timer,
bool converged)
{OstreamManager clout(std::cout,”getResults”);
SuperVTMwriter2D<T> vtkWriter(“rectangle2″);
SuperLatticePhysVelocity2D<T, NSDESCRIPTOR> velocity(NSlattice, converter);
SuperLatticePhysPressure2D<T, NSDESCRIPTOR> presure(NSlattice, converter);
SuperLatticePhysTemperature2D<T, NSDESCRIPTOR, TDESCRIPTOR> temperature(ADlattice, converter);
SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry1(NSlattice, superGeometry);
vtkWriter.addFunctor( presure );
vtkWriter.addFunctor( velocity );
vtkWriter.addFunctor( temperature );
vtkWriter.addFunctor( geometry1 );const int saveIter = converter.getLatticeTime(0.005);
if (iT == 0) {
SuperLatticeGeometry2D<T, NSDESCRIPTOR> geometry(NSlattice, superGeometry);
SuperLatticeCuboid2D<T, NSDESCRIPTOR> cuboid(NSlattice);
SuperLatticeRank2D<T, NSDESCRIPTOR> rank(NSlattice);
vtkWriter.write(geometry);
vtkWriter.write(cuboid);
vtkWriter.write(rank);
vtkWriter.createMasterFile();
}
/// Writes the VTK files and prints statistics
if (iT%saveIter == 0 || converged) {
/// Timer console output
timer.update(iT);
timer.printStep();/// Lattice statistics console output
NSlattice.getStatistics().print(iT,converter.getPhysTime(iT));
vtkWriter.write(iT);
BlockReduction2D2D<T> planeReduction( temperature, 600, BlockDataSyncMode::ReduceOnly );
// // write output as JPEG
heatmap::plotParam<T> jpeg_Param;
jpeg_Param.maxValue = Thot;
jpeg_Param.minValue = Tcold;
heatmap::write(planeReduction, iT,jpeg_Param);SuperEuklidNorm2D<T, NSDESCRIPTOR> normVel( velocity );
BlockReduction2D2D<T> planeReduction1( normVel, 600, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
heatmap::plotParam<T> jpeg_Param2;
jpeg_Param2.maxValue = converter.getCharPhysVelocity()*1.5;
heatmap::write(planeReduction1, iT, jpeg_Param2);
}}
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) ly/N, // physDeltaX
(T) 4.63e-7,
(T) ly, // charPhysLength
(T) 0.005, // charPhysVelocity
(T) 1.47e-5, // physViscosity
(T) 1.225, // physDensity
(T) 0.025105, // physThermalConductivity
(T) 1006.3, // physSpecificHeatCapacity
(T) 0.00425,//Ra * 1e-5 * 1e-5 / Pr / 9.81 / (Thot – Tcold) / pow(ly, 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 = 1;//7–>1
#endif
CuboidGeometry2D<T> cuboidGeometry(cuboid, converter.getPhysDeltaX(), noOfCuboids);//cuboidGeometry.setPeriodicity(true, false);
HeuristicLoadBalancer<T> loadBalancer(cuboidGeometry);
SuperGeometry2D<T> superGeometry(cuboidGeometry, loadBalancer, 2);
prepareGeometry(superGeometry, converter);
/// === 3rd Step: Prepare Lattice ===
SuperLattice2D<T, TDESCRIPTOR> ADlattice(superGeometry);
//maybe we can specific the geometry domain,eg inlet
SuperLattice2D<T, NSDESCRIPTOR> NSlattice(superGeometry);ForcedBGKdynamics<T, NSDESCRIPTOR> NSbulkDynamics(
converter.getLatticeRelaxationFrequency(),
instances::getBulkMomenta<T,NSDESCRIPTOR>());AdvectionDiffusionBGKdynamics<T, TDESCRIPTOR> TbulkDynamics (
converter.getLatticeThermalRelaxationFrequency(),
instances::getAdvectionDiffusionBulkMomenta<T,TDESCRIPTOR>()
);
SourcedAdvectionDiffusionBGKdynamics<T, TDESCRIPTOR> SinDynamics (
converter.getLatticeThermalRelaxationFrequency(),
instances::getAdvectionDiffusionBulkMomenta<T,TDESCRIPTOR>()
);// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//
// This coupling must be necessarily be put on the Navier-Stokes lattice!!
// !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!//std::vector<T> dir{0.0, 1.0};//1–>5
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, 1, 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);
prepareLattice(converter,
NSlattice, ADlattice,
NSbulkDynamics, TbulkDynamics,SinDynamics,
superGeometry );
//xinjia sboundarycondition/// === 4th Step: Main Loop with Timer ===
Timer<T> timer(converter.getLatticeTime(maxPhysT), superGeometry.getStatistics().getNvoxel() );
timer.start();util::ValueTracer<T> converge(converter.getLatticeTime(0.005),epsilon);//50–>10
for (std::size_t iT = 0; iT < converter.getLatticeTime(maxPhysT); ++iT) {if (converge.hasConverged()) {
clout << “Simulation converged.” << 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();
}postprocess-image:temperature
I have two questions.The first is how to make the color of fin the same as its temperature.The other one is why the sin-fin shape is incomplete.I would appreciate it if anyone could help me solve this!Yours sincerely,
MikeDecember 26, 2021 at 4:07 pm #6223mathiasKeymasterDear Mike,
I recomment using Paraview for more involved psotprocessing. What Do you mean by “incomplete”?
Best
MathiasDecember 29, 2021 at 5:43 am #6226MikeParticipantHello Mathias,
Thanks for your apply! By “incomplete” I mean the fin shape showed in the post-process image is not what I imagine, it seems part of the fin is missing due to resolution setting or something else.
The image is listed blow,I intend to design a sin-shaped fin, but it came out as a triangle like.
temperature image
Best
Mike- This reply was modified 3 years, 1 month ago by Mike.
December 29, 2021 at 10:26 am #6228AdrianKeymasterDoes the shape move towards what you expect if you increase the resolution (Set as deltaX in the unit converter)? Note that the problem could also be the discretization used by the heatmap writer.
How does the raw VTK output look?
December 30, 2021 at 4:12 am #6230MikeParticipantDear Adrian,
Thanks for your reply! Yes, the fin-shape get finer if I increase the resolution, which does not make sense and prolong the simulation time.Also, I have tried to find if I can solve this problem by change some parameters in heatmap namespace,but I can not.What’s more, I notice that the fin-shape in paraview is similar to what it looks like in openlb getresults,which persuade me it may has nothing to do with heatmap, I don’t know if I get this right……
As for VTK output, I know little about it.And it is too long and unreadable that I hesitate to put it here.But if it helps, I would be glad to do this!
Best
Mike- This reply was modified 3 years, 1 month ago by Mike.
December 30, 2021 at 11:03 am #6232AdrianKeymasterI just tried your code and increasing the resolution seems to work correctly. Due to the use of regular grids in LBM, increasing the resolution is the primary way of resolving geometry details. The corresponding increase in (sequential) simulation time is handled by the very good parallelizability. As your example currently only uses ~20k cells you can significantly increase the resolution even on a single compute node.
w.r.t. VTK: There is no need to look at the raw files for this issue, you already checked the output in ParaView (which visualizes the VTK files)
December 30, 2021 at 3:49 pm #6233MikeParticipantDear Adrian,
Thank you again for your reply! Earlier I was thinking something goes wrong within my code about IndicatorSin2D, or there exist a parameter I don’t find deep inside the Indicator like a geometry resolution or something, which can be adjusted to refine the Indicator shape.
Also I’m curious about what lead to the problem that the fin-shape is not the same as I suppose.Is it because that I created in a wrong way or is it about the display problem or something else.What’s more, when I use ug or freecad to create 3D geometry model and import it to openlb, the fin-shape is pretty nice.What is the mechanism behind this.I would appreciate it if you can help!
Best
MikeJanuary 4, 2022 at 2:52 pm #6240AdrianKeymasterAs you implemented a custom indicator for your fin shape (which is the correct way to go for non-primitive shapes that are not easily reproduced using indicator arithmetic) there is likely something wrong in the custom
operator()
. Can you tell us more about the specific “fin” geometry you want to model?The reason why this is not an issue for CAD-imported geometries is that the STL reader indicator is designed to correctly import arbitrary STLs.
-
AuthorPosts
- You must be logged in to reply to this topic.