MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
LiquidMigrationWilletInteraction.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 
30 #include "InteractionHandler.h"
31 #include "DPMBase.h"
32 #include <iomanip>
33 #include <fstream>
34 
41  unsigned timeStamp)
42  : BaseInteraction(P, I, timeStamp)
43 {
44  liquidBridgeVolume_ = 0.0;
45  wasInContact_ = false;
46 #ifdef DEBUG_CONSTRUCTOR
47  std::cout << "LiquidMigrationWilletInteraction::LiquidMigrationWilletInteraction() finished" << std::endl;
48 #endif
49 }
50 
51 //used for mpi
53  : BaseInteraction()
54 {
55  liquidBridgeVolume_ = 0.0;
56  wasInContact_ = false;
57 #ifdef DEBUG_CONSTRUCTOR
58  std::cout << "LiquidMigrationWilletInteraction::LiquidMigrationWilletInteraction() finished" << std::endl;
59 #endif
60 }
61 
66  : BaseInteraction(p)
67 {
70 #ifdef DEBUG_CONSTRUCTOR
71  std::cout << "LiquidMigrationWilletInteraction::LiquidMigrationWilletInteraction(const LiquidMigrationWilletInteraction &p finished" << std::endl;
72 #endif
73 }
74 
79 {
80 #ifdef DEBUG_DESTRUCTOR
81  std::cout << "LiquidMigrationWilletInteraction::~LiquidMigrationWilletInteraction() finished" << std::endl;
82 #endif
83 }
84 
86 {
87  rupture();
88 };
89 
93 void LiquidMigrationWilletInteraction::write(std::ostream& os) const
94 {
95  os
96  << " wasInContact " << wasInContact_
97  << " liquidBridgeVolume " << liquidBridgeVolume_;
98 }
99 
104 {
105  std::string dummy;
106  is >> dummy >> wasInContact_;
107  is >> dummy >> liquidBridgeVolume_;
108 }
109 
114 {
115  // Adding no capillary force for liquid bridge volume = 0
116  if (getLiquidBridgeVolume() == 0) return;
117 
118  if (getOverlap() >= 0)
119  {
120  // if particles are in contact
121  const LiquidMigrationWilletSpecies* species = getSpecies();
122  const Mdouble effectiveRadius = 2.0 * getEffectiveRadius();
123  const Mdouble fdotn = -2.0 * constants::pi * effectiveRadius * species->getSurfaceTension()
124  * std::cos(species->getContactAngle());
125  addForce(getNormal() * fdotn);
126  }
127  else if (wasInContact_)
128  {
129  // if particles are not in contact, but within their interaction distance
130  const LiquidMigrationWilletSpecies* species = getSpecies();
131  const Mdouble effectiveRadius = 2.0 * getEffectiveRadius();
132  const Mdouble s_c = -getOverlap() * std::sqrt(effectiveRadius / getLiquidBridgeVolume());
133  const Mdouble fdotn = -2.0 * constants::pi * effectiveRadius * species->getSurfaceTension()
134  * std::cos(species->getContactAngle()) / (1 + (1.05 + 2.5 * s_c) * s_c);
135  addForce(getNormal() * fdotn);
136  }
137 }
138 
141 {
142  if (wasInContact_)
143  {
144  if (-getOverlap() >= getRuptureDistance())
145  {
146  rupture();
147  }
148  }
149  else
150  {
151  if (getOverlap() >= 0)
152  {
153  form();
154  }
155  }
156 }
157 
159 {
160  //form a bridge
161  //todo: extend to neighbours
162 
163  wasInContact_ = true;
164  const LiquidMigrationWilletSpecies* species = getSpecies();
165  LiquidFilmParticle* IParticle = dynamic_cast<LiquidFilmParticle*>(getI());
166  LiquidFilmParticle* PParticle = dynamic_cast<LiquidFilmParticle*>(getP());
167  if (IParticle == nullptr) //if I is a wall
168  {
169  //do not form bridge if the volume is below minimum
170  if (PParticle->getLiquidVolume() < species->getLiquidBridgeVolumeMin())
171  {
172  return;
173  }
174  //if below max bridge volume, move all liquid from film to volume
175  else if (PParticle->getLiquidVolume() <= species->getLiquidBridgeVolumeMax())
176  {
177  liquidBridgeVolume_ = PParticle->getLiquidVolume();
178  PParticle->setLiquidVolume(0.0);
179  }
180  //if above max bridge volume, fill the liquid bridge and keep the rest of the liquid in the particle
181  else
182  {
184  PParticle->setLiquidVolume(PParticle->getLiquidVolume() - species->getLiquidBridgeVolumeMax());
185  }
186 // if (liquidBridgeVolume_) logger(INFO,"Forming liquid bridge of volume % between particles % and wall %",liquidBridgeVolume_,getP()->getId(),getI()->getId());
187  }
188  else if (PParticle == nullptr) //if P is a wall
189  {
190  logger(ERROR,"Should not happen");
191  }
192  else //if P and I are particles
193  {
194  //if I is a ghost particle, apply volume change only to real particles (this removes the possibility that contacts are established only on one side of the periodic boundary; the same problem could exist for ruptures!)
195  LiquidFilmParticle* IParticleReal;
196  if (IParticle->getPeriodicFromParticle())
197  {
198  IParticleReal = dynamic_cast<LiquidFilmParticle*>(IParticle->getPeriodicFromParticle());
199  }
200  else
201  {
202  IParticleReal = IParticle;
203  }
204 
205  Mdouble distributableLiquidVolume =
206  PParticle->getLiquidVolume() + IParticleReal->getLiquidVolume();
207  //assign all liquid of the contacting particles to the bridge,
208  //if the total volume does not exceed LiquidBridgeVolumeMax
210  if (distributableLiquidVolume <= species->getLiquidBridgeVolumeMin())
211  {
212  return;
213  }
214  else if (distributableLiquidVolume <= species->getLiquidBridgeVolumeMax())
215  {
216  liquidBridgeVolume_ = distributableLiquidVolume;
217  if (!IParticle->getPeriodicFromParticle() ||
218  PParticle->getIndex() < IParticle->getPeriodicFromParticle()->getIndex())
219  {
220  PParticle->setLiquidVolume(0.0);
221  IParticleReal->setLiquidVolume(0.0);
222  }
223  }
224  else //if the total volume exceeds LiquidBridgeVolumeMax, only distribute the max value
225  {
227  Mdouble pFraction =
228  PParticle->getLiquidVolume() / distributableLiquidVolume;
229  if (!IParticle->getPeriodicFromParticle() ||
230  PParticle->getIndex() < IParticle->getPeriodicFromParticle()->getIndex())
231  {
232  PParticle->addLiquidVolume(-pFraction * species->getLiquidBridgeVolumeMax());
233  IParticleReal->addLiquidVolume(-(1.0 - pFraction) * species->getLiquidBridgeVolumeMax());
234  }
235  }
236 // if (liquidBridgeVolume_) logger(INFO,"Forming liquid bridge of volume % between particles % and % (MPI %%, V % %, overlap %)",liquidBridgeVolume_,getP()->getId(),getI()->getId(),PParticle->isMPIParticle(),IParticle->isMPIParticle(),PParticle->getLiquidVolume(),IParticle->getLiquidVolume(),getOverlap());
237  }
238 }
239 
241 {
242  int numContacts = 0;
243  for (auto i : interactable->getInteractions())
244  {
246  LiquidFilmParticle* jIParticle = dynamic_cast<LiquidFilmParticle*>(i->getI());
247  if (j != this && jIParticle != nullptr && j->getLiquidBridgeVolume() != 0.0)
248  numContacts++;
249  }
250  return numContacts;
251 }
252 
254 {
255  // remove the contact history
256  wasInContact_ = false;
257 
258  //if the bridge is already empty, do nothing
259  if (getLiquidBridgeVolume() == 0.0)
260  return;
261 
262  //else rupture a bridge
263  const LiquidMigrationWilletSpecies* species = getSpecies();
264  LiquidFilmParticle* IParticle = dynamic_cast<LiquidFilmParticle*>(getI());
265  LiquidFilmParticle* PParticle = dynamic_cast<LiquidFilmParticle*>(getP());
266  if (IParticle == nullptr) //if I is a wall
267  {
268  int numContactsP = 0;
269  for (auto i : getP()->getInteractions())
270  {
272  if (j != this && j != nullptr && j->getLiquidBridgeVolume() != 0.0)
273  {
274  numContactsP++;
275  }
276 
277  }
278 // logger(INFO,"Rupturing liquid bridge of volume % between particles % and wall % (numContacts %)",liquidBridgeVolume_,getP()->getId(),getI()->getId(),numContactsP);
279  if (numContactsP > 0)
280  {
281  // Updating new volume with distribution less than critical volume
282  Mdouble ExcessBridgeVolume = 0.0;
283  Mdouble newVolume = liquidBridgeVolume_ * species->getDistributionCoefficient() / (numContactsP);
284 
285  for (auto i : getP()->getInteractions())
286  {
288  dynamic_cast<LiquidMigrationWilletInteraction*>(i);
289  if (j != this && j != nullptr && j->getLiquidBridgeVolume() != 0.0)
290  {
292  j->getLiquidBridgeVolume() + newVolume);
293  if (j->getLiquidBridgeVolume() >=
294  species->getLiquidBridgeVolumeMax())
295  {
296  ExcessBridgeVolume +=
298  - species->getLiquidBridgeVolumeMax();
300  species->getLiquidBridgeVolumeMax());
301  }
302  }
303  }
304  PParticle->setLiquidVolume(
305  PParticle->getLiquidVolume() + ExcessBridgeVolume +
307 
309  }
310  else
311  {
312  PParticle->setLiquidVolume(PParticle->getLiquidVolume() + liquidBridgeVolume_);
313  }
314  liquidBridgeVolume_ = 0.0;
315  for (auto i : *getHandler())
316  {
318  dynamic_cast<LiquidMigrationWilletInteraction*>(i);
319  }
320  for (auto i : getHandler()->getDPMBase()->particleHandler)
321  {
322  LiquidFilmParticle* j = dynamic_cast<LiquidFilmParticle*>(i);
323  }
324  }
325  else if (PParticle == nullptr) //if P is a wall
326  {
327  logger(ERROR,"this should not happen");
328  }
329  else //if P and I are particles
330  {
331  //count interaction partners of p (this contact, ghosts and interaction without liquid bridge dont count)
332  int numContactsP = getNumberOfContacts(getP());
333 // logger(INFO,"Rupturing liquid bridge of volume % between particles % and % (numContacts %, MPI %%, overlap %)",liquidBridgeVolume_,getP()->getId(),getI()->getId(),numContactsP,PParticle->isMPIParticle(),IParticle->isMPIParticle(),getOverlap());
334  if (numContactsP < 1) //if P has only one contact (the one that gets ruptured), pass the fluid into it)
335  {
336  PParticle->addLiquidVolume(0.5 * liquidBridgeVolume_);
337  }
338  else //if P has only multiple contacts pass the fluid into it)
339  {
340  // move fluid to neighbouring contacts
341  Mdouble perContactVolume = 0.5 * liquidBridgeVolume_ * species->getDistributionCoefficient() / numContactsP;
342  for (auto i : getP()->getInteractions())
343  {
345  dynamic_cast<LiquidMigrationWilletInteraction*>(i);
346  LiquidFilmParticle* jIParticle =
347  dynamic_cast<LiquidFilmParticle*>(i->getI());
348  if (j != this && jIParticle != nullptr && j->getLiquidBridgeVolume() != 0.0)
349  {
350  Mdouble excessVolume =
351  perContactVolume - (species->getLiquidBridgeVolumeMax() - j->getLiquidBridgeVolume());
352  if (excessVolume < 0.0)
353  {
354  j->addLiquidBridgeVolume(perContactVolume);
355  PParticle->addLiquidVolume(
356  0.5 * liquidBridgeVolume_ * (1 - species->getDistributionCoefficient()) / numContactsP);
357  }
358  else
359  {
361  PParticle->addLiquidVolume(excessVolume);
362  PParticle->addLiquidVolume(
363  0.5 * liquidBridgeVolume_ * (1 - species->getDistributionCoefficient()) / numContactsP);
364  }
365  }
366  }
367  Mdouble PParticle_fin = PParticle->getLiquidVolume();
368 
369  }
370 
371  //count interaction partners of i (ghosts and interaction without liquid bridge dont count)
372  int numContactsI = getNumberOfContacts(getI());
373  if (numContactsI < 1)
374  {
375  IParticle->addLiquidVolume(0.5 * liquidBridgeVolume_);
376  //std::cout << "added " << 0.5 * liquidBridgeVolume_
377  // << " to particle #" << IParticle->getIndex()
378  //<< ", V=" << IParticle->getLiquidVolume() << std::endl;
379  }
380  else
381  {
382  // move fluid to neighbouring contacts
383  Mdouble perContactVolume =
384  0.5 * liquidBridgeVolume_ * species->getDistributionCoefficient() / (numContactsI);
385  for (BaseInteraction* i : getI()->getInteractions())
386  {
388  LiquidFilmParticle* jIParticle = dynamic_cast<LiquidFilmParticle*>(i->getI());
389  if (j != this && jIParticle != nullptr && j->getLiquidBridgeVolume() != 0.0)
390  {
391  Mdouble excessVolume =
392  perContactVolume
393  - species->getLiquidBridgeVolumeMax()
394  + j->getLiquidBridgeVolume();
395  if (excessVolume < 0.0)
396  {
397  j->addLiquidBridgeVolume(perContactVolume);
398  IParticle->addLiquidVolume(
399  0.5 * liquidBridgeVolume_ * (1 - species->getDistributionCoefficient()) /
400  (numContactsI));
401  }
402  else
403  {
404  //std::cout << "excess i-Volume " << excessVolume << std::endl;
406  IParticle->addLiquidVolume(excessVolume);
407  IParticle->addLiquidVolume(
408  0.5 * liquidBridgeVolume_ * (1 - species->getDistributionCoefficient()) /
409  (numContactsI));
410  }
411  //~ std::cout << "added " << perContactVolume
412  //~ << " to i-contact #" << j->getIndex()
413  //~ << ", V=" << j->getLiquidBridgeVolume() << std::endl;
414  }
415  }
416  }
417 
418  //~ std::cout << "ruptured liquid bridge #" << getId()
419  //~ << " between particles #"
420  //~ << PParticle->getIndex()
421  //~ << " and #"
422  //~ << IParticle->getIndex()
423  //~ << ": V" << liquidBridgeVolume_
424  //~ << " -> Vp" << PParticle->getLiquidVolume()
425  //~ << " Vi" << IParticle->getLiquidVolume()
426  //~ << " Np" << numContactsP
427  //~ << " Ni" << numContactsI
428  //~ << std::endl;
429 
430  //to balance the added volume, remove the liquid from the bridge
431  liquidBridgeVolume_ = 0.0;
432 
433  }
434 }
435 
440 {
442  return 0.0;
443 }
444 
449 {
450  return static_cast<const LiquidMigrationWilletSpecies*>(getBaseSpecies()->getAdhesiveForce());
451 ;
452 }
453 
458 {
459  return "LiquidMigrationWillet";
460 }
461 
463 {
464  return liquidBridgeVolume_;
465 }
466 
468 {
469  liquidBridgeVolume_ = liquidBridgeVolume;
470 }
471 
473 {
474  liquidBridgeVolume_ += liquidBridgeVolume;
475 }
476 
478 {
479  return wasInContact_;
480 }
481 
483 {
484  wasInContact_ = wasInContact;
485 }
486 
488 {
489  const LiquidMigrationWilletSpecies* species = getSpecies();
490  return (1.0 + 0.5 * species->getContactAngle()) * cbrt(liquidBridgeVolume_);
491 }
492 
494 {
495  return 1;
496 }
497 
499 {
500  return "Float32";
501 }
502 
504 {
505  return "liquidBridgeRadius";
506 
507 }
508 
509 std::vector<Mdouble> LiquidMigrationWilletInteraction::getFieldVTK(unsigned i) const
510 {
511  return std::vector<Mdouble>(1, cbrt(liquidBridgeVolume_));
512 }
513 
515 {
516  Mdouble volume = 0;
517  for (auto p : h)
518  {
519  auto l = dynamic_cast<LiquidFilmParticle*>(p);
520  if (l)
521  {
522  volume += l->getLiquidVolume();
523  }
524  }
525  return volume;
526 }
527 
529 {
530  Mdouble volume = 0;
531  for (auto i : h)
532  {
533  auto l = dynamic_cast<LiquidMigrationWilletInteraction*>(i);
534  if (l)
535  {
536  volume += l->getLiquidBridgeVolume();
537  }
538  }
539  return volume;
540 }
std::string getTypeVTK(unsigned i) const override
unsigned int getIndex() const
Returns the index of the object in the handler.
Definition: BaseObject.h:118
static Mdouble getTotalLiquidBridgeVolume(InteractionHandler &)
unsigned getNumberOfFieldsVTK() const override
writes extra information to the VTK output
LiquidMigrationWilletSpecies contains the parameters used to describe a short-range force caused by l...
Mdouble getEffectiveRadius() const
Returns a Mdouble to the effective radius of the interaction. (Not corrected for the overlap) ...
int getNumberOfContacts(BaseInteractable *interactable)
Defines the liquid bridge willet interaction between two particles or walls.
static Mdouble getTotalLiquidFilmVolume(ParticleHandler &)
void addLiquidVolume(Mdouble liquidVolume)
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
InteractionHandler * getHandler() const
Gets a point to the interaction handlers to which this interaction belongs.
double Mdouble
Definition: GeneralDefine.h:34
Mdouble getLiquidBridgeVolumeMin() const
used to access the Volume of the liquid bridge.
BaseAdhesiveForce * getAdhesiveForce() const
Definition: BaseSpecies.h:152
Mdouble getElasticEnergy() const override
Returns the amount of Elastic energy involved in an interaction. Basically used in case you want to w...
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
BaseParticle * getPeriodicFromParticle() const
Returns the 'original' particle this one's a periodic copy of.
Definition: BaseParticle.h:338
BaseInteractable * getI()
Returns a pointer to the second object involved in the interaction (often a wall or a particle)...
Mdouble getContactAngle() const
used to access the contact angle between particle and liquid bridge surface.
Stores information about interactions between two interactable objects; often particles but could be ...
Mdouble getLiquidBridgeVolumeMax() const
used to access the Volume of the liquid bridge.
const Vec3D & getNormal() const
Gets the normal vector between the two interacting objects.
std::vector< Mdouble > getFieldVTK(unsigned i) const override
void write(std::ostream &os) const override
Interaction print function, which accepts an std::ostream as input.
bool wasInContact_
A history parameter to store if the particles were in contact or not. Useful to compute adhesive forc...
Mdouble cos(Mdouble x)
Definition: ExtendedMath.cc:64
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
Mdouble getDistributionCoefficient() const
used to access the surface tension of the liquid.
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
void actionsAfterTimeStep() override
test if particle needs to be ruptured
const Mdouble pi
Definition: ExtendedMath.h:45
void actionsOnErase() override
If an interaction needs to do something before it gets erased, add it here. E.g. Liquid bridges ruptu...
ParticleHandler particleHandler
An object of the class ParticleHandler, contains the pointers to all the particles created...
Definition: DPMBase.h:1329
void setLiquidVolume(Mdouble liquidVolume)
Container to store Interaction objects.
std::string getBaseName() const
Returns the name of the interaction, see Interaction.h.
Mdouble getOverlap() const
Returns a Mdouble with the current overlap between the two interacting objects.
void addForce(Vec3D force)
add an force increment to the total force.
std::string getNameVTK(unsigned i) const override
Container to store all BaseParticle.
void addLiquidBridgeVolume(Mdouble liquidBridgeVolume)
Mdouble getLiquidVolume() const
void setLiquidBridgeVolume(Mdouble liquidBridgeVolume)
const std::vector< BaseInteraction * > & getInteractions() const
Returns a list of interactions which belong to this interactable.
Defines the basic properties that a interactable object can have.
void computeAdhesionForce()
Computes the adhesive forces for liquid bridge Willet type of interaction.
DPMBase * getDPMBase()
Gets the problem that is solved using this handler.
Definition: BaseHandler.h:725
Mdouble getSurfaceTension() const
used to access the surface tension of the liquid.
const LiquidMigrationWilletSpecies * getSpecies() const
Returns a pointer to the adhesive force species LiquidMigrationWilletSpecies.