yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 4119
author: Nolan Dyck <ndyck@xxxxxx>
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Mon 2014-07-28 08:24:36 +0200
message:
  Update BubbleMat.hpp and BubbleMat.cpp
added:
  examples/bouncingbubble.py
modified:
  doc/references.bib
  pkg/dem/BubbleMat.cpp
  pkg/dem/BubbleMat.hpp


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'doc/references.bib'
--- doc/references.bib	2014-05-08 05:57:04 +0000
+++ doc/references.bib	2014-07-28 06:24:36 +0000
@@ -34,6 +34,16 @@
 	year = "2005"
 }
 
+@Article{ Chan2011,
+	title = "Film drainage and coalescence between deformable drops and bubbles.",
+	author = "D. Chan and E. Klaseboer and R. Manica",
+	journal = "Soft Matter",
+	pages = "2235-2264",
+	volume = "7",
+	number = "6",
+	year = "2011"
+}
+
 @Article{ Chareyre2002a,
 	title = "Theoretical versus experimental modeling of the anchorage capacity of geotextiles in trenches.",
 	author = "B. Chareyre and L. Briancon and P. Villard",

=== added file 'examples/bouncingbubble.py'
--- examples/bouncingbubble.py	1970-01-01 00:00:00 +0000
+++ examples/bouncingbubble.py	2014-07-28 06:24:36 +0000
@@ -0,0 +1,18 @@
+
+rad = 5e-1
+#O.materials.append(FrictMat(young=1e3,density=1000))
+O.materials.append(BubbleMat(density=1))
+O.bodies.append([
+   utils.sphere(center=(0,0,0),radius=rad,fixed=True),
+   utils.sphere((0,0,2*rad*1.1),rad)
+])
+O.dt = 1e-6
+
+O.engines=[
+   ForceResetter(),
+   InsertionSortCollider([Bo1_Sphere_Aabb()]),
+   InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_BubbleMat_BubbleMat_BubblePhys()],[Law2_ScGeom_BubblePhys_Bubble(surfaceTension=0.035)]),
+#   InteractionLoop([Ig2_Sphere_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),
+   NewtonIntegrator(damping=0.1,gravity=(0,0,-9.81e-2))
+]
+O.saveTmp()

=== modified file 'pkg/dem/BubbleMat.cpp'
--- pkg/dem/BubbleMat.cpp	2014-07-25 16:12:46 +0000
+++ pkg/dem/BubbleMat.cpp	2014-07-28 06:24:36 +0000
@@ -13,40 +13,49 @@
 
 	shared_ptr<BubblePhys> phys(new BubblePhys());
 	interaction->phys = phys;
-	BubbleMat* mat1 = YADE_CAST<BubbleMat*>(m1.get());
-	BubbleMat* mat2 = YADE_CAST<BubbleMat*>(m2.get());
-
-	// averaging over both materials
-	phys->surfaceTension = .5*(mat1->surfaceTension + mat2->surfaceTension);
 }
 
-
 /********************** BubblePhys ****************************/
 CREATE_LOGGER(BubblePhys);
-Real BubblePhys::computeForce(Real penetrationDepth, Real surfaceTension, Real rAvg, int newtonIter, Real newtonTol) {
-	if (penetrationDepth <= 0.0) { return 0.0; }
-
-	Real f,df,ll,retOld,residual;
-	Real c1 = 1./(4*Mathr::PI*surfaceTension);
-	Real c2 = 1./(8*Mathr::PI*surfaceTension*rAvg);
-	Real ret=1./c2;
-	int i = 0;
-	do {
-		retOld = ret;
-		ll = log( ret*c2 );
-		f = penetrationDepth - ret*c1*ll;
-		df = -c1*(ll + 1);
-		ret -= f/df;
-		residual = std::abs(ret - retOld)/retOld;
-		if (i++ > newtonIter) {
-			throw runtime_error("BubblePhys::computeForce: no convergence\n");
-		}
-	} while (residual > newtonTol);
-	return ret;
-}
-
-
-
+void BubblePhys::computeCoeffs(Real pctMaxForce, Real surfaceTension, Real c1)
+{
+	    Real Fmax = pctMaxForce*c1*rAvg;
+	    Real logPct = log(pctMaxForce/4);
+	    Dmax = pctMaxForce*rAvg*logPct;
+	    Real dfdd = c1/(logPct+1);
+	    coeffB = dfdd/Fmax;
+	    coeffA = Fmax/exp(coeffB*Dmax);
+	    fN = 0.1*Fmax;
+}
+
+Real BubblePhys::computeForce(Real separation, Real surfaceTension, Real rAvg, int newtonIter, Real newtonTol, Real c1, Real fN, BubblePhys* phys) {
+
+	if (separation >= phys->Dmax) {
+	  Real f,df,g,retOld,residual;
+	  Real c2 = 1./(4*c1*rAvg);
+	  Real ret = fN;
+	  int i = 0;
+	  do {				//[Chan2011], Loop solves modified form of equation (25) using Newton-Raphson method
+		  retOld = ret;
+		  g = log(ret*c2);
+		  f = separation*c1 - ret*g;
+		  df = g+1;
+		  ret += f/df;
+		  if(ret <= 0.0){	//Need to make sure ret > 0, otherwise the next iteration will try to evaluate the natural logarithm of a negative number, which results in NaN
+		    ret = 0.9*fabs(ret);
+		    residual = newtonTol*2;	//Also, if, by chance, retOld = 0.9*ret it would cause the loop to exit based on the boolean residual > newtonTol, so we force the next iteration by setting residual = newtonTol*2
+		  }	
+		  else {residual = fabs(ret - retOld)/retOld;}
+		  if (i++ > newtonIter) {
+			  throw runtime_error("BubblePhys::computeForce: no convergence\n");
+		  }
+	  } while (residual > newtonTol);
+	  return ret;
+	}
+	else {				//Artificial Extension of [Chan2011] equation 25 to approximiate behaviour outside the valid regime (large penetration cases)
+	  return phys->coeffA*exp(phys->coeffB*separation);
+	}
+}
 
 /********************** Law2_ScGeom_BubblePhys_Bubble ****************************/
 CREATE_LOGGER(Law2_ScGeom_BubblePhys_Bubble);
