yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3492
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Mon 2014-10-20 15:42:59 +0200
message:
  Reintroduce viscosity in SPH-modell.
modified:
  doc/references.bib
  pkg/common/SPHEngine.cpp
  pkg/dem/ViscoelasticPM.cpp
  pkg/dem/ViscoelasticPM.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-10-17 17:14:29 +0000
+++ doc/references.bib	2014-10-20 13:42:59 +0000
@@ -826,3 +826,16 @@
   pages = {543-574},
   doi = {10.1146/annurev.aa.30.090192.002551}
 }
+
+@article{Morris1997,
+  title = "Modeling Low Reynolds Number Incompressible Flows Using \{SPH\} ",
+  journal = "Journal of Computational Physics ",
+  volume = "136",
+  number = "1",
+  pages = "214 - 226",
+  year = "1997",
+  note = "",
+  issn = "0021-9991",
+  doi = "http://dx.doi.org/10.1006/jcph.1997.5776";,
+  url = "http://www.sciencedirect.com/science/article/pii/S0021999197957764";,
+}

=== modified file 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp	2014-10-17 17:14:29 +0000
+++ pkg/common/SPHEngine.cpp	2014-10-20 13:42:59 +0000
@@ -112,8 +112,24 @@
   
   const BodyContainer& bodies = *scene->bodies;
   
+  //////////////////////////////////////////////////////////////////
+  // Copy-paste
+  // Handle periodicity.
+  const Vector3r shift2 = scene->isPeriodic ? scene->cell->intrShiftPos(I->cellDist): Vector3r::Zero(); 
+  const Vector3r shiftVel = scene->isPeriodic ? scene->cell->intrShiftVel(I->cellDist): Vector3r::Zero(); 
+  
   const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
   const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
+
+  const Vector3r c1x = (geom.contactPoint - de1.pos);
+  const Vector3r c2x = (geom.contactPoint - de2.pos - shift2);
+  
+  const Vector3r relativeVelocity = (de1.vel+de1.angVel.cross(c1x)) - (de2.vel+de2.angVel.cross(c2x)) + shiftVel;
+  const Real normalVelocity	= geom.normal.dot(relativeVelocity);
+  
+  // Copy-paste  
+  //////////////////////////////////////////////////////////////////
+  
   
   const Real Mass1 = bodies[id1]->state->mass;
   const Real Mass2 = bodies[id2]->state->mass;
@@ -125,7 +141,7 @@
   
   if (xixj.norm() < phys.h) {
     Real fpressure = 0.0;
-    if (bodies[id1]->rho!=0.0 and bodies[id2]->rho!=0.0) {
+    if (Rho1!=0.0 and Rho2!=0.0) {
       // from [Monaghan1992], (3.3), multiply by Mass2, because we need a force, not du/dt
       fpressure = - Mass1 * Mass2 * (
                   bodies[id1]->press/(Rho1*Rho1) + 
@@ -134,7 +150,16 @@
                   * phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
     }
     
-    force = fpressure*geom.normal;
+    Vector3r fvisc = Vector3r::Zero();
+    if (Rho1!=0.0 and Rho2!=0.0) {
+      // from [Morris1997], (22), multiply by Mass2, because we need a force, not du/dt
+      fvisc = phys.mu * Mass1 * Mass2 *
+      (-normalVelocity*geom.normal)/(Rho1*Rho1) *
+      1 / (xixj.norm()) *
+      phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
+      //phys.kernelFunctionCurrentVisco(xixj.norm(), phys.h);
+    }
+    force = fpressure*geom.normal + fvisc;
     return true;
   } else {
     return false;

=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2014-10-17 17:14:29 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2014-10-20 13:42:59 +0000
@@ -255,7 +255,7 @@
 #ifdef YADE_SPH
 	if (mat1->SPHmode and mat2->SPHmode)  {
 		phys->SPHmode=true;
-		phys->mu=(mat1->mu+mat2->mu)/2.0;
+		phys->mu=(mat1->mu+mat2->mu);
 		phys->h=(mat1->h+mat2->h)/2.0;
 	}
 	

=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp	2014-10-17 17:14:29 +0000
+++ pkg/dem/ViscoelasticPM.hpp	2014-10-20 13:42:59 +0000
@@ -34,7 +34,7 @@
 		((Real,mR,0.0,,"Rolling resistance, see [Zhou1999536]_."))
 #ifdef YADE_SPH
 		((bool,SPHmode,false,,"True, if SPH-mode is enabled."))
-		((Real,mu,-1,, "Viscosity. See Mueller [Mueller2003]_ ."))                                              // [Mueller2003], (14)
+		((Real,mu,-1,, "Viscosity. See Mueller [Morris1997]_ ."))                                              // [Mueller2003], (14)
 		((Real,h,-1,,  "Core radius. See Mueller [Mueller2003]_ ."))                                            // [Mueller2003], (1)
 		((int,KernFunctionPressure,Lucy,, "Kernel function for pressure calculation (by default - Lucy). The following kernel functions are available: Lucy=1."))
 		((int,KernFunctionVisco,   Lucy,, "Kernel function for viscosity calculation (by default - Lucy). The following kernel functions are available: Lucy=1."))