Hi,
I am trying to set a flow in a 3D circular pipe based one the aneurysm bounceback case (especially for the boundary conditions). I have a poiseuille flow at the inlet and bounce back boundary conditions for the walls.
I set-up a lattice with a resolution of 160 (160 nodes along a diameter of the inlet) which leads to a pretty big 21.5 M M nodes lattice. With this lattice, only computations with Reynolds lower that 160 converge. Higher Reynolds procude velocity spots close to the wall which lead to velocity divergence.
The mesh requirements seem extremely important for such a low Reynolds number. Does anyone has any experience for the lattices requirements for this kind of wall driven flow ?
Another explanation may be that there is something wrong with my code…
My code is devided into 3 files: One with the main program, one with additional functions and one with the instruction to read the .XML input file. I put them at the end of the post.
Thanks in advance
Julien
Here is the main program file:
[code=“cpp”]
#include “palabos3D.h”
#include “palabos3D.hh”
#include “param_class.h”
#include “subroutines.h”
#include <sys/time.h>
using namespace plb;
using namespace std;
typedef double T;
typedef Array<T,3> Velocity;
#define PLB_DEBUG
#define DESCRIPTOR descriptors::D3Q19Descriptor
plint extraLayer = 1; // Make the bounding box larger; for visualization purposes
// only. For the simulation, it is OK to have extraLayer=0.
const plint blockSize = 50; // Zero means: no sparse representation.
const plint envelopeWidth = 1; // For standard BGK dynamics.
const plint extendedEnvelopeWidth = 0; // Set to 2 for Guo off lattice boundary condition
plint refDirection = 0;
plint openingSortDirection = 2;
plint margin = 3; // Extra margin of allocated cells around the obstacle.
plint borderWidth = 1; // Because the Guo boundary condition acts in a one-cell layer.
// Requirement: margin>=borderWidth.
// plint referenceDirection = 0;
bool poiseuilleInlet = true;
bool performOutput = true;
TriangleSet* triangleSet = 0;
template
struct Opening {
bool inlet;
Array<T,3> center;
T innerRadius;
};
std::vector<Opening > openings;
int main(int argc, char* argv[]) {
plbInit(&argc, &argv);
global::directories().setOutputDir("./");
global::IOpolicy().activateParallelIO(true);
// Time variables declaration
struct timeval start, end;
gettimeofday(&start, NULL);
T delta_T;
/// -------------------------------------------- ///
/// ---------------- Parameters ---------------- ///
/// -------------------------------------------- ///
// Reading of the XML input file
std::string xmlFileName;
try {
global::argv(1).read(xmlFileName);
}
catch (PlbIOException& exception) {
pcout << “Wrong parameters; the syntax is: "
<< (std::string)global::argv(0) << " input-file.xml” << std::endl;
return -1;
}
// Read the XML file and print of the log file
Param_Pipe param;
param.read_XML(xmlFileName);
param.writeLog(param);
string GeoName = “Geometry/Pipe.stl”;
/// -------------------------------------------- ///
/// --------------- Mesh reading --------------- ///
/// -------------------------------------------- ///
// At this part, the surface geometry is read into a data structure
// comprised by a set of triangles. The DBL constant means that double
// precision accuracy will be used (generally the recommended choice).
triangleSet = new TriangleSet(param.get_GeoName(), DBL);
DEFscaledMesh* defMesh =
new DEFscaledMesh(*triangleSet, param.get_resolution()*2, refDirection, margin, extraLayer);
TriangleBoundary3D boundary(*defMesh);
delete defMesh;
boundary.getMesh().inflate();
// Verification of the imposed dx value
T dx_relative = fabs(boundary.getDx() - param.get_dx()) / boundary.getDx();
if (dx_relative > 1.e-6) {
pcout << " Problem with dx calculation: "
<< " Dx from mesh = " << boundary.getDx() << " and "
<< " Dx from param = " << param.get_dx() << endl;
return -1;
}
Array<T,3> location(boundary.getPhysicalLocation());
// Localization of the STL holes to chose inlets or outlets
std::vectorstd::string openingType;
openingType = param.get_OpeningType();
openings.resize(openingType.size());
for (pluint i=0; i<openingType.size(); ++i) {
std::string next_opening = util::tolower(openingType[i]);
if (next_opening==“inlet”) {
openings[i].inlet = true;
}
else if (next_opening==“outlet”) {
openings[i].inlet = false;
}
else {
plbIOError(“Unknown opening type.”);
}
}
T dx = param.get_dx();
for (pluint i=0; i<openings.size(); ++i) {
Opening<T>& opening = openings[i];
opening.center = computeBaryCenter (
boundary.getMesh(),
boundary.getInletOutlet(openingSortDirection)[i] );
opening.innerRadius = computeInnerRadius (
boundary.getMesh(),
boundary.getInletOutlet(openingSortDirection)[i] );
// T innerRadius = computeInnerRadius (boundary.getMesh(),
// boundary.getInletOutlet(openingSortDirection)[i]);
Array<T,3> centre = computeBaryCenter (boundary.getMesh(),
boundary.getInletOutlet(openingSortDirection)[i]);
}
// Voxalization of the geometry
const int flowType = voxelFlag::inside;
VoxelizedDomain3D voxelizedDomain (
boundary, flowType, extraLayer, borderWidth, extendedEnvelopeWidth, blockSize );
if (performOutput) {
pcout << getMultiBlockInfo(voxelizedDomain.getVoxelMatrix()) << std::endl;
}
MultiScalarField3D<int> flagMatrix((MultiBlock3D&)voxelizedDomain.getVoxelMatrix());
setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
voxelFlag::inside, flagMatrix.getBoundingBox(), 1);
setToConstant(flagMatrix, voxelizedDomain.getVoxelMatrix(),
voxelFlag::innerBorder, flagMatrix.getBoundingBox(), 1);
pcout << "Number of fluid cells: " << computeSum(flagMatrix) << std::endl;
/// ---------------------------------------------------- ///
/// --------------- Lattice and Dynamics --------------- ///
/// --------------------------------------------------- ///
// Instantioation of the BGK dynamics
Dynamics<T,DESCRIPTOR>* dynamics = 0;
dynamics = new SmagorinskyBGKdynamics<T,DESCRIPTOR>(param.get_Omega(),param.get_cSmago());
// dynamics = new BGKdynamics<T,DESCRIPTOR>(param.get_Omega());
// Creation of the lattice based on the Voxelized domain
std::auto_ptr<MultiBlockLattice3D<T,DESCRIPTOR> > lattice
= generateMultiBlockLattice<T,DESCRIPTOR> (
voxelizedDomain.getVoxelMatrix(), envelopeWidth, dynamics );
lattice->toggleInternalStatistics(false);
// Create the sponge zones
createSpongeZones(*lattice, param);
/// ---------------------------------------------------- ///
/// ---------------- Boundary condition ---------------- ///
/// --------------------------------------------------- ///
OnLatticeBoundaryCondition3D<T,DESCRIPTOR>* boundaryCondition =
createLocalBoundaryCondition3D<T,DESCRIPTOR>();
T dt = param.get_dt();
plint nx = lattice->getNx();
plint ny = lattice->getNy();
plint nz = lattice->getNz();
T diameterReal = 0.080;
plint diameter = util::roundToInt(diameterReal/dx);
Array<T,3> inletCenter (0., 0., 0.);
inletCenter = openings[0].center;
Box3D inletDomain(util::roundToInt(inletCenter[0])+1-diameter,
util::roundToInt(inletCenter[0])+1+diameter,
util::roundToInt(inletCenter[1])+1-diameter,
util::roundToInt(inletCenter[1])+1+diameter,
util::roundToInt(inletCenter[2])+1,
util::roundToInt(inletCenter[2])+1);
diameterReal = 0.080;
diameter = util::roundToInt(diameterReal/dx);
Array<T,3> outletCenter (0., 0., 0.);
outletCenter = openings[1].center;
Box3D outletDomain(util::roundToInt(outletCenter[0])+1-diameter,
util::roundToInt(outletCenter[0])+1+diameter,
util::roundToInt(outletCenter[1])+1-diameter,
util::roundToInt(outletCenter[1])+1+diameter,
util::roundToInt(outletCenter[2])-1,
util::roundToInt(outletCenter[2])-1);
Box3D behindInlet(inletDomain.x0, inletDomain.x1,
inletDomain.y0, inletDomain.y1,
inletDomain.z0-diameter, inletDomain.z0-1);
Box3D behindOutlet(outletDomain.x0, outletDomain.x1,
outletDomain.y0, outletDomain.y1,
outletDomain.z0+1, outletDomain.z0+diameter);
boundaryCondition->addVelocityBoundary2N(inletDomain, *lattice);
boundaryCondition->addPressureBoundary2P(outletDomain, *lattice);
// setBoundaryVelocity(*lattice, inletDomain, PoiseuilleVelocity(param) );
setBoundaryDensity(*lattice, outletDomain, (T)1.);
defineDynamics(*lattice, flagMatrix, lattice->getBoundingBox(), new BounceBack<T,DESCRIPTOR>(1.),0);
defineDynamics(*lattice, behindInlet, new BounceBack<T,DESCRIPTOR>(1.));
defineDynamics(*lattice, behindOutlet, new BounceBack<T,DESCRIPTOR>(1.));
iniLattice(*lattice, voxelizedDomain, param, inletCenter);
/// ----------------------------------------------------- ///
/// ---------------- Collide and Stream ----------------- ///
/// ----------------------------------------------------- ///
// Time steps loop
plb_ofstream PressureFile(“PROBE/Head_losses.dat”);
plb_ofstream EnergyFile(“PROBE/Energy.dat”);
plb_ofstream MaxVelFile(“PROBE/MaxVelocity.dat”);
plb_ofstream VelInletFile(“PROBE/InletVelocity.dat”);
plb_ofstream VelOutletFile(“PROBE/OutletVelocity.dat”);
plb_ofstream TimeFile(“PROBE/Time.dat”);
for (plint iT=0; iTdt<param.get_MaxT(); ++iT) {
if (iTdt < param.get_RampeTime()) {
T velFactor = iT*dt / param.get_RampeTime();
setBoundaryVelocity(*lattice, inletDomain, PoiseuilleVelocity(param,velFactor,
inletCenter) );
// setBoundaryVelocity(lattice, inletDomain, Array<T,3> (0., 0., velFactorparam.get_u_LB() ));
}
if (iT%25 == 0) {
// Calcul de la vitesse max
std::auto_ptr<MultiScalarField3D<T> > velocityNorm = computeVelocityNorm (*lattice);
T maxVel = computeMax(*velocityNorm)*dx/dt;
// Calcul de la vitesse moyenne en entrée et en sortie
T AvVel_inlet = computeAverage(*velocityNorm, inletDomain)*dx/dt;
T AvVel_outlet = computeAverage(*velocityNorm, outletDomain)*dx/dt;
// Calcul de l'energie moyenne sur le domaine
T AvEnergy = computeAverageEnergy(*lattice, lattice->getBoundingBox());
// Calcul du DP entre l'entrée et la sortie
T Pinlet = computeAverageDensity<T>(*lattice, inletDomain);
T Poutlet = computeAverageDensity<T>(*lattice, outletDomain);
T pressureDrop = (Pinlet - Poutlet) * DESCRIPTOR<T>::invCs2 * dx*dx / (dt*dt);
// Temps de calcul
gettimeofday(&end, NULL);
delta_T = ((end.tv_sec - start.tv_sec) * 1.e6 + end.tv_usec - start.tv_usec) / 1.e6;
// Ecriture des résultats dans le fichier output
pcout << " T= " << iT*dt << "; "
<< " Av energy: " << AvEnergy
<< " Max Velocity: " << maxVel
<< " Pressure inlet : " << Pinlet
<< endl;
// Ecriture des résultats dans des fichiers
PressureFile << iT*dt << " " << pressureDrop << endl;
EnergyFile << iT*dt << " " << AvEnergy << endl;
MaxVelFile << iT*dt << " " << maxVel << endl;
VelInletFile << iT*dt << " " << AvVel_inlet << endl;
VelOutletFile << iT*dt << " " << AvVel_outlet << endl;
TimeFile << iT*dt << " " << delta_T << endl;
// Vérification de la présence de NaN
if (AvEnergy != AvEnergy) {
pcout << "==> Computation diverged !!!" << endl;
pcout << "==> NaN internal energy" << endl;
return -1;
}
}
if (iT%plint(param.get_gifTime()/dt) == 0) {
Box3D slice (plint(nx/2), plint(nx/2), 0, ny-1, 0, nz-1);
T valMax = param.get_u_Physical();
writeVelImages(*lattice, inletDomain, "tmp/inlet", iT, 0., valMax, dx, dt);
writeVelImages(*lattice, outletDomain, "tmp/outlet", iT, 0., valMax, dx, dt);
writeVelImages(*lattice, slice, "tmp/slice", iT, 0., valMax, dx, dt);
}
if ( (iT%plint(param.get_vtiTime()/dt) == 0) && (iT != 0) ){
pcout << " \n Writing VTK postprocessing files \n" << endl;
writeVTK( *lattice, flagMatrix, lattice->getBoundingBox(),
"tmp/full_domain_"+util::val2str(iT), location, dx, dt);
}
lattice->collideAndStream();
}
PressureFile.close();
EnergyFile.close();
MaxVelFile.close();
VelInletFile.close();
VelOutletFile.close();
TimeFile.close();
}
Here is the file with the subroutines used by the main program:
[code="cpp"]
#include "palabos3D.h"
#include "palabos3D.hh"
#include <string>
#include <fstream>
#include <cmath>
using namespace plb;
using namespace std;
typedef double T;
#define DESCRIPTOR descriptors::D3Q19Descriptor
/// ------------------------------------------------- ///
/// - - - Calculation of the Poiseuille profile - - - ///
/// ------------------------------------------------- ///
/// Velocity on the parabolic Poiseuille profile
T poiseuilleVelocity(plint iX, plint iY, Param_Pipe<T> const& parameters,
Array<T,3> center) {
T dx = parameters.get_dx();
T radius = sqrt(pow((util::roundToInt(center[0])+1-iX),2.)
+ pow((util::roundToInt(center[1])+1-iY),2.));
T radiusRef = ((parameters.get_InletDiam())-dx)/dx;
if (radius < radiusRef) {
T factor = ( pow(radiusRef,2.)-pow(radius,2.) ) / pow(radiusRef,2.);
return parameters.get_u_LB()*factor;
}
else return 0.;
}
/// Class for the Poiseuille velocity profile
template<typename T>
class PoiseuilleVelocity {
public:
PoiseuilleVelocity(Param_Pipe<T> parameters_, T factor_, Array<T,3> center_)
: parameters(parameters_), factor(factor_), center(center_)
{ }
void operator()(plint iX, plint iY, plint iZ, Array<T,3>& u) const {
u[0] = T();
u[1] = T();
u[2] = factor*poiseuilleVelocity(iX,iY, parameters, center);
}
private:
Param_Pipe<T> parameters;
T factor;
Array<T,3> center;
};
/// Class for the Poiseuille velocity profile
template<typename T>
class PoiseuilleVelocityAndDensity {
public:
PoiseuilleVelocityAndDensity(Param_Pipe<T> parameters_, Array<T,3> center_)
: parameters(parameters_), center(center_)
{ }
void operator()(plint iX, plint iY, plint iZ, T& rho, Array<T,3>& u) const {
rho = 1.;
u[0] = T();
u[1] = T();
u[2] = poiseuilleVelocity(iX,iY, parameters, center);
}
private:
Param_Pipe<T> parameters;
Array<T,3> center;
};
/// ------------------------------------------ ///
/// - - - - - Lattice initialization - - - - - ///
/// ------------------------------------------ ///
void iniLattice( MultiBlockLattice3D<T,DESCRIPTOR>& lattice,
VoxelizedDomain3D<T>& voxelizedDomain,
Param_Pipe<T> param, Array<T,3> inletCenter )
{
// Switch all remaining outer cells to no-dynamics, except the outer
// boundary layer, and keep the rest as BGKdynamics.
defineDynamics(lattice, voxelizedDomain.getVoxelMatrix(), lattice.getBoundingBox(),
new NoDynamics<T,DESCRIPTOR>, voxelFlag::outside);
// At the beginning the fluid is considered at rest
initializeAtEquilibrium(lattice, lattice.getBoundingBox(), (T) 1., Array<T,3>((T) 0.,(T) 0.,(T) 0.));
// initializeAtEquilibrium(lattice, lattice.getBoundingBox(),
// PoiseuilleVelocityAndDensity<T>(param, inletCenter) );
lattice.initialize();
}
/// --------------------------------------------- ///
/// - - - - - Post-processing functions - - - - - ///
/// --------------------------------------------- ///
/// subroutines to write .VTK files for paraview
void writeVTK (
MultiBlockLattice3D<T,DESCRIPTOR> lattice,
MultiScalarField3D<int>& flagMatrix,
Box3D const& vtkDomain, std::string fname, Array<T,3> location, T dx, T dt)
{
VtkImageOutput3D<T> vtkOut(fname, dx, location);
vtkOut.writeData<3,float>(*computeVelocity(lattice), "u", dx/dt);
vtkOut.writeData<float>(*computeDensity(lattice), "p", util::sqr(dx/dt)*998.2);
vtkOut.writeData<float>(*copyConvert<int,T>(*extractSubDomain(flagMatrix,vtkDomain)), "flag", 1.);
}
/// subroutines to write .gif images
void writeVelImages ( MultiBlockLattice3D<T,DESCRIPTOR> lattice,
Box3D const& domainImage, std::string keyword, plint it, T valMin, T valMax,
T dx, T dt ) {
// On adimensionalise ces valeurs
valMin = valMin * dt/dx;
valMax = valMax * dt/dx;
string fname = keyword+"_"+util::val2str(it);
ImageWriter<T> imageWriter("leeloo");
imageWriter.writeGif(fname, *computeVelocityNorm(lattice, domainImage), valMin, valMax);
}
/// --------------------------------------------- ///
/// - - - - - - Sponge zones creation - - - - - - ///
/// --------------------------------------------- ///
void createSpongeZones(MultiBlockLattice3D<T,DESCRIPTOR> lattice, Param_Pipe<T> param) {
if (param.get_spongeLayer() < 0) {
pcout << "Generating an outlet viscosity sponge zone." << std::endl;
T bulkValue = param.get_Omega();
Array<plint,6> numSpongeCells;
numSpongeCells[0] = 0;
numSpongeCells[1] = 0;
numSpongeCells[2] = 0;
numSpongeCells[3] = 0;
numSpongeCells[4] = 0;
numSpongeCells[5] = param.get_spongeLayer();
std::vector<MultiBlock3D*> args;
args.push_back(&lattice);
plint nx = lattice.getNx();
plint ny = lattice.getNy();
plint nz = lattice.getNz();
applyProcessingFunctional(new ViscositySpongeZone<T,DESCRIPTOR>(
nx, ny, nz, bulkValue,numSpongeCells),
lattice.getBoundingBox(), args);
}
}
Here is the code for the class Param_Pipe that reads the XML input file and computes the LB parameters.
[code=“cpp”]
#include “palabos3D.h”
#include “palabos3D.hh”
#include
#include
#include
using namespace plb;
using namespace std;
typedef double T;
#define DESCRIPTOR descriptors::D3Q19Descriptor
template
class Param_Pipe {
public:
/// Function to read the XML file
void read_XML(std::string XMLfile) {
XMLreader document(XMLfile);
document["geometry"]["mesh"].read(geoFile);
document["geometry"]["inletDiam"].read(inletDiam);
document["geometry"]["u_Physical"].read(u_Physical);
document["geometry"]["u_LB"].read(u_LB);
document["geometry"]["openings"]["sortDirection"].read(openingSortDirection);
document["geometry"]["openings"]["type"].read(openingType);
document["fluid"]["kinematicViscosity"].read(nu_Physical);
document["fluid"]["density"].read(rho_Physical);
document["numerics"]["cellSize"].read(cellSize);
document["numerics"]["cSmago"].read(cSmago);
document["numerics"]["spongeLayer"].read(spongeLayer);
document["numerics"]["cSmagoTarget"].read(cSmagoTarget);
document["simulation"]["simTime"].read(simTime);
document["simulation"]["rampeTime"].read(rampeTime);
document["simulation"]["gifTime"].read(gifTime);
document["simulation"]["vtiTime"].read(vtiTime);
document["simulation"]["useAllDirections"].read(useAllDirections);
document["simulation"]["useRegularizedWall"].read(useRegularizedWall);
document["simulation"]["poiseuilleInlet"].read(poiseuilleInlet);
/// Functions returning geometry type parameters
// Name of the geometry STL file name
std::string get_GeoName() const { return geoFile; }
// Sort direction for the STL holes
plint get_SortDir() const { return openingSortDirection; }
// List of inlets and outlets according to SortDirection
std::vector<std::string> get_OpeningType() const { return openingType; }
/// Function returning physical parameters
// Inlet diameter (m)
T get_InletDiam() const { return inletDiam; }
// Average inlet velocity (m/s) to be imposed
T get_u_Physical() const { return u_Physical; }
// return the physical cinematic viscosity (m2/s)
T get_nuPhysical() const { return nu_Physical; }
// return the density (kg/m3)
T get_rhoPhysical() const { return rho_Physical; }
// Return the maximal simulation time
T get_MaxT() const { return simTime; }
// Return the time for increasing inlet velocity (s)
T get_RampeTime() const { return rampeTime; }
// Reynolds number
T get_Re() const { return u_Physical*inletDiam / nu_Physical; }
/// Function returning numerical parameters
// return the cell size (m)
T get_dx() const { return inletDiam/get_resolution(); }
// Return the resolution
plint get_resolution() const { return (int)(inletDiam/cellSize); }
// Return LB velocity (stability criterion)
T get_u_LB() const { return u_LB; }
// Return de time step (s)
T get_dt() const { return get_dx() * (u_LB/u_Physical); }
// Kinematic viscosity in LB units
T get_nu_LB() const { return nu_Physical*get_dt()/get_dx()/get_dx(); }
// Fluid relaxation time in LB
T get_Tau() const { return DESCRIPTOR<T>::invCs2*get_nu_LB()+(T)0.5; }
// Fluid relaxation frequency in LB
T get_Omega() const { return (T)1. / get_Tau(); }
// Returning the cSmago value
T get_cSmago() const { return cSmago; }
// Returning the number of sponge layer cells
plint get_spongeLayer() const { return spongeLayer; }
// Returning the cSmago target value for sponge layer area
T get_cSmagoTarget() const { return cSmagoTarget; }
// Time to process a GIF picture
T get_gifTime() const { return gifTime; }
// Time to process a VTI picture
T get_vtiTime() const { return vtiTime; }
/// Function returning booleans for chosing options
// return the used direction for the extrapolation scheme
bool get_useAllDirection() const { return useAllDirections; }
// return the use Regularized Wall
bool get_useRegularizedWall() const { return useRegularizedWall; }
// return if Poisseuille inlet is used
bool get_poiseuilleInlet() const { return poiseuilleInlet; }
/// Function to write the .log summary file
void writeLog(Param_Pipe const ¶m) {
std::string fullName = global::directories().getLogOutDir() + "log.dat";
std::ofstream ofile(fullName.c_str());
T dx = param.get_dx();
T dt = get_dt();
ofile << "Palabos Elbow pipe summary " << "\n\n";
ofile << " Physical values Number in LB units Number in physical units" << "\n \n";
ofile << "- Reynolds number: " << param.get_Re() << " "
<< param.get_Re() << " - \n";
ofile << "- Kinematic viscosity: " << param.get_nu_LB() << " "
<< param.get_nu_LB() * (dx*dx/dt) << " m2/s \n";
ofile << "- Lattice resolution: " << param.get_resolution() << " "
<< param.get_dx() << "\n";
ofile << "- Time step deltaT: " << "1.0 "
<< param.get_dt() << " s \n";
ofile << "- Solvent omega fluid: " << param.get_Omega() << "\n";
ofile << "- Caracteristic vel: " << param.get_u_LB() << " "
<< param.get_u_Physical() << " m/s \n";
}
private:
std::vectorstd::string openingType;
std::string geoFile;
T inletDiam;
T u_Physical;
T u_LB;
T nu_Physical;
T rho_Physical;
T cellSize;
T simTime;
T rampeTime;
T gifTime;
T vtiTime;
T cSmago;
T cSmagoTarget;
plint openingSortDirection;
plint spongeLayer;
bool useAllDirections;
bool useRegularizedWall;
bool poiseuilleInlet;
};