[atanuchaudhury1994@gmail.com]

But what is the problem with the given file. Why I am not able to edit it in Solidworks? I changed it to solid body still cannot do it.

• This reply was modified 1 year, 9 months ago by atanuchaudhury1994@gmail.com.

[atanuchaudhury1994@gmail.com]

I tried to make a new file but the simulation terminates like this:

[ThreadPool] Sucessfully initialized, numThreads=1
[Directories] Directory ./tmp/ created.
[Directories] Directory ./tmp/imageData/ created.
[Directories] Directory ./tmp/imageData/data/ created.
[Directories] Directory ./tmp/vtkData/ created.
[Directories] Directory ./tmp/vtkData/data/ created.
[Directories] Directory ./tmp/gnuplotData/ created.
[Directories] Directory ./tmp/gnuplotData/data/ created.
[UnitConverter] —————– UnitConverter information —————–
[UnitConverter] — Parameters:
[UnitConverter] Resolution: N= 20
[UnitConverter] Lattice velocity: latticeU= 0.005
[UnitConverter] Lattice relaxation frequency: omega= 1.88679
[UnitConverter] Lattice relaxation time: tau= 0.53
[UnitConverter] Characteristical length(m): charL= 0.1
[UnitConverter] Characteristical speed(m/s): charU= 0.2
[UnitConverter] Phys. kinematic viscosity(m^2/s): charNu= 0.002
[UnitConverter] Phys. density(kg/m^d): charRho= 1
[UnitConverter] Characteristical pressure(N/m^2): charPressure= 0
[UnitConverter] Mach number: machNumber= 0.00866025
[UnitConverter] Reynolds number: reynoldsNumber= 10
[UnitConverter] Knudsen number: knudsenNumber= 0.000866025
[UnitConverter]
[UnitConverter] — Conversion factors:
[UnitConverter] Voxel length(m): physDeltaX= 0.005
[UnitConverter] Time step(s): physDeltaT= 0.000125
[UnitConverter] Velocity factor(m/s): physVelocity= 40
[UnitConverter] Density factor(kg/m^3): physDensity= 1
[UnitConverter] Mass factor(kg): physMass= 1.25e-07
[UnitConverter] Viscosity factor(m^2/s): physViscosity= 0.2
[UnitConverter] Force factor(N): physForce= 0.0400001
[UnitConverter] Pressure factor(N/m^2): physPressure= 1600
[UnitConverter] ————————————————————-
[prepareGeometry] Prepare Geometry …
[SuperGeometry3D] cleaned 0 outer boundary voxel(s)
cylinder3d: ../../../src/functors/analytical/indicator/indicatorF3D.hh:643: olb::IndicatorCuboid3D<T>::IndicatorCuboid3D(olb::Vector<T, 3>, olb::Vector<T, 3>) [with S = float]: Assertion `_xLength>0 && _yLength>0 && _zLength>0′ failed.
Aborted

[atanuchaudhury1994@gmail.com]

And also when I open the .stl file from SolidWorks I donot get the option of Mesh Files because as per all the youtube videos I saw has this option.

[atanuchaudhury1994@gmail.com]

Dear Bukreev,
I have converted to solid body and is still cannot edit it in Solidworks. When I click at the “Imported1” and the Edit Feature, it tells me to open another file. Can you please suggest me on how to do it?

Thank you

[atanuchaudhury1994@gmail.com]

I tried to edit the stl file with Solidworks. But it seems tobe not editable. Can you suggest any other software for this please.

Thank you

[atanuchaudhury1994@gmail.com]

Dear Bukreev
I want to know one thing. Is the .stl file provided by OpenLB is not editable right?

Thank you

[atanuchaudhury1994@gmail.com]

Thank you. I will do that.

[atanuchaudhury1994@gmail.com]

ok then I will do it withb CAD. I think in Solidworks can also be done right? And I just need to replace the .stl file in the folder. Am I right?

[atanuchaudhury1994@gmail.com]

Dear Bukreev,
I tried to change the model creating a new geometry using Indicators like in cylinder2d but it didnot worked out. Do I need to preceed with the making of .stl file?

I changed this part in the prepreGeometry section:

// // Set material number for cylinder
// origin[0] = superGeometry.getStatistics().getMinPhysR( 2 )[0]+converter.getConversionFactorLength();
// extend[0] = ( superGeometry.getStatistics().getMaxPhysR( 2 )[0]-superGeometry.getStatistics().getMinPhysR( 2 )[0] )/2.;
// std::shared_ptr<IndicatorF3D<T>> cylinder = std::make_shared<IndicatorCuboid3D<T>>( extend, origin );
// superGeometry.rename( 2,5, cylinder );
// Set material number for cylinder (new geometry along Y-axis)
origin[1] = superGeometry.getStatistics().getMinPhysR(2)[1] + extend[1] / 4.; // Centered
extend[1] = extend[1] / 2.; // Increase along Y-axis
std::shared_ptr<IndicatorF3D<T>> cylinder = std::make_shared<IndicatorCuboid3D<T>>(extend, origin);
superGeometry.rename(2, 5, cylinder);

[atanuchaudhury1994@gmail.com]

Thank you for your reply. I will try them.

[atanuchaudhury1994@gmail.com]

ok the 2nd problem is solved for now. Please tell me any solution for the 1st problem I tried like this in the code:
void PrepareGeometry //
// existing code
//Added
// Increase the Y-axis length by 50%
extend[1] = 1.5 * extend[1];
// Removes all not needed boundary voxels outside the surface
superGeometry.clean();
superGeometry.checkForErrors();

superGeometry.print();

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

It didn’t worked out. Can you please tell me how to change the dimenssions? Do I have to edit the code or I have to change the stl file?

[atanuchaudhury1994@gmail.com]

Thank you for the suggestion. I am provoding a code snippet below. Is this the way I need to implement it? Is it correct? Please reply,

void prepareGeometry(UnitConverter<T, DESCRIPTOR> const& converter, IndicatorF3D<T>& indicator,
STLreader<T>& stlReader, SuperGeometry<T, 3>& superGeometry)
{
// … (existing code)

// Set material number for additional cylinders
origin[0] = superGeometry.getStatistics().getMinPhysR(2)[0] + 2 * converter.getConversionFactorLength();
extend[0] = (superGeometry.getStatistics().getMaxPhysR(2)[0] – superGeometry.getStatistics().getMinPhysR(2)[0]) / 2.;

// Create additional cylinder indicators
std::shared_ptr<IndicatorF3D<T>> cylinder2 = std::make_shared<IndicatorCuboid3D<T>>(extend, origin);
superGeometry.rename(2, 6, cylinder2); // Material number 6

std::shared_ptr<IndicatorF3D<T>> cylinder3 = std::make_shared<IndicatorCuboid3D<T>>(extend, origin);
superGeometry.rename(2, 7, cylinder3); // Material number 7

// Removes all not needed boundary voxels outside the surface
superGeometry.clean();
superGeometry.checkForErrors();

superGeometry.print();

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

[atanuchaudhury1994@gmail.com]

The output that I am getting are the VTK files, drag file and the info that I have provided before.

Thank you

[atanuchaudhury1994@gmail.com]

This is the code, the changes that I have made is adding Smagorinsky model, change in some parameters, put SetSlip boundary condition in upper and lower wall and made uniform flow velocity at inlet.

#include “olb2D.h”
#include “olb2D.hh”

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

using T = FLOATING_POINT_TYPE;
using DESCRIPTOR = D2Q9<>;

//Changed
/// Selectable bulk models in this case
enum class BulkModel {
RLB,
Smagorinsky,
ShearSmagorinsky,
ConsistentStrainSmagorinsky,
Krause
};

/// Returns BulkModel from CLI argument
BulkModel bulkModelFromString(std::string name) {
if (name == “RLB”) {
return BulkModel::RLB;
}
if (name == “Smagorinsky”) {
return BulkModel::Smagorinsky;
}
if (name == “ShearSmagorinsky”) {
return BulkModel::ShearSmagorinsky;
}
if (name == “ConsistentStrainSmagorinsky”) {
return BulkModel::ConsistentStrainSmagorinsky;
}
if (name == “Krause”) {
return BulkModel::Krause;
}
throw std::runtime_error(name + ” is not a valid BulkModel”);
}

#define BOUZIDI

// Parameters for the simulation setup
const int N = 200; // resolution of the model
const T Re = 200.; // Reynolds number
const T maxPhysT = 20; // max. simulation time in s, SI unit
const T L = 0.1/N; // latticeL
const T lengthX = 2.2;
//const T lengthY = .41+L;
//const T lengthY = 20 * 2 * radiusCylinder + L; // making the domain height to cylinder diameter ratio
const T lengthY = 2.0+L;
const T centerCylinderX = 0.2;
//const T centerCylinderY = 0.2+L/2.;
const T centerCylinderY = lengthY / 2; // centering the cylinder in the domain
const T radiusCylinder = 0.05;

// Stores geometry information in form of material numbers
void prepareGeometry( UnitConverter<T, DESCRIPTOR> const& converter,
SuperGeometry<T,2>& superGeometry,
std::shared_ptr<IndicatorF2D<T>> circle)
{
OstreamManager clout( std::cout,”prepareGeometry” );
clout << “Prepare Geometry …” << std::endl;

Vector<T,2> extend( lengthX,lengthY );
Vector<T,2> origin;

superGeometry.rename( 0,2 );

superGeometry.rename( 2,1,{1,1} );

// Set material number for inflow
extend[0] = 2.*L;
origin[0] = -L;
IndicatorCuboid2D<T> inflow( extend, origin );
superGeometry.rename( 2,3,1,inflow );
// Set material number for outflow
origin[0] = lengthX-L;
IndicatorCuboid2D<T> outflow( extend, origin );
superGeometry.rename( 2,4,1,outflow );
// Set material number for cylinder
superGeometry.rename( 1,5, circle );

// Removes all not needed boundary voxels outside the surface
superGeometry.clean();
superGeometry.checkForErrors();

superGeometry.print();

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

// Set up the geometry of the simulation
void prepareLattice( SuperLattice<T,DESCRIPTOR>& sLattice,
UnitConverter<T, DESCRIPTOR> const& converter,
SuperGeometry<T,2>& superGeometry,
std::shared_ptr<IndicatorF2D<T>> circle, BulkModel bulkModel) //Changed
{
OstreamManager clout( std::cout,”prepareLattice” );
clout << “Prepare Lattice …” << std::endl;

const T omega = converter.getLatticeRelaxationFrequency();

// Material=1 –>bulk dynamics
auto bulkIndicator = superGeometry.getMaterialIndicator({1});
//sLattice.defineDynamics<BGKdynamics>(bulkIndicator);

// Material=0 –>do nothing
sLattice.defineDynamics<NoDynamics>(superGeometry, 0);

// Material=1 –>bulk dynamics
// Material=3 –>bulk dynamics (inflow)
// Material=4 –>bulk dynamics (outflow)
//auto bulkIndicator = superGeometry.getMaterialIndicator({1, 3, 4});
//Changed
switch (bulkModel) {
case BulkModel::RLB:
sLattice.defineDynamics<RLBdynamics>(bulkIndicator);
break;
case BulkModel::ShearSmagorinsky:
sLattice.defineDynamics<ShearSmagorinskyBGKdynamics>(bulkIndicator);
sLattice.setParameter<collision::LES::Smagorinsky>(0.15);
break;
case BulkModel::Krause:
sLattice.defineDynamics<KrauseBGKdynamics>(bulkIndicator);
sLattice.setParameter<collision::LES::Smagorinsky>(0.15);
case BulkModel::ConsistentStrainSmagorinsky:
sLattice.defineDynamics<ConStrainSmagorinskyBGKdynamics>(bulkIndicator);
sLattice.setParameter<collision::LES::Smagorinsky>(0.15);
break;
default:
sLattice.defineDynamics<SmagorinskyBGKdynamics>(bulkIndicator);
sLattice.setParameter<collision::LES::Smagorinsky>(T(0.15));
break;
}

// Material=2 –>bounce back
//setBounceBackBoundary(sLattice, superGeometry, 2);
setSlipBoundary(sLattice, superGeometry, 2);

// Setting of the boundary conditions

//if boundary conditions are chosen to be local
//setLocalVelocityBoundary(sLattice, omega, superGeometry, 3);
//setLocalPressureBoundary(sLattice, omega, superGeometry, 4);

//if boundary conditions are chosen to be interpolated
setInterpolatedVelocityBoundary(sLattice, omega, superGeometry, 3);
setInterpolatedPressureBoundary(sLattice, omega, superGeometry, 4);

// Material=5 –>bouzidi / bounce back
#ifdef BOUZIDI
setBouzidiBoundary(sLattice, superGeometry, 5, *circle);
#else
setBounceBackBoundary(sLattice, superGeometry, 5);
#endif

// Initial conditions
AnalyticalConst2D<T,T> rhoF( 1 );
std::vector<T> velocity( 2,T( 0 ) );
AnalyticalConst2D<T,T> uF( velocity );

// Initialize all values of distribution functions to their local equilibrium
sLattice.defineRhoU( bulkIndicator, rhoF, uF );
sLattice.iniEquilibrium( bulkIndicator, rhoF, uF );

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

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

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

// Generates a slowly increasing inflow for the first iTMaxStart timesteps
void setBoundaryValues( SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T, DESCRIPTOR> const& converter, int iT,
SuperGeometry<T,2>& superGeometry )
{

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

// No of time steps for smooth start-up
int iTmaxStart = converter.getLatticeTime( maxPhysT*0.4 );
//int iTupdate = 5;
int iTupdate = 5;

if ( iT%iTupdate==0 && iT<= iTmaxStart ) {
// Smooth start curve, sinus
// SinusStartScale<T,int> StartScale(iTmaxStart, T(1));

// Smooth start curve, polynomial
//PolynomialStartScale<T,T> StartScale( iTmaxStart, T( 1 ) );

// // Creates and sets the Poiseuille inflow profile using functors
// T iTvec[1] = {T( iT )};
// T frac[1] = {};
// StartScale( frac,iTvec );
// T maxVelocity = converter.getCharLatticeVelocity()*3./2.*frac[0];
// T distance2Wall = L/2.;
// Poiseuille2D<T> poiseuilleU( superGeometry, 3, maxVelocity, distance2Wall );

// sLattice.defineU( superGeometry, 3, poiseuilleU );

// Creates and sets the uniform inflow profile
Vector<T,2> inletVelocity(3.,0.) ;
inletVelocity[0] = converter.getCharLatticeVelocity();
AnalyticalConst<2,T,T> uniformVelocity(inletVelocity);
sLattice.defineU(superGeometry,3,uniformVelocity);

//clout << “step=” << iT << “; maxVel=” << inletVelocity[0] << std::endl;

sLattice.setProcessingContext<Array<momenta::FixedVelocityMomentumGeneric::VELOCITY>>(
ProcessingContext::Simulation);
}

}

// Computes the pressure drop between the voxels before and after the cylinder
void getResults( SuperLattice<T, DESCRIPTOR>& sLattice,
UnitConverter<T, DESCRIPTOR> const& converter, std::size_t iT,
SuperGeometry<T,2>& superGeometry, util::Timer<T>& timer )
{
OstreamManager clout( std::cout,”getResults” );

// Gnuplot constructor (must be static!)
// for real-time plotting: gplot(“name”, true) // experimental!
static Gnuplot<T> gplot( “drag” );

SuperVTMwriter2D<T> vtmWriter( “cylinder2d” );
SuperLatticePhysVelocity2D<T, DESCRIPTOR> velocity( sLattice, converter );
SuperLatticePhysPressure2D<T, DESCRIPTOR> pressure( sLattice, converter );
SuperLatticeGeometry2D<T, DESCRIPTOR> materials( sLattice, superGeometry );
SuperLatticeRefinementMetricKnudsen2D<T, DESCRIPTOR> quality( sLattice, converter);
SuperRoundingF2D<T> roundedQuality ( quality, RoundingMode::NearestInteger);
SuperDiscretizationF2D<T> discretization ( roundedQuality, 0., 2. );

vtmWriter.addFunctor( materials );
vtmWriter.addFunctor( quality );
vtmWriter.addFunctor( velocity );
vtmWriter.addFunctor( pressure );

const int vtkIter = converter.getLatticeTime( .3 );
const int statIter = converter.getLatticeTime( .1 );

T point[2] = {};
point[0] = centerCylinderX + 3*radiusCylinder;
point[1] = centerCylinderY;
AnalyticalFfromSuperF2D<T> intpolateP( pressure, true );
T p;
intpolateP( &p,point );

if ( iT == 0 ) {
// Writes the geometry, cuboid no. and rank no. as vti file for visualization
SuperLatticeGeometry2D<T, DESCRIPTOR> geometry( sLattice, superGeometry );
SuperLatticeCuboid2D<T, DESCRIPTOR> cuboid( sLattice );
SuperLatticeRank2D<T, DESCRIPTOR> rank( sLattice );
vtmWriter.write( geometry );
vtmWriter.write( cuboid );
vtmWriter.write( rank );

vtmWriter.createMasterFile();
}

if ( iT%statIter == 0 ) {
sLattice.setProcessingContext(ProcessingContext::Evaluation);

// Timer console output
timer.update( iT );
timer.printStep();

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

// Drag, lift, pressure drop
AnalyticalFfromSuperF2D<T> intpolatePressure( pressure, true );
SuperLatticePhysDrag2D<T,DESCRIPTOR> drag( sLattice, superGeometry, 5, converter );

T point1[2] = {};
T point2[2] = {};

point1[0] = centerCylinderX – radiusCylinder;
point1[1] = centerCylinderY;

point2[0] = centerCylinderX + radiusCylinder;
point2[1] = centerCylinderY;

T p1, p2;
intpolatePressure( &p1,point1 );
intpolatePressure( &p2,point2 );

clout << “pressure1=” << p1;
clout << “; pressure2=” << p2;

T pressureDrop = p1-p2;
clout << “; pressureDrop=” << pressureDrop;

int input[3] = {};
T _drag[drag.getTargetDim()];
drag( _drag,input );
clout << “; drag=” << _drag[0] << “; lift=” << _drag[1] << std::endl;

// set data for gnuplot: input={xValue, yValue(s), names (optional), position of key (optional)}
gplot.setData( converter.getPhysTime( iT ), {_drag[0], 5.58}, {“drag(openLB)”, “drag(schaeferTurek)”}, “bottom right”, {‘l’,’l’} );

// every (iT%vtkIter) write an png of the plot
if ( iT%( vtkIter ) == 0 ) {
// writes pngs: input={name of the files (optional), x range for the plot (optional)}
gplot.writePNG( iT, maxPhysT );
}
}

// Writes the vtk files
if ( iT%vtkIter == 0 && iT > 0 ) {
vtmWriter.write( iT );

{
SuperEuklidNorm2D<T, DESCRIPTOR> normVel( velocity );
BlockReduction2D2D<T> planeReduction( normVel, 600, BlockDataSyncMode::ReduceOnly );
// write output as JPEG
heatmap::write(planeReduction, iT);
}
{
BlockReduction2D2D<T> planeReduction( discretization, 600, BlockDataSyncMode::ReduceOnly );
heatmap::plotParam<T> jpeg_scale;
jpeg_scale.name = “quality”;
jpeg_scale.colour = “blackbody”;
heatmap::write( planeReduction, iT, jpeg_scale );
}
}

// write pdf at last time step
if ( iT == converter.getLatticeTime( maxPhysT )-1 ) {
// writes pdf
gplot.writePDF();
}
}

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

UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR> const converter(
int {N}, // resolution: number of voxels per charPhysL
(T) 0.560, // latticeRelaxationTime: relaxation time, have to be greater than 0.5!
(T) 2.0*radiusCylinder, // charPhysLength: reference length of simulation geometry
//(T) 0.1, // charPhysLength: reference length of simulation geometry
(T) 0.2, // charPhysVelocity: maximal/highest expected velocity during simulation in __m / s__
(T) 0.2*2.*radiusCylinder/Re, // physViscosity: physical kinematic viscosity in __m^2 / s__
//(T) 0.00000179, // physViscosity: physical kinematic viscosity in __m^2 / s__
//(T) 1.0 // physDensity: physical density in __kg / m^3__
(T) 997.0 // physDensity: physical density in __kg / m^3__
);
// Prints the converter log as console output
converter.print();
// Writes the converter log in a file
converter.write(“cylinder2d”);

// === 2rd Step: Prepare Geometry ===
Vector<T,2> extend( lengthX,lengthY );
Vector<T,2> origin;
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, L, noOfCuboids );

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

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

Vector<T,2> center( centerCylinderX,centerCylinderY );
std::shared_ptr<IndicatorF2D<T>> circle = std::make_shared<IndicatorCircle2D<T>>( center, radiusCylinder );

prepareGeometry( converter, superGeometry, circle );

// === 3rd Step: Prepare Lattice ===
SuperLattice<T, DESCRIPTOR> sLattice( superGeometry );

//Added
// Set the dynamics model to Smagorinsky
BulkModel bulkModel = bulkModelFromString(“Smagorinsky”);

//prepareLattice and set boundaryConditions
//Changed
prepareLattice( sLattice, converter, superGeometry, circle, bulkModel);

// === 4th Step: Main Loop with Timer ===
clout << “starting simulation…” << std::endl;
util::Timer<T> timer( converter.getLatticeTime( maxPhysT ), superGeometry.getStatistics().getNvoxel() );
timer.start();

for ( std::size_t iT = 0; iT < converter.getLatticeTime( maxPhysT ); ++iT ) {
// === 5th Step: Definition of Initial and Boundary Conditions ===
setBoundaryValues( sLattice, converter, iT, superGeometry );

// === 6th Step: Collide and Stream Execution ===
sLattice.collideAndStream();

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

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

[atanuchaudhury1994@gmail.com]

I am getting these output:
—————– UnitConverter information —————–
— Parameters:
Resolution: N= 200
Lattice velocity: latticeU= 0.02
Lattice relaxation frequency: omega= 1.78571
Lattice relaxation time: tau= 0.56
Characteristical length(m): charL= 0.1
Characteristical speed(m/s): charU= 0.2
Phys. kinematic viscosity(m^2/s): charNu= 0.0001
Phys. density(kg/m^d): charRho= 997
Characteristical pressure(N/m^2): charPressure= 0
Mach number: machNumber= 0.034641
Reynolds number: reynoldsNumber= 200
Knudsen number: knudsenNumber= 0.000173205

— Conversion factors:
Voxel length(m): physDeltaX= 0.0005
Time step(s): physDeltaT= 5e-05
Velocity factor(m/s): physVelocity= 10
Density factor(kg/m^3): physDensity= 997
Mass factor(kg): physMass= 1.24625e-07
Viscosity factor(m^2/s): physViscosity= 0.005
Force factor(N): physForce= 0.024925
Pressure factor(N/m^2): physPressure= 99700
————————————————————-
Can you please specify where I can find them?

[Author]

Posts