I am attempting to model a liquid drop on a solid surface and modify the contact angles by varying the solid node densities. There was a post by jlatt which indicated that the solid node densities could be modified by using the following command (solid density of 0.5 for example):
new BounceBack<T,Descriptor>(0.5);
I have also seen the following command used to set up the boundary node solid densities (on all sides and using 0.8 for example) at another post:
setBoundaryDensity(lattice, lattice.getBoundingBox(), 0.8 );
I have attempted to enter this into my code but have not been able to do so correctly. I get compilation errors because I am not sure which approach to use and where to insert it correctly. I am using the Shan-Chen processor and external moment BGK dynamics.
Perhaps it has not been implemented correctly. Any suggestions as to how to properly implement this would be greatly appreciated! The code is provided below for reference.
Thanks,
Dan
#include "palabos2D.h"
#include "modified_code/shanChenProcessor2D.hh" //Include modified .hh files before palabos2D.hh
#include "palabos2D.hh"
#include <cstdlib>
#include <iostream>
#include <iomanip>
using namespace plb;
using namespace std;
typedef double T;
///////////////////////////////////////////User input parameters//////////////////////////////////////////
// Lattice geometry
#define DESCRIPTOR descriptors::ForcedShanChenD2Q9Descriptor
const plint nx = 200; // x lattice dimension
const plint ny = 200; // y lattice dimension
//Constant accel due to an external body force (gravity) acting on the fluid
const T accel = 0.0000005; //this value can be set to 0 to turn this off.
//LBM parameters
const T omega = 1.; // Choice of the relaxation parameter
//Interparticle potential
// For the choice of the parameters G, rho0, and psi0, we refer to the book
// Michael C. Sukop and Daniel T. Thorne (2006),
// Lattice Boltzmann Modeling; an Introduction for Geoscientists and Engineers.
// Springer-Verlag Berlin/Heidelberg.
const T G = -120.0;
const T rho0 = 200.0;
const T psi0 = 4.0;
//Initial conditions
const T RhoAll = 86; // All cells have a uniform initially density to represent a vapor or a liquid ...
const T RhoCenter = 525; //Density in the liquid blob
const plint radius = 80;
const plint centerX = (nx/2); //centered in the x-direction
const plint centerY = radius; //touching the lower boundary
Array<T,2> u0(0,0); //Initial zero velocity; to be applied at all lattice cells
//Output directory
char opDir[]="/bessemer/dlr28/palabos/LiqBlob/tmp5000/";
///////////////////////////// End of user input; code starts here //////////////////////////////////////////////
// Initialize the lattice at zero velocity and constant density, except
// for a liquid blob on a circular sub-domain.
void initializeBlob(plint iX, plint iY, T& rho, Array<T,2>& u) {
u = u0;
if( (iX-centerX)*(iX-centerX) + (iY-centerY)*(iY-centerY) < radius*radius) {
rho = RhoCenter;
}
//use this when doing the liquid drop on a surface
else{
rho=RhoAll;
}
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Write a VTK file which can be post-processed for example in ParaView.
void writeVTK(MultiBlockLattice2D<T,DESCRIPTOR>& lattice,plint iter)
{
VtkImageOutput2D<T> vtkOut(createFileName("vtk", iter, 7));
vtkOut.writeData<float>(*computeDensity(lattice), "density", 1.);
}
///////////////////////////// Main program ///////////////////////////////////////////////////////////////////
int main(int argc, char *argv[])
{
plbInit(&argc, &argv);
plint iniT, maxIter, saveIter, statIter, checkPointIter, solidRho;
if (argc != 7) {
pcout << "Error; the syntax is \"" << argv[0] << " start-iter end-iter save-iter stat-iter checkpoint-iter solid-node-density\"," << endl;
}
stringstream iniTstr, endTstr, savestr, statstr, checkPointstr, solidRhostr;
iniTstr << argv[1]; iniTstr >> iniT;// start iteration here - 0 for begining >0 for loading prior simulation
endTstr << argv[2]; endTstr >> maxIter;// Iterate this many steps.
savestr << argv[3]; savestr >> saveIter;// start iteration here
statstr << argv[4]; statstr >> statIter;// Iterate this many steps.
checkPointstr << argv[5]; checkPointstr >> checkPointIter; //create saved iteration state for future loading
solidRhostr << argv[6]; solidRhostr >> solidRho;// Solid node density to control contact angles
//Output directory setup
//======================
global::directories().setOutputDir(opDir);
// Lattice creation
// ================
MultiBlockLattice2D<T, DESCRIPTOR> lattice (
nx,ny, new ExternalMomentBGKdynamics<T, DESCRIPTOR>(omega));
//Setting the solid boundary node density
setBoundaryDensity(lattice, Box2D(0,nx-1, 0,0), solidRho);
//Boundary conditions
//===================
//Periodic horizontal boundary condition
lattice.periodicity().toggle(0,true);
// Initialization
//===============
// Initialize the f's to be feq's at the initial density and velocity fields
initializeAtEquilibrium (lattice, lattice.getBoundingBox(), initializeBlob);
// Let's have gravity acting on the fluid.
setExternalVector(lattice, lattice.getBoundingBox(), DESCRIPTOR<T>::ExternalField::forceBeginsAt, Array<T,2>(0.,-accel));
// Add the data processor which implements the Shan/Chen interaction potential.
plint processorLevel = 1;
integrateProcessingFunctional (
new ShanChenSingleComponentProcessor2D<T,DESCRIPTOR> (
G, new interparticlePotential::PsiShanChen94<T>(psi0,rho0) ),
lattice.getBoundingBox(), lattice, processorLevel );
lattice.initialize();
// Load saved data from a previous simulation, if the initial time step
// is larger than zero.
if (iniT>0) {
//loadRawMultiBlock(lattice, "checkpoint.dat");
loadBinaryBlock(lattice, "/bessemer/dlr28/palabos/LiqBlob/checkpoint.dat");
}
// Run the simulation
// ==================
pcout << "Starting simulation" << endl;
pcout << "Timestep Density ( Avg, Min, Max)"<<endl; //AJ Changed output format
for (int iT=0; iT<maxIter; ++iT) {
if (iT%statIter==0) {
auto_ptr<MultiScalarField2D<T> > rho( computeDensity(lattice) );
pcout << setw(8)<< iT <<setw(12)<< computeAverage(*rho) << setw(12) << computeMin(*rho) <<setw(12)<< computeMax(*rho) << endl;
}
if (iT%saveIter==0) {
ImageWriter<T>("leeloo.map").writeScaledGif (
createFileName("rho", iT, 6), *computeDensity(lattice) );
writeVTK(lattice,iT);
}
if (iT%checkPointIter==0 && iT>=iniT) {
pcout << "Saving the state of the simulation ..." << endl;
//saveRawMultiBlock(lattice, "checkpoint.dat");
saveBinaryBlock(lattice, "/bessemer/dlr28/palabos/LiqBlob/tmp5000/checkpoint.dat");
}
lattice.collideAndStream();
}
}