yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10792
[Branch ~yade-pkg/yade/git-trunk] Rev 3942: Implement Liquid Migration model (experimental).
------------------------------------------------------------
revno: 3942
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Thu 2014-05-08 07:57:04 +0200
message:
Implement Liquid Migration model (experimental).
According to [Mani2013]. All relevant code is screened by
YADE_LIQMIGRATION proprocessor directieves, so it will not
hurt productivity of the builds, where this feature
will be disabled.
The feature is disabled by default.
added:
examples/capillary/liquidmigration/
examples/capillary/liquidmigration/showcase.py
modified:
CMakeLists.txt
core/Body.hpp
core/Scene.cpp
core/Scene.hpp
doc/references.bib
pkg/dem/VTKRecorder.cpp
pkg/dem/VTKRecorder.hpp
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 'CMakeLists.txt'
--- CMakeLists.txt 2014-05-05 18:38:31 +0000
+++ CMakeLists.txt 2014-05-08 05:57:04 +0000
@@ -16,6 +16,7 @@
# ENABLE_PFVFLOW: enable PFVFLOW-option, FlowEngine (ON by default)
# ENABLE_SPH: enable SPH-option, Smoothed Particle Hydrodynamics (OFF by default, experimental)
# ENABLE_LBMFLOW: enable LBMFLOW-option, LBM_ENGINE (OFF by default)
+# ENABLE_LIQMIGRATION: enable LIQCONTROL-option, see [Mani2013] for details (OFF by default)
# runtimePREFIX: used for packaging, when install directory is not the same is runtime directory (/usr/local by default)
# CHUNKSIZE: set >1, if you want several sources to be compiled at once. Increases compilation speed and RAM-consumption during it (1 by default).
@@ -120,6 +121,7 @@
OPTION(ENABLE_PFVFLOW "Enable flow engine (experimental)" ${DEFAULT})
OPTION(ENABLE_SPH "Enable SPH" OFF)
OPTION(ENABLE_LBMFLOW "Enable LBM engine (very experimental)" OFF)
+OPTION(ENABLE_LIQMIGRATION "Enable liquid control (very experimental), see [Mani2013] for details" OFF)
#===========================================================
# Use Eigen3 by default
@@ -297,6 +299,13 @@
SET(DISABLED_FEATS "${DISABLED_FEATS} SPH")
ENDIF(ENABLE_SPH)
#===============================================
+IF(ENABLE_LIQMIGRATION)
+ ADD_DEFINITIONS("-DYADE_LIQMIGRATION")
+ SET(CONFIGURED_FEATS "${CONFIGURED_FEATS} LIQCONTROL")
+ELSE(ENABLE_LIQMIGRATION)
+ SET(DISABLED_FEATS "${DISABLED_FEATS} LIQCONTROL")
+ENDIF(ENABLE_LIQMIGRATION)
+#===============================================
IF(ENABLE_GL2PS)
FIND_PACKAGE(GL2PS)
IF(GL2PS_FOUND)
=== modified file 'core/Body.hpp'
--- core/Body.hpp 2014-04-09 14:03:16 +0000
+++ core/Body.hpp 2014-05-08 05:57:04 +0000
@@ -94,6 +94,10 @@
((Real,press,0.0,, "Pressure (only for SPH-model)")) // [Mueller2003], (12)
((Real,Cs,0.0,, "Color field (only for SPH-model)")) // [Mueller2003], (15)
#endif
+#ifdef YADE_LIQMIGRATION
+ ((Real,Vf, -1.0,, "Individual amount of liquid"))
+ ((Real,Vmin,-1.0,, "Minimal amount of liquid"))
+#endif
,
/* ctor */,
/* py */
=== modified file 'core/Scene.cpp'
--- core/Scene.cpp 2013-08-23 15:21:20 +0000
+++ core/Scene.cpp 2014-05-08 05:57:04 +0000
@@ -203,5 +203,6 @@
}
}
bound->min=mn;
- bound->max=mx;
+ bound->max=mx;
}
+
=== modified file 'core/Scene.hpp'
--- core/Scene.hpp 2013-10-02 12:47:38 +0000
+++ core/Scene.hpp 2014-05-08 05:57:04 +0000
@@ -23,13 +23,23 @@
#ifndef HOST_NAME_MAX
#define HOST_NAME_MAX 255
#endif
-
+#ifdef YADE_OPENMP
+ #include<omp.h>
+#endif
class Bound;
#ifdef YADE_OPENGL
class OpenGLRenderer;
#endif
+#ifdef YADE_LIQMIGRATION
+struct intReal {
+ public:
+ id_t id;
+ Real Vol;
+};
+#endif
+
class Scene: public Serializable{
public:
//! Adds material to Scene::materials. It also sets id of the material accordingly and returns it.
@@ -60,6 +70,11 @@
shared_ptr<Engine> engineByName(const string& s);
+ #ifdef YADE_LIQMIGRATION
+ OpenMPVector<Interaction* > addIntrs; //Array of added interactions, needed for liquid migration.
+ OpenMPVector<intReal > delIntrs; //Array of deleted interactions, needed for liquid migration.
+ #endif
+
#ifdef YADE_OPENGL
shared_ptr<OpenGLRenderer> renderer;
#endif
=== modified file 'doc/references.bib'
--- doc/references.bib 2014-04-15 13:55:10 +0000
+++ doc/references.bib 2014-05-08 05:57:04 +0000
@@ -709,3 +709,19 @@
keywords = {liquid bridge, polydisperse, cohesion, capillary force, humid granular media, DEM simulation},
year = {2006},
}
+
+@article{Mani2013,
+ year={2013},
+ issn={1434-5021},
+ journal={Granular Matter},
+ volume={15},
+ number={4},
+ doi={10.1007/s10035-012-0387-3},
+ title={Liquid migration in sheared unsaturated granular media},
+ url={http://dx.doi.org/10.1007/s10035-012-0387-3},
+ publisher={Springer Berlin Heidelberg},
+ keywords={Wet granular matter; Contact dynamics simulations; Liquid bridge; Cohesion; Liquid migration},
+ author={Mani, Roman and Kadau, Dirk and Herrmann, HansJ.},
+ pages={447-454},
+ language={English}
+}
=== added directory 'examples/capillary/liquidmigration'
=== added file 'examples/capillary/liquidmigration/showcase.py'
--- examples/capillary/liquidmigration/showcase.py 1970-01-01 00:00:00 +0000
+++ examples/capillary/liquidmigration/showcase.py 2014-05-08 05:57:04 +0000
@@ -0,0 +1,126 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+# The script implements the show case in this article [Mani2013]
+from yade import utils, plot
+o = Omega()
+fr = 0.5;rho=2000
+tc = 0.001; en = 0.7; et = 0.7;
+o.dt = 1.0
+
+
+r1 = 1.0
+r2 = 1.0
+Gamma = 20.6*1e-3
+Theta = 0
+VB3 = 74.2*1e-12
+
+
+CapillarType = "Lambert"
+
+mat1 = O.materials.append(ViscElCapMat(frictionAngle=fr,density=rho,Vb=VB3,gamma=Gamma,theta=Theta,Capillar=True,CapillarType=CapillarType,tc=tc,en=en,et=et))
+
+d = 0.9999
+
+id1 = O.bodies.append(sphere(center=[0,0,0], radius=r2,material=mat1,fixed=True, color=[0,1,0]))
+id2 = O.bodies.append(sphere(center=[0,(r1+r2)*d,0], radius=r2,material=mat1,fixed=True, color=[0,1,0]))
+id3 = O.bodies.append(sphere(center=[(r1+r2)*d,(r1+r2)*d,0], radius=r2,material=mat1,fixed=True, color=[0,1,0]))
+id4 = O.bodies.append(sphere(center=[(r1+r2)*d,(r1+r2)*d*2,0], radius=r2,material=mat1,fixed=True, color=[0,1,0]))
+id5 = O.bodies.append(sphere(center=[(r1+r2)*d*2,(r1+r2)*d,0], radius=r2,material=mat1,fixed=True, color=[0,1,0]))
+
+
+Vf = 0.5e-1
+Vfmin = 0.1e-1
+
+O.bodies[id1].Vf = Vf
+O.bodies[id1].Vmin = Vfmin
+
+O.bodies[id2].Vf = Vf
+O.bodies[id2].Vmin = Vfmin
+
+O.bodies[id3].Vf = Vf
+O.bodies[id3].Vmin = Vfmin
+
+O.bodies[id4].Vf = Vf
+O.bodies[id4].Vmin = Vfmin
+
+O.bodies[id5].Vf = Vf
+O.bodies[id5].Vmin = Vfmin
+
+vel = 0.0
+O.bodies[id1].state.vel=[0,0,vel]
+O.bodies[id2].state.vel=[0,0,-vel]
+O.bodies[id3].state.vel=[vel,0,0]
+O.bodies[id4].state.vel=[-vel,0,0]
+
+o.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Sphere_Aabb()],verletDist=(r1+r2)*5.0),
+ InteractionLoop(
+ [Ig2_Sphere_Sphere_ScGeom()],
+ [Ip2_ViscElCapMat_ViscElCapMat_ViscElCapPhys()],
+ [Law2_ScGeom_ViscElCapPhys_Basic()],
+ ),
+ LiqControl(),
+ NewtonIntegrator(damping=0,gravity=[0,0,0]),
+ PyRunner(command='showData()',iterPeriod=1),
+]
+
+
+def showData():
+ print "Step %d"%O.iter
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id1, O.bodies[id1].Vf, O.bodies[id1].Vmin)
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id2, O.bodies[id2].Vf, O.bodies[id2].Vmin)
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id3, O.bodies[id3].Vf, O.bodies[id3].Vmin)
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id4, O.bodies[id4].Vf, O.bodies[id4].Vmin)
+ print "idB=%d, Vf=%s, Vmin=%s;"%(id5, O.bodies[id5].Vf, O.bodies[id5].Vmin)
+
+ try:
+ print "Interaction[1, 2].Vb=%s"%(O.interactions[id1,id2].phys.Vb)
+ except:
+ pass
+
+ try:
+ print "Interaction[2, 3].Vb=%s"%(O.interactions[id2,id3].phys.Vb)
+ except:
+ pass
+
+ try:
+ print "Interaction[3, 4].Vb=%s"%(O.interactions[id3,id4].phys.Vb)
+ except:
+ pass
+
+ try:
+ print "Interaction[3, 5].Vb=%s"%(O.interactions[id3,id5].phys.Vb)
+ except:
+ pass
+ print
+
+
+showData()
+
+O.run(1, True)
+
+for i in range(5):
+ O.bodies[i].Vf = 0
+ O.bodies[i].Vmin = 0
+
+O.interactions[id1,id2].phys.Vb = 1.0
+O.interactions[id1,id2].phys.Vmax = 5.0
+O.interactions[id2,id3].phys.Vb = 1.0
+O.interactions[id2,id3].phys.Vmax = 5.0
+O.interactions[id3,id4].phys.Vb = 1.0
+O.interactions[id3,id4].phys.Vmax = 5.0
+O.interactions[id3,id5].phys.Vb = 1.0
+O.interactions[id3,id5].phys.Vmax = 5.0
+
+O.run(1, True)
+
+
+vel = 1.0
+O.bodies[id3].state.vel=[vel,0,0]
+O.bodies[id4].state.vel=[vel,0,0]
+O.bodies[id5].state.vel=[vel,0,0]
+
+from yade import qt
+qt.View()
=== modified file 'pkg/dem/VTKRecorder.cpp'
--- pkg/dem/VTKRecorder.cpp 2014-04-09 14:03:16 +0000
+++ pkg/dem/VTKRecorder.cpp 2014-05-08 05:57:04 +0000
@@ -30,7 +30,9 @@
#include<yade/pkg/dem/WirePM.hpp>
#include<yade/pkg/dem/JointedCohesiveFrictionalPM.hpp>
#include<yade/pkg/dem/Shop.hpp>
-
+#ifdef YADE_LIQMIGRATION
+ #include<yade/pkg/dem/ViscoelasticCapillarPM.hpp>
+#endif
YADE_PLUGIN((VTKRecorder));
CREATE_LOGGER(VTKRecorder);
@@ -70,6 +72,7 @@
else if(rec=="jcfpm") recActive[REC_JCFPM]=true;
else if(rec=="cracks") recActive[REC_CRACKS]=true;
else if(rec=="pericell" && scene->isPeriodic) recActive[REC_PERICELL]=true;
+ else if(rec=="liquidcontrol") recActive[REC_LIQ]=true;
else LOG_ERROR("Unknown recorder named `"<<rec<<"' (supported are: all, spheres, velocity, facets, boxes, color, stress, cpm, wpm, intr, id, clumpId, materialId, jcfpm, cracks, pericell). Ignored.");
}
// cpm needs interactions
@@ -80,7 +83,10 @@
// wpm needs interactions
if(recActive[REC_WPM]) recActive[REC_INTR]=true;
-
+
+ // liquid control needs interactions
+ if(recActive[REC_LIQ]) recActive[REC_INTR]=true;
+
// spheres
vtkSmartPointer<vtkPoints> spheresPos = vtkSmartPointer<vtkPoints>::New();
@@ -116,6 +122,20 @@
spheresCoordNumbSPH->SetName("SPH_Neigh");
#endif
+#ifdef YADE_LIQMIGRATION
+ vtkSmartPointer<vtkFloatArray> spheresLiqVol = vtkSmartPointer<vtkFloatArray>::New();
+ spheresLiqVol->SetNumberOfComponents(1);
+ spheresLiqVol->SetName("Liq_Vol");
+
+ vtkSmartPointer<vtkFloatArray> spheresLiqVolIter = vtkSmartPointer<vtkFloatArray>::New();
+ spheresLiqVolIter->SetNumberOfComponents(1);
+ spheresLiqVolIter->SetName("Liq_VolIter");
+
+ vtkSmartPointer<vtkFloatArray> spheresLiqVolTotal = vtkSmartPointer<vtkFloatArray>::New();
+ spheresLiqVolTotal->SetNumberOfComponents(1);
+ spheresLiqVolTotal->SetName("Liq_VolTotal");
+#endif
+
vtkSmartPointer<vtkFloatArray> spheresMask = vtkSmartPointer<vtkFloatArray>::New();
spheresMask->SetNumberOfComponents(1);
spheresMask->SetName("mask");
@@ -258,7 +278,16 @@
vtkSmartPointer<vtkFloatArray> crackNorm = vtkSmartPointer<vtkFloatArray>::New();
crackNorm->SetNumberOfComponents(3);
crackNorm->SetName("crackNorm");
-
+
+#ifdef YADE_LIQMIGRATION
+ vtkSmartPointer<vtkFloatArray> liqVol = vtkSmartPointer<vtkFloatArray>::New();
+ liqVol->SetNumberOfComponents(1);
+ liqVol->SetName("liqVol");
+
+ vtkSmartPointer<vtkFloatArray> liqVolNorm = vtkSmartPointer<vtkFloatArray>::New();
+ liqVolNorm->SetNumberOfComponents(1);
+ liqVolNorm->SetName("liqVolNorm");
+#endif
// the same for newly created cracks
// vtkSmartPointer<vtkPoints> crackPosNew = vtkSmartPointer<vtkPoints>::New();
// vtkSmartPointer<vtkCellArray> crackCellsNew = vtkSmartPointer<vtkCellArray>::New();
@@ -365,11 +394,16 @@
intrIsCohesive->InsertNextValue(jcfpmphys->isCohesive);
intrIsOnJoint->InsertNextValue(jcfpmphys->isOnJoint);
intrForceN->InsertNextValue(fn);
- }
-
- else {
+ } else {
intrForceN->InsertNextValue(fn);
}
+#ifdef YADE_LIQMIGRATION
+ if (recActive[REC_LIQ]) {
+ const ViscElCapPhys* capphys = YADE_CAST<ViscElCapPhys*>(I->phys.get());
+ liqVol->InsertNextValue(capphys->Vb);
+ liqVolNorm->InsertNextValue(capphys->Vb/capphys->Vmax);
+ }
+#endif
}
}
}
@@ -438,6 +472,12 @@
spheresPressSPH->InsertNextValue(b->press);
spheresCoordNumbSPH->InsertNextValue(b->coordNumber());
#endif
+#ifdef YADE_LIQMIGRATION
+ spheresLiqVol->InsertNextValue(b->Vf);
+ const Real tmpVolIter = liqVolIterBody(b);
+ spheresLiqVolIter->InsertNextValue(tmpVolIter);
+ spheresLiqVolTotal->InsertNextValue(tmpVolIter + b->Vf);
+#endif
if (recActive[REC_MATERIALID]) spheresMaterialId->InsertNextValue(b->material->id);
continue;
}
@@ -579,6 +619,11 @@
spheresUg->GetPointData()->AddArray(spheresPressSPH);
spheresUg->GetPointData()->AddArray(spheresCoordNumbSPH);
#endif
+#ifdef YADE_LIQMIGRATION
+ spheresUg->GetPointData()->AddArray(spheresLiqVol);
+ spheresUg->GetPointData()->AddArray(spheresLiqVolIter);
+ spheresUg->GetPointData()->AddArray(spheresLiqVolTotal);
+#endif
if (recActive[REC_STRESS]){
spheresUg->GetPointData()->AddArray(spheresNormalStressVec);
spheresUg->GetPointData()->AddArray(spheresShearStressVec);
@@ -592,6 +637,7 @@
if (recActive[REC_JCFPM]) {
spheresUg->GetPointData()->AddArray(damage);
}
+
if (recActive[REC_MATERIALID]) spheresUg->GetPointData()->AddArray(spheresMaterialId);
#ifdef YADE_VTK_MULTIBLOCK
@@ -673,6 +719,12 @@
intrPd->SetLines(intrCells);
intrPd->GetCellData()->AddArray(intrForceN);
intrPd->GetCellData()->AddArray(intrAbsForceT);
+#ifdef YADE_LIQMIGRATION
+ if (recActive[REC_LIQ]) {
+ intrPd->GetCellData()->AddArray(liqVol);
+ intrPd->GetCellData()->AddArray(liqVolNorm);
+ }
+#endif
if (recActive[REC_JCFPM]) {
intrPd->GetCellData()->AddArray(intrIsCohesive);
intrPd->GetCellData()->AddArray(intrIsOnJoint);
=== modified file 'pkg/dem/VTKRecorder.hpp'
--- pkg/dem/VTKRecorder.hpp 2014-01-16 14:57:09 +0000
+++ pkg/dem/VTKRecorder.hpp 2014-05-08 05:57:04 +0000
@@ -11,7 +11,7 @@
class VTKRecorder: public PeriodicEngine {
public:
- enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_JCFPM,REC_CRACKS,REC_WPM,REC_PERICELL};
+ enum {REC_SPHERES=0,REC_FACETS,REC_BOXES,REC_COLORS,REC_MASS,REC_CPM,REC_INTR,REC_VELOCITY,REC_ID,REC_CLUMPID,REC_SENTINEL,REC_MATERIALID,REC_STRESS,REC_MASK,REC_RPM,REC_JCFPM,REC_CRACKS,REC_WPM,REC_PERICELL,REC_LIQ};
virtual void action();
void addWallVTK (vtkSmartPointer<vtkQuad>& boxes, vtkSmartPointer<vtkPoints>& boxesPos, Vector3r& W1, Vector3r& W2, Vector3r& W3, Vector3r& W4);
YADE_CLASS_BASE_DOC_ATTRS_CTOR(VTKRecorder,PeriodicEngine,"Engine recording snapshots of simulation into series of \\*.vtu files, readable by VTK-based postprocessing programs such as Paraview. Both bodies (spheres and facets) and interactions can be recorded, with various vector/scalar quantities that are defined on them.\n\n:yref:`PeriodicEngine.initRun` is initialized to ``True`` automatically.",
=== modified file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 2014-04-14 14:44:08 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-05-08 05:57:04 +0000
@@ -83,6 +83,9 @@
if (not(phys.liqBridgeCreated) and phys.Capillar and geom.penetrationDepth>=0) {
phys.liqBridgeCreated = true;
phys.liqBridgeActive = false;
+ #ifdef YADE_LIQMIGRATION
+ scene->addIntrs.push_back(I);
+ #endif
Sphere* s1=dynamic_cast<Sphere*>(bodies[id1]->shape.get());
Sphere* s2=dynamic_cast<Sphere*>(bodies[id2]->shape.get());
if (s1 and s2) {
@@ -114,6 +117,12 @@
VLiqBridg -= phys.Vb;
NLiqBridg -= 1;
}
+ #ifdef YADE_LIQMIGRATION
+ const intReal B1={id1, phys.Vb/2.0};
+ const intReal B2={id2, phys.Vb/2.0};
+ scene->delIntrs.push_back(B1);
+ scene->delIntrs.push_back(B2);
+ #endif
scene->interactions->requestErase(I);
return;
};
@@ -319,14 +328,13 @@
* Please, use this model only for testing purposes.
*
*/
-
+
const Real R = phys.R;
const Real Gamma = phys.gamma;
const Real D = -geom.penetrationDepth;
const Real V = phys.Vb;
const Real Theta = phys.theta;
-
const Real a = -1.1*pow((V/(R*R*R)), -0.53);
const Real b = (-0.148*log(V/(R*R*R)) - 0.96)*Theta*Theta -0.0082*log(V/(R*R*R)) + 0.48;
const Real c = 0.0018*log(V/(R*R*R)) + 0.078;
@@ -339,3 +347,199 @@
Real Law2_ScGeom_ViscElCapPhys_Basic::None_f(const ScGeom& geom, ViscElCapPhys& phys) {
return 0;
}
+
+#ifdef YADE_LIQMIGRATION
+YADE_PLUGIN((LiqControl));
+void LiqControl::action(){
+ // This function implements liquid migration model, introduced here [Mani2013]
+ mapBodyInt bI;
+ mapBodyInt bodyNeedUpdate;
+
+ // Calculate, how much new contacts will be at each body
+ for (unsigned int i=0; i<scene->addIntrs.size(); i++) {
+ addBodyMapInt( bI, scene->addIntrs[i]->getId1() );
+ addBodyMapInt( bI, scene->addIntrs[i]->getId2() );
+ }
+
+ // Update volume water at each deleted interaction for each body
+ for (unsigned int i=0; i<scene->delIntrs.size(); i++) {
+ shared_ptr<Body> b = Body::byId(scene->delIntrs[i].id,scene);
+ b->Vf += scene->delIntrs[i].Vol;
+ addBodyMapInt(bodyNeedUpdate, scene->delIntrs[i].id);
+ liqVolRup += scene->delIntrs[i].Vol;
+ }
+ scene->delIntrs.clear();
+
+ // Update volume bridge at each new added interaction
+ mapBodyReal bodyUpdateLiquid;
+ for (unsigned int i=0; i<scene->addIntrs.size(); i++) {
+ shared_ptr<Body> b1 = Body::byId(scene->addIntrs[i]->getId1(),scene);
+ shared_ptr<Body> b2 = Body::byId(scene->addIntrs[i]->getId2(),scene);
+
+ const id_t id1 = b1->id;
+ const id_t id2 = b2->id;
+
+ ViscElCapPhys* Vb=dynamic_cast<ViscElCapPhys*>(scene->addIntrs[i]->phys.get());
+ const Real Vmax = vMax(b1, b2);
+ Vb->Vmax = Vmax;
+
+ Real Vf1 = 0.0;
+ Real Vf2 = 0.0;
+
+ if ((b1->Vmin)<b1->Vf) {
+ Vf1 = (b1->Vf - b1->Vmin)/bI[id1];
+ }
+
+ if ((b2->Vmin)<b2->Vf) {
+ Vf2 = (b2->Vf - b2->Vmin)/bI[id2];
+ }
+
+ Real Vrup = Vf1+Vf2;
+
+ if(mask!=0 && ((b1->groupMask & b2->groupMask & mask)==0)) {
+ Vf1 = 0;
+ Vf2 = 0;
+ Vrup = 0;
+ } else if (Vrup > Vmax) {
+ Vf1 *= Vmax/Vrup;
+ Vf2 *= Vmax/Vrup;
+ Vrup = Vf1 + Vf2;
+ }
+
+ liqVolShr += Vrup;
+ addBodyMapReal(bodyUpdateLiquid, id1, -Vf1);
+ addBodyMapReal(bodyUpdateLiquid, id2, -Vf2);
+
+ Vb->Vb = Vrup;
+ }
+
+ scene->addIntrs.clear();
+
+ // Update water volume in body
+ for (mapBodyReal::const_iterator it = bodyUpdateLiquid.begin(); it != bodyUpdateLiquid.end(); ++it) {
+ Body::byId(it->first)->Vf += it->second;
+ }
+
+ // Update contacts around body
+ for (mapBodyInt::const_iterator it = bodyNeedUpdate.begin(); it != bodyNeedUpdate.end(); ++it) {
+ updateLiquid(Body::byId(it->first));
+ }
+
+}
+
+void LiqControl::updateLiquid(shared_ptr<Body> b){
+ if (b->Vf<=b->Vmin) {
+ return;
+ } else {
+ // How much liquid can body share
+ const Real LiqCanBeShared = b->Vf - b->Vmin;
+
+ // Check how much liquid can accept contacts
+ Real LiqContactsAccept = 0.0;
+ unsigned int contactN = 0;
+ for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
+ if(!((*it).second) or !(((*it).second)->isReal())) continue;
+ ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
+ if (physT->Vb<physT->Vmax) {
+ LiqContactsAccept+=physT->Vmax-physT->Vb;
+ contactN++;
+ }
+ }
+
+ if (contactN>0) {
+ //There are some contacts, which can be filled
+ Real FillLevel = 0.0;
+ if (LiqContactsAccept > LiqCanBeShared) { // Share all available liquid from body to contacts
+ const Real LiquidWillBeShared = b->Vf - b->Vmin;
+ b->Vf = b->Vmin;
+ FillLevel = LiquidWillBeShared/LiqContactsAccept;
+ } else { // Not all available liquid from body can be shared
+ b->Vf -= LiqContactsAccept;
+ FillLevel = 1.0;
+ }
+
+ for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
+ if(!((*it).second) or !(((*it).second)->isReal())) continue;
+ ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
+ if (physT->Vb<physT->Vmax) {
+ liqVolShr += (physT->Vmax - physT->Vb)*FillLevel;
+ physT->Vb += (physT->Vmax - physT->Vb)*FillLevel;
+ }
+ }
+ return;
+ } else {
+ return;
+ }
+ }
+}
+
+void LiqControl::addBodyMapInt( mapBodyInt & m, Body::id_t b ){
+ mapBodyInt::const_iterator got;
+ got = m.find (b);
+ if ( got == m.end() ) {
+ m.insert (mapBodyInt::value_type(b,1));
+ } else {
+ m[b] += 1;
+ }
+}
+
+void LiqControl::addBodyMapReal( mapBodyReal & m, Body::id_t b, Real addV ) {
+ mapBodyReal::const_iterator got;
+ got = m.find (b);
+ if ( got == m.end() ) {
+ m.insert (mapBodyReal::value_type(b, addV));
+ } else {
+ m[b] += addV;
+ }
+}
+
+Real LiqControl::vMax(shared_ptr<Body> const b1, shared_ptr<Body> const b2) {
+ Sphere* s1=dynamic_cast<Sphere*>(b1->shape.get());
+ Sphere* s2=dynamic_cast<Sphere*>(b2->shape.get());
+ Real minR = 0.0;
+ if (s1 and s2) {
+ minR = std::min (s1->radius, s2->radius);
+ } else if (s1 and not(s2)) {
+ minR = s1->radius;
+ } else {
+ minR = s2->radius;
+ }
+ return vMaxCoef*minR*minR*minR;
+}
+
+Real liqVolIterBody (shared_ptr<Body> b) {
+ Real LiqVol = 0.0;
+ for(Body::MapId2IntrT::iterator it=b->intrs.begin(),end=b->intrs.end(); it!=end; ++it) {
+ if(!((*it).second) or !(((*it).second)->isReal())) continue;
+ ViscElCapPhys* physT=dynamic_cast<ViscElCapPhys*>(((*it).second)->phys.get());
+ if (physT->Vb>0) {
+ LiqVol += physT->Vb/2.0;
+ }
+ }
+ return LiqVol;
+}
+
+Real LiqControl::liqVolBody (id_t id) const {
+ Scene* scene=Omega::instance().getScene().get();
+ const BodyContainer& bodies = *scene->bodies;
+ if (id >=0 and bodies[id]) {
+ if (bodies[id]->Vf > 0) {
+ return bodies[id]->Vf + liqVolIterBody(bodies[id]);
+ } else {
+ return liqVolIterBody(bodies[id]);
+ }
+ }
+ else return -1;
+}
+
+Real LiqControl::totalLiqVol(int mask=0) const{
+ Scene* scene=Omega::instance().getScene().get();
+ Real totalLiqVol = 0.0;
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if((mask>0 && (b->groupMask & mask)==0) or (!b)) continue;
+ totalLiqVol += liqVolIterBody(b);
+ if (b->Vf > 0) {totalLiqVol +=b->Vf;}
+ }
+ return totalLiqVol;
+}
+#endif
=== modified file 'pkg/dem/ViscoelasticCapillarPM.hpp'
--- pkg/dem/ViscoelasticCapillarPM.hpp 2014-04-14 14:11:53 +0000
+++ pkg/dem/ViscoelasticCapillarPM.hpp 2014-05-08 05:57:04 +0000
@@ -1,4 +1,5 @@
#include"ViscoelasticPM.hpp"
+#include <boost/unordered_map.hpp>
class ViscElCapMat : public ViscElMat {
public:
@@ -31,7 +32,11 @@
((Real,Vb,NaN,,"Liquid bridge volume [m^3]"))
((Real,gamma,NaN,,"Surface tension [N/m]"))
((Real,theta,NaN,,"Contact angle [rad]"))
- ((CapType,CapillarType,None_Capillar,,"Different types of capillar interaction: Willett_numeric, Willett_analytic, Weigert, Rabinovich, Lambert, Soulie")),
+ ((CapType,CapillarType,None_Capillar,,"Different types of capillar interaction: Willett_numeric, Willett_analytic, Weigert, Rabinovich, Lambert, Soulie"))
+#ifdef YADE_LIQMIGRATION
+ ((Real,Vmax,-1,,"Maximal liquid bridge volume [m^3]"))
+#endif
+ ,
createIndex();
)
REGISTER_CLASS_INDEX(ViscElCapPhys,ViscElPhys);
@@ -72,3 +77,31 @@
DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(Law2_ScGeom_ViscElCapPhys_Basic);
+
+#ifdef YADE_LIQMIGRATION
+typedef boost::unordered_map<Body::id_t, int> mapBodyInt;
+typedef boost::unordered_map<Body::id_t, Real> mapBodyReal;
+class LiqControl: public PartialEngine{
+ public:
+ virtual void action();
+ void addBodyMapInt( mapBodyInt & m, Body::id_t b );
+ void addBodyMapReal( mapBodyReal & m, Body::id_t b, Real addV );
+ Real vMax(shared_ptr<Body> b1, shared_ptr<Body> b2);
+ Real totalLiqVol(int mask) const;
+ Real liqVolBody(id_t id) const;
+ void updateLiquid(shared_ptr<Body> b);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(LiqControl,PartialEngine,"This engine implements liquid migration model, introduced here [Mani2013]_ . ",
+ ((int,mask,0,, "Bitmask for liquid creation."))
+ ((Real,liqVolRup,0.,, "Liquid volume (integral value), which has been freed after rupture occured, [m^3]."))
+ ((Real,liqVolShr,0.,, "Liquid volume (integral value), which has been shared among of contacts, [m^3]."))
+ ((Real,vMaxCoef,0.03,, "Coefficient for vMax, [-]."))
+ ,/* ctor */
+ ,/* py */
+ .def("totalLiq",&LiqControl::totalLiqVol,(python::arg("mask")=0),"Return total volume of water in simulation.")
+ .def("liqBody",&LiqControl::liqVolBody,(python::arg("id")=-1),"Return total volume of water in body.")
+ );
+};
+
+Real liqVolIterBody (shared_ptr<Body> b);
+REGISTER_SERIALIZABLE(LiqControl);
+#endif