@@ -54,13 +63,21 @@
 bool Law2_ScGeom_BubblePhys_Bubble::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
 	ScGeom* geom=static_cast<ScGeom*>(_geom.get());
 	BubblePhys* phys=static_cast<BubblePhys*>(_phys.get());
-
+	
+	if(geom->penetrationDepth <= 0.0) {
+		scene->interactions->requestErase(I);
+		return;
+	}
+	
 	if (I->isFresh(scene)) {
+		c1 = 2*Mathr::PI*surfaceTension;
 		phys->rAvg = .5*(geom->refR1 + geom->refR2);
+		phys->computeCoeffs(pctMaxForce,surfaceTension,c1);
 	}
 
 	Real &fN = phys->fN;
-	fN = phys->computeForce(geom->penetrationDepth, phys->surfaceTension, phys->rAvg, phys->newtonIter, phys->newtonTol);
+	fN = phys->computeForce(-geom->penetrationDepth, surfaceTension, phys->rAvg, phys->newtonIter, phys->newtonTol, c1, fN, phys);
+	phys->fN = fN;
 	Vector3r &normalForce = phys->normalForce;
 	normalForce = fN*geom->normal;
 

=== modified file 'pkg/dem/BubbleMat.hpp'
--- pkg/dem/BubbleMat.hpp	2014-07-18 18:18:50 +0000
+++ pkg/dem/BubbleMat.hpp	2014-07-28 06:24:36 +0000
@@ -15,7 +15,7 @@
 		((Real,surfaceTension,0.07197,,"The surface tension in the fluid surrounding the bubbles. The default value is that of water at 25 degrees Celcius."))
 		,
 		createIndex();
-		density=1; // TODO density default value
+		density=1000;
 	);
 	REGISTER_CLASS_INDEX(BubbleMat,Material);
 };
@@ -24,9 +24,14 @@
 
 /********************** BubblePhys ****************************/
 class BubblePhys : public IPhys {
+	private:
+
+	Real coeffA,coeffB; //Coefficents for artificial curve
+  
 	public:
 
-	static Real computeForce(Real penetrationDepth, Real surfaceTension, Real rAvg, int newtonIter, Real newtonTol);
+	void computeCoeffs(Real pctMaxForce,Real surfaceTension, Real c1);
+	static Real computeForce(Real separation, Real surfaceTension, Real rAvg, int newtonIter, Real newtonTol, Real c1, Real fN, BubblePhys* phys);
 
 	virtual ~BubblePhys(){};
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(BubblePhys,IPhys,"Physics of bubble-bubble interactions, for use with BubbleMat",
@@ -34,6 +39,7 @@
 		((Real,surfaceTension,NaN,,"Surface tension of the surrounding liquid"))
 		((Real,fN,NaN,,"Contact normal force"))
 		((Real,rAvg,NaN,,"Average radius of the two interacting bubbles"))
+		((Real,Dmax,NaN,,"Maximum penetrationDepth of the bubbles before the force displacement curve changes to an artificial exponential curve. Setting this value will have no effect. See Law2_ScGeom_BubblePhys_Bubble::pctMaxForce for more information"))
 		((Real,newtonIter,50,,"Maximum number of force iterations allowed"))
 		((Real,newtonTol,1e-6,,"Convergence criteria for force iterations"))
 		,
@@ -63,10 +69,16 @@
 
 /********************** Law2_ScGeom_BubblePhys_Bubble ****************************/
 class Law2_ScGeom_BubblePhys_Bubble : public LawFunctor{
+	private:
+	  
+	  Real c1; //Coeff used for many contacts
+  
 	public:
 	bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* interaction);
 	FUNCTOR2D(GenericSpheresContact,BubblePhys);
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_BubblePhys_Bubble,LawFunctor,"Constitutive law for Bubble model.",
+		((Real,pctMaxForce,0.1,,"Chan[2011] states the contact law is valid only for small interferences; therefore an exponential force-displacement curve models the contact stiffness outside that regime (large penetration). This artificial stiffening ensures that bubbles will not pass through eachother or completely overlap during the simulation. The maximum force is Fmax = (2*pi*surfaceTension*rAvg). pctMaxForce is the percentage of the maximum force dictates the separation threshold, Dmax, for each contact. Penetrations less than Dmax calculate the reaction force from the derived contact law, while penetrations equal to or greater than Dmax calculate the reaction force from the artificial exponential curve."))
+		((Real,surfaceTension,0.07197,,"The surface tension in the liquid surrounding the bubbles. The default value is that of water at 25 degrees Celcius."))
 		,
 		/*ctor*/,
 	);