MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ChargedBondedInteraction.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 
29 #include "Particles/BaseParticle.h"
30 #include "InteractionHandler.h"
31 #include <iomanip>
32 //#include <cassert>
33 
36 
43  : BaseInteraction(P, I, timeStamp)
44 {
45  //ensuring that, by default, particles are not 'bonded'
46  //i.e. they will not unintentionally 'stick' to any overlapping particles!
47  bonded_ = false;
48 #ifdef DEBUG_CONSTRUCTOR
49  std::cout<<"ChargedBondedInteraction::ChargedBondedInteraction() finished"<<std::endl;
50 #endif
51 }
52 
57  : BaseInteraction(p)
58 {
59  //carrying the history parameter over for copied particles to ensure that any bonded particles
60  //remain bonded!
61  bonded_ = p.bonded_;
62 #ifdef DEBUG_CONSTRUCTOR
63  std::cout<<"ChargedBondedInteraction::ChargedBondedInteraction(const ChargedBondedInteraction &p finished"<<std::endl;
64 #endif
65 }
66 
68 {
69 #ifdef MERCURY_USE_MPI
70  logger(FATAL,"ChargedBondedInteractions are currently not implemented in parallel MercuryDPM");
71 #endif
72 }
73 
78 {
79 #ifdef DEBUG_DESTRUCTOR
80  std::cout<<"ChargedBondedInteraction::ChargedBondedInteractionaction() finished"<<std::endl;
81 #endif
82 }
83 
87 void ChargedBondedInteraction::write(std::ostream& os UNUSED) const
88 {
89  os << " bonded " << bonded_;
90 }
91 
95 void ChargedBondedInteraction::read(std::istream& is UNUSED)
96 {
97  std::string dummy;
98  //logger(INFO,"ChargedBondedSpecies %",dummy);
99  is >> dummy >> bonded_;
100 }
101 
106 {
107 
108  const ChargedBondedSpecies* species = getSpecies();
109  //std::cout << getSpecies()->getCharge() << std::endl;
110 
111  //creating local parameters to store the charges of both particles
112  //involved in the interaction to allow for quick calculation
113  const auto pSpecies = dynamic_cast<const ChargedBondedSpecies*>(getP()->getSpecies());
114  const auto iSpecies = dynamic_cast<const ChargedBondedSpecies*>(getI()->getSpecies());
115  logger.assert(pSpecies,"No ChargedBondedSpecies");
116  logger.assert(iSpecies,"No ChargedBondedSpecies");
117  const int pCharge = pSpecies->getCharge();
118  const int iCharge = iSpecies->getCharge();
119 
120  //similarly, creating local parameters to store the relevant stiffness
121  //and max force values
122  const Mdouble k = species->getAdhesionStiffness();
123  const Mdouble fMax = species->getAdhesionForceMax();
124 
125  const Mdouble kWaals = species->getVanDerWaalsStiffness();
126  const Mdouble fMaxWaals = species->getVanDerWaalsForceMax();
127  const Mdouble rWaals = fMaxWaals / kWaals;
128 
129 
130  //First, adding bonded force if applicable
131  if (bonded_ && getOverlap() >= 0)
132  {
133  addForce(getNormal() * (-species->getBondForceMax()
135  return;
136  }
137 
138 
139  //determining which of the three possible cases for force based on charge -
140  //repulsive (like charges), attractive (unlike charges) or none (1 or more uncharged) -
141  //is relevant for the current combination of particle charges...
142  //(Note that the charge set function contains a safety check that means charge can only be
143  // +/- 1, i.e. the expressions used below should not produce errors!
144 
145  //case 1 - 1 or more particles has no charge, i.e. no EM force between them
146  if ((pCharge == 0) or (iCharge == 0))
147  {
148  //No need to write anything here as nothing needs to be returned!
149  //std::cout << "no charge" << std::endl;
150  }
151  //case 2: unlike charges --> attractive force
152  else if (pCharge == -iCharge)
153  {
154  //std::cout << "dissimilar charge" << std::endl;
155  //in the case of particles with opposing charges, simply implements the
156  //standard attractive force used for adhesive particle interactions
157  if (getOverlap() >= 0)
158  {
159  addForce(getNormal() * -fMax);
160  addForce(getNormal() * -fMaxWaals);
161  }
162  else if (getOverlap() >= -rWaals)
163  {
164  addForce(getNormal() * (-kWaals * getOverlap() - fMaxWaals));
165  addForce(getNormal() * (-k * getOverlap() - fMax));
166  }
167  else
168  {
169  addForce(getNormal() * (-k * getOverlap() - fMax));
170  }
171  }
172  //case 3: like charges --> repulsive force
173  else if (pCharge == iCharge)
174  {
175  //std::cout << "similar charge" << std::endl;
176  //in the case of particles with like charges, reverse the direction of the force applied
177  //such that particles repel one another
178  if (getOverlap() >= 0)
179  {
180  addForce(getNormal() * +fMax);
181  addForce(getNormal() * -fMaxWaals);
183  }
184  else if (getOverlap() >= -rWaals)
185  {
186  addForce(getNormal() * (-kWaals * getOverlap() - fMaxWaals));
187  addForce(getNormal() * (+k * getOverlap() + fMax));
188  //std::cout << "Waals = " << getNormal() * (-kWaals * getOverlap() - fMaxWaals) << std::endl;
189  }
190  else
191  {
192  addForce(getNormal() * (+k * getOverlap() + fMax));
193  }
194  }
195  //if none of the above are satisfied, something must have gone very wrong!
196  else
197  {
198  std::cerr << "Error: Particle charge has erroneous value" << std::endl;
199  exit(-1);
200  }
201 }
202 
209 {
210  const ChargedBondedSpecies* species = getSpecies();
211  const auto pSpecies = static_cast<const ChargedBondedSpecies*>(getP()->getSpecies()->getAdhesiveForce());
212  const auto iSpecies = static_cast<const ChargedBondedSpecies*>(getI()->getSpecies()->getAdhesiveForce());
213  logger.assert(pSpecies,"No ChargedBondedSpecies");
214  logger.assert(iSpecies,"No ChargedBondedSpecies");
215  const int pCharge = pSpecies->getCharge();
216  const int iCharge = iSpecies->getCharge();
217 
218  const Mdouble k = species->getAdhesionStiffness();
219  const Mdouble fMax = species->getAdhesionForceMax();
221 
222  const Mdouble kWaals = species->getVanDerWaalsStiffness();
223  const Mdouble fMaxWaals = species->getVanDerWaalsForceMax();
224  const Mdouble rWaals = (fMaxWaals == 0) ? 0 : (fMaxWaals / kWaals);
225 
226 
227  //First, adding bonded force if applicable
228  if (bonded_ && getOverlap() >= 0)
229  {
230  //comment to ignore BondForce
231  Mdouble elasticEnergyAtEquilibrium = getElasticEnergyAtEquilibrium(species->getBondForceMax());
232  return -species->getBondForceMax() * getOverlap() + elasticEnergyAtEquilibrium;
233  }
234 
235  Mdouble elasticEnergy = 0.0;
236  if ((pCharge != 0) && (iCharge != 0))
237  {
238  if (pCharge == -iCharge)
239  {
240  if (getOverlap() >= 0)
241  {
242  elasticEnergy -= (0.5 * rWaals + getOverlap()) * fMaxWaals
243  + (0.5 * r + getOverlap()) * fMax;
244  }
245  else if (getOverlap() >= -rWaals)
246  {
247  elasticEnergy -= (0.5 * kWaals * mathsFunc::square(getOverlap() + rWaals)) +
248  (0.5 * k * mathsFunc::square(getOverlap() + r));
249  }
250  else
251  {
252  elasticEnergy -= (0.5 * k * mathsFunc::square(getOverlap() + r));
253  }
254  }
255  //case 3: like charges --> repulsive force
256  else if (pCharge == iCharge)
257  {
258  if (getOverlap() >= 0)
259  {
260  elasticEnergy -= (0.5 * rWaals + getOverlap()) * fMaxWaals
261  - (0.5 * r + getOverlap()) * fMax;
262  }
263  else if (getOverlap() >= -rWaals)
264  {
265  elasticEnergy -= (0.5 * kWaals * mathsFunc::square(getOverlap() + rWaals)) -
266  (0.5 * k * mathsFunc::square(getOverlap() + r));
267  }
268  else
269  {
270  elasticEnergy += (0.5 * k * mathsFunc::square(getOverlap() + r));
271  }
272  }
273  else
274  {
275  logger(ERROR, "Particle charge has erroneous value");
276  }
277  }
278  return elasticEnergy;
279 }
280 
285 {
286  return static_cast<const ChargedBondedSpecies*> (getBaseSpecies()->getAdhesiveForce()); //downcast
287 }
288 
293 {
294  return "ChargedBonded";
295 }
296 
302 {
303  bonded_ = true;
304 }
305 
312 {
313  bonded_ = false;
314 }
Mdouble getVanDerWaalsForceMax() const
Logger< MERCURY_LOGLEVEL > logger("MercuryKernel")
double Mdouble
Definition: GeneralDefine.h:34
const ParticleSpecies * getSpecies() const
Returns a pointer to the species of this BaseInteractable.
BaseAdhesiveForce * getAdhesiveForce() const
Definition: BaseSpecies.h:152
std::string getBaseName() const
Returns the name of the interaction, see Interaction.h.
Mdouble getAdhesionForceMax() const
Allows the spring constant to be accessed.
void write(std::ostream &os) const override
Interaction print function, which accepts an std::ostream as input.
BaseInteractable * getI()
Returns a pointer to the second object involved in the interaction (often a wall or a particle)...
Mdouble getInteractionDistance() const
returns the largest separation distance at which adhesive short-range forces can occur.
Definition: BaseSpecies.h:146
Stores information about interactions between two interactable objects; often particles but could be ...
ChargedBondedSpecies contains the parameters used to describe a linear reversible short-range force...
const Vec3D & getNormal() const
Gets the normal vector between the two interacting objects.
Mdouble getBondDissipation() const
Allows the additional dissipation used to damp oscillations between bondd particles to be accessed...
BaseInteractable * getP()
Returns a pointer to first object involved in the interaction (normally a particle).
Mdouble getVanDerWaalsStiffness() const
Mdouble getNormalRelativeVelocity() const
Returns a double which is the norm (length) of the relative velocity vector.
const BaseSpecies * getBaseSpecies() const
Return a constant point to BaseSpecies of the interaction.
void read(std::istream &is) override
Interaction read function, which accepts an std::istream as input.
#define UNUSED
Definition: GeneralDefine.h:39
Mdouble getAdhesionStiffness() const
Allows the spring constant to be accessed.
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.
Mdouble getElasticEnergy() const override
Returns the amount of Elastic energy involved in an interaction. Basically used in case you want to w...
const ChargedBondedSpecies * getSpecies() const
Returns a pointer to the adhesive force species ChargedBondedSpecies.
~ChargedBondedInteraction() override
Destructor.
Mdouble getBondForceMax() const
Allows the maximal force for 'bonding' particles together to be accessed.
virtual Mdouble getElasticEnergyAtEquilibrium(Mdouble adhesiveForce) const
void bond()
A pair of functions which can be used to fix or unfix a pair of overlapping particles.
Defines the basic properties that a interactable object can have.
void computeAdhesionForce()
Computes the adhesive forces.
T square(const T val)
squares a number
Definition: ExtendedMath.h:104
bool bonded_
A history parameter to store if the particles were in contact or not. Useful to compute adhesive forc...