yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3505
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Fri 2014-10-24 10:39:43 +0200
message:
  Change sign in calculation of viscosity in SPH
modified:
  pkg/common/SPHEngine.cpp


--
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 'pkg/common/SPHEngine.cpp'
--- pkg/common/SPHEngine.cpp	2014-10-23 13:46:20 +0000
+++ pkg/common/SPHEngine.cpp	2014-10-24 08:39:43 +0000
@@ -197,31 +197,27 @@
   const Vector3r xixj = de2.pos - de1.pos;
   
   if ( phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h)) {
-    if (xixj.norm() < phys.h) {
-      Real fpressure = 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]->state->press/(Rho1*Rho1) + 
-                    bodies[id2]->state->press/(Rho2*Rho2) 
-                    )
-                    * phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
-      }
-      
-      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*Rho2) *
-        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;
-    }
+    Real fpressure = 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]->state->press/(Rho1*Rho1) + 
+                  bodies[id2]->state->press/(Rho2*Rho2) 
+                  )
+                  * phys.kernelFunctionCurrentPressure(xixj.norm(), phys.h);
+    }
+    
+    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*Rho2) *
+      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;
   }