yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3688
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Fri 2015-06-19 15:01:53 +0200
message:
  Implement viscous damping for capillary phase.
modified:
  pkg/dem/ViscoelasticCapillarPM.cpp
  pkg/dem/ViscoelasticCapillarPM.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 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp	2015-04-29 17:12:57 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp	2015-06-19 13:01:53 +0000
@@ -48,6 +48,13 @@
       throw runtime_error("Theta should be equal for both particles!.");
     }
     
+  
+    if (mat1->dcap == mat2->dcap) {
+      phys->dcap = mat1->dcap;
+    } else {
+      throw runtime_error("Theta should be equal for both particles!.");
+    }
+    
     if (mat1->CapillarType == mat2->CapillarType and mat2->CapillarType != ""){
       
       if      (mat1->CapillarType == "Willett_numeric")  {phys->CapillarType = Willett_numeric;  phys->CapFunct = Law2_ScGeom_ViscElCapPhys_Basic::Willett_numeric_f;}
@@ -123,7 +130,32 @@
         VLiqBridg += phys.Vb;
         NLiqBridg += 1;
       }
-      phys.normalForce = -phys.CapFunct(geom, phys)*geom.normal;
+      
+      const auto normalCapForceScalar = phys.CapFunct(geom, phys);
+      Real dampCapForceScalar = 0.0;
+      
+      if (phys.dcap) {
+        const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
+        const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
+      
+        Vector3r& shearForce = phys.shearForce;
+        if (I->isFresh(scene)) shearForce=Vector3r(0,0,0);
+        shearForce = geom.rotate(shearForce);
+        
+        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 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 auto normalVelocity	= geom.normal.dot(relativeVelocity);
+        
+        dampCapForceScalar = -phys.dcap * normalVelocity;
+      }
+      
+      phys.normalForce = -(normalCapForceScalar + dampCapForceScalar)*geom.normal;
+      
       if (I->isActive) {
         addForce (id1,-phys.normalForce,scene);
         addForce (id2, phys.normalForce,scene);

=== modified file 'pkg/dem/ViscoelasticCapillarPM.hpp'
--- pkg/dem/ViscoelasticCapillarPM.hpp	2014-12-10 10:52:48 +0000
+++ pkg/dem/ViscoelasticCapillarPM.hpp	2015-06-19 13:01:53 +0000
@@ -14,6 +14,7 @@
 		((Real,Vb,0.0,,"Liquid bridge volume [m^3]"))
 		((Real,gamma,0.0,,"Surface tension [N/m]"))
 		((Real,theta,0.0,,"Contact angle [°]"))
+		((Real,dcap,0.0,,"Damping coefficient for the capillary phase [-]"))
 		((std::string,CapillarType,"",,"Different types of capillar interaction: Willett_numeric, Willett_analytic [Willett2000]_ , Weigert [Weigert1999]_ , Rabinovich [Rabinov2005]_ , Lambert (simplified, corrected Rabinovich model) [Lambert2008]_ ")),
 		createIndex();
 	);
@@ -38,6 +39,7 @@
 		((Real,gamma,0.0,,"Surface tension [N/m]"))
 		((Real,theta,0.0,,"Contact angle [rad]"))
 		((CapType,CapillarType,None_Capillar,,"Different types of capillar interaction: Willett_numeric, Willett_analytic, Weigert, Rabinovich, Lambert, Soulie"))
+		((Real,dcap,0.0,,"Damping coefficient for the capillary phase [-]"))
 #ifdef YADE_LIQMIGRATION
 		((bool,LiqMigrEnabled,,,"True, if liquid migration mechanism is needed."))
 		((Real,Vmax,0.0,,"Maximal liquid bridge volume [m^3]"))