Hi all,
I need to implement a DEM simulation involving capillary force between particles. I use a hooke/stiffness contact model.
The capillary force model is the Willett model:
Willett, C. D., Adams, M. J., Johnson, S. A., & Seville, J. P. K. (2000). Capillary bridges between two spherical bodies. Langmuir, 16 (24), 9396–9405.
The capillary force is given by a liquid bridge that is created when two particles come into contact and persist until a rupture distance is reached.
The first step (of a long path) is:
assuming that at each contact a liquid bridge is created and that each bridge has a fixed volume of liquid (no volume of liquid conservation/transfer),
assuming monodisperse system
in this way, Willett's model is defined by three constants: contact angle, surface tension and volume of liquid bridge.
while (not in contact) F_lb=0
while (in contact) F_lb=constant value
if (after contant) and while (distance < rapture distance) F_lb= function of the separation distance
if (distance > rapture distance) F_lb=0
In the manual, I don't find a way to do something like this. The only thing I can imagine is to create liquid bridge as bond and define a force depending on the distance between particles.
while (not in contact) F_lb=0
while (in contact) F_lb=constant value ------------- create bond (fix bond/create??)
if (after contant) and while (distance < rapture distance) F_lb= function of the separation distance ----- model bonding force as function of sep dist
if (distance > rapture distance) ----- delete bond
Can I do something like this?
Does anyone implemented something like this before? Any suggestion?
Thank you.
Giuseppe
aaigner | Wed, 01/28/2015 - 13:34
Add a new cohesion model
Hi Giuseppe,
I would add a new cohesion model: First of all .. check the file cohesion_model_sjkr.h.
You can add your force F_lb during contact and set the cdata.touch flag.
The non-contact force (after an inital contact) can be added in the noCollision function. In addition you set the cdata.touch flag to false only if distance > rapture distance.
That should work.
Cheers
Andreas
giuraso | Sat, 01/31/2015 - 15:37
Thank you Andreas
As you suggested, I added a modified sjkr cohesion model, named liquidbridge. I am able to make something happen only after the collision, I have tested it for two particles so far. The problem now is: in noCollision the force i_forces.delta_F and j_forces.delta_F are non-zero inside the function but the particles experience 0 force. Moreover, cdata.en is 0 0 0 inside noCollision, so I set for fx a value without en[0]
const double fx = Fn_coh;
but the force on the particle is still 0 (I can see it in the dump file).
Following, the noCollision function, (all parameters of the model are declared directly inside the function)
void noCollision(ContactData & cdata, ForceData & i_forces, ForceData & j_forces)
{
const double Rp=2.0*1.0*(1+0.0/2)*(pow(0.040,1.0/3.0)+pow(0.040,2.0/3.0)/10)+cdata.radi+cdata.radj; //Rp=2*S+ri+rj
double r=sqrt(cdata.rsq);
std::cout<<"hi ";
if ((rcg_active())
error->cg(FLERR,"cohesion model liquidbridge");
}
then in global_properties.cpp
ScalarProperty* createNewParam(PropertyRegistry & registry, const char * caller, bool)
{
return createScalarProperty(registry, NewParam, caller);
}
in global_properties.h
ScalarProperty* createNewParam(PropertyRegistry & registry, const char * caller, bool sanity_checks);
and in the input script
fix id all property/global NewParam value
is it correct?
Thank you!
Giuseppe
giuraso | Sun, 02/01/2015 - 11:16
EDIT
Some lines of my previous post are missing.
noCollision function:
void noCollision(ContactData & cdata, ForceData & i_forces, ForceData & j_forces)
{
const double Rp=2.0*1.0*(1+0.0/2)*(pow(0.040,1.0/3.0)+pow(0.040,2.0/3.0)/10)+cdata.radi+cdata.radj; //Rp=2*S+ri+rj
double r=sqrt(cdata.rsq);
std::cout<<"hi ";
if ((r<Rp) && (*cdata.touch!=0)){
*cdata.touch=1;
const double L=sqrt(0.040/1.0); //Vp=40 nl *10^9
double Splus=(r-cdata.radj-cdata.radi)/2.0*1/L;
const double Fn_coh = -2.0*3.14*1.0/(1+2.10*Splus+10.0*pow(Splus,2.0));
std::cout<<"hi Rp "<<Rp<<" Fn "<<Fn_coh<<"en0 "<<cdata.en[0]<<"en1 "<<cdata.en[1]<<"en2 "<<cdata.en[2];
const double fx = Fn_coh;
const double fy = Fn_coh * cdata.en[1];
const double fz = Fn_coh * cdata.en[2];
i_forces.delta_F[0] += fx;
i_forces.delta_F[1] += fy;
i_forces.delta_F[2] += fz;
j_forces.delta_F[0] -= fx;
j_forces.delta_F[1] -= fy;
j_forces.delta_F[2] -= fz;
}
else{
if(cdata.touch) *cdata.touch &= ~TOUCH_COHESION_MODEL;
std::cout<<" rsq "<< cdata.rsq<<" ";
*cdata.touch=0;
std::cout<<" bey";
}
std::cout<<" touch "<< *cdata.touch<<" ";
}
Second question: to add a new parameter and set its value in the input script, can I do something like this:
CohesionModel(LAMMPS * lmp, IContactHistorySetup*) : Pointers(lmp), cohEnergyDens(NULL), newParam(NULL)
{
}
void registerSettings(Settings&) {}
void connectToProperties(PropertyRegistry & registry) {
registry.registerProperty("cohEnergyDens", &MODEL_PARAMS::createCohesionEnergyDensity);
registry.registerProperty("newParam", &MODEL_PARAMS::createNewParam);
registry.connect("cohEnergyDens", cohEnergyDens,"cohesion_model liquidbridge");
registry.connect("newParam", newParam,"cohesion_model liquidbridge");
if(force->cg_active())
error->cg(FLERR,"cohesion model liquidbridge");
}
then in global_properties.cpp
ScalarProperty* createNewParam(PropertyRegistry & registry, const char * caller, bool)
{
return createScalarProperty(registry, NewParam, caller);
}
in global_properties.h
ScalarProperty* createNewParam(PropertyRegistry & registry, const char * caller, bool sanity_checks);
and in the input script
fix id all property/global NewParam value
is it correct?
Thank you!
Giuseppe
ckloss | Mon, 02/09/2015 - 13:44
Hi Giuseppe,
Hi Giuseppe,
collision data is not avalailble in noCollision(), you have to calculate it if you need it
Also, there is a flag which defines if noCollision() adds a force, if you set it to true the force will be updated
Best wishes
Christoph